Three-dimensional Shape Explains Star Formation Mystery of California and Orion A

Plane-of-the-sky studies have provided valuable information about the star-forming interstellar medium (ISM). However, our knowledge of the true, three-dimensional (3D) structure of molecular clouds remains limited, which hampers our understanding of how the clouds evolve and form stars. But now, thanks to the unprecedented astrometry from the Gaia mission (Gaia Collaboration et al. 2016), studying the 3D structure of molecular clouds has become possible (e.g., Großschedl et al. 2018; Leike et al. 2020; Rezaei et al. 2020; Zucker et al. 2020; Kainulainen et al. 2022). In this Letter, we compare the 3D shapes of two molecular clouds and demonstrate how that knowledge has a fundamental effect in understanding their star formation activities.

The California and Orion A molecular clouds are two of the most massive giant molecular clouds (GMCs) within 500 pc from the Sun. The two clouds are located relatively near one another and show comparable kinematics and similar filamentary shapes and sizes in the plane of the sky (Lada et al. 2009). However, their star formation activities differ significantly, with Orion A having an order of magnitude higher star formation rate than California (Lada et al. 2009). The difference triggers an immediate question about what regulates star formation and what role the individual cloud's properties like shape and size play in setting its star formation activity.

To deepen our understanding of the contrasting star formation rates in California and Orion A, we investigate their 3D shape by mapping the dust distributions toward the two clouds. The detailed 3D dust distribution toward Orion A was already presented in Rezaei et al. (2020) where we revealed new information on its 3D shape, discovered a foreground dust component, and demonstrated the extending tail of the cloud to further distances as suggested by Großschedl et al. (2018). Here, we look into the 3D distribution of the dust toward California to map its 3D substructures that remain hidden in the plane-of-the-sky data.

While there are a couple of local dust maps of the Milky Way, our current work toward California and Orion A is a unique map of its kind. Leike et al. (2020) provide a 3D dust map of the solar neighborhood using Gaia DR2 and metric Gaussian variational inference that is then used by Zucker et al. (2021) to define cloud boundaries. However, Leike et al. (2020) provides data for all lines of sight (LOSs) only up to 370 pc and for some LOSs close to 500 pc, depending on the coordinates. Their results beyond 370 pc need to be treated with extreme caution, as it reaches the boundaries of the map that affects the density, shape, and distance to the clouds (Zucker et al. 2021). Apart from the coverage limits of Leike et al. (2020), our 3D dust mapping technique has the advantage of taking into account distance uncertainties; as a result, it can exploit a full data set like Gaia without the necessity to cut on noisy data. Those stars are typically located in dusty regions, which are important for mapping the dense parts of the molecular clouds. Moreover, owing to the analytical solution for our posterior calculation, approximations are minimal; therefore, our method provides accurate and reliable estimates (see Rezaei et al. 2017, 2018b, 2020 for more details).

Our 3D dust mapping technique has been extensively explained in Rezaei et al. (2017, 2018b). In addition, new developments have been introduced in Rezaei et al. (2020). Here we briefly summarize the main aspects of our method.

Our technique consists of a nonparametric method that uses 3D positions of stars together with their LOS attenuation as the input data and predicts dust densities for arbitrary points in the same 3D space. We divide the LOS toward stars into 1D cells in order to model the integrated attenuation to each star as sum of the dust densities along its LOS by connecting all 1D cells using a Gaussian process (GP) prior. GP takes into account the neighboring correlations using a correlation length, λ. Another hyperparameter of the model is the variance, θ, which sets the amplitude of density variations. All hyperparameters are fixed based on the input data as explained in Rezaei et al. (2017, 2018b). The final resolution of the map is set by the typical separation between input stars. Apart from considering the 3D neighboring correlation, our method has the advantage of taking into account both distance and extinction uncertainties toward individual stars that result in producing robust estimates of the distance and 3D shapes of the molecular clouds. Despite common artifacts in most of the 3D dust maps, our results are devoid of discontinuity and "fingers-of-God" artifact.

The early installment of the third Gaia data release (Gaia EDR3) consists of the full astrometric solution for around 1.5 billion sources, and a significant advantage over Gaia DR2 parallaxes with 30% increase in precision (Gaia Collaboration et al. 2021). We use the 3D positions of stars from Gaia EDR3; only limiting our sample to parallax uncertainty of less than 100%. Since we are dealing with large parallax uncertainties, we use the geometric distance estimates from Bailer-Jones et al. (2021) to map the full shape of the cloud.

