© A. Bracco et al. 2019

1 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 are ideal to probe the two-point statistics in polarization across the sky.

In the case of CMB, the B mode power wouldpartly 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. (2017, 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 large-scale 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 inSect. 6. Two appendices (Appendices A and B) clarify our data analysis.

2 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.

2.1 Planck polarization data

We used publicly available Planck PR3 data¹ at 353 GHz (Planck Collaboration III 2019) in HEALPix² 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).

2.2 Column density map

We considered the total gas column-density map, N_H, derived from the dust optical depth at 353 GHz, τ₃₅₃. The τ₃₅₃ 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 τ₃₅₃ map at its nominal resolution of 5′. To scale from τ₃₅₃ to N_H we adopted the same convention as in Planck Collaboration Int. XXXV (2016),

Fig. 1 Filtered column density map, δNH, derived from Planck data of interstellar dust emission. The Gould Belt is qualitatively represented by the two white dashed lines. Within the Gould Belt some bright and close-by molecular clouds can be identified (<600 pc from the Sun): Taurus, Perseus, and California in the extreme east (A); Cepheus and Polaris in the northeast (B); Ophiuchus above the Galactic center (C); Musca and Chamaeleon in the southwest (D); Orion in the extreme west (E). A Galactic coordinates grid centered in (l, b) = (0°,0°) is added with steps of 30° both in longitude and latitude.

$${\tau }_{353}/N_{\text{H}}=1.2\times {10}^{-26}{\text{cm}}^2.$$(1)

Variations in dust opacity are present even in the diffuse ISM and the opacity increases systematically by a factor of 2 from the diffuse to the denser ISM (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 gradient, we filtered N_H. The filtered N_H map is $\delta N_{\text{H}}=N_{\text{H}}-{\mathcal{N}}_{\text{H}}$, where ${\mathcal{N}}_{\text{H}}$ is the column-density map smoothed toa 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).

2.3 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², 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 $\delta N_{\text{H}}^{\star }=\langle \delta N_{\text{H}}(>95\%)\rangle$, 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%.

Fig. 2 Mask showing the 552 circular sky patches (gray spots) considered in the present analysis. We masked the sky for Galactic latitudes |b| < 5° and |b| > 60°. The black area is masked. A Galactic coordinates grid centered in (l, b) = (0°, 0°) is added with steps of 30° both in longitude and in latitude.

3 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.

3.1 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 $${(Q+iU)}^{\prime }(n)=e^{\mp 2i\beta }(Q+iU)(n),$$(2)

where n is the position in the sky and β is the rotation of the plane-of-the-sky reference (e₁, e₂) in $e_1^{\prime }=\cos \beta {\text{e}}_1+\sin \beta {\text{e}}_2$ and $e_2^{\prime }=-\sin \beta {\text{e}}_1+\cos \beta {\text{e}}_2$. We note that in the following Stokes I is alternativelyreferred 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 $$\begin{array}{l}T(n)={\sum }_{\mathcal{l}m}{a}_{T,\mathcal{l}m}{Y}_{\mathcal{l}m}(n),\\ (Q+iU)(n)={\sum }_{\mathcal{l}m}{a}_{2,\mathcal{l}{m}_2}{Y}_{\mathcal{l}m}(n),\\ (Q-iU)(n)={\sum }_{\mathcal{l}m}{a}_{-2,\mathcal{l}{m}_{-2}}{Y}_{\mathcal{l}m}(n),\end{array}$$(3)

and used thespin-raising (lowering) operators, ð₊ (ð₋), to get two rotationally invariant quantities $$\begin{array}{ll}{ð}_{-}^2(Q+iU)(n)\hfill & ={\sum }_{\mathcal{l}m}{\left[\frac{(\mathcal{l}+2)!}{(\mathcal{l}-2)!}\right]}^{1/2}{a}_{2,\mathcal{l}m}{Y}_{\mathcal{l}m}(n),\hfill \\ {ð}_{+}^2(Q-iU)(n)\hfill & ={\sum }_{\mathcal{l}m}{\left[\frac{(\mathcal{l}+2)!}{(\mathcal{l}-2)!}\right]}^{1/2}{a}_{-2,\mathcal{l}m}{Y}_{\mathcal{l}m}(n).\hfill \end{array}$$(4)

