# A 3D model of polarized dust emission in the Milky Way

A 3D model of polarized dust emission in the Milky Way Abstract We present a three-dimensional model of polarized 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 polarization fluctuations. The model is constrained to match observed dust emission on large scales, and match on smaller scales extrapolations of observed intensity and polarization power spectra. This model can be used to investigate the impact of plausible complexity of the polarized dust foreground emission on the analysis and interpretation of future cosmic microwave background polarization observations. polarization, dust, extinction, cosmic background radiation, cosmology: observations, diffuse radiation, submillimetre: ISM 1 INTRODUCTION Since the discovery of the cosmic microwave background (CMB) in 1965 (Penzias & Wilson 1965), significant efforts have been devoted to precise characterization of its emission, and to understanding the cosmological implications of its tiny temperature and polarization anisotropies, detected first with COBE-DMR (Smoot et al. 1992) for temperature, and with DASI for polarization (Kovac et al. 2002). Many experiments have gradually improved the measurement of CMB temperature and polarization 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 ≃5 arcmin angular scale, providing a wealth of information on the cosmological scenario. The Planck Collaboration XIII (2016c) has shown that both the CMB temperature and E-mode polarization power spectra were remarkably consistent with a spatially flat cosmology specified by six parameters, the so-called Λ cold dark matter model, with cosmic structures seeded at very early times by quantum fluctuations of space–time during an epoch of cosmic inflation. The accurate measurement of CMB polarization, 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 polarization B modes on large angular scale (see Kamionkowski & Kovetz 2016, for a review). CMB polarization also offers the opportunity to map the dark matter in the Universe that is responsible of slight distortions in polarization 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 & Planck Collaborations 2015) showed that the signal was mostly due to contamination of the observed map by polarized 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 2016 October 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 ∼ 0.001 (CORE Collaboration 2016). This requires subtracting at least 99 per cent of dust emission from the maps, or modelling 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 optimized 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; Betoule et al. 2009; Dunkley et al. 2009; Efstathiou, Gratton & Paci 2009; Bonaldi & Ricciardi 2011; Errard & Stompor 2012; Bonaldi, Ricciardi & Brown 2014; Remazeilles et al. 2016; Stompor, Errard & Poletti 2016; Remazeilles et al. 2017), using component separation methods mostly developed in the context of the analysis of WMAP and Planck intensity and polarization 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 (3D) model of the polarized dust emission, constrained by observations, that considers the spatial variation of the spectral index and of the temperature along the line of sight (LOS), 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 3D model of dust emission, which cannot be obtained without additional observations of the 3D dust, but a plausible 3D 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 polarization 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 polarization from that of dust with future multifrequency observations of polarized emission at millimetre wavelengths. This paper is organized as follows. In Section 2, we justify the need for 3D 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 3D dust data cube in temperature and polarization 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 polarization 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 3D model in Section 7. We conclude in Section 8. 2 WHY A 3D MODEL? Previous authors such as, e.g. Fauvet et al. (2011), O'Dea et al. (2012), and Vansyngel et al. (2017) have considered a 3D model of dust distribution and of the Galactic magnetic field (GMF) to model the spatial structure of dust polarization. Ghosh et al. (2017) complement this with an analysis of correlations of the direction of the GMF with the orientation of dust filaments, as traced by H i 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 3D 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 polarization (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 submillimetre 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 blackbody (MBB) of the form   $$I_{\nu }=\tau (\nu _0) \left(\frac{\nu }{\nu _0}\right)^{\beta } B_{\nu }(T),$$ (1)where Iν is the specific intensity at frequency ν and Bν(T) is the Planck blackbody function for dust at temperature T. In the frequency range we are considering, the optical depth τ(ν) scales as (ν/ν0)β, where β is a spectral index that depends on the chemical composition and structure of dust grains. Here, ν0 is a reference frequency at which a reference optical depth τ(ν0) is estimated (we use ν0 = 353 GHz throughout this paper). Using dust template observations in the Planck 353, 545, and 857 GHz channels and the IRAS 100 μm map, it is possible to fit for τ(ν0), T and β in each pixel. This fit, performed by the Planck Collaboration XI (2014), shows clear evidence for a variation across the sky of the best-fitting temperature and spectral index, with T mostly ranging from about 15 K to about 27 K and β 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 3D. Indeed, if dust properties vary across the sky, they must also vary along the LOS. This means that even if one single MBB is (empirically) a good fit to the average emission coming from a given region of the 3D 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 β(r) (in fact, a continuum, weighted by a local elementary optical depth dτ(r, ν0)):   $$I_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, \frac{\mathrm{d}\tau (r,\nu _0)}{\mathrm{d}r} \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \, B_\nu (T(r)),$$ (2)where r is the distance along the LOS and where, again, τ(r, ν0) is an optical depth at frequency ν0, T(r) is a temperature, and β(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 LOS with lower β 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 ν, which precludes fits of dust emission performed at high frequency to be valid at lower frequencies. To properly account for such LOS inhomogeneities, a 3D model of dust emission, with dust emission law variations both across and along the LOS, is needed. This 3D mixture of inhomogeneous emission would also naturally impact the polarized emission of galactic dust. The preferential alignment of elongated dust grains perpendicularly to the local magnetic field $$\boldsymbol B$$ results in a net sky polarization that is, on the plane of the sky, orthogonal to the component $$\boldsymbol B_{\perp }$$ of $$\boldsymbol B$$ that is perpendicular to the LOS. 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 polarized emission proportional to an intrinsic local polarization fraction p(r). Linear polarization Stokes parameters Q and U can be written as   $$Q_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, p(r) \frac{\mathrm{d}\tau }{\mathrm{d}r} B_\nu (T(r)) \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \cos 2\psi (r) \sin ^k \alpha (r)$$ (3)and   $$U_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, p(r) \frac{\mathrm{d}\tau }{\mathrm{d}r} B_\nu (T(r)) \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \sin 2\psi (r) \sin ^k \alpha (r),$$ (4)where, in the healpix CMB polarization convention,   $$\cos 2\psi = \frac{B_{\theta }^2-B_{\varphi }^2}{B_{\perp }^2}, \, \sin 2\psi = \frac{2B_\theta B_\varphi }{B_{\perp }^2}, \, \sin \alpha = {\frac{B_{\perp }}{B}},$$ (5)and where k is an exponent that takes into account depolarization 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, θ, and φ are spherical heliocentric coordinates. In equations (3) and (4), we recognize an overall intensity term (equal to the integrand in equation 2), multiplied by a polarization fraction p(r), an orientation term cos 2ψ(r) or sin 2ψ(r), and a geometrical term sin kα(r) that depends on the direction of the magnetic field with respect to the LOS. 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 work1 as it does not change the polarization angle on the sky, while the polarization maps will ultimately be re-normalized to match the total observed dust polarization at 353 GHz. This re-normalization somewhat corrects for possible inadequacy or inaccuracy of the assumption made for the geometrical term. Since all parameters (p, τ, T, β, ψ, and α) vary along the LOS, the total polarized emission is a superposition of emissions with different polarization angles and different emission laws. As a consequence, the polarization fraction will change with frequency (i.e. intensity and polarization have different emission laws); in addition, the polarization will rotate as a function of frequency, depending on the relative level of emission of various regions along the LOS. This polarization rotation effect would also naturally generate decorrelation of polarization emission at various frequencies. Such an effect that has been reported in Planck observations (Planck Collaboration L 2017), 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 submillimetre 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, Delabrouille & Cardoso 2011) is a component separation method that extracts from noisy multifrequency 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, Petrushev & Ward 2006; Faÿ et al. 2008; Marinucci 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-yr temperature data (Delabrouille et al. 2009). They were further used on the 7-yr and 9-yr temperature and polarization maps (Basak & Delabrouille 2012, 2013). GNILC has been used by the Planck collaboration to separate dust emission from Cosmic Infrared Background (Planck Collaboration XLVIII 2016e). 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 polarization), and of reduced contamination by the cosmic infrared background fluctuations (for intensity). However, different templates of dust emission in intensity and polarization 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 subdominant 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 ‘2D’ model of the form of equation (1) uses Planck maps of τ(ν0) and β that are obtained from a fit of the GNILC dust maps between 353 and 3000 GHz, obtained as described in Planck Collaboration XLVIII (2016e). For polarization, we apply independently GNILC on Planck 30–353 GHz E and B polarization maps (Data Release 2; Planck Collaboration I 2016b) to obtain polarized galactic emission maps in the seven polarized Planck channels. These maps are specifically produced for the present analysis, and are not part of the Planck archive. Dust-dominated E and B polarization maps at ν = 143, 217, and 353 GHz are shown in Fig. 1. The polarization maps with best dust signal-to-noise ratio are at ν = 353 GHz. The other polarization maps are not used further in our model.2 Here, the GNILC processing is mostly used as a pixel-dependent de-noising of the 353 GHz polarization map. A model that fully exploits the multifrequency information in the Planck data is postponed to future work. The needlet decomposition extends down to 5 arcmin angular resolution for intensity and 1° for polarization. Figure 1. View largeDownload slide T , E, and B maps at 353, 217, and 143 GHz, obtained with a generalized needlet ILC analysis of Planck HFI public data products. Figure 1. View largeDownload slide T , E, and B maps at 353, 217, and 143 GHz, obtained with a generalized needlet ILC analysis of Planck HFI public data products. 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–13.7 arcmin . These maps will be used to infer some information about the distribution of dust along the LOS, which will be used to generate our 3D model of polarized dust emission. 4 MULTILAYER MODELLING STRATEGY We approximate the continuous integrals of equations (2)–(4) as discrete sums over independent layers of emission, indexed by i, so that we have, for the intensity,   $$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)).$$ (6)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 Ti(p) and a spectral index β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 pixelized maps) for I, Q, and U at some reference frequency ν0, as well as scaling parameter maps T and β, all such that the total emission matches the observed sky. By ‘layer’ we mean a component, loosely associated with different distances from us, but which could equally well be a component associated with 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 six 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 Ti(p) and βi(p), i.e. a total of 5N maps, where N is the number of layers used in the model. For any N ≥ 2, we need additional data or constraints. We thus use the 3D 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 Ii coming from different layers i. We group the dust extinction maps in six ‘layers’ (shown in Fig. 2) 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 with 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. Figure 2. View largeDownload slide Maps of starlight extinction tracing the interstellar dust optical depth in shells at different distance from the Sun (maps obtained from the 3D maps of Green et al. 2015). Grey areas correspond to regions that have not been observed. Figure 2. View largeDownload slide Maps of starlight extinction tracing the interstellar dust optical depth in shells at different distance from the Sun (maps obtained from the 3D maps of Green et al. 2015). Grey areas correspond to regions that have not been observed. We then further use a 3D model of the GMF 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 3D density distribution of dust grains in the Galaxy is not. Simple models consider a galactocentric radius and height function:   $$n_d(R,z) = n_0 \exp (-R/h_R)\, {\rm sech}^2(z/h_z),$$ (7)where (R, z) are cylindrical coordinates centred at the Galactic centre, and where hR = 3 kpc and hz = 0.1 kpc. Such models cannot reproduce the observed intermediate and small-scale structure of dust emission.3 On the other hand, the maps of Green et al. (2015) trace the dust density distribution, and are directly proportional to the optical depth τ 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 fi(p) of the total opacity associated with each layer (so that ∀p, ∑i fi(p) = 1). We then construct the opacity map for each layer as the product τi(ν0) = fi τ(ν0), where τ(ν0) is the opacity at 353 GHz obtained in the Planck MBB fit. For our 3D model, we must face the practical difficulty that the maps of Green et al. (2015) do not cover the full sky (Fig. 2). 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 fi 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 τ(ν0) at 353 GHz and obtain the six maps of extinction shown in Fig. 3 . Figure 3. View largeDownload slide Full sky optical depth layers at 353 GHz, scaled to match the total 353 GHz extinction map of Planck Collaboration XLVIII (2016e). The fraction of optical depth in each layer is obtained from the maps of Green et al. (2015) where missing sky pieces in the 3D model are filled-in using symmetry arguments. Figure 3. View largeDownload slide Full sky optical depth layers at 353 GHz, scaled to match the total 353 GHz extinction map of Planck Collaboration XLVIII (2016e). The fraction of optical depth in each layer is obtained from the maps of Green et al. (2015) where missing sky pieces in the 3D model are filled-in using symmetry arguments. We then compute the corresponding brightness in a given layer by multiplying by the Planck function together with the spectral index correction (equation 8), using for this an average temperature and spectral index for each layer.4 We get, for each layer i, an initial estimate of the intensity   $$\widetilde{I}^i_{\nu } = f_i \, \tau (\nu _0) \left(\frac{\nu }{\nu _0}\right)^{\beta _i} \, B_\nu (T_i).$$ (8)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 ν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   $${I^i_{\nu _0} = \widetilde{I}^i_{\nu _0} + f_i(I_{\nu _0}- \widetilde{I}_{\nu _0})},$$ (9)and by construction we now have $${I}_{\nu _0} = \sum _i {I}^i_{\nu _0}$$. The full model across frequencies is   $$I^i_{\nu } = I^i_{\nu _0} \left(\frac{\nu }{\nu _0}\right)^{\beta _i} \, \frac{B_\nu (T_i)}{B_{\nu _0}(T_i)},$$ (10)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. (2017), 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 fi of emission in all the layers is constant over the sky. This is not compatible with a truly 3D 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 LOS 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 ≃ 30°, left of the Galactic centre, 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 with 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 with Orion is at a distance of 150–400 pc. As discussed by Rezaei et al. (2017), future analyses of the Gaia satellite data are likely to drastically improve the 3D reconstruction of Galactic dust. For this work, we use the maps of Fig. 3, 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 with 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 XLIV 2016d, 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 this work. 4.2 Polarization layers We model polarization using equations (3) and (4). Geometric terms depending on ψ and α are computed using a simple large-scale model of the 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, Mollerach & Roulet 1999; Tinyakov & Tkachev 2002), defines the radial and azimuthal field components (in Galactocentric cylindrical coordinates) as   $$B_r=B(r,\theta ,z)\sin q, \, \, \, \, \, \, B_{\theta }=-B(r,\theta ,z)\cos q,$$ (11)where q is the pitch angle of the logarithmic spiral, and where the function B(r, θ, z) is defined as   $$B(r,\theta ,z) = - B_0(r)\cos \left(\theta + \beta \log \frac{r}{r_0} \right) \exp (-|z|/z_0),$$ (12)where β = 1/tan q. We model the regular magnetic field using such a BSS parametrization, in which we consider the z-component of the GMF 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 ≤ 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° to −55° depending on the tracer used to determine it, with the most commonly cited value being around −11.5°. A possible explanation for the wide range of pitch angles determined from different data sets 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 polarized dust emission on large scales, and at high galactic latitude in particular. The simple large-scale density model of equation (7) together with the BSS large -scale magnetic field from equations (11) and (12) can be integrated following equations (2)–(4) to provide a first guess of dust intensity and polarization distribution for each layer ($$I^i_m,Q^i_m,U^i_m$$). We initially assume that the intrinsic local polarization fraction p(r) in equations (3) and (4) is constant and equal to 20 per cent. Since we already have layers of intensity emission ($$I^i_{353}$$), the polarized emission in each layer i can be generated as   $$\widetilde{Q}^{i}_{353}=\left(\frac{Q^{i}_m}{I_m^i}\right)I^i_{353}, \, \, \, \, \, \, \, \, \, \, \, \, \widetilde{U}^{i}_{353}=\left(\frac{U^{i}_m}{I_m^i}\right)I^i_{353},$$ (13)The best-fitting pitch angle q can be found minimizing some function of the difference between the simple polarization model obtained from equation (13) and the observations. We minimize the L1 norm of the difference in Q and U, summed for all the pixels at high galactic latitude:   $$G(q)= \sum _{p} \left( \left|\widetilde{Q}^{{\rm model}}_{353}(q)- Q^{{\rm obs}}_{353}\right| + \left|\widetilde{U}^{{\rm model}}_{353}(q) -U^{{\rm obs}}_{353}\right| \right),$$ (14)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 equation (13):   $$\widetilde{Q}^{{\rm model}}_{353} = \sum _{i=1}^{N} \widetilde{Q}^{i}_{353}$$ (15)and similarly for U. We find that a pitch angle of −33° provides the best fit of the GNILC maps by the BSS model at galactic latitude |b| ≥ 15°, 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°) weighted with some pixel-dependent weights Fi:   $$Q^{i}_{353}=\left(\frac{Q^{i}_m}{I_m^i}\right)I^i_{353}+F_i \left[Q^{{\rm obs}}_{353}-\sum _{j=1}^{N}\left(\frac{Q_{m}^j}{I_m^j}\right)I^j_{353}\right],$$ (16)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 Fi must be such that the polarization fraction after the redistribution of residuals does not exceed some maximum value pmax, which is a free parameter of our model, and which we pick to be 25 per cent. We fix the value of Fi as Fi = Pi/∑jPj, i.e. proportionally to the polarized dust emission fraction in each layer, unless the resulting polarization fraction exceeds pmax. When this happens, we redistribute the polarization excess in neighbouring layers. The first term in the sum on the right hand side of equation (16) is the predicted polarization of layer i, based on a polarization 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 polarized layers $$Q_m^i$$ and $$U_m^i$$ given by the large-scale model of the magnetic field while Fig. 5 shows the polarized layers after redistributing the residuals all over the former layers. After adding the small-scale 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 polarization 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 polarized dust emission on intermediate scales. In addition, the amplitude of the modelled polarized 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. Figure 4. View largeDownload slide U and Q dust emission layers obtained using a model of dust fraction in each layer based on a simple model of dust density distribution in the Galaxy, and a large-scale bi-symmetric spiral model of GMF to infer thermal dust polarization emission from dust intensity maps at 353 GHz. Figure 4. View largeDownload slide U and Q dust emission layers obtained using a model of dust fraction in each layer based on a simple model of dust density distribution in the Galaxy, and a large-scale bi-symmetric spiral model of GMF to infer thermal dust polarization emission from dust intensity maps at 353 GHz. Figure 5. View largeDownload slide U and Q dust emission layers after renormalization of the sum to match the observed sky (Planck HFI GNILC dust polarization maps at 353 GHz). Figure 5. View largeDownload slide U and Q dust emission layers after renormalization of the sum to match the observed sky (Planck HFI GNILC dust polarization maps at 353 GHz). Figure 6. View largeDownload slide U and Q layers after matching with the observed sky (as in Fig. 5), after adding random small-scale fluctuations at a level matching an extrapolation of the temperature and polarization angular power spectra and cross-spectra. Figure 6. View largeDownload slide U and Q layers after matching with the observed sky (as in Fig. 5), after adding random small-scale fluctuations at a level matching an extrapolation of the temperature and polarization angular power spectra and cross-spectra. 5 SMALL SCALES The polarization maps we have generated are normalized to match the observed dust polarization in the GNILC 353 GHz maps obtained as described in Section 3. However, the polarization GNILC maps are produced at 1° angular resolution. In the galactic plane, where the polarized signal is strong, this is the actual resolution of the GNILC map. At high galactic latitude however, the amount of polarized 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 polarized emission so far: everywhere on the sky on scales smaller than 1° (because the GNILC maps are produced at 1° 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 polarized 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, etc.). This provides us with only partial information, making it extremely difficult to generate an accurate 3D 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 μ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 equations (11) and (12) to get a first guess of the layer-dependent dust polarization, but we do not attempt to directly model the 3D 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 polarization 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 polarized dust on scales where those observations are reliable, and extend the power spectra of our maps at high l in polarization, independently for each layer, to empirically model the effect of a small-scale turbulent component of the GMF, on scales missing or noisy in the GNILC 353 GHz map. To do so, we add small-scale fluctuations independently in each layer of our model, both for intensity and for polarization. In the case of intensity, we simply fit the power spectrum of the original map in the multipole interval 30 ≤ l ≤ 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 lognormal 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 XXX (2016a), we assume a power-law dependence for EE and BB power spectra at high l, of the form Cl = A(l/lfit)α with α = −2.42 for both E and B. We use a Gaussian distribution, instead of lognormal, for polarization fields. For each layer, we fix the amplitude A and lfit to match the power spectrum of the large-scale map for that layer in the range 30 ≤ l ≤ 100. The amplitude of the small-scale fluctuations is scaled by the polarized intensity map in each layer. The randomly generated T and E harmonic coefficients are drawn with 30 per cent 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ℓ for temperature and hℓ for polarization, i.e.   $$a_{\ell m}^{T, \rm {obs}} = b_\ell \, a_{\ell m}^{T, \rm {sky}}; \quad a_{\ell m}^{E, \rm {obs}} = h_\ell \, a_{\ell m}^{E, \rm {sky}}$$ (17)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ℓm by forming   $$a_{\ell m}^{T, \rm {sim}} = a_{\ell m}^{T, \rm {obs}} + \sqrt{1-b_\ell ^2} \, a_{\ell m}^{T, \rm {rnd}}$$ (18)and similarly   $$a_{\ell m}^{E, \rm {sim}} = a_{\ell m}^{E, \rm {obs}} + \sqrt{1-h_\ell ^2} \, a_{\ell m}^{E, \rm {rnd}},$$ (19)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 arcmin for all layers in intensity, 2.5° for polarization layers 1, 2, and 3 (emission mostly at high galactic latitude), and of 2° for polarization layers 4, 5, and 6 (emission mostly near to the Galactic plane). These simulated sets of aℓm have correct $$C_\ell ^{TT}$$ and $$C_\ell ^{EE}$$, but not cross-spectrum $$C_\ell ^{TE}$$. Indeed   $$C_\ell ^{TE, {\rm sim}} = C_\ell ^{TE} \left[ b_\ell h_\ell + \sqrt{1-b_\ell ^2}\sqrt{1-h_\ell ^2} \right].$$ (20)we obtain final simulated aℓm as   $$a_{\ell m}^{\rm {final}} = \left[ C_\ell \right]^{1/2} [ C_\ell ^{\rm sim} ]^{-1/2} \, a_{\ell m}^{\rm {sim}},$$ (21)where for each ℓ, Cℓ, and $$C_\ell ^{\rm sim}$$ are 2 × 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 polarized 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 polarization fraction decreases with the polarization fraction, as seen in Fig. 9. Figure 7. View largeDownload slide First row: Full sky Q and U maps given by the BSS model. Second row: Q and U GNILC maps. Third row: Q and U total simulated maps after matching the GNILC maps on large scales and adding random small-scale fluctuations. The BSS model provides only a crude approximation of the observed dust emission. Figure 7. View largeDownload slide First row: Full sky Q and U maps given by the BSS model. Second row: Q and U GNILC maps. Third row: Q and U total simulated maps after matching the GNILC maps on large scales and adding random small-scale fluctuations. The BSS model provides only a crude approximation of the observed dust emission. Figure 8. View largeDownload slide Layers of polarized intensity ($$P = \sqrt{Q^2+U^2}$$), as modelled in our work. Figure 8. View largeDownload slide Layers of polarized intensity ($$P = \sqrt{Q^2+U^2}$$), as modelled in our work. Figure 9. View largeDownload slide Histograms of polarization fraction for each layer. We only use pixels where the polarization fraction is well defined, i.e. I(p) ≠ 0. This excludes high galactic latitude pixels for the most distant layers. Figure 9. View largeDownload slide Histograms of polarization fraction for each layer. We only use pixels where the polarization fraction is well defined, i.e. I(p) ≠ 0. This excludes high galactic latitude pixels for the most distant layers. 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-scale 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°, 50°) 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. Figure 10. View largeDownload slide T, E, B power spectra for each layer. The first three rows also display the power spectra for 75 per cent and 25 per cent of the sky. Figure 10. View largeDownload slide T, E, B power spectra for each layer. The first three rows also display the power spectra for 75 per cent and 25 per cent of the sky. Figure 11. View largeDownload slide TT, EE, BB, TE power spectra of both GNILC maps and of simulated maps including small-scale fluctuations. Figure 11. View largeDownload slide TT, EE, BB, TE power spectra of both GNILC maps and of simulated maps including small-scale fluctuations. Figure 12. View largeDownload slide Modelled E and B modes maps at 353 GHz, after adding small-scale fluctuations, adding-up six layers of emission (see the text). Figure 12. View largeDownload slide Modelled E and B modes maps at 353 GHz, after adding small-scale fluctuations, adding-up six layers of emission (see the text). Figure 13. View largeDownload slide Observed and modelled E and B modes maps at 353 GHz – detail around (l, b) = (0°, 50°). Top row: T, E, and B modes, observed with Planck after GNILC processing; Bottom row: Modelled T, E, and B modes at Nside= 512, after adding small-scale fluctuations, adding-up six layers of emission. Figure 13. View largeDownload slide Observed and modelled E and B modes maps at 353 GHz – detail around (l, b) = (0°, 50°). Top row: T, E, and B modes, observed with Planck after GNILC processing; Bottom row: Modelled T, E, and B modes at Nside= 512, after adding small-scale fluctuations, adding-up six layers of emission. Figure 14. View largeDownload slide E and B power spectra of both GNILC maps and GNILC maps + small-scale fluctuations at 143 and 217 GHz. Figure 14. View largeDownload slide E and B power spectra of both GNILC maps and GNILC maps + small-scale fluctuations at 143 and 217 GHz. 6 SCALING LAWS We now need a prescription for scaling the 353 GHz polarized 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 MBB with pixel-dependent temperature and spectral index), but now, we must define as many templates of T(p) and β(p) as there are layers in our model, i.e. six maps of T and six maps of β. A complete description of the temperature and the spectral index distribution in 3D 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 polarization maps, scaled to other frequencies than 353 GHz, should match the observed polarization 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–β scatter plot. With only 353 GHz polarization 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 XI (2014), which we assume to hold for all of the I, Q, and U Stokes parameters. We compute, for each layer i:   \begin{eqnarray} T_{{\rm avg}}^i&=&\sum _{p} w^i_T(p) T_d(p), \nonumber\\ \beta _{{\rm avg}}^i&=&\sum _{p} w^i_\beta (p) \beta _d(p), \end{eqnarray} (22)where Td(p) and βd(p) are the best-fitting 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_{{\rm avg}}^i$$ and spectral index $$\beta _{{\rm avg}}^i$$,   $$w^i_T(p) = w^i_\beta (p) = f_i(p),$$ (23)i.e. we empirically weight the maps by the pixel-dependent fraction fi(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_{{\rm avg}}^i$$ and $$\beta _{{\rm avg}}^i$$ are constant across the sky in a given layer. This however implements only a variability of the physical parameters T and β 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 β (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_{{\rm avg}}^i$$ and $$\beta _{{\rm avg}}^i$$ following the auto and cross-spectra of the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). To take into account the non-Gaussianity of the distribution of T and β, 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 β. 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. Figure 15. View largeDownload slide Top left: Power spectra used to draw random realizations of temperature and spectral index maps (note the negative sign of the Tβ cross-spectrum). Bottom left: Scatter plot of T and β for a pair of random maps (right), showing an overall anticorrelation and the same general behaviour as observed by Planck Collaboration XI (2014) on Planck observations (see their fig. 16). Right: Maps of randomly generated temperature and spectral index for the first layer, with Tavg = 19.10, σT = 2.059, βavg = 1.627, σβ = 0.209. Figure 15. View largeDownload slide Top left: Power spectra used to draw random realizations of temperature and spectral index maps (note the negative sign of the Tβ cross-spectrum). Bottom left: Scatter plot of T and β for a pair of random maps (right), showing an overall anticorrelation and the same general behaviour as observed by Planck Collaboration XI (2014) on Planck observations (see their fig. 16). Right: Maps of randomly generated temperature and spectral index for the first layer, with Tavg = 19.10, σT = 2.059, βavg = 1.627, σβ = 0.209. 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 ν = 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 β 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 β 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 β in each layer for one single realization of the simulation, compared with those from the Planck MBB fit and the fit performed on this realization. The average values from several simulations are in good agreement with those of the Planck MMB fit. Table 1. Averages and standard deviation values of temperature and spectral index in each layer, for a simulation with 6.87 arcmin 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  Tavg  19.10  18.96  18.98  19.35  19.23  20.05  σT  2.059  2.100  2.022  2.076  2.117  2.069  βavg  1.627  1.628  1.598  1.538  1.513  1.689  σβ  0.209  0.210  0.207  0.208  0.202  0.204    $$T_{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{T}^{{\rm {MMB}}}$$  $$\beta _{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{\beta }^{{\rm {MMB}}}$$      Planck fit  19.396  1.247  1.598  0.126      Simul. fit  19.389  1.253  1.598  0.135      Layer  1  2  3  4  5  6  Tavg  19.10  18.96  18.98  19.35  19.23  20.05  σT  2.059  2.100  2.022  2.076  2.117  2.069  βavg  1.627  1.628  1.598  1.538  1.513  1.689  σβ  0.209  0.210  0.207  0.208  0.202  0.204    $$T_{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{T}^{{\rm {MMB}}}$$  $$\beta _{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{\beta }^{{\rm {MMB}}}$$      Planck fit  19.396  1.247  1.598  0.126      Simul. fit  19.389  1.253  1.598  0.135      View Large 7 VALIDATION AND PREDICTIONS We use our model to generate maps of polarized dust emission at 143, 217 GHz, and compare them to Planck observations in polarization and intensity (Fig. 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 polarization 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 16. View largeDownload slide GNILC maps both in intensity and polarization are shown in the first and the fourth row (subindex G), while maps obtained using our 3D model are shown in the second and the fifth row (subindex m). The differences between them are also shown in the third and sixth row. Note the different colour scales for difference maps. Figure 16. View largeDownload slide GNILC maps both in intensity and polarization are shown in the first and the fourth row (subindex G), while maps obtained using our 3D model are shown in the second and the fifth row (subindex m). The differences between them are also shown in the third and sixth row. Note the different colour scales for difference maps. Fig. 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. Fig. 18 shows cross-correlation between various modelling options, showing that those models differ only at a subdominant level. These correlations are computed for maps smoothed to 2° angular resolution over 70 per cent 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 3D model in which the temperature and spectral index are constant in each layer, using the average values from Table 1; A 2D 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 equation 1); A 3D model in which each layer has a different pixel-dependent map of T and β (the main model developed in this paper). Figure 17. View largeDownload slide Cross-correlation between simulations and observations for T, E, and B power spectra at 143 and 217 GHz 70 per cent of sky. We show in blue and green the correlation in intensity and polarization between maps generated with our model and the observations. While blue curves are computed using one single value of temperature and spectral index per layer, green curves consider one template per layer, with fluctuations of the temperature and the spectral index (model b). Red curves show the correlation between the observed polarized sky maps and maps obtained from a 2D model, i.e. one single template for temperature and spectral index from the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). Figure 17. View largeDownload slide Cross-correlation between simulations and observations for T, E, and B power spectra at 143 and 217 GHz 70 per cent of sky. We show in blue and green the correlation in intensity and polarization between maps generated with our model and the observations. While blue curves are computed using one single value of temperature and spectral index per layer, green curves consider one template per layer, with fluctuations of the temperature and the spectral index (model b). Red curves show the correlation between the observed polarized sky maps and maps obtained from a 2D model, i.e. one single template for temperature and spectral index from the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). Figure 18. View largeDownload slide Cross-correlations of T, E, and B between various modelling options. Differences between these models in polarization are at the level of a few per cent at most for ℓ ≤ 100. Figure 18. View largeDownload slide Cross-correlations of T, E, and B between various modelling options. Differences between these models in polarization are at the level of a few per cent at most for ℓ ≤ 100. We see an excellent correlation overall in all cases, of more than 96 per cent at 143 GHz, and more than 99 per cent at 217 GHz for polarization, 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 polarization 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 ℓ. 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 polarization, 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 β 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 2D model with one single MBB per pixel. The same effect is observed both in intensity and in polarization. Figure 19. View largeDownload slide Left: Average total sky emission in intensity for our 3D multi-MBB model as compared to a ‘2D’ model with one single MBB per pixel. Both have the same intensity at 353 GHz by construction. The 3D model has flatter emission law at low frequency, an effect that originates from the increasing importance at low frequency of components with flatter spectral index that may be subdominant at higher frequency where the emission is dominated by hotter components. Right: Ratio of the average emission law of the 3D model and the 2D model, for both intensity and polarized intensity. Figure 19. View largeDownload slide Left: Average total sky emission in intensity for our 3D multi-MBB model as compared to a ‘2D’ model with one single MBB per pixel. Both have the same intensity at 353 GHz by construction. The 3D model has flatter emission law at low frequency, an effect that originates from the increasing importance at low frequency of components with flatter spectral index that may be subdominant at higher frequency where the emission is dominated by hotter components. Right: Ratio of the average emission law of the 3D model and the 2D model, for both intensity and polarized intensity. Finally, we can compute the level of decorrelation between polarization 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, de Oliveira-Costa & Hamilton 2003; Eriksen et al. 2004; Delabrouille et al. 2009), CCA (Bonaldi et al. 2006), or SMICA (Delabrouille, Cardoso & Patanchon 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 polarization maps (both E and B) at 143 or 217, and 353 GHz (see Fig. 20) for our 3D model. The correlations obtained in both cases are ranging from 97 per cent on small scale to close to 100 per cent on large scales, which is larger than what is observed on real Planck maps (Planck Collaboration L 2017). This shows that even if our multilayer model adds a level of complexity to dust emission modelling, it cannot produce a decorrelation between frequencies as strong as originally claimed in the first analysis of Planck polarization maps (Planck Collaboration L 2017). Our model, however, is compatible with the lack of evidence for such decorrelation between 217 and 353 GHz at the 0.4 per cent level for 55 ≤ ℓ ≤ 90 claimed in Sheehy & Slosar (2017), and predicts increased decorrelation (of the order of 1–2 per cent) between 143 and 353 GHz over the same range of ℓ. More multifrequency observations of polarized dust emission are necessary to better model dust polarization 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. Figure 20. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over 70 per cent of sky. Figure 20. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over 70 per cent of sky. Figure 21. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over different sky fractions. Figure 21. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over different sky fractions. 8 CONCLUSION We have developed a 3D model of polarized Galactic dust emission that is consistent with the large-scale Planck HFI polarization 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 polarization 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 arcmin . For polarization, small-scale fluctuations of emission originating from the turbulence of the GMF are randomly generated on scales smaller than 2° or 2.5°, 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 per cent 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 LOS. We hence assign to each layer of emission a different, pixel-dependent, scaling law in the form of a MBB emission characterized, 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 polarization forbids such an ambition. None the less, the model provides a means to simulate a dust component that features some of the plausible complexity of the polarized dust component, while being compatible with the observed large-scale polarized emission at 353 GHz and with most of the observed statistical properties of dust (temperature and polarization power spectra, amplitude and correlation of temperature and spectral index of the best-fitting MBB emission). However, this model fails to predict the strong decorrelation of dust polarization between frequency channels on small angular scales seen in Planck Collaboration L (2017), a limitation that must be addressed in the future if that decorrelation is confirmed. In the meantime, we expect these simulated maps to be useful to investigate the component separation problem for future CMB polarization 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 pixelization scheme (Górski et al. 2005), available from the healpix webpage,5 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. Footnotes 1 Nor does the specific analytic form of the depolarization function. 2 After 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. 3 However, 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 polarization fraction of dust emission in each layer. 4 We postpone to Section 6 the discussion of temperature and spectral index maps. In equation (8), we use average values given in Table 1. 5 http://healpix.sourceforge.net REFERENCES Abazajian K. N. et al.  , 2016, preprint (arXiv:1610.02743) André P. et al.  , 2014, J. Cosmol. Astropart. Phys. , 2, 006 CrossRef Search ADS   Basak S., Delabrouille J., 2012, MNRAS , 419, 1163 CrossRef Search ADS   Basak S., Delabrouille J., 2013, MNRAS , 435, 18 CrossRef Search ADS   Bennett C. L. et al.  , 2013, ApJS , 208, 20 CrossRef Search ADS   Benoit A. et al.  , 2002, Astropart. Phys. , 17, 101 CrossRef Search ADS   Betoule M., Pierpaoli E., Delabrouille J., Le Jeune M., Cardoso J.-F., 2009, A&A , 503, 691 CrossRef Search ADS   BICEP2 Collaboration, 2014, Phys. Rev. Lett. , 112, 241101 CrossRef Search ADS PubMed  BICEP2/Keck & Planck Collaborations, 2015, Phys. Rev. Lett. , 114, 101301 CrossRef Search ADS PubMed  Bonaldi A., Ricciardi S., 2011, MNRAS , 414, 615 CrossRef Search ADS   Bonaldi A., Bedini L., Salerno E., Baccigalupi C., de Zotti G., 2006, MNRAS , 373, 271 CrossRef Search ADS   Bonaldi A., Ricciardi S., Brown M. L., 2014, MNRAS , 444, 1034 CrossRef Search ADS   Cardoso J.-F., Le Jeune M., Delabrouille J., Betoule M., Patanchon G., 2008, IEEE J. Sel. Top. Signal Process. , 2, 735 CrossRef Search ADS   Challinor A. et al.  , 2017, preprint (arXiv:1707.02259) CORE Collaboration, 2016, preprint (arXiv:1612.08270) Das S. et al.  , 2014, J. Cosmol. Astropart. Phys. , 4, 014 CrossRef Search ADS   de Bernardis P. et al.  , 2000, Nature , 404, 955 CrossRef Search ADS PubMed  Delabrouille J., Cardoso J.-F., 2009, in Martínez V. J., Saar E., Martínez-González E., Pons-Bordería M.-J., eds, Lecture Notes in Physics, Vol. 665, Data Analysis in Cosmology . Springer-Verlag, Berlin, p. 159 Delabrouille J., Cardoso J.-F., Patanchon G., 2003, MNRAS , 346, 1089 CrossRef Search ADS   Delabrouille J., Cardoso J.-F., Le Jeune M., Betoule M., Fay G., Guilloux F., 2009, A&A , 493, 835 CrossRef Search ADS   Delabrouille J. et al.  , 2013, A&A , 553, A96 CrossRef Search ADS   Delabrouille J. et al.  , 2017, preprint (arXiv:1706.04516) Dunkley J. et al.  , 2009, in Dodelson S. et al.  , eds, AIP Conf. Ser. Vol. 1141, CMB Polarization Workshop: Theory and Foregrounds: CMBPol Mission Concept Study . Am. Inst. Phys., New York, p. 222 Efstathiou G., Gratton S., Paci F., 2009, MNRAS , 397, 1355 CrossRef Search ADS   Eriksen H. K., Banday A. J., Górski K. M., Lilje P. B., 2004, ApJ , 612, 633 CrossRef Search ADS   Errard J., Stompor R., 2012, Phys. Rev. D , 85, 083006 CrossRef Search ADS   Fauvet L. et al.  , 2011, A&A , 526, A145 CrossRef Search ADS   Faÿ G., Guilloux F., Betoule M., Cardoso J.-F., Delabrouille J., Le Jeune M., 2008, Phys. Rev. D , 78, 083013 CrossRef Search ADS   Ghosh T. et al.  , 2017, A&A , 601, A71 CrossRef Search ADS   Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ , 622, 759 CrossRef Search ADS   Green G. M. et al.  , 2015, ApJ , 810, 25 CrossRef Search ADS   Han J. L., 2017, ARA&A , 55, 111 CrossRef Search ADS   Han J. L., Qiao G. J., 1993, Acta Astrophys. Sin. , 13, 385 Hanany S. et al.  , 2000, ApJ , 545, L5 CrossRef Search ADS   Harari D., Mollerach S., Roulet E., 1999, J. High Energy Phys. , 8, 022 CrossRef Search ADS   Haverkorn M., 2015, in Lazarian A., de Gouveia Dal Pino E. M., Melioli C., eds, Astrophysics and Space Science Library, Vol. 407, Magnetic Fields in Diffuse Media . Springer-Verlag Berlin, p. 483 Jaffe T. R., Leahy J. P., Banday A. J., Leach S. M., Lowe S. R., Wilkinson A., 2010, MNRAS , 401, 1013 CrossRef Search ADS   Jones W. C. et al.  , 2006, ApJ , 647, 823 CrossRef Search ADS   Kamionkowski M., Kovetz E. D., 2016, ARA&A , 54, 227 CrossRef Search ADS   Kogut A. et al.  , 2011, J. Cosmol. Astropart. Phys. , 7, 025 CrossRef Search ADS   Kovac J. M., Leitch E. M., Pryke C., Carlstrom J. E., Halverson N. W., Holzapfel W. L., 2002, Nature , 420, 772 CrossRef Search ADS PubMed  Lallement R., Snowden S., Kuntz K. D., Dame T. M., Koutroumpa D., Grenier I., Casandjian J. M., 2016, A&A , 595, A131 CrossRef Search ADS   Leach S. M. et al.  , 2008, A&A , 491, 597 CrossRef Search ADS   Lewis A., Challinor A., 2006, Phys. Rep. , 429, 1 CrossRef Search ADS   Marinucci D. et al.  , 2008, MNRAS , 383, 539 CrossRef Search ADS   Matsumura T. et al.  , 2014, J. Low Temp. Phys. , 176, 733 CrossRef Search ADS   Narcowich F., Petrushev P., Ward J., 2006, SIAM J. Math. Anal. , 38, 574 CrossRef Search ADS   O'Dea D. T., Clark C. N., Contaldi C. R., MacTavish C. J., 2012, MNRAS , 419, 1795 CrossRef Search ADS   Penzias A. A., Wilson R. W., 1965, ApJ , 142, 419 CrossRef Search ADS   Planck Collaboration XI, 2014, A&A , 571, A11 CrossRef Search ADS   Planck Collaboration XXX, 2016a, A&A , 586, A133 CrossRef Search ADS   Planck Collaboration I, 2016b, A&A , 594, A1 CrossRef Search ADS   Planck Collaboration XIII, 2016c, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration XLIV, 2016d, A&A , 596, A105 CrossRef Search ADS   Planck Collaboration XLVIII, 2016e, A&A , 596, A109 CrossRef Search ADS   Planck Collaboration L, 2017, A&A , 599, A51 CrossRef Search ADS   Reichardt C. L. et al.  , 2009, ApJ , 694, 1200 CrossRef Search ADS   Remazeilles M., Delabrouille J., Cardoso J.-F., 2011, MNRAS , 418, 467 CrossRef Search ADS   Remazeilles M., Dickinson C., Eriksen H. K. K., Wehus I. K., 2016, MNRAS , 458, 2032 CrossRef Search ADS   Remazeilles M. et al.  , 2017, preprint (arXiv:1704.04501) Rezaei Kh. S., Bailer-Jones C. A. L., Hanson R. J., Fouesneau M., 2017, A&A , 598, A125 CrossRef Search ADS   Sheehy C., Slosar A., 2017, preprint (arXiv:1709.09729) Smoot G. F. et al.  , 1992, ApJ , 396, L1 CrossRef Search ADS   Sofue Y., Fujimoto M., 1983, ApJ , 265, 722 CrossRef Search ADS   Stanev T., 1997, ApJ , 479, 290 CrossRef Search ADS   Stompor R., Errard J., Poletti D., 2016, Phys. Rev. D , 94, 083526 CrossRef Search ADS   Story K. T. et al.  , 2013, ApJ , 779, 86 CrossRef Search ADS   Tassis K., Pavlidou V., 2015, MNRAS , 451, L90 CrossRef Search ADS   Tauber J. A. et al.  , 2010, A&A , 520, A1 CrossRef Search ADS   Tegmark M., de Oliveira-Costa A., Hamilton A. J., 2003, Phys. Rev. D , 68, 123523 CrossRef Search ADS   The COrE Collaboration, 2011, preprint (arXiv:1102.2181) Tinyakov P. G., Tkachev I. I., 2002, Astropart. Phys. , 18, 165 CrossRef Search ADS   Tucci M., Martínez-González E., Vielva P., Delabrouille J., 2005, MNRAS , 360, 935 CrossRef Search ADS   Vansyngel F. et al.  , 2017, A&A , 603, A62 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# A 3D model of polarized dust emission in the Milky Way