Similar to Rezaei et al. (2020), we use the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) and the Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010) photometry to calculate extinctions to individual stars using the Rayleigh–Jeans Color Excess (RJCE; Majewski et al. 2011) method, which provides us with the extinction in the K_S band. The crossmatch between 2MASS and WISE catalogs with Gaia EDR3 sources are provided on Gaia Archive. ³

After calculating extinctions, we select our final sample based on the position of stars on the dereddened color–magnitude diagram in order to remove the outliers (see Rezaei et al. 2018a, 2020). For our current work toward the California cloud, 155° < l < 170° and −14° < b < −6°, we use around 160,000 stars as our final input sample. The hyperparameters of the method for the aforementioned data are cell size'= 5 pc, λ = 20 pc, and θ = 4 × 10⁻⁸ pc⁻². It is important to note that the choice of hyperparameters is not completely arbitrary. As discussed in detail in Rezaei et al. (2017, 2018b), the hyperparameters are calculated according to input data. The final resolution of the map is not set by these parameters, but rather the input data (see Appendix B for an example).

Figure 1 shows the 3D distribution of dust toward the California region from two viewing angles. Rather than being a filament perpendicular to the LOS, California is an extended flat, sheet-like structure, with a bubble at one side, where higher-density structures form within the sheet. From the face-on view, the longest length of the sheet is approximately 120 pc, with the shortest length of about 80 pc, while from the edge-on view, the width of the sheet appears very narrow (as expected from the plane-of-the-sky observations). Figure 2 demonstrates various components of the cloud at multiple distance slices: while the bubble centered at l = 158°, b = −10° appears clearly at 455 pc and expands over only 25 pc in radius, some other parts of the cloud extend along the sheet for over 100 pc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Having uncertainties on the dust density predictions is crucial for evaluating the significance and validity of different 3D substructures, as explained in Rezaei et al. (2020). Figure 3 shows dust density predictions with their uncertainties as a function of distance for three LOSs toward different parts of the cloud (also marked in Figure 2 in the plane of the sky). The red curve that is toward the lowest part of the bubble has a dominant, symmetrical overdensity at 445 pc, and there is a shallower overdensity in the foreground at 290 pc. The foreground overdensity is in agreement with an increase in reddening in the Green et al. (2019) map at the same distance, and an increase in density in Leike et al. (2020). The blue and gray curves, however, recover a broader distribution that represents the extent of the sheet along the LOS. While the highest peak of dust density along (l, b) = (161, −9) (blue curve) is around 450 pc (as estimated by Zucker et al. 2020), it captures the most elongated part of the sheet extending from 410 to 530 pc, giving it the length of about 120 pc. There is also a shallow overdensity at closer distances toward this LOS that is related to the Taurus cloud on the foreground. As seen from the gray curve, (l, b) = (165, −8.5), the peak density of the eastern part of California is at a farther distance of 514 pc.

Zoom In Zoom Out Reset image size

It is important to note that the extension of the clouds along various LOSs are driven by the input data; in fact, as explained in our previous works, our method is capable of capturing multiscale density substructures that are smaller or larger than the input correlation length. In addition, the peak densities presented here are not representative of the high-density small-scale molecular gas but an average density within the resolution of the map; therefore, much higher density cores can be located within our dust clumps. We demonstrate these points with a mock data set in Appendix B.

We can reconstruct the 2D extinction map toward California by integrating our predicted 3D dust densities along various LOSs. As can be seen from Figure 7, the final 2D projected map reaches extinction of ∼0.5 mag in the K_s band (∼5 mag in A_V ), and there is a clear resemblance to higher-resolution maps (e.g., CO-based maps or NICEST/NICER extinction maps).

Even though very similar in the plane of the sky (Lada et al. 2009), California's 3D structure appears very different from that of Orion A. Within our resolution limit, only one apparent feedback-driven substructure is seen in California (a bubble on the western side), while Orion A has several such structures. Beyond this substructure, California is distributed along a flat sheet, while Orion A is more filamentary also in 3D. The two clouds are shown together in Figure 4. The substructures in Orion A are more complex than in California, likely caused by feedback processes and multiple episodes of star formation in the region (e.g., Schlafly et al. 2015; Rezaei et al. 2020). In addition, substructures in Orion A have higher densities than in California (see Figure 4). Since the shallow overdensities in California are spread along the flat sheet, they add up in the plane of the sky to represent extinction properties similar to those of Orion A. However, the 3D view reveals the dramatic differences between the clouds.

Zoom In Zoom Out Reset image size

