© ESO, 2016

1. Introduction

The formation and evolution of molecular clouds (MCs) and their substructures, from filaments to cores and eventually to stars, is the product of the interaction between turbulence, magnetic fields, and gravity (Bergin & Tafalla 2007; McKee & Ostriker 2007). The study of the relative importance of these dynamical processes is limited by the observational techniques used to evaluate them. These limitations have been particularly critical when integrating magnetic fields into the general picture of MC dynamics (Elmegreen & Scalo 2004; Crutcher 2012; Heiles & Haverkorn 2012; Hennebelle & Falgarone 2012).

There are two primary methods of measuring magnetic fields in the dense interstellar medium (ISM). First, observation of the Zeeman effect in molecular lines provides the line-of-sight component of the field B_∥ (Crutcher 2005). Second, polarization maps – in extinction from background stars and emission from dust – reveal the orientation of the field averaged along the line of sight and projected on the plane of the sky (Hiltner 1949; Davis & Greenstein 1951; Hildebrand 1988; Planck Collaboration Int. XXI 2015).

Analysis of the Zeeman effect observations presented by Crutcher et al. (2010) shows that in the diffuse ISM sampled by Hi lines (n_H< 300 cm^-3), the maximum magnetic field strength B_max does not scale with density. This is interpreted as the effect of diffuse clouds assembled by flows along magnetic field lines, which would increase the density but not the magnetic field strength. In the denser regions (n_H> 300 cm^-3), probed by OH and CN spectral lines, the same study reports a scaling of the maximum magnetic field strength ${\mathit{B}}_{\max }\mathrm{\propto }{\mathit{n}}_{\mathrm{H}}^{\mathrm{0.65}}$. The latter observation can be interpreted as the effect of isotropic contraction of gas too weakly magnetized for the magnetic field to affect the morphology of the collapse. However, given that the observations are restricted to pencil-like lines of sight and the molecular tracers are not homogeneously distributed, the Zeeman effect measurements alone are not sufficient to determine the relative importance of the magnetic field at the multiple scales within MCs.

The observation of starlight polarization provides an estimate of the projected magnetic field orientation in particular lines of sight. Starlight polarization observations show coherent magnetic fields around density structures in MCs (Pereyra & Magalhães 2004; Franco et al. 2010; Sugitani et al. 2011; Chapman et al. 2011; Santos et al. 2014). The coherent polarization morphology can be interpreted as the result of dynamically important magnetic fields. However, these observations alone are not sufficient to map even the projected magnetic field morphology fully and in particular do not tightly constrain the role of magnetic fields in the formation of structure inside MCs.

The study of magnetic field orientation within the MCs is possible through the observation of polarized thermal emission from dust. Far-infrared and submillimetre polarimetric observations have been limited to small regions up to hundreds of square arcminutes within clouds (Li et al. 2006; Matthews et al. 2014) or to large sections of the Galactic plane at a resolution of several degrees (Benoît et al. 2004; Bierman et al. 2011). On the scale of prestellar cores and cloud segments, these observations reveal both significant levels of polarized emission and coherent field morphologies (Ward-Thompson et al. 2000; Dotson et al. 2000; Matthews et al. 2009).

The strength of the magnetic field projected on the plane of the sky (B_⊥) can be estimated from polarization maps using the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951; Chandrasekhar & Fermi 1953). As discussed in Appendix D, it is assumed that the dispersion in polarization angle ς_ψ¹ is entirely due to incompressible and isotropic turbulence. Turbulence also affects the motion of the gas and so broadens profiles of emission and absorption lines, as quantified by dispersion σ_v∥. In the DCF interpretation B_⊥ is proportional to the ratio σ_v∥/ς_ψ. Application of the DCF method to subregions of the Taurus MC gives estimates of B_⊥ ≈ 10 μG in low-density regions and ≈ 25 to ≈ 42 μG inside filamentary structures (Chapman et al. 2011). Values of B_⊥ ≈ 760 μG have been found in dense parts of the Orion MC region (Houde et al. 2009). Because of the experimental difficulties involved in producing large polarization maps, a complete statistical study of the magnetic field variation across multiple scales is not yet available.

Additional information on the effects of the magnetic field on the cloud structure is found by studying the magnetic field orientation inferred from polarization observations relative to the orientation of the column density structures. Patterns of relative orientation have been described qualitatively in simulations of magnetohydrodynamic (MHD) turbulence with different degrees of magnetization. This is quantified as half the ratio of the gas pressure to the mean-field magnetic pressure (Ostriker et al. 2001; Heitsch et al. 2001), with the resulting turbulence ranging from sub-Alfvénic to super-Alfvénic. Quantitative analysis of simulation cubes, where the orientation of B is available directly, reveals a preferred orientation relative to density structures that depends on the initial magnetization of the cloud (Hennebelle 2013; Soler et al. 2013). Using simple models of dust grain alignment and polarization efficiency to produce synthetic observations of the simulations, Soler et al. (2013) showed that the preferred relative orientation and its systematic dependence on the degree of magnetization are preserved.

Observational studies of relative orientation have mostly relied on visual inspection of polarization maps (e.g., Myers & Goodman 1991; Dotson 1996). This is adequate for evaluating general trends in the orientation of the field. However, it is limited ultimately by the need to represent the field orientation with pseudo-vectors, because when a large polarization map is to be overlaid on a scalar-field map, such as intensity or column density, only a selection of pseudo-vectors can be plotted. On the one hand, if the plotted pseudo-vectors are the result of averaging the Stokes parameters over a region, then the combined visualization illustrates different scales in the polarization and in the scalar field. On the other hand, if the plotted pseudo-vectors correspond to the polarization in a particular pixel, then the illustrated pattern is influenced by small-scale fluctuations that might not be significant in evaluating any trend in relative orientation.

Tassis et al. (2009) present a statistical study of relative orientation between structures in the intensity and the inferred magnetic field from polarization measured at 350μm towards 32 Galactic clouds in maps of a few arcminutes in size. Comparing the mean direction of the field to the semi-major axis of each cloud, they find that the field is mostly perpendicular to that axis. Similarly, Li et al. (2013) compared the relative orientation in 13 clouds in the Gould Belt, calculating the main cloud orientation from the extinction map and the mean orientation of the intercloud magnetic field from starlight polarization. That study reported a bimodal distribution of relative cloud and field orientations; that is, some MCs are oriented perpendicular and some parallel to the mean orientation of the intercloud field. In both studies each cloud constitutes one independent observation of relative orientation, so that the statistical significance of each study depends on the total number of clouds observed. In a few regions of smaller scales, roughly a few tenths of a parsec, Koch et al. (2013) report a preferred orientation of the magnetic field, inferred from polarized dust emission, parallel to the gradient of the emission intensity.

By measuring the intensity and polarization of thermal emission from Galactic dust over the whole sky and down to scales that probe the interiors of nearby MCs, Planck² provides an unprecedented data set from a single instrument and with a common calibration scheme, for studying the morphology of the magnetic field in MCs and the surrounding ISM, as illustrated for the Taurus region in Fig. 1. We present a quantitative analysis of the relative orientation in a set of nearby (d< 450 pc) well-known MCs to quantify the role of the magnetic field in the formation of density structures on physical scales ranging from tens of parsecs to approximately one parsec in the nearest clouds.

The present work is an extension of previous findings, as reported by the Planck collaboration, on their study of the polarized thermal emission from Galactic dust. Previous studies include an overview of this emission (Planck Collaboration Int. XIX 2015), which reported dust polarization percentages up to 20% at low N_H, decreasing systematically with increasing N_H to a low plateau for regions with N_H> 10²² cm^-2. Planck Collaboration Int. XX (2015) presented a comparison of the polarized thermal emission from Galactic dust with results from simulations of MHD turbulence, focusing on the statistics of the polarization fractions and angles. Synthetic observations were made of the simulations under the simple assumption of homogeneous dust grain alignment efficiency. Both studies reported that the largest polarization fractions are reached in the most diffuse regions. Additionally, there is an anti-correlation between the polarization percentage and the dispersion of the polarization angle. This anti-correlation is reproduced well by the synthetic observations, indicating that it is essentially caused by the turbulent structure of the magnetic field.

Over most of the sky Planck Collaboration Int. XXXII (2016) analysed the relative orientation between density structures, which is characterized by the Hessian matrix, and polarization, revealing that most of the elongated structures (filaments or ridges) have counterparts in the Stokes Q and U maps. This implies that in these structures, the magnetic field has a well-defined mean direction on the scales probed by Planck. Furthermore, the ridges are predominantly aligned with the magnetic field measured on the structures. This statistical trend becomes more striking for decreasing column density and, as expected from the potential effects of projection, for increasing polarization fraction. There is no alignment for the highest column density ridges in the N_H ≳ 10²² cm^-2 sample. Planck Collaboration Int. XXXIII (2016) studied the polarization properties of three nearby filaments, showing by geometrical modelling that the magnetic field in those representative regions has a well-defined mean direction that is different from the field orientation in the surroundings.

In the present work, we quantitatively evaluate the relative orientation of the magnetic field inferred from the Planck polarization observations with respect to the gas column density structures, using the histogram of relative orientations (HRO, Soler et al. 2013). The HRO is a novel statistical tool that quantifies the relative orientation of each polarization measurement with respect to the column density gradient, making use of the unprecedented statistics provided by the Planck polarization observations. The HRO can also be evaluated in both 3D simulation data cubes and synthetic observations, thereby providing a direct comparison between observations and the physical conditions included in MHD simulations. We compare the results of the HRO applied to the Planck observations with the results of the same analysis applied to synthetic observations of MHD simulations of super-Alfvénic, Alfvénic, and sub-Alfvénic turbulence.

Thus by comparison with numerical simulations of MHD turbulence, the HRO provides estimates of the magnetic field strength without any of the assumptions involved in the DCF method. For comparison, we estimate B_⊥ using the DCF method and the related method described by Hildebrand et al. (2009; DCF+SF, for DCF plus structure function) and provide a critical assessment of their applicability.

This paper is organized as follows. Section 2 introduces the Planck353 GHz polarization observations, the gas column density maps, and the CO line observations used to derive the velocity information. The particular regions where we evaluate the relative orientation between the magnetic field and the column density structures are presented in Sect. 3. Section 4 describes the statistical tools used for the study of these relative orientations. In Sect. 5 we discuss our results and their implications in the general picture of cloud formation. Finally, Sect. 6 summarizes the main results. Additional information on the selection of the polarization data, the estimation of uncertainties affecting the statistical method, and the statistical significance of the relative orientation studies can be found in Appendices A−C, respectively. Appendix D presents alternative estimates of the magnetic field strength in each region.

2. Data

