Dust polarization modeling at large-scale over the Northern Galactic cap using EBHIS and Planck data

TL;DR

This paper develops a detailed dust polarization model using Planck and HI survey data to improve the separation of Galactic dust emission from the CMB B-mode polarization signal, aiding future cosmological measurements.

Contribution

It introduces a seven-parameter dust model that accurately reproduces observed dust polarization statistics over the northern Galactic cap.

Findings

Model reproduces observed polarization statistics

Uses HI data as a proxy for dust intensity

Assists in separating Galactic dust from CMB signals

Abstract

The primary source of systematic uncertainty in the quest for the B-mode polarization of the Cosmic Microwave Background (CMB) introduced by primordial gravitational waves is polarized thermal emission from Galactic dust. Therefore, accurate characterization and separation of the polarized thermal dust emission is an essential step in distinguishing such a faint CMB B-mode signal. We provide a modelling framework to simulate polarized thermal dust emission based on the model described in Ghosh et al. (2017), making use of both the Planck dust and Effelsberg-Bonn HI surveys over the northern Galactic cap. Our seven-parameter dust model, incorporating both HI gas in three different column density templates as a proxy for spatially variable dust intensity and a phenomenological model of Galactic magnetic field, is able to reproduce both 1- and 2-point statistics of the observed dust polarization maps seen by Planck at 353 GHz over a selected low-column density region in the northern Galactic cap. This work has important applications in assessing the accuracy of component separation methods and in quantifying the confidence level of separating polarized Galactic emission and the CMB B-mode signal, as is needed for ongoing and future CMB missions.

Tables

Table 1. Table 1: Mean polarization fraction at 353 GHz and orientation of the large-scale GMF using different subsets of the Planck  data.

Dataset mean polarization subsets large-scale GMF fraction ($p_d$) $l_0$ $b_0$ Half Missions 8.8 % HM1 61∘  $-{17}^{\circ }$ HM2 64∘  $-{19}^{\circ }$ Years 8.8 % YR1 61∘  $-{17}^{\circ }$ YR2 64∘  $-{19}^{\circ }$ Odd-Even 8.8 % odd 63∘  $-{18}^{\circ }$ even 62∘  $-{18}^{\circ }$ Mean – – 62∘  $-{18}^{\circ }$

Table 2. Table 2: Fitted dust power spectra of the Planck  data at 353 GHz and of the dust model over mask65.

Parameter      Planck 353 GHz data      Dust model ${\alpha }_{EE}$      $-2.34\pm 0.16$      $-2.56\pm 0.18$ ${\alpha }_{BB}$      $-2.51\pm 0.21$      $-2.56\pm 0.25$ ${\alpha }_{TE}$      $-2.43\pm 0.38$      $-2.75\pm 0.16$ ${\chi }_{EE}^2\left(N_{d.o.f.}=4\right)$      2.2      1.3 ${\chi }_{BB}^2\left(N_{d.o.f.}=4\right)$      1.8      1.3 ${\chi }_{TE}^2\left(N_{d.o.f.}=4\right)$      2.3      8.3 $A_{EE}$ [$\mu K_{\mathrm{CMB}}^2$]      $39.1\pm 2.2$      $39.9\pm 2.7$ $A_{BB}$ [$\mu K_{\mathrm{CMB}}^2$]      $22.5\pm 1.6$      $22.0\pm 1.4$ $A_{TE}$ [$\mu K_{\mathrm{CMB}}^2$]      $74.2\pm 10.9$      $65.6\pm 4.0$ $⟨A_{BB}/A_{EE}⟩$      $0.58\pm 0.05$      $0.57\pm 0.06$ $⟨A_{TE}/A_{EE}⟩$      $1.90\pm 0.30$      $1.85\pm 0.17$

Equations

$$N_{\text{mathsc}Hi}$$

$$T_{\mathrm{b}}(v)=\sum _i{T}_0^i\exp \left[-\frac{1}{2}{\left(\frac{v_{LSR}-v_{\mathrm{c}}^i}{{\sigma }^i}\right)}^2\right],$$

$$p_{\mathrm{d}}^2=⟨\frac{Q_{\mathrm{d}}^{{\mathrm{s}}_1}{Q}_{\mathrm{d}}^{{\mathrm{s}}_2}+U_{\mathrm{d}}^{{\mathrm{s}}_1}{U}_{\mathrm{d}}^{{\mathrm{s}}_2}}{I_{G353}^2}⟩,$$

$$\displaystyle\left[\begin{array}[]{c}Q_{\rm d}^{\rm R}\\ U_{\rm d}^{\rm 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_{\rm d}\\ U_{\rm d}\end{array}\right]\ ,$$

$$I_{\mathrm{m}}(\widehat{\text{@vec}n})$$

$$Q_{\mathrm{m}}(\widehat{\text{@vec}n})$$

$$U_{\mathrm{m}}(\widehat{\text{@vec}n})$$

$${\cos }^2\gamma =1-|\widehat{\text{@vec}B}.\widehat{\text{@vec}r}{|}^2,$$

$${\chi }_{XX}^2=\sum _{{\ell }_{\min }}^{{\ell }_{\max }}{\left[\frac{{\mathcal{D}}_{\ell }^{XX}-{\mathcal{M}}_{\ell }^{XX}({ϵ}_{353},p_0,f_{\mathrm{M}}^{u/w})}{{\sigma }_{\ell }^{XX}}\right]}^2,$$

$${\chi }^2=\sum _{i=1}^{N_{\mathrm{pix}}}{\left[I_{G353}-s{I}_{\mathrm{m}}({ϵ}_{353},p_0,f_{\mathrm{M}}^{u/w})-o\right]}^2.$$

Full text

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 2018), also produce a $B$-mode signal at microwave frequencies (Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2018). 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 ($\@vec{B}_{\rm 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 $\@vec{B}_{\rm 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, $\psi$ (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 $\psi$ 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 (2018).

The Planck maps were also used to measure power spectra of dust polarization. At intermediate and high Galactic latitudes, in the multipole range, $40<\ell<600$, the ratio of dust $EE$ to $BB$ power (${\cal D}_{\ell}^{EE}/{\cal D}_{\ell}^{BB}$) is found to be around 2 (Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2018). 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) and Kandel et al. (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 $\mathsc{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 $\mathsc{Hi}$ emission features along the LOS, later used by Clark & Hensley (2019) to model dust polarization. Clearly, $\mathsc{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 $\mathsc{Hi}$ emission data into three distinct maps of $\mathsc{Hi}$ column density referred to as ‘$\mathsc{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\leq-30^{\circ}$) where dust and $\mathsc{Hi}$ emission are well correlated. In this paper, we use this framework to fit the Planck PR3 data (Planck Collaboration I 2018) over a fraction of the northern Galactic cap at $b\geq 30^{\circ}$. 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 $\mathsc{Hi}$ Survey (EBHIS). Section 3 describes the choice of sky region and $\mathsc{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 $\mathsc{Hi}$ templates.

2 Datasets used

2.1 Planck dust polarization maps

In this paper, we use publicly available Planck PR3 data111http://www.cosmos.esa.int/web/planck/pla (Planck Collaboration I 2018) 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 (${\rm K}_{\rm CMB}$, Planck Collaboration III 2018). 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 2018). 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 ${\mu\rm K}_{\rm CMB}$  at 353 GHz, which needs to be subtracted (Planck Collaboration III 2018). To correct for the Galactic $\mathsc{Hi}$ offset and the contribution from dust emission associated with $\mathsc{Hii}$ emission (Planck Collaboration XII 2018), a ‘fiducial’ offset of 63 ${\mu\rm K}_{\rm CMB}$ is added back to the CIB-subtracted $\tt GNILC$ dust intensity map (hereafter, we refer to this final map as $I_{G353}$). The fiducial offset has an uncertainty of 40 ${\mu\rm K}_{\rm CMB}$, with an associated ‘low’ and ‘high’ offset of 23 ${\mu\rm K}_{\rm CMB}$ and 103 ${\mu\rm K}_{\rm CMB}$, respectively.

The Planck Stokes $Q$ and $U$ maps at 353 GHz have a beam resolution of 4$\aas@@fstack{\prime}$82 (FWHM), and the $I_{G353}$ map has the beam resolution of 5*′* (FWHM). To increase the signal-to-noise ratio, we smooth the dust Stokes $I$, $Q$, and $U$ maps to a common resolution of 60*′* or 80*′* (FWHM) and reproject on the $\tt HEALPix$ grid222https://healpix.jpl.nasa.gov (Górski et al. 2005) with $N_{\rm 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 $\mathsc{Hi}$ data

We use $\mathsc{Hi}$ 21-cm line spectroscopic data from EBHIS333EBHIS data are available at http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/585/A41 (Kerp et al. 2011; Winkel, B. 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$\aas@@fstack{\prime}$8 (FWHM), spectral resolution $\delta v=1.44$  km s*-1*(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 $\tt HEALPix$ grid with $N_{\rm side}=1024$ (${\theta}_{\rm pix}$ = 3$\aas@@fstack{\prime}$4). The EBHIS data also form the northern part of the all-sky product HI4PI (HI4PI Collaboration et al. 2016).

For optically thin emission, the total $\mathsc{Hi}$ column density ($N_{\mathsc{Hi}}$) can be obtained by integrating the brightness temperature ($T_{\rm b}$) over velocity channels (Wilson et al. 2009):

[TABLE]

3 Exploiting the $\mathsc{Hi}$  data

3.1 Region selection

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

Ultraviolet observations of molecular hydrogen, H2 (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 H2 becomes significant for sight lines where the total $N_{\mathsc{H}}$ exceeds $4\times 10^{20}\thinspace$ cm*-2* (Arendt et al. 1998). Therefore, we restrict ourselves to the sky region in which $N_{\mathsc{Hi},50}$ is below a threshold $3.8\times 10^{20}\thinspace$ cm*-2*, where $N_{\mathsc{Hi},50}$ from Eq. (1) is evaluated over the restricted velocity interval $|v|\leq 50$  km s*-1*. We explored extending the range up to $80$  km s*-1* and because there is little additional gas in this range our results below from the modelling analysis (Sects. 6 to 8) are robust. We also get the same model parameters for regions restricted by $N_{\mathsc{Hi},50}$ thresholds in the range $3.4-4.0\times 10^{20}\thinspace$ cm*-2*. The sky fraction increases beneficially with threshold, but not significantly beyond the adopted $3.8\times 10^{20}\thinspace$ cm*-2*. However, the model parameters change significantly for thresholds above $4\times 10^{20}\thinspace$ cm*-2*, 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_{\mathsc{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 ($\sigma_{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\thinspace\sigma_{cg}$. A stable mask is obtained after five iterations, a region of $5900\thinspace\text{deg}^{2}$ comprising 65 % of the northern Galactic cap ($b\geq 30^{\circ}$), 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^{\circ}$ (FWHM). After apodization, we have an effective sky fraction, $f_{\rm sky}^{\rm 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_{\ell}$ from the pseudo-spectrum, these isolated islands introduce correlations from very high $\ell$ to low $\ell$ via the mode-mode coupling matrix (Hivon et al. 2002). In our analysis, we eliminate isolated islands containing less than 20 pixels using the $\tt Process\char 95\relax mask$ routine of $\tt HEALPix$.

Over mask65 the mean $\mathsc{Hi}$ column density over the restricted velocity interval is $\left\langle N_{\mathsc{Hi},50}\right\rangle=1.85\times 10^{20}$ cm*-2* and the mean dust intensity is $\left\langle I_{G353}\right\rangle=293\thinspace{\mu\rm K}_{\rm CMB}$ at 353 GHz. The linear correlation between $I_{G353}$ and total $N_{\mathsc{Hi}}$  yields a slope or emissivity $\epsilon_{353,50}$ of $137\thinspace{\mu\rm K}_{\rm CMB}(10^{20}\thinspace$ cm${}^{-2})^{-1}$, where again the 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\thinspace{\mu\rm K}_{\rm CMB}$ is satisfactory, within the systematic uncertainty of the zero level (Sect. 2.1).

3.2 $\mathsc{Hi}$ templates

The diffuse ISM is a complex turbulent multiphase medium. The $\mathsc{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.

$\mathsc{Hi}$ spectra can be used to build maps of $\mathsc{Hi}$ emission associated notionally with the CNM, UNM, and WNM. First, the observed brightness temperature profiles of $\mathsc{Hi}$ spectra are decomposed into Gaussian components (Haud 2000, 2013) under the assumption that the random velocities in any $\mathsc{Hi}$ ‘cloud’ have a Gaussian distribution. Thus,

[TABLE]

where the summation is over all Gaussian components needed to describe the profile along a LOS, $v_{\rm LSR}$ is the velocity of the $\mathsc{Hi}$ gas, and $T_{0}^{i}$, $v_{\rm c}^{i}$, and ${\sigma}^{i}$ are the peak brightness temperature (expressed in Kelvin), central peak velocity, and 1 $\sigma$ standard deviation of Gaussian component $i$, respectively.

Following the detailed physically motivated prescription in section 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$\aas@@fstack{\prime}$8 (FWHM), which we have smoothed to a common resolution of 60*′*(FWHM) at $N_{\rm side}$= 128. The resulting maps, to be used as $\mathsc{Hi}$ templates below, are shown in Fig. 3. In Appendix. A, we present an analysis that uses alternative $\mathsc{Hi}$ templates based on a decomposition by Kalberla & Haud (2018).

4 Observed statistics from ${\it 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_{\rm d}^{2}$) combining independent subsets of the Planck data at 353 GHz:

[TABLE]

where ‘d’ stands for the Planck data, $(s_{1},s_{2})$ stands for different subsets of Planck data at 353 GHz: \big{\{}(HM1, HM2); (Y1, Y2); (odd, even)\big{\}}, and $\left<...\right>$ implies an average over different data subsets. In the left panel of Fig. 4, we present the histogram of $p_{\rm d}^{2}$. The uncertainties (1 $\sigma$) in $p^{2}_{\rm d}$ 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_{\rm d}$ over mask65 are listed in Table 1 for the different data subsets.

4.2 Polarization angle

We computed the Stokes parameter ratios, $q_{\rm d}=Q_{\rm d}/I_{G353}$ and $u_{\rm d}=U_{\rm d}/I_{G353}$ from Planck data at 353 GHz. As shown in the top panel of Fig. 5, maps of $q_{\rm d}$ and $u_{\rm d}$ centred at the northern Galactic pole reveal a ‘butterfly’ pattern associated with an ordered large-scale GMF.

To determine $l_{0}$ and $b_{0}$ 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_{0}$ and $b_{0}$ for different data subsets are noted in Table 1 and the lower panel of Fig. 5 presents the best-fit model. These values, $(l_{0},b_{0})=(62^{\circ},-18^{\circ})$, give the ‘orientation’ of the GMF, but not the direction; i.e. the $\pi$ 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_{0},b_{0})=(74^{\circ},+24^{\circ})$. The change of sign of $b_{0}$ 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_{0}$ and $b_{0}$ are then used to calculate a polarization angle map, ${\psi}_{0}$, for the mean GMF, where ${\psi}_{0}$ 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 ${\psi}_{0}$ at respective pixels using the relation,

[TABLE]

where $Q_{\rm d}^{\rm R}$ and $U_{\rm d}^{\rm R}$ are Stokes parameters in the rotated frame. The polarization angle $\psi_{\rm d}^{\rm R}$ is derived using the relation $\psi_{\rm d}^{\rm R}=0.5\thinspace\mathrm{atan2}(-U_{\rm d}^{\rm R},Q_{\rm d}^{\rm R})$,

where the two-argument function $\mathrm{atan2}(-U,Q)$ is used in place of $\mathrm{atan}(-U/Q)$ to avoid the $\pi$-ambiguity. The minus sign is introduced to compute $\psi_{\rm d}^{\rm R}$ in the IAU convention (Hamaker & Bregman 1996) from $\tt HEALPix$ maps given in the $\tt COSMO$  convention444https://healpix.jpl.nasa.gov/html/intronode12.htm used by Planck.

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

4.3 Dust power spectra

We analyse the dust polarization power spectra at 353 GHz over mask65,555${\cal D}^{XX}_{\ell}\equiv\ell(\ell+1)C^{XX}_{\ell}/(2\pi$), where $XX\thinspace=EE,BB,TE$. ${\cal D}_{\ell}^{EE}$, ${\cal D}_{\ell}^{BB}$, and ${\cal D}_{\ell}^{TE}$, in the multipole range $40\leq\ell\leq 160$. These were computed using $\tt Xpol$ (Tristram et al. 2005) from cross spectra of three datasets: half-missions (HM1 $\times$ HM2), years (Y1 $\times$ Y2), and odd-even surveys (Odd $\times$Even). In Fig. 6, we present the mean of the three cross-power spectra in six multipole bins. The dust ${\cal D}_{\ell}^{EE}$ and ${\cal D}_{\ell}^{BB}$ spectra are corrected for CMB contributions in harmonic space using the Planck best-fit power spectra (Planck Collaboration VI 2018), whereas we keep ${\cal D}_{\ell}^{TE}$ unaltered. Uncertainties (1 $\sigma$) on the binned dust power spectra are the quadratic sum of statistical noise computed by $\tt 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 $\tt SMICA$ and $\tt SEVEM$ maps (Planck Collaboration IV 2018) 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, ${\cal D}_{\ell}^{XX}=A_{XX}(\ell/80)^{\alpha_{XX}+2}$, where $A_{XX}$ is the best-fit amplitude at $\ell=80$, $\alpha_{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 $\alpha_{XX}$ and the respective values of $\chi^{2}$ 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.

5 Multiphase model of polarized dust emission

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

We use the fact that dust emission and $N_{\mathsc{Hi}}$ are correlated so that templates of $N_{\mathsc{Hi}}$ can be used as proxies for spatially variable dust intensity. For optically thin dust emission, the model Stokes parameters $I_{\rm m}$, $Q_{\rm m}$, and $U_{\rm m}$ at 353 GHz can be written as

[TABLE]

where $N$ is the number of distinct templates , $\hat{\@vec{n}}$ is the direction vector, $p_{0}$ is a parameter related to the dust grain properties (Lee & Draine 1985; Draine & Fraisse 2009; Planck Collaboration Int. XX 2015), $\gamma$ is the angle made by the local magnetic field with the plane of the sky, $\psi$ is the polarization angle measured counterclockwise from Galactic north,666IAU convention, the minus sign producing $U$ in the $\tt COSMO$ convention. $\epsilon_{353}$ is the dust emissivity at 353 GHz for each $\mathsc{Hi}$ template, and $N_{\mathsc{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 by a polarization model that is inherently minimalist and not dependent on fine-tuning. For simplicity, $\epsilon_{353}$ 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 $\gamma$ and $\psi$ is developed as in Planck Collaboration XLIV (2016). Following Jaffe et al. (2010), the GMF, $\@vec{B}$, is expressed as a vector sum of an ordered component ($\@vec{B}_{\rm ord}$) and a fluctuating component ($\@vec{B}_{\rm 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, $\@vec{B}_{\rm 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 $\@vec{B}_{\rm ord}$, oriented toward the $l_{0}$ and $b_{0}$ that are already determined. However, the components of $\@vec{B}_{\rm turb}$ are different in each template and are taken from independent realisations of a Gaussian random field that has a power spectrum with $C_{\ell}\propto\ell^{\alpha_{\rm M}}$ for $\ell\geq 2$. The strength of $|\@vec{B}_{\rm turb}|$ relative to $|\@vec{B}_{\rm ord}|$ is parameterised by $f_{\rm 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 $\mathsc{Hi}$ emission maps (see section 4.3 in TG17), underlying the potential for $TE$ correlation and $E$-$B$ power asymmetry.

We can use the proxy $N_{\mathsc{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 $\@vec{B}_{\rm POS}$ and dust structures represented by the CNM template and model the polarization angles $\psi^{\rm c}$ as above. However, for the UNM and WNM templates we assume no such alignment777In this case, there is little direct evidence for or against alignment, and again this is a minimalist approach. In the actual turbulent ISM the evolutionary relationship between gas and fields in thermal phases is complex, whether time-dependent or steady-state, and probably varies with the position because of its particular history and the ambient pressure, radiative environment, etc. and so we apply a different procedure in which the polarization angles $\psi^{\rm u}$ and $\psi^{\rm w}$ follow from random Gaussian realisations of the components of $\@vec{B}_{\rm turb}$ (Planck Collaboration Int. XLIV 2016). We simulate the $\cos^{2}\thinspace\gamma$ factor in Eqs. (11) along each LOS for each of three templates using the relation (Planck Collaboration Int. XLIV 2016)

[TABLE]

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

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

• •

$\epsilon_{353}$ is the mean dust emissivity at 353 GHz (all model Stokes parameters scale with this value).

• •

$p_{0}$ is a polarization parameter that combines the polarization degree of interstellar dust grains and their alignment efficiency.

• •

Galactic coordinates $l_{0}$ and $b_{0}$ determine the direction of the ordered component of the GMF.

• •

$f_{\rm M}^{\rm c}$ parameterises the relative strength of $|\@vec{B}_{\rm turb}|$ and $|\@vec{B}_{\rm ord}|$ in the CNM template. The corresponding parameter is taken to be the same in the UNM and WNM templates and is denoted $f_{\rm M}^{\rm u/w}$.

• •

$\alpha_{\rm 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_{\rm m}$, $Q_{\rm m}$, and $U_{\rm 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 2018) 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_{\rm 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 $\@vec{B}_{\rm ord}$ is $l_{0}=62^{\circ}$ and $b_{0}=-18^{\circ}$ (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 $N_{\mathsc{Hi}}^{\rm c}$ template  we compute the polarization angles using the method described in Sect. 5. In particular, we compute maps of $\cos 2\psi^{\rm c}$. In the CNM the assumed alignment of the local GMF with respect to the gas structures constrains the pair of parameters, $f_{\rm M}^{\rm c}$ and $\alpha_{\rm M}$, that characterise $\@vec{B}_{\rm turb}$.

To evaluate these constraints, following TG17 we simulate $\cos 2\psi$ for different values of the pairs assuming the adopted orientation of $\@vec{B}_{\rm ord}$, using algebra described in section 4.1 of Planck Collaboration Int. XLIV (2016). For each pair, we compute means of power spectrum amplitudes of $\cos 2\psi$ from 100 Monte-Carlo realisations over mask65 and fit these with a power-law model ${\cal D}_{\ell}\propto\ell^{\alpha_{\rm M}+2}$ within multipole range $40\leq\ell\leq 160$ (Fig. 7). Comparing the power spectra of the $\cos 2\psi^{\rm c}$ map with the simulated spectra of $\cos 2\psi$ for different pairs of $f_{\rm M}$ and $\alpha_{\rm M}$, we adopt the following pair of parameter values for the CNM template: $f_{\rm M}^{\rm c}=0.5^{+0.05}_{-0.04}$ and $\alpha_{\rm M}=-2.4\pm 0.1$.

Applying the same procedure to the UNM and WNM templates results in a value of $f_{\rm M}^{\rm u/w}$ much higher than $f_{\rm M}^{\rm 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_{\rm M}^{\rm u/w}=0.1$.

6.2 Additional model parameters

We determine the three remaining model parameters, $\epsilon_{353}$, $p_{0}$, and $f_{\rm M}^{\rm u/w}$, by jointly minimising the following two expressions for $\chi^{2}$. One is in harmonic space (Fig. 6, Table 2):

[TABLE]

where $XX=\{EE,BB,TE\}$, ${\cal D}_{\ell}^{XX}$ and ${\cal M}_{\ell}^{XX}$ are binned power spectra of dust within the multipole range $40\leq\ell\leq 160$ for the data and model, respectively, and $\sigma_{\ell}^{XX}$ are the corresponding standard deviations. The other is in pixel space (Fig. 8):

[TABLE]

Through combining statistics of the polarized dust observations in both harmonic and pixel space, we break the degeneracy between the two model parameters $\epsilon_{353}$ and $p_{0}$ in Eq. (11). The $\chi^{2}$ of the T$-$T correlation between $I_{G353}$ and $I_{\rm 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 $\epsilon_{353}$ during the iterative solution of the two equations.

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

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_{\rm m}^{2}$, from 100 realisations using Eq. (11). 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_{\rm 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 compute the polarization angle, $\psi^{\rm R}_{\rm m}$, 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. (10). The polarization angle $\psi^{\rm R}_{i}$ is then derived from the rotated $Q_{i}$ and $U_{i}$ maps. The mean polarization angle, $\psi^{\rm R}_{\rm m}$, for the model is computed by taking the mean of the histograms of $\psi^{\rm R}_{i}$ 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^{\rm R}_{\rm d}$. The dispersion of $\psi^{\rm R}_{\rm m}$ derived from a Gaussian fit is $17$ .${}^{\circ}5\pm 0$ .${}^{\circ}4$, which is slightly higher than that for the data $\psi^{\rm R}_{\rm d}$: $16$ .${}^{\circ}6\pm 0$ .${}^{\circ}2$. The difference might come from the fact that the model is not fitted to the Planck data at $\ell<40$, which cannot be accomplished accurately over the limited area of mask65.

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: ${\cal M}_{\ell}^{XX}=A_{XX}(\ell/80)^{\alpha_{XX}+2}$ within the multipole range $40\leq\ell\leq 160$. Error bars in each bin are the standard deviations computed from 100 realisations. Best-fit amplitudes, $A_{XX}$ and exponents, $\alpha_{XX}$ along with the corresponding $\chi^{2}$ values, are listed in the right column of Table 2. We compare the results with the Planck dust power spectra in Fig. 6. The ${\cal D}_{\ell}^{EE}$ and ${\cal D}_{\ell}^{BB}$ amplitudes and their ratio match the values from the data accurately. The ${\cal D}_{\ell}^{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_{\rm 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\thinspace{\mu\rm K}_{\rm CMB}$, again within the uncertainty of the zero level in the intensity, $40\thinspace{\mu\rm K}_{\rm 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 $\delta$. An inverse relationship between $\mathcal{S}$ and $p_{\rm MAS}$ was found, where $p_{\rm 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 $\delta=30{{}^{\scriptstyle\prime}}$. Planck Collaboration XII (2018) extended this work over the full sky using $\tt GNILC$ processed intensity and polarization maps at a resolution of 160*′* (FWHM) with $\delta=80{{}^{\scriptstyle\prime}}$ and found that $\mathcal{S}\propto p_{\rm MAS}^{-1}$.

Here, we study the relationship between $\mathcal{S}$ and $p_{\rm MAS}$ over mask65 using 353 GHz maps at a resolution of 80*′* (FWHM) with $\delta=40{{}^{\scriptstyle\prime}}$. The noise bias parameter needed for $p_{\rm MAS}$ (Montier et al. 2015a, b) is estimated from the smoothed noise covariance matrices, ${\sigma}_{II},{\sigma}_{IQ},{\sigma}_{IU},{\sigma}_{QQ},{\sigma}_{QU},{\sigma}_{UU}$ at a resolution of 80*′* (FWHM) and $N_{\rm side}=128$.

In the top panel of Fig. 9, we present the two-dimensional joint distribution of $\mathcal{S}$ and $p_{\rm MAS}$ for the Planck data over mask65. The running mean of $\mathcal{S}$ in each bin of $p_{\rm MAS}$ follows an inverse relationship, $\mathcal{S}=(0$ .${}^{\circ}43\pm 0$ .${}^{\circ}04)/p_{\rm MAS}$. In the bottom panel of Fig. 9, we show the corresponding $\mathcal{S}$ - $p_{\rm MAS}$ distribution for our model. This too has an inverse relationship and the slope [math] .${}^{\circ}50\pm 0$ .${}^{\circ}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.

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^{\circ}$ (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 $\mathsc{Hi}$ data: in particular, three independent templates notionally representing the 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 $\psi$ 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_{\rm d}$ and $u_{\rm 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_{0},b_{0})=(62^{\circ},-18^{\circ})$ (Sect. 4.2). The best-fit value of $l_{0}$ 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 $\@vec{B}_{\rm POS}$ and the $\mathsc{Hi}$ structure in the CNM template (Sect. 5).

• •

The observed distributions of $p$ and $\psi$ over mask65 show a scatter that is comparable to that of the distributions reported for the whole sky in Planck Collaboration XII (2018). 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 $\mathsc{Hi}$ templates (Sect. 3.2).

• •

The best-fit value of the parameter, $p_{0}$, 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 $\epsilon_{353}=146\thinspace{\mu\rm K}_{\rm CMB}(10^{20}\thinspace$ 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_{\rm M}^{\rm 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_{\rm M}^{\rm c}=0.5$ characterising the CNM (Sect. 6.1).

• •

The spectral index of the turbulent component of the GMF, $\alpha_{\rm 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_{\rm 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 (2018) that the inverse relationship of $\mathcal{S}$ and $p_{\rm 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 $\mathsc{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. DA 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 $\tt 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 $\mathsc{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 $\mathsc{Hi}$ templates

Here we describe the development of an alternative set of $\mathsc{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 $\mathsc{Hi}$ emission spectra from an intermediate product of the HI4PI all-sky survey (HI4PI Collaboration et al. 2016), applying the algorithm described in Haud (2000) with some modifications (see their section 2.2). Based on line widths, KH18 defined $\mathsc{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_{\rm D}$), which depend on $|v|$ (see, e.g. their figures 3 and 7). For example, for velocities $|v|\thinspace\leq 8$  km s*-1* over the full sky, the mean line widths are 3.6, 9.6, and 23.3  km s*-1* 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_{\rm c}^{i}$ lies within the velocity range of $|v_{\rm c}|\thinspace\leq 24$  km s*-1*, which encompasses most of the local $\mathsc{Hi}$ gas (Kalberla & Haud 2018) but ignores most intermediate velocity gas (IVC). This decomposition was done at the native resolution in maps with $N_{\rm side}$= 1024. As in Kalberla & Haud (2018), maps of the integrated emission of $\mathsc{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_{\rm side}$= 128. The smoothing reduces uncertainties within each map that arise from separating the gas into components. These maps are shown in Fig. 10. Hereafter, we refer to these as the KH18 templates.

We refer to the sum of the column densities of the three templates as $N_{\mathsc{Hi},24}$. The correlation of dust emission with $N_{\mathsc{Hi},24}$ on the selected region mask65 is not satisfactory, having an emissivity $\epsilon_{353,24}$ of $121\thinspace{\mu\rm K}_{\rm CMB}(10^{20}\thinspace$ cm${}^{-2})^{-1}$, lower than $\epsilon_{353,50}$, and also an unacceptably large offset, $115\thinspace{\mu\rm K}_{\rm CMB}$. This arises because of the restricted velocity range adopted from KH18: $N_{\mathsc{Hi},24}\leq N_{\mathsc{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 $\mathsc{Hi}$. The fourth template is simply the column density map $N_{\mathsc{Hi},50}$ - $N_{\mathsc{Hi},24}$, also shown in Fig. 10.

The sum of the KH18 CNM and LNM templates is related to our CNM template used in the main paper, as shown in Fig. 11. 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 $\@vec{B}_{\rm POS}$ of the model GMF is aligned with $\mathsc{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_{\rm M}^{\rm c}=0.5$ and $\alpha_{\rm M}=-2.4$ (see Fig. 12). Because the error bars on $\cos 2\psi$ estimated from simulated maps are significant, on the basis of Fig. 12 we adopt the same pair of parameter values for the LNM. In the WNM, we adopt $f_{\rm M}^{\rm 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 $\epsilon_{353}=146\thinspace{\mu\rm K}_{\rm CMB}(10^{20}\thinspace$ cm${}^{-2})^{-1}$ and $p_{0}=19$ %.

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. 13. For $\ell$ bins between 40 and 160, the alternative set of $\mathsc{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 $\ell$-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.