© ESO 2020

1 Introduction

Gravitational waves, possibly generated during the inflation epoch introduce B-mode polarization in the cosmic microwave background (CMB; Starobinskij 1979; Fabbri & Pollock 1983). The measurement of this primordial CMB B-mode signal is the intense focus of ongoing and proposed CMB experiments. Galactic polarized foregrounds, especially polarized thermal dust emission (Martin 2007; Draine & Fraisse 2009; Vaillancourt & Matthews 2012; Planck Collaboration Int. XXII 2015; Ashton et al. 2018; Planck Collaboration XII 2020), also produce a B-mode signal at microwave frequencies (Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2020). Accurate modelling and separation of polarized dust emission is therefore a necessary step in the detection of a primordial CMB B-mode signal. Residual foregrounds due to imperfect component separation can be misinterpreted as detection of primordial CMB B-modes (Remazeilles et al. 2016).

Polarized thermal dust emission is expected from aspherical grains aligned with respect to the local Galactic magnetic field (GMF; Stein 1966). Many distinct mechanisms have been elaborated to account for the alignment of dust grains, starting from the first quantitative approach by Davis & Greenstein (1951); for a review, see Lazarian (2007) and Andersson et al. (2015).

Empirically, the all-sky map of dust polarization from Planck at 353 GHz (Planck Collaboration I 2016) reveals a connection between the dust intensity structures and the local orientation of the GMF projected on the plane of the sky (B_POS) (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXV 2016). In the low-column-density or diffuse interstellar medium (ISM), the orientation of the dust intensity structures is observed to be preferentially parallel to B_POS (Planck Collaboration Int. XXXII 2016). Going to high-column-density regions, including star-forming molecular clouds, the preferred orientation of the dust intensity structures changes from parallel to perpendicular (Planck Collaboration Int. XXXV 2016; Soler et al. 2017; Jow et al. 2018). This change in relative orientation was predicted using sub- and trans-Alfvenic magnetohydrodynamic (MHD) simulations (Soler et al. 2013), and highlights the role of the magnetic field during the formation of the molecular clouds.

Away from the Galactic plane, Planck polarization maps reveal a large scatter in the distributions of the polarization fraction, p, and the polarization angle, ψ (Planck Collaboration Int. XIX 2015). A comparison with maps computed from a simulation of MHD turbulence in Planck Collaboration Int. XX (2015) shows that the large scatter of p is associated with variations in the orientation of the GMF along the line of sight (LOS), causing a depolarization effect. The simulation also reproduces the inverse relationship between p and polarization angle dispersion function, $\mathcal{S}$ (Planck Collaboration Int. XIX 2015). Subsequently, Planck Collaboration Int. XLIV (2016) connected the distributions of p and ψ with the amplitude of turbulent magnetic field using a phenomenological model. This model is discussed further and compared with the Planck data in Planck Collaboration XII (2020).

The Planck maps were also used to measure power spectra of dust polarization. At intermediate and high Galactic latitudes, in the multipole range, 40 < ℓ < 600, the ratio of dust EE to BB power ($D_{\mathcal{l}}^{EE}/D_{\mathcal{l}}^{BB}$) is found to be around 2 (Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2020). Planck Collaboration Int. XXXVIII (2016) show that both the E–B power asymmetry and the correlation between dust temperature and E mode polarization, i.e. the dust TE correlation, can be accounted for by the alignment between the orientation of the filamentary structure of interstellar matter in the diffuse ISM, as traced by total dust emission, and the orientation of the GMF, as inferred from dust polarization (Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XXXVIII 2016). Caldwell et al. (2017) and Kandel et al. (2018) investigated whether these observed properties could be related to turbulence in the magnetised ISM, considering the contributions from slow, fast, and Alfven MHD modes. Caldwell et al. (2017) conclude that the E–B power asymmetry and positive TE correlation cannot both be accounted for with their model. Using the same theoretical framework, Kandel et al. (2017, 2018) reached the opposite conclusion, assuming that any correlation between the gas density and the magnetic field is negligible. However, this assumption is challenged by the above-mentioned observed alignment between the filamentary structure and the magnetic field. Filamentary structures identified in HI 21-cm line spectroscopic data cubes, which were shown to trace density structure in the cold neutral medium (CNM; Clark et al. 2019), are also found to be aligned with the GMF (Clark et al. 2014, 2015; Martin et al. 2015; Kalberla et al. 2016). Furthermore, Clark (2018) reported a correlation between the polarization fraction, p, and the degree of coherence of the orientation of HI emission features along the LOS, later used by Clark & Hensley (2019) to model dust polarization. Clearly, HI data contain valuable information on the structure of the GMF in the neutral atomic ISM.

Ghosh et al. (2017, hereafter TG17) present a phenomenological model of dust polarization. This TG17 model combines the framework introduced by Planck Collaboration Int. XLIV (2016) with a decomposition of HI emission data into three distinct maps of HI column density referred to as “HI templates”. By adjusting just a few parameters, the TG17 model reproduces the one- and two-point statistical properties of dust polarization over that fraction of the southern Galactic cap (defined by the region of b ≤ −30°) where dust and HI emission are well correlated. In this paper, we use this framework to fit the Planck PR3 data (Planck Collaboration I 2020) over a fraction of the northern Galactic cap at b ≥ 30°. The main goal is to extend the sky area available to fit and test the TG17 model.

This paper is organised as follows. In Sect. 2, we describe the data used from Planck and the Effelsberg-Bonn HI Survey (EBHIS). Section 3 describes the choice of sky region and HI templates used in our analysis. We summarise the statistical properties of the dust polarization over the selected region in Sect. 4. In Sect. 5, we briefly describe the phenomenological modelling framework. Section 6 describes how the model parameters were evaluated. We present separate aspects of the results in Sects. 7 and 8. Finally, we summarise our results in Sect. 9. In Appendix A, we explore a model based on alternative HI templates.

2 Datasets used

2.1 Planck dust polarization maps

In this paper, we use publicly available Planck PR3 data¹ (Planck Collaboration I 2020) at 353 GHz to study the statistical properties of the dust polarization. These maps are produced using only the polarization-sensitive bolometers (PSBs) and are expressed in thermodynamic temperature units (K_CMB, Planck Collaboration III 2020). We also use various subsets of the Planck polarization data at 353 GHz, namely, the half-mission maps (HM1 and HM2), yearly surveys (Y1 and Y2), and odd and even surveys (O and E) to debias the effect of instrumental noise (Sect. 4).

For dust intensity, we use the generalised needlet internal linear combination (GNILC) processed Stokes I map (Planck Collaboration IV 2020). Planck does not have the ability to measure absolute emission, and so it is necessary to correct for the zero level of the dust intensity map. By construction, the GNILC dust intensity map has a cosmic infrared background (CIB) monopole contribution of 452 μK_CMB at 353 GHz, which needs to be subtracted (Planck Collaboration III 2020). To correct for the Galactic HI offset and the contributionfrom dust emission associated with HII emission (Planck Collaboration XII 2020), a “fiducial” offset of 63 μK_CMB is added back to the CIB-subtracted GNILC dust intensity map (hereafter, we refer to this final map as I_G353). The fiducial offset has an uncertainty of 40 μK_CMB, with an associated “low” and “high” offset of 23 μK_CMB and 103 μK_CMB, respectively.

The Planck Stokes Q and U maps at 353 GHz have a beam resolution of4.′82 (FWHM), and the I_G353 map has the beam resolution of5′ (FWHM). To increase the signal-to-noise ratio (S/N), we smooth the dust Stokes I, Q, and U maps to a common resolution of 60′ or 80′ (FWHM) and reproject on the HEALPix grid² (Górski et al. 2005) with N_side = 128. To compare our dust model with the Planck data, we work with the 60′ (FWHM) smoothed maps. However, to study the inverse relationship between $\mathcal{S}$ and p (Sect. 8) we work with 80′ (FWHM) smoothed Planck data.

2.2 EBHIS and HI4PI H I data

We use HI 21-cm line spectroscopic data from EBHIS³ (Kerp et al. 2011; Winkel et al. 2016), which mapped the Milky Way gas in the northern Galactic sky with the 100 m telescope at Effelsberg. The survey has an angular resolution of 10.′8 (FWHM), spectral resolution δv = 1.44 km s⁻¹ (FWHM), and rms brightness temperature uncertainty of 90 mK. Velocities are with respect to the local standard of rest (LSR). The EBHIS data are projected on a HEALPix grid with N_side = 1024 (θ_pix = 3.′4). The EBHIS data also form the northern part of the all-sky product HI4PI (HI4PI Collaboration 2016).

For optically thin emission, the total HI column density (N_HI) can be obtained by integrating the brightness temperature (T_b) over velocity channels (Wilson et al. 2009): $$\begin{array}{ll}{\text{N}}_{\text{HI}}\hfill & =1.82\times {10}^{18}\times {\displaystyle \int T_b}\text{d}v{\text{cm}}^{-2}.\hfill \end{array}$$(1)

Fig. 1 Histogram of residuals at different iterations. The blue curve is the Gaussian fitted to residuals at the final iteration. Gray shadedregions are 1 σcg, 2 σcg, and 3 σcg regions of the fitted Gaussian.

3 Exploiting the HI data

3.1 Region selection

We describe the procedure to select the low-column-density region in the northern Galactic cap in which HI and dust emission are highly correlated, making HI a good proxy for dust.

Ultraviolet observations of molecular hydrogen, H₂ (Savage et al. 1977; Gillmon et al. 2006)and early results from Planck (Planck Collaboration XXIV 2011) show that dust emission associated with gas in the form of H₂ becomes significant for sight lines where the total N_H exceeds 4 × 10²⁰ cm⁻² (Arendt et al. 1998). Therefore, we restrict ourselves to the sky region in which N_HI,50 is below a threshold 3.8 × 10²⁰ cm⁻², where N_HI,50 from Eq. (1) is evaluated over the restricted velocity interval |v| ≤ 50 km s⁻¹. We explored extending the range up to 80 km s⁻¹ and because there is little additional gas in this range our results below from the modelling analysis (Sects. 6–8) are robust. We also get the same model parameters for regions restricted by N_HI,50 thresholds in the range 3.4 − 4.0 × 10²⁰ cm⁻². The sky fraction increases beneficially with threshold, but not significantly beyond the adopted 3.8 × 10²⁰ cm⁻². However, the model parameters change significantly for thresholds above 4 × 10²⁰ cm⁻², presumably because of the increased and unaccounted molecular fraction.

We build our mask using an iterative correlation analysis method as described in Planck Collaboration Int. XVII (2014). The initial mask contains unmasked pixels for which both the EBHIS and Planck data are available over the northern Galactic cap. At each iteration, we compute the linear correlation between I_G353 and total N_HI over unmasked pixels of the binary mask produced in the previous iteration. We compute residuals by subtracting the fit to the correlation from the observed dust emission for each unmasked pixel and find the standard deviation (σ_cg) characterising the Gaussian core of these residuals (Fig. 1).

Next, the binary mask is updated by masking pixels for which the absolute value of the residual is greater than 3 σ_cg. A stable mask is obtained after five iterations, a region of 5900 deg² comprising 65% of the northern Galactic cap (b ≥ 30°), as shown in Fig. 2.

Hereafter we refer to this region as mask65. For later analysis, we apodize this binary mask by convolving it with a Gaussian window function of 2° (FWHM). After apodization, we have an effective sky fraction, $f_{\text{sky}}^{\text{eff}}$, of 0.143 (14.3%). In the stable mask, there are some isolated islands with very few pixels. When we derive the full-sky angular power spectrum C_ℓ from the pseudo-spectrum, these isolated islands introduce correlations from very high ℓ to low ℓ via the mode-mode coupling matrix (Hivon et al. 2002). In our analysis, we eliminate isolated islands containing less than 20 pixels using the Process_mask routine of HEALPix.

Over mask65 the mean HI column density over the restricted velocity interval is $\langle N_{\text{HI},50}\rangle =1.85\times {10}^{20}{\text{cm}}^{-2}$ and the mean dust intensity is $\langle I_{\text{G353}}\rangle =$293 μK_CMB at 353 GHz. The linear correlation between I_G353 and total N_HI yields a slope or emissivity ɛ_353,50 of 137 μK_CMB ${({10}^{20}{\text{cm}}^{-2})}^{-1}$, where againthe subscript 50 denotes the restricted velocity interval. This value is about 2% higher than that found for the diffuse sky studied by Planck Collaboration XI (2014). The offset of 40 μK_CMB is satisfactory, within the systematic uncertainty of the zero level (Sect. 2.1).

Fig. 2 Binary mask defining mask65, the dark regions selecting 65% of the northern Galactic cap (b ≥ 30°).

3.2 HI templates

The diffuse ISM is a complex turbulent multiphase medium. The HI gas comprises two thermally stable phases, the cold and warm neutral medium (CNM and WNM, respectively), along with an additional thermally unstable phase at an intermediate temperature, hereafter referred to as the unstable neutral medium (UNM) following TG17.

HI spectra can be used to build maps of HI emission associated notionally with the CNM, UNM, and WNM. First, the observed brightness temperature profiles of HI spectra are decomposed into Gaussian components (Haud 2000, 2013) under the assumption that the random velocities in any HI “cloud” have a Gaussian distribution. Thus, $$T_{\text{b}}(v)={\sum }_i{T}_0^i\exp \left[-\frac{1}{2}{\left(\frac{v_{LSR}-v_{\text{c}}^i}{{\sigma }^i}\right)}^2\right],$$(2)

where the summation is over all Gaussian components needed to describe the profile along a LOS, v_LSR is the velocity of the HI gas, and $T_0^i$, $v_{\text{c}}^i$, and σⁱ are the peak brightness temperature (expressed in Kelvin), central peak velocity, and 1 σ standard deviation of Gaussian component i, respectively.

Following the detailed physically motivated prescription in Sect. 2.2 of TG17, we produce CNM, UNM, and WNM column density maps. This involves using weight factors applied in summing the contributions from the Gaussian components. The weights involve two parameters of the polarization model of TG17 that constrain the range of centroid velocities of the Gaussians that contribute to each of the three maps. Here we simply adopt the TG17 values for these two parameters and the same weighting scheme, which should not be materially different in the north as compared to the south. TG17 noted the beneficial effect of having considerable column density in the three maps and that although the CMN contains some higher velocity dispersion gas, subdivision into two has little effect on the model. As we show below, this prescription does indeed lead to an acceptable polarization model for mask65.

These three maps have an original resolution of 10.′8 (FWHM), which we have smoothed to a common resolution of 60′ (FWHM) at N_side = 128. The resulting maps, to be used as HI templates below, are shown in Fig. 3. In Appendix A, we present an analysis that uses alternative HI templates based on a decomposition by Kalberla & Haud (2018).

Fig. 3 Northern Galactic orthographic projection of integrated HI emission maps of the three HI templates called CNM (left panel), UNM (middle panel), and WNM (right panel). Over our selected region mask65 (Sect. 3.1), the mean column densities of the three templates are $\langle N_{\text{HI}}^{\text{c}}\rangle =0.90\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{u}}\rangle =0.32\times {10}^{20}{\text{cm}}^{-2}$, and $\langle N_{\text{HI}}^{\text{w}}\rangle =0.56\times {10}^{20}{\text{cm}}^{-2}$, respectively.

4 Observed statistics from Planck dust observations over the selected region

In this section, we analyse the statistical properties of dust polarization over mask65 in both pixel and harmonic space.

4.1 Polarization fraction

To avoid a bias from data noise, we compute the square of the polarization fraction ($p_{\text{d}}^2$) combining independent subsets of the Planck data at 353 GHz: $$p_{\text{d}}^2=<\frac{Q_{\text{d}}^{{\text{s}}_1}{Q}_{\text{d}}^{{\text{s}}_2}+U_{\text{d}}^{{\text{s}}_1}{U}_{\text{d}}^{{\text{s}}_2}}{I_{\text{G353}}{}^2}>,$$(3)

where “d” stands for the Planck data, (s₁, s₂) stands for different subsets of Planck data at 353 GHz: {(HM1, HM2); (Y1, Y2); (odd, even)}, and $<\mathrm{...}>$ implies an average over different data subsets. In the left panel of Fig. 4, we present the histogram of $p_{\text{d}}^2$. The uncertainties (1 σ) in $p_{\text{d}}^2$ are the quadratic sum of data systematic errors estimated from the three data subsets and the uncertainty in I_G353 associated with the determination of the zero level of Galactic emission. The mean values of p_d over mask65 are listed in Table 1 for the different data subsets.

4.2 Polarization angle

We computed the Stokes parameter ratios, q_d = Q_d∕I_G353 and u_d = U_d∕I_G353 from Planck data at 353 GHz. As shown in the top panel of Fig. 5, maps of q_d and u_d centred at the northern Galactic pole reveal a “butterfly” pattern associated with an ordered large-scale GMF.

To determine l₀ and b₀ for this large-scale GMF, we fit these data with a normalised model of Stokes parameters as described in “step A” in Planck Collaboration Int. XLIV (2016). The best-fit values of l₀ and b₀ for different data subsets are noted in Table 1 and the lower panel of Fig. 5 presents the best-fit model. These values, (l₀, b₀) = (62°, −18°), give the “orientation” of the GMF, but not the direction; i.e. the π ambiguity cannot be resolved using these data, though it can be resolved using rotation measures (Xu & Han 2019). The corresponding values found in TG17 for the southern Galactic cap are (l₀, b₀) = (74°, +24°). The change of sign of b₀ from north to south reflects the deformation of the large-scale GMF by the local bubble (Alves et al. 2018).

The mean best-fit values of l₀ and b₀ are then used to calculate a polarization angle map, ψ₀, for the mean GMF, where ψ₀ is the polarization angle of the model without the turbulent component of the GMF. At each sky pixel, we then rotate Q and U values with respect to the new reference angle ψ₀ at respective pixels using the relation, $$\left[\begin{array}{c}{Q}_{\text{d}}^{\text{R}}\\ U_{\text{d}}^{\text{R}}\end{array}\right]=\left(\begin{array}{cc}\cos 2{\psi }_0& -\sin 2{\psi }_0\\ \sin 2{\psi }_0& \cos 2{\psi }_0\end{array}\right)\left[\begin{array}{c}{Q}_{\text{d}}\\ U_{\text{d}}\end{array}\right],$$(4)

where $Q_{\text{d}}^{\text{R}}$ and $U_{\text{d}}^{\text{R}}$ are Stokes parameters in the rotated frame. The polarization angle ${\psi }_{\text{d}}^{\text{R}}$ is derived using the relation ${\psi }_{\text{d}}^{\text{R}}=0.5\text{atan2}(-U_{\text{d}}^{\text{R}},Q_{\text{d}}^{\text{R}})$, where the two-argument function atan2(−U, Q) is used in place of atan(−U∕Q) to avoid the π-ambiguity. The minus sign is introduced to compute ${\psi }_{\text{d}}^{\text{R}}$ in the IAU convention (Hamaker & Bregman 1996) from HEALPix maps given in the COSMO convention⁴ used by Planck.

Figure 4 (right panel) presents the distribution of ${\psi }_{\text{d}}^{\text{R}}$ over mask65. The 1 σ error bars on the distribution of ${\psi }_{\text{d}}^{\text{R}}$ are computed from the standard deviation of the results from the six independent subsets of Planck data in Table 1. The 1 σ dispersion of ${\psi }_{\text{d}}^{\text{R}}$ derived from a Gaussian fit to the distribution is 16. °6 ± 0. °2.

Fig. 4 Left panel: histogram of the square of the polarization fraction, p2, at 353 GHz. Right panel: histogram of ψR, the difference between the measured polarization angle at 353 GHz and that computed for the mean GMF orientation over mask65. Green triangles are mean $p_{\text{d}}^2$ (${\psi }_{\text{d}}^{\text{R}}$) estimated from Planck data subsets, and black squares are mean $p_{\text{m}}^2$ (${\psi }_{\text{m}}^{\text{R}}$) estimated from 100 realisations of model maps. The vertical dashed line in the left panel corresponds to the model parameter p0 = 15.5% ($p_0^2=0.024$), while the vertical dashed line in the right panel corresponds to a model without a turbulent component of the GMF.

4.3 Dust power spectra

We analyse the dust polarization power spectra at 353 GHz over mask65⁵, $D_{\mathcal{l}}^{EE}$, $D_{\mathcal{l}}^{BB}$, and $D_{\mathcal{l}}^{TE}$, in the multipole range 40 ≤ ℓ ≤ 160. These were computed using Xpol (Tristram et al. 2005) from cross spectra of three datasets: half-missions (HM1 × HM2), years (Y1 × Y2), and odd-even surveys (Odd × Even). In Fig. 6, we present the mean of the three cross-power spectra in six multipole bins. The dust $D_{\mathcal{l}}^{EE}$ and $D_{\mathcal{l}}^{BB}$ spectra are corrected for CMB contributions in harmonic space using the Planck best-fit power spectra (Planck Collaboration VI 2020), whereas we keep $D_{\mathcal{l}}^{TE}$ unaltered. Uncertainties (1 σ) on the binned dust power spectra are the quadratic sum of statistical noise computed by Xpol analytically, and systematic uncertainties computed from the standard deviation from the three Planck data subsets.

We also checked the CMB correction at the map level. We subtracted the half-mission and odd-even component-separated CMB SMICA and SEVEM maps (Planck Collaboration IV 2020) from the total Planck maps and then computed cross-power spectra using these four subsets of the maps. Subtraction of CMB does not introduce any noticeable changes in the dust polarization power spectra.

Binned dust power spectra are well described by a power-law model, $D_{\mathcal{l}}^{XX}=A_{XX}{(\mathcal{l}/80)}^{{\alpha }_{XX}+2}$, where A_XX is the best-fit amplitude at ℓ = 80, α_XX is the best-fit spectral index, and again XX = {EE, BB, TE} (Planck Collaboration Int. XXX 2016).The best-fit values of A_XX and α_XX and the respective values of χ² are listed in the middle column of Table 2.

The ratio of BB to EE power is about 0.6 over mask65, consistent with the E–B asymmetry result of Planck Collaboration Int. XXX (2016). We detect a significant positive TE correlation over mask65, as in TG17 for the southern Galactic cap.

Fig. 5 Upper panel: northern Galactic orthographic projection of qd (left) and ud (right) for the Planck data. Lower panel: same for the best-fit model where the large-scale GMF is directed toward (l0, b0) = (62°, −18°) and scaled with a mean polarization fraction of 8.8%.

Fig. 6 Left panel: dust EE, BB, and TE cross-power spectra computed from the subsets of Planck data at 353 GHz over mask65. Dashed lines represent the best-fit power laws. Right panel: similar plots as left panel computed from 100 realisations of the dust model maps. Error bars are 1 σ uncertainties as explained in the main text. The filled areas represent the Planck dust power spectra measurements over mask65.

5 Multiphase model of polarized dust emission

We use the modelling framework from TG17, which is based on a decomposition of HI emission-line data and incorporates the phenomenological magnetic field model described in Planck Collaboration Int. XLIV (2016). We briefly describe the salient concepts.

We use the fact that dust emission and N_HI are correlated so that templates of N_HI can be used as proxies for spatially variable dust intensity. For optically thin dust emission, the model Stokes parameters I_m, Q_m, and U_m at 353 GHz can be written as $$\begin{array}{l}{I}_{\text{m}}(\widehat{n})= {\sum }_{i=1}^N\left[1-p_0\left({\cos }^2{\gamma }^i(\widehat{n})-\frac{2}{3}\right)\right]{\varepsilon }_{353}{\text{N}}_{\text{HI}}^i(\widehat{n})\\ Q_{\text{m}}(\widehat{n})= {\sum }_{i=1}^N{p}_0{\cos }^2{\gamma }^i(\widehat{n})\cos 2{\psi }^i(\widehat{n}){\varepsilon }_{353}{\text{N}}_{\text{HI}}^i(\widehat{n})\\ U_{\text{m}}(\widehat{n})=-{\sum }_{i=1}^N{p}_0{\cos }^2{\gamma }^i(\widehat{n})\sin 2{\psi }^i(\widehat{n}){\varepsilon }_{353}{\text{N}}_{\text{HI}}^i(\widehat{n}),\end{array}$$(5)

where N is the number of distinct templates, $\widehat{n}$ is the direction vector, p₀ is a parameter related to the dust grain properties (Lee & Draine 1985; Draine & Fraisse 2009; Planck Collaboration Int. XX 2015), γ is the angle made by the local magnetic field with the plane of the sky, ψ is the polarization angle measured counterclockwise from Galactic north⁶, ɛ₃₅₃ is the dust emissivity at 353 GHz for each HI template, and N_HI is the column density within the particular template. In this work, N = 3, which is represented by the three distinct templates for CNM, UNM, and WNM (Sect. 3.2). This is a remarkably small number to describe the ISM. However, our focus is on demonstrating what essentials can nevertheless be captured bya polarization model that is inherently minimalist and not dependent on fine-tuning. For simplicity, ɛ₃₅₃ is assumed to be the same for all three templates.

The connection between the dust polarization and the structure of the GMF via the angles γ and ψ is developed as in Planck Collaboration Int. XLIV (2016). Following Jaffe et al. (2010), the GMF, B, is expressed as a vector sum of an ordered component (B_ord) and a fluctuating component (B_turb). Because most of the dust emission comes from the thin Galactic disk, Planck Collaboration Int. XLIV (2016) assume that there is an ordered large-scale magnetic field, B_ord, in the solar neighbourhood (approximately 200 parsec). The butterfly pattern discussed in Sect. 4.2 supports the simplifying assumption that each template has the same B_ord, oriented toward the l₀ and b₀ that are already determined. However, the components of B_turb are different in each template and are taken from independent realisations of a Gaussian random field that has a power spectrum with C_ℓ ∝ ℓ^α_M for ℓ ≥ 2. The strength of |B_turb| relative to |B_ord| is parameterised by $f_{\text{M}}^i$. This phenomenological model captures schematically the association between the structure of the multiphase ISM and that of the GMF, an interdependence that we consider to be essential for modelling the dust polarization. This interdependence is usually ignored in 3D models of the GMF, for example those with Gaussian random magnetic fields as in Levrier et al. (2018) and Wang et al. (2020). However, our approach does not allow us to include the divergence-free constraint because it is not a 3D model.

Here, this key interdependence is introduced explicitly by aligning the model GMF with the structures in the HI emission maps (see Sect. 4.3 in TG17), underlying the potential for TE correlation and E–B power asymmetry.

We can use the proxy ${\text{N}}_{\text{HI}}{}^i$ templates as pure E-mode polarization maps, transforming them to Q and U maps and then use these latter to compute maps of the polarization angle describing the orientation. Given the evidence (Sect. 1), as in TG17, we assume alignment between the local B_POS and dust structures represented by the CNM template and model the polarization angles ψ^c as above. However, for the UNM and WNM templates we assume no such alignment⁷ and so we apply a different procedure in which the polarization angles ψ^u and ψ^w follow from random Gaussian realisations of the components of B_turb (Planck Collaboration Int. XLIV 2016). We simulate the cos² γ factor in Eqs. (5) along each LOS for each of three templates using the relation (Planck Collaboration Int. XLIV 2016) $${\cos }^2\gamma =1-|\widehat{B}.\widehat{r}{|}^2,$$(6)