, Volume 476 (1) – May 1, 2018
21 pages

/lp/ou_press/a-3d-model-of-polarized-dust-emission-in-the-milky-way-00Xuq0UtRl
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty204
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present a three-dimensional model of polarized 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 polarization fluctuations. The model is constrained to match observed dust emission on large scales, and match on smaller scales extrapolations of observed intensity and polarization power spectra. This model can be used to investigate the impact of plausible complexity of the polarized dust foreground emission on the analysis and interpretation of future cosmic microwave background polarization observations. polarization, dust, extinction, cosmic background radiation, cosmology: observations, diffuse radiation, submillimetre: ISM 1 INTRODUCTION Since the discovery of the cosmic microwave background (CMB) in 1965 (Penzias & Wilson 1965), significant efforts have been devoted to precise characterization of its emission, and to understanding the cosmological implications of its tiny temperature and polarization anisotropies, detected first with COBE-DMR (Smoot et al. 1992) for temperature, and with DASI for polarization (Kovac et al. 2002). Many experiments have gradually improved the measurement of CMB temperature and polarization 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 ≃5 arcmin angular scale, providing a wealth of information on the cosmological scenario. The Planck Collaboration XIII (2016c) has shown that both the CMB temperature and E-mode polarization power spectra were remarkably consistent with a spatially flat cosmology specified by six parameters, the so-called Λ cold dark matter model, with cosmic structures seeded at very early times by quantum fluctuations of space–time during an epoch of cosmic inflation. The accurate measurement of CMB polarization, 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 polarization B modes on large angular scale (see Kamionkowski & Kovetz 2016, for a review). CMB polarization also offers the opportunity to map the dark matter in the Universe that is responsible of slight distortions in polarization 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 & Planck Collaborations 2015) showed that the signal was mostly due to contamination of the observed map by polarized 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 2016 October 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 ∼ 0.001 (CORE Collaboration 2016). This requires subtracting at least 99 per cent of dust emission from the maps, or modelling 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 optimized 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; Betoule et al. 2009; Dunkley et al. 2009; Efstathiou, Gratton & Paci 2009; Bonaldi & Ricciardi 2011; Errard & Stompor 2012; Bonaldi, Ricciardi & Brown 2014; Remazeilles et al. 2016; Stompor, Errard & Poletti 2016; Remazeilles et al. 2017), using component separation methods mostly developed in the context of the analysis of WMAP and Planck intensity and polarization 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 (3D) model of the polarized dust emission, constrained by observations, that considers the spatial variation of the spectral index and of the temperature along the line of sight (LOS), 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 3D model of dust emission, which cannot be obtained without additional observations of the 3D dust, but a plausible 3D 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 polarization 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 polarization from that of dust with future multifrequency observations of polarized emission at millimetre wavelengths. This paper is organized as follows. In Section 2, we justify the need for 3D 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 3D dust data cube in temperature and polarization 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 polarization 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 3D model in Section 7. We conclude in Section 8. 2 WHY A 3D MODEL? Previous authors such as, e.g. Fauvet et al. (2011), O'Dea et al. (2012), and Vansyngel et al. (2017) have considered a 3D model of dust distribution and of the Galactic magnetic field (GMF) to model the spatial structure of dust polarization. Ghosh et al. (2017) complement this with an analysis of correlations of the direction of the GMF with the orientation of dust filaments, as traced by H i 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 3D 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 polarization (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 submillimetre 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 blackbody (MBB) of the form   $$I_{\nu }=\tau (\nu _0) \left(\frac{\nu }{\nu _0}\right)^{\beta } B_{\nu }(T),$$ (1)where Iν is the specific intensity at frequency ν and Bν(T) is the Planck blackbody function for dust at temperature T. In the frequency range we are considering, the optical depth τ(ν) scales as (ν/ν0)β, where β is a spectral index that depends on the chemical composition and structure of dust grains. Here, ν0 is a reference frequency at which a reference optical depth τ(ν0) is estimated (we use ν0 = 353 GHz throughout this paper). Using dust template observations in the Planck 353, 545, and 857 GHz channels and the IRAS 100 μm map, it is possible to fit for τ(ν0), T and β in each pixel. This fit, performed by the Planck Collaboration XI (2014), shows clear evidence for a variation across the sky of the best-fitting temperature and spectral index, with T mostly ranging from about 15 K to about 27 K and β 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 3D. Indeed, if dust properties vary across the sky, they must also vary along the LOS. This means that even if one single MBB is (empirically) a good fit to the average emission coming from a given region of the 3D 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 β(r) (in fact, a continuum, weighted by a local elementary optical depth dτ(r, ν0)):   $$I_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, \frac{\mathrm{d}\tau (r,\nu _0)}{\mathrm{d}r} \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \, B_\nu (T(r)),$$ (2)where r is the distance along the LOS and where, again, τ(r, ν0) is an optical depth at frequency ν0, T(r) is a temperature, and β(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 LOS with lower β 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 ν, which precludes fits of dust emission performed at high frequency to be valid at lower frequencies. To properly account for such LOS inhomogeneities, a 3D model of dust emission, with dust emission law variations both across and along the LOS, is needed. This 3D mixture of inhomogeneous emission would also naturally impact the polarized emission of galactic dust. The preferential alignment of elongated dust grains perpendicularly to the local magnetic field $$\boldsymbol B$$ results in a net sky polarization that is, on the plane of the sky, orthogonal to the component $$\boldsymbol B_{\perp }$$ of $$\boldsymbol B$$ that is perpendicular to the LOS. 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 polarized emission proportional to an intrinsic local polarization fraction p(r). Linear polarization Stokes parameters Q and U can be written as   $$Q_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, p(r) \frac{\mathrm{d}\tau }{\mathrm{d}r} B_\nu (T(r)) \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \cos 2\psi (r) \sin ^k \alpha (r)$$ (3)and   $$U_{\nu } = \int _{0}^{\infty } \mathrm{d}r \, p(r) \frac{\mathrm{d}\tau }{\mathrm{d}r} B_\nu (T(r)) \left(\frac{\nu }{\nu _0}\right)^{\beta (r)} \sin 2\psi (r) \sin ^k \alpha (r),$$ (4)where, in the healpix CMB polarization convention,   $$\cos 2\psi = \frac{B_{\theta }^2-B_{\varphi }^2}{B_{\perp }^2}, \, \sin 2\psi = \frac{2B_\theta B_\varphi }{B_{\perp }^2}, \, \sin \alpha = {\frac{B_{\perp }}{B}},$$ (5)and where k is an exponent that takes into account depolarization 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, θ, and φ are spherical heliocentric coordinates. In equations (3) and (4), we recognize an overall intensity term (equal to the integrand in equation 2), multiplied by a polarization fraction p(r), an orientation term cos 2ψ(r) or sin 2ψ(r), and a geometrical term sin kα(r) that depends on the direction of the magnetic field with respect to the LOS. 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 work1 as it does not change the polarization angle on the sky, while the polarization maps will ultimately be re-normalized to match the total observed dust polarization at 353 GHz. This re-normalization somewhat corrects for possible inadequacy or inaccuracy of the assumption made for the geometrical term. Since all parameters (p, τ, T, β, ψ, and α) vary along the LOS, the total polarized emission is a superposition of emissions with different polarization angles and different emission laws. As a consequence, the polarization fraction will change with frequency (i.e. intensity and polarization have different emission laws); in addition, the polarization will rotate as a function of frequency, depending on the relative level of emission of various regions along the LOS. This polarization rotation effect would also naturally generate decorrelation of polarization emission at various frequencies. Such an effect that has been reported in Planck observations (Planck Collaboration L 2017), 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 submillimetre 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, Delabrouille & Cardoso 2011) is a component separation method that extracts from noisy multifrequency 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, Petrushev & Ward 2006; Faÿ et al. 2008; Marinucci 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-yr temperature data (Delabrouille et al. 2009). They were further used on the 7-yr and 9-yr temperature and polarization maps (Basak & Delabrouille 2012, 2013). GNILC has been used by the Planck collaboration to separate dust emission from Cosmic Infrared Background (Planck Collaboration XLVIII 2016e). 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 polarization), and of reduced contamination by the cosmic infrared background fluctuations (for intensity). However, different templates of dust emission in intensity and polarization 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 subdominant 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 ‘2D’ model of the form of equation (1) uses Planck maps of τ(ν0) and β that are obtained from a fit of the GNILC dust maps between 353 and 3000 GHz, obtained as described in Planck Collaboration XLVIII (2016e). For polarization, we apply independently GNILC on Planck 30–353 GHz E and B polarization maps (Data Release 2; Planck Collaboration I 2016b) to obtain polarized galactic emission maps in the seven polarized Planck channels. These maps are specifically produced for the present analysis, and are not part of the Planck archive. Dust-dominated E and B polarization maps at ν = 143, 217, and 353 GHz are shown in Fig. 1. The polarization maps with best dust signal-to-noise ratio are at ν = 353 GHz. The other polarization maps are not used further in our model.2 Here, the GNILC processing is mostly used as a pixel-dependent de-noising of the 353 GHz polarization map. A model that fully exploits the multifrequency information in the Planck data is postponed to future work. The needlet decomposition extends down to 5 arcmin angular resolution for intensity and 1° for polarization. Figure 1. View largeDownload slide T , E, and B maps at 353, 217, and 143 GHz, obtained with a generalized needlet ILC analysis of Planck HFI public data products. Figure 1. View largeDownload slide T , E, and B maps at 353, 217, and 143 GHz, obtained with a generalized needlet ILC analysis of Planck HFI public data products. 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–13.7 arcmin . These maps will be used to infer some information about the distribution of dust along the LOS, which will be used to generate our 3D model of polarized dust emission. 4 MULTILAYER MODELLING STRATEGY We approximate the continuous integrals of equations (2)–(4) as discrete sums over independent layers of emission, indexed by i, so that we have, for the intensity,   $$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)).$$ (6)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 Ti(p) and a spectral index β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 pixelized maps) for I, Q, and U at some reference frequency ν0, as well as scaling parameter maps T and β, all such that the total emission matches the observed sky. By ‘layer’ we mean a component, loosely associated with different distances from us, but which could equally well be a component associated with 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 six 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 Ti(p) and βi(p), i.e. a total of 5N maps, where N is the number of layers used in the model. For any N ≥ 2, we need additional data or constraints. We thus use the 3D 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 Ii coming from different layers i. We group the dust extinction maps in six ‘layers’ (shown in Fig. 2) 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 with 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. Figure 2. View largeDownload slide Maps of starlight extinction tracing the interstellar dust optical depth in shells at different distance from the Sun (maps obtained from the 3D maps of Green et al. 2015). Grey areas correspond to regions that have not been observed. Figure 2. View largeDownload slide Maps of starlight extinction tracing the interstellar dust optical depth in shells at different distance from the Sun (maps obtained from the 3D maps of Green et al. 2015). Grey areas correspond to regions that have not been observed. We then further use a 3D model of the GMF 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 3D density distribution of dust grains in the Galaxy is not. Simple models consider a galactocentric radius and height function:   $$n_d(R,z) = n_0 \exp (-R/h_R)\, {\rm sech}^2(z/h_z),$$ (7)where (R, z) are cylindrical coordinates centred at the Galactic centre, and where hR = 3 kpc and hz = 0.1 kpc. Such models cannot reproduce the observed intermediate and small-scale structure of dust emission.3 On the other hand, the maps of Green et al. (2015) trace the dust density distribution, and are directly proportional to the optical depth τ 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 fi(p) of the total opacity associated with each layer (so that ∀p, ∑i fi(p) = 1). We then construct the opacity map for each layer as the product τi(ν0) = fi τ(ν0), where τ(ν0) is the opacity at 353 GHz obtained in the Planck MBB fit. For our 3D model, we must face the practical difficulty that the maps of Green et al. (2015) do not cover the full sky (Fig. 2). 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 fi 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 τ(ν0) at 353 GHz and obtain the six maps of extinction shown in Fig. 3 . Figure 3. View largeDownload slide Full sky optical depth layers at 353 GHz, scaled to match the total 353 GHz extinction map of Planck Collaboration XLVIII (2016e). The fraction of optical depth in each layer is obtained from the maps of Green et al. (2015) where missing sky pieces in the 3D model are filled-in using symmetry arguments. Figure 3. View largeDownload slide Full sky optical depth layers at 353 GHz, scaled to match the total 353 GHz extinction map of Planck Collaboration XLVIII (2016e). The fraction of optical depth in each layer is obtained from the maps of Green et al. (2015) where missing sky pieces in the 3D model are filled-in using symmetry arguments. We then compute the corresponding brightness in a given layer by multiplying by the Planck function together with the spectral index correction (equation 8), using for this an average temperature and spectral index for each layer.4 We get, for each layer i, an initial estimate of the intensity   $$\widetilde{I}^i_{\nu } = f_i \, \tau (\nu _0) \left(\frac{\nu }{\nu _0}\right)^{\beta _i} \, B_\nu (T_i).$$ (8)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 ν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   $${I^i_{\nu _0} = \widetilde{I}^i_{\nu _0} + f_i(I_{\nu _0}- \widetilde{I}_{\nu _0})},$$ (9)and by construction we now have $${I}_{\nu _0} = \sum _i {I}^i_{\nu _0}$$. The full model across frequencies is   $$I^i_{\nu } = I^i_{\nu _0} \left(\frac{\nu }{\nu _0}\right)^{\beta _i} \, \frac{B_\nu (T_i)}{B_{\nu _0}(T_i)},$$ (10)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. (2017), 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 fi of emission in all the layers is constant over the sky. This is not compatible with a truly 3D 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 LOS 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 ≃ 30°, left of the Galactic centre, 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 with 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 with Orion is at a distance of 150–400 pc. As discussed by Rezaei et al. (2017), future analyses of the Gaia satellite data are likely to drastically improve the 3D reconstruction of Galactic dust. For this work, we use the maps of Fig. 3, 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 with 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 XLIV 2016d, 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 this work. 4.2 Polarization layers We model polarization using equations (3) and (4). Geometric terms depending on ψ and α are computed using a simple large-scale model of the 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, Mollerach & Roulet 1999; Tinyakov & Tkachev 2002), defines the radial and azimuthal field components (in Galactocentric cylindrical coordinates) as   $$B_r=B(r,\theta ,z)\sin q, \, \, \, \, \, \, B_{\theta }=-B(r,\theta ,z)\cos q,$$ (11)where q is the pitch angle of the logarithmic spiral, and where the function B(r, θ, z) is defined as   $$B(r,\theta ,z) = - B_0(r)\cos \left(\theta + \beta \log \frac{r}{r_0} \right) \exp (-|z|/z_0),$$ (12)where β = 1/tan q. We model the regular magnetic field using such a BSS parametrization, in which we consider the z-component of the GMF 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 ≤ 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° to −55° depending on the tracer used to determine it, with the most commonly cited value being around −11.5°. A possible explanation for the wide range of pitch angles determined from different data sets 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 polarized dust emission on large scales, and at high galactic latitude in particular. The simple large-scale density model of equation (7) together with the BSS large -scale magnetic field from equations (11) and (12) can be integrated following equations (2)–(4) to provide a first guess of dust intensity and polarization distribution for each layer ($$I^i_m,Q^i_m,U^i_m$$). We initially assume that the intrinsic local polarization fraction p(r) in equations (3) and (4) is constant and equal to 20 per cent. Since we already have layers of intensity emission ($$I^i_{353}$$), the polarized emission in each layer i can be generated as   $$\widetilde{Q}^{i}_{353}=\left(\frac{Q^{i}_m}{I_m^i}\right)I^i_{353}, \, \, \, \, \, \, \, \, \, \, \, \, \widetilde{U}^{i}_{353}=\left(\frac{U^{i}_m}{I_m^i}\right)I^i_{353},$$ (13)The best-fitting pitch angle q can be found minimizing some function of the difference between the simple polarization model obtained from equation (13) and the observations. We minimize the L1 norm of the difference in Q and U, summed for all the pixels at high galactic latitude:   $$G(q)= \sum _{p} \left( \left|\widetilde{Q}^{{\rm model}}_{353}(q)- Q^{{\rm obs}}_{353}\right| + \left|\widetilde{U}^{{\rm model}}_{353}(q) -U^{{\rm obs}}_{353}\right| \right),$$ (14)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 equation (13):   $$\widetilde{Q}^{{\rm model}}_{353} = \sum _{i=1}^{N} \widetilde{Q}^{i}_{353}$$ (15)and similarly for U. We find that a pitch angle of −33° provides the best fit of the GNILC maps by the BSS model at galactic latitude |b| ≥ 15°, 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°) weighted with some pixel-dependent weights Fi:   $$Q^{i}_{353}=\left(\frac{Q^{i}_m}{I_m^i}\right)I^i_{353}+F_i \left[Q^{{\rm obs}}_{353}-\sum _{j=1}^{N}\left(\frac{Q_{m}^j}{I_m^j}\right)I^j_{353}\right],$$ (16)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 Fi must be such that the polarization fraction after the redistribution of residuals does not exceed some maximum value pmax, which is a free parameter of our model, and which we pick to be 25 per cent. We fix the value of Fi as Fi = Pi/∑jPj, i.e. proportionally to the polarized dust emission fraction in each layer, unless the resulting polarization fraction exceeds pmax. When this happens, we redistribute the polarization excess in neighbouring layers. The first term in the sum on the right hand side of equation (16) is the predicted polarization of layer i, based on a polarization 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 polarized layers $$Q_m^i$$ and $$U_m^i$$ given by the large-scale model of the magnetic field while Fig. 5 shows the polarized layers after redistributing the residuals all over the former layers. After adding the small-scale 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 polarization 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 polarized dust emission on intermediate scales. In addition, the amplitude of the modelled polarized 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. Figure 4. View largeDownload slide U and Q dust emission layers obtained using a model of dust fraction in each layer based on a simple model of dust density distribution in the Galaxy, and a large-scale bi-symmetric spiral model of GMF to infer thermal dust polarization emission from dust intensity maps at 353 GHz. Figure 4. View largeDownload slide U and Q dust emission layers obtained using a model of dust fraction in each layer based on a simple model of dust density distribution in the Galaxy, and a large-scale bi-symmetric spiral model of GMF to infer thermal dust polarization emission from dust intensity maps at 353 GHz. Figure 5. View largeDownload slide U and Q dust emission layers after renormalization of the sum to match the observed sky (Planck HFI GNILC dust polarization maps at 353 GHz). Figure 5. View largeDownload slide U and Q dust emission layers after renormalization of the sum to match the observed sky (Planck HFI GNILC dust polarization maps at 353 GHz). Figure 6. View largeDownload slide U and Q layers after matching with the observed sky (as in Fig. 5), after adding random small-scale fluctuations at a level matching an extrapolation of the temperature and polarization angular power spectra and cross-spectra. Figure 6. View largeDownload slide U and Q layers after matching with the observed sky (as in Fig. 5), after adding random small-scale fluctuations at a level matching an extrapolation of the temperature and polarization angular power spectra and cross-spectra. 5 SMALL SCALES The polarization maps we have generated are normalized to match the observed dust polarization in the GNILC 353 GHz maps obtained as described in Section 3. However, the polarization GNILC maps are produced at 1° angular resolution. In the galactic plane, where the polarized signal is strong, this is the actual resolution of the GNILC map. At high galactic latitude however, the amount of polarized 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 polarized emission so far: everywhere on the sky on scales smaller than 1° (because the GNILC maps are produced at 1° 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 polarized 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, etc.). This provides us with only partial information, making it extremely difficult to generate an accurate 3D 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 μ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 equations (11) and (12) to get a first guess of the layer-dependent dust polarization, but we do not attempt to directly model the 3D 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 polarization 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 polarized dust on scales where those observations are reliable, and extend the power spectra of our maps at high l in polarization, independently for each layer, to empirically model the effect of a small-scale turbulent component of the GMF, on scales missing or noisy in the GNILC 353 GHz map. To do so, we add small-scale fluctuations independently in each layer of our model, both for intensity and for polarization. In the case of intensity, we simply fit the power spectrum of the original map in the multipole interval 30 ≤ l ≤ 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 lognormal 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 XXX (2016a), we assume a power-law dependence for EE and BB power spectra at high l, of the form Cl = A(l/lfit)α with α = −2.42 for both E and B. We use a Gaussian distribution, instead of lognormal, for polarization fields. For each layer, we fix the amplitude A and lfit to match the power spectrum of the large-scale map for that layer in the range 30 ≤ l ≤ 100. The amplitude of the small-scale fluctuations is scaled by the polarized intensity map in each layer. The randomly generated T and E harmonic coefficients are drawn with 30 per cent 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ℓ for temperature and hℓ for polarization, i.e.   $$a_{\ell m}^{T, \rm {obs}} = b_\ell \, a_{\ell m}^{T, \rm {sky}}; \quad a_{\ell m}^{E, \rm {obs}} = h_\ell \, a_{\ell m}^{E, \rm {sky}}$$ (17)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ℓm by forming   $$a_{\ell m}^{T, \rm {sim}} = a_{\ell m}^{T, \rm {obs}} + \sqrt{1-b_\ell ^2} \, a_{\ell m}^{T, \rm {rnd}}$$ (18)and similarly   $$a_{\ell m}^{E, \rm {sim}} = a_{\ell m}^{E, \rm {obs}} + \sqrt{1-h_\ell ^2} \, a_{\ell m}^{E, \rm {rnd}},$$ (19)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 arcmin for all layers in intensity, 2.5° for polarization layers 1, 2, and 3 (emission mostly at high galactic latitude), and of 2° for polarization layers 4, 5, and 6 (emission mostly near to the Galactic plane). These simulated sets of aℓm have correct $$C_\ell ^{TT}$$ and $$C_\ell ^{EE}$$, but not cross-spectrum $$C_\ell ^{TE}$$. Indeed   $$C_\ell ^{TE, {\rm sim}} = C_\ell ^{TE} \left[ b_\ell h_\ell + \sqrt{1-b_\ell ^2}\sqrt{1-h_\ell ^2} \right].$$ (20)we obtain final simulated aℓm as   $$a_{\ell m}^{\rm {final}} = \left[ C_\ell \right]^{1/2} [ C_\ell ^{\rm sim} ]^{-1/2} \, a_{\ell m}^{\rm {sim}},$$ (21)where for each ℓ, Cℓ, and $$C_\ell ^{\rm sim}$$ are 2 × 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 polarized 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 polarization fraction decreases with the polarization fraction, as seen in Fig. 9. Figure 7. View largeDownload slide First row: Full sky Q and U maps given by the BSS model. Second row: Q and U GNILC maps. Third row: Q and U total simulated maps after matching the GNILC maps on large scales and adding random small-scale fluctuations. The BSS model provides only a crude approximation of the observed dust emission. Figure 7. View largeDownload slide First row: Full sky Q and U maps given by the BSS model. Second row: Q and U GNILC maps. Third row: Q and U total simulated maps after matching the GNILC maps on large scales and adding random small-scale fluctuations. The BSS model provides only a crude approximation of the observed dust emission. Figure 8. View largeDownload slide Layers of polarized intensity ($$P = \sqrt{Q^2+U^2}$$), as modelled in our work. Figure 8. View largeDownload slide Layers of polarized intensity ($$P = \sqrt{Q^2+U^2}$$), as modelled in our work. Figure 9. View largeDownload slide Histograms of polarization fraction for each layer. We only use pixels where the polarization fraction is well defined, i.e. I(p) ≠ 0. This excludes high galactic latitude pixels for the most distant layers. Figure 9. View largeDownload slide Histograms of polarization fraction for each layer. We only use pixels where the polarization fraction is well defined, i.e. I(p) ≠ 0. This excludes high galactic latitude pixels for the most distant layers. 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-scale 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°, 50°) 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. Figure 10. View largeDownload slide T, E, B power spectra for each layer. The first three rows also display the power spectra for 75 per cent and 25 per cent of the sky. Figure 10. View largeDownload slide T, E, B power spectra for each layer. The first three rows also display the power spectra for 75 per cent and 25 per cent of the sky. Figure 11. View largeDownload slide TT, EE, BB, TE power spectra of both GNILC maps and of simulated maps including small-scale fluctuations. Figure 11. View largeDownload slide TT, EE, BB, TE power spectra of both GNILC maps and of simulated maps including small-scale fluctuations. Figure 12. View largeDownload slide Modelled E and B modes maps at 353 GHz, after adding small-scale fluctuations, adding-up six layers of emission (see the text). Figure 12. View largeDownload slide Modelled E and B modes maps at 353 GHz, after adding small-scale fluctuations, adding-up six layers of emission (see the text). Figure 13. View largeDownload slide Observed and modelled E and B modes maps at 353 GHz – detail around (l, b) = (0°, 50°). Top row: T, E, and B modes, observed with Planck after GNILC processing; Bottom row: Modelled T, E, and B modes at Nside= 512, after adding small-scale fluctuations, adding-up six layers of emission. Figure 13. View largeDownload slide Observed and modelled E and B modes maps at 353 GHz – detail around (l, b) = (0°, 50°). Top row: T, E, and B modes, observed with Planck after GNILC processing; Bottom row: Modelled T, E, and B modes at Nside= 512, after adding small-scale fluctuations, adding-up six layers of emission. Figure 14. View largeDownload slide E and B power spectra of both GNILC maps and GNILC maps + small-scale fluctuations at 143 and 217 GHz. Figure 14. View largeDownload slide E and B power spectra of both GNILC maps and GNILC maps + small-scale fluctuations at 143 and 217 GHz. 6 SCALING LAWS We now need a prescription for scaling the 353 GHz polarized 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 MBB with pixel-dependent temperature and spectral index), but now, we must define as many templates of T(p) and β(p) as there are layers in our model, i.e. six maps of T and six maps of β. A complete description of the temperature and the spectral index distribution in 3D 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 polarization maps, scaled to other frequencies than 353 GHz, should match the observed polarization 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–β scatter plot. With only 353 GHz polarization 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 XI (2014), which we assume to hold for all of the I, Q, and U Stokes parameters. We compute, for each layer i:   \begin{eqnarray} T_{{\rm avg}}^i&=&\sum _{p} w^i_T(p) T_d(p), \nonumber\\ \beta _{{\rm avg}}^i&=&\sum _{p} w^i_\beta (p) \beta _d(p), \end{eqnarray} (22)where Td(p) and βd(p) are the best-fitting 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_{{\rm avg}}^i$$ and spectral index $$\beta _{{\rm avg}}^i$$,   $$w^i_T(p) = w^i_\beta (p) = f_i(p),$$ (23)i.e. we empirically weight the maps by the pixel-dependent fraction fi(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_{{\rm avg}}^i$$ and $$\beta _{{\rm avg}}^i$$ are constant across the sky in a given layer. This however implements only a variability of the physical parameters T and β 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 β (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_{{\rm avg}}^i$$ and $$\beta _{{\rm avg}}^i$$ following the auto and cross-spectra of the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). To take into account the non-Gaussianity of the distribution of T and β, 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 β. 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. Figure 15. View largeDownload slide Top left: Power spectra used to draw random realizations of temperature and spectral index maps (note the negative sign of the Tβ cross-spectrum). Bottom left: Scatter plot of T and β for a pair of random maps (right), showing an overall anticorrelation and the same general behaviour as observed by Planck Collaboration XI (2014) on Planck observations (see their fig. 16). Right: Maps of randomly generated temperature and spectral index for the first layer, with Tavg = 19.10, σT = 2.059, βavg = 1.627, σβ = 0.209. Figure 15. View largeDownload slide Top left: Power spectra used to draw random realizations of temperature and spectral index maps (note the negative sign of the Tβ cross-spectrum). Bottom left: Scatter plot of T and β for a pair of random maps (right), showing an overall anticorrelation and the same general behaviour as observed by Planck Collaboration XI (2014) on Planck observations (see their fig. 16). Right: Maps of randomly generated temperature and spectral index for the first layer, with Tavg = 19.10, σT = 2.059, βavg = 1.627, σβ = 0.209. 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 ν = 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 β 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 β 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 β in each layer for one single realization of the simulation, compared with those from the Planck MBB fit and the fit performed on this realization. The average values from several simulations are in good agreement with those of the Planck MMB fit. Table 1. Averages and standard deviation values of temperature and spectral index in each layer, for a simulation with 6.87 arcmin 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  Tavg  19.10  18.96  18.98  19.35  19.23  20.05  σT  2.059  2.100  2.022  2.076  2.117  2.069  βavg  1.627  1.628  1.598  1.538  1.513  1.689  σβ  0.209  0.210  0.207  0.208  0.202  0.204    $$T_{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{T}^{{\rm {MMB}}}$$  $$\beta _{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{\beta }^{{\rm {MMB}}}$$      Planck fit  19.396  1.247  1.598  0.126      Simul. fit  19.389  1.253  1.598  0.135      Layer  1  2  3  4  5  6  Tavg  19.10  18.96  18.98  19.35  19.23  20.05  σT  2.059  2.100  2.022  2.076  2.117  2.069  βavg  1.627  1.628  1.598  1.538  1.513  1.689  σβ  0.209  0.210  0.207  0.208  0.202  0.204    $$T_{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{T}^{{\rm {MMB}}}$$  $$\beta _{{\rm avg}}^{{\rm {MMB}}}$$  $$\sigma _{\beta }^{{\rm {MMB}}}$$      Planck fit  19.396  1.247  1.598  0.126      Simul. fit  19.389  1.253  1.598  0.135      View Large 7 VALIDATION AND PREDICTIONS We use our model to generate maps of polarized dust emission at 143, 217 GHz, and compare them to Planck observations in polarization and intensity (Fig. 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 polarization 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 16. View largeDownload slide GNILC maps both in intensity and polarization are shown in the first and the fourth row (subindex G), while maps obtained using our 3D model are shown in the second and the fifth row (subindex m). The differences between them are also shown in the third and sixth row. Note the different colour scales for difference maps. Figure 16. View largeDownload slide GNILC maps both in intensity and polarization are shown in the first and the fourth row (subindex G), while maps obtained using our 3D model are shown in the second and the fifth row (subindex m). The differences between them are also shown in the third and sixth row. Note the different colour scales for difference maps. Fig. 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. Fig. 18 shows cross-correlation between various modelling options, showing that those models differ only at a subdominant level. These correlations are computed for maps smoothed to 2° angular resolution over 70 per cent 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 3D model in which the temperature and spectral index are constant in each layer, using the average values from Table 1; A 2D 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 equation 1); A 3D model in which each layer has a different pixel-dependent map of T and β (the main model developed in this paper). Figure 17. View largeDownload slide Cross-correlation between simulations and observations for T, E, and B power spectra at 143 and 217 GHz 70 per cent of sky. We show in blue and green the correlation in intensity and polarization between maps generated with our model and the observations. While blue curves are computed using one single value of temperature and spectral index per layer, green curves consider one template per layer, with fluctuations of the temperature and the spectral index (model b). Red curves show the correlation between the observed polarized sky maps and maps obtained from a 2D model, i.e. one single template for temperature and spectral index from the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). Figure 17. View largeDownload slide Cross-correlation between simulations and observations for T, E, and B power spectra at 143 and 217 GHz 70 per cent of sky. We show in blue and green the correlation in intensity and polarization between maps generated with our model and the observations. While blue curves are computed using one single value of temperature and spectral index per layer, green curves consider one template per layer, with fluctuations of the temperature and the spectral index (model b). Red curves show the correlation between the observed polarized sky maps and maps obtained from a 2D model, i.e. one single template for temperature and spectral index from the MBB fit obtained on the observed Planck dust maps (Planck Collaboration XLVIII 2016e). Figure 18. View largeDownload slide Cross-correlations of T, E, and B between various modelling options. Differences between these models in polarization are at the level of a few per cent at most for ℓ ≤ 100. Figure 18. View largeDownload slide Cross-correlations of T, E, and B between various modelling options. Differences between these models in polarization are at the level of a few per cent at most for ℓ ≤ 100. We see an excellent correlation overall in all cases, of more than 96 per cent at 143 GHz, and more than 99 per cent at 217 GHz for polarization, 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 polarization 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 ℓ. 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 polarization, 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 β 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 2D model with one single MBB per pixel. The same effect is observed both in intensity and in polarization. Figure 19. View largeDownload slide Left: Average total sky emission in intensity for our 3D multi-MBB model as compared to a ‘2D’ model with one single MBB per pixel. Both have the same intensity at 353 GHz by construction. The 3D model has flatter emission law at low frequency, an effect that originates from the increasing importance at low frequency of components with flatter spectral index that may be subdominant at higher frequency where the emission is dominated by hotter components. Right: Ratio of the average emission law of the 3D model and the 2D model, for both intensity and polarized intensity. Figure 19. View largeDownload slide Left: Average total sky emission in intensity for our 3D multi-MBB model as compared to a ‘2D’ model with one single MBB per pixel. Both have the same intensity at 353 GHz by construction. The 3D model has flatter emission law at low frequency, an effect that originates from the increasing importance at low frequency of components with flatter spectral index that may be subdominant at higher frequency where the emission is dominated by hotter components. Right: Ratio of the average emission law of the 3D model and the 2D model, for both intensity and polarized intensity. Finally, we can compute the level of decorrelation between polarization 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, de Oliveira-Costa & Hamilton 2003; Eriksen et al. 2004; Delabrouille et al. 2009), CCA (Bonaldi et al. 2006), or SMICA (Delabrouille, Cardoso & Patanchon 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 polarization maps (both E and B) at 143 or 217, and 353 GHz (see Fig. 20) for our 3D model. The correlations obtained in both cases are ranging from 97 per cent on small scale to close to 100 per cent on large scales, which is larger than what is observed on real Planck maps (Planck Collaboration L 2017). This shows that even if our multilayer model adds a level of complexity to dust emission modelling, it cannot produce a decorrelation between frequencies as strong as originally claimed in the first analysis of Planck polarization maps (Planck Collaboration L 2017). Our model, however, is compatible with the lack of evidence for such decorrelation between 217 and 353 GHz at the 0.4 per cent level for 55 ≤ ℓ ≤ 90 claimed in Sheehy & Slosar (2017), and predicts increased decorrelation (of the order of 1–2 per cent) between 143 and 353 GHz over the same range of ℓ. More multifrequency observations of polarized dust emission are necessary to better model dust polarization 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. Figure 20. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over 70 per cent of sky. Figure 20. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over 70 per cent of sky. Figure 21. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over different sky fractions. Figure 21. View largeDownload slide Correlation between maps at different frequencies obtained with our 3D model, computed over different sky fractions. 8 CONCLUSION We have developed a 3D model of polarized Galactic dust emission that is consistent with the large-scale Planck HFI polarization 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 polarization 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 arcmin . For polarization, small-scale fluctuations of emission originating from the turbulence of the GMF are randomly generated on scales smaller than 2° or 2.5°, 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 per cent 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 LOS. We hence assign to each layer of emission a different, pixel-dependent, scaling law in the form of a MBB emission characterized, 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 polarization forbids such an ambition. None the less, the model provides a means to simulate a dust component that features some of the plausible complexity of the polarized dust component, while being compatible with the observed large-scale polarized emission at 353 GHz and with most of the observed statistical properties of dust (temperature and polarization power spectra, amplitude and correlation of temperature and spectral index of the best-fitting MBB emission). However, this model fails to predict the strong decorrelation of dust polarization between frequency channels on small angular scales seen in Planck Collaboration L (2017), a limitation that must be addressed in the future if that decorrelation is confirmed. In the meantime, we expect these simulated maps to be useful to investigate the component separation problem for future CMB polarization 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 pixelization scheme (Górski et al. 2005), available from the healpix webpage,5 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. Footnotes 1 Nor does the specific analytic form of the depolarization function. 2 After 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. 3 However, 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 polarization fraction of dust emission in each layer. 4 We postpone to Section 6 the discussion of temperature and spectral index maps. In equation (8), we use average values given in Table 1. 5 http://healpix.sourceforge.net REFERENCES Abazajian K. N. et al.  , 2016, preprint (arXiv:1610.02743) André P. et al.  , 2014, J. Cosmol. Astropart. Phys. , 2, 006 CrossRef Search ADS   Basak S., Delabrouille J., 2012, MNRAS , 419, 1163 CrossRef Search ADS   Basak S., Delabrouille J., 2013, MNRAS , 435, 18 CrossRef Search ADS   Bennett C. L. et al.  , 2013, ApJS , 208, 20 CrossRef Search ADS   Benoit A. et al.  , 2002, Astropart. Phys. , 17, 101 CrossRef Search ADS   Betoule M., Pierpaoli E., Delabrouille J., Le Jeune M., Cardoso J.-F., 2009, A&A , 503, 691 CrossRef Search ADS   BICEP2 Collaboration, 2014, Phys. Rev. Lett. , 112, 241101 CrossRef Search ADS PubMed  BICEP2/Keck & Planck Collaborations, 2015, Phys. Rev. Lett. , 114, 101301 CrossRef Search ADS PubMed  Bonaldi A., Ricciardi S., 2011, MNRAS , 414, 615 CrossRef Search ADS   Bonaldi A., Bedini L., Salerno E., Baccigalupi C., de Zotti G., 2006, MNRAS , 373, 271 CrossRef Search ADS   Bonaldi A., Ricciardi S., Brown M. L., 2014, MNRAS , 444, 1034 CrossRef Search ADS   Cardoso J.-F., Le Jeune M., Delabrouille J., Betoule M., Patanchon G., 2008, IEEE J. Sel. Top. Signal Process. , 2, 735 CrossRef Search ADS   Challinor A. et al.  , 2017, preprint (arXiv:1707.02259) CORE Collaboration, 2016, preprint (arXiv:1612.08270) Das S. et al.  , 2014, J. Cosmol. Astropart. Phys. , 4, 014 CrossRef Search ADS   de Bernardis P. et al.  , 2000, Nature , 404, 955 CrossRef Search ADS PubMed  Delabrouille J., Cardoso J.-F., 2009, in Martínez V. J., Saar E., Martínez-González E., Pons-Bordería M.-J., eds, Lecture Notes in Physics, Vol. 665, Data Analysis in Cosmology . Springer-Verlag, Berlin, p. 159 Delabrouille J., Cardoso J.-F., Patanchon G., 2003, MNRAS , 346, 1089 CrossRef Search ADS   Delabrouille J., Cardoso J.-F., Le Jeune M., Betoule M., Fay G., Guilloux F., 2009, A&A , 493, 835 CrossRef Search ADS   Delabrouille J. et al.  , 2013, A&A , 553, A96 CrossRef Search ADS   Delabrouille J. et al.  , 2017, preprint (arXiv:1706.04516) Dunkley J. et al.  , 2009, in Dodelson S. et al.  , eds, AIP Conf. Ser. Vol. 1141, CMB Polarization Workshop: Theory and Foregrounds: CMBPol Mission Concept Study . Am. Inst. Phys., New York, p. 222 Efstathiou G., Gratton S., Paci F., 2009, MNRAS , 397, 1355 CrossRef Search ADS   Eriksen H. K., Banday A. J., Górski K. M., Lilje P. B., 2004, ApJ , 612, 633 CrossRef Search ADS   Errard J., Stompor R., 2012, Phys. Rev. D , 85, 083006 CrossRef Search ADS   Fauvet L. et al.  , 2011, A&A , 526, A145 CrossRef Search ADS   Faÿ G., Guilloux F., Betoule M., Cardoso J.-F., Delabrouille J., Le Jeune M., 2008, Phys. Rev. D , 78, 083013 CrossRef Search ADS   Ghosh T. et al.  , 2017, A&A , 601, A71 CrossRef Search ADS   Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ , 622, 759 CrossRef Search ADS   Green G. M. et al.  , 2015, ApJ , 810, 25 CrossRef Search ADS   Han J. L., 2017, ARA&A , 55, 111 CrossRef Search ADS   Han J. L., Qiao G. J., 1993, Acta Astrophys. Sin. , 13, 385 Hanany S. et al.  , 2000, ApJ , 545, L5 CrossRef Search ADS   Harari D., Mollerach S., Roulet E., 1999, J. High Energy Phys. , 8, 022 CrossRef Search ADS   Haverkorn M., 2015, in Lazarian A., de Gouveia Dal Pino E. M., Melioli C., eds, Astrophysics and Space Science Library, Vol. 407, Magnetic Fields in Diffuse Media . Springer-Verlag Berlin, p. 483 Jaffe T. R., Leahy J. P., Banday A. J., Leach S. M., Lowe S. R., Wilkinson A., 2010, MNRAS , 401, 1013 CrossRef Search ADS   Jones W. C. et al.  , 2006, ApJ , 647, 823 CrossRef Search ADS   Kamionkowski M., Kovetz E. D., 2016, ARA&A , 54, 227 CrossRef Search ADS   Kogut A. et al.  , 2011, J. Cosmol. Astropart. Phys. , 7, 025 CrossRef Search ADS   Kovac J. M., Leitch E. M., Pryke C., Carlstrom J. E., Halverson N. W., Holzapfel W. L., 2002, Nature , 420, 772 CrossRef Search ADS PubMed  Lallement R., Snowden S., Kuntz K. D., Dame T. M., Koutroumpa D., Grenier I., Casandjian J. M., 2016, A&A , 595, A131 CrossRef Search ADS   Leach S. M. et al.  , 2008, A&A , 491, 597 CrossRef Search ADS   Lewis A., Challinor A., 2006, Phys. Rep. , 429, 1 CrossRef Search ADS   Marinucci D. et al.  , 2008, MNRAS , 383, 539 CrossRef Search ADS   Matsumura T. et al.  , 2014, J. Low Temp. Phys. , 176, 733 CrossRef Search ADS   Narcowich F., Petrushev P., Ward J., 2006, SIAM J. Math. Anal. , 38, 574 CrossRef Search ADS   O'Dea D. T., Clark C. N., Contaldi C. R., MacTavish C. J., 2012, MNRAS , 419, 1795 CrossRef Search ADS   Penzias A. A., Wilson R. W., 1965, ApJ , 142, 419 CrossRef Search ADS   Planck Collaboration XI, 2014, A&A , 571, A11 CrossRef Search ADS   Planck Collaboration XXX, 2016a, A&A , 586, A133 CrossRef Search ADS   Planck Collaboration I, 2016b, A&A , 594, A1 CrossRef Search ADS   Planck Collaboration XIII, 2016c, A&A , 594, A13 CrossRef Search ADS   Planck Collaboration XLIV, 2016d, A&A , 596, A105 CrossRef Search ADS   Planck Collaboration XLVIII, 2016e, A&A , 596, A109 CrossRef Search ADS   Planck Collaboration L, 2017, A&A , 599, A51 CrossRef Search ADS   Reichardt C. L. et al.  , 2009, ApJ , 694, 1200 CrossRef Search ADS   Remazeilles M., Delabrouille J., Cardoso J.-F., 2011, MNRAS , 418, 467 CrossRef Search ADS   Remazeilles M., Dickinson C., Eriksen H. K. K., Wehus I. K., 2016, MNRAS , 458, 2032 CrossRef Search ADS   Remazeilles M. et al.  , 2017, preprint (arXiv:1704.04501) Rezaei Kh. S., Bailer-Jones C. A. L., Hanson R. J., Fouesneau M., 2017, A&A , 598, A125 CrossRef Search ADS   Sheehy C., Slosar A., 2017, preprint (arXiv:1709.09729) Smoot G. F. et al.  , 1992, ApJ , 396, L1 CrossRef Search ADS   Sofue Y., Fujimoto M., 1983, ApJ , 265, 722 CrossRef Search ADS   Stanev T., 1997, ApJ , 479, 290 CrossRef Search ADS   Stompor R., Errard J., Poletti D., 2016, Phys. Rev. D , 94, 083526 CrossRef Search ADS   Story K. T. et al.  , 2013, ApJ , 779, 86 CrossRef Search ADS   Tassis K., Pavlidou V., 2015, MNRAS , 451, L90 CrossRef Search ADS   Tauber J. A. et al.  , 2010, A&A , 520, A1 CrossRef Search ADS   Tegmark M., de Oliveira-Costa A., Hamilton A. J., 2003, Phys. Rev. D , 68, 123523 CrossRef Search ADS   The COrE Collaboration, 2011, preprint (arXiv:1102.2181) Tinyakov P. G., Tkachev I. I., 2002, Astropart. Phys. , 18, 165 CrossRef Search ADS   Tucci M., Martínez-González E., Vielva P., Delabrouille J., 2005, MNRAS , 360, 935 CrossRef Search ADS   Vansyngel F. et al.  , 2017, A&A , 603, A62 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations