A 3-D model of polarised dust emission in the Milky Way

TL;DR

This paper introduces a 3D model of polarized dust emission in the Milky Way, accounting for variations along the line of sight and small-scale fluctuations, to aid in CMB polarization studies.

Contribution

The paper presents a novel 3D model of polarized dust emission that incorporates line-of-sight variations and small-scale fluctuations constrained by observations.

Findings

Model matches large-scale dust emission observations.

Extrapolates observed intensity and polarization power spectra to smaller scales.

Provides a tool for assessing foreground impact on CMB polarization analysis.

Abstract

We present a three-dimensional model of polarised galactic dust emission that takes into account the variation of the dust density, spectral index and temperature along the line of sight, and contains randomly generated small scale polarisation fluctuations. The model is constrained to match observed dust emission on large scales, and match on smaller scales extrapolations of observed intensity and polarisation power spectra. This model can be used to investigate the impact of plausible complexity of the polarised dust foreground emission on the analysis and interpretation of future CMB polarisation observations.

Tables

Table 1. Table 1: Averages and standard deviation values of temperature and spectral index in each layer, for a simulation with 6.87 ′ pixels HEALPix pixels at Nside =512. The average and standard deviation of the resulting temperature and spectral index, as obtained from an MBB fit on the total intensity maps at 353, 545, 857 and 3000 GHz, is compared to what is obtained on Planck observations.

Layer	1	2	3	4	5	6
$T_{\text{avg}}$	19.10	18.96	18.98	19.35	19.23	20.05
${\sigma }_T$	2.059	2.100	2.022	2.076	2.117	2.069
${\beta }_{\text{avg}}$	1.627	1.628	1.598	1.538	1.513	1.689
${\sigma }_{\beta }$	0.209	0.210	0.207	0.208	0.202	0.204
	$T_{\text{avg}}^{\text{MMB}}$	${\sigma }_T^{\text{MMB}}$	${\beta }_{\text{avg}}^{\text{MMB}}$	${\sigma }_{\beta }^{\text{MMB}}$
Planck fit	19.396	1.247	1.598	0.126
Simul. fit	19.389	1.253	1.598	0.135

Equations

$$I_{\nu }=\tau ({\nu }_0){\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }{B}_{\nu }(T)$$

$$I_{\nu }={\int }_0^{\infty }dr\frac{d\tau (r,{\nu }_0)}{dr}{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta (r)}{B}_{\nu }(T(r)),$$

$$Q_{\nu }={\int }_0^{\infty }drp(r)\frac{d\tau }{dr}{B}_{\nu }(T(r)){\left(\frac{\nu }{{\nu }_0}\right)}^{\beta (r)}\cos 2\psi (r){\sin }^k\alpha (r)$$

$$U_{\nu }={\int }_0^{\infty }drp(r)\frac{d\tau }{dr}{B}_{\nu }(T(r)){\left(\frac{\nu }{{\nu }_0}\right)}^{\beta (r)}\sin 2\psi (r){\sin }^k\alpha (r)$$

$$\cos 2\psi =\frac{B_{\theta }^2-B_{\phi }^2}{B_{\perp }^2},\sin 2\psi =\frac{2B_{\theta }{B}_{\phi }}{B_{\perp }^2},\sin \alpha =\frac{B_{\perp }}{B},$$

$$I_{\nu }(p)=\sum _1^N{I}_{\nu }^i(p)=\sum _1^N{\tau }_i({\nu }_0){\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_i(p)}{B}_{\nu }(T_i(p)).$$

$$n_d(R,z)=n_0\exp (-R/h_R){\text{sech}}^2(z/h_z)$$

$${\tilde{I}}_{\nu }^i=f_i\tau ({\nu }_0){\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_i}{B}_{\nu }(T_i).$$

$$I_{{\nu }_0}^i={\tilde{I}}_{{\nu }_0}^i+f_i(I_{{\nu }_0}-{\tilde{I}}_{{\nu }_0}),$$

$$I_{\nu }^i=I_{{\nu }_0}^i{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_i}\frac{B_{\nu }(T_i)}{B_{{\nu }_0}(T_i)},$$

$$B_r=B(r,\theta ,z)\sin q,B_{\theta }=-B(r,\theta ,z)\cos q$$

$$B(r,\theta ,z)=-B_0(r)\cos \left(\theta +\beta \log \frac{r}{r_0}\right)\exp (-|z|/z_0),$$

$${\tilde{Q}}_{353}^i=\left(\frac{Q_m^i}{I_m^i}\right)I_{353}^i,{\tilde{U}}_{353}^i=\left(\frac{U_m^i}{I_m^i}\right)I_{353}^i,$$

$$G(q)=\sum _p\left(\left|{\tilde{Q}}_{353}^{\text{model}}(q)-Q_{353}^{\text{obs}}\right|+\left|{\tilde{U}}_{353}^{\text{model}}(q)-U_{353}^{\text{obs}}\right|\right)$$

$${\tilde{Q}}_{353}^{\text{model}}=\sum _{i=1}^N{\tilde{Q}}_{353}^i$$

$$Q_{353}^i=\left(\frac{Q_m^i}{I_m^i}\right)I_{353}^i+F_i\left[Q_{353}^{\text{obs}}-\sum _{j=1}^N\left(\frac{Q_m^j}{I_m^j}\right)I_{353}^j\right],$$

$$\begin{gathered} a_{\ell m}^{T,\mathrm{obs}}=b_{\ell }{a}_{\ell m}^{T,\mathrm{sky}}; \\ {a}_{\ell m}^{E,\mathrm{obs}}=h_{\ell }{a}_{\ell m}^{E,\mathrm{sky}} \end{gathered}$$

$$a_{\ell m}^{T,\mathrm{sim}}=a_{\ell m}^{T,\mathrm{obs}}+\sqrt{1-b_{\ell }^2}{a}_{\ell m}^{T,\mathrm{rnd}}$$

$$a_{\ell m}^{E,\mathrm{sim}}=a_{\ell m}^{E,\mathrm{obs}}+\sqrt{1-h_{\ell }^2}{a}_{\ell m}^{E,\mathrm{rnd}},$$

$$C_{\ell }^{TE,\mathrm{sim}}=C_{\ell }^{TE}\left[b_{\ell }{h}_{\ell }+\sqrt{1-b_{\ell }^2}\sqrt{1-h_{\ell }^2}\right].$$

$$a_{\ell m}^{\mathrm{final}}={\left[C_{\ell }\right]}^{1/2}[C_{\ell }^{\mathrm{sim}}{]}^{-1/2}{a}_{\ell m}^{\mathrm{sim}}$$

$$T_{\text{avg}}^i=\sum _p{w}_T^i(p)T_d(p)$$

$${\beta }_{\text{avg}}^i=\sum _p{w}_{\beta }^i(p){\beta }_d(p)$$

$$w_T^i(p)=w_{\beta }^i(p)=f_i(p),$$

Full text

††pubyear: 2016††pagerange: A 3-D model of polarised dust emission in the Milky Way–A 3-D model of polarised dust emission in the Milky Way

1 Introduction

Since the discovery of the cosmic microwave background (CMB) in 1965 (Penzias & Wilson, 1965), significant efforts have been devoted to precise characterisation of its emission, and to understanding the cosmological implications of its tiny temperature and polarisation anisotropies, detected first with COBE-DMR (Smoot et al., 1992) for temperature, and with DASI for polarisation (Kovac et al., 2002). Many experiments have gradually improved the measurement of CMB temperature and polarisation power spectra. Experiments on stratospheric balloons, notably Boomerang (de Bernardis et al., 2000; Jones et al., 2006), Maxima (Hanany et al., 2000), and Archeops (Benoit et al., 2002), detected with high significance the first acoustic peak in the CMB temperature power spectrum, and made the first measurements of the temperature power spectrum over a large range of angular scales. The WMAP satellite (Bennett et al., 2013) produced the first high signal-to-noise ratio full-sky CMB map and power spectrum from the largest scales to the third acoustic peak, opening the path to precision cosmology with the CMB. These observations have been completed by power spectra measurements from many ground based experiments, for instance ACBAR (Reichardt et al., 2009) and more recently ACT (Das et al., 2014) and SPT (Story et al., 2013) on scales smaller than observed with the balloons and space missions.

Planck, the latest space mission to-date, launched by ESA in 2009 (Tauber et al., 2010), has mapped CMB anisotropies with extraordinary precision down to $\simeq 5^{\prime}$ angular scale, providing a wealth of information on the cosmological scenario. The Planck Collaboration (2016) has shown that both the CMB temperature and E-mode polarisation power spectra were remarkably consistent with a spatially flat cosmology specified by six parameters, the so-called $\Lambda$CDM model, with cosmic structures seeded at very early times by quantum fluctuations of spacetime during an epoch of cosmic inflation.

The accurate measurement of cosmic microwave background polarisation, including inflationary and lensing $B$ modes, is the next objective of CMB observations. Such a measurement offers a unique opportunity to confirm the inflationary scenario, through the detection of the imprint of primordial inflationary gravitational waves on CMB polarisation $B$ modes on large angular scale (see Kamionkowski & Kovetz, 2016, for a review). CMB polarisation also offers the opportunity to map the dark matter in the universe that is responsible of slight distortions in polarisation patterns by the process of gravitational lensing of the background CMB (Lewis & Challinor, 2006; Challinor et al., 2017).

In 2014, the BICEP2 collaboration claimed evidence for primordial CMB $B$ modes with a tensor-to-scalar ratio $r=0.2$ (BICEP2 Collaboration, 2014). However, a joint analysis with Planck mission data (BICEP2/Keck and Planck Collaborations, 2015) showed that the signal was mostly due to contamination of the observed map by polarised dust emission from the Milky Way rather than gravitational waves from inflation. Future space missions such as COrE (The COrE Collaboration, 2011) and its more recent version, CORE (with a capital “R”), proposed to ESA in October 2016 in answer to the “M5” call for a medium-size mission (Delabrouille et al., 2017), PIXIE (Kogut et al., 2011), PRISM (André et al., 2014), LiteBIRD (Matsumura et al., 2014), and ground-based experiments such as CMB-S4 (Abazajian et al., 2016), plan to reach a sensitivity in $r$ as low as $r\sim 0.001$ (CORE Collaboration et al., 2016). This requires subtracting at least 99% of dust emission from the maps, or modeling the contribution of dust to the measured CMB B-mode angular power spectrum at the level of $10^{-4}$ precision or better. The feasibility of such dust-cleaning critically depends on the (unknown) complexity of dust emission down to that relative level, and on the number and central frequencies of frequency channels used in the observation (to be optimised in the design phase of future CMB experiments).

Investigations of the feasibility of measuring CMB $B$ modes in the presence of foreground astrophysical emission have been pursued by a number of authors (Tucci et al., 2005; Efstathiou et al., 2009; Betoule et al., 2009; Dunkley et al., 2009; Bonaldi & Ricciardi, 2011; Errard & Stompor, 2012; Bonaldi et al., 2014; Remazeilles et al., 2016; Stompor et al., 2016; Remazeilles et al., 2017), using component separation methods mostly developed in the context of the analysis of WMAP and Planck intensity and polarisation observations (see, e.g., Leach et al., 2008; Delabrouille & Cardoso, 2009, for reviews and comparisons of component separation methods). Conclusions on the achievable limit on $r$ drastically depend on the assumed complexity of the foreground emission model (see Delabrouille et al., 2013, for a widely used sky modelling tool), the number of components included, and on whether the component separation method that is used is or is not perfectly matched to the model used in the simulations.

In this paper we present a three-dimensional (3-D) model of the polarised dust emission, constrained by observations, that considers the spatial variation of the spectral index and of the temperature along the line of sight, and can help give insight on the feasibility and complexity of dust-cleaning in future CMB observations in the presence of a model of dust emission more complex and more realistic than what has been used in previous work. The objective is not an accurate 3-D model of dust emission, which cannot be obtained without additional observations of the 3-D dust, but a plausible 3-D model that is compatible with observed dust emission and its spatial variations, and at the same time implements a complexity which, although not strictly necessary yet to fit current observations, is likely to be detectable in future sensitive CMB polarisation surveys. This model can be used to infer properties such as decorrelation between frequencies and flattening of the spectral index at low frequencies, and also to test the possibility to separate CMB polarisation from that of dust with future multifrequency observations of polarised emission at millimeter wavelengths.

This paper is organized as follows. In section 2, we justify the need for 3-D modelling and discuss plausible consequences on the properties of dust maps across scales and frequencies. Section 3 presents the observations that are used in the construction of our dust model. In section 4, we present the strategy that is used to make a 3-D dust data cube in temperature and polarisation using the (incomplete) observations at hand. As these available observations have limited angular resolution, we describe in section 5 how to extend the model to smaller scales, in preparation for future high-resolution sensitive polarisation experiments. Section 6 describes our prescription for scaling the dust emission across frequencies. We compare simulated maps with existing observations and discuss implications of the 3-D model in section 7. We conclude in section 8.

2 Why a 3-D model?

Previous authors such as, e.g., Fauvet et al. (2011), O’Dea et al. (2012), and Vansyngel et al. (2016) have considered a 3-D model of dust distribution and of the Galactic magnetic field to model the spatial structure of dust polarisation. Ghosh et al. (2016) complement this with an analysis of correlations of the direction of the Galactic magnetic field with the orientation of dust filaments, as traced by Hi data. However, all of these approaches produce single templates of dust emission at a specific frequency but do not attempt at the same time to model the 3-D dependence of the dust emission law. This misses one of the key aspects of dust emission that is crucial to disentangling its emission from that of CMB polarisation (see Tassis & Pavlidou, 2015).