where $\widehat{B}$ is the unit vector of total GMF and $\widehat{r}$ is the unit vector for a given LOS.

In summary, the seven parameters in the dust model are as follows:

• ɛ₃₅₃ is the mean dust emissivity at 353 GHz (all model Stokes parameters scale with this value);

• p₀ is a polarization parameter that combines the polarization degree of interstellar dust grains and their alignment efficiency;

• galactic coordinates l₀ and b₀ determine the direction of the ordered component of the GMF;

• $f_{\text{M}}^{\text{c}}$ parameterises the relative strength of |B_turb| and |B_ord| in the CNM template. The corresponding parameter is taken to be the same in the UNM and WNM templates and is denoted $f_{\text{M}}^{\text{u/w}}$;

• α_M is the exponent of the power spectrum of the turbulent component of the GMF, assumed to be the same in each template.

We simulate a set of 100 Monte-Carlo model realisations of Stokes I_m, Q_m, and U_m maps at 353 GHz for a set of parameters. Because the Planck data have noise, we add two independent half-mission end-to-end noise realisations provided by Planck (Planck Collaboration III 2020) to each noise-free dust model polarization map to produce two maps (Q and U) with independent noise. We do not add noise to the I_m maps because the I_G353 map that we are using has negligible noise at 353 GHz.

6 Method of constraining model parameters

The best-fit orientation of B_ord is l₀ = 62° and b₀ = −18° (see Table 1 in Sect. 4.2). Here, we describe how we fit the other five parameters of the model.

6.1 Turbulent magnetic field in the CNM template

From the ${\text{N}}_{\text{HI}}{}^{\text{c}}$ template wecompute the polarization angles using the method described in Sect. 5. In particular, we compute maps of cos 2ψ^c. In the CNM the assumed alignment of the local GMF with respect to the gas structures constrains the pair of parameters, $f_{\text{M}}^{\text{c}}$ and α_M, that characterise B_turb.

To evaluate these constraints, following TG17 we simulate cos 2ψ for different values of the pairs assuming the adopted orientation of B_ord, using algebra described in Sect. 4.1 of Planck Collaboration Int. XLIV (2016). For each pair, we compute means of power spectrum amplitudes of cos2ψ from 100 Monte-Carlo realisations over mask65 and fit these with a power-law model $D_{\mathcal{l}}\propto {\mathcal{l}}^{{\alpha }_{\text{M}}+2}$ within multipole range 40 ≤ ℓ ≤ 160 (Fig. 7). Comparing the power spectra of the cos2ψ^c map with the simulated spectra of cos2ψ for different pairs of f_M and α_M, we adopt the following pair of parameter values for the CNM template: $f_{\text{M}}^{\text{c}}={0.5}_{-0.04}^{+0.05}$ and α_M = −2.4 ± 0.1.

Applying the same procedure to the UNM and WNM templates results in a value of $f_{\text{M}}^{\text{u/w}}$ much higher than $f_{\text{M}}^{\text{c}}$ for the CNM template. However, as seen next, the model fit does not support such high turbulence in the UNM or WNM. As discussed in Sect. 6.2, the model is in favour of a very low turbulence in both the UNM and WNM, with $f_{\text{M}}^{\text{u/w}}=0.1$.

6.2 Additional model parameters

We determine the three remaining model parameters, ɛ₃₅₃, p₀, and $f_{\text{M}}^{\text{u/w}}$, by jointly minimising the following two expressions for χ². One is in harmonic space (Fig. 6, Table 2): $${\chi }_{XX}^2={\sum }_{{\mathcal{l}}_{\text{min}}}^{{\mathcal{l}}_{\text{max}}}{\left[\frac{D_{\mathcal{l}}^{XX}-M_{\mathcal{l}}^{XX}({\varepsilon }_{353},p_0,f_{\text{M}}^{\text{u/w}})}{{\sigma }_{\mathcal{l}}^{XX}}\right]}^2,$$(7)

where XX = {EE, BB, TE}, $D_{\mathcal{l}}^{XX}$ and $M_{\mathcal{l}}^{XX}$ are binned power spectra of dust within the multipole range 40 ≤ ℓ ≤ 160 for the data and model, respectively, and ${\sigma }_{\mathcal{l}}^{XX}$ are the corresponding standard deviations. The other is in pixel space (Fig. 8): $${\chi }^2={\sum }_{i=1}^{N_{\text{pix}}}{\left[I_{\text{G353}}-s{I}_{\text{m}}({\varepsilon }_{353},p_0,f_{\text{M}}^{\text{u/w}})-o\right]}^2.$$(8)

Through combining statistics of the polarized dust observations in both harmonic and pixel space, we break the degeneracy between the two model parameters ɛ₃₅₃ and p₀ in Eq. (5). The χ² of the T−T correlation between I_G353 and I_m is minimised in pixel space over mask65. To match the observed dust amplitude at 353 GHz, the value of s should be close to 1 and is kept so by adjusting ɛ₃₅₃ during the iterative solution of the two equations.

The best-fit values are ${\varepsilon }_{353}=146\mu {\text{K}}_{\text{CMB}}{({10}^{20}{\text{cm}}^{-2})}^{-1}$, p₀ = 15.5%, and $f_{\text{M}}^{\text{u/w}}=0.1$. Compared to the values of ɛ₃₅₃ and p₀ found by TG17 in their analysis of the southern Galactic cap, our values differ by factors of 1.2 and 0.84, respectively.

Fig. 7 Means of power spectrum amplitudes from simulated model maps of cos2ψ within the multipole range 40 ≤ ℓ ≤ 160 for different values of fM and αM and the best-fit power laws (dashed lines). The error bars (1 σ) are standard deviations computed from 100 realisations. Black circles are power spectrum amplitudes from the map of cos 2ψc computed from the CNM template.

7 Comparison of the model with Planck observations

Using the methods detailed in Sect. 6, we determined parameters of the Stokes dust model such that several statistical properties computed from the model match the Planck dust observations at 353 GHz over mask65 in both pixel and harmonic space. The quality of these matches to constraining the Planck polarization data is demonstrated below. Furthermore, we show in Sect. 8 that the dust model is able to reproduce the observed inverse relationship between the polarization angle dispersion and the polarization fraction over mask65, even though these statistics were not used in determining the model parameters.

7.1 Polarization fraction

We compute the mean square of the polarization fraction, $p_{\text{m}}^2$, from 100 realisations using Eq. (5). The error bars on the mean histogram are the standard deviations estimated over the 100 realisations. As shown in Fig. 4 (left), our model reproduces the statistical distribution of $p_{\text{d}}^2$ quite accurately, including the negative values of $p_d^2$ and the largest values beyond $p_0^2$.

7.2 Polarization angle

