The link between E-B polarization modes and gas column density from interstellar dust emission

Context. The analysis of the Planck polarization E and B mode power spectra of interstellar dust emission at 353 GHz recently raised new questions on the impact of Galactic foregrounds to the detection of the polarization of the Cosmic Microwave Background (CMB) and on the physical properties of the interstellar medium (ISM). In the diffuse ISM at high-latitude a clear E − B asymmetry is observed, with twice as much power in E modes than in B modes; as well as a positive correlation between the total power, T , and both E and B modes, presently interpreted in terms of the link between the structure of interstellar matter and that of the Galactic magnetic field. Aims. In this paper we aim at extending the Planck analysis of the high-latitude sky to low Galactic latitude, investigating the correlation between the T -E-B autoand cross-correlation power spectra with the gas column density from the diffuse ISM to molecular clouds. Methods. We divide the sky between Galactic latitude |b| > 5◦ and |b| < 60◦ in 552 circular patches, with an area of ∼400◦2, and we study the cross-correlations between the T -E-B power spectra and the column density of each patch using the latest release of the Planck polarization data. Results. We find that the B-to-E power ratio (DBB ` /DEE ` ) and the T E correlation ratio (rT E) depend on column density. While the former increases going from the diffuse ISM to molecular clouds in the Gould Belt, the latter decreases. This systematic variation must be related to actual changes in ISM properties. The data show significant scatter about this mean trend. The variations of DBB ` /DEE ` and rT E are observed to be anti-correlated for all column densities. In the diffuse ISM, the variance of these two ratios is consistent with a stochastic non-Gaussian model in which the values ofDBB ` /DEE ` and rT E are fixed. We finally discuss the dependencies of T B and EB with column density, which are however hampered by instrumental noise. Conclusions. For the first time, this work shows significant variations of the T -E-B power spectra of dust polarized emission across a large fraction of the Galaxy. Their dependence on multipole and gas column density is key for accurate forecasts of next generation CMB experiments and for constraining present models of ISM physics (i.e., dust properties and interstellar turbulence) that are considered responsible for the observed T -E-B signals.