From Eq. (4), the expansion coefficients are written as $$\begin{array}{l}{a}_{T,\mathcal{l}m}={\displaystyle \int Y_{\mathcal{l}m}^{\ast }}(n)T(n)\text{d}\Omega ,\\ a_{2,\mathcal{l}m}={\left[\frac{(\mathcal{l}+2)!}{(\mathcal{l}-2)!}\right]}^{-1/2}{\displaystyle \int Y_{\mathcal{l}m}^{\ast }}(n){ð}_{-}^2(Q+iU)(n),\\ a_{-2,\mathcal{l}m}={\left[\frac{(\mathcal{l}+2)!}{(\mathcal{l}-2)!}\right]}^{-1/2}{\displaystyle \int Y_{\mathcal{l}m}^{\ast }}(n){ð}_{+}^2(Q-iU)(n),\end{array}$$(5)

which can be linearly combined into $$\begin{array}{ll}{a}_{E,\mathcal{l}m}\hfill & =-(a_{2,\mathcal{l}m}+a_{-2,\mathcal{l}m})/2,\hfill \\ a_{B,\mathcal{l}m}\hfill & =i(a_{2,\mathcal{l}m}-a_{-2,\mathcal{l}m})/2.\hfill \end{array}$$(6)

The E and B modes, i.e., the scalar and pseudo-scalar fields, respectively, are defined as $$\begin{array}{ll}E(n)\hfill & ={\sum }_{\mathcal{l}m}{a}_{E,\mathcal{l}m}{Y}_{\mathcal{l}m}(n)\hfill \\ B(n)\hfill & ={\sum }_{\mathcal{l}m}{a}_{B,\mathcal{l}m}{Y}_{\mathcal{l}m}(n).\hfill \end{array}$$(7)

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 ℓ, $$C_{\mathcal{l}}^{XY}=\frac{1}{2\mathcal{l}+1}{\sum }_m\langle a_{X,\mathcal{l}m}^{\ast }{a}_{Y,\mathcal{l}m}\rangle ,$$(8)

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 $${\mathcal{D}}_{\mathcal{l}}^{XY}=\mathcal{l}(\mathcal{l}+1)C_{\mathcal{l}}^{XY}/(2\pi ).$$(9)

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: $$r^{XY}=\frac{C_{\mathcal{l}}^{XY}}{\sqrt{C_{\mathcal{l}}^{XX}\times C_{\mathcal{l}}^{YY}}},$$(10)

so that in case of perfect positive (negative) correlation r^XY = 1 (−1), and in case of absence of correlation r^XY = 0.

Fig. 3 Median values of ${\mathcal{D}}_{\mathcal{l}}$ in the five bins in multipole (see text) for E (red) and B (blue) modes over the 552 circular regions in Fig. 2. The dashed lines correspond to the best-fit power-law spectra to the observed Planck data. The slopes of about − 2.4 are consistent those what presented in P16XXX.

3.2 Power-spectra analysis

We computed the T–E–B power spectra in Eq. (8) for each circular sky patch using the XPOL³ 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 ${\mathcal{D}}_{\mathcal{l}}^{EE}$ and ${\mathcal{D}}_{\mathcal{l}}^{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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ ratios for each multipole-bin 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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$). As shown in Fig. B.1, neglecting bin 0 allowed us to ensure on the other multipole-bins a level of cosmic variance ($\Delta {\mathcal{D}}_{\mathcal{l}}/{\mathcal{D}}_{\mathcal{l}}$ in the figure) within our 12° circular patches (θ_max in the figure) below 20%.

Fig. 4 Histograms of the ratio ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ for the multipole bins considered in this work (see text). Centers shown in: ℓ = 35 (bin 0), ℓ= 80 (bin 1), ℓ= 200 (bin 2), ℓ= 550 (bin 3), ℓ= 1150 (bin 4); the corresponding median values are overlaid in dashed vertical lines.