We computethe polarization angle, ${\psi }_{\text{m}}^{\text{R}}$, from 100 realisations of the model, following the method described in Sect. 4.2 and used for the data. For each sky realisation i, we fit the ordered component of GMF using $l_0^i$ and $b_0^i$ and from the map of the mean polarization angle of the GMF ${\psi }_0^i$. At each sky pixel in the realisation of the Stokes model maps we then rotate Q_i and U_i with respect to the new reference angle ${\psi }_0^i$ at respective pixels using Eq. (4). The polarization angle ${\psi }_i^{\text{R}}$ is then derived from the rotated Q_i and U_i maps. The mean polarization angle, ${\psi }_{\text{m}}^{\text{R}}$, for the model is computed by taking the mean of the histograms of ${\psi }_i^{\text{R}}$ distributions for each bin. Similarly, the uncertainties for each bin of the histogram are computed from the standard deviation of 100 realisations. As shown in Fig. 4 (right), our model reproduces the statistical distribution of ${\psi }_{\text{d}}^{\text{R}}$. The dispersion of ${\psi }_{\text{m}}^{\text{R}}$ derived from a Gaussian fit is 17. °5 ± 0. °4, which is slightly higher than that for the data ${\psi }_{\text{d}}^{\text{R}}$: 16. °6 ± 0. °2. The difference might come from the fact that the model is not fitted to the Planck data at ℓ < 40, which cannot be accomplished accurately over the limited area of mask65.

Fig. 8 Correlation plot between Im and IG353 at 353 GHz. The black dashed line is the result of the joint optimisation of Eqs. (7) and (8).

7.3 Dust polarization spectra

We compute mean cross-power spectra from 100 sets of two independent realisations of the Stokes model maps. We fit the binned power spectra with power laws: $M_{\mathcal{l}}^{XX}=A_{XX}{(\mathcal{l}/80)}^{{\alpha }_{XX}+2}$ within the multipole range 40 ≤ ℓ ≤ 160. Error bars in each bin are the standard deviations computed from 100 realisations. Best-fit amplitudes, A_XX and exponents, α_XX along with the corresponding χ² values, are listed in the right column of Table 2. We compare the results with the Planck dust power spectra in Fig. 6. The $D_{\mathcal{l}}^{EE}$ and $D_{\mathcal{l}}^{BB}$ amplitudes and their ratio match the values from the data accurately. The $D_{\mathcal{l}}^{TE}$ amplitude for the dust model is slightly lower, but consistent with the Planck data within its uncertainty.

7.4 Total intensity

Figure 8 shows the tight correlation between I_m and I_G353 at 353 GHz, which demonstrates that statistically, the model dust intensity reproduces the observed Planck dust intensity. The best-fit dashed line shown has slope s = 1.0 and offset o = −35 μK_CMB, again within the uncertainty of the zero level in the intensity, 40 μK_CMB (Sect. 2.1).

8 Inverse relationship between the polarization angle dispersion and the polarization fraction

The polarization angle dispersion function, $\mathcal{S}$, was introduced in Planck Collaboration Int. XIX (2015) to quantify the local dispersion of the dust polarization angle for a given angular resolution and lag δ. An inverse relationship between $\mathcal{S}$ and p_MAS was found, where p_MAS is the modified asymptotic estimator of the polarization fraction (Plaszczynski et al. 2014). That analysis was done over 42% of the sky at low and intermediate Galactic latitudes at a resolution of 1° (FWHM) with a lag δ = 30′. Planck Collaboration XII (2020) extended this work over the full sky using GNILC processed intensity and polarization maps at a resolution of 160′ (FWHM) with δ = 80′ and found that $\mathcal{S}\propto p_{\text{MAS}}{}^{-1}$.

Here, westudy the relationship between $\mathcal{S}$ and p_MAS over mask65 using 353 GHz maps at a resolution of 80′ (FWHM) with δ = 40′. The noise bias parameter needed for p_MAS (Montier et al. 2015a,b) is estimated from the smoothed noise covariance matrices, σ_II, σ_IQ, σ_IU, σ_QQ, σ_QU, σ_UU at a resolution of 80′ (FWHM) and N_side = 128.

In the top panel of Fig. 9, we present the two-dimensional joint distribution of $\mathcal{S}$ and p_MAS for the Planck data over mask65. The running mean of $\mathcal{S}$ in each bin of p_MAS follows an inverse relationship, $\mathcal{S}=$(0. °43 ± 0. °04)∕p_MAS. In the bottom panel of Fig. 9, we show the corresponding $\mathcal{S}$–p_MAS distribution for our model. This too has an inverse relationship and the slope 0. °50 ± 0. °05 is very close to that for the Planck data. This is remarkable given that this statistical property of dust polarization has not been exploited in the fit of our model parameters.

Fig. 9 Top panel: two-dimensional histogram of the joint distribution of $\mathcal{S}$ and pMAS over mask65 for the Planck data at a resolution of 80′ with lag δ = 40′. The black circles show the mean of $\mathcal{S}$ in bins of pMAS containing the same number of map pixels. Error bars are the standard deviation of $\mathcal{S}$ in each bin. The black dotted line is a fit of the running mean. Bottom panel: same as top panel, but for the model.

9 Summary and discussion

In this paper, we analyse the statistical properties of the Planck dust polarization maps at 353 GHz at a 60′ resolution (FWHM) over a low-column-density region that accounts for 65% of the northern Galactic cap at latitudes larger than 30° (mask65). We make use of the multiphase dust polarization model described in TG17. The model shows how the dust polarization across the sky can be approximated on the basis of proxy HI data: in particular, three independent templates notionally representingthe contributions from the CNM, UNM, and WNM. The model has seven adjustable parameters whose values are determined by reproducing Planck observations at 353 GHz: in particular, the one-point statistics of p and ψ in pixel space, and the EE, BB, and TE power spectra in harmonic space. Our main results can be summarised as follows.

• The butterfly pattern seen in q_d and u_d maps (Fig. 5) around the pole is associated with an ordered GMF, which we find has a mean orientation in Galactic coordinates toward (l₀, b₀) = (62°, −18°) (Sect. 4.2). The best-fit value of l₀ is roughly consistent with the earlier values derived from starlight polarization (Ellis & Axon 1978; Heiles 1996).

• From the Planck data, we find a BB∕EE power ratio of 0.58 and significant positive TE correlation over mask65 (Sect. 4.3). The key property in our model that allows the model to reproduce these observations (Sect. 7.3) is the alignment of the local magnetic field B_POS and the HI structure in the CNM template (Sect. 5).

• The observed distributions of p and ψ over mask65 show a scatter that is comparable to that of the distributions reported for the whole sky in Planck Collaboration XII (2020). Our model successfully reproduced these one-point statistics associated with LOS depolarization (Sects. 7.1 and 7.2) by introducing fluctuations in the GMF orientation that are uncorrelated between the three independent HI templates (Sect. 3.2).

• The best-fit value of the parameter, p₀, that measures the grain alignment efficiency combined with the intrinsic polarization fraction of interstellar dust emission at 353 GHz, is 15.5%. The best-fit value of the dust emissivity is ${\varepsilon }_{353}=146\mu {\text{K}}_{\text{CMB}}{({10}^{20}{\text{cm}}^{-2})}^{-1}$ (Sect. 6.2).

• To match the observed EE, BB, and TE amplitudes, the model fit yields a low value of the parameter $f_{\text{M}}^{\text{u/w}}$ specifying the relative amplitude of the turbulent component of the GMF in the UNM and WNM. This value, 0.1, is significantly smaller than the value $f_{\text{M}}^{\text{c}}=0.5$ characterising the CNM (Sect. 6.1).

• The spectral index of the turbulent component of the GMF, α_M = −2.4, fitted over mask65 (Sect. 6.1), is consistent with the value reported in TG17. These two complementary analyses reveal that the spectral index of turbulent magnetic field closely matches the power-law index of dust polarization power spectra, as expected.

• Our model also reproduces the inverse relationship between the polarization angle dispersion, $\mathcal{S}$, and the polarization fraction, p_MAS, present in the Planck data, despite the fact that we do not utilise this phenomenon in fitting the model parameters (Sect. 8). Our work reinforces the conclusion in Planck Collaboration XII (2020) that the inverse relationship of $\mathcal{S}$ and p_MAS is a generic feature associated with the GMF structure.

This work demonstrates that the phenomenological model introduced in TG17 for the southern Galactic cap also works well in the northern hemisphere. The next step in our modelling work will be to extend this framework to multiple frequencies by incorporating spectral energy distributions for the dust emission associated with the three HI templates thus introducing the potential for frequency decorrelation of the dust polarization. This is a necessary step towards investigating the utility of this framework for evaluating component separation methods for future CMB missions.

Acknowledgements

We gratefully acknowledge the use of the Aquila cluster at NISER, Bhubaneswar. D.A. acknowledges the University Grants Commission India for providing financial support as Senior Research Fellow. This work was supported by the Science and Engineering Research Board, Department of Science and Technology, Govt. of India grant number SERB/ECR/2018/000826 and the Natural Sciences and Engineering Research Council of Canada. Some of the results in this paper have been derived using the HEALPix package. The Planck Legacy Archive (PLA) contains all public products originating from the Planck mission, and we take the opportunity to thank ESA/Planck and the Planck Collaboration for the same. This work has made use of HI data of the EBHIS survey headed by the Argelander-Institut für Astronomie (PI: J. Kerp) in collaboration with the Max-Planck-Institut für Radioastronomie and funded by the Deutsche Forschungsgemeinschaft (grants KE757/7-1 to 7-3). The Gaussian decomposition of the EBHIS data was supported by the Estonian Research Council grant IUT26-2, and by the European Regional Development Fund (TK133).

Appendix A Analysis with an alternative set of HI templates

Here we describe the development of an alternative set of HI templates and then the results from a modelling analysis similar to that presented in the main paper.

Kalberla & Haud (2018; hereafter KH18) performed a Gaussian decomposition of EBHIS HI emission spectra from an intermediate product of the HI4PI all-sky survey (HI4PI Collaboration 2016), applying the algorithm described in Haud (2000) with some modifications (see their Sect. 2.2). Based on line widths, KH18 defined HI components called the cold, luke-warm, and warm neutral medium (CNM, LNM, and WNM, respectively). The separation was done using intervals of line width or the corresponding Doppler temperature (T_D), which depend on |v| (see, e.g. their Figs. 3 and 7). For example, for velocities |v| ≤ 8 km s⁻¹ over the full sky, the mean line widths are 3.6, 9.6, and 23.3 km s⁻¹ FWHM for the CNM, LNM, and WNM components, respectively.

Making use of this Gaussian decomposition, we consider here Gaussian components whose central peak velocity $v_{\text{c}}^i$ lies within the velocity range of |v_c| ≤ 24 km s⁻¹, which encompasses most of the local HI gas (Kalberla & Haud 2018) but ignores most intermediate velocity gas (IVC). This decomposition was done at the native resolution in maps with N_side = 1024. As in Kalberla & Haud (2018), maps of the integrated emission of HI were made for the CNM, LNM, and WNM components at an angular resolution of 30′ (FWHM). For our alternative analysis, these column density maps were further smoothed to 60′ (FWHM) and downgraded to N_side = 128. The smoothing reduces uncertainties within each map that arise from separating the gas into components. These maps are shown in Fig. A.1. Hereafter, we refer to these as the KH18 templates.

We refer to the sum of the column densities of the three templates as N_HI,24. The correlation of dust emission with N_HI,24 on the selected region mask65 is not satisfactory, having an emissivity ɛ_353,24 of $121\mu {\text{K}}_{\text{CMB}}{({10}^{20}{\text{cm}}^{-2})}^{-1}$, lower than ɛ_353,50, and also an unacceptably large offset, 115 μK_CMB. This arises because of the restricted velocity range adopted from KH18: N_HI,24 ≤ N_HI,50 and is not an adequate proxy for all of the dust. Therefore, to make the model more rigorous we introduced an additional fourth template to compensate for the missing HI. The fourth template is simply the column density map N_HI,50 –N_HI,24, also shown in Fig. A.1.

The sum of the KH18 CNM and LNM templates is related to our CNM template used in the main paper, as shown in Fig. A.2. The slope of the pixel by pixel correlation between the sum of the KH18 CNM and LNM templates and our CNM template is 1.02 over mask65, which confirms that CNM and LNM from KH18 can be considered embedded in the gas in our CNM template. Therefore, in building this alternative model we assume that the local B_POS of the model GMF is aligned with HI structures in both of the CNM and LNM templates of KH18 (Sect. 5).

Applying the technique discussed in Sect. 6.1, we find $f_{\text{M}}^{\text{c}}=0.5$ and α_M = −2.4 (see Fig. A.3). Because the error bars on cos2ψ estimated from simulated maps are significant, on the basis of Fig. A.3 we adopt the same pair of parameter values for the LNM. In the WNM, we adopt $f_{\text{M}}^{\text{w}}$ a lower value, 0.1. The fourth template is assumed to be unpolarized. For these alternative templates, the other three parameters optimised (Sect. 6.2) are ${\varepsilon }_{353}=146\mu {\text{K}}_{\text{CMB}}{({10}^{20}{\text{cm}}^{-2})}^{-1}$ and p₀ = 19%.

Fig. A.1 Northern Galactic orthographic projection of NHI of the three KH18 templates, CNM (upper left), LNM (upper right), and WNM (lower left), and the additional fourth template (lower right) over our selected region mask65 (Sect. 3.1). The mean column densities of the four templates are $\langle N_{\text{HI}}^{\text{c}}\rangle =0.19\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{l}}\rangle =0.52\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{w}}\rangle =0.75\times {10}^{20}{\text{cm}}^{-2}$, and $\langle N_{\text{HI}}^{\text{4th}}\rangle =0.37\times {10}^{20}{\text{cm}}^{-2}$, respectively.