Introduction
The Galactic polarized light emitted by interstellar dust grains is considered a major foreground for detecting primordial B modes of the cosmic microwave background (CMB; Planck Collaboration Int. XXX 2016, hereafter P16XXX). The E-B mode decomposition was introduced by Zaldarriaga & Seljak (1997) to characterize the polarization of the CMB as it allows for building an orthogonal base for linear polarization that is invariant under rotation, in contrast to the Stokes parameters Q and U, and separates the CMB polarization in components of different physical origins. More generally, as E and B modes are scalar (parity-even) and pseudo-scalar (parity-odd) quantities, respectively, their auto-and cross-correlation power spectra Based on observations obtained with Planck (http://www.esa. int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. are ideal to probe the two-point statistics in polarization across the sky.
In the case of CMB, the B mode power would partly be the result of tensor perturbations in the early Universe generated by primordial gravitational waves during the epoch of cosmic inflation (Kamionkowski et al. 1997). Such a detection would represent an indirect proof of the paradigm of cosmic inflation after the Big Bang. Until now a wealth of experiments from the ground, balloons, and satellites have reached the required sensitivity to perform accurate measurements of the CMB anisotropies both in intensity and in polarization (e.g., DASI (Carlstrom & DASI Collaboration 2000), ACT (Marriage & Atacama Cosmology Telescope Team 2009), POLARBEAR (Kermish et al. 2012), BICEP1/2 (Pryke & BICEP2 and Keck-Array Collaborations 2013), BOOMERanG (de Bernardis et al. 2000), SPIDER (Fraisse et al. 2013), WMAP (Bennett et al. 2013), Planck (Planck Collaboration I 2016. However, the extraction of the cosmological signal is still limited by the ability of controlling instrumental systematics and subtracting foreground contamination that add to the primordial radiation. Above 100 GHz the most important CMB foreground is interstellar dust emission. Thanks to the first full-sky maps in polarization at 353 GHz obtained with the Planck satellite (Planck Collaboration Int. XIX 2015), it has been possible to quantify the levels of E and B modes from Galactic dust. Focusing on the high and intermediate Galactic latitude sky (|b| > 35 • ) P16XXX, first, and more recently Planck Collaboration XI (2019, hereafter P18XI), showed that on average (i) the dusty Milky Way produces twice as much power in E modes than in the B modes (also referred to as E-B asymmetry); (ii) a positive correlation exists over a wide range of angular scales (for multipoles > 5) between E modes and the total intensity, Stokes I, alternatively referred to as T ; (iii) a hint of a positive correlation at large angular scales (for multipoles < 100) between T and B modes is present as well.
The origin of these observational results has yet to be established. More work is therefore needed to model these results as CMB foregrounds. They are the consequence of the physical processes in the interstellar medium (ISM) that generate and affect dust polarization. Dust grains aligned with the interstellar magnetic field (i.e., Chandrasekhar & Fermi 1953;Davis & Greenstein 1951;Lazarian & Hoang 2007;Hoang & Lazarian 2016;Hoang et al. 2018) and mixed with interstellar gas emit thermal radiation with a polarization vector preferentially perpendicular to the local orientation of the magnetic field. Hence, dust polarization observations are a suitable probe of the physical coupling between the gas dynamics and the magnetic-field structure, giving insight into magnetohydrodynamical (MHD) turbulence in the ISM (e.g., Brandenburg & Lazarian 2013).
The possibility that the cross-correlations between dust polarization power spectra are related to MHD turbulence in the ISM has been recently investigated by several authors, although no general agreement has been achieved yet. Kritsuk et al. (2018) and Kandel et al. (2017Kandel et al. ( , 2018 suggested that sub-Alfvénic turbulence at high Galactic latitude (with Alfvén Mach number M A < 0.5) may reproduce the E-B asymmetry and the positive T -E correlation at > 50. Caldwell et al. (2017) on the contrary concluded that only a narrow range of theoretical parameters in MHD simulations would account for the observations, suggesting that Planck results may likely connect to the large-scale driving of ISM turbulence. The E-B asymmetry was also found to be produced by inhomogeneous helical turbulence in Brandenburg et al. (2019), investigating the role of magnetic helicity in the emergence of parity-odd or -even quantities in interstellar polarized emission. The variety and complexity of simulated scenarios that are able to reproduce the E-B decomposition from Planck is described as well in Kim et al. (2019). The authors presented a first statistical analysis of all-sky synthetic maps of dust polarization at 353 GHz produced with the TIGRESS MHD simulations. Displacing the view point within a kpc-scale shearing box, they found large fluctuations of E-B asymmetry and T -E correlation depending both on the observer's position and on temporal fluctuations of ISM properties due to bursts of star formation. The observer's environment, and the role of the largescale Galactic magnetic field in the solar neighborhood, were also considered in Bracco et al. (2019) as a possible explanation for the positive T -E and T -B correlations at very low multipoles ( < 50) via a left-handed helical component.
For multipoles > 50, sub/trans-Alfvénic turbulence in the diffuse ISM was independently suggested by additional observational evidence. Sub/trans-Alfvénic turbulence would explain the overall alignment of the magnetic-field morphology with the distribution of filamentary matter-density structures observed with dust emission at high Galactic latitude (Planck Collaboration Int. XXXII 2016;Planck Collaboration Int. XXXVIII 2016;Soler & Hennebelle 2017). The alignment between density structures and magnetic fields, as suggested by Zaldarriaga (2001, hereafter Z01), would generate more E mode power compared to the B modes and naturally explain the positive correlation between T and E, at least on angular scales typical of interstellar filaments (for multipoles > 50).
An analysis of the histograms of relative orientation (HROs) between magnetic-field and density structures showed a change in trend from the diffuse ISM to dense molecular clouds in the Galaxy, where the magnetic field appears to be mostly perpendicular to the densest matter structures (Planck Collaboration Int. XXXV 2016). Such perpendicular configuration would produce a negative T -E correlation (see Z01). Going from the diffuse ISM to dense molecular clouds there would be a transition producing more random orientations between the magnetic-field and the density structures, reducing the E-B asymmetry. Thus, if the interpretation of the dust polarization power spectra in terms of the correlation between magnetic fields and filamentary density structures is right, we expect a density dependence of the E-B mode decomposition as well.
In this paper we present an observational work, in which we extend the Planck analysis reported in P16XXX to low Galactic latitude to investigate the dependence between the gas column density derived from the Planck dust emission data and the E and B mode power of dust polarization at 353 GHz. The paper is organized as follows: in Sect. 2.1 we describe the Planck data used in the analysis; Sect. 3 presents the E and B decomposition and the power spectra at intermediate and low Galactic latitudes; in Sect. 4 we show the correlation between the dust polarization power spectra and the gas column density; in Sect. 5 we discuss of our results. A summary is presented in Sect. 6. Two appendices (Appendices A and B) clarify our data analysis.

Data description
In this section we provide a description of the Planck polarization data and the column density map, and we describe how we divide the intermediate and low Galactic latitude sky to define the regions of interest for this analysis.

Planck polarization data
We used publicly available Planck PR3 data 1 at 353 GHz (Planck Collaboration III 2019) in HEALPix 2 format. These maps are produced from polarization sensitive bolometers and expressed in thermodynamic temperature units (K CMB , Planck Collaboration III 2019). We also used subsets of the Planck polarization data at 353 GHz, namely, the half-mission maps (HM1 and HM2), to debias the effect of instrumental noise in the autocorrelation power spectra. We used the raw Stokes IQU maps at 353 GHz at their nominal beam resolution of 4.82 full width half maximum (FWHM).

Column density map
We considered the total gas column-density map, N H , derived from the dust optical depth at 353 GHz, τ 353 . The τ 353 map (Planck Collaboration XI 2014) was obtained from the all-sky Planck intensity observations at 353, 545, and 857 GHz, and the IRAS observations at 100 µm, which were fitted using a modified blackbody spectrum. We used the τ 353 map at its nominal resolution of 5 . To scale from τ 353 to N H we adopted the same convention as in Planck Collaboration Int. XXXV (2016), (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 (Martin et al. 2012;Planck Collaboration XI 2014).
In this work, similar to what was done in Planck Collaboration Int. XXXV (2016), we want to analyze the column density of local molecular clouds around the Sun. Thus, to focus on these dense clouds, and to reduce the contribution to the total N H coming from the large-scale Galactic density gra- where N H is the column-density map smoothed to a FWHM of 12 • . The choice of this scale for the background column density is clarified in Sect. 2.3.
As shown in Fig. 1, the densest regions in δN H correspond to well-known molecular clouds in the Gould Belt: Taurus, Perseus, and California in the extreme east (labeled as A); Cepheus and Polaris in the northeast (labeled as B); Ophiuchus above the Galactic center (labeled as C); Musca and Chamaeleon in the southwest (labeled as D); and Orion in the extreme west (labeled as E).

Selected sky regions
To study the variations of the E-B mode power spectra across the sky, we divided it at intermediate and low Galactic latitudes (|b| < 60 • ) in circular patches of 12 • radius (with an area of 400 deg 2 , or a sky fraction of f sky ∼ 1%, see Appendix B) using a HEALPix grid at N side = 8 to get the central pixel of each patch. We chose this radius to be consistent with the analysis presented in P16XXX. It also explains our choice of filtering N H (see Sect. 2.2). To avoid strong depolarization caused by long lines of sight across the Galaxy, we masked the thin Galactic disk for |b| < 5 • (Planck Collaboration Int. XIX 2015). Hence, we generated a sample of 552 sky patches (see Fig. 2), within which we estimated average gas column density and dust polarization power spectra. For each circular patch, the column density value that we considered is represented by the parameter δN H = δN H (>95%) , where the brackets refer to the average over the 5% densest pixels within each patch. This choice allowed us to keep a high dynamic range in column density among the different patches. Results do not significantly change considering 10% instead of 5%. A17, page 3 of 10 A&A 632, A17 (2019)

E-B mode decomposition: methods
This section describes the formalism used to build the E and B mode power spectra from the observed Stokes Q and U parameters. We also show their values within the 552 sky patches introduced in Sect. 2.3.

E-B mode formalism
Computing angular power spectra of Stokes parameters requires some discussion. Stokes I is a scalar quantity that is invariant under rotation. The Stokes Q and U are not. Following Zaldarriaga & Seljak (1997) they transform as where n is the position in the sky and β is the rotation of the plane-of-the-sky reference (e 1 , e 2 ) in e 1 = cos β e 1 + sin β e 2 and e 2 = − sin β e 1 + cos β e 2 . We note that in the following Stokes I is alternatively referred to as T (n) for consistency with previous works. The authors of the aforementioned paper expanded these quantities in the appropriate spin-weighted basis (spherical harmonics) as and used the spin-raising (lowering) operators, ð + ( ð − ), to get two rotationally invariant quantities From Eq. (4), the expansion coefficients are written as which can be linearly combined into The E and B modes, i.e., the scalar and pseudo-scalar fields, respectively, are defined as These two quantities are rotationally invariant and they differ for parity symmetries (i.e., changing the sign of the x-axis only). Since Q (n ) = Q(n) and U (n ) = −U(n), from Eqs. (5) and (6), we can show that E (n ) = E(n) while B (n ) = −B(n). Thereby, E and B modes are even and odd quantities, respectively, under parity transformations.
The usual statistical description of the three scalar and pseudo-scalar quantities defined above (T, E, and B) is based on their power spectra as a function of the multipole , where X and Y may refer to T , E, or B. Power spectra are named auto-power spectra when X = Y and cross-power spectra when X Y. Alternatively we can use the quantity In this work we also used the normalized parameter, r XY , to quantify the correlation among the power spectra and already shown in P18XI. It is defined as follows: so that in case of perfect positive (negative) correlation r XY = 1 (−1), and in case of absence of correlation r XY = 0.

Power-spectra analysis
We computed the T -E-B power spectra in Eq. (8) for each circular sky patch using the XPOL 3 code, which is the generalization to polarization of XSPECT (Tristram et al. 2005). The XSPECT code corrects for incomplete sky coverage, pixel and beam window functions. In order not to correlate noise in the autocorrelated power spectra (i.e., X = Y) we always cross-correlate the HM1 and HM2 independent subsets of the data. We binned the power spectra in five principal multipole-bins centered in =35 (hereafter, bin 0), 80 (bin 1), 200 (bin 2), 550 (bin 3), and 1150 (bin 4), respectively. The corresponding widths are 15, 40, 200, 500, and 1200 from bin 0 to bin 4.
In Fig. 3 we show the median values, and the corresponding standard deviations over the full sample of 552 circular patches  of D EE and D BB for each selected bin in multipole. On average, the E-B power spectra at these intermediate and low Galactic latitudes are consistent with those presented in P16XXX at high latitude. The histograms of the D BB /D EE ratios for each multipolebin are shown in Fig. 4. These distributions enabled us to choose a specific selection of bins. In the rest of the analysis we considered neither bin 0 nor bin 4, as bin 0 is highly affected by cosmic variance in small sky patches, and bin 4 is contaminated by noise at full Planck resolution (see the corresponding negative tail in D BB /D EE ). As shown in Fig. B.1, neglecting bin 0 allowed us to ensure on the other multipole-bins a level of cosmic variance (∆D /D in the figure) within our 12 • circular patches (θ max in the figure) below 20%.

E-B mode power spectra versus δ N H
Based on the methodology described above we are now able to study variations of D BB /D EE , r T E , r T B , and r EB as a function of δN H for the 552 circular patches at intermediate and low Galactic latitudes.
In Fig. 5 we show two scatter plots: D BB /D EE versus δN H (left panel) and r T E versus δN H (right panel). In the former plot a change in the E-B asymmetry with column density can be clearly seen. In the diffuse ISM, or δN H < 3 × 10 20 cm −2 , D BB / D EE is consistent with the value of 0.52 reported in P16XXX. In denser circular patches the ratio tends to increase toward unity; that is, the amount of power in E and B modes for the densest regions is almost the same. As shown by the inset, this effect appears more important for bin 2 and bin 3, which depart from the value of P16XXX of at least 3σ for δN H > 10 21 cm −2 . However, owing to the larger sample variance of bin 1, this apparent scale dependence must be considered with caution.
In the right panel of the figure an anticorrelation between r T E and δN H can be viewed. As in the left panel, r T E is compatible with diffuse ISM value of 0.36 presented in P18XI for δN H < 3 × 10 20 cm −2 . However, as shown by the linear fit of r T E for δN H > 10 21 cm −2 , for denser regions r T E decreases with column density. The solid line corresponding to the fit could be used to infer the behavior of r T E if data at higher angular resolution were available. A finer angular resolution would allow us to access to larger column densities otherwise smoothed by the Planck beam. This may suggest that r T E would be significantly negative for δN H > 10 22 cm −2 .
Gray shaded areas in both plots define regions in which data noise and systematic effects dominate the signal. These are the only causes producing negative values of the D BB /D EE and values of r T E larger than unity. We want to stress that the overall scatter of the correlations is not primarily caused by noise, as explained in more detail in Appendix A. It is mostly related to sample variance of a non-Gaussian signal, such as that of interstellar dust polarization, in small sky patches across the sky. In the same appendix we also present the 2D probability density function of D BB /D EE and r T E (see Fig. A.2), which shows an intrinsic anticorrelation between the two parameters. Figure 6 shows the dependence of r T B and r EB on δN H . These two parameters are noisier. No dependence on column density can be seen. We also find that, in spite of the large scatter, the median values of r T B and r EB for δN H > 3 × 10 20 cm −2 (see dashed horizontal lines) are systematically larger, and non-zero, at large scales (bin 1 and bin 2) rather than at small scales (bin 3). As explained in Appendix A, this effect is not due to data noise or to the data analysis.
For the less noisy parameters, we also produced N side = 8 maps of D BB /D EE and r T E (see Fig. 7). These show how their variations appear correlated with δN H , with organized, nonrandom, patterns over the intermediate and low latitude sky. As quantified by the scatter plots, the largest values of δN H correspond to largest values of D BB /D EE and to the lowest values of r T E .

Discussion
Our work extends the Planck analysis on the E and B modes of dust polarization at 353 GHz from the diffuse ISM (see P16XXX and P18XI) to denser regions in molecular clouds of the Gould Belt at low Galactic latitude. This study is important both for a better understanding of how interstellar dust affects CMB polarization and for establishing a link between the E-B mode decomposition of dust polarized emission and the ISM physics.
We focused on the link between the variations of the E and B power spectra, and their cross-correlation coefficients (r XY , where X and Y are equal to T , E, or B), with δN H (see Sect. 2.2). We confirmed the average values of the B-to-E power ratio, D BB /D EE , and r T E in the diffuse ISM (δN H < 3 × 10 20 cm −2 ) reported in P16XXX and P18XI. However, for denser regions (δN H > 3 × 10 20 cm −2 ) we found clear departures from these mean values with signs of correlation between D BB /D EE and δN H , and anticorrelation between r T E and δN H . These statistical trends quantified by the scatter plots in Fig. 5 are also shown in the maps presented in Fig. 7. We found an intrinsic anticorrelation between D BB /D EE and r T E as well.
These results strengthened the interpretation of the E-B asymmetry, and the positive T E correlation, in terms of the alignment between the magnetic field orientation and the density filamentary structures in the ISM, as already claimed in Planck Collaboration Int. XXXVIII (2016) for the diffuse medium. A positive r T E and a D BB /D EE less than unity, would be both naturally produced by filamentary structures aligned with the orientation of the interstellar magnetic field (see Z01), which was proved true in the diffuse ISM at high latitude by Planck Collaboration Int. XXXII (2016) and Planck Collaboration Int. XXXVIII (2016).
The same alignment was also observed in the diffuse surrounding of molecular clouds in the Gould Belt. However, these lower latitude regions also present a gradual change in relative orientation, or a smooth transition from parallel to perpendicular, for denser and denser matter structures with respect to the magnetic field (Planck Collaboration Int. XXXV 2016; Soler & Hennebelle 2017). This change in relative orientation is considered representative of the dynamical properties of molecular clouds. Based on the comparison between data and numerical simulations, the change in relative orientation with increasing matter density is indicative of molecular clouds dominated by their self-gravity in sub/trans-Alfvénic MHD turbulent media (Soler et al. 2013). Always following Z01, a relative perpendicular orientation between filamentary density structures and the magnetic field would still produce D BB /D EE < 1, but would have values of r T E < 0. The extrapolation of r T E with δN H in the right panel of Fig. 5 indeed shows that, for δN H > 10 22 cm −2 , r T E may gradually change and become negative. This value of column density is also very close to that quoted in Planck Collaboration Int. XXXV (2016), which corresponds to the change in relative orientation. The smooth change in relative orientation between density filamentary structures and magneticfield orientation would produce a transition in the values of D BB /D EE , which, as shown in the left panel of Fig. 5, would first increase toward unity and decrease again once most of the dense structures would be perpendicular to the magnetic-field orientation. However, at the angular resolution of Planck we do not have access to enough statistics to trace the densest filamentary structures in molecular clouds (Planck Collaboration Int. XXXIII 2016). If the dependence of r T E and D BB /D EE on δN H was indeed related to the change in relative orientation between density and magnetic-field structures, we would also expect a dependence of r T E and D BB /D EE on the scale, as the smaller the scale in molecular clouds the denser the region. In Fig. 5 we show a possible hint of this effect since the increase of D BB /D EE with δN H appears stronger at small scales for bin 2 and bin 3. The angular scales probed by the three multipole bins are 2.25 • , Fig. 7. Maps at N side = 8 of the column-density parameter, δN H (top-center), D BB /D EE (left), and r T E (right) for bins 1, 2, and 3. Nonrandom patterns changing with multipole, and related to the morphology of δN H , can be seen across the sky. 54 arcmin, and 19 arcmin for bin 1, 2, and 3, respectively. For close-by molecular clouds (∼300 pc from us) these sizes would correspond to ∼11 pc, ∼5 pc, and ∼2 pc, which represent largeto medium-size scales in molecular clouds. It would be ideal to reach smaller scales to probe even denser regions. However, going to smaller angular scales (larger multipoles) significantly increases the noise contribution, as shown by bin 4 in Fig. 4. Thus, in order to confirm our interpretation, higher angular resolution polarization surveys designed to probe interstellar dust emission would be necessary (e.g., BFORE, Bryan et al. 2018).
Another result that extends the recent finding of P18XI is that r T B , and maybe r EB , may indeed differ from zero, with a stronger positive signal at large scale compared to small scales. Bracco et al. (2019) suggested that the T B positive correlation at very large scale (for multipoles l < 50) may be principally caused by the Galactic magnetic-field structure in the solar neighborhood, which would leave an imprint of a left-handed helical component on the T B correlation on scales of a few hundred parsecs. However, at the angular scales probed in this work, other processes may be at play, since for the closest Gould Belt clouds we would be probing physical scales of a few parsecs. Further investigation is needed to understand what kind of mechanisms may generate the T B correlation in molecular clouds.
From previous works it is worth noting that most of the effort has been put toward understanding the level of E-B asymmetry.
Our analysis shows that, although such value is on average true in the diffuse ISM, large variations are found across the sky. These variations have organized patterns at intermediate and low Galactic latitudes (Fig. 7), and must be related to intrinsic changes in ISM physics and interstellar dust properties. Kim et al. (2019) used MHD simulations to produce allsky synthetic observations to study the E-B asymmetry. They concluded that the observed power spectra strongly fluctuate depending both on the position of the observer and on temporal fluctuations of ISM properties owing to variations of the star formation process. For the first time, our work shows that the level of E-B asymmetry in real observational data may indeed significantly vary depending on the sky position. However, comparing observational data and all-sky non-Gaussian stochastic models of dust polarization, we showed that most of the variations of E-B modes in the diffuse ISM are likely due to sample variance across the sky rather than to intrinsic physical differences among the sky patches. This is not true in the dense ISM, where the E-B decomposition depends on the value of the gas column density, thus likely on the physics of the observed ISM region. This is important for modeling the impact of dust polarization in CMB studies and for assessing the link between E-B modes and ISM physics.

Summary
We have presented a novel analysis of the Planck polarization data at 353 GHz that extends the study of the T -E-B mode power spectra of interstellar dust to low Galactic latitudes (|b| < 60 • and |b| > 5 • ). We investigated the correlation between these power spectra and the gas column density, which, in the selected sky, is dominated by the emission of molecular clouds in the Gould Belt. Our analysis is relevant to better characterize the statistical properties of dust polarization, both to model Galactic foreground emission to the CMB polarization and to study the dynamical properties of the ISM.
We divided the selected sky in 552 identical circular patches within which we could define mean values of column density, δN H , and of T -E-B power spectra for multipoles between 80 < < 550. We thus studied the respective auto-and crosscorrelations (r XY , with X and Y equal to T, E, B).
We found that the B-to-E power ratio, D BB /D EE , correlates with column density, δN H . While for δN H < 3 × 10 20 cm −2 the values of D BB /D EE are consistent with what was already found in the diffuse ISM (D BB /D EE ≈ 0.5, P16XXX, P18XI), for larger column density the ratio increases approaching unity.
We found that the positive T E correlation observed in the diffuse ISM (r T E ≈ 0.36, P18XI) is on average compatible with our results for δN H < 3 × 10 20 cm −2 . However, for denser regions we found a clear anticorrelation between r T E and δN H , with r T E approaching zero for our densest sample of column density in molecular clouds of the Gould Belt. This trend suggests that r T E could become negative for δN H > 10 22 cm −2 , corresponding to a perpendicular relative orientation between density structures and magnetic field in molecular clouds (see Z01). This would be consistent with the analysis of HROs in dense molecular clouds (i.e., Planck Collaboration Int. XXXV 2016). Only high-resolution polarization surveys of dust emission will allow us to confirm this interpretation.
We found an anticorrelation between r T E and D BB /D EE . We confirmed that, as shown in P18XI, the median value of the r T B may be positive and non-zero at large scale (for multipoles l ≈ 80). We did not find any dependence between δN H and r T B , or r EB , however, this may result from the low signal-to-noise in T B and EB.
We found that the E-B mode dust power spectra show strong variations compared to the mean values reported in previous works. These variations, seen correlated on the sky, are not due to noise. In the diffuse ISM they are mainly caused by small sample variance of a highly non-Gaussian signal such as interstellar dust polarization. In the dense ISM, however, they appear to be correlated with the column density suggesting that we may effectively trace changes of ISM physical properties (i.e., Galactic magnetic-field structure and interstellar turbulence). This is both relevant to model the impact of dust polarization as a CMB foreground and for understanding the link between the E-B mode decomposition and ISM physics.  In this appendix we test the methodology described in Sect. 3.2, performing the same analysis on non-Gaussian simulations of the polarized sky, which have the property of reproducing the 1-and 2-point statistics of the Planck polarization data at high Galactic latitude (Planck Collaboration Int. XLIV 2016; Vansyngel et al. 2017). These simulations are stochastic models of polarized dust emission on the sphere. The method builds on the understanding of Galactic polarization in terms of the structure of the Galactic magnetic field and its coupling with interstellar matter and turbulence through a handful of parameters. The simulated maps do not correspond to Gaussian random fields as shown in Fig. 5 of Vansyngel et al. (2017). We generate two sets of simulations, with and without a noise realization from Planck (including systematic effects), in which the input values of r T E and D BB /D EE are fixed to 0.36 and 0.52, respectively. Results can be seen in Fig. A.1 for the simulated r T E parameter against the observed δN H . The effect of noise does not significantly increase the variance in the correlation plots, which is dominated by the intrinsic variance among the different sky patches in the simulations, as is detailed in the following. Moreover, as expected, no dependence exists between the observed δN H and the simulated r T E . We also notice that, regardless of the multipole bin, the input values for the medians of r T E (solid horizontal lines) are obtained in output. The same is found for the simulated r T B and r EB parameters where the input values are set to 0. These two parameters do not show any systematic decrease in the median values with scale as observed in the Planck data. Thus, we conclude that the decrease in the median values of r T B and r EB observed in the data, from large to small scale, cannot be caused by noise or by our methodology. We suggest that, unless