In addition, other studies have estimated the age of California to be around 1–3 Myr (e.g., Covey et al. 2010; Wolk et al. 2010; Imara et al. 2017), while Orion A seems to be in an older evolutionary stage with older stellar populations (5–10 Myr; e.g., Bouy et al. 2014; Zari et al. 2019) associated with its discovered foreground bubble (Rezaei et al. 2020), which seems to have triggered the next generation of star formation in the main Orion A filament in the background.

We also demonstrate that the knowledge of clouds' 3D shape is an important factor in interpreting their star formation activities and in considering star formation relations such as the Kennicutt–Schmidt (KS) relation. Without observational biases, the mean surface density of the California cloud changes drastically depending on the viewing angle: the total mass of the cloud is the same irrespective of the viewing angle, but the surface area is an order of magnitude larger when viewed from the direction perpendicular to the Galactic plane (perpendicular to the cloud's sheet-like morphology). The change in area is less dramatic for Orion A that is more filamentary than sheet-like. This effect has two potentially important consequences:

First, the location of California in the KS relation depends significantly on the viewing angle. The star formation rate (SFR) in a molecular cloud in the Milky Way is calculated from the number of young stellar objects (YSOs) in the region, assuming an average mass and age for the YSOs. Since the number of observed YSOs does not change with the angle, then SFR will remain constant. As the cloud area changes, but the mass and SFR remain constant, California can move in the relation along a line with the slope of unity (see Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

How exactly the different structures of the clouds affect their location in the observed KS relation depends on the details of the observational techniques used to derive the masses and areas of the clouds. In Kainulainen et al. (2022), we showed that in particular the column density threshold used to define the clouds has a strong effect on this. In particular, using plane-of-the-sky data, Lada et al. (2013) found that there is no KS relation between the clouds in the solar neighborhood. Our result opens a door to speculate that the relation may well exist, but, depending on the size and morphology of the clouds, is not recovered from the plane-of-the-sky data. This could be either due to systematic effects in the cloud orientations, coincidence in the sampling of random orientation angles, or the details of the cloud definition process (especially thresholding; Kainulainen et al. 2022).

Second, our result also suggests a link between the 3D morphology of molecular clouds and their star formation activity. The plane-of-the-sky works indicate that the California and Orion A clouds are morphologically similar; the origin of the differing star formation activities has been attributed to differing amounts of dense gas (e.g., Kainulainen et al. 2009; Lada et al. 2013; Kainulainen & Federrath 2017) and/or environmental effects (e.g., Schlafly et al. 2015). Our work shows that the two clouds are not morphologically similar, which immediately raises the question about the importance of morphology. Again, providing 3D shape information for a larger sample of clouds will be vital to establish this.

Finally, we provide the most up-to-date distance estimates to various LOSs toward California and Orion A using Gaia EDR3 (Table 1). We note that the geometric distance estimates from GEDR3 parallaxes by Bailer-Jones et al. (2021) are overestimated because their prior's scale length is larger than our targeted cloud distance closer than 500 pc. Therefore, after we map the full 3D shape of the cloud using all stars that match our input criteria (see Section 2), in order to get accurate distances to individual dust components along various LOSs, we only use stars from Gaia EDR3 with parallax estimates better than 20% and directly invert them to get their distances.

There is a perfect agreement between our distance estimates and distance estimates of Zucker et al. (2020) for LOSs toward where there is a prominent peak. However, where there are multiple peaks near each other along an LOS (like the foreground cloud in Orion A) or where there are no dominant peaks (e.g., elongations in California), our distance estimates differ. Although distances in Zucker et al. (2020) are predicted on higher spatial resolution than that of our work, it is possible that due to the nature of the Zucker et al. (2020) technique, multiple nearby peaks and the physical shape of the clouds remain undetected in their work. It is important to note that, as explained in Rezaei et al. (2020), our predictions are made for fixed points in 3D; therefore, there is no direct uncertainty associated with our distances. However, our method provides the full probability distribution function for density predictions at each point that can then be propagated into distance uncertainties. Having done so, the typical uncertainty in our predicted distances is 10 pc.

We have determined the 3D structure of the California molecular cloud using our advanced mapping technique that takes into account distance and extinction uncertainties to individual stars and considers correlations between neighboring points in 3D. Having the full probability density function for each predicted point, i.e., mean and standard deviation of the predicted density, allows us to verify the results; therefore, our final products are robust 3D dust maps of the molecular clouds.

We demonstrated that California is a sheet-like structure extended around 120 pc along the LOS. It consists of low-density substructures that align with the LOS resulting in higher column density in the plane of the sky, comparable to that of Orion A. The surface density of California, however, is much lower than Orion A from the viewing angle perpendicular to the Galactic plane. In addition, the substructures of the two clouds show substantial differences with Orion A being a filamentary structure with several feedback bubbles, while California has only one isolated bubble on one side and largely resembles an unperturbed sheet. The two clouds, even though similar in the plane of the sky, are remarkably different in 3D. The dramatic differences likely contribute to the different star formation efficiencies of the clouds. Our results demonstrate the importance of 3D information in understanding the star formation activities of molecular clouds.

We also showed that the 3D information plays a dramatic role in setting the mean surface density of the clouds, which can then significantly affect their appearance in scaling relations such as their location in the KS relation. This further indicates that the column density thresholds are not necessarily reliable tools to determine core star formation activities in molecular clouds. We have presented the results for California in this Letter, but plan to expand this study to more clouds in order to have a better understanding of star formation in the solar neighborhood and its connection to extragalactic studies.

We would like to sincerely thank the anonymous referee for their thoughtful and constructive comments that improved the quality of the paper. S.R. thanks Henrik Beuther for fruitful discussions. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. J.K. has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 639459 (PROMISE).

Figure 6 shows the 3D distribution of California from different viewing angles which better demonstrates the full structure of the cloud. Figure 7 shows the reconstructed extinctions (${A}_{{K}_{s}}$) using the predicted 3D densities. For each LOS the predicted densities are summed up and multiplied by the cell sizes in order to get the projected extinctions. Despite the low resolution, the structures resemble that of the higher-resolution maps (Lada et al. 2009). Contours of ${A}_{{K}_{s}}=0.3$ mag from Green et al. (2019) are shown in the figure.

Figure 8 shows the reconstructed extinctions (${A}_{{K}_{s}}$) using the predicted 3D densities as a function of distance for three LOSs shown in Figure 3, together with those from the Green et al. (2019) map. The slope of the reddening curves represents the amount of dust density in the corresponding distance. Therefore, an increase in the reddening in Figure 8 is equivalent to a peak in our LOS density profile in Figure 3. Overall, there is a fairly good agreement between our predicted distances to the clouds and those from Green et al. (2019); the only difference is for the red curve, where the secondary peak in our map is not recovered by Green et al. (2019).

Apart from the predicted distances, predicted extinction in Green et al. (2019) seems to appear in steps; allocating almost all predicted dust densities toward a cloud to a certain distance, while an increase in extinction in our map is more gradual. This could be the result of different techniques used in each work: the method used in Green et al. (2019) is not optimal for predicting multiple overdensities near one another or recovering an extended structure. Our technique, on the other hand, tends to smooth the predictions, therefore underestimating the maximum predicted density for each peak. It is important to note, however, that the amount of elongation caused by the smoothing in our method is less than the 10 pc uncertainties; thus, the predicted elongated structures in our map are real effects of the data rather than an artifact caused by the technique (see Appendix B for an example using mock data).

Here we demonstrate the ability of the method to capture the locations and widths of the clouds quite precisely. Also, we demonstrate how the peak density predictions in the map are affected by its resolution.

We simulate mock observations using two dust clouds with different sizes and densities, and stars with similar spatial density distribution as that of Gaia. Here are the true properties of the simulated clouds:

1. Cloud 1 at 450 pc, 5 pc width, density of 2 × 10⁻³ mag pc⁻¹,

2. Cloud 2 at 650 pc, 70 pc width, density of 1 × 10⁻³ mag pc⁻¹.

We then calculate the LOS extinction to stars in this region and add 20% noise to both distances and extinctions in order to have a more realistic simulation. We use these as the input data and use the model to find the underlying 3D dust densities.

As can be seen in Figure 9, the model can predict the location of the clouds and their widths quite precisely. The predicted widths (elongation) of the clouds at the half-maximum are 15 pc for the first one and 85 pc for the second one, which is consistent with their true widths within ±10 pc uncertainties. It is worth mentioning that even though the correlation length of the model was 20 pc, it recovered clouds with completely different lengths from 5 to 70 pc; indicating that the elongation of the cloud is mainly determined by the input data. Another important point about the predictions are the maximum densities: since the size of the first cloud is smaller than the scale lengths of the model, the maximum predicted density is less than half of the "true" input density. These points indicate the following: (1) Our predicted distances are quite precise. (2) The extension of the clouds is reliable and determined by the input stars (with the predicted uncertainties). (3) The maximum densities probed by the map is an average density within the resolution of the map. The black line is the mean and the blue shaded area is the standard deviation of the predictions.