2.1. Thermal dust polarization

Over the whole sky Planck observed the linear polarization (Stokes Q and U) in seven frequency bands from 30 to 353 GHz (Planck Collaboration I 2014). In this study, we used data from the High Frequency Instrument (HFI, Lamarre et al. 2010) at 353 GHz, the highest frequency band that is sensitive to polarization. Towards MCs the contribution of the cosmic microwave background (CMB) polarized emission is negligible at 353 GHz, making this the Planck map that is best suited to studying the spatial structure of the dust polarization (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015).

We used the Stokes Q and U maps and the associated noise maps made from five independent consecutive sky surveys of the Planck cryogenic mission, which together correspond to the DR3 (delta-DX11d) internal data release. We refer to previous Planck publications for the data processing, map making, photometric calibration, and photometric uncertainties (Planck Collaboration II 2014; Planck Collaboration V 2014; Planck Collaboration VI 2014; Planck Collaboration VIII 2014). As in the first Planck polarization papers, we used the International Astronomical Union (IAU) conventions for the polarization angle, measured from the local direction to the north Galactic pole with positive values increasing towards the east.

The maps of Q, U, their respective variances ${\mathit{\sigma }}_{\mathit{Q}}^{\mathrm{2}}$, ${\mathit{\sigma }}_{\mathit{U}}^{\mathrm{2}}$, and their covariance σ_QU are initially at resolution in HEALPix format³ with a pixelization at N_side = 2048, which corresponds to an effective pixel size of . To increase the signal-to-noise ratio (S/N) of extended emission, we smoothed all the maps to 10′ resolution using a Gaussian approximation to the Planck beam and the covariance smoothing procedures described in Planck Collaboration Int. XIX (2015).

The maps of the individual regions are projected and resampled onto a Cartesian grid by using the gnomonic projection procedure described in Paradis et al. (2012). The HRO analysis is performed on these projected maps.

2.2. Column density

We used the dust optical depth at 353 GHz (τ₃₅₃) as a proxy for the gas column density (N_H). The τ₃₅₃ map (Planck Collaboration XI 2014) was derived from the all-sky Planck intensity observations at 353, 545, and 857 GHz, and the IRAS observations at 100μm, which were fitted using a modified black body spectrum. Other parameters obtained from this fit are the temperature and the spectral index of the dust opacity. The τ₃₅₃ map, computed initially at 5′ resolution, was smoothed to 10′ to match the polarization maps. The errors resulting from smoothing the product τ₃₅₃ map, rather than the underlying data, are negligible compared to the uncertainties in the dust opacity and do not significantly affect the results of this study.

To scale from τ₃₅₃ to N_H, following Planck Collaboration XI (2014), we adopted the dust opacity found using Galactic extinction measurements of quasars, ${\mathit{\tau }}_{\mathrm{353}}\mathit{/}{\mathit{N}}_{\mathrm{H}}\mathrm{=}\mathrm{1.2}\mathrm{\times }{\mathrm{10}}^{-26}\mathrm{cm}\mathrm{2.}$(1)

Variations in dust opacity are present even in the diffuse ISM and the opacity increases systematically by a factor of 2 from the diffuse to the denser ISM (Planck Collaboration XXIV 2011; Martin et al. 2012; Planck Collaboration XI 2014), but our results do not critically depend on this calibration.

3. Analysed regions

The selected regions, shown in Fig. 2, correspond to nearby (d< 450 pc) MCs, whose characteristics are well studied and can be used for cloud-to-cloud comparison (Poppel 1997; Reipurth 2008). Their properties are summarized in Table 1, which includes: Galactic longitude l and latitude b at the centre of the field; field size Δl × Δb; estimate of distance; mean and maximum total column densities from dust, ⟨N_H⟩ and $\max {}^{\mathrm{(}}\mathit{N}_{\mathrm{H}}{}^{\mathrm{)}}$, respectively; and mean H₂ column density from CO.

In the table the regions are organized from the nearest to the farthest in three groups: (a) regions located at d ≈ 150 pc, namely Taurus, Ophiuchus, Lupus, Chamaeleon-Musca, and Corona Australis (CrA); (b) regions located at d ≈ 300 pc, Aquila Rift and Perseus; and (c) regions located at d ≈ 450 pc, IC 5146, Cepheus, and Orion.

Among the clouds in the first group (all shown in Fig. 3, left column) are Ophiuchus and Lupus, which are two regions with different star-forming activities but are close neighbours within an environment disturbed by the Sco-Cen OB association (Wilking et al. 2008; Comerón 2008). Chamaeleon-Musca is a region evolving in isolation,and it is relatively unperturbed (Luhman 2008). Taurus (see also Fig. 1), a cloud with low-mass star formation, appears to be formed by the material swept up by an ancient superbubble centred on the Cas-Tau group (Kenyon et al. 2008). Finally, CrA is one of the nearest regions with recent intermediate- and low-mass star formation, possibly formed by a high-velocity cloud impact on the Galactic plane (Neuhäuser & Forbrich 2008).

In the second group we consider Aquila Rift and Perseus, shown in Fig. 4 (left column). Aquila Rift is a large complex of dark clouds where star formation proceeds in isolated pockets (Eiroa et al. 2008; Prato et al. 2008). The Perseus MC is the most active site of on-going star formation within 300 pc of the Sun. It features a large velocity gradient and is located close to hot stars that might have impacted its structure (Bally et al. 2008).

In the third group are IC 5146, Cepheus, and Orion, shown in Fig. 5 (left column). IC 5146 is an MC complex in Cygnus. It includes an open cluster surrounded by a bright optical nebulosity called the Cocoon nebula, and a region of embedded lower-mass star formation known as the IC 5146 Northern⁴ Streamer (Harvey et al. 2008). The Cepheus Flare, called simply Cepheus in this study, is a large complex of dark clouds that seems to belong to an even larger expanding shell from an old supernova remnant (Kun et al. 2008). Orion is a dark cloud complex with on-going high and low mass star formation, whose structure appears to be affected by multiple nearby hot stars (Bally 2008).

When taking background/foreground emission and noise within these regions into account, pixels are selected for analysis according to criteria for the gradient of the column density (Appendix A.2) and the polarization (Appendix A.2).

4. Statistical study of the relative orientation of the magnetic field and column density structure

4.1. Methodology

4.1.1. Histogram of relative orientations

We quantify the relative orientation of the magnetic field with respect to the column density structures using the HRO (Soler et al. 2013). The column density structures are characterized by their gradients, which are by definition perpendicular to the iso-column density curves (see calculation in Appendix B.1). The gradient constitutes a vector field that we compare pixel by pixel to the magnetic field orientation inferred from the polarization maps.

In practice we use τ₃₅₃ as a proxy for N_H (Sect. 2.2). The angle φ between B_⊥ and the tangent to the τ₃₅₃ contours is evaluated using⁵$\mathit{\phi }\mathrm{=}\arctan \left(\hspace{0.17em}\mathrm{|}\mathrm{\nabla }\hspace{0.17em}{\mathit{\tau }}_{\mathrm{353}}\mathrm{\times }{Ê}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{,}\hspace{0.17em}\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}\mathrm{·}{Ê}\hspace{0.17em}\right)\mathit{,}$(2)where, as illustrated in Fig. 6, ∇ τ₃₅₃ is perpendicular to the tangent of the iso-τ₃₅₃ curves, the orientation of the unit polarization pseudo-vector Ê, perpendicular to B_⊥, is characterized by the polarization angle $\mathit{\psi }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \mathrm{(}\mathrm{-}\mathit{U,Q}\mathrm{)}\mathit{,}$(3)and in Eq. (2), as implemented, the norm actually carries a sign when the range used for φ is between − 90° and 90°.

The uncertainties in φ due to the variance of the τ₃₅₃ map and the noise properties of Stokes Q and U at each pixel are characterized in Appendix B.

The gradient technique is one of multiple methods for characterizing the orientation of structures in a scalar field. Other methods, which include the Hessian matrix analysis (Molinari et al. 2011; Planck Collaboration Int. XXXII 2016) and the inertia matrix (Hennebelle 2013), are appropriate for measuring the orientation of ridges, i.e., the central regions of filamentary structures. The gradient technique is sensitive to contours and in that sense it is better suited to characterizing changes in the relative orientation in extended regions, not just on the crests of structures (Soler et al. 2013; Planck Collaboration Int. XXXII 2016). Additionally, the gradient technique can sample multiple scales by increasing the size of the vicinity of pixels used for its calculation (derivative kernel; see Appendix B.1). Previous studies that assign an average orientation of the cloud (Tassis et al. 2009; Li et al. 2013) are equivalent to studying the relative orientation using a derivative kernel close to the cloud size.

The selected pixels belong to the regions of each map where the magnitude of the gradient | ∇τ₃₅₃ | is greater than in a diffuse reference field (Appendix A). This selection criterion aims at separating the structure of the cloud from the structure of the background using the reference field as a proxy. For each region the selected reference field is the region with the same size and Galactic latitude that has the lowest average N_H (see Appendix A.1).

In addition to the selection on | ∇τ₃₅₃ |, we only consider pixels where the norms of the Stokes Q and U are larger than in the diffuse reference field, therefore minimizing the effect of background/foreground polarization external to the cloud. The relative orientation angle, φ, is computed by using polarization measurements with a high S/N in Stokes Q and U, i.e., only considering pixels with | Q | /σ_Q or | U | /σ_U> 3. This selection allows the unambiguous definition of Ê by constraining the uncertainty in the polarization angle (see Appendix A.2).

Once we have produced a map of relative orientations for selected pixels following Eq. (2), we divide the map into bins of N_H containing an equal number of pixels and generate a histogram of φ for each bin. The shape of the histogram is used to evaluate the preferred relative orientation in each bin directly. A concave histogram, peaking at 0°, corresponds to the preferred alignment of B_⊥ with the N_H contours. A convex histogram, peaking at 90° and/or − 90°, corresponds to the preferred orientation of B_⊥ perpendicular to the N_H contours.

The HROs in each region are computed in 25 N_H bins having equal numbers of selected pixels (10 bins in two regions with fewer pixels, CrA and IC 5146). The number of N_H bins is 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 from each histogram. The typical number of pixels per bin of N_H ranges from approximately 600 in CrA to around 4000 in Chamaeleon-Musca. We use 12 angle bins of width 15°.