Within our multiphase dust modelling framework (Sect. 5), these seven parameters and four templates can reproduce (Sect. 7) all of the observed statistical properties of dust (Sect. 4) quite well, including the power spectra as shown in Fig. A.4. For ℓ bins between 40 and 160, the alternative set of HI templates is able to produce dust EE and BB power spectra with similar strength as the Planck data. The overall BB power is at the right level and the BB∕EE ratio is about 0.57.

However, the predicted TE is an exception. There are two concerns. First, the overall level of the TE power across most scales is low relative to the Planck data (shown originally in Fig. 6 left). Second, the pattern of ℓ-dependent deviations from a smooth power-law dependence is more pronounced compared to our model in Fig. 6. These deviations could come from excess small-scale power in the KH18 CNM or LNM templates that corrupts the TE (and EE) power spectra at small scales.

Fig. A.2 Northern Galactic orthographic projection of our CNM template (left panel), the KH18 CNM plus LNM templates (middle panel), and the difference between the two (right panel).

Fig. A.3 Mean of simulated power spectra of cos2ψ within the multipole range 40 ≤ ℓ ≤ 160 for different values of fM and αM and best fit power-law (dashed lines). The error-bars (1 σ) are standard deviations computed from 100 realisations.Black circles and black diamonds are power spectra of cos 2ψ computed (Sect. 3.2) from the CNM and LNM templates of KH18.

Fig. A.4 Dust EE, BB, and TE cross-power spectra in units of $\mu {\text{K}}_{\text{CMB}}{}^2$ computed from the subsets of 100 realisations of the dust model maps over mask65 using the CNM, LNM, and WNM of KH18 and the fourth template. The filled areas represent the Planck dust power spectra measurements over mask65 computed in Sect. 4.3.

References

All Tables

All Figures

	Fig. 1 Histogram of residuals at different iterations. The blue curve is the Gaussian fitted to residuals at the final iteration. Gray shadedregions are 1 σcg, 2 σcg, and 3 σcg regions of the fitted Gaussian.
In the text

	Fig. 2 Binary mask defining mask65, the dark regions selecting 65% of the northern Galactic cap (b ≥ 30°).
In the text

Fig. 3 Northern Galactic orthographic projection of integrated HI emission maps of the three HI templates called CNM (left panel), UNM (middle panel), and WNM (right panel). Over our selected region mask65 (Sect. 3.1), the mean column densities of the three templates are $\langle N_{\text{HI}}^{\text{c}}\rangle =0.90\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{u}}\rangle =0.32\times {10}^{20}{\text{cm}}^{-2}$, and $\langle N_{\text{HI}}^{\text{w}}\rangle =0.56\times {10}^{20}{\text{cm}}^{-2}$, respectively.

In the text

Fig. 4 Left panel: histogram of the square of the polarization fraction, p2, at 353 GHz. Right panel: histogram of ψR, the difference between the measured polarization angle at 353 GHz and that computed for the mean GMF orientation over mask65. Green triangles are mean $p_{\text{d}}^2$ (${\psi }_{\text{d}}^{\text{R}}$) estimated from Planck data subsets, and black squares are mean $p_{\text{m}}^2$ (${\psi }_{\text{m}}^{\text{R}}$) estimated from 100 realisations of model maps. The vertical dashed line in the left panel corresponds to the model parameter p0 = 15.5% ($p_0^2=0.024$), while the vertical dashed line in the right panel corresponds to a model without a turbulent component of the GMF.

In the text

	Fig. 5 Upper panel: northern Galactic orthographic projection of qd (left) and ud (right) for the Planck data. Lower panel: same for the best-fit model where the large-scale GMF is directed toward (l0, b0) = (62°, −18°) and scaled with a mean polarization fraction of 8.8%.
In the text

Fig. 6 Left panel: dust EE, BB, and TE cross-power spectra computed from the subsets of Planck data at 353 GHz over mask65. Dashed lines represent the best-fit power laws. Right panel: similar plots as left panel computed from 100 realisations of the dust model maps. Error bars are 1 σ uncertainties as explained in the main text. The filled areas represent the Planck dust power spectra measurements over mask65.

In the text

	Fig. 7 Means of power spectrum amplitudes from simulated model maps of cos2ψ within the multipole range 40 ≤ ℓ ≤ 160 for different values of fM and αM and the best-fit power laws (dashed lines). The error bars (1 σ) are standard deviations computed from 100 realisations. Black circles are power spectrum amplitudes from the map of cos 2ψc computed from the CNM template.
In the text

	Fig. 8 Correlation plot between Im and IG353 at 353 GHz. The black dashed line is the result of the joint optimisation of Eqs. (7) and (8).
In the text

Fig. 9 Top panel: two-dimensional histogram of the joint distribution of $\mathcal{S}$ and pMAS over mask65 for the Planck data at a resolution of 80′ with lag δ = 40′. The black circles show the mean of $\mathcal{S}$ in bins of pMAS containing the same number of map pixels. Error bars are the standard deviation of $\mathcal{S}$ in each bin. The black dotted line is a fit of the running mean. Bottom panel: same as top panel, but for the model.

In the text

Fig. A.1 Northern Galactic orthographic projection of NHI of the three KH18 templates, CNM (upper left), LNM (upper right), and WNM (lower left), and the additional fourth template (lower right) over our selected region mask65 (Sect. 3.1). The mean column densities of the four templates are $\langle N_{\text{HI}}^{\text{c}}\rangle =0.19\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{l}}\rangle =0.52\times {10}^{20}{\text{cm}}^{-2}$, $\langle N_{\text{HI}}^{\text{w}}\rangle =0.75\times {10}^{20}{\text{cm}}^{-2}$, and $\langle N_{\text{HI}}^{\text{4th}}\rangle =0.37\times {10}^{20}{\text{cm}}^{-2}$, respectively.

In the text

	Fig. A.2 Northern Galactic orthographic projection of our CNM template (left panel), the KH18 CNM plus LNM templates (middle panel), and the difference between the two (right panel).
In the text

	Fig. A.3 Mean of simulated power spectra of cos2ψ within the multipole range 40 ≤ ℓ ≤ 160 for different values of fM and αM and best fit power-law (dashed lines). The error-bars (1 σ) are standard deviations computed from 100 realisations.Black circles and black diamonds are power spectra of cos 2ψ computed (Sect. 3.2) from the CNM and LNM templates of KH18.
In the text

	Fig. A.4 Dust EE, BB, and TE cross-power spectra in units of $\mu {\text{K}}_{\text{CMB}}{}^2$ computed from the subsets of 100 realisations of the dust model maps over mask65 using the CNM, LNM, and WNM of KH18 and the fourth template. The filled areas represent the Planck dust power spectra measurements over mask65 computed in Sect. 4.3.
In the text