Using Herschel and Planck observations to delineate the role of magnetic ﬁelds in molecular cloud structure (Corrigendum)

We present a study of the relative orientation between the magnetic field projected onto the plane of sky (B⊥) on scales down to 0.4 pc, inferred from the polarized thermal emission of Galactic dust observed by Planck at 353 GHz, and the distribution of gas column density (NH) structures on scales down to 0.026 pc, derived from the observations by Herschel in submillimeter wavelengths, toward ten nearby (d < 450 pc) molecular clouds. Using the histogram of relative orientation technique in combination with tools from circular statistics, we found that the mean relative orientation between NH and B⊥ toward these regions increases progressively from 0°, where the NH structures lie mostly parallel to B⊥, with increasing NH, in many cases reaching 90°, where the NH structures lie mostly perpendicular to B⊥. We also compared the relative orientation between NH and B⊥ and the distribution of NH, which is characterized by the slope of the tail of the NH probability density functions (PDFs). We found that the slopes of the NH PDF tail are steepest in regions where NH and B⊥ are close to perpendicular. This coupling between the NH distribution and the magnetic field suggests that the magnetic fields play a significant role in structuring the interstellar medium in and around molecular clouds. However, we found no evident correlation between the star formation rates, estimated from the counts of young stellar objects, and the relative orientation between NH and B⊥ in these regions.


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 transor 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 . 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 , 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 3 cm −3 ). 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 1 .

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.

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-tonoise ratio (S/N), we use the maps smoothed to a resolution of 10 FWHM. The maps of B ⊥ 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.

Column density
We used the publicly available column density maps presented in Abreu-Vicente et al. (2017), which are based on the combination of 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 allsky 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 . 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.

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 where the pseudo-vectorê 353 corresponds to the orientation of the Planck 353 GHz linear polarization defined by 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 Planck-XXXV 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 Planck-XXXV. 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.

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 where w k corresponds to the statistical weights assigned to each angle φ k . We selected the reference angle φ 0 = 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 cos 2φ 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 σ 2 = 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φ k / √ 1/2, which can be estimated as We consider this value to be the uncertainty in the values of V. Jow et al. (2018) used Monte Carlo simulations to show A96, page 3 of 26 A&A 629, A96 (2019) 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.

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 where Its associated standard deviation is defined as where the mean resultant length, r, which is itself a measure of concentration in the angle distribution, is defined as 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.

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 where ∆l = 14 is the pixel size of the N H maps and FWHM Pla = 10 corresponds to the Planck beam.

Relative orientation between N H and B ⊥
The general results 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, φ , and the projected Rayleigh statistic, V, introduced in Sect. 3. Because some of the subregions correspond to previously studied objects presented in the 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.

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 nearinfrared 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 10 (N H /cm −2 ) ≈ 21.7 cm −2 . This behavior is similar to that found in the L1489 dark cloud subregion, where the smooth increase of φ from 0 • is interrupted around log 10 (N H /cm −2 ) ≈ 21.5 cm −2 . 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.

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 10 (N H /cm −2 ) > 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 .

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 10 (N H /cm −2 ) < 21.8, and it is only V < 0 in the highest N H bin, where φ is approximately 60 • , as shown in Fig. 1. A96, page 4 of 26 Top: plane-of-the-sky magnetic field (B ⊥ ) and column density (N H ) 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., N H contours mostly parallel to B ⊥ . Values of ξ < 0 correspond to histograms with many counts around 90 • , i.e., N H 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., N H contours mostly parallel to B ⊥ . Values of V < 0 correspond to φ k clustered around 90 • , i.e., N H contours mostly perpendicular to B ⊥ . Fourth row: mean relative orientation angle ( φ ) defined in Eq. (5).

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 10 (N H /cm −2 ) > 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 10 (N H /cm −2 ) ≈ 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 . A96, page 6 of 26

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 10 (N H /cm −2 ) > 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 10 (N H /cm −2 ) < 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 10 (N H /cm −2 ) = 21.03.

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 ) 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 independent measurements 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.

Perseus
The HRO study of Perseus , 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 10 (N H /cm −2 ) < 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.

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 10 (N H /cm −2 ) > 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 .

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) 13 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 10 (N H /cm −2 ) > 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).

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 A96, page 7 of 26 A&A 629, A96 (2019) Fig. 3. Same as Fig. 1 for Perseus, IC 5146, and Cepheus.
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 10 (N H /cm −2 ) 21.8 and it extends to the highest N H bin, although the maximum 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 A96, page 8 of 26 J. D. Soler: Using Herschel and Planck to delineate the role of magnetic fields in molecular cloud structure 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 . φ ≈ 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 α.

Relative orientation of N H and B ⊥ and column density PDFs
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, ψ (CC) and ψ (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 | ψ (CC) − ψ (off) |. This excludes these MCs and minimizes the effect of the background. A selection of MCs based on | ψ (CC) − ψ (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. A96, page 10 of 26 Fig. 6. Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), and mean relative orientation angle ( φ , bottom) as a function of the slope of the N H distribution (α, Eq. (10)) within the last closed contour of the studied regions before (left) and after (right) the selection criterion | ψ (in) − ψ (out) | > 20 • , where ψ (in) and ψ (out) are the mean B ⊥ orientation angles inside and outside the largest closed N H contour, respectively.
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.

Relative orientation and SFRs
We used the reported 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 A96, page 11 of 26 A&A 629, A96 (2019) Fig. 7. Relative orientation parameter (ξ, top), projected Rayleigh statistic (V, middle), and mean 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 | ψ (in) − ψ (out) | > 20 • and < 20 • , respectively. 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.

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 A96, page 12 of 26 J. D. Soler: Using Herschel and Planck to delineate the role of magnetic fields in molecular cloud structure 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 , 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.

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 to B ⊥ 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 HI observations (Heiles & Crutcher 2005) and the anisotropy of HI structures with respect to B ⊥ (Clark et al. 2014(Clark et al. , 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 A96, page 13 of 26 A&A 629, A96 (2019) 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.

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 be those 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, suggests that 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.

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 A96, page 14 of 26 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. of nearby MCs and the magnetic fields in more distant MCs are necessary, however, to find conclusive evidence of these effects.

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 Planck-XXXV 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.).

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 and Serpens 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 A96, page 15 of 26 A&A 629, A96 (2019) 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.