# Planck intermediate results XXXV. Probing the role of the magnetic field in the formation of structure in molecular clouds

Planck Collaboration, P. A. R. Ade, N. Aghanim, M. I. R. Alves, M. Arnaud, D. Arzoumanian, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Planck Collaboration, P. A. R. Ade, N. Aghanim, M. I. R. Alves, M. Arnaud, et al., Planck intermediate results XXXV. Probing the role of the magnetic field in the formation of structure in molecular clouds. Astronomy and Astrophysics - A&A, 2016, 586, pp.A138. 10.1051/0004-6361/201525896. [hal-02137621](https://hal.science/hal-02137621)

HAL Id: hal-02137621

https://hal.science/hal-02137621

Submitted on 19 Jan 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**Planck intermediate results**

XXXV. Probing the role of the magnetic field
in the formation of structure in molecular clouds

Planck Collaboration: P. A. R. Ade<sup>83</sup>, N. Aghanim<sup>56</sup>, M. I. R. Alves<sup>56</sup>, M. Arnaud<sup>70</sup>, D. Arzoumanian<sup>56</sup>, M. Ashdown<sup>66,5</sup>, J. Aumont<sup>56</sup>, C. Baccigalupi<sup>81</sup>, A. J. Banday<sup>89,8</sup>, R. B. Barreiro<sup>61</sup>, N. Bartolo<sup>27,62</sup>, E. Battaner<sup>90,91</sup>, K. Benabed<sup>57,88</sup>, A. Benoît<sup>54</sup>, A. Benoît-Lévy<sup>21,57,88</sup>, J.-P. Bernard<sup>89,8</sup>, M. Bersanelli<sup>30,46</sup>, P. Bielewicz<sup>78,8,81</sup>, J. J. Bock<sup>63,9</sup>, L. Bonavera<sup>61</sup>, J. R. Bond<sup>7</sup>, J. Borrill<sup>11,85</sup>, F. R. Bouchet<sup>57,84</sup>, F. Boulanger<sup>56</sup>, A. Bracco<sup>56</sup>, C. Burigana<sup>45,28,47</sup>, E. Calabrese<sup>87</sup>, J.-F. Cardoso<sup>71,1,57</sup>, A. Catalano<sup>72,69</sup>, H. C. Chiang<sup>24,6</sup>, P. R. Christensen<sup>79,33</sup>, L. P. L. Colombo<sup>20,63</sup>, C. Combet<sup>72</sup>, F. Couchot<sup>68</sup>, B. P. Crill<sup>63,9</sup>, A. Curto<sup>61,5,66</sup>, F. Cuttaia<sup>45</sup>, L. Danese<sup>81</sup>, R. D. Davies<sup>64</sup>, R. J. Davis<sup>64</sup>, P. de Bernardis<sup>29</sup>, A. de Rosa<sup>45</sup>, G. de Zotti<sup>42,81</sup>, J. Delabrouille<sup>1</sup>, C. Dickinson<sup>64</sup>, J. M. Diego<sup>61</sup>, H. Dole<sup>56,55</sup>, S. Donzelli<sup>46</sup>, O. Doré<sup>63,9</sup>, M. Douspis<sup>56</sup>, A. Ducout<sup>57,52</sup>, X. Dupac<sup>35</sup>, G. Efstathiou<sup>58</sup>, F. Elsner<sup>21,57,88</sup>, T. A. Enßlin<sup>76</sup>, H. K. Eriksen<sup>59</sup>, D. Falceta-Gonçalves<sup>82,34</sup>, E. Falgarone<sup>69</sup>, F. Ferrière<sup>89,8</sup>, F. Finelli<sup>45,47</sup>, O. Forni<sup>89,8</sup>, M. Frailis<sup>44</sup>, A. A. Fraisse<sup>24</sup>, E. Franceschi<sup>45</sup>, A. Frejsel<sup>79</sup>, S. Galeotta<sup>44</sup>, S. Galli<sup>65</sup>, K. Ganga<sup>1</sup>, T. Ghosh<sup>56</sup>, M. Giard<sup>89,8</sup>, E. Gjerløw<sup>59</sup>, J. González-Nuevo<sup>16,61</sup>, K. M. Górski<sup>63,92</sup>, A. Gregorio<sup>31,44,50</sup>, A. Gruppuso<sup>45</sup>, J. E. Gudmundsson<sup>24</sup>, V. Guillet<sup>56</sup>, D. L. Harrison<sup>58,66</sup>, G. Helou<sup>9</sup>, P. Hennebelle<sup>70</sup>, S. Henrot-Versillé<sup>68</sup>, C. Hernández-Monteagudo<sup>10,76</sup>, D. Herranz<sup>61</sup>, S. R. Hildebrandt<sup>63,9</sup>, E. Hivon<sup>57,88</sup>, W. A. Holmes<sup>63</sup>, A. Hornstrup<sup>13</sup>, K. M. Huffenberger<sup>22</sup>, G. Hurier<sup>56</sup>, A. H. Jaffe<sup>52</sup>, T. R. Jaffe<sup>89,8</sup>, W. C. Jones<sup>24</sup>, M. Juvela<sup>23</sup>, E. Keihänen<sup>23</sup>, R. Keskitalo<sup>11</sup>, T. S. Kisner<sup>74</sup>, J. Knoche<sup>76</sup>, M. Kunz<sup>14,56,2</sup>, H. Kurki-Suonio<sup>23,40</sup>, G. Lagache<sup>4,56</sup>, J.-M. Lamarre<sup>69</sup>, A. Lasenby<sup>5,66</sup>, M. Lattanzi<sup>28</sup>, C. R. Lawrence<sup>63</sup>, R. Leonardi<sup>35</sup>, F. Levrier<sup>69</sup>, M. Liguori<sup>27,62</sup>, P. B. Lilje<sup>59</sup>, M. Linden-Vørnle<sup>13</sup>, M. López-Caniego<sup>35,61</sup>, P. M. Lubin<sup>25</sup>, J. F. Macías-Pérez<sup>72</sup>, D. Maino<sup>30,46</sup>, N. Mandolesi<sup>45,28</sup>, A. Mangilli<sup>56,68</sup>, M. Maris<sup>44</sup>, P. G. Martin<sup>7</sup>, E. Martínez-González<sup>61</sup>, S. Masi<sup>29</sup>, S. Matarrese<sup>27,62,38</sup>, A. Melchiorri<sup>29,48</sup>, L. Mendes<sup>35</sup>, A. Mennella<sup>30,46</sup>, M. Migliaccio<sup>58,66</sup>, M.-A. Miville-Deschênes<sup>56,7</sup>, A. Moneti<sup>57</sup>, L. Montier<sup>89,8</sup>, G. Morgante<sup>45</sup>, D. Mortlock<sup>52</sup>, D. Munshi<sup>83</sup>, J. A. Murphy<sup>77</sup>, P. Naselsky<sup>79,33</sup>, F. Nati<sup>24</sup>, C. B. Netterfield<sup>17</sup>, F. Noviello<sup>64</sup>, D. Novikov<sup>75</sup>, I. Novikov<sup>79,75</sup>, N. Oppermann<sup>7</sup>, C. A. Oxborrow<sup>13</sup>, L. Pagano<sup>29,48</sup>, F. Pajot<sup>56</sup>, R. Paladini<sup>53</sup>, D. Paoletti<sup>45,47</sup>, F. Pasian<sup>44</sup>, L. Perotto<sup>72</sup>, V. Pettorino<sup>39</sup>, F. Piacentini<sup>29</sup>, M. Piat<sup>1</sup>, E. Pierpaoli<sup>20</sup>, D. Pietrobon<sup>63</sup>, S. Plaszczynski<sup>68</sup>, E. Pointecouteau<sup>89,8</sup>, G. Polenta<sup>3,43</sup>, N. Ponthieu<sup>56,51</sup>, G. W. Pratt<sup>70</sup>, S. Prunet<sup>57,88</sup>, J.-L. Puget<sup>56</sup>, J. P. Rachen<sup>18,76</sup>, M. Reinecke<sup>76</sup>, M. Remazeilles<sup>64,56,1</sup>, C. Renault<sup>72</sup>, A. Renzi<sup>32,49</sup>, I. Ristorcelli<sup>89,8</sup>, G. Rocha<sup>63,9</sup>, M. Rossetti<sup>30,46</sup>, G. Roudier<sup>1,69,63</sup>, J. A. Rubiño-Martín<sup>60,15</sup>, B. Rushholme<sup>53</sup>, M. Sandri<sup>45</sup>, D. Santos<sup>72</sup>, M. Savelainen<sup>23,40</sup>, G. Savini<sup>80</sup>, D. Scott<sup>19</sup>, J. D. Soler<sup>56,*</sup>, V. Stolyarov<sup>5,86,67</sup>, R. Sudiwala<sup>83</sup>, D. Sutton<sup>58,66</sup>, A.-S. Suur-Uski<sup>23,40</sup>, J.-F. Sygnet<sup>57</sup>, J. A. Tauber<sup>36</sup>, L. Terenzi<sup>37,45</sup>, L. Toffolatti<sup>16,61,45</sup>, M. Tomasi<sup>30,46</sup>, M. Tristram<sup>68</sup>, M. Tucci<sup>14</sup>, G. Umana<sup>41</sup>, L. Valenziano<sup>45</sup>, J. Valiviita<sup>23,40</sup>, B. Van Tent<sup>73</sup>, P. Vielva<sup>61</sup>, F. Villa<sup>45</sup>, L. A. Wade<sup>63</sup>, B. D. Wandelt<sup>57,88,26</sup>, I. K. Wehus<sup>63</sup>, N. Ysard<sup>23</sup>, D. Yvon<sup>12</sup>, and A. Zonca<sup>25</sup>

(Affiliations can be found after the references)

Received 13 February 2015 / Accepted 25 May 2015

ABSTRACT

Within ten nearby (d < 450 pc) Gould belt molecular clouds we evaluate statistically the relative orientation between the magnetic field projected on the plane of sky, inferred from the polarized thermal emission of Galactic dust observed by *Planck* at 353 GHz, and the gas column density structures, quantified by the gradient of the column density, $N_H$. The selected regions, covering several degrees in size, are analysed at an effective angular resolution of 10' FWHM, thus sampling physical scales from 0.4 to 40 pc in the nearest cloud. The column densities in the selected regions range from $N_H$ ≈ 10<sup>21</sup> to 10<sup>23</sup> cm<sup>-2</sup>, and hence they correspond to the bulk of the molecular clouds. The relative orientation is evaluated pixel by pixel and analysed in bins of column density using the novel statistical tool called "histogram of relative orientations". 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, independent of the local environment. Within most clouds we find that the relative orientation changes progressively with increasing $N_H$, from mostly parallel or having no preferred orientation to mostly perpendicular. In simulations of magnetohydrodynamic turbulence in molecular clouds this trend in relative orientation is a signature of Alfvénic or sub-Alfvénic turbulence, implying that the magnetic field is significant for the gas dynamics at the scales probed by *Planck*. We compare the deduced magnetic field strength with estimates we obtain from other methods and discuss the implications of the *Planck* observations for the general picture of molecular cloud formation and evolution.

**Key words.** ISM: general – ISM: magnetic fields – ISM: clouds – dust, extinction – submillimeter: ISM – infrared: ISM

**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

<footnote>* Corresponding author: Juan D. Soler, e-mail: jsolerpu@ias.u-psud.fr</footnote>
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 \text{ 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 \text{ cm}^{-3}$), probed by OH and CN spectral lines, the same study reports a scaling of the maximum magnetic field strength $B_{max} \propto n_H^{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 $\zeta_\psi$<fnref n="1" /> 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 $\sigma_{v_\parallel}$. In the DCF interpretation $B_⊥$ is proportional to the ratio $\sigma_{v_\parallel}/\zeta_\psi$. Application of the DCF method to subregions of the Taurus MC gives estimates of $B_⊥ \approx 10 \mu\text{G}$ in low-density regions and $\approx 25$ to $\approx 42 \mu\text{G}$ inside filamentary structures (Chapman et al. 2011). Values of $B_⊥ \approx 760 \mu\text{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 \mu\text{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<fnref n="2" /> provides an unprecedented data set from a single instrument and with

<footnote n="1">¹ We use $\zeta_\psi$ to avoid confusion with $\sigma_\psi$, which was introduced in Planck Collaboration Int. XIX (2015) as the uncertainty in the polarization angle $\psi$.</footnote>

<footnote n="2">² Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries France and Italy), with contributions from NASA (USA) and telescope reflectors provided by a collaboration between ESA and a scientific consortium led and funded by Denmark.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Magnetic field and column density measured by Planck towards the Taurus MC. The colours represent column density. The “drapery” pattern, produced using the line integral convolution method (LIC, Cabral & Leedom 1993), indicates the orientation of magnetic field lines, orthogonal to the orientation of the submillimetre polarization.</figcaption></figure>

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^{22} \text{ 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 \ge 10^{22} \text{ 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 Planck 353 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 $\sigma_Q^2$, $\sigma_U^2$, and their covariance $\sigma_{QU}$ are initially at 4'8 resolution in HEALPix format<fnref n="3" /> with a pixelization at $N_{side} = 2048$, which corresponds to an effective pixel size of 1'7. 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 ($\tau_{353}$) as a proxy for the gas column density ($N_H$). The $\tau_{353}$ 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 $\mu$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 $\tau_{353}$ map, computed initially at 5' resolution, was smoothed to 10' to match the polarization maps. The errors resulting from smoothing the product $\tau_{353}$ 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 $\tau_{353}$ to $N_H$, following Planck Collaboration XI (2014), we adopted the dust opacity found using Galactic extinction measurements of quasars,

$$ \tau_{353}/N_H = 1.2 \times 10^{-26} \text{ cm}^2. \qquad (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; Reiputh 2008). Their properties are summarized in Table 1, which includes: Galactic longitude $l$ and latitude $b$ at the centre of the field; field size $\Delta l \times \Delta b$; estimate of distance; mean and maximum total column densities from dust, $\langle N_H \rangle$ and max ($N_H$), respectively; and mean $H_2$ 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 \approx 150$ pc, namely Taurus, Ophiuchus, Lupus, Chamaeleon-Musca, and Corona Australis (CrA); (b) regions located at $d \approx 300$ pc, Aquila Rift and Perseus; and (c) regions located at $d \approx 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).

<footnote n="3">³ Górski et al. (2005), http://healpix.sf.net</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Locations and sizes of the regions selected for analysis. The background map is the gas column density, $N_H$, derived from the dust optical depth at 353 GHz (Planck Collaboration XI 2014).</figcaption></figure>

<table><caption>**Table 1.** Locations and properties of the selected regions.</caption><thead><tr><th>Region</th><th><i>l</i><br>[°]</th><th><i>b</i><br>[°]</th><th>Δ<i>l</i><br>[°]</th><th>Δ<b<br>[°]</th><th>Distance<sup>a</sup><br>[pc]</th><th>⟨N<sub>H</sub>⟩<sup>b</sup><br>[10<sup>21</sup> cm<sup>-2</sup>]</th><th>Max (N<sub>H</sub>)<sup>b</sup><br>[10<sup>21</sup> cm<sup>-2</sup>]</th><th>⟨N<sub>H2</sub>⟩<sup>c</sup><br>[10<sup>21</sup> cm<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Taurus</td><td>172.5</td><td>-14.5</td><td>15.0</td><td>15.0</td><td>140</td><td>5.4</td><td>51.9</td><td>1.6</td></tr><tr><td>Ophiuchus</td><td>354.0</td><td>17.0</td><td>13.0</td><td>13.0</td><td>140</td><td>4.4</td><td>103.3</td><td>1.1</td></tr><tr><td>Lupus</td><td>340.0</td><td>12.7</td><td>12.0</td><td>12.0</td><td>140</td><td>3.8</td><td>30.8</td><td>1.2</td></tr><tr><td>Chamaeleon-Musca</td><td>300.0</td><td>-15.0</td><td>16.0</td><td>16.0</td><td>160</td><td>2.3</td><td>29.7</td><td>1.3</td></tr><tr><td>Corona Australis (CrA)</td><td>0.0</td><td>-22.0</td><td>12.0</td><td>12.0</td><td>170</td><td>1.1</td><td>40.5</td><td>1.2</td></tr><tr><td>Aquila Rift</td><td>27.0</td><td>8.0</td><td>12.0</td><td>12.0</td><td>260</td><td>9.3</td><td>58.7</td><td>1.9</td></tr><tr><td>Perseus</td><td>159.0</td><td>-20.0</td><td>9.0</td><td>9.0</td><td>300</td><td>3.9</td><td>94.8</td><td>2.6</td></tr><tr><td>IC 5146</td><td>94.0</td><td>-5.5</td><td>5.0</td><td>5.0</td><td>400</td><td>3.7</td><td>22.6</td><td>1.0</td></tr><tr><td>Cepheus</td><td>110.0</td><td>15.0</td><td>16.0</td><td>16.0</td><td>440</td><td>4.2</td><td>21.3</td><td>1.2</td></tr><tr><td>Orion</td><td>212.0</td><td>-16.0</td><td>16.0</td><td>16.0</td><td>450</td><td>5.0</td><td>93.6</td><td>2.2</td></tr></tbody></table>

**Notes.** (a) The estimates of distances are from: Schlafly et al. (2014) for Taurus, Ophiuchus, Perseus, IC 5146, Cepheus, and Orion; Knude & Hog (1998) for Lupus and CrA; Whittet et al. (1997) for Chamaeleon-Musca; and Straizys et al. (2003) for Aquila Rift. (b) Estimated from $\tau_{353}$ using Eq. (1) for the selected pixels defined in Appendix A. (c) Using the line integral $W_{CO}$ over $-10 < v_{\parallel}/(\text{km s}^{-1}) < 10$ from the Dame et al. (2001) survey and $X_{CO} = 1.8 \times 10^{20} \text{ cm}^{-2} \text{ K}^{-1} \text{ km}^{-1} \text{ s}$.

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<fnref n="4" /> 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

<footnote n="4">⁴ In equatorial coordinates.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Left: column density map, log₁₀(NH/cm⁻²), overlaid with magnetic field pseudo-vectors whose orientations are inferred from the Planck 353 GHz polarization observations. The length of the pseudo-vectors is normalized so does not reflect the polarization fraction. In this first group, the regions analysed are, from top to bottom, Taurus, Ophiuchus, Lupus, Chamaeleon-Musca, and CrA. Right: HROs for the lowest, an intermediate, and the highest NH bin (black, blue, and red, respectively). For a given region, bins have equal numbers of selected pixels (see Sect. 4.1.1 and Appendix A) within the NH ranges labelled. The intermediate bin corresponds to selected pixels near the blue contours in the column density images. The horizontal dashed line corresponds to the average per angle bin of 15°. The widths of the shaded areas for each histogram correspond to the ±1σ uncertainties related to the histogram binning operation. Histograms peaking at 0° correspond to B⊥ predominantly aligned with iso-NH contours. Histograms peaking at 90° and/or -90° correspond to B⊥ predominantly perpendicular to iso-NH contours.</figcaption></figure>

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<fnref n="5" />

$$ \phi = \arctan(|\nabla \tau_{353} \times \hat{E}|, \nabla \tau_{353} \cdot \hat{E}), \quad (2) $$

where, as illustrated in Fig. 6, ∇τ₃₅₃ is perpendicular to the tangent of the iso-τ₃₅₃ curves, the orientation of the unit polarization

<footnote n="5">⁵ In this paper we use the version of arctan with two signed arguments to resolve the π ambiguity in the orientation of pseudo-vectors (Planck Collaboration Int. XIX 2015).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3. continued.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. Same as Fig. 3 for the second group, Aquila Rift and Perseus.</figcaption></figure>

pseudo-vector $\hat{E}$, perpendicular to $B_⊥$, is characterized by the polarization angle

$$ \psi = \frac{1}{2} \arctan(-U, Q), \qquad (3) $$

and in Eq. (2), as implemented, the norm actually carries a sign when the range used for $\phi$ is between -90° and 90°.

The uncertainties in $\phi$ due to the variance of the $\tau_{353}$ 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 $|\nabla\tau_{353}|$ 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 $|\nabla\tau_{353}|$, 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, $\phi$, is computed by using polarization measurements with a high S/N in Stokes Q and U, i.e., only considering pixels with $|Q|/\sigma_Q$ or $|U|/\sigma_U > 3$. This selection allows
<figure><img src="image_8.png" /><figcaption>Fig. 5. Same as Fig. 3 for the third group, IC 5146, Cepheus, and Orion.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Schematic of the vectors involved in the calculation of the relative orientation angle $\phi$.</figcaption></figure>

the unambiguous definition of $\hat{E}$ 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 $\phi$ 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 $\mathbf{B}_{\perp}$ with the $N_H$ contours. A convex histogram, peaking at 90° and/or -90°, corresponds to the preferred orientation of $\mathbf{B}_{\perp}$ 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 \approx 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_H$ value. 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\sigma$ uncertainties related to the histogram binning operation, which are greater than the uncertainties produced by the variances of $Q$, $U$, and $\tau_{353}$ (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 $\mathbf{B}_{\perp}$ 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 \approx 300$ pc, Aquila Rift and Perseus, are shown in the right-hand column of Fig. 4. They

<table><caption>**Table 2.** Fit of $\xi$ vs. $\log_{10}(N_H/cm^{-2})$.</caption><thead><tr><th>Region</th><th>C<sub>HRO</sub></th><th>X<sub>HRO</sub></th></tr></thead><tbody><tr><td>Taurus</td><td>-0.53</td><td>21.84</td></tr><tr><td>Ophiuchus</td><td>-0.22</td><td>22.70</td></tr><tr><td>Lupus</td><td>-0.28</td><td>21.72</td></tr><tr><td>Chamaeleon-Musca</td><td>-0.51</td><td>21.67</td></tr><tr><td>Corona Australia (CrA)</td><td>-0.11</td><td>24.14</td></tr><tr><td>Aquila Rift</td><td>-0.60</td><td>22.23</td></tr><tr><td>Perseus</td><td>-0.45</td><td>21.76</td></tr><tr><td>IC 5146</td><td>-0.68</td><td>21.79</td></tr><tr><td>Cepheus</td><td>-0.44</td><td>21.90</td></tr><tr><td>Orion</td><td>-0.28</td><td>21.88</td></tr></tbody></table>

**Notes.** See Eq. (6) and Fig. 7.

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 \approx 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 $\xi$

The changes in the HROs are quantified using the histogram shape parameter $\xi$, defined as

$$ \xi = \frac{A_c - A_e}{A_c + A_e}, \qquad (4) $$

where $A_c$ is the area in the centre of the histogram ($-22.5° < \phi < 22.5°$) and $A_e$ the area in the extremes of the histogram ($-90.0° < \phi < -67.5°$ and $67.5° < \phi < 90.0°$). The value of $\xi$, 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 $\mathbf{B}_{\perp}$ mostly aligned with $N_H$ contours would have $\xi > 0$. A convex histogram corresponding to $\mathbf{B}_{\perp}$ mostly perpendicular to $N_H$ contours would have $\xi < 0$. A flat histogram corresponding to no preferred relative orientation would have $\xi \approx 0$.

The uncertainty in $\xi$, $\sigma_{\xi}$, is obtained from

$$ \sigma_{\xi}^{2} = \frac{4 (A_{c}^{2} \sigma_{A_{c}}^{2} + A_{e}^{2} \sigma_{A_{e}}^{2})}{(A_{c} + A_{e})^{4}}. \qquad (5) $$

The variances of the areas, $\sigma_{A_e}^2$ and $\sigma_{A_c}^2$, characterize the jitter of the histograms. If the jitter is large, $\sigma_{\xi}$ is large compared to $|\xi|$ and the relative orientation is indeterminate. The jitter depends on the number of bins in the histogram, but $\xi$ does not.

Figure 7 illustrates the change in $\xi$ as a function of $\log_{10}(N_H/cm^{-2})$ of each bin. For most of the clouds, $\xi$ 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 $\xi$ from positive to negative are seen in Taurus, Chamaeleon-Musca, Aquila Rift, Perseus, and IC 5146.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Histogram shape parameter $\xi$ (Eqs. (4) and (5)) calculated for the different $N_H$ bins in each region. The cases $\xi > 0$ and $\xi < 0$ correspond to the magnetic field oriented mostly parallel and perpendicular to the structure contours, respectively. For $\xi \approx 0$ there is no preferred orientation. The displayed values of $C_{HRO}$ and $X_{HRO}$ were calculated from Eq. (6) and correspond to the grey dashed line in each plot.</figcaption></figure>

The trend in $\xi$ vs. $\log_{10}(N_H/cm^{-2})$ can be fit roughly by a linear relation

$$ \xi = C_{HRO} [\log_{10}(N_H/cm^{-2}) - X_{HRO}] \quad (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_{10}(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 $\xi$ 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 $\xi$ is seen in CrA, where $\xi$ 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_{10}(N_H/cm^{-2}) < 21.6$; see also Fig. 11 below) located directly south and directly east of the
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Same as Fig. 3 for the two test regions located directly east (ChamEast) and south of Chamaeleon-Musca (ChamSouth).</figcaption></figure>

Chamaeleon-Musca region. As in CrA, we find that $B_{\perp}$ is mostly parallel to the $N_H$ contours, a fairly constant $\xi$, and no indication of predominantly perpendicular relative orientation.

## 4.2. Comparisons with previous studies

The above trends in relative orientation between $B_{\perp}$ and the $N_H$ contours in targeted MCs, where $B_{\perp}$ tend to become perpendicular to the $N_H$ contours at high $\bar{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_{\perp}$ is mostly aligned with the lowest $N_H$ contours ($20.8 < \log_{10}(N_H/\text{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
Planck Collaboration: Probing the role of the magnetic field in the formation of structure in molecular clouds

<figure><img src="image_12.png" /></figure>

**Fig. 9.** Maps of the absolute value of the relative orientation angle, $|\phi|$, in the Taurus region. These maps are produced after smoothing the input maps to beam FWHMs of 10', 15', 30', and 60' and then resampling the grid to sample each beam FWHM with the same number of pixels. The regions in red correspond to $B_\perp$ close to perpendicular to $N_H$ structures. The regions in blue correspond to $B_\perp$ close to parallel to $N_H$ structures. The black contour, corresponding to the $N_H$ value of the intermediate contour introduced in Fig. 3, provides a visual reference to the cloud structure.

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_{10}(N_H/cm^{-2}) \ge 21.5$.

# 5. Discussion

## 5.1. The relative orientation between $B_\perp$ 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 $\nabla\tau_{353}$. 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).
<figure><img src="image_13.png" /><figcaption>Fig. 10. HROs of the Taurus region after smoothing the input maps to beam FWHMs of 15', 30', and 60', shown from left to right, respectively.</figcaption></figure>

### 5.1.2. Statistical significance of the HRO signal

Our results reveal a systemic change of $\xi$ 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 $\xi$ 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 $\epsilon = 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 isodensity 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, $\xi$, show that the simulation with super-Alfvénic turbulence does not undergo a transition in relative orientation from parallel to perpendicular for $\log_{10}(N_H/cm^{-2}) < 23$. In contrast, most of the observed clouds show a decrease in $\xi$ with increasing $N_H$, close to the trends seen from the simulations with Alfvénic or sub-Alfvénic turbulence for $\log_{10}(N_H/cm^{-2}) < 23$. Furthermore, $X_{HRO}$, the value of $\log_{10}(N_H/cm^{-2})$ where $\xi$ 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 ($\sigma_v = 2.0$ km s$^{-1}$ and $n = 500$ cm$^{-3}$) are typical of those in the selected regions ($\sigma_{v\parallel}$ is given in Table D.1), the similarities in the dependence of $\xi$ 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,
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Histogram shape parameter $\xi$ (Eqs. (4) and (5)) calculated for the different $N_H$ bins in each region. Top: relative orientation in synthetic observations of simulations with super-Alfvénic (blue), Alfvénic (green), and sub-Alfvénic (red) turbulence, as detailed in Soler et al. (2013). Middle: relative orientation in the regions selected from the Planck all-sky observations, from Fig. 7. The blue data points correspond to the lowest $N_H$ regions (CrA and the test regions in Fig. 8, ChamSouth and ChamEast) and the orange correspond to the rest of the clouds. Bottom: comparison between the trends in the synthetic observations (in colours) and the regions studied (grey). The observed smooth transition from mostly parallel ($\xi > 0$) to perpendicular ($\xi < 0$) is similar to the transition in the simulations for which the turbulence is Alfvénic or sub-Alfvéic.</figcaption></figure>

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_\perp$ 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 $\xi$ 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 sub-structures, 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_{10}(N_H/cm^{-2}) \approx 21.7$ above which the preferential orientation of $B_\perp$ 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_{\perp}$ 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_{\parallel}$ (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_{\perp}$ and $N_H$ structures and to model the potential imprint in $B_{\parallel}$ observations and in $B_{\perp}$ 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 ($\epsilon$) 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^{22} \text{ cm}^{-2}$ (reported in Planck Collaboration Int. XIX 2015), when seen at small scales, might be interpreted in terms of a decrease of $\epsilon$ 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 $\epsilon$ 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 $\epsilon$ 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_{\perp}$ 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 $\epsilon$ 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_{\perp}$), 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_{\perp}$ 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_{10}(N_H/\text{cm}^{-2}) \approx 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 \mu\text{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.

*Acknowledgements.* 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

Andersson, B.-G. 2015, in Astrophys. Space Sci. Libr. 407, eds. A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, 59

Arzoumanian, D., André, P., Peretto, N., & Könyves, V. 2013, A&A, 553, A119

Audit, E., & Hennebelle, P. 2005, A&A, 433, 1

Ballesteros-Paredes, J., Hartmann, L., & Vázquez-Semadeni, E. 1999, ApJ, 527, 285

Bally, J. 2008, Overview of the Orion Complex, ed. B. Reipurth, 4, 459

Bally, J., Walawender, J., Johnstone, D., Kirk, H., & Goodman, A. 2008, The Perseus Cloud, ed. B. Reipurth, 308

Banerjee, R., Vázquez-Semadeni, E., Hennebelle, P., & Klessen, R. S. 2009, MNRAS, 398, 1082

Benoît, A., Ade, P., Amblard, A., et al. 2004, A&A, 424, 571

Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339

Bierman, E. M., Matsumura, T., Dowell, C. D., et al. 2011, ApJ, 741, 81

Burkhart, B., Lazarian, A., Leão, I. C., de Medeiros, J. R., & Esquivel, A. 2014, ApJ, 790, 130

Cabral, B., & Leedom, L. C. 1993, in Special Interest Group on GRAPHics and Interactive Techniques Proceedings, 263

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113

Chapman, N. L., Goldsmith, P. F., Pineda, J. L., et al. 2011, ApJ, 741, 21

Comerón, F. 2008, The Lupus Clouds, ed. B. Reipurth, 5, 295

Crutcher, R. M. 2005, in Magnetic Fields in the Universe: From Laboratory and Stars to Primordial Structures, eds. E. M. de Gouveia dal Pino, G. Lugones, & A. Lazarian, AIP Conf. Ser., 784, 129

Crutcher, R. M. 2012, ARA&A, 50, 29

Crutcher, R. M., Nutter, D. J., Ward-Thompson, D., & Kirk, J. M. 2004, ApJ, 600, 279

Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792

Davis, L. 1951, Phys. Rev., 81, 890

Davis, Jr., L., & Greenstein, J. L. 1951, ApJ, 114, 206

Dotson, J. L. 1996, ApJ, 470, 566

Dotson, J. L., Davidson, J., Dowell, C. D., Schleuning, D. A., & Hildebrand, R. H. 2000, ApJS, 128, 335

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)

Eiroa, C., Djupvik, A. A., & Casali, M. M. 2008, The Serpens Molecular Cloud, ed. B. Reipurth, 693

Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211

Falceta-Gonçalves, D., Lazarian, A., & Kowal, G. 2008, ApJ, 679, 537

Franco, G. A. P., Alves, F. O., & Girart, J. M. 2010, ApJ, 723, 146

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852

Harvey, P. M., Huard, T. L., Jørgensen, J. K., et al. 2008, ApJ, 680, 495

Heiles, C., & Haverkorn, M. 2012, Space Sci. Rev., 166, 293

Heitsch, F., Zweibel, E. G., Mac Low, M.-M., Li, P., & Norman, M. L. 2001, ApJ, 561, 800

Heitsch, F., Slyz, A. D., Devriendt, J. E. G., Hartmann, L. W., & Burkert, A. 2006, ApJ, 648, 1052

Heitsch, F., Stone, J. M., & Hartmann, L. W. 2009, ApJ, 695, 248

Hennebelle, P. 2013, A&A, 556, A153

Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55

Heyer, M., Gong, H., Ostriker, E., & Brunt, C. 2008, ApJ, 680, 420

Hildebrand, R. H. 1988, QJRAS, 29, 327

Hildebrand, R. H., Kirby, L., Dotson, J. L., Houde, M., & Vaillancourt, J. E. 2009, ApJ, 696, 567

Hiltner, W. A. 1949, Science, 109, 165

Hoang, T., & Lazarian, A. 2009, ApJ, 697, 1316

Houde, M., Vaillancourt, J. E., Hildebrand, R. H., Chitsazzadeh, S., & Kirby, L. 2009, ApJ, 706, 1504

Houde, M., Rao, R., Vaillancourt, J. E., & Hildebrand, R. H. 2011, ApJ, 733, 109

Houde, M., Fletcher, A., Beck, R., et al. 2013, ApJ, 766, 49

Inoue, T., & Inutsuka, S.-i. 2009, ApJ, 704, 161

Jones, T. J., Bagley, M., Krejny, M., Andersson, B.-G., & Bastien, P. 2015, AJ, 149, 31

Kenyon, S. J., Gómez, M., & Whitney, B. A. 2008, Low Mass Star Formation in the Taurus-Auriga Clouds, ed. B. Reipurth, 4, 405

Knude, J., & Hog, E. 1998, A&A, 338, 897

Kobulinicky, H. A., Molnar, L. A., & Jones, T. J. 1994, AJ, 107, 1433

Koch, P. M., Tang, Y.-W., & Ho, P. T. P. 2013, ApJ, 775, 77

Koyama, H., & Inutsuka, S.-i. 2002, ApJ, 564, L97

Kun, M., Kiss, Z. T., & Balog, Z. 2008, Star Forming Regions in Cepheus, ed. B. Reipurth, 136

Lamarre, J., Puget, J., Ade, P. A. R., et al. 2010, A&A, 520, A9

Lazarian, A., & Hoang, T. 2007, MNRAS, 378, 910

Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700

Li, H.-B., & Houde, M. 2008, ApJ, 677, 1151

Li, H., Griffin, G. S., Krejny, M., et al. 2006, ApJ, 648, 340

Li, H.-B., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707

Luhman, K. L. 2008, Chamaeleon, ed. B. Reipurth, 5, 169

Martin, P. G., Roy, A., Bontemps, S., et al. 2012, ApJ, 751, 28

Matthaeus, W. H., Pouquet, A., Mininni, P. D., Dmitruk, P., & Breech, B. 2008, Phys. Rev. Lett., 100, 085003

Matthews, B. C., Wilson, C. D., & Fiege, J. D. 2001, ApJ, 562, 400

Matthews, B. C., McPhee, C. A., Fissel, L. M., & Curran, R. L. 2009, ApJS, 182, 143

Matthews, T. G., Ade, P. A. R., Angilè, F. E., et al. 2014, ApJ, 784, 116

Mckee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Mckee, C. F., Zweibel, E. G., Goodman, A. A., & Heiles, C. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 327

Molinari, S., Schisano, E., Faustini, F., et al. 2011, A&A, 530, A133

Moneti, A., Pipher, J. L., Helfer, H. L., McMillan, R. S., & Perry, M. L. 1984, ApJ, 282, 508

Montier, L., Plaszczynski, S., Levrier, F., et al. 2015, A&A, 574, A135

Mouschovias, T. C. 1976a, ApJ, 206, 753

Mouschovias, T. C. 1976b, ApJ, 207, 141

Mouschovias, T. C. 1991, ApJ, 373, 169

Myers, P. C., & Goodman, A. A. 1991, ApJ, 373, 509

Nakano, T. 1998, ApJ, 494, 587

Nakano, T., & Nakamura, T. 1978, PASJ, 30, 671

Neuhäuser, R., & Forbrich, J. 2008, The Corona Australis Star Forming Region, ed. B. Reipurth, 735

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38

Paradis, D., Dobashi, K., Shimoikura, T., et al. 2012, A&A, 543, A103

Pereyra, A., & Magalhães, A. M. 2004, ApJ, 603, 584

Planck Collaboration XXIV. 2011, A&A, 536, A24

Planck Collaboration I. 2014, A&A, 571, A1

Planck Collaboration II. 2014, A&A, 571, A2

Planck Collaboration V. 2014, A&A, 571, A5

Planck Collaboration VI. 2014, A&A, 571, A6

Planck Collaboration VIII. 2014, A&A, 571, A8

Planck Collaboration XI. 2014, A&A, 571, A11

Planck Collaboration Int. XIX. 2015, A&A, 576, A104

Planck Collaboration Int. XX. 2015, A&A, 576, A105

Planck Collaboration Int. XXI. 2015, A&A, 576, A106

Planck Collaboration Int. XXXII. 2016, A&A, 586, A135

Planck Collaboration Int. XXXIII. 2016, A&A, 586, A136

Poppel, W. 1997, Fund. Cosmic Phys., 18, 1

Prato, L., Rice, E. L., & Dame, T. M. 2008, Where are all the Young Stars in Aquila?, ed. B. Reipurth, 4, 18

Reipurth, B. 2008, Handbook of Star Forming Regions, Volume I: The Northern Sky, 4

Santos, F. P., Franco, G. A. P., Roman-Lopes, A., Reis, W., & Román-Zúñiga, C. G. 2014, ApJ, 783, 1

Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, ApJ, 747, 21

Schlaflay, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29

Serkowski, K. 1958, Acta Astron., 8, 135

Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128

Spitzer, Jr., L. 1968, Interscience Tracts on Physics and Astronomy, 28, 262

Sridhar, S., & Goldreich, P. 1994, ApJ, 432, 612

Straižys, V., Černis, K., & Bartašiūtė, S. 2003, A&A, 405, 585

Sugitani, K., Nakamura, F., Watanabe, M., et al. 2011, ApJ, 734, 63

Tassis, K., Dowell, C. D., Hildebrand, R. H., Kirby, L., & Vaillancourt, J. E. 2009, MNRAS, 399, 1681

Ward-Thompson, D., Kirk, J. M., Crutcher, R. M., et al. 2000, ApJ, 537, L135

Whittet, D. C. B., Prusti, T., Franco, G. A. P., et al. 1997, A&A, 327, 1194

Whittet, D. C. B., Hough, J. H., Lazarian, A., & Hoang, T. 2008, ApJ, 674, 304
Wilking, B. A., Gagné, M., & Allen, L. E. 2008, Star Formation in the ρ Ophiuchi Molecular Cloud, ed. B. Reipurth, 5, 351

Zweibel, E. G. 1996, in Polarimetry of the Interstellar Medium, ed. W. G. Roberge, & D. C. B. Whittet, ASP Conf. Ser., 97, 486

¹ APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

² African Institute for Mathematical Sciences, 6-8 Melrose Road, Muizenberg, Cape Town, South Africa

³ Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy

⁴ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France

⁵ Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK

⁶ Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Westville Campus, Private Bag X54001, Durban 4000, South Africa

⁷ CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada

⁸ CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France

⁹ California Institute of Technology, Pasadena, California, USA

¹⁰ Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001 Teruel, Spain

¹¹ Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, USA

¹² DSM/Irfu/SPP, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France

¹³ DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark

¹⁴ Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland

¹⁵ Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain

¹⁶ Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33007 Oviedo, Spain

¹⁷ Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, ON M5S Ontario, Canada

¹⁸ Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands

¹⁹ Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1 British Columbia, Canada

²⁰ Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA

²¹ Department of Physics and Astronomy, University College London, London WC1E 6BT, UK

²² Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, FL 32306, USA

²³ Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00014 Helsinki, Finland

²⁴ Department of Physics, Princeton University, Princeton, NJ 08544-0708, USA

²⁵ Department of Physics, University of California, Santa Barbara, CA 93106, USA

²⁶ Department of Physics, University of Illinois at Urbana-Champaign, 1110 West Green Street, Urbana, IL 61801, USA

²⁷ Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy

²⁸ Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy

²⁹ Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, 00185 Roma, Italy

³⁰ Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16, 20133 Milano, Italy

³¹ Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, 34128 Trieste, Italy

³² Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00133 Roma, Italy

³³ Discovery Center, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark

³⁴ Escola de Artes, Ciências e Humanidades, Universidade de São Paulo, Rua Arlindo Bettio 1000, CEP 03828-000, São Paulo, Brazil

³⁵ European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28691 Madrid, Spain

³⁶ European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

³⁷ Facoltà di Ingegneria, Università degli Studi e-Campus, via Isimbardi 10, Novedrate (CO), 22060, Italy

³⁸ Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L'Aquila, Italy

³⁹ HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany

⁴⁰ Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00100 Helsinki, Finland

⁴¹ INAF-Osservatorio Astrofisico di Catania, via S. Sofia 78, 95123 Catania, Italy

⁴² INAF-Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

⁴³ INAF-Osservatorio Astronomico di Roma, via di Frascati 33, 00078 Monte Porzio Catone, Italy

⁴⁴ INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, Trieste, Italy

⁴⁵ INAF/IASF Bologna, via Gobetti 101, 40129 Bologna, Italy

⁴⁶ INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy

⁴⁷ INFN, Sezione di Bologna, via Irnerio 46, 40126 Bologna, Italy

⁴⁸ INFN, Sezione di Roma 1, Università di Roma Sapienza, P.le Aldo Moro 2, 00185 Roma, Italy

⁴⁹ INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1 Roma, Italy

⁵⁰ INFN/National Institute for Nuclear Physics, via Valerio 2, 34127 Trieste, Italy

⁵¹ IPAG: Institut de Planétologie et d'Astrophysique de Grenoble, Université Grenoble Alpes, IPAG, 38000 Grenoble, France, CNRS, IPAG, 38000 Grenoble, France

⁵² Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK

⁵³ Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

⁵⁴ Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, 38042 Grenoble, France

⁵⁵ Institut Universitaire de France, 103 bd Saint-Michel, 75005 Paris, France

⁵⁶ Institut d'Astrophysique Spatiale, CNRS (UMR 8617) Université Paris-Sud 11, Bâtiment 121, 91400 Orsay, France

⁵⁷ Institut d'Astrophysique de Paris, CNRS (UMR 7095), 98bis boulevard Arago, 75014 Paris, France

⁵⁸ Institut of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁵⁹ Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0315 Oslo, Norway

⁶⁰ Instituto de Astrofísica de Canarias, C/vía Láctea s/n, La Laguna 38205 Tenerife, Spain

⁶¹ Instituto de Física de Cantabria (CSIC-Universidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain

⁶² Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy

⁶³ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA

⁶⁴ Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK

⁶⁵ Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
⁶⁶ Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

⁶⁷ Kazan Federal University, 18 Kremlyovskaya St., 420008 Kazan, Russia

⁶⁸ LAL, Université Paris-Sud, CNRS/IN2P3, 91405 Orsay, France

⁶⁹ LERMA, CNRS, Observatoire de Paris, 61 avenue de l'Observatoire, 75014 Paris, France

⁷⁰ Laboratoire AIM, IRFU/Service d'Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France

⁷¹ Laboratoire Traitement et Communication de l'Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault 75634 Paris Cedex 13, France

⁷² Laboratoire de Physique Subatomique et Cosmologie, Université Grenoble-Alpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France

⁷³ Laboratoire de Physique Théorique, Université Paris-Sud 11 & CNRS, Bâtiment 210, 91405 Orsay, France

⁷⁴ Lawrence Berkeley National Laboratory, Berkeley, California, USA

⁷⁵ Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., Moscow, GSP-7, 117997, Russia

⁷⁶ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany

⁷⁷ National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland

⁷⁸ Nicolaus Copernicus Astronomical Center, Bartycka 18, 00-716 Warsaw, Poland

⁷⁹ Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark

⁸⁰ Optical Science Laboratory, University College London, Gower Street, London WC1E 6BT, UK

⁸¹ SISSA, Astrophysics Sector, via Bonomea 265, 34136, Trieste, Italy

⁸² SUPA, Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

⁸³ School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK

⁸⁴ Sorbonne Université-UPMC, UMR 7095, Institut d'Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France

⁸⁵ Space Sciences Laboratory, University of California, Berkeley, CA 94720, USA

⁸⁶ Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167 Karachai-Cherkessian Republic, Russia

⁸⁷ Sub-Department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK

⁸⁸ UPMC Univ Paris 06, UMR 7095, 98bis boulevard Arago, 75014 Paris, France

⁸⁹ Université de Toulouse, UPS-OMP, IRAP, 31028 Toulouse Cedex 4, France

⁹⁰ University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain

⁹¹ University of Granada, Instituto Carlos I de Física Teórica y Computacional, 18071 Granada, Spain

⁹² Warsaw University Observatory, Aleje Ujazdowskie 4, 00-478 Warszawa, Poland
# 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.

## A.1. Gradient mask

The dust optical depth, $\tau_{353}$, observed in each region, can be interpreted as

$$ \tau_{353}^{\text{OBS}} = \tau_{353}^{\text{MC}} + \tau_{353}^{\text{BG}} + \delta_{\tau_{353}}, \quad (A.1) $$

where $\tau_{353}^{\text{MC}}$ is the optical depth of the MC, $\tau_{353}^{\text{BG}}$ is the optical depth of the diffuse regions behind and/or in front of the cloud (background/foreground), and $\delta_{\tau_{353}}$ the noise in the optical depth map with variance $\sigma_{\tau_{353}}^2 = \delta_{\tau_{353}}^2$.

The gradient of the optical depth can be then written as

$$ \nabla \tau_{353}^{\text{OBS}} = \nabla \tau_{353}^{\text{MC}} + \nabla(\tau_{353}^{\text{BG}} + \delta_{\tau_{353}}). \quad (A.2) $$

We quantify the contribution of the background/foreground and the noise, $\nabla(\tau_{353}^{\text{BG}} + \delta_{\tau_{353}})$, by evaluating $\nabla\tau_{353}$ 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, $\langle|\nabla\tau_{353}^{\text{REF}}|\rangle$, and use this value as a threshold for selecting the regions of the map where $\nabla\tau_{353}$ carries significant information about the structure of the cloud. We note that this threshold includes a contribution from the noise $\nabla\delta_{\tau_{353}}$. The HROs presented in this study correspond to regions in each field where $|\nabla\tau_{353}| > \langle|\nabla\tau_{353}^{\text{REF}}|\rangle$.

## A.2. Polarization mask

The total Stokes parameters $Q$ and $U$ measured in each region can be interpreted as

$$ Q^{\text{OBS}} = Q^{\text{MC}} + Q^{\text{BG}} + \delta_Q, \quad U^{\text{OBS}} = U^{\text{MC}} + U^{\text{BG}} + \delta_U, \quad (A.3) $$

where $Q^{\text{MC}}$ and $U^{\text{MC}}$ correspond to the polarized emission from the MC, $Q^{\text{BG}}$ and $U^{\text{BG}}$ correspond to the polarized emission from the diffuse background/foreground, and $\delta_Q$ and $\delta_U$ are the noise contributions to the observations, such that the variances $\sigma_Q^2 = \delta_Q^2$ and $\sigma_U^2 = \delta_U^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, $Q_{\text{rms}}^{\text{REF}}$ and $U_{\text{rms}}^{\text{REF}}$. The HROs presented in this study correspond to pixels in each region where $|Q| > 2Q_{\text{rms}}^{\text{REF}}$ or $|U| > 2U_{\text{rms}}^{\text{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{Q^2 + U^2} > 2\sqrt{(Q_{\text{rms}}^{\text{REF}})^2 + (U_{\text{rms}}^{\text{REF}})^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\sigma_Q$ or $|U| > 3\sigma_U$. This aims to select pixels where the uncertainty in the polarization angle is smaller than the size of the angle bins used for

<table><caption>Table A.1. Selection of data.</caption><thead><tr><th>Region</th><th>⟨N<sup>BG</sup><sub>H</sub>⟩<br>[10<sup>20</sup> cm<sup>-2</sup>]</th><th>f<sub>v</sub><br>[%]</th><th>f<sub>pol</sub><br>[%]</th><th>f<sub>tot</sub><br>[%]</th></tr></thead><tbody><tr><td>Taurus</td><td>6.2</td><td>66</td><td>78</td><td>28</td></tr><tr><td>Ophiuchus</td><td>5.6</td><td>65</td><td>82</td><td>31</td></tr><tr><td>Lupus</td><td>9.6</td><td>65</td><td>67</td><td>24</td></tr><tr><td>Chamaeleon-Musca</td><td>6.3</td><td>49</td><td>83</td><td>31</td></tr><tr><td>Corona Australia (CrA)</td><td>5.4</td><td>34</td><td>40</td><td>4</td></tr><tr><td>Aquila Rift</td><td>18.4</td><td>48</td><td>96</td><td>32</td></tr><tr><td>Perseus</td><td>5.3</td><td>60</td><td>59</td><td>16</td></tr><tr><td>IC 5146</td><td>26.4</td><td>38</td><td>93</td><td>29</td></tr><tr><td>Cepheus</td><td>7.5</td><td>66</td><td>80</td><td>36</td></tr><tr><td>Orion</td><td>6.1</td><td>63</td><td>67</td><td>24</td></tr><tr><td>ChamEast</td><td>6.3</td><td>33</td><td>38</td><td>10</td></tr><tr><td>ChamSouth</td><td>4.9</td><td>36</td><td>42</td><td>13</td></tr></tbody></table>

**Notes.** Mean column density of the background/foreground for each region $\langle N_H^{\text{BG}} \rangle$ estimated from a reference field at the same Galactic latitude; percentage $f_v$ of all pixels where $|\nabla\tau_{353}| > \langle|\nabla\tau_{353}^{\text{REF}}|\rangle$; percentage $f_{\text{pol}}$ of all pixels where $|Q| > 2Q_{\text{rms}}^{\text{REF}}$ or $|U| > 2U_{\text{rms}}^{\text{REF}}$ (first polarization criterion) and where $|Q| > 3\sigma_Q$ or $|U| > 3\sigma_U$ (second polarization criterion); and percentage $f_{\text{tot}}$ of all pixels used for the HRO analysis.

the constructions of the HRO (Serkowski 1958; Montier et al. 2015). In terms of the total polarized intensity, $P = \sqrt{Q^2 + U^2}$, and following Eqs. (B.4) and (B.5) in Planck Collaboration Int. XIX (2015), the second criterion corresponds to $P/\sigma_P > 3$ and uncertainties in the polarized orientation angle $\sigma_\psi < 10^\circ$.

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 $\langle N_H^{\text{BG}} \rangle$. 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^\circ$ (see Sect. 4.1.1).

## B.1. Calculation of $\nabla\tau_{353}$ and the uncertainty of its orientation

The optical depth gradient ($\nabla\tau_{353}$) is calculated by convolving the $\tau_{353}$ map with a Gaussian derivative kernel (Soler et al. 2013), such that $\nabla\tau_{353}$ corresponds to

$$ \nabla\tau_{353} = (G_x \otimes \tau_{353})\hat{i} + (G_y \otimes \tau_{353})\hat{j} = g_x\hat{i} + g_y\hat{j}, \quad (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-$\tau_{353}$ contour is calculated from
the components of the gradient vector,

$$ \theta = \arctan(-g_x, g_y). \qquad (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

$$ \delta_{g_x} = G_x \otimes \delta_{\tau_{353}}; \quad \delta_{g_y} = G_y \otimes \delta_{\tau_{353}}, \qquad (B.3) $$

from which we obtain the $\sigma_{g_x}^2 = \overline{\delta_{g_x}^2}$ and $\sigma_{g_y}^2 = \overline{\delta_{g_y}^2}$.

The standard deviation of the angle $\theta$ can be written as

$$ \sigma_{\theta} = \sqrt{\left(\frac{\partial \theta}{\partial g_x}\right)^2 \sigma_{g_x}^2 + \left(\frac{\partial \theta}{\partial g_y}\right)^2 \sigma_{g_y}^2}, \qquad (B.4) $$

which corresponds to

$$ \sigma_{\theta} = \frac{1}{g_x^2 + g_y^2} \sqrt{g_y^2 \sigma_{g_x}^2 + g_x^2 \sigma_{g_y}^2}. \qquad (B.5) $$

In the application discussed here, the standard deviations in the $\tau_{353}$ 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.

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

### 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 $\sigma_Q$ and $\sigma_U$ (Planck Collaboration II 2014; Planck Collaboration VI 2014; Planck Collaboration V 2014; Planck Collaboration VIII 2014). Given that $\sigma_{QU}$ is smaller than $\sigma_Q^2$ and $\sigma_U^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 $\tau_{353}$ 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 $\xi$. Together, the low noise in the maps of $\tau_{353}$ 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

$$ \sigma_k^2 = h_k \left(1 - \frac{h_k}{h_{\text{tot}}}\right), \qquad (B.6) $$

where $h_k$ is the number of samples in the kth bin, and $h_{\text{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.

<figure><img src="image_15.png" /><figcaption>Fig. B.1. HROs in the Taurus region that correspond to the indicated $N_H$ bins. The plotted values are obtained using the original $\tau_{353}$ map at 10' resolution and maps of the Stokes parameters $Q_r$ and $U_r$, which correspond to 1000 random noise realizations. Each realization is generated using a Gaussian probability density function centred on the measured values Q and U with variances $\sigma_Q^2$ and $\sigma_U^2$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. B.2. Histogram shape parameter, $\xi$, as a function of $\log_{10}(N_H/cm^{-2})$ in the Taurus region. The values are obtained using the $\tau_{353}$ map at 10' resolution and maps of the Stokes parameters $Q_r$ and $U_r$ that correspond to 1000 random-noise realizations. Each realization is generated using a Gaussian probability density function centred on the measured values Q and U with variances $\sigma_Q^2$ and $\sigma_U^2$. By joining the values at each $N_H$ bin, we find a trend very close to the black line in Fig. 7, with little dispersion from the noise; much larger are the uncertainties in evaluating $\xi$ at each $N_H$ bin, as given in Fig. 7 but not shown here.</figcaption></figure>

### B.2.2. Uncertainties in the histogram shape parameter $\xi$

As in the case of the HRO, the uncertainty in $\xi$, as defined in Eq. (4), can be quantified using the random realizations introduced in the previous section. Figure B.2 shows the dependence of $\xi$ on $\log_{10}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 $\xi$ 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 $\xi$ in Fig. B.2.
<figure><img src="image_17.png" /><figcaption>**Fig. C.1.** Histogram shape parameter, $\xi$, as a function of $\log_{10}(N_{\text{H}}/\text{cm}^{-2})$ in the Taurus region. Blue lines join the $\xi$ values obtained using the $\tau_{353}$ map at 10' resolution and each of 1000 random realizations of $Q$ and $U$ maps corresponding to a uniform distribution of $\psi$. Results in black and grey are from the analysis of the Planck data, as reported in Fig. 7.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. C.2.** Like Fig. C.1 but for values of $\xi$ (red) obtained using 1000 random realizations of $Q$ and $U$ maps corresponding to a Gaussian distribution of $\psi$ centred on $\psi_0 = 0^{\circ}$ and with standard deviation $\sigma_{\psi} = 45^{\circ}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. C.3.** Like Fig. C.1 but for values of $\xi$ (green) obtained using 1000 random realizations of $Q$ and $U$ maps corresponding to a power spectrum $P(k) \propto k^{\alpha_M}$, with $\alpha_M = -1.5, -2.5$, and $-3.5$.</figcaption></figure>

# Appendix C: Statistical significance of the HRO signal

## C.1. A product of chance?

To investigate various potential sources of the signal found in the HROs we use the $\tau_{353}$ 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 $\psi$ 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_{\text{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 $\psi$ with a Gaussian distribution and unit-length polarization pseudo-vectors. The polarization angle distribution is centred at $\psi_0 = 0^{\circ}$ with a standard deviation $\sigma_{\psi} = 45^{\circ}$. 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) \propto k^{\alpha_M}$. For the spectral indices we adopted, $\alpha_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_{\text{H}}$ found in the HROs and $\xi$ do not arise by chance from these potential sources.

## 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 $\xi$. However, they are bound to produce very different statistics compared to those in the real observations. In reality the orientations of $B_{\perp}$ 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_{\perp}$ and of the random realizations of $B_{\perp}$.

To evaluate whether the signal present in the HROs arises from an actual physical relationship between $B_{\perp}$ and $\tau_{353}$, we introduce randomness by shifting the Stokes $Q$ and $U$ maps with respect to the $\tau_{353}$ map and then calculate the corresponding HROs and $\xi$ as a function of $\log_{10}(N_{\text{H}}/\text{cm}^{-2})$. The intrinsic statistical properties of the two maps are unchanged because the two maps are unchanged, only shifted.

If the trend in $\xi$ as a function of $\log_{10}(N_{\text{H}}/\text{cm}^{-2})$ were arising randomly, these trends would be unchanged even for significant shifts. Instead the results of this experiment, illustrated in
<figure><img src="image_20.png" /><figcaption>**Fig. C.4.** Histogram shape parameter, $\xi$, as a function of $\log_{10}(N_{\text{H}}/\text{cm}^{-2})$ calculated from the $\tau_{353}$ map at 10' resolution in the Taurus region and the Stokes Q and U maps shifted in Galactic longitude and latitude by the values indicated, for $\Delta I$ and $\Delta b$, respectively.</figcaption></figure>

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.

### C.3. Projection effects

We evaluate the statistical significance of the relative orientation between $B_{\perp}$ and $\nabla\tau_{353}$ 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 $\hat{n}$ can be written as

$$u = v - (v \cdot \hat{n}) \hat{n}. \qquad (C.1)$$

<figure><img src="image_21.png" /><figcaption>**Fig. C.5.** Normalized distributions of relative orientations of random vectors in 3D (top) and of their projections in 2D (bottom). The curves correspond to mostly perpendicular vectors (red), uniform distribution of relative orientations (black), and mostly parallel vectors (green). The solid, dashed, and dotted lines correspond to dispersions of the relative orientation angle $\zeta_{\alpha} = 18^{\circ}, 36^{\circ},$ and $72^{\circ}$, respectively.</figcaption></figure>

Consider two unit vectors in 3D, $v_1$ and $v_2$, separated by an angle $\alpha$, such that

$$\cos \alpha = v_1 \cdot v_2. \qquad (C.2)$$

The angle $\beta$ between the projections of these two vectors onto a plane normal to $\hat{n}$, $u_1$, and $u_2$, can be written as

$$\begin{aligned} \cos \beta &= \frac{u_1 \cdot u_2}{|u_1||u_2|} \\ &= \frac{(v_1 \cdot v_2 - (v_1 \cdot \hat{n})(v_2 \cdot \hat{n}))}{(v_1 \cdot v_1 - (v_1 \cdot \hat{n})^2)^{1/2} (v_2 \cdot v_2 - (v_2 \cdot \hat{n})^2)^{1/2}}. \end{aligned} \qquad (C.3)$$

Given a particular distribution of angles between the vectors $v_1$ and $v_2$, this expression, which is solved numerically, is useful for evaluating the resulting distribution of angles between the projected vectors $u_1$ and $u_2$. Without any loss of generality, we can assume that $v_1$ is oriented along the axis of a spherical coordinate system, such that $v_1 = \hat{k}$, and that $v_2$ is oriented at polar angle $\theta$ and azimuth $\phi$, such that $v_2 = \cos \phi \sin \theta \hat{i} + \sin \phi \sin \theta \hat{j} + \cos \theta \hat{k}$. Then, we can simulate a distribution of $\cos \alpha$ by simulating a distribution of $\cos \theta$ and generating $\phi$ with a uniform distribution.

We thus generate a set of vectors $v_2$ that follow a particular distribution of $\cos \alpha$. We choose three examples: a uniform distribution of relative orientation between $v_2$ and $v_1$, $v_2$ vectors that are mostly parallel to $v_1$, and $v_2$ vectors that are mostly perpendicular to $v_1$. The last two are Gaussian distributions centred at $\cos \alpha = 0$ (for mostly parallel) or $\cos \alpha = \pm 1$ (for mostly perpendicular, given the periodicity of $\alpha$) with a dispersion $\zeta_{\alpha}$. These distributions are shown in the top panel of Fig. C.5 for three values of $\zeta_{\alpha}$.
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.

### D.1. Davis-Chandrasekhar-Fermi method

The DCF method estimates the strength of $B_⊥$ in a region using the dispersion of the polarization angles $\zeta_ψ$. Assuming that the magnetic field is frozen into the gas and that the dispersion of the $B_⊥$ orientation angles (or equivalently $\zeta_ψ$) is due to transverse incompressible Alfvén waves, then

$$B_{\perp}^{\text{DCF}} = \sqrt{4\pi\rho} \frac{\sigma_{v_{\parallel}}}{\zeta_{\psi}}, \qquad (\text{D.1})$$

where $\sigma_{v_\parallel}$ is the dispersion of the radial velocity of the gas (Appendix D.3) and $\rho$ the gas mass density.

We calculate $\zeta_ψ$ directly from Stokes $Q$ and $U$ using

$$\zeta_{\psi} = \sqrt{\langle (\Delta\psi)^2 \rangle} \qquad (\text{D.2})$$

and

$$\Delta\psi = \frac{1}{2} \arctan (Q \langle U \rangle - \langle Q \rangle U, Q \langle Q \rangle + \langle U \rangle U), \qquad (\text{D.3})$$

where $\langle \dots \rangle$ denotes an average over the selected pixels in each map (Planck Collaboration Int. XIX 2015).

### 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, $\Delta\psi(\ell) = \psi(\mathbf{x}) - \psi(\mathbf{x} + \ell)$, between pairs of $B_⊥$ vectors separated by displacements $\ell$ in the plane of the sky. Assuming that the angle differences are statistically isotropic (i.e., they depend only on $\ell = |\ell|$ and not on the orientation of $\ell$), they can be binned by

distance, $\ell$. From the $N(\ell)$ pairs of $B_⊥$ vectors for that bin, the square of the second-order structure function is

$$S_2^2(\ell) = \langle [\Delta\psi(\mathbf{x}, \ell)]^2 \rangle_x = \left\langle \frac{1}{N(\ell)} \sum_{i=1}^{N(\ell)} (\Delta\psi_{x,i})^2 \right\rangle_x, \qquad (\text{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

$$\Delta\psi_{x,i} = \frac{1}{2} \arctan (Q_i U_x - Q_x U_i, Q_i Q_x + U_i U_x), \qquad (\text{D.5})$$

where the subscripts $x$ and $i$ represent the central and displaced positions, respectively.

Hildebrand et al. (2009) assume that $B(\mathbf{x})$ is composed of a large-scale structured field, $B_0(\mathbf{x})$, and a random component, $B_r(\mathbf{x})$, which are statistically independent of each other. Assuming that $B_0(\mathbf{x})$ is a smoothly varying quantity, its contribution to $S_2^2(\ell)$ should increase in proportion to $\ell^2$ for distances that are much smaller than the scales at which $B_0$ 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 $\ell$, Hildebrand et al. (2009) derived the approximation

$$S_2^2(\ell) = b^2 + m^2 \ell^2, \qquad (\text{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^2$ 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{\langle B_{r,\perp}^2 \rangle^{1/2}}{B_{0,\perp}} = \frac{b}{\sqrt{2-b^2}}, \qquad (\text{D.7})$$

where $\langle B_{r,\perp}^2 \rangle^{1/2}$ stands for the root mean square (rms) variations about the large-scale magnetic field, $B_{0,\perp}$. 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{\langle B_{r,\perp}^2 \rangle^{1/2}}{B_{0,\perp}} = \frac{\sigma_{v_{\parallel}}}{V_{A,\perp}}, \qquad (\text{D.8})$$

where

$$V_{A,\perp} = \frac{B_{0,\perp}}{\sqrt{4\pi\rho}}, \qquad (\text{D.9})$$

is the speed of the transverse incompressible Alfvén waves. Combining Eqs. (D.7)–(D.9) results in

$$B_{\perp}^{\text{DCF+SF}} \equiv B_{0,\perp}^{\text{DCF+SF}} = \sqrt{4\pi\rho} \frac{\sigma_{v_{\parallel}} \sqrt{2-b^2}}{b}. \qquad (\text{D.10})$$
<figure><img src="image_22.png" /><figcaption>**Fig. D.1.** Line-of-sight centroid velocity $v_\parallel$ maps inferred from CO emission (Dame et al. 2001). Areas shown have sufficient S/N in the polarization observations (second polarization criterion in Appendix A.2) for the DCF and DCF+SF analyses and lie in the range $-10 < v_\parallel/(\text{km s}^{-1}) < 10$.</figcaption></figure>
<table><caption>**Table D.1.** Magnetic properties of the selected regions.</caption><thead><tr><th>Region</th><th>σ<sub>v||</sub><br>[km s<sup>-1</sup>]</th><th>ζ<sub>ψ</sub><br>[deg]</th><th>b<br>[deg]</th><th>B<sub>⊥</sub><sup>DCF</sup><br>[μG]</th><th>B<sub>⊥</sub><sup>DCF+SF</sup><br>[μG]</th><th>λ<sub>obs</sub><sup>DCF</sup></th><th>λ<sub>obs</sub><sup>DCF+SF</sup></th></tr></thead><tbody><tr><td>Taurus</td><td>1.2 ± 0.5</td><td>43 ± 0.1</td><td>23 ± 0.05</td><td>0.39 ± 0.01</td><td>13 ± 5</td><td>32 ± 13</td><td>0.4 ± 0.4</td><td>0.2 ± 0.1</td></tr><tr><td>Ophiuchus</td><td>0.9 ± 0.4</td><td>29 ± 0.3</td><td>20 ± 0.04</td><td>0.36 ± 0.01</td><td>13 ± 6</td><td>25 ± 11</td><td>0.4 ± 0.4</td><td>0.2 ± 0.2</td></tr><tr><td>Lupus</td><td>1.5 ± 0.6</td><td>46 ± 0.7</td><td>30 ± 0.06</td><td>0.52 ± 0.01</td><td>14 ± 5</td><td>29 ± 11</td><td>0.3 ± 0.2</td><td>0.2 ± 0.1</td></tr><tr><td>Chamaeleon-Musca</td><td>1.0 ± 0.4</td><td>36 ± 0.3</td><td>23 ± 0.05</td><td>0.40 ± 0.01</td><td>12 ± 5</td><td>27 ± 11</td><td>0.4 ± 0.3</td><td>0.2 ± 0.2</td></tr><tr><td>Corona Australis (CrA)</td><td>0.6 ± 0.2</td><td>59 ± 0.1</td><td>30 ± 0.07</td><td>0.52 ± 0.01</td><td>5 ± 2</td><td>12 ± 5</td><td>0.9 ± 0.9</td><td>0.3 ± 0.3</td></tr><tr><td>Aquila Rift</td><td>1.9 ± 0.6</td><td>43 ± 0.5</td><td>23 ± 0.09</td><td>0.40 ± 0.01</td><td>20 ± 6</td><td>50 ± 15</td><td>0.3 ± 0.2</td><td>0.1 ± 0.1</td></tr><tr><td>Perseus</td><td>1.5 ± 0.6</td><td>38 ± 0.3</td><td>29 ± 0.11</td><td>0.50 ± 0.01</td><td>17 ± 7</td><td>30 ± 11</td><td>0.3 ± 0.3</td><td>0.2 ± 0.2</td></tr><tr><td>IC 5146</td><td>1.7 ± 0.6</td><td>69 ± 0.1</td><td>49 ± 0.11</td><td>0.85 ± 0.01</td><td>11 ± 4</td><td>18 ± 6</td><td>0.5 ± 0.3</td><td>0.3 ± 0.2</td></tr><tr><td>Cepheus</td><td>1.6 ± 0.6</td><td>43 ± 0.2</td><td>20 ± 0.04</td><td>0.35 ± 0.01</td><td>16 ± 6</td><td>47 ± 18</td><td>0.3 ± 0.1</td><td>0.1 ± 0.0</td></tr><tr><td>Orion</td><td>1.7 ± 0.6</td><td>36 ± 0.1</td><td>26 ± 0.06</td><td>0.45 ± 0.01</td><td>20 ± 7</td><td>38 ± 14</td><td>0.3 ± 0.3</td><td>0.2 ± 0.2</td></tr></tbody></table>

**Notes.** The tabulated values are: the velocity dispersion, $\sigma_{v||}$; the dispersion of the polarization orientation angle, $\zeta_{\psi}$; the turbulent contribution to the angular dispersion, $b$, calculated from a linear fit to Eq. (D.6); the magnetic field strengths, $B_{\perp}^{\text{DCF}}$ and $B_{\perp}^{\text{DCF+SF}}$, calculated with the DCF method, of Eq. (D.1), and the DCF+SF method, of Eq. (D.10), with $b$ values obtained from a fit to Eq. (D.6) in the range $50' \le l \le 200'$; and their corresponding observed mass-to-flux ratios, $\lambda_{\text{obs}}^{\text{DCF}}$ and $\lambda_{\text{obs}}^{\text{DCF+SF}}$. The reported uncertainties are from the appropriate propagation of errors.

### D.3. Calculation

The estimates of velocity dispersion $\sigma_{v||}$ and mass density $\rho$ are the same in the two methods (DCF and DCF+SF).

We obtain $\sigma_{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_{||}/(\text{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 $\sigma_{v||}$. This is tabulated along with other properties in Table D.1.

The mass density, $\rho = \mu n m_p$, is the product of the proton mass, $m_p$, the mean molecular weight per hydrogen molecule, $\mu = 2.8$, and the mean number density, $n$. We require a value of $n$, which for the discussions here we approximate to be $100 \text{ 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 $\zeta_{\psi}$ corresponds to the centred second moment of the $\psi$ distribution evaluated in pixels within the selected region. Following the analysis presented in Planck Collaboration Int. XIX (2015) for $S(\ell)$, the variance on $\zeta_{\psi}$ can be expressed as

$$ \sigma_{\zeta_{\psi}}^2 = \frac{1}{N^2 \zeta_{\psi}^2} \left[ \left( \sum_{i=1}^{N} \zeta_i \right)^2 \sigma_{\psi}^2 + \sum_{i=1}^{N} \zeta_i^2 \sigma_{\psi_i}^2 \right], \quad (\text{D.11}) $$

where $\psi_i$ and $\sigma_{\psi_i}$ are the orientation angle and uncertainty for each polarization pseudo-vector, $N$ is the number of $\psi$ measurements, $\zeta_i = \psi_i - \langle\psi\rangle$, and $\langle\psi\rangle$ is the average orientation angle in

each region. Like other quadratic functions, $\sigma_{\zeta_{\psi}}^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 $\sigma_{v||}$ are considerably larger, the bias correction does not have a significant effect on our estimates. We apply Eq. (D.1) to estimate $B_{\perp}^{\text{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 $S_2^2(\ell)$ for each of the regions. The resolution of the data used is 10'. We evaluate $S_2^2(\ell)$ in steps of 3.44 for lags $0' < \ell < 34.4$ and in steps of 34.4 for lags $40' < \ell < 200'$. For each considered lag $\ell$, $\Delta\psi_{x,i}$ is evaluated pixel by pixel, considering all the pixels located in an annulus with radius $\ell$. The term $\Delta\psi_{x,i}$ is only considered in the calculation of $S_2^2(\ell)$ 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 $\ell$ above the resolution of the data and between 50' and 200', where the behaviour of $S_2^2(\ell)$ is approximately linear in $\ell^2$. We then make a linear fit of $S_2^2(\ell)$ using $\ell^2$ as the independent variable and estimate the value of $b$ and its uncertainty, $\sigma_b$. The results are plotted in Fig. D.2.

Using the calculated $b$ values$^6$, we apply Eq. (D.10) to estimate $B_{\perp}^{\text{DCF+SF}}$, and the results are listed in Table D.1. As expected, the DCF+SF method produces estimates of $B_{\perp}$ 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 $\ell$-range and so propagates into a different value of $B_{\perp}^{\text{DCF+SF}}$.

### D.4. Mass-to-flux ratios

The critical value for the mass that can be supported against gravity by a magnetic flux $\Phi$ can be estimated to first order for a uniform disk from $(M/\Phi)_{\text{crit}} \equiv 1/(2\pi G^{1/2})$

<footnote>$^6$ For Eq. (D.10) $b$ is required to be in radians, and so those are also listed in Table D.1.</footnote>
<figure><img src="image_23.png" /><figcaption>Fig. D.2. Structure function $S_2(l)$ calculated from the Q and U maps of the selected regions following Eq. (D.4). The black lines indicate the fits given by Eq. (D.6). The vertical dashed line marks the common 10' resolution of the data used in the analysis; the data are correlated for low values of $l$ causing the drop in $S_2(l)$.</figcaption></figure>

(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

$$ \lambda = \frac{(M/\Phi)}{(M/\Phi)_{\text{crit}}} = 7.6 \times 10^{-21} \frac{(N_{\text{H}_2}^{\text{mft}}/\text{cm}^{-2})}{(B^{\text{mft}}/\mu\text{G})}, \quad (\text{D.12}) $$

where $N_{\text{H}_2}^{\text{mft}}$ and $B^{\text{mft}}$ are the H$_2$ 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 $\lambda > 1$; otherwise, when $\lambda < 1$, the cloud is sub-critical, magnetically supported against gravitational collapse.

What is observable is $\lambda_{\text{obs}}$, in which $N_{\text{H}_2}^{\text{mft}}/B^{\text{mft}}$ is replaced by $N_{\text{H}_2}/B_{\perp}$. We evaluate $\lambda_{\text{obs}}$ by combining the value of $\langle N_{\text{H}_2} \rangle$ computed from the integrated CO line emission and the conversion factor $X_{\text{CO}} = (1.8 \pm 0.3) \times 10^{20} \text{ cm}^{-2} \text{ K}^{-1} \text{ km}^{-1} \text{ s}$ (Dame et al. 2001) and $B_{\perp}$ estimated with the DCF and DCF+SF methods. The $X_{\text{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_{\perp}$.
The calculated values of $\lambda_{\text{obs}}$ are listed in Table D.1. They are consistent with being less than unity.

From this we might obtain $\lambda$ by judicious deprojection, since

$$ \lambda = (N_{H_2}^{\text{mft}}/N_{H_2})(B_{\perp}/B^{\text{mft}}) \times \lambda_{\text{obs}} \equiv f_{\text{dp}} \lambda_{\text{obs}}. \quad (\text{D.13}) $$

Here, $B_{\perp}$ is always less than $B^{\text{mft}}$, pushing $f_{\text{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, $N_{H_2}^{\text{mft}} < N_{H_2}$, again lowering $f_{\text{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_{\text{dp}} = 1/3$, yielding

$$ \overline{M/\Phi} = \int_{0}^{\pi/2} \frac{M \cos \beta}{\Phi/\sin \beta} \sin \beta \, d\beta = \frac{1}{3}(M/\Phi)_{\text{obs}}, \quad (\text{D.14}) $$

where $\beta$ is the inclination angle with respect to the line of sight, and the $\sin \beta$ dependence in the flux comes from $B = B_{\perp} \sin \beta$, (Crutcher et al. 2004). For a sphere there is no $\cos \beta$ dependence and $f_{\text{dp}} = 1/2$. For a prolate spheroid elongated along the orientation of the field, the mass is multiplied by $\sin \beta$ instead of $\cos \beta$, resulting in $f_{\text{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.

## 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 $B_{\perp}^{\text{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 $\zeta_{\psi}$ and to low values of $B_{\perp}$. 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 $B_{\perp}^{\text{DCF+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_0$, is to cause the square of the second-order structure function, $S_2^2(\ell)$, to increase as $\ell^2$. This corresponds to a very specific correlation function for $B_0$. 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_{\perp}$. 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 ($\delta$) 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 $\delta$ and the corrections would have to rely on rough estimates of this value and the depth of integration ($\Delta$). Following Eq. (29) in Houde et al. (2009), rough estimates $\delta \approx 0.2$ pc and $\Delta \approx 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 $B_{\perp}^{\text{DCF+SF}}$ close to those of $B_{\perp}^{\text{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 $C \approx 0.5$ provides a good approximation to the actual magnetic field strength in cases where $\zeta_{\psi} < 25^{\circ}$. However, the effect of nonlinear amplitudes is uncertain (Zweibel 1996), and the method fails for values of $\zeta_{\psi} > 25^{\circ}$, 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 $\zeta_{\psi}$, based on a fit to the $B_{\perp}^{\text{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 $\sigma_{v_{\parallel}}$ and $\zeta_{\psi}$ 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 $\sigma_{v_{\parallel}}$ is based on the emission-line profile $v_{\parallel}$ 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_{\perp}$ are projected and integrated along the line of sight as pseudo-vectors.

Finally, in the estimates of $B_{\perp}^{\text{DCF}}$ and $B_{\perp}^{\text{DCF+SF}}$, a common mean density $\rho$ has been adopted, while in practice $\rho$ is different from cloud to cloud. Direct estimation of the values of $\rho$ from the measured column densities $\langle N_{\text{H}} \rangle$ and $\langle N_{\text{H}_2} \rangle$ 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_{\perp}$.

Given the line-of-sight integration and the fact that $\zeta_{\psi} > 25^{\circ}$ uniformly, the calculated values could be considered as upper limits to the actual $B_{\perp}$. 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.