Dust is made of grains of different size and chemical composition absorbing and scattering light in the ultraviolet, optical and near-infrared, and re-radiating it in the mid- to far-infrared. Being made of structured baryonic matter (atoms, molecules, grains), dust interacts with the radiation field through many different processes. Empirically, at millimetre and sub-millimetre wavelengths, the observed emission in broad frequency bands is dominated by thermal emission at a temperature $T$, well fit in the optically thin limit by a modified black body (MBB) of the form

[TABLE]

where $I_{\nu}$ is the specific intensity at frequency $\nu$ and $B_{\nu}(T)$ is the Planck blackbody function for dust at temperature $T$. In the frequency range we are considering, the optical depth $\tau(\nu)$ scales as $(\nu/\nu_{0})^{\beta}$, where $\beta$ is a spectral index that depends on the chemical composition and structure of dust grains. Here $\nu_{0}$ is a reference frequency at which a reference optical depth $\tau(\nu_{0})$ is estimated (we use $\nu_{0}=353\,$GHz throughout this paper).

Using dust template observations in the Planck 353, 545, and 857 GHz channels and the IRAS 100$\,\mu$m map, it is possible to fit for $\tau(\nu_{0})$, $T$ and $\beta$ in each pixel. This fit, performed by the Planck Collaboration (2014), shows clear evidence for a variation across the sky of the best-fit temperature and spectral index, with $T$ mostly ranging from about 15 K to about 27 K and $\beta$ ranging from about 1.2 to about 2.2. Such variations are expected by reason of variations of dust chemical composition and size, and of variations of the stellar radiation field, as a function of local physical and environmental conditions.

In this paper, we propose to revisit this model to make it three-dimensional. Indeed, if dust properties vary across the sky, they must also vary along the line of sight (LOS). This means that even if one single modified blackbody is (empirically) a good fit to the average emission coming from a given region of the 3-dimensional Milky Way as observed with the best current signal-to-noise ratio, the integrated emission in a given LOS must be a superposition of several such MBB emissions with varying $T(r)$ and $\beta(r)$ (in fact, a continuum, weighted by a local elementary optical depth $d\tau(r,\nu_{0})$):

[TABLE]

where $r$ is the distance along the line of sight and where, again, $\tau(r,\nu_{0})$ is an optical depth at frequency $\nu_{0}$, $T(r)$ is a temperature, and $\beta(r)$ a spectral index, now all dependent on the distance $r$ from the observer.

As a sum of MBBs is not a MBB, this mixture of dust emissions is at best only approximately a MBB. For instance, regions along the line of sight with lower $\beta$ contribute relatively more at low frequency than at high frequency. This would then naturally generate an effect of flattening of the observed dust spectral index at low $\nu$, which precludes fits of dust emission performed at high frequency to be valid at lower frequencies. To properly account for such LOS inhomogeneities, a 3-D model of dust emission, with dust emission law variations both across and along the line of sight, is needed.

This 3-D mixture of inhomogeneous emission would also naturally impact the polarised emission of galactic dust. The preferential alignment of elongated dust grains perpendicularly to the local magnetic field $\vec{B}$ results in a net sky polarisation that is, on the plane of the sky, orthogonal to the component $\vec{B}_{\perp}$ of $\vec{B}$ that is perpendicular to the line of sight. The efficiency of grain alignment depends on the local physical properties of the interstellar medium (density, which impacts the collisions between grains; irradiation). Each region emits polarised emission proportional to an intrinsic local polarisation fraction $p(r)$. Linear polarisation Stokes parameters $Q$ and $U$ can be written as

[TABLE]

and

[TABLE]

where, in the HEALPix CMB polarisation convention,

[TABLE]

and where $k$ is an exponent that takes into account de-polarisation and projection effects linked to the local geometry and the alignment of grains. In these equations, $r$ is the distance to the observer, i.e. $r$, $\theta$ and $\varphi$ are spherical heliocentric coordinates. In Eqs. 3 and 4, we recognise an overall intensity term (equal to the integrand in Eq. 2), multiplied by a polarisation fraction $p(r)$, an orientation term $\cos 2\psi(r)$ or $\sin 2\psi(r)$, and a geometrical term $\sin^{k}\alpha(r)$ that depends on the direction of the magnetic field with respect to the line of sight. In the absence of strong theoretical or observational constraints on the value of $k$, we follow Fauvet et al. (2011) and assume $k=3$. This choice, although arguably somewhat arbitrary, does not impact much the rest of this work111Nor does the specific analytic form of the de-polarisation function. as it does not change the polarisation angle on the sky, while the polarisation maps will ultimately be re-normalised to match the total observed dust polarisation at 353 GHz. This re-normalisation somewhat corrects for possible inadequacy or inaccuracy of the assumption made for the geometrical term.

Since all parameters ($p$, $\tau$, $T$, $\beta$, $\psi$, and $\alpha$) vary along the line of sight, the total polarised emission is a superposition of emissions with different polarisation angles and different emission laws. As a consequence, the polarisation fraction will change with frequency (i.e. intensity and polarisation have different emission laws); in addition, the polarisation will rotate as a function of frequency, depending on the relative level of emission of various regions along the line of sight. This polarisation rotation effect would also naturally generate decorrelation of polarisation emission at various frequencies. Such an effect that has been reported in Planck observations (Planck Collaboration et al., 2016a), but is the object of debate following a subsequent analysis that does not confirm the statistical significance of the observed decorrelation (Sheehy & Slosar, 2017).

3 Observations

Full-sky (or near-full-sky) dust emission is observed at sub-millimetre wavelength by Planck and IRAS. We process the Planck 2015 data release maps with a Generalized Needlet Internal Linear Combination (GNILC) method to separate dust emission from other astrophysical emissions and to reduce noise contamination. GNILC (Remazeilles et al., 2011) is a component separation method that extracts from noisy multi-frequency observations a multiscale model of significant emissions, based on the comparison of auto and cross-spectra with the level of noise locally in needlet space. Needlets (Narcowich et al., 2006; Marinucci et al., 2008; Faÿ et al., 2008) are a tight frame of space-frequency functions (which serve as a redundant decomposition basis). The use of needlets for component separation by Needlet Internal Linear Combination (NILC) was introduced in the analysis of WMAP 5-year temperature data (Delabrouille et al., 2009). They were further used on the 7-year and 9-year temperature and polarisation maps (Basak & Delabrouille, 2012, 2013).

GNILC has been used by the Planck collaboration to separate dust emission from Cosmic Infrared Background (Planck Collaboration et al., 2016b). We use the corresponding dust maps to constrain our model of dust emission in intensity. GNILC maps offer the advantage of reduced noise level (for both intensity and polarisation), and of reduced contamination by the cosmic infrared background fluctuations (for intensity). However, different templates of dust emission in intensity and polarisation could have been used instead, as long as those maps are not too noisy, nor contaminated by systematic effects such that, for instance, the intensity map is negative in some pixels, or that dust is a sub-dominant component in some pixels or at some angular scales (problems of that sort are usually present in maps that have not been processed to avoid these issues specifically).

From now on, the single greybody ‘2-D’ model of the form of Eq. 1 uses Planck maps of $\tau(\nu_{0})$ and $\beta$ that are obtained from a fit of the GNILC dust maps between 353 and 3000 GHz, obtained as described in Planck Collaboration et al. (2016b). For polarisation, we apply independently GNILC on Planck 30-353 GHz $E$ and $B$ polarisation maps (Data Release 2, Planck Collaboration et al., 2016d) to obtain polarised galactic emission maps in the seven polarised Planck channels. These maps are specifically produced for the present analysis, and are not part of the Planck archive. Dust-dominated $E$ and $B$ polarisation maps at $\nu=143$, 217 and 353 GHz are shown in Fig. 3. The polarisation maps with best dust signal-to-noise ratio are at $\nu=353\,$GHz. The other polarisation maps are not used further in our model.222After the GNILC process to de-noise the observations, these maps bring only limited additional information: considering their noise level, their dust component over most of the sky is obtained largely by GNILC from their correlation with the 353 GHz map locally in needlet space. Here, the GNILC processing is mostly used as a pixel-dependent de-noising of the 353 GHz polarisation map. A model that fully exploits the multi-frequency information in the Planck data is postponed to future work. The needlet decomposition extends down to $5^{\prime}$ angular resolution for intensity and $1^{\circ}$ for polarisation.

Three-dimensional maps of interstellar dust optical depth, as traced by starlight extinction, have been derived by Green et al. (2015) based on the reddening of 800 million stars detected by PanSTARRS 1 and 2MASS, covering three-quarters of the sky. The maps are grouped in 31 bins out to distance moduli from 4 to 19 (corresponding to distances from 63 pc to 63 kpc) and have a hybrid angular resolution, with most of maps at an angular resolution of $3.4^{\prime}$ to $13.7^{\prime}$. These maps will be used to infer some information about the distribution of dust along the line of sight, which will be used to generate our three-dimensional model of polarised dust emission.

4 Multilayer modeling strategy

We approximate the continuous integrals of Eqs. 2, 3 and 4 as discrete sums over independent layers of emission, indexed by $i$, so that we have, for the intensity,

[TABLE]

Each layer is then characterized by maps of Stokes parameters $I_{\nu}^{i}(p)$, $Q_{\nu}^{i}(p)$ and $U_{\nu}^{i}(p)$, with a frequency scaling, for each sky pixel $p$, in the form a single MBB emission law with a temperature $T_{i}(p)$ and a spectral index $\beta_{i}(p)$ (both assumed to be the same for all three Stokes parameters). We want to find a way to assign to each such layer plausible templates (full-sky pixellised maps) for $I$, $Q$ and $U$ at some reference frequency $\nu_{0}$, as well as scaling parameter maps $T$ and $\beta$, all such that the total emission matches the observed sky. By ‘layer’ we mean a component, loosely associated to different distances from us, but which could equally well be a component associated to a specific population of dust grains.

The problem is clearly degenerate. Starting from only four dust-dominated maps of $I$ (Planck and IRAS maps from 353 to 3000 GHz obtained after the GNILC analysis to remove CIB contamination), and one map of each of $Q$ and $U$ (both at 353 GHz), for a total of 6 maps, we propose to model dust emission with $3N$ maps of Stokes parameters $I_{\nu_{0}}^{i}(p)$, $Q_{\nu_{0}}^{i}(p)$ and $U_{\nu_{0}}^{i}(p)$ and $2N$ maps of emission law parameters $T_{i}(p)$ and $\beta_{i}(p)$, i.e. a total of $5N$ maps, where $N$ is the number of layers used in the model.

For any $N\geq 2$, we need additional data or constraints. We thus use the 3-D maps of dust extinction from Green et al. (2015) to decompose the observed intensity map $I$ at some reference frequency as a sum of intensity maps $I_{i}$ coming from different layers $i$. We group the dust extinction maps in 6 ‘layers’ (shown in Fig. 1) by simple coaddition of the corresponding optical depths. Six layers are sufficient for our purpose and provide a better estimate of the optical thickness associated to each layer than if we tried to use more. Three of these layers map the dust emission at high galactic latitude, while three map most of the emission close to the galactic plane. We choose the smallest possible homogeneous pixel size, corresponding to HEALPix Nside$=64$. These choices could be revisited in the future, in particular when more data become available.

We then further use a 3-D model of the Galactic magnetic field to generate $Q$ and $U$ maps for each layer. Finally, the total emission from all layers is readjusted so that the sum matches the observed sky at the reference frequency. We detail each of these steps in the following subsections.

4.1 Intensity Layers

Although the general shape and density distribution of the Galaxy is known, the exact 3-D density distribution of dust grains in the Galaxy is not. Simple models consider a galacto-centric radius and height function:

[TABLE]

where $(R,z)$ are cylindrical coordinates centred at the Galactic center, and where $h_{R}=3\,$kpc and $h_{z}=0.1\,$kpc. Such models cannot reproduce the observed intermediate and small scale structure of dust emission.333However, they can be used to get an initial estimate of the dust density on very large scale. We will make use of this in the next section for an initial guess of the polarisation fraction of dust emission in each layer.

On the other hand, the maps of Green et al. (2015) trace the dust density distribution, and are directly proportional to the optical depth $\tau$ at visible wavelength. We select six primary shells within distance moduli of 4, 7, 10, 13, 16, and 19 (corresponding to distances of 63, 251, 1000, 3881, 15849, and 63096 pc from the Sun), and use those maps to compute, in each pixel, an estimate of the fraction $f_{i}(p)$ of the total opacity associated to each layer (so that $\forall p$, $\sum_{i}f_{i}(p)=1$). We then construct the opacity map for each layer as the product $\tau_{i}(\nu_{0})=f_{i}\,\tau(\nu_{0})$, where $\tau(\nu_{0})$ is the opacity at 353 GHz obtained in the Planck MBB fit.

For our 3-D model, we must face the practical difficulty that the maps of Green et al. (2015) do not cover the full sky (Fig. 1). For a full-sky model, the missing sky regions must be filled-in with a simulation or a best-guess estimate. We use the maps where they are defined to evaluate the relative fraction $f_{i}$ of dust in each shell $i$. For each pixel where the layers are not defined, we use symmetry arguments and copy the average fraction from regions centred on pixels at the same absolute Galactic latitude and longitude. This gives us a plausible dust fraction in the region not covered in the decomposition of Green et al. (2015). We then use these fractions of emission to decompose the total map of optical depth $\tau(\nu_{0})$ at 353 GHz and obtain the six maps of extinction shown in Fig. 2.

We then compute the corresponding brightness in a given layer by multiplying by the Planck function together with the spectral index correction (Eq. 8), using for this an average temperature and spectral index for each layer.444We postpone to section 6 the discussion of temperature and spectral index maps. In equation 8, we use average values given in table 1. We get, for each layer $i$, an initial estimate of the intensity

[TABLE]

The sum $\widetilde{I}_{\nu_{0}}=\sum_{i}\widetilde{I}^{i}_{\nu_{0}}$ however does not exactly match the observed Planck map $I_{\nu_{0}}$ at $\nu_{0}=353$ GHz. We readjust the layers by redistributing the residual error in the various layers, with weights proportional to the fraction of dust in each layer, to get:

[TABLE]

and by construction we now have ${I}_{\nu_{0}}=\sum_{i}{I}^{i}_{\nu_{0}}$. The full model across frequencies is:

[TABLE]

with $I^{i}_{\nu_{0}}$ computed following equations 8 and 9. In this way, we have six different maps of dust intensity that add-up to the observed Planck dust intensity emission at 353 GHz.

We note that our model differs from that of Vansyngel et al. (2016), who instead make the simplifying assumption that the intensity template in all the layers they use is the same. The consequence of this approximation is that the fraction $f_{i}$ of emission in all the layers is constant over the sky. This is not compatible with a truly three-dimensional model: galactic structures cannot be expected to be spread over all layers of emission with a proportion that does not depend on the direction of observation.

The decomposition we implement in our model is just one of many possible ways to separate the total map of dust optical depth into several contributions. A close look at what we obtain shows several potential inaccuracies. For instance, some compact structures are clearly visible in more than one map, while it is not very likely that they all happen to be precisely at the edge between layers or elongated along the line of sight so that they extend over more than one layer. This ‘finger of God’ effect is likely to be due to errors in the determination of the distance or of the extinction of stars, which, as a result, spreads the estimated source of extinction over a large distance span. The north polar spur (extending from the galactic plane at $l\simeq 30^{\circ}$, left of the Galactic center, towards the north Galactic pole) is clearly visible both in the first two maps. According to Lallement et al. (2016), it should indeed extend over both layers. On the other hand, structures associated to the Orion-Eridanus bubble (right of the maps, below the Galactic plane) can be seen in all three first maps, from less than 60 pc to more than 250 pc, while most of the emission associated to Orion is at a distance of 150-400 pc. As discussed by Rezaei Kh. et al. (2016), future analyses of the Gaia satellite data are likely to drastically improve the 3-D reconstruction of Galactic dust. For the present work, we use the maps of Fig. 2, noting that for our purpose what really matters is not the actual distance of any structure, but whether such a structure is likely to emit with more than one single MBB emission law. Certainly, a complex region such as Orion cannot be expected to be in thermal equilibrium and constituted of homogeneous populations of dust grains, and thus modelling its emission with more than one map is in fact preferable for our purpose. The same holds for distant objects such as the large and small Magellanic clouds and associated tidal structures, wrongly associated to nearby layers of emission by the procedure we use to fill the missing sky regions. Hence, the ‘layers’ presented here should be understood as layers of emission with roughly one single MBB (per pixel), originating mostly from a given range of distances from the Earth (see also Planck Collaboration et al., 2016e, for a discussion of emission layers and their connection to spatial shells or different phases of the ISM). While this decomposition is not exact, it matches the purposes of the present work.

4.2 Polarisation layers

We model polarisation using Eq. 3 and 4. Geometric terms depending on $\psi$ and $\alpha$ are computed using a simple large scale model of the Galactic magnetic field (GMF). This regular magnetic field is assumed to roughly follow the spiral arms of the Milky Way. Several plausible configurations have been proposed, based on rotational symmetry around the Galactic centre, and on mirror symmetry with respect to the Galactic plane. A widely used parametrization, named in the literature as bisymmetric spiral (BSS, Sofue & Fujimoto, 1983; Han & Qiao, 1993; Stanev, 1997; Harari et al., 1999; Tinyakov & Tkachev, 2002), defines the radial and azimuthal field components (in Galactocentric cylindrical coordinates) as:

[TABLE]

where $q$ is the pitch angle of the logarithmic spiral, and where the function $B(r,\theta,z)$ is defined as:

[TABLE]

where $\beta=1/\tan q$. We model the regular magnetic field using such a BSS parameterization, in which we consider the $z$-component of the Galactic magnetic field to be zero. The model is restricted for $r>1$ kpc to avoid divergence of the field at small radius (and is hence assumed to vanish for $r\leq 1$ kpc).

The value of the pitch angle of the spiral arms in the Milky Way is still a matter of debate in the community. Estimates of this angle range from $-5^{\circ}$ to $-55^{\circ}$ depending on the tracer used to determine it, with the most commonly cited value being around $-11.5^{\circ}$. A possible explanation for the wide range of pitch angles determined from different datasets is that the pitch angle is not constant but varies with radius, meaning the spirals are not exactly logarithmic (e.g. slightly irregular).

In our case, the model should reproduce as well as possible the polarised dust emission on large scales, and at high galactic latitude in particular. The simple large scale density model of Eq. 7 together with the BSS large scale magnetic field from Eqs. 11 and 12 can be integrated following Eqs. 2, 3 and 4 to provide a first guess of dust intensity and polarisation distribution for each layer ($I^{i}_{m},Q^{i}_{m},U^{i}_{m}$). We initially assume that the intrinsic local polarisation fraction $p(r)$ in Eqs. 3 and 4 is constant and equal to $20\%$. Since we already have layers of intensity emission ($I^{i}_{353}$), the polarised emission in each layer $i$ can be generated as:

[TABLE]

The best fit pitch angle $q$ can be found minimizing some function of the difference between the simple polarisation model obtained from Eq. 13 and the observations. We minimize the L1 norm of the difference in $Q$ and $U$, summed the for all the pixels at high galactic latitude:

[TABLE]

where the dependence on the pitch angle $q$ has been specified for clarity, and where the total modelled $Q$ is the sum of the simple layer contributions from Eq. 13:

[TABLE]

and similarly for $U$. We find that a pitch angle of $-33^{\circ}$ provides the best fit of the GNILC maps by the BSS model at galactic latitude $|b|\geq 15^{\circ}$, which is the region of the sky with more interest for CMB observations.

Finally, to match the observations, we redistribute in the modelled layers of emission the residuals (observed emission minus modelled emission for $q=33^{\circ}$) weighted with some pixel-dependent weights $F_{i}$:

[TABLE]

This guarantees that the model matches the observation at 353 GHz on the angular scales that are observed with good signal-to-noise ratio by Planck. However, these weights $F_{i}$ must be such that the polarisation fraction after the redistribution of residuals does not exceed some maximum value $p_{\rm max}$, which is a free parameter of our model, and which we pick to be 25%. We fix the value of $F_{i}$ as $F_{i}=P_{i}/\sum_{j}P_{j}$, i.e., proportionally to the polarised dust emission fraction in each layer, unless the resulting polarisation fraction exceeds $p_{\rm max}$. When this happens, we redistribute the polarisation excess in neighbouring layers. The first term in the sum on the right hand side of Eq. 16 is the predicted polarisation of layer $i$, based on a polarisation fraction predicted by the BSS magnetic field applied to an intensity map for that layer. The second term is the correction that is applied to force the sum of all layers’s emissions to match the observed sky. The $U$ Stokes parameters is modelled in a similar way.

