© J. D. Soler 2019

1 Introduction

The physical processes that regulate the formation of stars are one of the main open subjects of research in contemporary astrophysics (for a review see McKee & Ostriker 2007; Dobbs et al. 2014; Molinari et al. 2014). Interstellar magnetic fields are a long-standing candidate for the control of the star formation process (see, e.g., Mouschovias et al. 2006). Observations indicate that in flux-freezing conditions, the interstellar magnetic fields, gas pressure, cosmic-ray pressure, and gravity are the most important forces on the diffuse gas (Heiles & Crutcher 2005). It is likely that magnetic fields are responsible for reducing the star formation rate in dense molecular clouds (MCs) by a factor of a few by strongly shaping the interstellar gas (Hennebelle & Inutsuka 2019, and references therein). However, observing the magnetic fields in and around MCs and determining their exact influence is a difficult task.

Observations of the Zeeman splitting in emission and absorption by species such as neutral hydrogen (H) and the hydroxyl (OH) and cyanide (CN) radicals provide the only direct detections of the magnetic field strength in the stellar medium (see Crutcher 2012, for a review). The requirements of both high sensitivity and control of systematic effects considerably limit the coverage of the observations of the Zeeman splitting in MCs, however (Bourke et al. 2001; Troland & Crutcher 2008). Further information on the magnetic fields in MCs is provided by observations of the linear polarization, in extinction from background stars and in emission from interstellar dust grains, reveal the orientation of the interstellar magnetic field averaged along the line of sight (LOS) and projected onto the plane of the sky, B_⊥, (Hiltner 1949; Davis & Greenstein 1951; Hildebrand 1988; Pattle & Fissel 2019).

The magnetic field strength can be assessed indirectly from the linear polarization observations using the Davis– Chandrasekhar–Fermi method (DCF, Davis 1951; Chandrasekhar & Fermi 1953), which combines estimations of the density, velocity dispersion, and polarization angle dispersion toward the studied regions (see, e.g., Houde et al. 2009; Pattle et al. 2017a). An alternative to directly estimating the magnetic field strength is to statistically study the magnetic field orientation and its correlation with the observed column density, N_H, distribution. Recent studies of numerical simulations of magnetohydrodynamic (MHD) turbulence indicate that there is a relation between the statistical trends of relative orientation between N_H and B_⊥ and the initial magnetization of the MCs (Soler et al. 2013; Chen et al. 2016; Soler & Hennebelle 2017). The observational study of this trend toward ten nearby MCs reveals based on the Planck observations that on spatial scales ranging from tens of parsecs to 0.2 pc, turbulence is trans- or sub-Alfvénic, that is, the magnetic fields are at least in equipartion with turbulence (Planck Collaboration Int. XXXV 2016, from here on PlanckXXXV).

In this work, we extend the relative orientation study presented in PlanckXXXV to the N_H structures observed by the Herschel satellite at higher resolution than that obtained by Planck (Pilbratt et al. 2010). Because of the difference in angular resolution between the linear polarization observations that are used to infer B_⊥ (10′ FWHM at 353 GHz) and the submillimeter observations that are used to determine N_H (36.′′0 FWHM), the scientific question we address is what is the relation between the N_H structures and the larger-scale magnetic field. This approach is not a substitute for the observation of polarized emission by dust and magnetic fields at higher angular resolutions, but rather uses the existing observations to explore the coupling between matter and magnetic fields across scales. This is central for our understanding of the formation of density structures in a magnetized medium

We here determine whether the relative orientation trends reported in PlanckXXXV are also present when we compare the cloud-scale B_⊥ observed by Planck and the N_H structures observed by Herschel. We also investigate whether the distribution of N_H is related to the relative orientation between N_H and B_⊥. Finally, we study whether the relative orientation between the Planck B_⊥ and the Herschel N_H structures is related to the star formation rates in the studied regions. These three questions can only be addressed with the higher resolution Herschel observations, which allow us to separate different regions within the MCs that were studied in PlanckXXXV and estimate their internal N_H distribution.

The first question addresses the anchoring of density structures by Galactic magnetic fields, which has been suggested as one of the observational proofs of the importance of the magnetic fields in MC formation (Li et al. 2014). Recent studies of MHD simulations indicate that a strong magnetization preserves not only the field direction, but also the orientation of density structures across scales by creating a strong anisotropy in the flows that form and structure MCs (see, e.g., Hull et al. 2017; Gómez et al. 2018; Mocz & Burkhart 2018). Recent observational studies of the relative orientation between the structures that were traced by nine molecular rotational emission lines and B_⊥ toward the Vela C region (d ≈700 pc) indicate a connection between the structure of dense gas on small scales and the cloud-scale magnetic field (Fissel et al. 2019). We here therefore aim to evaluate whether the orientation of the cloud-scale magnetic field is related to the orientation of the density structures at smaller scales in nearby MCs.

The second question focuses on the probability distribution function (PDF) of N_H, which is a frequently used tool for describing the structure of MCs (see, e.g., Lombardi et al. 2008; Kainulainen et al. 2009; Goodman et al. 2009; Schneider et al. 2011). Theoretical studies indicate that the N_H PDFs of MCs are characterized by a log-normal peak and a power-law tail toward high N_H. This is interpreted as the effect of turbulent motions and gravity, respectively (see, e.g., Vazquez-Semadeni 1994; Ballesteros-Paredes et al. 2011; Burkhart 2018). Recent studies based on MHD simulations indicate that magnetic fields can affect N_H PDFs. Specifically, MCs with subcritical mass-to-flux ratios show significantly steeper power tails than supercritical MCs (Auddy et al. 2018). A proof of concept of the study of the relative orientation of N_H and B_⊥ and its relation with the N_H PDFs in observations can be found in Soler et al. (2017), where the authors considered the Serpens Main 2 region in the Aquila Rift.

Finally, the third question concentrates on the relation between the relative orientation of N_H and B_⊥ and the star formation rates (SFRs). This relation has been proclaimed an observational test of the role of magnetic fields as a primary regulator of star formation (Li et al. 2017). We use the reported SFR values evaluated from counts of young stellar objects (YSOs; Lada et al. 2010; Evans et al. 2014) and directly compare them to the N_H and B_⊥ relative orientation trends toward the studied regions.

Throughout this paper we assume that B_⊥ is well represented by the emission from the magnetically aligned dust grains (for a review, see Andersson et al. 2015). This assumption is justified both by most recent observational evidence (Planck Collaboration Int. XX 2015) and by synthetic observations (Seifried et al. 2019)that indicate that the dust grains remain well aligned even at the highest column densities that are relevant for this study (n > 10³ cm⁻³).

This paper is organized as follows. We introduce the Planck and Herschel observations in Sect. 2. Section 3 describes how we implemented the method called the histogram of relative orientations (HRO), which we used for the systematic analysis of the relative orientation between N_H and B_⊥. We present the results of HRO analysis and study the relation between the relative orientations of N_H -B_⊥ and the shape of the N_H PDFs and SFRs in Sect. 4. We discuss the physical implications of our main results in Sect. 5. Finally, Sect. 6 presents our conclusions and future prospects. We reserve the technical details and additional information for a set of appendices. Appendix A presents complementary information on the HRO analysis. Appendix B presents the HRO analysis of the Herschel 250 μm maps toward the Musca and the IC 5146 regions. Appendix C presents some additional comparisons of the relative orientation of N_H and B_⊥ and its relation with the N_H PDFs. All the routines used in this paper, as well as other tools for the analysis of observations and simulations, are publicly available¹.

2 Observations

We present the study of ten nearby (d < 450 pc) MCs that have previously been studied in PlanckXXXV. These regions were partly covered in observations by ESA’s Herschel satellite. We study the subregions using the data products described below. For the sake of consistency, we use the same distances to the regions indicated in Table 1 of PlanckXXXV and use Galactic coordinates in the analysis. These selections do not affect our results.

2.1 Polarization

Planck observed linearly polarized emission (Stokes Q and U) in seven frequency bands from 30 to 353 GHz over the whole sky. In this study, we used the publicly available PR3 data from the High Frequency Instrument (HFI, Lamarre et al. 2010) at 353 GHz. Toward MCs, the contribution of the polarized emission from the cosmic microwave background (CMB) is negligible in this frequency band. This Planck map is therefore best suited for studying the polarized emission from interstellar dust grains (Planck Collaboration Int. XIX 2015).

The maps of Stokes Q, U, their respective variances σ_QQ, σ_UU, and their covariance σ_QU are initially at 4.′8 resolution in HEALPIX format with a pixelization at N_side =2048, which corresponds to an effective pixel size of 1.′72. For the sake of comparison with PlanckXXXV and to increase the signal-to-noise ratio (S/N), we use the maps smoothed to a resolution of 10′ FWHM. The maps ofB_⊥ inferred from these polarization observations are presented in the top panels of Figs. 1–4. Detailed maps of the Herschel-observed regions are presented in Appendix A.

Fig. 1 Top: plane-of-the-sky magnetic field (B⊥) and column density (NH) measured by Planck toward the ten nearby molecular clouds studied in PlanckXXXV. The colors represent the column density inferred from the dust opacity at 353 GHz. The “drapery” pattern, produced using the LIC, indicates the orientation of magnetic field lines, which is orthogonal to the orientation of the submillimeter polarization. The rectangular frames indicate the regions where Herschel observations are available. Detailed maps of these Herschel-observed regions are presented in Appendix A. Second row: relative orientation parameter (ξ), a measure of the trends in the HRO introduced in PlanckXXXV. Values of ξ > 0 correspond to histograms with many counts around 0°, i.e., NH contours mostly parallel to B⊥. Values of ξ < 0 correspond to histograms with many counts around 90°, i.e.,NH contours mostly perpendicular to B⊥. A selection of the HROs is presented in Appendix A. Third row: projected Rayleigh statistic (V), a statistical test of the nonuniformity and unimodal clustering in the distribution of relative orientation angles (ϕk), defined in Eq. (3). Values of V > 0 correspond to ϕk clustered around 0°, i.e., NH contours mostly parallel to B⊥. Values of V < 0 correspond to ϕk clustered around 90°, i.e., NH contours mostly perpendicular to B⊥. Fourth row: mean relative orientation angle ($<\varphi >$) defined in Eq. (5).

2.2 Column density

We used thepublicly available column density maps presented in Abreu-Vicente et al. (2017), which are based on the combinationof the Herschel data and Planck observations. The Herschel data correspond to the publicly available observations obtained with the Spectral and Photometric Imaging Receiver (SPIRE) in three bands centered on 250, 350, and 500 μm (17.′′6, 23.′′ 9, and 35.′′ 2 FWHM resolution, respectively) and with the Photodetector Array Camera and Spectrometer (PACS) in a band centered on 160 μm (13.′′6 FWHM resolution). The Planck data in the corresponding bands were derived from the Planck all-sky foreground dust emission model (Planck Collaboration XI 2014). The Herschel data and Planck observations in each band were cross-calibrated and combined in Fourier space using the feathering method (Csengeri et al. 2016). This method produces maps that improve the previous constant-offset corrections to the Herschel maps by accounting for variations in the background emission levels. The feathered data are convolved to the common 36.′′ 0 F4WHM resolution and are fit using a modified blackbody spectrum to infer the optical depth that is used as a proxy for N_H. Further details are described in Abreu-Vicente et al. (2017).

The Herschel observations toward the MCs studied in PlanckXXXV cover objects in the Lynds Catalog of Dark Nebulae (LDN, Lynds 1962) or simply split the MCs into smaller arbitrary portions that facilitated observing them, for example, North and South in Perseus. We separately analyzed each of the Herschel-observed segments and combined them only in the regions where these segments correspond to a single well-known region, such as CrA, Perseus, Orion A, and Orion B.

3 Histogram of relative orientations

We quantified the relative orientation between the column density structures and B_⊥ using the histogram of relative orientations (HRO, Soler et al. 2013) technique. In the HRO, the N_H structures are characterized by their gradients, ∇N_H, which are by definition perpendicular to the iso-N_H curves. The gradient constitutes a vector field that we compared pixel by pixel to the B_⊥ orientation inferred from the Planck 353-GHz polarization maps.

The angle ϕ between B_⊥ and the tangent to the N_H contours was evaluated using $$\varphi =\arctan \left(|\nabla N_{\text{H}}\times {\widehat{e}}_{353}|,\nabla N_{\text{H}}\cdot {\widehat{e}}_{353}\right),$$(1)

where the pseudo-vector ê₃₅₃ corresponds to the orientation of the Planck 353 GHz linear polarization defined by $$\psi =\frac{1}{2}\arctan \left(-U,Q\right).$$(2)

In Eq. (1) as implemented, the norm carries a sign when the range used for ϕ is between − 90° and 90°. We computed the gradients using the Gaussian derivatives introduced in PlanckXXXV with a kernel size equal to five pixels, thus establishing that the orientation of the N_H structures is characterized on a scale of approximately 1.′17, given the pixel size of the N_H maps Δl = 14′′.

In this implementation, we report the results of the HRO analysis in terms of the histogram-dependent relative orientation parameter, ξ, introduced in Eq. (4) of PlanckXXXV. We also bypass the histogram step, however, by directly applying the tools of circular statistics introduced in Jow et al. (2018), that is, the projected Rayleigh statistic, V, and the mean relative orientation angle, ⟨ϕ⟩ (see Batschelet 1981, for general references in circular statistics). For the sake of completeness, we present the HROs and the maps of the individual subregions in Appendix A.

As in PlanckXXXV, the HRO analysis was performed in N_H bins with equal numbers of pixels, which guarantees comparable statistics. The number of N_H bins was determined by requiring enough bins to resolve the highest N_H regions and at the same time maintaining enough pixels per N_H bin to obtain significant statistics. The measurement errors in the values of ξ, V, and ⟨ϕ⟩ were estimated using Monte Carlo realizations of the noise, quantified by the covariances of the intensity maps in the Herschel bands and the Planck Stokes parameters at 353 GHz, following the procedure described in Appendix B of PlanckXXXV. As in PlanckXXXV, we minimized the effect of noise by smoothing the Planck maps to a 10′ FWHM resolution and applying a 1.′17 FWHM derivative kernel to estimate the gradients. Thus, the reported error bars around the values of ξ, V, and ⟨ϕ⟩ correspond to the standard deviations, which are much larger than the estimated effect of the measurement noise.

Using ξ, V, and ⟨ϕ⟩ as metrics to characterize the distribution of relative orientation angles ϕ_k allows for a more complete interpretation of the HRO results than in PlanckXXXV. The values of ξ provide a description of the shape of the HROs. The values of V constitute an optimal statistical test independent of the HRO binning, and ⟨ϕ⟩ allows evaluating the prevalence of relative orientations different to 0 or 90°, as we detail below.

3.1 Projected Rayleigh statistic V

The projected Rayleigh statistic (V) is a test of the nonuniformity and unimodal clustering of the angle distribution (see, e.g., Chap. 6 in Mardia & Jupp 2009). For a set of N relative orientation vectors ϕ_k, it is defined as $$V\equiv 2\frac{{[{\sum }_k^N\cos (2({\varphi }_k-{\varphi }_0))]}^2}{{\sum }_k^N{w}_k},$$(3)

where w_k corresponds to the statistical weights assigned to each angle ϕ_k. We selected the reference angle ϕ₀ = 0° such that V > 0 corresponds to ϕ_k clustered around 0°, or N_H contours mostly parallel to B_⊥. In the same way, V < 0 corresponds to ϕ_k clustered around 90°, or N_H contours mostly perpendicular to B_⊥.

The projected Rayleigh (V) is primarily a hypothesis test, which means that if |V | is small, there is no evidence that N_H and B_⊥ are oriented toward the reference angles 0 or 90°. If |V | is relatively large, however, there must be some concentration around the reference angles. The larger the value of V, the better chance there is of rejecting the null hypothesis of randomness (Batschelet 1981). A useful analogy for the interpretation of V is the random walk problem. If we set all the weights w_k = 1, Eq. (3) is equivalent to the displacement from the origin in a random walk if steps of unit length were taken in the direction of each angle 2ϕ_k. If the distribution of angles is uniform, then the expectation value of V is 0.

We evaluated the uncertainty that a single V measurement reflects the dispersion of the relative angles following the description presented in Jow et al. (2018). Assuming that the angles ϕ_i are independently and uniformly distributed, it follows that every cos2ϕ_k in Eq. (4) is independently and identically distributed, with a mean of 1∕2. By the central limit theorem, it follows that in the limit of N →∞, V follows a normal distribution with μ = 0 and σ² =1∕2. The asymptotic limit of the V distribution is therefore the standard normal distribution.

For a general distribution of angles, the variance in V in the high N limit is the variance of each $\cos 2{\varphi }_k/\sqrt{1/2}$, which can be estimated as $${\sigma }_V\equiv {\left[\frac{2{\sum }_k^N\cos (2{({\varphi }_k)}^2-V^2}{{\sum }_k^N{w}_k}\right]}^{1/2}.$$(4)

We consider this value to be the uncertainty in the values of V. Jow et al. (2018) used Monte Carlo simulations to show that for uniformly distributed samples of ϕ_k the expectation value of V converges to 0 and σ_V →1. Thus, measurements of V ≫1 indicate a significant detection of parallel relative orientation, while measurements of V ≪− 1 indicate a significant detection of perpendicular relative orientation. Given our selection of a high S/N in the emission and polarization observations, σ_V is dominated by the measurement of a single V in a particular MC rather that by the uncertainties in the observations.

3.2 Mean relative orientation angle ⟨ϕ⟩

If there is a concentration of the ϕ_k values, as can be inferred from the values of V, the mean relative orientation angle is defined as $$<\varphi >\equiv \arctan \left(y,x\right),$$(5)

where $$x\equiv \frac{{\sum }_k^N{w}_k\cos {\varphi }_k}{{\sum }_k^N{w}_k}\text{ and }y\equiv \frac{{\sum }_k^N{w}_k\sin {\varphi }_k}{{\sum }_k^N{w}_k}.$$(6)

Its associated standard deviation is defined as $${\sigma }_{\varphi }\equiv {\left[-2\log (r)\right]}^{1/2},$$(7)

where the mean resultant length, r, which is itself a measure of concentration in the angle distribution, is defined as $$r\equiv {(x^2+y^2)}^{1/2}.$$(8)

If the ϕ_k values are tightly clustered, r is close to 1. If the ϕ_k values are broadly distributed, r is close to 0. A preferential relative orientation is only well defined if r is greater than 0 and σ_ϕ →∞ as r → 0.

3.3 Statistical weights

In this implementation, we accounted for the oversampling of ϕ values introduced by the pixel size of the maps in relation to the size of the beam using the statistical weights $$w_k={\left(\frac{\Delta l}{FWH{M}_{\text{Pla}}}\right)}^2,$$(9)

where Δl = 14′′ is the pixel size of the N_H maps and FWHM_Pla = 10′ corresponds to the Planck beam.

4 Results

4.1 Relative orientation between N_H and B_⊥

The generalresults of the HRO analysis of the ten MCs are presented in the three lower panels of Figs. 1–4. These panels correspond to the mean resultant length, r, the mean relative orientation angle, $<\varphi >$, and the projected Rayleigh statistic, V, introduced in Sect. 3. Because some of the subregions correspond to previously studied objects presented inthe literature, we comment on the HRO results for each of them. This division does not imply that these subregions correspond to a single astronomical object and is made for the sake of convenience in the discussion of our results.

4.1.1 Taurus

The largest subregion covered by the Herschel observations toward Taurus corresponds to the elongated Barnard Dark Nebulae B213 and the dark cloud L1495, which was studied in Palmeirim et al. (2013). Previous studies based on near-infrared polarization observations showed that B_⊥ is mostly perpendicular to B213 (Chapman et al. 2011). The HRO analysis toward this region indicates that the trend in the relative orientation between N_H and B_⊥ corresponds to the archetypical behavior identified in PlanckXXXV for the Taurus region, that is, N_H and B_⊥ are mostly parallel at the lowest N_H and mostly perpendicular at the highest N_H. This is evident in the transition from V > 0 to V < 0 and ⟨ϕ⟩ from 0° to 90° with increasing N_H.

The HRO analysis toward the L1551 region reveals that the smooth increment in ⟨ϕ⟩ with increasing N_H is interrupted around log₁₀(N_H∕cm⁻²) ≈ 21.7 cm⁻². This behavior is similar to that found in the L1489 dark cloud subregion, where the smooth increase of ⟨ϕ⟩ from 0° is interrupted around log₁₀(N_H∕cm⁻²) ≈ 21.5 cm⁻². In principle, the behavior of ⟨ϕ⟩ in L1551 and L1489 can be interpreted as the effect of the projection of the magnetic field and the density structure onto the plane of the sky. The relative orientation between the 2D projections of two vectors that are perpendicular in 3D is likely to be less than 90° for many of the possible configurations of the vectors and the plane of projection. However, the low values of V indicate that the relative orientation angle distributions are not significantly unimodal and centered on 0° or 90°, thus preventing us from drawing further conclusions from ⟨ϕ⟩ alone.

4.1.2 Ophiuchus

The largest subregion in Ophiuchus corresponds to the L1688 dark cloud, an exhaustively studied region of low-mass star formation (Wilking et al. 2008). This region shows a similar behavior to that observed in the B213/L1495 region in Taurus, that is, a clear transition of ⟨ϕ⟩ from 0° to 90° and V > 0 to V < 0 with increasing N_H. This transition from N_H and B_⊥ mostly parallel at the lowest N_H and mostly perpendicular at the log₁₀(N_H∕cm⁻²) > 21.8 is much clearer than that found in PlanckXXXV studies of the whole Ophiuchus region.

The other two subregions in Ophiuchus correspond to the dark clouds that are often called streamers. These trail L1688 to the southeast (L1712) and the east (L1704 and L1709). The HRO analysis toward the L1712 region reveals that ⟨ϕ⟩ is roughly between 20 and 40°, although ξ ≈ 0 and V ≈ 0. Toward L1704/L1709, the HRO analysis reveals that N_H and B_⊥ are mostly parallel, ⟨ϕ⟩ ≈ 0°. In sum, the region with lowest N_H values (L1704/L1709) shows N_H and B_⊥ mostly parallel, the region with intermediate N_H values (L1712) has a transition from mostly parallel to roughly 40°, and finally, the region with the highest N_H values (L1688) shows N_H and B_⊥ changing from mostly parallel to mostly perpendicular with increasing N_H.

4.1.3 Lupus

Lupus I and Lupus III show a clear transition of V > 0 to V < 0 and ⟨ϕ⟩ from 0° to 90° with increasing N_H. In the Lupus IV subregion, which is located at lower Galactic latitude than the other two subregions, V ≈ 0 in for log₁₀(N_H∕cm⁻²) < 21.8, and it is only V < 0 in the highest N_H bin, where ⟨ϕ⟩ is approximately 60°, as shown in Fig. 1.

4.1.4 Chamaeleon-Musca

Musca is the quintessential example of a well-defined elongated molecular cloud with a homogeneous B_⊥ perpendicular to it, as first revealed by starlight polarization observations (Pereyra & Magalhães 2004). The HRO analysis confirms this visual assessment with clearly defined negative values of ξ and V and ⟨ϕ⟩ very close to 90°, as shown in Fig. 2. However, the HOG analysis does not reveal a significant trend of N_H structures being parallel to B_⊥ for log₁₀(N_H∕cm⁻²) > 21.0. This is most likely the result of the smoothing of the Herschel observations to a common resolution to produce the N_H maps, which washes out most of the striations that are clearly visible at 250 μm in its native resolution. Further HRO analysis of the Herschel 250 μm maps and B_⊥, presented in Appendix B, shows a clear transition of V > 0 to V < 0 and ⟨ϕ⟩ from 0° to 90° with increasing N_H. This quantitatively confirms the relative orientation assessment reported in Cox et al. (2016).

The two Chamaeleon clouds covered by the Herschel observations, for which we used the designations in Luhman (2008), show dissimilar behaviors. Chamaeleon I (Cham I) shows a clear transition from V > 0 to V < 0 and ⟨ϕ⟩ from 0° to 90° with increasing N_H, very much like the one found throughout the entire Chamaeleon-Musca region in PlanckXXXV. In contrast, Chamaeleon II (Cham II) presents values of ξ and V close to 0, which indicates that the distribution of relative orientation angles does not clearly peak on either 0° or 90°, as can be inferred from the ⟨ϕ⟩ values shown in Fig. 2.

Cham I and Cham II have similar sizes and masses, but the SFR, based on counting YSOs identified by their infrared excess, is higher in Cham I by roughly a factor of six (Evans et al. 2014). Furthermore, Alves de Oliveira et al. (2014) identified clear differences between the two Cham I and Cham II clouds based on the Herschel observations. On the one hand, Cham I has formed stars, but it seems to have arrived at the end of its star formation. Morphologically, it appears as a ridge that is embedded in low-density gas with striations. On the other hand, Cham II has a clumpy structure with ongoing star formation, and faint striations are only marginally observed.

Toward Cham II the values of V ≈ 0 indicate that the distribution of relative orientation angle is neither unimodal nor centered on 0° or 90°. Still, the values of ⟨ϕ⟩ are close to 0° in the lowest N_H bins, increase up to 60° around log₁₀(N_H∕cm⁻²) ≈ 21.6, and then decrease to roughly 20°. This behavior is one of the two exceptions that we found to the general increment of ⟨ϕ⟩ with increasing N_H.

Fig. 2 Same as Fig. 1 for Chamaeleon-Musca, CrA, and the Aquila Rift.

4.1.5 Corona Australis

The HRO study of the Corona Australis (CrA), presented in Fig. 2, indicates that throughout the entire region, the N_H structures are predominantly parallel to B_⊥ over most N_H values and are only perpendicular in the last N_H bin, which corresponds to log₁₀(N_H∕cm⁻²) > 22. However, the significance of the relative orientation in terms of ξ and V toward CrA is not very high. It is mostly due the correlation introduced by the Planck beam and not by the contrast in the histograms shown in Fig. A.5.

To facilitate the discussion, we divided the CrA region into two portions. CrA North includes the Cr-A, Cr-B, Cr-C, and Cr-D, designations introduced in Nutter et al. (2005). These are the densest regions within CrA and contain most of the gravitationally bound cores (Bresnahan et al. 2018). CrA South includes N_H structures with cometary appearance with fewer starless cores, candidate prestellar cores, and confirmed prestellar cores than CrA North. The HRO analysis reveals significant differences between the two regions. On the one hand, CrA North shows a transition from ξ and V > 0 to ξ and V < 0 with increasing N_H and with ⟨ϕ⟩ around 0° at roughly log₁₀(N_H∕cm⁻²) < 21.7 and then increasing to approximately 70°. Although the statistical significance of the ξ and V < 0 are low, the HROs presented in Fig. A.5 further confirm that they represent a substantial change in the distribution of ϕ toward that portion of the CrA region. On the other hand, CrA South shows ξ and V > 0 and ⟨ϕ⟩ around 0° for all N_H, which only reaches a maximum value log₁₀(N_H∕cm⁻²) = 21.03.

4.1.6 Aquila Rift

The HRO analysis of Serpens Main 2 (Bontemps et al. 2010; Könyves et al. 2015), which has been presented as a proof of concept in (Soler et al. 2017) and is shown here in Fig. 2, reveals a clear transition from N_H and B_⊥ mostly parallel at the lowest N_H and mostly perpendicular at the highest N_H, which is evident in the change of ⟨ϕ⟩ from 0° to 90° and V > 0 to V < 0 with increasing N_H in the rightmost panels of Fig. 2. This confirms the trends found by estimating B_⊥ based on near-infrared polarimetry and its qualitative comparison with the N_H structure in the region (Sugitani et al. 2011).

A similar behavior is found in the Serpens dark cloud (Eiroa et al. 2008), where the significance in the values of V is lower, a smaller cloud has less independentmeasurements of ϕ, but the values of ⟨ϕ⟩ are clearly close to 90° in the highest N_H bin. Toward the Serpens Main 1 and the Serpens West regions, the trends in ξ, V, and ⟨ϕ⟩ reveal that N_H and B_⊥ are mostly parallel, although the values of V ≈ 0 indicate that the significance of these values is low, mostly due the correlation introduced by the Planck beam and not by the contrast in the histograms shown in Fig. A.8.

4.1.7 Perseus

The HRO study of Perseus (Bally et al. 2008), presented in Fig. 3, indicates that over the whole region the relative orientation trends are consistent with those found in PlanckXXXV, that is, from ξ > 0 and V > 0 to ξ < 0 and V < 0 with increasing N_H. To facilitate discussion, we divided Perseus into two portions, following the two Herschel maps toward Perseus. Perseus North includes B5 (L1471), which is one of the most extensively studied dark clouds in the sky, and IC 348. Perseus South includes B1, the dark nebulae L1448 and 1455, and NGC 1333, which is currently the most active region of star formation in Perseus (Sadavoy et al. 2014).

The main difference between Perseus North and Perseus South is the significance in the values of ξ and V, which are roughly a factor for two higher in Perseus North. This difference is not simply the result of the larger number of independent ϕ observations, as shown by the normalized quantity ξ, but it reflects a significant difference in the distribution of its values in each N_H bin, as further illustrated in Fig. A.6. Perseus North and South both show the transition from ξ > 0 and V > 0 to ξ < 0 and V < 0, which corresponds to N_H and B_⊥ being mostly parallel at the lowest N_H and mostly perpendicular at the highest N_H. However, in the highest N_H bin of Perseus South, corresponding to 22.2 < log₁₀(N_H∕cm⁻²) < 23.5, ξ and V return to 0 and ⟨ϕ⟩ is roughly 50°. The values of ξ and V indicate that the distribution of ϕ is neither unimodal nor centered on 0 or 90°. This behavior is the second of the exceptions that we found to the general trends discussed in Sect. 5.

Fig. 3 Same as Fig. 1 for Perseus, IC 5146, and Cepheus.

4.1.8 IC 5146

The IC 5146 region is composed of a bulbous dark cloud at the end of a long dark streamer extending to the northwest, where a dark cloud complex lies (Herbig & Reipurth 2008). The HRO toward that region, presented in Fig. 3, shows clearly defined negative values of ξ and V and ⟨ϕ⟩ very close to 90°, particularly at log₁₀(N_H∕cm⁻²) > 21.6. As in the case of Musca, we repeated the HRO analysis of IC 5146 using the Herschel 250 μm map in Appendix B, which revealed that the lower resolution in the N_H maps washes out structures that are clearly parallel to B_⊥ and would complete the clear transition from ξ > 0 and V > 0 to ξ < 0 and V < 0 with increasing N_H.

4.1.9 Cepheus

The region of the Cepheus flare (Kun et al. 2008) studied in PlanckXXXV corresponds to the northeast portion of the region covered in the Yonekura et al. (1997) ¹³CO survey. Although in PlanckXXXV this region was analyzed without considering individual clouds, the Herschel observations toward Cepheus cover four regions around dark cloud associations: L1241 (Cepheus 2 in Evans et al. 2014), L1172/74 (Cepheus 3), L1157/1152/1155/1147 (Cepheus 4), and L1228 (Cepheus 5). We therefore present the HRO analysis for each of these separate regions in Fig. 3.

The HRO analysis toward L1228 shows a clear transition from V > 0 and ⟨ϕ⟩ ≈ 0° to V < 0 and ⟨ϕ⟩ ≈ 90° with increasing N_H, that is, from N_H and B_⊥ being mostly parallel at the lowest N_H and mostly perpendicular at the highest N_H. The analysis of L1241 shows a similar trend, with less significance and mostly at log₁₀(N_H∕cm⁻²) > 22.0. In contrast with L1228, which is the region with the highest SFR in Cepheus, L1241 shows very little evidence of current star formation (Kirk et al. 2009).

Toward L1172, N_H and B_⊥ appear to be mostly parallel across all N_H values, but the trends show that ξ and V tend to 0 and ⟨ϕ⟩ increasing progressively from 0° to approximately 30°. In contrast, toward L1157, N_H and B_⊥ appear to be perpendicular across all N_H values. The SFRs, based on counting YSOs, indicate that L1157 is currently forming eight times more stars than L1172 (Evans et al. 2014).

4.1.10 Orion

The HRO study of the Orion A, Orion B, and L1622, presented in Fig. 4, indicates that the relative orientation trends are consistent with those found in PlanckXXXV, that is, from ξ > 0 and V > 0 to ξ < 0 and V < 0, or N_H and B_⊥ changing from mostly parallel to mostly perpendicular, with increasing N_H. The high values of ξ and V indicate that these relative orientation trends are very significant, although the values of ⟨ϕ⟩ only reach approximately 70° in the highest N_H bins of the three regions. The similarity in the relative orientation in the three regions is clear in the values of ξ, V, and ⟨ϕ⟩ in the range 21.0 ≲ log₁₀(N_H∕cm⁻²) ≲ 21.8 and it extends to the highest N_H bin, although themaximum N_H values are lower for L1622 and Orion B than for Orion A.

The relative orientation in the highest N_H bins in Orion A, which indicate that N_H and B_⊥ are mostly perpendicular, is consistent with the relative orientation reported at smaller scales toward the eastern end of the region, in an area of roughly 6′ × 6′ around the Kleinmann-Low nebula (KL, Kleinmann & Low 1967) and the Becklin-Neugebauer object (BN, Becklin & Neugebauer 1967), using observations by the SCUBA-2 Polarimeter (POL-2) on the James Clerk Maxwell Telescope (JCMT) at 850 μm (14.′′1 resolution, Pattle et al. 2017b) and the High-Resolution Airborne Wideband Camera-Plus (HAWC+) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA) at 53, 89, 154, and 214 μm (4.′′9, 7.′′ 8, 13.′′ 6, and 18.′′ 2 resolution, respectively, Chuss et al. 2019). This indicates that the hourglass configuration of B_⊥ reported in Pattle et al. (2017b) and Chuss et al. (2019) in arcminute scales is not decoupled from the large-scale B_⊥ configuration, which is perpendicular in degree scales to the integral-shaped filament that hosts the KL-BN region.

The size of the region allows us to further divide the HRO analysis of Orion A into three subregions that we present in the rightmost panels of Fig. 4. The western region contains the KL-BN region and most of the YSOs in Orion A (Großschedl et al. 2019), the center region contains most of the L1641 dark cloud (Polychroni et al. 2013), and finally, the eastern region contains the L1647 dark cloud. Although the western and eastern regions have known differences in terms of the slope of their N_H distribution and number of Class 0 protostars (Stutz & Kainulainen 2015), their behavior in terms of the relative orientation between N_H and B_⊥ is very similar. It changes from mostly parallel to mostly perpendicular with increasing N_H. In contrast, the center region does not show a very clear preferential relative orientation, as inferred from ξ and V ≈ 0, although the ⟨ϕ⟩ values indicate a tendency for the relative orientation to increase from roughly 0° to 50° with increasing N_H.

Fig. 4 Same as Fig. 1 for Orion.

4.2 Relative orientation of N_H and B_⊥ and column density PDFs

Following the considerations presented in Ossenkopf-Okada et al. (2016) and Alves et al. (2017), we focused our analysis of the N_H PDF on the area within the last closed contour in each of the subregions. With this selection, the N_H PDF can be reasonably modeled as $${\log }_{10}(\text{PDF})=C-\alpha {\log }_{10}(N_H),$$(10)

as illustrated in Fig. 5. Although there are deviations from this linear behavior in some regions, the values of the slope α provide a general characterization of the N_H that we can compare to the trends in the relative orientation of N_H -B_⊥.

The left-hand side of Fig. 6 shows a comparison of the values of ξ, V, and ⟨ϕ⟩ evaluated within the last closed contours, and values of α in all of the studied MC subregions. At first glimpse, this figure reveals a general lack of correlation between the values α and the metrics of the relative orientation. The MCs with the steepest N_H PDFs (α >3.2) show values of ξ, V, and ⟨ϕ⟩ that are consistent with N_H and B_⊥ being mostly parallel. Most of the MCs show slopes 3.2 >α <1.8, and within this range lie clear examples of clouds that are predominantly perpendicular, such as Musca (α ≈2.9, ⟨ϕ⟩ ≈ 79° ± 20°) and CrA South (α ≈2.9, ⟨ϕ⟩ ≈ 1° ± 26°). The MCs with the shallowest N_H PDFs (1.8 >α) also show values of ξ, V, and ⟨ϕ⟩ that are consistent with N_H and B_⊥ being mostly parallel. However, the MCs with α <3.2 suggest a general trend: the values of ξ and V appear to progressively decrease with increasing α and change from mostly positive to mostly negative, which is reflected in the increase in ⟨ϕ⟩ from 0 to 90° with increasing α.

We further investigated the correlation between the relative orientation of N_H -B_⊥ and the N_H PDFs by applying a selection criterion that accounts for the potential effect of the polarization background in the B_⊥ that we observed toward the MCs. For this purpose, we evaluated the mean B_⊥ orientation angle within and outside the largest closed N_H contour, $<\psi {>}^{(\text{CC})}$ and $<\psi {>}^{(\text{off})}$, respectively. In principle, MCs where B_⊥ is dominated by the polarization background, here defined as the polarization signal outside of the largest closed N_H contour, have relatively low values of $|<\psi {>}^{(\text{CC})}-<\psi {>}^{(\text{off})}|$. This excludes these MCs and minimizes the effect of the background. A selection of MCs based on $|<\psi {>}^{(\text{CC})}-<\psi {>}^{(\text{off})}|$ automatically excludes MCs where B_⊥ within and outside the largest closed N_H contour are possibly parallel, such as Musca, but it is a sensible initial standard that does not require any further information or modeling of the spatial distribution of dust along the line of sight. We note that this criterion is more effective and simple than a similar approach based on the evaluation of the polarization fraction (p =P∕I), which is more sensitive to the degeneracy between the orientation of the mean magnetic field with respect to the LOS and the amplitude of the turbulent component of the field.

We selected MCs where $|<\psi {>}^{(\text{CC})}-<\psi {>}^{(\text{off})}|$ > 20° and report their values of ξ, V, and ⟨ϕ⟩ in comparison with the values of α on the right-hand side of Fig. 6. With the exception of the L1241 in Cepheus, the selected MCs show a clear trend of decreasing ξ and V values with increasing N_H PDF slopes, that is, the MCs with the steepest N_H PDF tend to be those with N_H and B_⊥ close to perpendicular and those with the shallowest tend to be those with N_H and B_⊥ close to parallel. This tendency is well represented by the L1489 (α ≈1.5, ⟨ϕ⟩ ≈ 28° ± 26°) and the B213/L1495 (α ≈3.1, ⟨ϕ⟩ ≈ 79° ± 20°) regions in Taurus, the South (α ≈1.9, ⟨ϕ⟩ ≈ 29° ± 28°) and North (α ≈2.6, ⟨ϕ⟩ ≈ 28° ± 26°) regions in Perseus, and Orion A (α ≈2.0, ⟨ϕ⟩ ≈ 74° ± 25°), L1622 (α ≈2.4, ⟨ϕ⟩ ≈ 58° ± 27°), and Orion B (α ≈2.7, ⟨ϕ⟩ ≈ 44° ± 28°) regions in Orion.

The results in Fig. 6 indicate that each cloud as a whole, if we define the cloud as the area within the last closed contour, is not exclusively parallel or perpendicular to B_⊥ but its found in a range of ⟨ϕ⟩ between 0 to 90°. This constitutes a critical difference with respect to the results of the studies of the cloud elongation, characterized by its autocorrelation, and the mean field direction (Li et al. 2014). Given that the gradient method for characterizing the relative orientations includes the information from the N_H orientation and B_⊥ within the last closed contour, it is less sensitive to that particular selection, which is an arbitrary boundary determined by the edges of the observed area.

Fig. 5 Probability density functions of NH. The black lines correspond to the fits of Eq. (10) to the NH PDFs within the last closed contour of the subregions that are indicated by the vertical lines.

Fig. 6 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle ($<\varphi >$, bottom) as a function of the slope of the NH distribution (α, Eq. (10)) within the last closed contour of the studied regions before (left) and after (right) the selection criterion $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20°, where $<\psi {>}^{(\text{in})}$ and $<\psi {>}^{(\text{out})}$ are the meanB⊥ orientation angles insideand outside the largest closed NH contour, respectively.

4.3 Relative orientation and SFRs

We used thereported values of the SFRs evaluated from YSO counts (Lada et al. 2010; Evans et al. 2014) and directly compared them to the values ofξ, V and ⟨ϕ⟩ in the regions we studied, as presented in Fig. 7. We considered two portions of each region to estimate the relative orientation: (i) the last closed contour, and (ii) the contour defined by the 95th percentile of N_H. The first selection aims to evaluate the N_H -B_⊥ relative orientation in the largest well-delineated object within the limits of the Herschel map. The second selection aims to evaluate the N_H -B_⊥ relative orientation in the highest N_H portion of each Herschel map, which in principle corresponds to the regions that have accumulated matter more effectively and potentially are more prone to star formation.

When the largest closed contours are considered, we found that the region with the highest SFR per unit of mass (Orion A) shows ⟨ϕ⟩ close to 90°, that is, the N_H structures are close to perpendicular to B_⊥. Regions such as Taurus, Ophiuchus, Perseus, and Orion B show relatively similar SFRs, but different relative orientations: N_H and B_⊥ are preferentially parallel (Ophiuchus), are preferentially perpendicular (Perseus), or have other angles (Taurus and Orion B). We also found that some of the regions with the lowest SFRs (Cham II, Musca, and Lupus) show ⟨ϕ⟩ values close to 90°, but another region with low SFR (Lupus IV) presents ⟨ϕ⟩ ≈ 45°. In sum, however, we do not find an evident correlation between the relative orientation between N_H and B_⊥ and the SFRs in the regions we studied.

The lack of correlation between the relative orientation of N_H -B_⊥ and the SFRs is also present in the highest N_H portions of the studied clouds, which we selected using the 95th percentile of N_H. With this selection, most of the regions show ⟨ϕ⟩ values that are consistent with preferentially perpendicular N_H and B_⊥ in clouds with both high and low SFRs. As expected from the results of the HRO analysis, Cham II shows a change in ⟨ϕ⟩ from 90° to 0°, which is an exception to the general trend of ⟨ϕ⟩ changing from 0° to 90° with increasing N_H.

To facilitate comparison with Li et al. (2017), we also compared the ratio of the SFR to the cloud mass with the values of ξ, V, and ⟨ϕ⟩ in the regions we studied, as presented in Fig. 8. Our comparison is based exclusively on the regions reported in Lada et al. (2010), where the masses are all calculated within the A_V > 8 contours. This selection does not imply any loss of generality and is made to maintain a consistent definition of the masses using published values. Figure 8 shows that both for the last closed contours and for the 95th percentile of N_H, there is no evidence of a clear relation between the relative orientation of N_H -B_⊥ and the ratio of the SFR to cloud mass.

5 Discussion

5.1 General trends in the relative orientation of N_H and B_⊥

If we can point to a common denominator in this HRO study, it is that there is a general trend for ξ and V to decrease and ⟨ϕ⟩ increase from 0° with increasing N_H, illustrated in Figs. 1–4. There are two exceptions that are worth mentioning, however: Cham II (leftmost panel of Fig. 2) and Perseus South (leftmost panel of Fig. 3), where the highest N_H bins show values of V that are less negative and ⟨ϕ⟩ lower than in the preceding N_H bins. There are observational reasons, for example, the orientation with respect to the LOS, and physical conditions, for example, the tangling of the magnetic field at higher densities, that can be used to explain both of these exceptions, as we further discuss in Sect. 5. Nevertheless, it is still remarkable that such a heterogeneous set of MCs for the most part show avsimilar behavior in the relative orientation between the 36′′ -scale N_H structures and the 10′-scale B_⊥.

The change in V from positive to negative and ⟨ϕ⟩ from 0 to 90° is particularly clear in some of the objects that have previously been identified as archetypes of the paradigm of less dense and magnetically-aligned filaments feeding dense filaments perpendicular to the magnetic field (André et al. 2014), for example, B213 in Taurus, Musca, and Serpens Main 2 in Aquila Rift. Here, we have identified that this trend is generally followed in the studied MCs using a method that does not require the identification of filaments. The column densities where the transition from mostly parallel to mostly perpendicular is found are close to those where the change of the slopes in the N_H distributions (Kainulainen et al. 2009) and the scaling of the LOS magnetic field (B_||) with increasing N_H, inferred from the Zeeman observations, have been reported previously (Crutcher 2012). In clouds where the mostly parallel portion of the general trend is missing, for example, Musca and IC 5146, we have shown that the increase in angular resolution obtained when the N_H maps are replaced with the higher resolution 250 μm Herschel observations recovers the mostly parallel to mostly perpendicular N_H and B_⊥ trend, as shown in Appendix B.