The HROs of the first group of regions, the nearest at d ≈ 150 pc, are shown in the right-hand column of Fig. 3. For the sake of clarity, we only present the histograms that correspond to three bins, namely the lowest and highest N_H and an intermediate N_Hvalue. The intermediate bin is the 12th (sixth in two regions with fewer pixels, CrA and IC 5146), and it corresponds to pixels near the blue contour in the image in the left-hand column of Fig. 3. The widths of the shaded areas for each histogram correspond to the 1σ uncertainties related to the histogram binning operation, which are greater than the uncertainties produced by the variances of Q, U, and τ₃₅₃ (Appendix B). The sharp and narrow features (“jitter”) in the HROs are independent of these variances. They are the product of sampling the spatial correlations in the magnetic field over a finite region of the sky together with the histogram binning; these features average out when evaluating the relative orientation over larger portions of the sky (Planck Collaboration Int. XXXII 2016).

Although often asymmetric, most histograms reveal a change in the preferred relative orientation across N_H bins. The most significant feature in the HROs of Taurus, Ophiuchus, and Chamaeleon-Musca is the drastic change in relative orientation from parallel in the lowest N_H bin to perpendicular in the highest N_H bin. In Lupus the behaviour at low N_H is not clear, but at high N_H it is clearly perpendicular. In contrast, CrA tends to show B_⊥ as parallel in the intermediate N_H bin, but no preferred orientation in the other N_H bins.

The HROs of the clouds located at d ≈ 300 pc, Aquila Rift and Perseus, are shown in the right-hand column of Fig. 4. They indicate that the relative orientation is usually parallel in the lowest N_H bins and perpendicular in the highest N_H bins.

The HROs of the third group, located at d ≈ 400−450 pc, IC 5146, Cepheus, and Orion, are presented in the right-hand column of Fig. 5. In both IC 5146 and Orion the HROs for the highest N_H bins reveal a preferred orientation of the field perpendicular to the N_H contours (Orion is quite asymmetric), whereas the HROs corresponding to the low and intermediate N_H bins reveal a preferred alignment of the field with N_H structures. This trend is also present, but less pronounced, in the Cepheus region.

4.1.2. Histogram shape parameter ξ

The changes in the HROs are quantified using the histogram shape parameter ξ, defined as $\mathit{\xi }\mathrm{=}\frac{{\mathit{A}}_{\mathrm{c}}\mathrm{-}{\mathit{A}}_{\mathrm{e}}}{{\mathit{A}}_{\mathrm{c}}\mathrm{+}{\mathit{A}}_{\mathrm{e}}}\mathit{,}$(4)where A_c is the area in the centre of the histogram () and A_e the area in the extremes of the histogram ( and ). The value of ξ, the result of the integration of the histogram over 45° ranges, is independent of the number of bins selected to represent the histogram if the bin widths are smaller than the integration range.

A concave histogram corresponding to B_⊥ mostly aligned with N_H contours would have ξ> 0. A convex histogram corresponding to B_⊥ mostly perpendicular to N_H contours would have ξ< 0. A flat histogram corresponding to no preferred relative orientation would have ξ ≈ 0.

The uncertainty in ξ, σ_ξ, is obtained from ${\mathit{\sigma }}_{\mathit{\xi }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}\hspace{0.17em}\mathrm{(}{\mathit{A}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{A}}_{\mathrm{c}}}^{\mathrm{2}}\mathrm{+}{\mathit{A}}_{\mathrm{c}}^{\mathrm{2}}{{\mathit{\sigma }}_{{\mathit{A}}_{\mathrm{e}}}^{\mathrm{2}}}^{\mathrm{)}}}{{\mathrm{(}{\mathit{A}}_{\mathrm{c}}\mathrm{+}{\mathit{A}}_{\mathrm{e}}\mathrm{)}}^{\mathrm{4}}}\mathrm{·}$(5)The variances of the areas, ${\mathit{\sigma }}_{{\mathit{A}}_{\mathrm{e}}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{{\mathit{A}}_{\mathrm{c}}}^{\mathrm{2}}$, characterize the jitter of the histograms. If the jitter is large, σ_ξ is large compared to | ξ | and the relative orientation is indeterminate. The jitter depends on the number of bins in the histogram, but ξ does not.

Figure 7 illustrates the change in ξ as a function of log ₁₀(N_H/ cm^-2) of each bin. For most of the clouds, ξ is positive in the lowest and intermediate N_H bins and negative or close to zero in the highest bins. The most pronounced changes in ξ from positive to negative are seen in Taurus, Chamaeleon-Musca, Aquila Rift, Perseus, and IC 5146.

The trend in ξ vs. log ₁₀(N_H/ cm^-2) can be fit roughly by a linear relation $\mathit{\xi }\mathrm{=}{\mathit{C}}_{\mathrm{HRO}}\hspace{0.17em}\mathrm{[}{\log }_{\mathrm{10}}\mathrm{(}{\mathit{N}}_{\mathrm{H}}\mathit{/}{\mathrm{cm}}^{-2}\mathrm{)}\mathrm{-}{\mathit{X}}_{\mathrm{HRO}}\mathrm{]}\mathit{.}$(6)The values of C_HRO and X_HRO in the regions analysed are summarized in Table 2. For the clouds with a pronounced change in relative orientation the slope C_HRO is steeper than about − 0.5, and the value X_HRO for the log ₁₀(N_H/ cm^-2) at which the relative orientation changes from parallel to perpendicular is greater than about 21.7. Ophiuchus, Lupus, Cepheus, and Orion are intermediate cases, where ξ definitely does not go negative in the data, but seems to do so by extrapolation; these tend to have a shallower C_HRO and/or a higher X_HRO.

The least pronounced change in ξ is seen in CrA, where ξ is consistently positive in all bins, and the slope is very flat. We applied the HRO analysis to a pair of test regions (Fig. 8) with even lower N_H values (log ₁₀(N_H/ cm^-2) < 21.6; see also Fig. 11 below) located directly south and directly east of the Chamaeleon-Musca region. As in CrA, we find that B_⊥ is mostly parallel to the N_H contours, a fairly constant ξ, and no indication of predominantly perpendicular relative orientation.

4.2. Comparisons with previous studies

The above trends in relative orientation between B_⊥ and the N_H contours in targeted MCs, where B_⊥ tend to become perpendicular to the N_H contours at high N_H, agree with the results of the Hessian matrix analysis applied to Planck observations over the whole sky, as reported in Fig. 15 of Planck Collaboration Int. XXXII (2016).

Evidence for preferential orientations in sections of some regions included in this study has been reported previously. The Taurus region has been the target of many studies. Moneti et al. (1984) and Chapman et al. (2011) find evidence using infrared polarization of background stars for a homogeneous magnetic field perpendicular to the embedded dense filamentary structure. High-resolution submillimetre observations of intensity with Herschel find evidence of faint filamentary structures (“striations”), which are well correlated with the magnetic field orientation inferred from starlight polarization (Palmeirim et al. 2013) and perpendicular to the filament B211. Heyer et al. (2008) report a velocity anisotropy aligned with the magnetic field, which can be interpreted as evidence of the channeling effect of the magnetic fields. But the magnetic field in B211 and the dense filamentary structures are not measured directly. As described above, using Planck polarization we find that B_⊥ is mostly aligned with the lowest N_H contours (20.8 < log ₁₀(N_H/ cm^-2) < 21.3; see also Planck Collaboration Int. XXXIII 2016), although the aligned structures do not correspond to the striations, which are not resolved at 10′ resolution. However, we also find that at higher N_H the relative orientation becomes perpendicular.

Similar studies in other regions have found evidence of striations correlated with the starlight-inferred magnetic field orientation and perpendicular to the densest filamentary structures. The regions studied include Serpens South (Sugitani et al. 2011), which is part of Aquila Rift in this study, Musca (Pereyra & Magalhães 2004), and the Northern Lupus cloud (Matthews et al. 2014). By studying Planck polarization in larger regions around the targets of these previous observations, we show that a systematic change in relative orientation is the prevailing statistical trend in clouds that reach log ₁₀(N_H/ cm^-2) ≳ 21.5.

5. Discussion

5.1. The relative orientation between B_⊥ and N_H structures

5.1.1. Spatial distribution of the HRO signal

The maps obtained using Eq. (2)characterize the relative orientation in each region, without assuming an organization of the N_H structures in ridges or filaments; HROs basically just sample the orientation of N_H contours. However, the resulting maps of the relative orientation angle, shown for the Taurus region in Fig. 9, reveal that the regions that are mostly oriented parallel or perpendicular to the field form continuous patches, indicating that the HRO signal is not only coming from variations in the field or the N_H contours at the smallest scales in the map.

HRO analysis on larger scales in a map can be achieved by considering a larger vicinity of pixels for calculating the gradient ∇τ₃₅₃. This operation is equivalent to calculating the next-neighbour gradient on a map first smoothed to the scale of interest. The results for relative orientation after smoothing to resolutions of 15′, 30′, and 60′ are illustrated for the Taurus region in Fig. 9. Figure 10 shows that the corresponding HROs have a similar behaviour for the three representative N_H values. These results confirm that the preferred relative orientation is not particular to the smallest scales in the map, but corresponds to coherent structures in N_H.

A study of the preferred orientation for the whole cloud would be possible by smoothing the column density and polarization maps to a scale comparable to the cloud size. The statistical significance of such a study would be limited to the number of clouds in the sample and would not be directly comparable to previous studies of relative orientation of clouds, where elongated structures were selected to characterize the mean orientation of each cloud (Li et al. 2013).

5.1.2. Statistical significance of the HRO signal

Our results reveal a systemic change of ξ with N_H, suggesting a systematic transition from magnetic field mostly parallel to N_H contours in the lowest N_H bins to mostly perpendicular in the highest N_H bins of the clouds studied. The statistical significance of this change can in principle be evaluated by considering the geometrical effects that influence this distribution. In Appendix C.1, using simulations of Q and U maps, we eliminate the possibility that this arises from random magnetic fields, random spatial correlations in the field, or the large scale structure of the field. In Appendix C.2 we simply displace the Q and U maps spatially and repeat the analysis, showing that the systemic trend of ξ vs. N_H disappears for displacements greater than 1°.

Using a set of Gaussian models, Planck Collaboration Int. XXXII (2016) estimated the statistical significance of this transition in terms of the relative orientation between two vectors in 3D and their projection in 2D. As these authors emphasized, two vectors that are close to parallel in 3D would be projected as parallel in 2D for almost all viewing angles for which the projections of both vectors have a non-negligible length, but on the other hand, the situation is more ambiguous seen in projection for two vectors that are perpendicular in 3D, because they can be projected as parallel in 2D depending on the angle of viewing. The quantitative effects are illustrated by the simulations in Appendix C.3, where we consider distributions of vectors in 3D that are mostly parallel, mostly perpendicular, or have no preferred orientation. The projection tends to make vector pairs look more parallel in 2D, but the distribution of relative orientations in 2D is quite similar though not identical to the distribution in 3D. In particular, the signal of perpendicular orientation is not erased and we can conclude that two projected vectors with non-negligible lengths that are close to perpendicular in 2D must also be perpendicular in 3D. This would apply to the mostly perpendicular orientation for the highest bin of N_H in the Taurus region.

5.2. Comparison with simulations of MHD turbulence

As a complement to observations, MHD simulations can be used to directly probe the actual 3D orientation of the magnetic field B with respect to the density structures. The change in the relative orientation with N_H was previously studied in MHD simulations using the inertia matrix and the HRO analysis (Hennebelle 2013; Soler et al. 2013). Soler et al. (2013) showed that in 3D, the change in the relative orientation is related to the degree of magnetization. If the magnetic energy is above or comparable to the kinetic energy (turbulence that is sub-Alfvénic or close to equipartition), the less dense structures tend to be aligned with the magnetic field and the orientation progressively changes from parallel to perpendicular with increasing density. In the super-Alfvénic regime, where the magnetic energy is relatively low, there appears to be no change in relative orientation with increasing density, with B and density structures being mostly parallel.

Soler et al. (2013) describe 2D synthetic observations of the MHD simulations. The synthetic observations are produced by integrating the simulation cubes along a direction perpendicular to the mean magnetic field and assuming a homogeneous dust grain alignment efficiency ϵ = 1.0. The angular resolution of the simulation is obtained by assuming a distance d = 150 pc and convolving the projected map with a Gaussian beam of 10′ full width at half maximum (FWHM). The trends in the relative orientation with N_H seen in 3D are also seen using the these 2D synthetic observations. Given that sub-Alfvénic or close to Alfvénic turbulence does not significantly disturb the well-ordered mean magnetic field, the orientation of B perpendicular to the iso-density contours is projected well for lines of sight that are not close to the mean magnetic field orientation. In contrast, the projected relative orientation produced by super-Alfvénic turbulence does not necessarily reflect the relative orientation in 3D as a result of the unorganized field structure.

The direct comparison between the HROs of the regions in this study and of the synthetic observations is presented in Fig. 11. The trends in the relative orientation parameter, ξ, show that the simulation with super-Alfvénic turbulence does not undergo a transition in relative orientation from parallel to perpendicular for log ₁₀(N_H/ cm^-2) < 23. In contrast, most of the observed clouds show a decrease in ξ with increasing N_H, close to the trends seen from the simulations with Alfvénic or sub-Alfvénic turbulence for log ₁₀(N_H/ cm^-2) < 23. Furthermore, X_HRO, the value of log ₁₀(N_H/ cm^-2) where ξ goes through zero, is near 21.7, which is consistent with the behaviour seen in the simulations with super-Alfvénic or Alfvénic turbulence. Given that the physical conditions in the simulations (σ_v = 2.0 km s^-1 and n = 500 cm^-3) are typical of those in the selected regions (σ_v∥ is given in Table D.1), the similarities in the dependence of ξ on N_H suggest that the strength of the magnetic field in most of the regions analysed would be about the same as the mean magnetic fields in the Alfvénic and sub-Alfvénic turbulence simulations, which are 3.5 and 11 μG, respectively. However, more precise estimates of the magnetic field strength coming directly from the HROs would require further sampling of the magnetization in the MHD simulations and detailed modelling of the effects of the line-of-sight integration.

Indirectly, the presence of a N_H threshold in the switch in preferred relative orientation between B_⊥ and the N_H structures hints that gravity plays a significant role. By contrast, for the regions with low average N_H (i.e., CrA and the two test regions; Fig. 11), there is little change in ξ and certainly no switch in the preferred relative orientation to perpendicular.

5.3. Physics of the relative orientation

The finding of dense structures mostly perpendicular to the magnetic field (and the small mass-to-flux ratios discussed in Appendix D.4) suggests that the magnetic field in most of the observed regions is significant for the structure and dynamics. However, discerning the underlying geometry is not obvious. As one guide, a frozen-in and strong interstellar magnetic field would naturally cause a self-gravitating, static cloud to become oblate, with its major axis perpendicular to the field lines, because gravitational collapse would be restricted to occurring along field lines (Mouschovias 1976a,b). In the case of less dense structures that are not self-gravitating, the velocity shear can stretch matter and field lines in the same direction, thereby producing aligned structures, as discussed in Hennebelle (2013) and Planck Collaboration Int. XXXII (2016).

If the MCs are isolated entities and the magnetic field is strong enough to set a preferred direction for the gravitational collapse, the condensations embedded in the cloud are not very likely to have higher column densities than their surroundings (Nakano 1998). This means that the formation of dense substructures, such as prestellar cores and stars, by gravitational collapse would be possible only if the matter decouples from the magnetic field. This is possible through the decoupling between neutral and ionized species (ambipolar diffusion, Mouschovias 1991; Li & Houde 2008) or through removal of magnetic flux from clouds via turbulent reconnection (Lazarian & Vishniac 1999; Santos-Lima et al. 2012).

Alternatively, if we regard MCs not as isolated entities but as the result of an accumulation of gas by large-scale flows (Ballesteros-Paredes et al. 1999; Hartmann et al. 2001; Koyama & Inutsuka 2002; Audit & Hennebelle 2005; Heitsch et al. 2006), the material swept up by colliding flows may eventually form a self-gravitating cloud. If the magnetic field is strong the accumulation of material is favoured along the magnetic field lines, thus producing dense structures that are mostly perpendicular to the magnetic field. The inflow of material might eventually increase the gravitational energy in parts of the cloud, thereby producing supercritical structures such as prestellar cores.

For supersonic turbulence in the ISM and MCs, density structures can be formed by gas compression in shocks. If the turbulence is strong with respect to the magnetic field (super-Alfvénic), gas compression by shocks is approximately isotropic; because magnetic flux is frozen into matter, field lines are dragged along with the gas, forming structures that tend to be aligned with the field. If the turbulence is weak with respect to the magnetic field (sub-Alfvénic), the fields produce a clear anisotropy in MHD turbulence (Sridhar & Goldreich 1994; Goldreich & Sridhar 1995; Matthaeus et al. 2008; Banerjee et al. 2009) and compression by shocks that is favored to occur along the magnetic field lines, creating structures perpendicular to the field. The cold phase gas that constitutes the cloud receives no information about the original flow direction because the magnetic field redistributes the kinetic energy of the inflows (Heitsch et al. 2009; Inoue & Inutsuka 2009; Burkhart et al. 2014). This seems to be the case in most of the observed regions, where the mostly perpendicular relative orientation between the magnetic field and the high column density structures is an indication of the anisotropy produced by the field.

The threshold of log ₁₀(N_H/ cm^-2) ≈ 21.7 above which the preferential orientation of B_⊥ switches to being perpendicular to the N_H contours is intriguing. Is there a universal threshold column density that is independent of the particular MC environment and relevant in the context of star formation? In principle, this threshold might be related to the column density of filaments at which substructure forms, as reported in an analysis of Herschel observations (Arzoumanian et al. 2013), but the Planck polarization observations leading to B_⊥ do not fully resolve such filamentary structures. In principle, this threshold might also be related to the column density at which the magnetic field starts scaling with density, according to the Zeeman effect observations of B_∥ (Fig. 7 in Crutcher 2012). However, establishing such relationships requires further studies with MHD simulations to identify what densities and scales influence the change in relative orientation between B_⊥ and N_H structures and to model the potential imprint in B_∥ observations and in B_⊥ observations to be carried out at higher resolution.

5.4. Effect of dust grain alignment

Throughout this study we assume that the polarized emission observed by Planck at 353 GHz is representative of the projected morphology of the magnetic field in each region; i.e., we assume a constant dust grain alignment efficiency (ϵ) that is independent of the local environment. Indeed, observations and MHD simulations under this assumption (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015) indicate that depolarization effects at large and intermediate scales in MCs might arise from the random component of the magnetic field along the line of sight. On the other hand, the sharp drop in the polarization fraction at N_H> 10²² cm^-2 (reported in Planck Collaboration Int. XIX 2015), when seen at small scales, might be interpreted in terms of a decrease of ϵ with increasing column density (Matthews et al. 2001; Whittet et al. 2008).

A leading theory for the process of dust grain alignment involves radiative torques by the incident radiation (Lazarian & Hoang 2007; Hoang & Lazarian 2009; Andersson 2015). A critical parameter for this mechanism is the ratio between the dust grain size and the radiation wavelength. As the dust column density increases, only the longer wavelength radiation penetrates the cloud and the alignment decreases. Grains within a cloud (without embedded sources) should have lower ϵ than those at the periphery of the same cloud. There is evidence for this from near-infrared interstellar polarization and submillimetre polarization along lines of sight through starless cores (Jones et al. 2015), albeit on smaller scales and higher column densities than considered here. If ϵ inside the cloud is very low, the observed polarized intensity would arise from the dust in the outer layers, tracing the magnetic field in the “skin” of the cloud. Then the observed orientation of B_⊥ is not necessarily correlated with the column density structure, which is seen in total intensity, or with the magnetic field deep in the cloud.

Soler et al. (2013) presented the results of HRO analysis on a series of synthetic observations produced using models of how ϵ might decrease with increasing density. They showed that with a steep decrease there is no visible correlation between the inferred magnetic field orientation and the high-N_H structure, corresponding to nearly flat HROs.

The HRO analysis of MCs carried out here reveals a correlation between the polarization orientation and the column density structure. This suggests that the dust polarized emission samples the magnetic field structure homogeneously on the scales being probed at the resolution of the Planck observations or, alternatively, that the field deep within high-N_H structures has the same orientation of the field in the skin.

6. Conclusions

We have presented a study of the relative orientation of the magnetic field projected on the plane of the sky (B_⊥), as inferred from the Planck dust polarized thermal emission, with respect to structures detected in gas column density (N_H). The relative orientation study was performed by using the histogram of relative orientations (HRO), a novel statistical tool for characterizing extended polarization maps. With the unprecedented statistics of polarization observations in extended maps obtained by Planck, we analyze the HRO in regions with different column densities within ten nearby molecular clouds (MCs) and two test fields.

In most of the regions analysed we find that the relative orientation between B_⊥ and N_H structures changes systematically with N_H  from being parallel in the lowest column density areas to perpendicular in the highest column density areas. The switch occurs at log ₁₀(N_H/ cm^-2) ≈ 21.7. This change in relative orientation is particularly significant given that projection tends to produce more parallel pseudo-vectors in 2D (the domain of observations) than exist in 3D.

The HROs in these MCs reveal that most of the high N_H structures in each cloud are mostly oriented perpendicular to the magnetic field, suggesting that they may have formed by material accumulation and gravitational collapse along the magnetic field lines. According to a similar study where the same method was applied to MHD simulations, this trend is only possible if the turbulence is Alfvénic or sub-Alfvénic. This implies that the magnetic field is significant for the gas dynamics on the scales sampled by Planck. The estimated mean magnetic field strength is about 4 and 12 μG for the case of Alfvénic and sub-Alfvénic turbulence, respectively.

We also estimate the magnetic field strength in the MCs studied using the DCF and DCF+SF methods. The estimates found seem consistent with the above values from the HRO analysis, but given the assumptions and systematic effects involved, we recommend that these rough estimates be treated with caution. According to these estimates the analysed regions appear to be magnetically sub-critical. This result is also consistent with the conclusions of the HRO analysis. Specific tools, such as the DCF and DCF+SF methods, are best suited to the scales and physical conditions in which their underlying assumptions are valid. The study of large polarization maps covering multiple scales calls for generic statistical tools, such as the HRO, for characterizing their properties and establishing a direct relation to the physical conditions included in MHD simulations.

The study of the structure on smaller scales is beyond the scope of this work, however, the presence of gravitationally bound structures within the MCs, such as prestellar cores and stars, suggests that the role of magnetic fields is changing on different scales. Even if the magnetic field is important in the accumulation of matter that leads to the formation of the cloud, effects such as matter decoupling from the magnetic field and the inflow of matter from the cloud environment lead to the formation of magnetically supercritical structures on smaller scales. Further studies will help to identify the dynamical processes that connect the MC structure with the process of star formation.

Appendix A: Selection of data

The HRO analysis is applied to each MC using common criteria for selecting the areas in which the relative orientation is to be assessed.

Appendix A.1: Gradient mask

The dust optical depth, τ₃₅₃, observed in each region, can be interpreted as ${\mathit{\tau }}_{\mathrm{353}}^{{obs}}\mathrm{=}{\mathit{\tau }}_{\mathrm{353}}^{{mc}}\mathrm{+}{\mathit{\tau }}_{\mathrm{353}}^{{bg}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}\mathit{,}$(A.1)where ${\mathit{\tau }}_{\mathrm{353}}^{{mc}}$ is the optical depth of the MC, ${\mathit{\tau }}_{\mathrm{353}}^{{bg}}$ is the optical depth of the diffuse regions behind and/or in front of the cloud (background/foreground), and δ_τ353 the noise in the optical depth map with variance ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{353}}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}^{\mathrm{2}}}$.

The gradient of the optical depth can be then written as $\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}^{{obs}}\mathrm{=}\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}^{{mc}}\mathrm{+}\mathrm{\nabla }\mathrm{(}{\mathit{\tau }}_{\mathrm{353}}^{{bg}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}\mathrm{)}\mathit{.}$(A.2)We quantify the contribution of the background/foreground and the noise, $\mathrm{\nabla }\mathrm{(}{\mathit{\tau }}_{\mathrm{353}}^{{bg}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}\mathrm{)}$, by evaluating ∇τ₃₅₃ in a reference field with lower submillimetre emission. Given that the dominant contribution to the background/foreground gradient would come from the gradient in emission from the Galactic plane, for each of the regions analysed we chose a reference field of the same size at the same Galactic latitude and with the lowest average N_H in the corresponding latitude band. We compute the average of the gradient norm in the reference field, ${\mathrm{⟨}}^{}\mathrm{|}\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}^{{ref}}{\mathrm{|}}^{\mathrm{⟩}}$, and use this value as a threshold for selecting the regions of the map where ∇τ₃₅₃ carries significant information about the structure of the cloud. We note that this threshold includes a contribution from the noise ∇δ_τ353. The HROs presented in this study correspond to regions in each field where $\mathrm{\left|}\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}\mathrm{\right|}\mathit{>}{}^{\mathrm{⟨}}\mathrm{|}\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}^{{ref}}{\mathrm{|}}^{\mathrm{⟩}}$.

Appendix A.2: Polarization mask

The total Stokes parameters Q and U measured in each region can be interpreted as ${\mathit{Q}}^{{obs}}\mathrm{=}{\mathit{Q}}^{{mc}}\mathrm{+}{\mathit{Q}}^{{bg}}\mathrm{+}{\mathit{\delta }}_{\mathrm{Q}}\hspace{0.17em}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{U}}^{{obs}}\mathrm{=}{\mathit{U}}^{{mc}}\mathrm{+}{\mathit{U}}^{{bg}}\mathrm{+}{\mathit{\delta }}_{\mathrm{U}}\mathit{,}$(A.3)where Q^mc and U^mc correspond to the polarized emission from the MC, Q^bg and U^bg correspond to the polarized emission from the diffuse background/foreground, and δ_Q and δ_U are the noise contributions to the observations, such that the variances ${\mathit{\sigma }}_{\mathit{Q}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{\delta }}_{\mathit{Q}}^{\mathrm{2}}}$ and ${\mathit{\sigma }}_{\mathit{U}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{\delta }}_{\mathit{U}}^{\mathrm{2}}}$.

As in the treatment of the gradient, we estimate the contributions of the background/foreground polarized emission and the noise using the rms of the Stokes parameters in the same reference field, ${\mathit{Q}}_{\mathrm{rms}}^{{ref}}$ and ${\mathit{U}}_{\mathrm{rms}}^{{ref}}$. The HROs presented in this study correspond to pixels in each region where $\mathrm{\left|}\mathit{Q}\mathrm{\right|}\mathit{>}\mathrm{2}{\mathit{Q}}_{\mathrm{rms}}^{{ref}}$ or $\mathrm{\left|}\mathit{U}\mathrm{\right|}\mathit{>}\mathrm{2}{\mathit{U}}_{\mathrm{rms}}^{{ref}}$. The “or” conditional avoids biasing the selected values of polarization. This first selection criterion provides a similar sample to the alternative coordinate-independent criterion $\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}}\mathit{>}\mathrm{2}\sqrt{\mathrm{(}{\mathit{Q}}_{\mathrm{rms}}^{{ref}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{U}}_{\mathrm{rms}}^{{ref}}{\mathrm{)}}^{\mathrm{2}}}$. This first criterion aims to distinguish between the polarized emission coming from the cloud and the polarized emission coming from the background/foreground estimated in the reference regions.

Additionally, as a second criterion, our sample is restricted to polarization measurements where | Q | > 3σ_Q or | U | > 3σ_U. This aims to select pixels where the uncertainty in the polarization angle is smaller than the size of the angle bins used for the constructions of the HRO (Serkowski 1958; Montier et al. 2015). In terms of the total polarized intensity, $\mathit{P}\mathrm{=}\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}}$, and following Eqs. (B.4) and (B.5) in Planck Collaboration Int. XIX (2015), the second criterion corresponds to P/σ_P> 3 and uncertainties in the polarized orientation angle σ_ψ< 10°.

The fractions of pixels considered in each region, after applying the selection criteria described above, are summarized in Table A.1. The largest masked portions of the regions correspond to the gradient mask, which selects mostly those areas of column density above the mean column density of the background/foreground ${}^{\mathrm{⟨}}{\mathit{N}}_{\mathrm{H}}^{{bg}}^{\mathrm{⟩}}$. The polarization mask provides an independent criterion that is less restrictive. The intersection of these two masks selects the fraction of pixels considered for the HRO analysis.

Appendix B: Construction of the histogram of relative orientations and related uncertainties

The HROs were calculated for 25 column density bins having equal numbers of selected pixels (10 bins in two regions with fewer pixels, CrA and IC 5146). For each of these HROs we use 12 angle bins of width 15° (see Sect. 4.1.1).

Appendix B.1: Calculation of ∇τ₃₅₃ and the uncertainty of its orientation

The optical depth gradient (∇τ₃₅₃) is calculated by convolving the τ₃₅₃ map with a Gaussian derivative kernel (Soler et al. 2013), such that ∇τ₃₅₃ corresponds to $\mathrm{\nabla }{\mathit{\tau }}_{\mathrm{353}}\mathrm{=}\mathrm{(}{\mathit{G}}_{\mathit{x}}\mathrm{\otimes }{\mathit{\tau }}_{\mathrm{353}}\mathrm{)}\hspace{0.17em}{î}\mathrm{+}\mathrm{(}{\mathit{G}}_{\mathit{y}}\mathrm{\otimes }{\mathit{\tau }}_{\mathrm{353}}\mathrm{)}\hspace{0.17em}{ĵ}\mathrm{=}{\mathit{g}}_{\mathit{x}}\hspace{0.17em}{î}\mathrm{+}{\mathit{g}}_{\mathit{y}}\hspace{0.17em}{ĵ}\mathit{,}$(B.1)where G_x and G_y are the kernels calculated using the x- and y-derivatives of a symmetric two-dimensional Gaussian function. The orientation of the iso-τ₃₅₃ contour is calculated from the components of the gradient vector, $\mathit{\theta }\mathrm{=}\arctan \mathrm{(}\mathrm{-}{\mathit{g}}_{\mathit{x}}\hspace{0.17em}\mathit{,}\hspace{0.17em}{\mathit{g}}_{\mathit{y}}\mathrm{)}\mathit{.}$(B.2)Because the calculation of the gradient through convolution is a linear operation, the associated uncertainties can be calculated using the same operation,so that ${\mathit{\delta }}_{{\mathit{g}}_{\mathit{x}}}\mathrm{=}{\mathit{G}}_{\mathit{x}}\mathrm{\otimes }{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}\hspace{0.17em}\mathrm{;}\hspace{0.17em}\hspace{0.17em}{\mathit{\delta }}_{{\mathit{g}}_{\mathit{y}}}\mathrm{=}{\mathit{G}}_{\mathit{y}}\mathrm{\otimes }{\mathit{\delta }}_{{\mathit{\tau }}_{\mathrm{353}}}\mathit{,}$(B.3)from which we obtain the ${\mathit{\sigma }}_{{\mathit{g}}_{\mathit{x}}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{\delta }}_{{\mathit{g}}_{\mathit{x}}}^{\mathrm{2}}}$ and ${\mathit{\sigma }}_{{\mathit{g}}_{\mathit{y}}}^{\mathrm{2}}\mathrm{=}\overline{{\mathit{\delta }}_{{\mathit{g}}_{\mathit{y}}}^{\mathrm{2}}}$.

The standard deviation of the angle θ can be written as ${\mathit{\sigma }}_{\mathit{\theta }}\mathrm{=}\sqrt{{\left(\frac{\mathit{\partial \theta }}{\mathit{\partial }{\mathit{g}}_{\mathit{x}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{g}}_{\mathit{x}}}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\partial \theta }}{\mathit{\partial }{\mathit{g}}_{\mathit{y}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{g}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{,}$(B.4)which corresponds to ${\mathit{\sigma }}_{\mathit{\theta }}\mathrm{=}\frac{\mathrm{1}}{{\mathit{g}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{g}}_{\mathit{y}}^{\mathrm{2}}}\hspace{0.17em}\sqrt{{\mathit{g}}_{\mathit{y}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{g}}_{\mathit{x}}}^{\mathrm{2}}\mathrm{+}{\mathit{g}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{g}}_{\mathit{y}}}^{\mathrm{2}}}\mathit{.}$(B.5)In the application discussed here, the standard deviations in the τ₃₅₃ map within the selected areas are much less than a few percentage points, so their effect on the estimate of the orientation of the gradient is negligible.

Appendix B.2: Uncertainties affecting the characterization of relative orientations within MCs

Appendix B.2.1: Uncertainties in the construction of the histogram

To estimate the uncertainty associated with the noise in Stokes Q and U, we produce 1000 noise realizations, Q_r and U_r, using Monte Carlo sampling. We assume that the errors are normally distributed and are centred on the measured values Q and U with dispersions σ_Q and σ_U (Planck Collaboration II 2014; Planck Collaboration VI 2014; Planck Collaboration V 2014; Planck Collaboration VIII 2014). Given that σ_QU is smaller than ${\mathit{\sigma }}_{\mathit{Q}}^{\mathrm{2}}$ and ${\mathit{\sigma }}_{\mathit{U}}^{\mathrm{2}}$, it is justified to generate Q_r and U_r independently of each other. We then introduce Q_r and U_r in the analysis pipeline and compute the HRO using the corresponding τ₃₅₃ map in each region. The results, presented in Fig. B.1 for the Taurus region, show that the noise in Q and U does not critically affect the shape of the HROs or the trend in ξ. Together, the low noise in the maps of τ₃₅₃ and the selection criteria for the polarization measurements ensure that the HRO is well determined.

Another source of uncertainty in the HRO resides in the histogram binning process. The variance in the kth histogram bin is given by ${\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}{\mathit{h}}_{\mathit{k}}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{h}}_{\mathit{k}}}{{\mathit{h}}_{\mathrm{tot}}}\right)\mathit{,}$(B.6)where h_k is the number of samples in the kth bin, and h_tot is the total number of samples.

Of the two independent sources of uncertainty above, we find that the largest contribution comes from the binning process, so that these are the ones shown as the shaded uncertainty ranges in all figures of HROs, for example for Taurus in Fig. 3. Because of the large number of samples in each histogram bin, the uncertainties in the HRO do not significantly affect the results of this study.

Appendix B.2.2: Uncertainties in the histogram shape parameter ξ

As in the case of the HRO, the uncertainty in ξ, as defined in Eq. (4), can be quantified using the random realizations introduced in the previous section. Figure B.2 shows the dependence of ξ on log ₁₀N_H obtained using Q_r and U_r for the Taurus region. The small variations around the trend line indicate that the uncertainties in Q and U do not significantly affect the trends discussed in this study, as expected from the behaviour of the histograms presented in Fig. B.1. The main source of uncertainty in the estimation of ξ is related to the histogram binning, characterized by the error bars calculated using Eq. (5). As seen in Fig. 7 and reproduced in Fig. C.1, these are much larger than the dispersion of the values of ξ in Fig. B.2.

Appendix C: Statistical significance of the HRO signal

Appendix C.1: A product of chance?

To investigate various potential sources of the signal found in the HROs we use the τ₃₅₃ map at 10′ resolution in each region in combination with Q and U maps that are produced with different recipes, each with 1000 realizations.

To eliminate random fields, we use Q and U maps produced with a random realization of ψ with a uniform distribution and with unit-length polarization pseudo-vectors. The results of this numerical experiment are shown in Fig. C.1 for the Taurus region. For each of the 1000 realizations, we join the values at each N_H bin to show the trend lines in order to compare with the lines in Fig. 7.

To eliminate the large-scale magnetic field as the source we use Q and U maps produced with random realizations of ψ with a Gaussian distribution and unit-length polarization pseudo-vectors. The polarization angle distribution is centred at ψ₀ = 0° with a standard deviation ς_ψ = 45°. The results of this numerical experiment are shown in Fig. C.2 for the Taurus region.

To eliminate random spatial correlations, we use Q and U maps produced from random realizations with a power spectrum P(k) ∝ k^α_m. For the spectral indices we adopted, α_m = −1.5, − 2.5, and − 3.5, correlations are introduced in the orientation of the polarization pseudo-vectors that are independent of the structure of matter. This test evaluates the statistical significance of the random sampling of spatial correlations in the magnetic field implied when we calculate the HROs in a finite region of the sky. The results of this numerical experiment are shown in Fig. C.3 for the Taurus region.

The results of these numerical experiments show that the trends with N_H found in the HROs and ξ do not arise by chance from these potential sources.

Appendix C.2: A product of random correlations in the polarization maps?

The random realizations of the magnetic field presented above are useful for characterizing the behaviour of ξ. However, they are bound to produce very different statistics compared to those in the real observations. In reality the orientations of B_⊥ have some non-trivial correlation structure in the map. In principle, the difference in the HROs might be due to the different correlations of the observed B_⊥ and of the random realizations of B_⊥.

To evaluate whether the signal present in the HROs arises from an actual physical relationship between B_⊥ and τ₃₅₃, we introduce randomness by shifting the Stokes Q and U maps with respect to the τ₃₅₃ map and then calculate the corresponding HROs and ξ as a function of log ₁₀(N_H/ cm^-2). The intrinsic statistical properties of the two maps are unchanged because the two maps are unchanged, only shifted.

If the trend in ξ as a function of log ₁₀(N_H/ cm^-2) were arising randomly, these trends would be unchanged even for significant shifts. Instead the results of this experiment, illustrated in Fig. C.4 for the Taurus region, show that the trends tend to disappear with increasing values of the size of the shift. For shifts of about 1°, the correlation between the magnetic field and the matter is still present, as expected from the results presented in Fig. 10, but for larger shifts the correlation is lost and the trend does not survive. Over the many MCs studied, the nature of the trends at large shifts appears to be random.

Appendix C.3: Projection effects

We evaluate the statistical significance of the relative orientation between B_⊥ and ∇τ₃₅₃ by considering the distribution of relative orientations between two vectors in 3D space compared to the distribution of relative orientation between their projections in 2D.

The projection of a vector v onto a plane normal to the unit vector can be written as ${u}\mathrm{=}{v}\mathrm{-}\mathrm{(}{v}\mathrm{·}{n̂}\mathrm{)}\hspace{0.17em}{n̂}\mathit{.}$(C.1)

Consider two unit vectors in 3D, v₁ and v₂, separated by an angle α, such that $\cos \mathit{\alpha }\mathrm{=}{{v}}_{\mathrm{1}}\mathrm{·}{{v}}_{\mathrm{2}}\mathit{.}$(C.2)The angle β between the projections of these two vectors onto a plane normal to , u₁, and u₂, can be written as $\begin{array}{ccc}\cos \mathit{\beta }& \mathrm{=}& \frac{{{u}}_{\mathrm{1}}\mathrm{·}{{u}}_{\mathrm{2}}}{\mathrm{\left|}{{u}}_{\mathrm{1}}\mathrm{\right|}\mathrm{\left|}{{u}}_{\mathrm{2}}\mathrm{\right|}}\\ & \mathrm{=}& \frac{\mathrm{(}\hspace{0.17em}{{v}}_{\mathrm{1}}\mathrm{·}{{v}}_{\mathrm{2}}\mathrm{-}\mathrm{(}{{v}}_{\mathrm{1}}\mathrm{·}{n̂}\mathrm{)}\mathrm{(}{{v}}_{\mathrm{2}}\mathrm{·}{n̂}\mathrm{)}{\hspace{0.17em}}^{\mathrm{)}}}{{\left(\hspace{0.17em}{{v}}_{\mathrm{1}}\mathrm{·}{{v}}_{\mathrm{1}}\mathrm{-}\mathrm{(}{{v}}_{\mathrm{1}}\mathrm{·}{n̂}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\left(\hspace{0.17em}{{v}}_{\mathrm{2}}\mathrm{·}{{v}}_{\mathrm{2}}\mathrm{-}\mathrm{(}{{v}}_{\mathrm{2}}\mathrm{·}{n̂}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$(C.3)Given a particular distribution of angles between the vectors v₁ and v₂, this expression, which is solved numerically, is useful for evaluating the resulting distribution of angles between the projected vectors u₁ and u₂. Without any loss of generality, we can assume that v₁ is oriented along the axis of a spherical coordinate system, such that ${{v}}_{\mathrm{1}}\mathrm{=}\hat{{k}}$, and that v₂ is oriented at polar angle θ and azimuth φ, such that ${{v}}_{\mathrm{2}}\mathrm{=}\cos \mathit{\phi }\sin \mathit{\theta }\hspace{0.17em}{î}\mathrm{+}\sin \mathit{\phi }\sin \mathit{\theta }\hspace{0.17em}{ĵ}\mathrm{+}\cos \mathit{\theta }\hspace{0.17em}\hat{{k}}$. Then, we can simulate a distribution of cosα by simulating a distribution of cosθ and generating φ with a uniform distribution.

We thus generate a set of vectors v₂ that follow a particular distribution of cosα. We choose three examples: a uniform distribution of relative orientation between v₂ and v₁, v₂ vectors that are mostly parallel to v₁, and v₂ vectors that are mostly perpendicular to v₁. The last two are Gaussian distributions centred at cosα = 0 (for mostly parallel) or cosα = ± 1 (for mostly perpendicular, given the periodicity of α) with a dispersion ς_α. These distributions are shown in the top panel of Fig. C.5 for three values of ς_α.

Using Eq. (C.3)we calculate the distribution of the angles between the projected vectors and the results are shown in the bottom panel of Fig. C.5. Projection generally tends to make vector pairs look more parallel in 2D than they are in 3D, regardless of the orientation in 3D, in agreement with the results of the Gaussian models presented in Planck Collaboration Int. XXXII (2016). Nevertheless, the distribution of relative orientations in 2D is qualitatively similar (though not identical) to the distribution in 3D. In particular, the signature of perpendicular orientation is not erased. We conclude that the observation in 2D of B_⊥ close to perpendicular to the column density structures is a direct indicator of the perpendicular configuration of B with respect to the density structures in 3D. In contrast, the parallel alignment of B_⊥ with the structure, or no signs of preferential orientation, while suggestive of parallel orientation in 3D, does not unambiguously rule out the presence of some close-to-perpendicular orientations within the distribution in 3D.

Appendix D: Alternative estimates of magnetic field strength

Given the historic importance of the DCF method (Davis 1951; Chandrasekhar & Fermi 1953) and the related DCF+SF method (Hildebrand et al. 2009), we have used them to estimate the magnetic field strength in each region. We also evaluated the mass-to-flux ratio, which is important in investigating the stability against gravitational collapse. We now provide a critical assessment of the applicability of these results.

Appendix D.1: Davis-Chandrasekhar-Fermi method

The DCF method estimates the strength of B_⊥ in a region using the dispersion of the polarization angles ς_ψ. Assuming that the magnetic field is frozen into the gas and that the dispersion of the B_⊥ orientation angles (or equivalently ς_ψ) is due to transverse incompressible Alfvén waves, then ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}\mathrm{=}\sqrt{\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }}\frac{{\mathit{\sigma }}_{{\mathit{v}}_{\mathrm{\parallel }}}}{{\mathit{\varsigma }}_{\mathit{\psi }}}\mathit{,}$(D.1)where σ_v∥ is the dispersion of the radial velocity of the gas (Appendix D.3) and ρ the gas mass density.

We calculate ς_ψ directly from Stokes Q and U using ${\mathit{\varsigma }}_{\mathit{\psi }}\mathrm{=}\sqrt{⟨\mathrm{(}\mathrm{\Delta }\mathit{\psi }{\mathrm{)}}^{\mathrm{2}}⟩}$(D.2)and $\mathrm{\Delta }\mathit{\psi }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \left(\mathit{Q}⟨\mathit{U}⟩\mathrm{-}⟨\mathit{Q}⟩\mathit{U}\hspace{0.17em}\mathit{,}\hspace{0.17em}\mathit{Q}⟨\mathit{Q}⟩\mathrm{+}⟨\mathit{U}⟩\mathit{U}\right)\mathit{,}$(D.3)where ⟨... ⟩ denotes an average over the selected pixels in each map (Planck Collaboration Int. XIX 2015).

Appendix D.2: Davis-Chandrasekhar-Fermi plus structure function method

As described by Hildebrand et al. (2009), the DCF+SF method characterizes the magnetic field dispersion about local structured fields by considering the difference in angle, Δψ(ℓ) = ψ(x)−ψ(x+ℓ), between pairs of B_⊥ vectors separated by displacements ℓ in the plane of the sky. Assuming that the angle differences are statistically isotropic (i.e., they depend only on ℓ = | ℓ | and not on the orientation of ℓ), they can binned by distance, ℓ. From the N(ℓ) pairs of B_⊥ vectors for that bin, the square of the second-order structure function is ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathrm{⟨}\mathrm{\left[}\mathrm{\Delta }\mathit{\psi }\mathrm{\right(}{x}\mathit{,\ell }\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\mathrm{⟩}}_{\mathit{x}}\mathrm{=}{⟨\frac{\mathrm{1}}{\mathit{N}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}}\mathrm{(}\mathrm{\Delta }{\mathit{\psi }}_{\mathit{x,i}}{\mathrm{)}}^{\mathrm{2}}⟩}_{\mathit{x}}\mathit{,}$(D.4)as introduced by Kobulnicky et al. (1994) and Falceta-Gonçalves et al. (2008). In terms of the Stokes parameters, each term in the sum can be written $\mathrm{\Delta }{\mathit{\psi }}_{\mathit{x,i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \left({\mathit{Q}}_{\mathit{i}}{\mathit{U}}_{\mathit{x}}\mathrm{-}{\mathit{Q}}_{\mathit{x}}{\mathit{U}}_{\mathit{i}}\hspace{0.17em}\mathit{,}\hspace{0.17em}{\mathit{Q}}_{\mathit{i}}{\mathit{Q}}_{\mathit{x}}\mathrm{+}{\mathit{U}}_{\mathit{i}}{\mathit{U}}_{\mathit{x}}\right)\mathit{,}$(D.5)where the subscripts x and i represent the central and displaced positions, respectively.

Hildebrand et al. (2009) assume that B(x) is composed of a large-scale structured field, B₀(x), and a random component, B_r(x), which are statistically independent of each other. Assuming that B₀(x) is a smoothly varying quantity, its contribution to ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ should increase in proportion to ℓ² for distances that are much smaller than the scales at which B₀ itself fluctuates. Additionally, assuming that turbulence occurs on scales that are small enough to completely decorrelate B_r for the range of scales probed by the displacements ℓ, Hildebrand et al. (2009) derived the approximation ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathit{b}}^{\mathrm{2}}\mathrm{+}{\mathit{m}}^{\mathrm{2}}{\mathit{\ell }}^{\mathrm{2}}\mathit{,}$(D.6)where the two terms on the right-hand side give the contributions from the random and large-scale magnetic fields, respectively. These assumptions are not necessarily valid for the range of scales in the Planck data.

Hildebrand et al. (2009) used Eq. (D.6)to estimate b² as the intercept and then related b to the ratio of the random to the large-scale magnetic field strength, both projected onto the plane of the sky, through $\frac{{\mathrm{⟨}{{\mathit{B}}_{\mathrm{r}\mathit{,}\mathrm{\perp }}^{\mathrm{2}}}^{\mathrm{⟩}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{\perp }}}\mathrm{=}\frac{\mathit{b}}{\sqrt{\mathrm{2}\mathrm{-}{\mathit{b}}^{\mathrm{2}}}}\mathit{,}$(D.7)where ${{}^{\mathrm{⟨}}{\mathit{B}}_{\mathrm{r}\mathit{,}\mathrm{\perp }}^{\mathrm{2}}^{\mathrm{⟩}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ stands for the root mean square (rms) variations about the large-scale magnetic field, B_{0, ⊥}. The same assumptions that result in Eq. (D.1)– i.e., considering only incompressible and isotropic turbulence, magnetic fields frozen into the gas, and dispersion of the B_⊥ orientation originating in transverse incompressible Alfvén waves – lead to $\frac{{\mathrm{⟨}{{\mathit{B}}_{\mathrm{r}\mathit{,}\mathrm{\perp }}^{\mathrm{2}}}^{\mathrm{⟩}}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{\perp }}}\mathrm{=}\frac{{\mathit{\sigma }}_{{\mathit{v}}_{\mathrm{\parallel }}}}{{\mathit{V}}_{{A}\mathit{,}\mathrm{\perp }}}\mathit{,}$(D.8)where ${\mathit{V}}_{{A}\mathit{,}\mathrm{\perp }}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{\perp }}}{\sqrt{\mathrm{4}\mathit{\pi \rho }}}$(D.9)is the speed of the transverse incompressible Alfvén waves. Combining Eqs. (D.7) – (D.9)results in ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}\mathrm{\equiv }{\mathit{B}}_{\mathrm{0}\mathit{,}\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}\mathrm{=}\sqrt{\mathrm{4}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}\mathit{\rho }}\hspace{0.17em}\frac{{\mathit{\sigma }}_{{\mathit{v}}_{\mathrm{\parallel }}}\hspace{0.17em}\sqrt{\mathrm{2}\mathrm{-}{\mathit{b}}^{\mathrm{2}}}}{\mathit{b}}\mathrm{·}$(D.10)

Appendix D.3: Calculation

The estimates of velocity dispersion σ_v∥ and mass density ρ are the same in the two methods (DCF and DCF+SF).

We obtain σ_v∥ from the most complete CO emission-line survey of the Milky Way currently available, that of Dame et al. (2001). This data set consists of 488 000 spectra that beam-sample (1/8°) a set of MCs and the Galactic plane over a strip from 4° to 10° wide in latitude. We find that material in the MCs has centroid velocities in the range − 10 <v_∥/ (km s^-1) < 10 in the first moment map. Furthermore, for the calculations below, we need to select pixels with sufficient S/N in the polarization observations and so use the second criterion described in Appendix A.2. This combination of cuts selects the areas illustrated in Fig. D.1. Over these areas we calculate the average of the velocity dispersion from the second moment map and use this value as σ_v∥. This is tabulated along with other properties in Table D.1.

The mass density, ρ = μnm_p, is the product of the proton mass, m_p, the mean molecular weight per hydrogen molecule, μ = 2.8, and the mean number density, n. We require a value of n, which for the discussions here we approximate to be 100 cm^-3. This is a typical value for MCs (Draine 2011) and in rough agreement with the column densities and cloud sizes presented in Table 1. In practice, n varies from one cloud to the next, and its estimation involves the assumption of a particular cloud geometry, resulting in additional uncertainties that are not considered in this study.

For the DCF method, the dispersion of the orientation angle ς_ψ corresponds to the centred second moment of the ψ distribution evaluated in pixels within the selected region. Following the analysis presented in Planck Collaboration Int. XIX (2015) for , the variance on ς_ψ can be expressed as ${\mathit{\sigma }}_{{\mathit{\varsigma }}_{\mathit{\psi }}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}^{\mathrm{2}}{\mathit{\varsigma }}_{\mathit{\psi }}^{\mathrm{2}}}\left[{\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\varsigma }}_{\mathit{i}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{\psi }}^{\mathrm{2}}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\varsigma }}_{\mathit{i}}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\psi }}_{\mathit{i}}}^{\mathrm{2}}\right]\mathit{,}$(D.11)where ψ_i and σ_ψi are the orientation angle and uncertainty for each polarization pseudo-vector, N is the number of ψ measurements, ς_i = ψ_i−⟨ψ⟩, and ⟨ψ⟩ is the average orientation angle in each region. Like other quadratic functions, ${\mathit{\sigma }}_{{\mathit{\varsigma }}_{\mathit{\psi }}}^{\mathrm{2}}$ is biased positively when noise is present, leading to an overestimation of this quantity. However, given that we limit our analysis to polarization measurements with high S/N and that the uncertainties in σ_v∥ are considerably larger, the bias correction does not have a significant effect on our estimates. We apply Eq. (D.1)to estimate ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}$ and propagate the errors to obtain the values listed in Table D.1.

For the DCF+SF method, to estimate the parameter b, we first calculate ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ for each of the regions. The resolution of the data used is 10′. We evaluate ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ in steps of for lags 0′  and in steps of for lags 40′ <ℓ< 200′. For each considered lag ℓ, Δψ_x,i is evaluated pixel by pixel, considering all the pixels located in an annulus with radius ℓ. The term Δψ_x,i is only considered in the calculation of ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ if there are at least three pixels in the annulus. As anticipated, the model from Eq. (D.6)does not agree with the data on all sampled scales, so the range of scales considered is limited to ℓ above the resolution of the data and between 50′ and 200′, where the behaviour of ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ is approximately linear in ℓ². We then make a linear fit of ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ using ℓ² as the independent variable and estimate the value of b and its uncertainty, σ_b. The results are plotted in Fig. D.2.