With this approach, we straightforwardly constrain the sum of emissions from all the layers to match the total observed emission for both $Q$ and $U$. Fig. 4 shows the polarised layers $Q_{m}^{i}$ and $U_{m}^{i}$ given by the large scale model of the magnetic field while Fig. 5 shows the polarised layers after redistributing the residuals all over the former layers. After adding the small scales features (next section), we get the maps displayed in Fig. 6. A visual comparison with Fig. 4 shows that while the regular BSS field model does a reasonable job at predicting the very large scale polarisation patterns (lowest modes of emission) at high galactic latitude (after picking the appropriate pitch angle), it fails at predicting most of the features of the observed polarised dust emission on intermediate scales. In addition, the amplitude of the modelled polarised emission at high Galactic latitude is seen to be too strong as compared to the observations. It is thus important, for the modelled emission to be reasonably consistent with Planck data, to enforce that the model match the observations, as we do, and not just rely on a simple regular model of the magnetic field, which does not exactly capture the observed features of the real emission.

5 Small scales

The polarisation maps we have generated are normalised to match the observed dust polarisation in the GNILC 353 GHz maps obtained as described in section 3. However, the polarisation GNILC maps are produced at $1^{\circ}$ angular resolution. In the galactic plane, where the polarised signal is strong, this is the actual resolution of the GNILC map. At high galactic latitude however, the amount of polarised sky emission power is low compared to noise even at intermediate scales. The GNILC processing then ‘filters’ the maps to subtract noise when no significant signal is locally detected. Hence, there is a general lack of small scales in the template $E$ and $B$ maps used to model polarised emission so far: everywhere on the sky on scales smaller than $1^{\circ}$ (because the GNILC maps are produced at $1^{\circ}$ angular resolution), but also on scales larger than that at high galactic latitude (because of the GNILC filtering). We must then complement the maps with small scales in a range of harmonic modes that depends on the layer considered, the first three layers covering most of the high galactic latitude sky, and the last three dominating the emission in the galactic plane and close to it, where the NILC filters less of the intermediate scales.

Small angular scale polarised emission arises from both small scale distribution of matter in three dimensions, but also from the fact that on small scales, the magnetic field becomes gradually more irregular, tangled and turbulent. Fully characterizing the strength, direction, and structure of the GMF in the entire Milky Way is a daunting task, involving measurements of very different observational tracers (see Han (2017) for a recent review). This field can be considered as a combination of a regular field as discussed above, complemented by a turbulent field that is caused by local phenomena such as supernova explosions and shock waves. The GMF is altered by gas dynamics, magnetic reconnection, turbulence effects. Observations constrain only one component of the magnetic field (e.g. strength or direction, parallel or perpendicular to the LOS) in one particular tracer (ionized gas, dense cold gas, dense dust, diffuse dust, cosmic-ray electrons…). This provides us with only partial information, making it extremely difficult to generate an accurate three-dimensional picture. The small scale magnetic field can be modelled with a combination of components that can be isotropic, or somewhat ordered with, e.g., a direction that does not vary on small scales while the sign of the $B$ vector does, as illustrated in Fig. 1 of Jaffe et al. (2010). The amplitude of these small scale fields depend on the turbulent energy density. In both the Milky Way and in other spiral galaxies, the fields have been found to be more turbulent within the material spiral arms than in between them (Jaffe et al., 2010). Different strategies to constrain the strength of the random magnetic fields (including or not both turbulent fields) estimate an amplitude of the turbulent field of about the same order of magnitude as that of the regular part, ranging however from 0.7 to 4 $\mu G$ for different estimates (Haverkorn, 2015). In a typical model, the power spectrum of the random magnetic field is assumed to follow a Kolmogorov spectrum (with spectral index $n=5/3$) with an outer scale of 100 pc.

In our work, we do model the large scale, regular magnetic field using the BSS model of Eqns. 11 and 12 to get a first guess of the layer-dependent dust polarisation, but we do not attempt to directly model the 3-D turbulent magnetic field. Indeed, it is not possible to implement a description of the real field down to those small scales, by lack of observations. The alternate strategy that consists in generating a random turbulent magnetic field, as in Fauvet et al. (2011), generates fluctuations with random phases and orientations, and dust polarisation fluctuations that cannot be expected to match those observed in the real sky. Hence – as we detail next – we propose instead to rely on the observed polarised dust on scales where those observations are reliable, and extend the power spectra of our maps at high $l$ in polarisation, independently for each layer, to empirically model the effect of a small scale turbulent component of the galactic magnetic field, on scales missing or noisy in the GNILC 353 GHz map.

To do so, we add small scales fluctuations independently in each layer of our model, both for intensity and for polarisation. In the case of intensity, we simply fit the power spectrum of the original map in the multipole interval $30\leq l\leq 300$, obtaining spectral indexes in harmonic space ranging from $-2.2$ to $-3.2$ as a function of the layer (steeper at further distances). We use these fitted spectral indexes to generate maps of fake intensity fluctuations, generated with a log-normal distribution in pixel space (so that the dust emission is never negative), with an amplitude proportional to the large scale intensity, and globally adjusted to match the level of the angular power spectrum. We use a similar prescription for $E$ and $B$, except that following the Planck results presented in Planck Collaboration et al. (2016c), we assume a power-law dependence for $EE$ and $BB$ power spectra at high $l$, of the form $C_{l}=A(l/l_{\text{fit}})^{\alpha}$ with $\alpha=-2.42$ for both $E$ and $B$. We use a Gaussian distribution, instead of log-normal, for polarisation fields. For each layer, we fix the amplitude $A$ and $l_{\text{fit}}$ to match the power spectrum of the large-scale map for that layer in the range $30\leq l\leq 100$. The amplitude of the small scale fluctuations is scaled by the polarised intensity map in each layer. The randomly-generated $T$ and $E$ harmonic coefficients are drawn with 30% correlation between the two, while $B$ is uncorrelated with both $T$ and $E$.

We then make combined maps which use large scales from the observations, and the smallest scales from the simulations, as follows. For each layer, we have an observed map, with a beam window function $b_{\ell}$ for temperature and $h_{\ell}$ for polarisation, i.e.

[TABLE]

and we have available $a_{\ell m}^{T,\rm{rnd}}$ and $a_{\ell m}^{E,\rm{rnd}}$ randomly generated following modelled statistics $C_{\ell}^{TT}$, $C_{\ell}^{EE}$ and $C_{\ell}^{TE}$, which we assume match the statistics of real sky emission. We complement the observed $a_{\ell m}$ by forming

[TABLE]

and similarly

[TABLE]

i.e., we make the transition between large and small scales in the harmonic domain using smooth harmonic windows, corresponding to that of a Gaussian beam of $5^{\prime}$ for all layers in intensity, $2.5^{\circ}$ for polarisation layers 1,2 and 3 (emission mostly at high galactic latitude), and of $2^{\circ}$ for polarisation layers 4, 5 and 6 (emission mostly near to the Galactic plane). These simulated sets of $a_{\ell m}$ have correct $C_{\ell}^{TT}$ and $C_{\ell}^{EE}$, but not cross spectrum $C_{\ell}^{TE}$. Indeed