Fig. 5 Left panel: ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ vs. the $\delta N_{\text{H}}^{\star }$ for all 552 circular regions in bin 1 (light blue), 2 (green), and 3 (purple). Horizontal dashed lines show: (i) the value of 1, corresponding to equal power in E and B modes; (ii)the value of 0.52, as reported in P16XXX. The inset shows the mean and standard deviation of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ for all multipole bins to highlight the increase of the ratio at high column density. The mean values are computed in bins of ${\log }_{10}\delta N_{\text{H}}^{\star }$ that contain 20 data points each. Right panel: scatter plot of $r_{\mathcal{l}}^{TE}$ vs. $\delta N_{\text{H}}^{\star }$. Colors are same as in the left panel. Horizontal dashed lines show: (i) the value of 0, corresponding to absence of correlation between intensity and E modes; (ii) the value of 0.36, as reported in P18XI. The solid line shows a fit of $r_{\mathcal{l}}^{TE}$ for $\delta N_{\text{H}}^{\star }>{10}^{21}{\text{cm}}^{-2}$. The gray-shaded areas in both panels refer to regions dominated by instrumental noise and data systematic effects.

4 E–B mode power spectra versus δ N ${}_{\text{H}}^{\star }$

Based on the methodology described above we are now able to study variations of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$, r^TE, r^TB, and r^EB as a function of $\delta N_{\text{H}}^{\star }$ for the 552 circular patches at intermediate and low Galactic latitudes.

In Fig. 5 we show two scatter plots: ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ versus $\delta N_{\text{H}}^{\star }$ (left panel) and r^TE versus $\delta N_{\text{H}}^{\star }$ (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 $\delta N_{\text{H}}^{\star }<3\times {10}^{20}$ cm⁻², ${\mathcal{D}}_{\mathcal{l}}^{BB}/$ ${\mathcal{D}}_{\mathcal{l}}^{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 $\delta N_{\text{H}}^{\star }>{10}^{21}$ cm⁻². 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^TE and $\delta N_{\text{H}}^{\star }$ can be viewed. As in the left panel, r^TE is compatible with diffuse ISM value of 0.36 presented in P18XI for $\delta N_{\text{H}}^{\star }<3\times {10}^{20}$ cm⁻². However, as shown by the linear fit of r^TE for $\delta N_{\text{H}}^{\star }>{10}^{21}{\text{cm}}^{-2}$, for denser regions r^TE decreases with column density. The solid line corresponding to the fit could be used to infer the behavior of r^TE 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^TE would be significantly negative for $\delta N_{\text{H}}^{\star }>{10}^{22}$ cm⁻².

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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and values of r^TE 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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and r^TE (see Fig. A.2), which shows an intrinsic anticorrelation between the two parameters.

Figure 6 shows the dependence of r^TB and r^EB on $\delta N_{\text{H}}^{\star }$. 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^TB and r^EB for $\delta N_{\text{H}}^{\star }>3\times {10}^{20}{\text{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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and r^TE (see Fig. 7). These show how their variations appear correlated with $\delta N_{\text{H}}^{\star }$, with organized, nonrandom, patterns over the intermediate and low latitude sky. As quantified by the scatter plots, the largest values of $\delta N_{\text{H}}^{\star }$ correspond to largest values of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and to the lowest values of r^TE.

Fig. 6 Top panel: rTB vs. $\delta N_{\text{H}}^{\star }$. Colors are same as in Fig. 5. Bottom panel: rEB vs. $\delta N_{\text{H}}^{\star }$. Colors aresame as in the top panel. Horizontal dashed lines in both panels show the median values of rTB and rEB for $\delta N_{\text{H}}^{\star }>3\times {10}^{20}{\text{cm}}^{-2}$ in the three multipole bins. Despite the large scatter, a systematic decrease of the median values with angular scale can be seen.

5 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 $\delta N_{\text{H}}^{\star }$ (see Sect. 2.2). We confirmed the average values of the B-to-E power ratio, ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$, and r^TE in the diffuse ISM ($\delta N_{\text{H}}^{\star }<3\times {10}^{20}$ cm⁻²) reported in P16XXX and P18XI. However, for denser regions ($\delta N_{\text{H}}^{\star }>3\times {10}^{20}$ cm⁻²) we found clear departures from these mean values with signs of correlation between ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and $\delta N_{\text{H}}^{\star }$, and anticorrelation between r^TE and $\delta N_{\text{H}}^{\star }$. These statistical trends quantified by the scatter plots in Fig. 5 are also shown in the mapspresented in Fig. 7. We found an intrinsic anticorrelation between ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ and r^TE as well.

These results strengthened the interpretation of the E–B asymmetry, and the positive TE 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^TE and a ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}<1$, but would have values of r^TE < 0. The extrapolation of r^TE with $\delta N_{\text{H}}^{\star }$ in the right panel of Fig. 5 indeed shows that, for $\delta N_{\text{H}}^{\star }>{10}^{22}$ cm⁻², r^TE 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 magnetic-field orientation would produce a transition in the values of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{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^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ on $\delta N_{\text{H}}^{\star }$ was indeed related to the change in relative orientation between density and magnetic-field structures, we would also expect a dependence of r^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{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 ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ with $\delta N_{\text{H}}^{\star }$ appears stronger at small scales for bin 2 and bin 3. The angular scales probed by the three multipole bins are 2.25°, 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 large- to 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^TB, 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 TB 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 TB 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 TB 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 all-sky 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.

Fig. 7 Maps at Nside = 8 of the column-density parameter, $\delta N_{\text{H}}^{\star }$ (top-center), ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ (left), and rTE (right) for bins 1, 2, and 3. Nonrandom patterns changing with multipole, and related to the morphology of $\delta N_{\text{H}}^{\star }$, can be seen across the sky.

6 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, $\delta N_{\text{H}}^{\star }$, and of T–E–B power spectra for multipoles between 80 < ℓ < 550. We thus studied the respective auto- and cross-correlations (r^XY, with X and Y equal to T, E, B).

We found that the B-to-E power ratio, ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$, correlates with column density, $\delta N_{\text{H}}^{\star }$. While for $\delta N_{\text{H}}^{\star }<3\times {10}^{20}$ cm⁻² the values of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ are consistent with what was already found in the diffuse ISM (${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}\approx 0.5$, P16XXX, P18XI), for larger column density the ratio increases approaching unity.

We found that the positive TE correlation observed in the diffuse ISM (r^TE ≈ 0.36, P18XI) is on average compatible with our results for $\delta N_{\text{H}}^{\star }<3\times {10}^{20}$ cm⁻². However, for denser regions we found a clear anticorrelation between r^TE and $\delta N_{\text{H}}^{\star }$, with r^TE approaching zero for our densest sample of column density in molecular clouds of the Gould Belt. This trend suggests that r^TE could become negative for $\delta N_{\text{H}}^{\star }>{10}^{22}$ cm⁻², 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^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$.

We confirmed that, as shown in P18XI, the median value of the r^TB may be positive and non-zero at large scale (for multipoles l ≈ 80). We did not find any dependence between $\delta N_{\text{H}}^{\star }$ and r^TB, or r^EB, however, this may result from the low signal-to-noise in TB 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.

Acknowledgements

We sincerely thank the anonymous referee for improving the paper with her/his valuable comments. We gratefully acknowledge the use of the Aquila cluster at NISER, Bhubaneswar. This research is partly supported by the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022). Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package.

Appendix A: Non-Gaussian simulations of dust polarization

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

Fig. A.1 Scatter plots of modeled rTE obtained with the simulations against the gas column density derived from the data, $\delta N_{\text{H}}^{\star }$, with (right panel) and without (left panel) Planck noise. No dependence with $\delta N_{\text{H}}^{\star }$ is expected. The variance observed in the scatter plots is not dominated by noise. We note that the median values in every multipole bin (see horizontal solid lines) correspond to the input values used to generate the simulations.

Fig. A.2 Probability density functions (pdf, 1D and 2D) of rTE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ only for bin 2 and bin 3 together with 1, 2, and 3σ contour levels for the observed data (left panel), the simulations without Planck noise (central panel), and the simulations with Planck noise added (right panel). The blue lines indicate the observed values in P16XXX and P18XI of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}=0.52$ and rTE = 0.36. On top of the 1D pdf the mean values and standard deviations (also represented by black vertical dashed lines) are shown. These plots were made with the PYTHON CORNER package (Foreman-Mackey 2016).

We generate two sets of simulations, with and without a noise realization from Planck (including systematic effects), in which the input values of r^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ are fixed to 0.36 and 0.52, respectively. Results can be seen in Fig. A.1 for the simulated r^TE parameter against the observed $\delta N_{\text{H}}^{\star }$. The effect ofnoise does not significantly increase the variance in the correlation plots, which is dominated by the intrinsic variance among the differentsky patches in the simulations, as is detailed in the following. Moreover, as expected, no dependence exists between the observed $\delta N_{\text{H}}^{\star }$ and the simulated r^TE. We also notice that, regardless of the multipole bin, the input values for the medians of r^TE (solid horizontal lines) are obtained in output. The same is found for the simulated r^TB 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^TB 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 residual (unknown) systematic effects in the data are present, the observed decrease may be true. However, due to the large scatter in the distributions of the observed r^TB and r^EB, it is not possible to quantify the significance of this statement.

In Fig. A.2 we show 1D and 2D probability density functions for r^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ for the observed and the simulated data, respectively, considering bin 2 and bin 3 together to increase the number statistics. The two parameters appear clearly anticorrelated both in the observations and in the simulations. From Eq. (5) in Vansyngel et al. (2017) the inverse dependence in the simulated spectra can be derived as ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}\approx -m{(r^{TE})}^2+q$, where m = 0.54 and q = 0.56. The effect of noise smooths the anticorrelation between r^TE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$, suggesting that the true anticorrelation in the Planck data is likely stronger. The spread about the mean values found in the noisy simulations allows us to statistically recover the observed data dispersion (see the standard deviations quoted in the figure), confirming that sample variance is a major responsible for the T–E–B power fluctuations from patch to patch over the sky at least in the diffuse ISM. This result validates the simulations presented in Vansyngel et al. (2017) for the statistical description of the polarized properties of the diffuse ISM even in small sky patches. However, as proved by our work, a significant dependence of the parameters with column density for $\delta N_{\text{H}}^{\star }>3\times {10}^{20}$ cm⁻² is observed. This is not captured yet by any existing model.

Appendix B: Cosmic variance per multiple bin

We show a figure that allows us to quantify the level of cosmic variance in each multipole bin used in the data analysis. Following Tegmark (1997), the cosmic variance can be estimated as $$\frac{\Delta {\mathcal{D}}_{\mathcal{l}}}{{\mathcal{D}}_{\mathcal{l}}}\approx \sqrt{\frac{2}{(2\mathcal{l}+1)L{f}_{\text{sky}}}},$$(B.1)

where f_sky is the sky fraction considered and related to the sky-patch size as f_sky = sin²(θ_max∕2), which in ourcase is θ_max = 12°; and L is the width of the ℓ-bin equal to 15, 40, 200, 500, and 1200 from bin 0 to bin 4, respectively. As shown in Fig. B.1, neglecting bin 0 enables us to limit the level of cosmic variance per bin below 20%. We note that this equation is not completely accurate for cross-spectra. In that case it would read $\Delta D_{\mathcal{l}}^{XY}=(1/{\nu }_{\mathcal{l}})[(D_{\mathcal{l}}^{XY})+D_{\mathcal{l}}^{XX}{D}_{\mathcal{l}}^{YY}]$. Moreover, these estimates are only valid in case of Gaussian random fields. The observed signal is not Gaussian, thus we expect a larger amount of variance per bin of a factor of a few (Vansyngel et al. 2017).

Fig. B.1 Cosmic variance per multipole bin vs. sky-patch size and sky fraction (Tegmark 1997). The central multipoles for each bin are as follows: ℓ = 35 (bin 0), ℓ = 80 (bin 1), ℓ = 200 (bin 2), ℓ = 550 (bin 3), and ℓ = 1150 (bin 4).

References

All Figures

Fig. 1 Filtered column density map, δNH, derived from Planck data of interstellar dust emission. The Gould Belt is qualitatively represented by the two white dashed lines. Within the Gould Belt some bright and close-by molecular clouds can be identified (<600 pc from the Sun): Taurus, Perseus, and California in the extreme east (A); Cepheus and Polaris in the northeast (B); Ophiuchus above the Galactic center (C); Musca and Chamaeleon in the southwest (D); Orion in the extreme west (E). A Galactic coordinates grid centered in (l, b) = (0°,0°) is added with steps of 30° both in longitude and latitude.

In the text

	Fig. 2 Mask showing the 552 circular sky patches (gray spots) considered in the present analysis. We masked the sky for Galactic latitudes |b| < 5° and |b| > 60°. The black area is masked. A Galactic coordinates grid centered in (l, b) = (0°, 0°) is added with steps of 30° both in longitude and in latitude.
In the text

	Fig. 3 Median values of ${\mathcal{D}}_{\mathcal{l}}$ in the five bins in multipole (see text) for E (red) and B (blue) modes over the 552 circular regions in Fig. 2. The dashed lines correspond to the best-fit power-law spectra to the observed Planck data. The slopes of about − 2.4 are consistent those what presented in P16XXX.
In the text

	Fig. 4 Histograms of the ratio ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ for the multipole bins considered in this work (see text). Centers shown in: ℓ = 35 (bin 0), ℓ= 80 (bin 1), ℓ= 200 (bin 2), ℓ= 550 (bin 3), ℓ= 1150 (bin 4); the corresponding median values are overlaid in dashed vertical lines.
In the text

Fig. 5 Left panel: ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ vs. the $\delta N_{\text{H}}^{\star }$ for all 552 circular regions in bin 1 (light blue), 2 (green), and 3 (purple). Horizontal dashed lines show: (i) the value of 1, corresponding to equal power in E and B modes; (ii)the value of 0.52, as reported in P16XXX. The inset shows the mean and standard deviation of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ for all multipole bins to highlight the increase of the ratio at high column density. The mean values are computed in bins of ${\log }_{10}\delta N_{\text{H}}^{\star }$ that contain 20 data points each. Right panel: scatter plot of $r_{\mathcal{l}}^{TE}$ vs. $\delta N_{\text{H}}^{\star }$. Colors are same as in the left panel. Horizontal dashed lines show: (i) the value of 0, corresponding to absence of correlation between intensity and E modes; (ii) the value of 0.36, as reported in P18XI. The solid line shows a fit of $r_{\mathcal{l}}^{TE}$ for $\delta N_{\text{H}}^{\star }>{10}^{21}{\text{cm}}^{-2}$. The gray-shaded areas in both panels refer to regions dominated by instrumental noise and data systematic effects.

In the text

Fig. 6 Top panel: rTB vs. $\delta N_{\text{H}}^{\star }$. Colors are same as in Fig. 5. Bottom panel: rEB vs. $\delta N_{\text{H}}^{\star }$. Colors aresame as in the top panel. Horizontal dashed lines in both panels show the median values of rTB and rEB for $\delta N_{\text{H}}^{\star }>3\times {10}^{20}{\text{cm}}^{-2}$ in the three multipole bins. Despite the large scatter, a systematic decrease of the median values with angular scale can be seen.

In the text

	Fig. 7 Maps at Nside = 8 of the column-density parameter, $\delta N_{\text{H}}^{\star }$ (top-center), ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ (left), and rTE (right) for bins 1, 2, and 3. Nonrandom patterns changing with multipole, and related to the morphology of $\delta N_{\text{H}}^{\star }$, can be seen across the sky.
In the text

Fig. A.1 Scatter plots of modeled rTE obtained with the simulations against the gas column density derived from the data, $\delta N_{\text{H}}^{\star }$, with (right panel) and without (left panel) Planck noise. No dependence with $\delta N_{\text{H}}^{\star }$ is expected. The variance observed in the scatter plots is not dominated by noise. We note that the median values in every multipole bin (see horizontal solid lines) correspond to the input values used to generate the simulations.

In the text

Fig. A.2 Probability density functions (pdf, 1D and 2D) of rTE and ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}$ only for bin 2 and bin 3 together with 1, 2, and 3σ contour levels for the observed data (left panel), the simulations without Planck noise (central panel), and the simulations with Planck noise added (right panel). The blue lines indicate the observed values in P16XXX and P18XI of ${\mathcal{D}}_{\mathcal{l}}^{BB}/{\mathcal{D}}_{\mathcal{l}}^{EE}=0.52$ and rTE = 0.36. On top of the 1D pdf the mean values and standard deviations (also represented by black vertical dashed lines) are shown. These plots were made with the PYTHON CORNER package (Foreman-Mackey 2016).

In the text

	Fig. B.1 Cosmic variance per multipole bin vs. sky-patch size and sky fraction (Tegmark 1997). The central multipoles for each bin are as follows: ℓ = 35 (bin 0), ℓ = 80 (bin 1), ℓ = 200 (bin 2), ℓ = 550 (bin 3), and ℓ = 1150 (bin 4).
In the text