Using the calculated b values⁶, we apply Eq. (D.10)to estimate ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}$, and the results are listed in Table D.1. As expected, the DCF+SF method produces estimates of B_⊥ that are about twice as large as those obtained with the DCF method. However, it is worth noting that because of the discrepancy between the model and the observations, the value of b depends on the selected ℓ-range and so propagates into a different value of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}$.

Appendix D.4: Mass-to-flux ratios

The critical value for the mass that can be supported against gravity by a magnetic flux Φ can be estimated to first order for a uniform disk from (M/ Φ)_crit≡1 / (2πG^{1 / 2}) (Nakano & Nakamura 1978). The precise value of the right-hand side changes for different geometries (e.g., Spitzer 1968; McKee et al. 1993). Stability can be assessed using $\mathit{\lambda }\mathrm{=}\frac{\mathrm{(}\mathit{M}\mathit{/}\mathrm{\Phi }\mathrm{)}}{\mathrm{(}\mathit{M}\mathit{/}\mathrm{\Phi }{\mathrm{)}}_{\mathrm{crit}}}\mathrm{=}\mathrm{7.6}\mathrm{\times }{\mathrm{10}}^{-21}\hspace{0.17em}\frac{\mathrm{(}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{mft}}\mathit{/}\mathrm{cm}-2\mathrm{)}}{\mathrm{(}{\mathit{B}}^{\mathrm{mft}}\mathit{/}\mathit{\mu }\mathrm{G}\mathrm{)}}\mathit{,}$(D.12)where ${\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{mft}}$ and B^mft are the H₂ column density and magnetic field strength along a magnetic flux tube (Crutcher et al. 2004). A cloud is supercritical and prone to collapse under its own gravity, when λ> 1; otherwise, when λ< 1, the cloud is sub-critical, magnetically supported against gravitational collapse.

What is observable is λ_obs, in which ${\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{mft}}\mathit{/}{\mathit{B}}^{\mathrm{mft}}$ is replaced by N_H2/B_⊥. We evaluate λ_obs by combining the value of ⟨N_H2⟩ computed from the integrated CO line emission and the conversion factor X_CO = (1.8 ± 0.3) × 10²⁰ cm^-2K^-1 km^-1 s (Dame et al. 2001) and B_⊥ estimated with the DCF and DCF+SF methods. The X_CO factor may show cloud-to-cloud or regional variations (Draine 2011), but we consider that these are not significant in comparison to the uncertainties involved in the estimation of B_⊥. The calculated values of λ_obs are listed in Table D.1. They are consistent with being less than unity.

From this we might obtain λ by judicious deprojection, since $\mathit{\lambda }\mathrm{=}\mathrm{(}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{mft}}\mathit{/}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{)}\hspace{0.17em}\mathrm{(}{\mathit{B}}_{\mathrm{\perp }}\mathit{/}{\mathit{B}}^{\mathrm{mft}}\mathrm{)}\mathrm{\times }{\mathit{\lambda }}_{\mathrm{obs}}\mathrm{\equiv }{\mathit{f}}_{\mathrm{dp}}\hspace{0.17em}{\mathit{\lambda }}_{\mathrm{obs}}\mathit{.}$(D.13)Here, B_⊥ is always less than B^mft, pushing f_dp below unity. The situation for column density depends on the geometry of the structure relative to the magnetic field. For a structure with the magnetic field along the short axis, ${\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{mft}}\mathit{<}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}}$, again lowering f_dp.

Statistically, the mean mass-to-flux ratio can be related to the observed value by assuming a particular geometry of the cloud, an ellipsoid with equatorial radius a and centre-to-pole distance c, and a magnetic field oriented along the polar axis of the ellipsoid. For an oblate spheroid, flattened perpendicular to the orientation of the magnetic field, f_dp = 1 / 3, yielding $\overline{\mathit{M}\mathit{/}\mathrm{\Phi }}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\frac{\mathit{M}\cos \mathit{\beta }}{\mathrm{\Phi }\mathit{/}\sin \mathit{\beta }}\hspace{0.17em}\sin \mathit{\beta }\hspace{0.17em}\mathrm{d}\mathit{\beta }\hspace{0.17em}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\mathrm{(}\mathit{M}\mathit{/}\mathrm{\Phi }{\mathrm{)}}_{\mathrm{obs}}\mathit{,}$(D.14)where β is the inclination angle with respect to the line of sight, and the sinβ dependence in the flux comes from B = B_⊥sinβ, (Crutcher et al. 2004). For a sphere there is no cosβ dependence and f_dp = 1 / 2. For a prolate spheroid elongated along the orientation of the field, the mass is multiplied by sinβ instead of cosβ, resulting in f_dp = 3 / 4. Investigating which geometry, if any, is most relevant to the actual MCs and the magnetized structures detected within them is obviously important before firm conclusions can be drawn regarding gravitational instability.

Appendix D.5: Discussion

Our estimates of the magnetic field strengths in the MCs analysed and the mass-to-flux ratios, presented in Table D.1, stem from the classic calculation presented by Chandrasekhar & Fermi (1953) and an updated interpretation presented by Hildebrand et al. (2009). Both of these methods assume very specific conditions, so (as we discuss below) this limits the conclusions that can be drawn from these estimates of the magnetic field strength. The deduced mass-to-flux ratios suggest that the clouds are potentially magnetically sub-critical, but we again need to be cautious about drawing conclusions from this application of the DCF and DCF+SF analyses alone.

In the case of the DCF method, the values of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}$ are obtained by assuming that the structure of the magnetic field is the product only of incompressible Alfvén waves, where the displacements are perpendicular to the direction of propagation. This is not the case for turbulence in MCs where the random component of the magnetic field can have any orientation. The dispersion measured about mean fields, assumed to be straight in DCF, may be much larger than should be attributed to MHD waves or turbulence, leading to an overestimation of ς_ψ and to low values of B_⊥. Moreover true interstellar turbulence in MCs involves not only incompressible Alfvén waves, but also compressible magneto-sonic waves, which do not satisfy Eq. (D.1). Furthermore, depending on the scales examined, the magnetic field may have structures due to effects such as differential rotation, gravitational collapse, or expanding Hii regions.

In the case of the DCF+SF method, the values of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}$ are obtained by assuming a very specific model of the magnetic field, which is in principle just a first-order approximation. First, this model assumes that the effect of the large-scale structured magnetic field, B₀, is to cause the square of the second-order structure function, ${\mathit{S}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$, to increase as ℓ². This corresponds to a very specific correlation function for B₀. Second, this model assumes that the dispersion of the random component of the field, B_r, is scale-independent, which is not realistic for the range of scales probed by Planck (Elmegreen & Scalo 2004; Hennebelle & Falgarone 2012).

In addition, the magnetic field orientation deduced from the polarization angle in a particular direction is not generally that of the field at a single point along the line of sight. The observed polarization is the average of various field pseudo-vectors weighted by local dust emission along the line of sight. The net effect of the integration of multiple uncorrelated components along the line of sight is an observed dispersion of the polarization angle that is smaller than the true 3D dispersion of the magnetic field orientation, thus leading to an overestimation of B_⊥. Myers & Goodman (1991) presented an analysis of this effect in terms of the number of correlation lengths of the magnetic field along the line of sight through a cloud, which they calculated empirically.

Houde et al. (2009) presented an extension of the DCF+SF method that includes the effect of signal integration along the line of sight and across the area subtended by the telescope beam. The extended method, also implemented in Houde et al. (2011, 2013), is based on the identification of the magnetized turbulence correlation length (δ) by means of the structure function of the polarization angles. In the case of the Planck 353 GHz observations, the angular resolution is not sufficient to identify δ and the corrections would have to rely on rough estimates of this value and the depth of integration (Δ). Following Eq. (29) in Houde et al. (2009), rough estimates δ ≈ 0.2 pc and Δ ≈ 10 pc result in a few correlation lengths across the beam, corresponding to correction factors around 0.4. Coincidentally, such correction factors lead to values of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}$ close to those of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}$ in Table D.1. But note that this correction relies on specific assumptions on the nature of the turbulence correlation function, it does not circumvent the necessity to observe the dust polarized emission with higher angular resolution to fully characterize the magnetic field.

MHD simulations provide a potentially useful guide to what modifications might be introduced into the DCF formula to allow for inhomogeneity and line-of-sight averaging. Using synthetic observations of MHD simulations, Ostriker et al. (2001) showed that correcting Eq. (D.1)by a factor provides a good approximation to the actual magnetic field strength in cases where ς_ψ< 25°. However, the effect of nonlinear amplitudes is uncertain (Zweibel 1996), and the method fails for values of ς_ψ> 25°, which is the case for all of the regions in this study (Table D.1). Falceta-Gonçalves et al. (2008) propose a method that is potentially valid for any value of ς_ψ, based on a fit to the ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}$ values obtained from maps at different resolutions, again concluding that the field should be lower than estimated with Eq. (D.1). However, this approach was not tested in MHD simulations that include gravity, which is the critical process that we aim to evaluate by using the DCF method.

Another strong assumption is that the behaviour of the velocity and the magnetic field are represented well by the observed quantities σ_v∥ and ς_ψ for a particular set of scales, which might not necessarily be the case. Even if the power spectra of v and B are comparable in 3D, the integration along the line of sight is different for the two quantities. The dispersion σ_v∥ is based on the emission-line profile v_∥ tracing a gas species, and while the emission is directly integrated along the line of sight, the line profile is possibly affected by radiative transfer and excitation effects. The tracer of the magnetic field is the optically-thin polarized submillimetre emission of dust, and both the polarization and B_⊥ are projected and integrated along the line of sight as pseudo-vectors.

Finally, in the estimates of ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}}$ and ${\mathit{B}}_{\mathrm{\perp }}^{\mathrm{DCF}\mathrm{+}\mathrm{SF}}$, a common mean density ρ has been adopted, while in practice ρ is different from cloud to cloud. Direct estimation of the values of ρ from the measured column densities ⟨N_H⟩ and ⟨N_H2⟩ relies heavily on the geometrical modelling of the cloud and the filling factors of each species, introducing uncertainties that affect the calculated values of B_⊥.

Given the line-of-sight integration and the fact that ς_ψ> 25° uniformly, the calculated values could be considered as upper limits to the actual B_⊥. However, other shortcomings, such as the assumptions about the correlation structure and the uncertainties in the determination of the density, do not necessarily bias the estimate of the magnetic field strength towards high values. In conclusion, the values presented in Table D.1 should be viewed only as a reference and only applied with caution, given the many assumptions in both methods at the scales considered.

---

Acknowledgments

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). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 267934.

References

All Tables

All Figures