5.2 Density structure within MCs anchored in the cloud-scale magnetic fields

We first observe about the relative orientation trends between the N_H structures and B_⊥, reported in Figs. 1–4 that the trend in the relative orientation between N_H and B_⊥ found in PlanckXXXV is also found in many of the subregions observed by Herschel. Notable examples are the B213 region in Taurus, L1688 in Ophiuchus, Lupus I and III, Chamaeleon I, Serpens Main 2, Perseus North, Orion A and B, and L1622, all of which show a progressive change in relative orientation with increasing N_H from mostly parallel to mostly perpendicular. However, the interpretation of these results is different from that in PlanckXXXV.

In this case, the relation corresponds to the coupling between different scales, where the orientation of the N_H structures bears a relation toB_⊥ on a physical scale that is an order of magnitude larger. Insofar as this correlation exists, there is a strong argument that the magnetic field plays a significant role in the distribution of the density structures inside MCs. It suggests that the strong magnetization in the diffuse gas, inferred from the Zeeman H I observations (Heiles & Crutcher 2005) and the anisotropy of H I structures with respect to B_⊥ (Clark et al. 2014, 2019), also affects the density structure of these MCs. If the matter and the magnetic fields are strongly coupled, as is most likely the case in the diffuse ISM, it is unlikely that the accumulation of the parcels of gas that eventually become MCs is made perpendicular to the magnetic field because of the effect of Lorentz force (see discussion in Soler & Hennebelle 2017). Numerical simulations of MC formation indeed indicate that the magnetic field delays the formation of dense molecular gas by imposing an anisotropy in the direction of the mass flows (Ntormousi et al. 2017; Girichidis et al. 2018). It is therefore not surprising that we find the densest parts of the MC perpendicular to the magnetic field.

The magnetic fields are not strong enough to fully determine the orientation of the N_H structures within some MCs, however. Inevitably, gravity may dominate the magnetic fields at least in portions of the MCs, as is evident by the presence of prestellar objects, either by accretion of matter from the larger scales or by the loss of magnetic flux by a diffusion process, such as ambipolar difussion or magnetic reconnection. It is tempting to predict that the trends of the relative orientation of N_H -B_⊥ we reported here extend to the structures at smaller scales, but we find at least two examples of regions where the N_H structures and B_⊥ are mostly perpendicular at high N_H but less perpendicular at the highest N_H: Cham II and the southern portion of Perseus. These examples further confirm that the relative orientation is not necessarily inherited from the largest scales and is subject to change depending on the local conditions, as observed in Zhang et al. (2014).

Numerical simulations of MHD turbulence indicate that only a very strong magnetic field can preserve the relative orientation with respect to the large-scale field from cloud scales (tens of parsecs) to disk scales and to an hourglass-shaped field below 1000 au scales (Hull et al. 2017). Just as the parsec-scale magnetic field is not directly inherited from the Galactic magnetic field, as is evident in the Orion-Eridanus superbubble (Soler et al. 2018), it is premature to extrapolate the relative orientation between N_H structures and B_⊥ at one particular scale to a smaller scale.

5.3 Relation of the relative orientation to the distribution of N_H within MCs

After the selection based on the polarization orientation angle, it is clear that the clouds with the shallowest N_H PDFs appear to be those where ⟨ϕ⟩ is closer to 0°, while the clouds with the steepest N_H PDF tails appear to bethose where ⟨ϕ⟩ is closer to 90°, as illustrated in Fig. 6. This trend is exactly opposite to what would be expected for the gravitational collapse of a magnetized MC. If the magnetic field is dynamically important, a shallow N_H PDF tail, which is interpreted as the effect of gravitational collapse, should appear after the accumulation of matter along the field lines, whose final product are density structures that are mostly perpendicular to the field. By the same token, a steep N_H PDF tail, which is associated with regions that are dominated by turbulence rather than by gravitational collapse, should correspond to relatively early stages of the accumulation of matter along the field lines, when density structures tend to be mostly parallel to the field, where the matter flows are restricted by the Lorentz force. However, these expectations are based on a laminar description of a fluid that is eminently turbulent.

Numerical simulations of MHD turbulence indicate that an increasing magnetic field strength results in narrower density PDFs with steeper tails (Collins et al. 2012; Molina et al. 2012). More recent studies of MHD turbulence in self-gravitating clouds indicate that the slope of the N_H PDF tail is significantly steeper in magnetically subcritical clouds than in supercritical clouds (Auddy et al. 2018). These numerical experiments suggest that the magnetic fields act as a density cushion in turbulent gas, shaping the cloud by preventing the gas from reaching very high densities. Thus, it is plausible that the regions where N_H is mostly perpendicular to B_⊥ and the N_H PDFs are the steepest correspond to magnetically subcritical clouds.

The progressive increase in the values of ⟨ϕ⟩ (decrease in ξ and V) with increasing N_H PDF tail slope is singular. This behavior, which is not exclusive of the last-closed-contour selection, as illustrated in Appendix C, suggeststhat the gravitational collapse and the turbulence are not exclusively responsible for the distribution of N_H and that the magnetic field plays an active role in structuring these MCs. Detailed studies of this trend, using MHD simulations of MCs with different magnetizations and turbulence levels at different evolutionary stages, are beyond the scope of this work, but are instrumental for characterizing the broad range of physical conditions that are responsible for the trends in Fig. 6.

5.4 Relation of the relative orientation to star formation

The correlation between the relative orientations of N_H -B_⊥ and the N_H distributions makes it tempting to extrapolate this trend to the scales of star formation. However, the lack of correlation between the relative orientations of N_H and B_⊥ and the SFRs, illustrated in Figs. 7 and 8, is not entirely unexpected. Without an unusually strong magnetic field, the coupling between the matter and the fields at 0.026 pc scales is on the one hand subject to a broad range of physical conditions that affect star formation at smaller scales, for example, protostellar outflows, turbulence induced by the gravitational collapse, or simply a change in the balance of magnetic to gravitational forces. On the other hand, the observed clouds are not all observed in the same period of the evolution and the SFR derived from the YSO counts only offer one snapshot of their evolution. It is possible that a cloud such as Musca or Lupus I has been formed by flows along the magnetic field and has effectively accumulated matter perpendicular to the field but is in an early stage of its evolution, thus showing only a few YSOs. It is also possible that one of these clouds was unable to accumulate enough matter from its surrounding to globally form as many stars as Orion A, for instance, and is in a late stage of its evolution, confined by pressure and gravity, but supported by the magnetic field against gravitational collapse. Even the eastern and western end of Orion A present the same relative orientation trends, but their star formation activity is very different, as illustrated in the rightmost panel of Fig. 4 or directly in Fig. 9.

These results call into question the interpretation of other analyses that are premised on the idea that the density structures are anchored in the large-scale magnetic field, which determines the elongation of the MCs, and sets their star formation activity. Based on this expectation, Li et al. (2017) proposed that magnetic fields are a primary regulator of the SFRs. We found, however, that the densest portions of the regions with highest SFRs are mostly perpendicular to B_⊥, as are the densest portions of some of the regions with the lowest SFRs, as illustrated in Fig. 7. This lack of correlation with the relative orientation of N_H -B_⊥ is also found in the ratios of the SFR to cloud mass presented in Fig. 8.

The main difference between our results and those reported in Li et al. (2017) lies in the definition of the relative orientation. In Li et al. (2017), the authors define a single structre elongation in N_H, estimated using the N_H correlation function, and a single mean B_⊥ by averaging the polarization observations. In the HRO, each cloud is characterized by the N_H contours, which are compared to the local B_⊥ orientation. This accounts for the internal N_H structure of the regions and the fact that they may not correspond to a single object anchored by the large-scale field. In this sense, the HRO provides greater statistical significance and a more complete description of the N_H structure.

Although the results of our analysis contradict the conclusions of Li et al. (2017), they do not discard the importance of the magnetic fields in the star formation process. Magnetic fields can globally reduce the SFRs in MCs by controlling the accumulation of the material that is available for gravitational collapse and by their influence on the stellar feedback that disperses this material back into the ISM (see, e.g., the discussions in Mac Low et al. 2017; Rahner et al. 2019). Further studies of the dynamics of nearby MCs and the magnetic fields in more distant MCs are necessary, however, to find conclusive evidence of these effects.

Fig. 7 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle (⟨ϕ⟩, bottom) within the largest close contours (left) and above the 95th percentile of column density (right) in each region plotted against the SFRs evaluated from the counts of YSOs reported in Lada et al. (2010) and Evans et al. (2014), represented by the squares and hexagons, respectively. The filled and open markers correspond to regions where $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20° and < 20°, respectively.

Fig. 8 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle (⟨ϕ⟩, bottom) within the largest close contours (left) and above the 95th percentile of column density (right) in each region plotted against the ratio of the SFR to cloud mass. The filled and open squares correspond to regions where $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20° and < 20°, respectively.

Fig. 9 Emission toward the Orion A region observed by Herschel at 160 (blue), 250 (green), and 500 μm (red) and B⊥ (drapery pattern)inferred from the Planck 353 GHz polarization observations. The gray disk represents the size of the Planck beam. The red dot represents the size of the Herschel 500 μm beam.

5.5 Projection onto the plane of the sky

The observed quantities N_H and B_⊥ are the result of the integration along the LOS and the projection onto the plane of the sky. The observed relative orientation may therefore not reflect the orientation of the density structures and the 3D magnetic field (B). One of the main consequences of projection is that two vectors that are parallel in 3D are almost always projected parallel in 2D, but two perpendicular vectors in 3D are only projected as perpendicular if they are in a plane that is almost parallel to the plane of the sky, as discussed in detail in Appendix C of PlanckXXXV and Sect. 4.1 of Zhang et al. (2014). This would imply that if the density and B have no preferred orientations in 3D, it is very unlikely to observe an elongated N_H structure that is almost perpendicular to B_⊥. And yet we see clear examples of this configuration in Taurus, Lupus I, Musca, IC 5146, and Orion A.

However, it is clear that the density and B are not entirely isotropic in the vicinity of the Sun, and toward the nearby clouds we may be observing the N_H and B_⊥ configuration in the walls of the cavity of hot ionized gas that extends all around the Sun, known as the Local Bubble (see Alves et al. 2018, and references therein), the expanding supernova remnant known as the Orion-Eridanus superbubble(see Soler et al. 2018, and references therein), and other nearby superbubbles (see, e.g., Bally 2008). These magnetized superbubbles would generate anisotropies in the mean direction of the magnetic field (⟨B⟩) that would favor the observations of B close to parallel to the plane of the sky. This would explain the prevalence of dense N_H structures perpendicular to B_⊥. Demonstrating this hypothesis is beyond the scope of this work, however.

The observations toward regions where ⟨B⟩ is close to parallel to the LOS should show higher dispersions of B_⊥ than regions where ⟨B⟩ is close to parallel to the plane of the sky (see, e.g., Chen et al. 2019). However, given that the dispersion of B_⊥ can also be the result of the tangling of the field caused by turbulence, it is difficult to determine whether the difference in the relative orientation trends between the clouds presented in this study are exclusively the result of the mean field orientation with respect to the LOS (see Planck Collaboration Int. XLIV 2016, for a discussion of this effect). Comparisons of B_⊥ inferred from starlight polarization and from the Planck 353 GHz observations suggest that in regions like Taurus and Musca, B is mostly in the plane of the sky, while in regions like Lupus I and the Pipe Nebula, ⟨B⟩ is either oriented close to the LOS or high levels of field tangling caused by turbulence (Soler et al. 2016). Breaking this degeneracy between turbulence and the orientation of ⟨B⟩ would require additional information, such as that provided by Zeeman observations and Faraday rotation (Tahani et al. 2018). This is the subject of forthcoming works (Tahani et al. 2019, Sullivan et al., in prep.).

6 Conclusions

We presented a study of the relative orientation between the column density structures, derived from the Herschel observations at 36.′′ 0 resolution, and the plane-of-the-sky magnetic field orientation, inferred from the Planck linear polarization observations at 10′ resolution, toward portions of ten nearby molecular clouds. This is an extension of the analysis presented in Planck Collaboration Int. XXXV (2016).

Our study indicates that the relative orientation trends found in PlanckXXXV, that is, that N_H and B_⊥ change progressively from mostly parallel to mostly perpendicular with increasing N_H, globally persist in the N_H structures that are observed by Herschel at higher angular resolution. The transition from mostly parallel to mostly perpendicular is particularly significant toward a set of subregions that includes B213 in Taurus, L1688 in Ophiuchus, Lupus I, Chamaeleon I, Musca, the northern portion of CrA, Serpens andSerpens Main 2 in the Aquila Rift, Perseus, IC 5146, Orion A, Orion B, and L1622. This common behavior is notable given that the clouds we considered represent a relatively heterogeneous set of physical conditions and environments in the vicinity of the Sun. This provides observational evidence of the coupling between the cloud-scale magnetic fields and the N_H structure within the MCs.

We also studied the relation between the relative orientation of N_H -B_⊥ and the distribution of N_H within the MCs, which is characterized by the N_H PDFs within the last closed N_H contours. We found that the regions with the steepest N_H PDF tails are those where N_H and B_⊥ are close to perpendicular. We also found that regions with shallower PDF tails show a broad range of global relative orientations of N_H and B_⊥, with the shallowest being those regions where N_H and B_⊥ are closer to parallel. These results provide observational evidence of the effects of the magnetic fields in the N_H distribution based on MHD simulations (see, e.g., Collins et al. 2012; Molina et al. 2012; Auddy et al. 2018). However, further numerical studies are necessary to determine the physical mechanisms that shape these nearby MCs and properly account for observational effects, such as the mean orientation of the field with respect to the LOS.

Finally,we investigated the relation between the relative orientation of N_H -B_⊥ and the SFRs and the SFR to cloud-mass ratios and found no significant trends. These results directly disagree with the studies that link the MC elongation and its orientation with respect to the mean orientation of B_⊥ with the aforemention quantities (Li et al. 2017). Our results do not exclude the magnetic field as an important agent in setting the SFRs, however, but call for further studies of their role in the dynamics of the gas reservoirs of these and other star-forming regions and its influence on the dispersion of that reservoir through stellar feedback.

We presented a study that uses the most recent estimates of N_H and B_⊥ toward nearby MCs that have become available through observations made with the ESA Planck and Herschel satellites. Still, our conclusions are limited by the difference in angular resolution between the intensity and polarization observations. This gap is expected to be bridged in the future by balloon-borne experiments, such as the Next-Generation Balloon-borne Large Aperture Submillimeter Telescope (BLAST-TNG, Dober et al. 2014), and proposed satellite missions, such as the Space Infrared Telescope for Cosmology and Astronomy (SPICA) satellite (André et al. 2019), which will yield extended polarization maps that will reveal further details of the magnetic fields within the most nearby regions of star formation.

Acknowledgements

J.D.S. acknowledges funding from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648 505. Part of the discussions that lead to this work took part under the program Milky-Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. J.D.S. thanks the various people who helped with their encouragement and conversation; notably, Henrik Beuther, Andrea Bracco, François Boulanger, Laura Fissel, Patrick Hennebelle, Mordecai MacLow, Peter Martin, Marc-Antoine Miville-Deschênes, and Douglas Scott. J.D.S. also thanks the anonymous reviewer for the thorough review and highly appreciates the comments and suggestions that significantly contributed to improving the quality of this work. The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). Herschel SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University Sussex (UK); Caltech, JPL, NHSC, University Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC (UK); and NASA (USA). Herschel PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KUL, CSL, IMEC (Belgium); CEA, OAMP (France); MPIA (Germany); IFSI, OAP/AOT, OAA/CAISMI, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI (Italy), and CICT/MCT (Spain).

Appendix A Histograms of the relative orientation

As its name indicates, the method histograms of relative orientation (HRO) characterizes the distribution of the relative orientation angles ϕ as defined in Eq. (1). For the sake of discussion, we have reported the trends in the HROs in the main text using the derived quantities that describe the shape of the distribution, as is the case of ξ and V, or its average values, as is the case of ⟨ϕ⟩. Here we present some representation HROs that correspond to the lowest, one intermediate, and to the highest column density (N_H) bins with an equal number of ϕ measurements. These N_H bins are derived from the inversion of the N_H distribution and aim to have comparable error bars.

For thesake of completeness, we also present maps of the gas column density (N_H) and the magnetic field projected onto the plane of the sky and integrated along the LOS (B_⊥) represented by pseudo-vectors and by the drapery pattern that is produced using the line integral convolution method (LIC, Cabral & Leedom 1993). Figures A.1–A.3 correspond to the Taurus, Ophiuchus, and Lupus regions, respectively. Figures A.4 and A.5 correspond to the Chamaeleon-Musca and the CrA regions. Figure A.6 corresponds to both the Perseus and IC 5146 regions. Finally, Figs. A.7–A.9 correspond to the Orion, the Aquila Rift, and the Cepheus regions, respectively.

Fig. A.1 Top: column density (NH) and magnetic field (B⊥) toward the Taurus subregions B213/L1495 (left), L1551 (center), and L1489 (right). The colors represent NH derived from the Herschel data and Planck observations. The drapery pattern represents B⊥, as inferred from the Planck 353 GHz observations. The black circle represents the size of the Planck beam. Bottom: histograms of relative orientations (HROs) for three representative NH bins.

Fig. A.2 Same as Fig. A.1 for the Ophiuchus subregions L1712 (left), L1704/L1409 (center), and L1688 (right).

Fig. A.3 Same as Fig. A.1 for the Lupus I (left), Lupus III (center), and Lupus IV (right) subregions.

Fig. A.4 Same as Fig. A.1 for the Chamaeleon I (left), Chamaeleon II (center), and Musca (right) subregions.

Fig. A.5 Same as Fig. A.1 for the subregions CrA North (left) and CrA South (right).

Fig. A.6 Same as Fig. A.1 for the Perseus North (left) and Perseus South (center), and IC 5146 (right) regions.

Fig. A.7 Same as Fig. A.1 for the Orion A (left), Orion B (center), and L1622 (right) regions.

Fig. A.8 Same as Fig. A.1 for the Aquila Rift subregions Serpens West, Serpens Main 1, Serpens Main 2, and Serpens(clockwise starting from the top left).

Fig. A.9 Same as Fig. A.1 for the Cepheus subregions L1175, L1172, L1228, and L1241 (clockwise starting from the topleft).

Appendix B Striationtest

Because the derivation of column density (N_H) from the combination of Herschel observations requires the degradation of the maps to the same common resolution, some of the striations or wispy elongated structures that appear abundantly around the dense filaments in the 250 μm maps are not present in the N_H maps. In order to study the relative orientation between the magnetic field (B_⊥) and the striations and establish a direct comparison with previous studies of these kind of structures (namely, Cox et al. 2016), in this appendix we quantify the relative orientation between the Herschel 250 μm (I₂₅₀) observationsand B_⊥ toward the Musca and IC 5146 regions using the HROs method.

The relative orientation parameter (ξ) and the projected Rayleigh statistic (V), shown in Fig. B.1, indicate that the structures sampled by the 250 μm observations show a clear transition from V > 0 and V < 0 with increasing I₂₅₀ toward both regions. Moreover, ⟨ϕ⟩ extends from less than 30° to roughly 90° with increasing I₂₅₀. The significance of the results toward Musca is greater than those toward IC 5146, as expected from the greater number of independent ϕ estimates toward that region. However, the HROs presented in Fig. B.2 reveal that the transition from N_H and B_⊥ being mostly parallel to mostly perpendicular is clearer in IC 5146.

Fig. B.1 Same as Fig. 1 for the comparison between the observations of Herschel observations at 250 μm and B⊥ inferred from Planck 353 GHz polarization toward the Musca and IC 5146 regions.

Fig. B.2 Top: emission at 250 μm and magnetic field toward the Musca (left) and IC 5146 (right) regions. The colors represent the emission observed in the Herschel 250 μm band. The draperypattern represents the magnetic field orientation averaged along the LOS and projected onto the plane of the sky, B⊥, as inferred from the Planck 353 GHz data. Bottom: histograms of relative orientations for three representative column density bins.

Appendix C Column density slopes and relative orientation

For the sake of completeness, we present the relation between the relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), and mean relative orientation angle ($<\varphi >$, bottom) as a function of the slope of the N_H distribution (α) using a different selection than that discussed in Sect. 5.3. Figure C.1 shows the result of this relation for the regions within the 80th (left) and 95th (right) percentile of N_H.

Fig. C.1 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle ($<\varphi >$, bottom) as a function of the slope of the NH distribution (α, Eq. (10)) within the 80th (left) and 95th (right) percentile of NH. All regions are selected using the criterion $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20°, where $<\psi {>}^{(\text{in})}$ and $<\psi {>}^{(\text{out})}$ are the meanB⊥ orientation angles insideand outside the largest closed NH contour, respectively.

References

All Figures

Fig. 1 Top: plane-of-the-sky magnetic field (B⊥) and column density (NH) measured by Planck toward the ten nearby molecular clouds studied in PlanckXXXV. The colors represent the column density inferred from the dust opacity at 353 GHz. The “drapery” pattern, produced using the LIC, indicates the orientation of magnetic field lines, which is orthogonal to the orientation of the submillimeter polarization. The rectangular frames indicate the regions where Herschel observations are available. Detailed maps of these Herschel-observed regions are presented in Appendix A. Second row: relative orientation parameter (ξ), a measure of the trends in the HRO introduced in PlanckXXXV. Values of ξ > 0 correspond to histograms with many counts around 0°, i.e., NH contours mostly parallel to B⊥. Values of ξ < 0 correspond to histograms with many counts around 90°, i.e.,NH contours mostly perpendicular to B⊥. A selection of the HROs is presented in Appendix A. Third row: projected Rayleigh statistic (V), a statistical test of the nonuniformity and unimodal clustering in the distribution of relative orientation angles (ϕk), defined in Eq. (3). Values of V > 0 correspond to ϕk clustered around 0°, i.e., NH contours mostly parallel to B⊥. Values of V < 0 correspond to ϕk clustered around 90°, i.e., NH contours mostly perpendicular to B⊥. Fourth row: mean relative orientation angle ($<\varphi >$) defined in Eq. (5).

In the text

	Fig. 2 Same as Fig. 1 for Chamaeleon-Musca, CrA, and the Aquila Rift.
In the text

	Fig. 3 Same as Fig. 1 for Perseus, IC 5146, and Cepheus.
In the text

	Fig. 4 Same as Fig. 1 for Orion.
In the text

	Fig. 5 Probability density functions of NH. The black lines correspond to the fits of Eq. (10) to the NH PDFs within the last closed contour of the subregions that are indicated by the vertical lines.
In the text

Fig. 6 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle ($<\varphi >$, bottom) as a function of the slope of the NH distribution (α, Eq. (10)) within the last closed contour of the studied regions before (left) and after (right) the selection criterion $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20°, where $<\psi {>}^{(\text{in})}$ and $<\psi {>}^{(\text{out})}$ are the meanB⊥ orientation angles insideand outside the largest closed NH contour, respectively.

In the text

Fig. 7 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle (⟨ϕ⟩, bottom) within the largest close contours (left) and above the 95th percentile of column density (right) in each region plotted against the SFRs evaluated from the counts of YSOs reported in Lada et al. (2010) and Evans et al. (2014), represented by the squares and hexagons, respectively. The filled and open markers correspond to regions where $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20° and < 20°, respectively.

In the text

Fig. 8 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle (⟨ϕ⟩, bottom) within the largest close contours (left) and above the 95th percentile of column density (right) in each region plotted against the ratio of the SFR to cloud mass. The filled and open squares correspond to regions where $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20° and < 20°, respectively.

In the text

	Fig. 9 Emission toward the Orion A region observed by Herschel at 160 (blue), 250 (green), and 500 μm (red) and B⊥ (drapery pattern)inferred from the Planck 353 GHz polarization observations. The gray disk represents the size of the Planck beam. The red dot represents the size of the Herschel 500 μm beam.
In the text

Fig. A.1 Top: column density (NH) and magnetic field (B⊥) toward the Taurus subregions B213/L1495 (left), L1551 (center), and L1489 (right). The colors represent NH derived from the Herschel data and Planck observations. The drapery pattern represents B⊥, as inferred from the Planck 353 GHz observations. The black circle represents the size of the Planck beam. Bottom: histograms of relative orientations (HROs) for three representative NH bins.

In the text

	Fig. A.2 Same as Fig. A.1 for the Ophiuchus subregions L1712 (left), L1704/L1409 (center), and L1688 (right).
In the text

	Fig. A.3 Same as Fig. A.1 for the Lupus I (left), Lupus III (center), and Lupus IV (right) subregions.
In the text

	Fig. A.4 Same as Fig. A.1 for the Chamaeleon I (left), Chamaeleon II (center), and Musca (right) subregions.
In the text

	Fig. A.5 Same as Fig. A.1 for the subregions CrA North (left) and CrA South (right).
In the text

	Fig. A.6 Same as Fig. A.1 for the Perseus North (left) and Perseus South (center), and IC 5146 (right) regions.
In the text

	Fig. A.7 Same as Fig. A.1 for the Orion A (left), Orion B (center), and L1622 (right) regions.
In the text

	Fig. A.8 Same as Fig. A.1 for the Aquila Rift subregions Serpens West, Serpens Main 1, Serpens Main 2, and Serpens(clockwise starting from the top left).
In the text

	Fig. A.9 Same as Fig. A.1 for the Cepheus subregions L1175, L1172, L1228, and L1241 (clockwise starting from the topleft).
In the text

	Fig. B.1 Same as Fig. 1 for the comparison between the observations of Herschel observations at 250 μm and B⊥ inferred from Planck 353 GHz polarization toward the Musca and IC 5146 regions.
In the text

Fig. B.2 Top: emission at 250 μm and magnetic field toward the Musca (left) and IC 5146 (right) regions. The colors represent the emission observed in the Herschel 250 μm band. The draperypattern represents the magnetic field orientation averaged along the LOS and projected onto the plane of the sky, B⊥, as inferred from the Planck 353 GHz data. Bottom: histograms of relative orientations for three representative column density bins.

In the text

Fig. C.1 Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), andmean relative orientation angle ($<\varphi >$, bottom) as a function of the slope of the NH distribution (α, Eq. (10)) within the 80th (left) and 95th (right) percentile of NH. All regions are selected using the criterion $|<\psi {>}^{(\text{in})}-<\psi {>}^{(\text{out})}|$ > 20°, where $<\psi {>}^{(\text{in})}$ and $<\psi {>}^{(\text{out})}$ are the meanB⊥ orientation angles insideand outside the largest closed NH contour, respectively.

In the text