[TABLE]

we obtain final simulated $a_{\ell m}$ as

[TABLE]

where for each $\ell$, $C_{\ell}$ and $C_{\ell}^{\rm sim}$ are $2\times 2$ matrices corresponding to the terms of the multivariate ($T$,$E$) power spectra of the model and of the simulated maps with small scales added.

Fig. 7 shows the maps of polarised emission after the various steps of our simulation process, summing up the contributions of all layers. Final maps of polarized intensity can be seen in Fig. 8. The percentage of polarized pixels with a given polarisation fraction decreases with the polarisation fraction, as seen in Fig. 9.

The power spectra of simulated maps in all the layers after this full process are shown in Fig. 10. The power spectra of the original GNILC maps with those resulting from the individual sum of the simulations with small scales fluctuations added in each layer is shown in Fig. 11: the missing power on small scales is complemented with fake, simulated small scale fluctuations. We show full-sky maps of $E$ and $B$ at 353 GHz in Fig. 12. A detail at $(l,b)=(0^{\circ},50^{\circ})$ is shown in Fig. 13. $E$ and $B$ power spectra of the original GNILC maps and the simulations at 143 and 217 GHz are shown in Fig. 14.

6 Scaling laws

We now need a prescription for scaling the 353 GHz polarised dust emission templates obtained above across the range of frequencies covered by Planck and future CMB experiments. We stick with the empirical form of dust emission laws (for each layer, a modified black body with pixel-dependent temperature and spectral index), but now, we must define as many templates of $T(p)$ and $\beta(p)$ as there are layers in our model, i.e. 6 maps of $T$ and 6 maps of $\beta$.

A complete description of the temperature and the spectral index distribution in 3-D would require observations of the intensity emission at different frequencies in each layer, which are not presently available. We cannot either use for each layer the same temperature and spectral index maps (otherwise there is no point using several layers to model the total emission). Finally, the scaling law we use for all the layers must be such that the final dust emission across frequencies should match the observations, i.e.: i- On average, the dust intensity scaled to other Planck frequencies (besides 353 GHz, at which matching the observations is enforced by construction) should be as close as possible to the actual Planck observed dust intensity; ii- Similarly, each of the dust Q and U polarisation maps, scaled to other frequencies than 353 GHz, should match the observed polarisation at those frequencies; iii- If we perform a MBB fit on our modelled dust intensity maps, the statistical distribution of temperature and spectral index should match those observed on the real sky: same angular power spectra and cross-spectra, similar non-stationary distribution of amplitudes of fluctuations across the sky, similar $T$-$\beta$ scatter plot.

With only 353 GHz polarisation maps with good signal-to-noise ratio, we construct our model for frequency scaling on intensity alone. In a first step, we make use of the fraction of dust assigned to each layer to compute the weighted mean of the spectral index and temperature maps for each layer, using the overall maps obtained from the MMB fit made in Planck Collaboration (2014), which we assume to hold for all of the $I$, $Q$ and $U$ Stokes parameters. We compute, for each layer $i$:

[TABLE]

where $T_{d}(p)$ and $\beta_{d}(p)$ are the best-fit values of the overall MMB fit of Planck dust emission in each pixel, and where $w^{i}_{T}(p)$ and $w^{i}_{\beta}(p)$ are some weights used for computing the average. We use the same weights both for temperature $T_{\text{avg}}^{i}$ and spectral index $\beta_{\text{avg}}^{i}$,

[TABLE]

i.e. we empirically weight the maps by the pixel-dependent fraction $f_{i}(p)$ of dust emission in layer $i$, to take into account the fact that we are mostly interested in the temperature and spectral index of the regions of sky where that layer contributes most to the total emission.

The simplest way to scale to other frequencies is to assume that $T_{\text{avg}}^{i}$ and $\beta_{\text{avg}}^{i}$ are constant across the sky in a given layer. This however implements only a variability of the physical parameters $T$ and $\beta$ along the LOS, and not across the sky anymore. It provides a (uniform) prediction of the scaling law in each layer that is informed by the observed emission law, but which does not reproduce the observed variability across the pixels of the globally fitted $T$ and $\beta$ (even if a global fit might find fluctuations because of the varying proportions of the various layers in the total emission as a function of sky pixel).

To generate fluctuations of the spectral index and temperature of dust emission in each layer, we first generate, for each layer, Gaussian random variations around $T_{\text{avg}}^{i}$ and $\beta_{\text{avg}}^{i}$ following the auto and cross spectra of the MBB fit obtained on the observed Planck dust maps (Planck Collaboration et al., 2016b).

To take into account the non-Gaussianity of the distribution of $T$ and $\beta$, we then re-map the fluctuations to match the observed probability distribution function in pixel space. This slightly changes the map spectra. As a final step, we thus re-filter the maps to match the observed auto and cross-spectra of $T$ and $\beta$. One such iteration yields simulated temperature and spectral index maps in good statistical agreement with the observations. In Fig. 15 we show a random realization of temperature and spectral index maps for the first layer, its power spectra and its scatter plot.

We then model the total emission at 353, 545, 857 GHz and 100 microns using those scaling laws, summing-up contributions from all six layers, and, in order to validate that the simulation is compatible with the observations, check with a MBB fit on the global map whether the distribution of the fitted parameters for the model is similar to that inferred on the real Planck observations.

We find two problems. First, the average temperature and spectral index fit on the total emission turn out to be slightly larger and smaller respectively than observed on the real sky. This is not surprising: as the emission in each pixel is a sum, the layer with the largest temperature and the smallest spectral index tends to dominate the total emission both at $\nu=3$ THz, pulling the temperature towards higher values, and at low frequency, pulling the spectral index towards lower values. We find that the average MBB fit temperature from the model matches the observations if we rescale the temperature in individual layers by a factor 0.982. Secondly, the standard deviations of the resulting fitted $T$ and $\beta$ are significantly smaller than those of the real sky, presumably because of averaging effects. We recover a global distribution of temperature and spectral index as fit on the total emission, if we rescale the amplitude of the temperature and spectral index fluctuations generated in each layer. We find that to match the observed $T$ and $\beta$ inhomogeneities of the MBB fit performed on GNILC Planck and IRAS maps, we need to multiply the amplitude of temperature fluctuations in each layer by 1.84 and the spectral index fluctuations by 1.94. With this re-scaling, we find a good match between the simulated and the observed temperature and spectral index distributions in the global MBB fit. Table 1 shows the standard deviation and the average values for $T$ and $\beta$ in each layer for one single realisation of the simulation, compared with those from the Planck MBB fit and the fit performed on this realisation. The average values from several simulations are in good agreement with those of the Planck MMB fit.

7 Validation and predictions

We use our model to generate maps of polarised dust emission at 143, 217 GHz, and compare them to Planck observations in polarisation and intensity (Figure 16). Even if our model is not specifically constrained to exactly match the observations at these other frequencies, we observe a reasonable overall agreement both for polarisation and intensity. Naturally, the discrepancies between model and observation become larger as we move further away from the reference frequency. Randomly drawn temperature and spectral index fluctuations are not expected to be those of the real microwave sky.

Figure 17 shows cross-correlation for E and B power spectra between Planck observation and the modelled dust maps when we use uniform temperature and spectral index map in each layer. Figure 18 shows cross-correlation between various modelling options, showing that those models differ only at a sub-dominant level. These correlations are computed for maps smoothed to $2^{\circ}$ angular resolution over $70\%$ of sky. Each figure compares the correlation as a function of angular scale between real-sky GNILC maps, as obtained from Planck data and modelled emission. Three models are considered: A 3-D model in which the temperature and spectral index are constant in each layer, using the average values from Table 1; A 2-D model in which the 353 GHz total maps of $E$ and $B$ are simply scaled using the temperature and spectral index from the fit on the intensity maps (from Eq. 1); A 3-D model in which each layer has a different pixel-dependent map of $T$ and $\beta$ (the main model developed in the present paper).

We see an excellent correlation overall in all cases, of more than 96% at 143 GHz, and more than 99% at 217 GHz for polarisation, slightly worse for intensity, a difference that might be due to the presence of other foreground emission in the Planck foreground intensity maps – free-free, point sources, and/or CO line contamination. This shows that the large scale polarisation maps are in excellent agreement with the observations across the frequency channels where there is the best sensitivity to the CMB. The correlation decreases at higher $\ell$. This is probably due to a combination of non-vanishing noise in the GNILC maps, residuals of small scale fluctuations in the template 353 GHz E and B maps that are used to model the total polarisation, and lack of small scales in the modelled scaling law of each layer. Because GNILC, in a way ‘selects’ modes that are correlated between channels, it may also be that the correlation of the model with the GNILC data is artificially high. We postpone further investigations of this possible effect to future work.

As expected, when random fluctuations of $T$ and $\beta$ are generated in each layer, the correlation with the real observations is reduced. We also compute the average intensity emission across frequencies (Fig. 19), and note that, as expected, our multilayer model has more power at low frequency than a 2-D model with one single MBB per pixel. The same effect is observed both in intensity and in polarisation.

Finally, we can compute the level of decorrelation between polarisation maps at different frequencies as predicted by our model. Understanding this decorrelation is essential for future component separation work to detect CMB $B$ modes with component separation methods that exploit correlations between foregrounds at different frequencies, such as variants of the ILC (Tegmark et al., 2003; Eriksen et al., 2004; Delabrouille et al., 2009), CCA (Bonaldi et al., 2006), or SMICA (Delabrouille et al., 2003; Cardoso et al., 2008; Betoule et al., 2009). We generate maps with small scales and with random fluctuations of temperature and spectral index in each layer. We compute the correlation between polarisation maps (both $E$ and $B$) at 143 GHz or 217 GHz, and 353 GHz see Fig. 20) for our 3D model. The correlations obtained in both cases are ranging from 97% on small scale to close to 100% on large scales, which is larger than what is observed on real Planck maps (Planck Collaboration et al., 2017). This shows that even if our multilayer model adds a level of complexity to dust emission modelling, it can not produce a decorrelation between frequencies as strong as originally claimed in the first analysis of Planck polarisation maps (Planck Collaboration et al., 2016a). Our model, however, is compatible with the lack of evidence for such decorrelation between 217 and 353 GHz at the 0.4% level for $55\leq\ell\leq 90$ claimed in Sheehy & Slosar (2017), and predicts increased decorrelation (of the order of 1 to 2%) between 143 and 353 GHz over the same range of $\ell$. More multifrequency observations of polarised dust emission are necessary to better model dust polarisation and refine these predictions. We also note that in our model, as shown in Fig. 21 the correlations do not significantly depend on the region of sky, as they remain similar for smaller sky fractions.

8 Conclusion

We have developed a 3-dimensional model of polarised Galactic dust emission that is consistent with the large scale Planck HFI polarisation observations at 143, 217 and 353 GHz. The model is composed of six layers of emission, loosely associated with different distance ranges from the solar system as estimated from stellar extinction data. Each of these layers is assigned an integrated intensity and polarisation emission at 353 GHz, adjusted so that the sum matches the Planck observation on large scales. Small scale fluctuations are randomly generated to model the emission on scales that have not been observed with sufficient signal to noise ratio with Planck. For intensity, these random small scales extend the dust template beyond the Planck resolution of about $5^{\prime}$. For polarisation, small scale fluctuations of emission originating from the turbulence of the galactic magnetic field are randomly generated on scales smaller than $2^{\circ}$ or $2.5^{\circ}$, depending on the layer of emission considered. The level and correlations of randomly generated fluctuations are adjusted to extend the observed multivariate spectrum of the $T$, $E$ and $B$ components of the observed dust emission, assuming a 30% correlation of $T$ and $E$.

One of the primary motivation of this work is the recognition of the fact that if the parameters that define the scaling of dust emission between frequencies of observation vary across the sky, they must also vary along the line of sight. We hence assign to each layer of emission a different, pixel-dependent, scaling law in the form of a modified blackbody emission characterised, for each pixel, by a temperature and an emissivity spectral index. Observational constraints to infer the real scaling law for each layer are lacking. We hence generate random scaling laws adjusted to match on average the observed global scaling, and with fluctuations of temperature and spectral index compatible with the observed distribution of these two parameters as fitted on the Planck HFI data.

The model developed here does not pretend to be exact. The lack of multifrequency high signal-to-noise dust observations in polarisation forbids such an ambition. Nonetheless, the model provides a means to simulate a dust component that features some of the plausible complexity of the polarised dust component, while being compatible with the observed large scale polarised emission at 353 GHz and with most of the observed statistical properties of dust (temperature and polarisation power spectra, amplitude and correlation of temperature and spectral index of the best-fit modified blackbody emission). However, this model fails to predict the strong decorrelation of dust polarisation between frequency channels on small angular scales seen in Planck Collaboration et al. (2016a), a limitation that must be addressed in the future if that decorrelation is confirmed. In the mean time, we expect these simulated maps to be useful to investigate the component separation problem for future CMB polarisation surveys such as CMB-S4, PIXIE, CORE, or LiteBIRD. Simulated maps at a set of observing frequencies can be made available by the authors upon request.

Acknowledgements

We thank François Boulanger, Jan Tauber, and Mathieu Remazeilles for useful discussions and valuable comments on the first draft of this paper. Extensive use of the HEALPix pixelisation scheme (Górski et al., 2005), available from the HEALPix webpage,555http://healpix.sourceforge.net was made for this research project. We thank Douglas Finkbeiner for pointing out a mistake in the assignment of distances to the various layers in the first preprint of this paper, and an anonymous referee for many useful comments and suggestions.