Issue 
A&A
Volume 603, July 2017



Article Number  A62  
Number of page(s)  14  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201629992  
Published online  10 July 2017 
Statistical simulations of the dust foreground to cosmic microwave background polarization
^{1} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{2} California Institute of Technology, Pasadena, California CA 91125, USA
^{3} Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98bis bd Arago, 75014 Paris, France
^{4} Sorbonne Universités, Institut Lagrange de Paris (ILP), 98bis Boulevard Arago, 75014 Paris, France
^{5} Laboratoire AIM, IRFU/Service d’Astrophysique, CEA/DSM, CNRS, Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{6} LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, École normale supérieure, 75005 Paris, France
^{7} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{8} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
Received: 1 November 2016
Accepted: 24 March 2017
The characterization of the dust polarization foreground to the cosmic microwave background (CMB) is a necessary step toward the detection of the Bmode signal associated with primordial gravitational waves. We present a method to simulate maps of polarized dust emission on the sphere that is similar to the approach used for CMB anisotropies. This method builds on the understanding of Galactic polarization stemming from the analysis of Planck data. It relates the dust polarization sky to the structure of the Galactic magnetic field and its coupling with interstellar matter and turbulence. The Galactic magnetic field is modeled as a superposition of a mean uniform field and a Gaussian random (turbulent) component with a powerlaw power spectrum of exponent α_{M}. The integration along the line of sight carried out to compute Stokes maps is approximated by a sum over a small number of emitting layers with different realizations of the random component of the magnetic field. The model parameters are constrained to fit the power spectra of dust polarization EE, BB, and TE measured using Planck data. We find that the slopes of the E and B power spectra of dust polarization are matched for α_{M} = −2.5, an exponent close to that measured for total dust intensity but larger than the Kolmogorov exponent − 11/3. The model allows us to compute multiple realizations of the Stokes Q and U maps for different realizations of the random component of the magnetic field, and to quantify the variance of dust polarization spectra for any given sky area outside of the Galactic plane. The simulations reproduce the scaling relation between the dust polarization power and the mean total dust intensity including the observed dispersion around the mean relation. We also propose a method to carry out multifrequency simulations, including the decorrelation measured recently by Planck, using a given covariance matrix of the polarization maps. These simulations are well suited to optimize component separation methods and to quantify the confidence with which the dust and CMB Bmodes can be separated in present and future experiments. We also provide an astrophysical perspective on our phenomenological modeling of the dust polarization spectra.
Key words: polarization / ISM: general / cosmic background radiation / submillimeter: ISM
© ESO, 2017
1. Introduction
An era of exponential expansion of the Universe, dubbed cosmic inflation, has been proposed to explain why the Universe is almost exactly Euclidean and nearly isotropic (Guth 1981; Linde 1982). One generic prediction of this theoretical paradigm is the existence of a background of gravitational waves, which produces a distinct, curllike, signature in the polarization of the cosmic microwave background (CMB), referred to as primordial Bmode polarization (Starobinskiǐ 1979). The detection of this signal would have a deep impact on cosmology and fundamental physics, motivating a number of experiments designed to measure the sky polarization at microwave frequencies. Current projects have achieved the sensitivity required to detect the CMB Bmode signal predicted by the simplest models of inflation (Abazajian et al. 2015; Kamionkowski & Kovetz 2016). Yet, any detection has relied on the proper removal of much brighter Galactic foregrounds.
Thermal emission from aspherical dust grains aligned with respect to the Galactic magnetic field (GMF) is the dominant polarized foreground for frequencies higher than about 70GHz (Dunkley et al. 2009; Planck Collaboration X 2016). From the analysis of the Planck^{1}353GHz polarization maps, we know that the primordial Bmode polarization of the CMB cannot be measured without subtracting the foreground emission, even in the faintest dustemitting regions at high Galactic latitude (Planck Collaboration Int. XXX 2016, hereafter PXXX). The observed correlation between the Bmode signal detected by BICEP2/Keck Array, on the one hand, and the Planck dust maps, on the other hand, has confirmed this conclusion (BICEP2/Keck Array and Planck Collaborations 2015).
To distinguish cosmological and Galactic foreground polarization signals, CMB experiments must rely on multifrequency observations. Component separation is a main challenge because the spatial structure of dust polarization is observed to vary with frequency (Planck Collaboration Int. L 2017). This introduces two questions that motivate our work. What design of CMB experiments and combination of groundbased, balloonborne, and space observations is best to achieve an optimal separation? How can confidence in the subtraction of foregrounds be quantified? To provide quantitative answers, we must be able to simulate observations of the polarized sky combining Galactic and CMB polarization. This paper presents a statistical model with a few parameters to simulate maps of dust polarization in a way similar to what is available for CMB anisotropies (Seljak & Zaldarriaga 1996).
The analysis of the Wilkinson Microwave Anisotropy Probe (WMAP) data and the preparation of the Planck project motivated a series of models of the polarized synchrotron and thermal dust emission at microwave frequencies (Page et al. 2007; MivilleDeschênes et al. 2008; Fauvet et al. 2011; O’Dea et al. 2012; Delabrouille et al. 2013). These early studies followed two distinct approaches. The first is to produce a sky that is as close as possible to the observed sky combining data templates and a spectral model. Prior to Planck, for dust polarization this was performed using stellar polarization data by Page et al. (2007), and WMAP observations of synchrotron polarization by Delabrouille et al. (2013). More recently, the simulations presented in Planck Collaboration XII (2016) use the Planck353GHz data to model dust polarization. This first approach is limited by the signaltonoise ratio of available data, which for dust polarization is low at high Galactic latitude even after smoothing to one degree angular resolution. The second approach is to simulate the polarization sky from a 3D model of the GMF and of the density structure of the interstellar medium (ISM), both its regular and turbulent components, as carried out by MivilleDeschênes et al. (2008), Fauvet et al. (2011) and O’Dea et al. (2012). This method connects the modeling of the microwave polarized sky to broader efforts to model the GMF (Waelkens et al. 2009; Jansson & Farrar 2012; Planck Collaboration Int. XLII 2016).
Planck polarization maps have been used to characterize the structure (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XX 2015) and the spectral energy distribution (SED) of polarized thermal emission from Galactic dust (Planck Collaboration Int. XXI 2015; Planck Collaboration Int. XXII 2015). Several studies have established the connection between the structure of the magnetic field and matter (Clark et al. 2014; Planck Collaboration Int. XX 2015; Planck Collaboration Int. XXXII 2016; Martin et al. 2015; Kalberla et al. 2016). The power spectra analysis presented in PXXX decomposes dust polarization into E (gradientlike) and B (curllike) modes (Zaldarriaga 2001; Caldwell et al. 2017). This analysis led to two unexpected results: a positive TE correlation and a ratio of about 2 between the E and B dust powers over the ℓ range 40 to 600. Clark et al. (2015) and Planck Collaboration Int. XXXVIII (2016) have showed that the observed TE correlation and asymmetry between E and Bmode power amplitudes for dust polarization could be both accounted for by the preferred alignment between the filamentary structure of the total intensity map and the orientation of the magnetic field inferred from the polarization angle.
The work presented here makes use of the model framework introduced in (Planck Collaboration Int. XLIV 2016; hereafter PXLIV). By analyzing Planck dust polarization maps toward the southern Galactic cap, the part of the sky used for CMB observations from Antartica and Atacama, PXLIV related the largescale patterns of the maps to the mean orientation of the magnetic field, and the scatter of the dust polarization angle and fraction (ψ and p) to the amplitude of its turbulent component. In this paper, we extend their work to produce Stokes maps that fit dust polarization power spectra including the TE correlation and the TT/EE and EE/BB power ratios at high and intermediate Galactic latitudes. In a companion paper Ghosh et al. (2017), the dust polarization of the southern sky region with the lowest dust column density is modeled using Hi observations and astrophysical insight to constrain their model parameters. In essence, our approach is more mathematical but it allows us to model dust polarization over a larger fraction of the sky. The two approaches are complementary and compared in this paper. We also present a mathematical process to introduce spatial decorrelation across microwave frequencies via the auto and cross spectra of dust polarization. By doing this, we obtain a model to compute independent realizations of dust polarization sky maps at one or multiple frequencies with a few parameters adjusted to fit the statistical properties inferred from the analysis of the Planck data away from the Galactic plane.
The paper is organized as follows. Sections 2 and 3 present the framework we use to model dust polarization in general terms. Our method is illustrated by producing simulated maps at 353GHz presented in Sect. 4. We show that these maps successfully match the statistical properties of dust polarization derived from the analysis of Planck data (Sect. 5). One method to compute dust polarization maps at multiple frequencies is presented in Sect. 6. We discuss the astrophysical implications of our work in Sect. 7. The main results of the paper are summarized in Sect. 8. Appendix A details how the simulated maps used in this study are computed. Appendix B shows how to compute the cross correlation between two frequency maps, when spectral differences about a mean SED may be parametrized with a spatially varying spectral index.
2. Astrophysical framework
To model dust polarization we used the framework introduced by PXLIV, which we briefly describe here. We refer to PXLIV for a detailed presentation and discussion of the astrophysical motivation and the simplifying assumptions of our modeling approach.
The polarization of thermal dust emission results from the alignment of aspherical grains with respect to the GMF (Stein 1966; Lee & Draine 1985; Planck Collaboration Int. XXI 2015). Within the hypothesis that grain polarization properties, including alignment, are homogeneous, the structure of the dust polarization sky reflects the structure of the magnetic field combined with that of matter. We assume that this hypothesis applies to the diffuse ISM where radiative torques provide a viable mechanism to align grains efficiently (Dolginov & Mitrofanov 1976; Andersson et al. 2015; Hoang & Lazarian 2016).
To compute the Stokes parameters I, Q, and U describing the linearly polarized thermal dust emission, we start from the integral equations in Sect. 3.2 and Appendix B of Planck Collaboration Int. XX (2015) for optically thin emission at frequency ν, i.e., (1)where S(ν) is the source function, τ_{ν} the optical depth, p_{0} a parameter related to dust polarization properties (the grain cross sections and the degree of alignment with the magnetic field), γ the angle that the local magnetic field makes with the plane of the sky, and φ the local polarization angle (see Fig. 14 in Planck Collaboration Int. XX 2015).
As in PXLIV, the integration along the line of sight is approximated by a sum over a finite number N of layers. This sum is written as (2)where S_{i}(ν) is the integral of the source function over layer i, and γ_{i} and φ_{i} define the magnetic field orientation within each layer. As discussed in PXLIV, the layers are a phenomenological means to model the density structure of the interstellar matter and the correlation length of the GMF. This approach accounts for both signatures of the turbulent magnetic field component in Galactic polarization maps: the depolarization resulting from the integration along the line of sight of emission with varying polarization orientations, and the scale invariant structure of the polarization maps across the sky reflecting the power spectrum of the turbulent component of the magnetic field (Cho & Lazarian 2002; Houde et al. 2009). It overcomes the difficulty of generating realizations of the turbulent component of the magnetic field in three dimensions over the celestial sphere. Ghosh et al. (2017) uses Hi data to associate the layers with different phases of the ISM, each of which provide a different intensity map. On the contrary, in the simulations presented in this paper, like in PXLIV, the term S_{i}(ν) in Eqs. (2) is a sky map assumed to be the same in each layer, i.e., it is independent of the index i. Thus we do not address the question of the physical meaning of the layers.
Through the angles γ_{i} and φ_{i}, the model relates the dust polarization to the structure of the GMF. The magnetic field B is expressed as the sum of its mean (ordered), B_{0}, and turbulent (random), B_{t}, components, (3)where and are unit vectors in the directions of B_{0} and B_{t}, and f_{M} a model parameter that sets the relative strength of the random component of the field. To simulate dust as a foreground to the CMB we need a description of the GMF within the solar neighborhood. We follow PXLIV in assuming that B_{0} has a fixed orientation in all layers. We ignore the structure of the GMF on galaxywide scales because the dust emission arises mainly from a thin disk with a relatively small scale height and we are interested in modeling dust polarization away from the Galactic plane. This scale height is not measured directly in the solar neighborhood but modeling of the dust emission from the Milky Way indicates that it is ~200pc at the solar distance from the Galactic center (Drimmel & Spergel 2001).
Each component of the vector field in 3D, in each layer, is obtained from independent Gaussian realizations of a powerlaw power spectrum, which is written as (4)Our modeling of is continuous over the celestial sphere and uncorrelated between layers. The coherence of the GMF orientation along the line of sight comes from the mean field and is controlled by the parameter f_{M}.
The model has six parameters: the Galactic longitude and latitude l_{0} and b_{0} defining the orientation of , the number of layers N, the spectral exponent α_{M}, and the effective polarization fraction of the dust emission p_{0}. The PXLIV authors used the same model to analyze the dust polarization measured by Planck at 353GHz over the southern Galactic cap (Galactic latitude b< −60°). They determined l_{0} = 70° ± 5° and b_{0} = 24° ± 5° by fitting the largescale pattern observed in the Stokes Q and U maps, and f_{M} = 0.9 ± 0.1, N = 7 ± 2 and p_{0} = 26 ± 3% by fitting the distribution function (onepoint statistics) of p^{2}, the square of the dust polarization fraction p, and of the polarization angle ψ, computed after removal of the regular pattern from the ordered component of the GMF.
Hereafter we label the Stokes maps computed from Eqs. (2) as I_{a},Q_{a}, and U_{a}. At this stage a, the power spectra of the model maps have equal EE and BB power, and no TE correlation at ℓ ≳ 30. This follows from the fact that our modeling does not include the alignment observed between the filamentary structure of the diffuse ISM and the GMF orientation. Some TE correlation is present at low ℓ because the mean GMF orientation is close to being within the Galactic disk, and we take into account the latitude dependence of the total dust intensity. In the next section, we explain how we modify the spherical harmonic decomposition of the stage a maps to introduce the TE correlation and the EB asymmetry, matching the Planck dust polarization power spectra in PXXX.
3. Introducing TE correlation and EB asymmetry
Our aim is to simulate maps that match given observables based on dust angular power spectra, namely the TE correlation, TT/EE and EE/BB ratios, and the BB spectrum without altering the statistics of p and ψ of the Stokes maps from stage a. We describe a generic process to construct such a set of Stokes maps (I_{b},Q_{b},U_{b}), later referred to as stage b maps. The process can be applied on a full, or a masked, sky.
We start with the Stokes maps (I_{a},Q_{a},U_{a}) obtained as described in Sect. 2. We compute the spherical harmonic coefficients of the stage b maps from those of the stage a maps as follows: (5)where and denote the coefficients of the X = T,E,B harmonic decomposition of stage a and b maps, respectively. The parameter ρ introduces the TE correlation and the factor f the EB asymmetry. The parameter t is a scaling factor for the intensity part and p_{0} is the polarization parameter introduced in Eqs. (1). These parameters control the amplitude of the TT, EE, BB, and TE power spectra of stage b maps. We note that Q_{a} and U_{a} scale linearly with p_{0} and thus that the two ratios and in Eqs. (5) are independent of p_{0}.
At this stage b, our modeling of the random component of the magnetic field is anisotropic, which is a fundamental characteristic of magnetohydrodynamical turbulence (Lazarian & Pogosyan 2012; Brandenburg & Lazarian 2013). The factors ρ and f introduce anisotropy in two ways. First, the T map, which is added to the polarization part through the parameter ρ, has a filamentary structure and thus is anisotropic. This amounts to adding an extra polarization layer that is perfectly aligned with the filamentary structure of the matter and is similar to what is carried out by Ghosh et al. (2017) for their cold neutral medium map. Second, the factor f breaks the symmetry between E and B, whereas the power is expected to be equally distributed between E and B modes in the case of isotropic turbulence (Caldwell et al. 2017). Through the parameter f, the random component of B is anisotropic in all layers and everywhere on the sky, unlike in Ghosh et al. (2017) where anisotropy is introduced in only one layer.
In the simplest case, ρ,f,t, and p_{0} are constants over the whole multipole range and in the most general case they are functions of ℓ and m. We find that the statistics of p and ψ found using the stage a maps are lost at stage b if f ≠ 1 or ρ ≠ 0 for very low multipoles. Thus, we look for a solution where the parameters t and p_{0} are constants but f and ρ depend on ℓ and tend toward 1 and 0 for very low ℓ values, respectively.
The power spectra of stage b maps are noted with X,Y = T,E,B and use the quantity . The t, p_{0}, ρ, and f coefficients in Eqs. (5) are chosen such that the power spectra of stage b maps match a given set of averaged ratios as follows: (6)where ℰ[·] is a given averaging process over multipoles. The absolute scaling is performed by matching the amplitude of one power spectrum. For this purpose, we use the BB spectrum because the main motivation of the simulations is to produce polarized dust skies for component separation of Bmodes. Thus, to Eqs. (6) we add the fourth constraint (7)where N_{B} is an overall factor that scales the BB power spectrum of stage a maps divided by p_{0} to the desired amplitude. The four parameters t, p_{0}, ρ, and f can be derived analytically from the four input parameters R_{TT}, R_{TE}, R_{BB}, and N_{B}. One can choose any values for R_{TT}, R_{TE}, R_{BB}, and N_{B}, as long as the normalization is positive and the ratios respect the condition forced by the positive definiteness of the power spectra covariance.
We construct the , , and according to their definitions in Eqs. (5). The final product is a triplet of Stokes maps (I_{b},Q_{b},U_{b}) that have the desired twopoint statistics.
4. Simulated maps
To illustrate our method, we apply the formalism presented in the previous sections and simulate dust polarization maps that fit the Planck power spectra. The input values for R_{TT}, R_{TE}, and R_{BB} in Eqs. (6) are derived from Planck data (Sect. 4.1). We introduce the simulated maps in Sect. 4.2. The method used to compute these maps is detailed in Appendix A.
4.1. Planck power spectra
The EE, BB, TE, and TB angular power spectra of dust polarization were measured using the Planck maps at 353GHz on the six large regions at high and intermediate Galactic latitude defined in PXXX. The effective sky fraction f_{sky}, after a 5° (FWHM) apodization, ranges from f_{sky} = 24% to f_{sky} = 72%. The regions are labeled LRxx, with xx the sky fraction in percent.
The EE and BB spectra reported in PXXX are well fitted by power laws with exponents , with no systematic dependence on the sky region. The amplitudes of the spectra at a reference multipole ℓ_{0} = 80, A^{EE,data}, were measured from powerlaw fits over the range 40 <ℓ< 600 with an index fixed to its mean value of − 2.42. These amplitudes are observed to increase with the mean total dust intensity in the mask, I_{dust}, following the law (X = E,B). We combine the amplitude A^{EE,data} and the EE to BB ratios listed in Table 1 of PXXX for their LR33 mask to compute the amplitude A^{BB,data} of the spectrum at ℓ = 80.
The values of the R_{TT} and R_{TE} ratios are not listed in Table 1 of PXXX. To determine these values, we combine the fit to the EE spectrum from PXXX, the TE spectrum plotted in Fig. B.1 of PXXX, and the TT spectrum we computed using the Planck dust map at 353 GHz obtained by Planck Collaboration Int. XLVIII (2016) after separation from the cosmic infrared background (CIB) anisotropies. The data points and error bars of the TE spectrum were provided to us by the contact author of PXXX. The spectra are binned between ℓ = 40 and ℓ = 600 with Δℓ = 20 and the binned spectra are noted (XY = TT,TE). We compute the ratios R_{XY} by comparing the measured power spectra with the powerlaw fit to the EE spectrum , minimizing the following chisquared: (8)where is the standard deviation error on output from the Xpol power spectrum estimator^{2}.
The values we use as input for the simulations are gathered in Table 1.
Input values for the simulations.
4.2. Simulated maps used in this study
Here and in Appendix A, we introduce the simulated maps and describe how we produce them.
We have analyzed the simulated maps over a larger sky area than in PXLIV. We have not, however, attempted to fit the PXLIV model of the mean field to the Planck data over a larger region. In particular, the adopted mean field direction is given by the same Galactic coordinates (l_{0},b_{0}) = (70°,24°). Although this specific choice affects the Q_{a} and U_{a} maps, it has no critical impact on the statistical results presented in the paper. Our fiducial set of values for N, f_{M}, and α_{M} is 4, 0.9, and − 2.5, respectively. To quantify the impact of these parameters on the model power spectra, we computed simulated maps for several combinations around the fiducial values within the constraints set by PXLIV. For N we considered two values 4 and 7, and for f_{M} the range 0.7 to 1.0. We explored a range of values of α_{M} from − 3.4 to − 2.2.
The method we followed to construct the stage a and b maps is described in Sects. A.1 and A.2. We produced our simulations at an angular resolution of 30′ on a HEALPix^{3} (Górski et al. 2005) grid with resolution parameter N_{side} = 256. Although the parameter p_{0} was computed at stage b, we needed an initial guess in order to compute the total intensity map of stage a maps (see Eqs. (2) for I(ν)). Based on PXLIV, we took p_{0} = 0.25. We used this value to compute stage b maps from Q_{a}/p_{0} and U_{a}/p_{0} that do not depend much on p_{0} (Sect. A.1 ).
The parameters t,p_{0},ρ, and f used to construct stage b maps were determined by the ratios R_{TT}, R_{TE}, and R_{BB} and the amplitude of the BB spectrum (Sect. A.2). We used the BB amplitude and the ratio values computed on the LR33 mask (Table 1). The corresponding values of the stage b parameters t, p_{0}, ρ and f are listed in Table 2 for our fiducial set of values for N, f_{M} and α_{M}. The value of p_{0}, 0.22 ± 0.05 agrees with that derived by PXLIV from their data fit, which we used to compute the stage a maps. Thus, it is not necessary to iterate the process. The scaling factor t of the Stokes I map is found to be unity within uncertainties.
Because the stage a maps have a high intensity contrast, the conversion from pixel space to spherical harmonic space induces leakage of power from the Galactic plane to high latitudes. In order to avoid this artifact, the brightest part of the Galactic plane must be masked before performing the transformation. The Planck collaboration provides eight Galactic masks for general purposes. They are derived from the 353 GHz intensity map by gradually thresholding the intensity after having subtracted the CMB. These masks are then apodized with a 2 degree Gaussian kernel and cover respectively 15, 33, 51, 62, 72, 81, 91, and 95% of the sky^{4}. The precise choice of the mask is not critical. We chose the mask corresponding to f_{sky} = 80%, which discards low Galactic latitude areas where our model with a uniform mean orientation of the field does not apply. The unmasked region is large enough to encompass all regions outside the Galactic plane that are relevant for CMB analyses.
As mentioned in Sect. 3, extending the EB asymmetry down to very low multipoles changes the onepoint statistics of fraction and angle of polarization. To prevent this artefact, we introduce the EB asymmetry and the TE correlation smoothly from low multipoles. In practice, the parameters ρ and f are functions of ℓ as follows: (9)Here w(ℓ) is a window function going smoothly from 0 to 1 around multipole ℓ_{c} and is defined as follows: (10)where we set ℓ_{c} = 30 and δℓ = 30. After this modification, the EB power ratio tends to 1 for ℓ<ℓ_{c} in agreement with the EE and BBPlanck353GHz power spectra presented in Fig. 20 of Planck Collaboration Int. XLVI (2016) at ℓ< 30. Figure 1 shows that the distributions (onepoint statistics) of p and ψ computed around the southern Galactic pole of the stage a and b maps are very similar.
Fig. 1 Probability distribution functions of p^{2} and ψ (top and bottom plots) for stage a and stage b maps (red and blue histograms). The maps were computed using the fiducial values of α_{M}, f_{M}, and N and the corresponding parameters t, p_{0}, ρ, and f introduced in Sect. 4.2. The distributions are computed on the southern Galactic polar cap (b ≤ −60°) as in PXLIV. The very close match between the corresponding histograms shows that the inclusion of the TE correlation and the EB asymmetry does not alter the onepoint statistics of the simulated maps. 
Fig. 2 Parameter p_{0} (top) and slopes of the EE (middle) and BB (bottom) power spectra of stage b maps vs. the slope α_{M} of the power spectrum of the turbulent component of the magnetic field (left) and vs. the relative strength of the turbulence f_{M} (right). In the left plots f_{M} = 0.8 and in the right plots α_{M} = −2.5. Red stars (blue squares) represent results for N = 4 (N = 7). The abscissae of the two sets of points are slightly shifted from their original values for a better visibility. The values observed in the data are represented by a gray shaded region for α_{EE} and α_{BB} (a dashed line for the mean, dark, and light gray for the 1 and 2σ uncertainties) and by a hatched regions for p_{0} (a dashed horizontal line for the mean, and a red 45° (resp. blue −45°) hatched region for 1σ uncertainty for N = 4 (resp. 7)). 
5. Model power spectra
In this section, we show that our simulated stage b maps reproduce the PlanckEE, BB, and TE dust spectra constraining the exponent α_{M} of the magnetic field power spectrum (Sect. 5.1), provide the statistical variance of the dust polarization power in a given ℓ bin (Sect. 5.2), and match the observed scaling between the spectra amplitude and the mean dust total intensity for both large and small sky regions (Sect. 5.3).
5.1. Matching Planck power spectra
To compare our model results directly with the analysis of the Planck data in PXXX, we compute power spectra of the simulated maps over the LR33 mask. The power spectra are computed using the PolSpice estimator (Chon et al. 2004) that corrects for multipoletomultipole coupling due to the masking. We checked that we obtain very similar results when the spectra are computed with the Xpol estimator.
For both values N = 4 and 7, we vary the parameters f_{M} and α_{M} as follows. First, we keep f_{M} fixed to 0.9 and let α_{M} vary from − 3.4 to − 2.0 in steps of 0.2 with the addition of − 2.5, then we keep fixed α_{M} to − 2.5 and let f_{M} vary from 0.7 to 1 in steps of 0.1. For each set of parameters, we compute a sample of 1000 realizations with the procedure described in Sect. 4.2 and Appendix A. The power spectra of stage b maps are binned from ℓ = 60 to 200 with a bin width of Δℓ = 20. We fit the model A^{XX}^{(}ℓ/ℓ_{0}^{)}^{αXX + 2} (X = E or B, ℓ_{0} = 80) to the sample mean spectrum . The weights used in the fit are the entries of the sample covariance matrix. For each pair of (f_{M},α_{M}) values, we can derive the mean and covariance of (A^{XX},α_{XX}) from the fit.
Figure 2 shows the changes in the parameter p_{0} and the spectral indices α_{EE} and α_{BB} when varying either f_{M} or α_{M}, for N = 4 and 7. The points are the sample means of the parameters p_{0}, α_{EE}, and α_{BB} and the error bars represent the sample standard deviation. The results are compared to the data values reported in PXLIV and in PXXX. In PXLIV, the authors constrain the value of p_{0} with onepoint statistics of the p^{2} and ψ around the south pole at a fixed number of layers N (see middle plot of Fig. 10 of PXLIV). Over the range of values we consider, p_{0}, α_{EE}, and α_{BB} are mostly sensitive to α_{M}. The comparison of the power spectra between simulations and data does not constraint f_{M} nor N. The parameter f_{M} affects both the dispersion of ψ and p through depolarization along the line of sight (PXLIV). These two effects modify the variance of the dust polarization in opposite directions. The fact that the parameter p_{0} is independent of f_{M} (see top right panel of Fig. 2) suggests that they compensate each other over the range of values we are considering.
The measured values of α_{EE} and α_{BB} constrain α_{M} to be − 2.5 within about 0.1. For steeper B_{t} spectra (α_{M} ≤ −2.8), α_{EE} and α_{BB} are roughly constant with mean values lower than the observed values. In this regime, turbulence is not significant over the ℓ range used in this analysis. The dust total intensity map and the changing orientation with respect to the line of sight of the mean magnetic field dominate the variance of the polarized maps. For α_{M} ≥ −2.6, α_{EE} and α_{BB} are roughly equal to α_{M} within a small positive offset of about 0.1. In other words, the exponents of the dust polarization spectra reproduce the exponent of the magnetic field power spectrum.
The parameter p_{0} may also be used to constrain α_{M}. If the p_{0} values from PXLIV for N = 4 and 7 hold for the LR33 region, we find that the model fit constrains α_{M} to be − 2.5 within an uncertainty of about 0.1 (top left panel of Fig. 2). The systematic dependence of p_{0} with α_{M} follows from dispersion of the B_{t} orientation on angular scales corresponding to multipoles ℓ> 40. For a given f_{M}, this dispersion decreases as the power spectrum of B_{t} steepens (i.e., toward low values of α_{M}). Hence, the observed amplitude of the BB spectrum is matched for increasing values of p_{0} when α_{M} decreases.
Fig. 3 EE (top curve, diamond symbols) and BB (bottom curve, star symbols) power spectra of the simulated maps and their fits for the LR33 sky region. The diamonds and the stars represent the mean value computed over 1000 realizations. The 1σ error bars are derived from the sample standard deviation of the power in each ℓ bin. The blue dashed lines represent the fits to the mean spectra and the blue dotted lines the 1σ error on the fits. The red shade areas represent the powerlaw fit and the 1σ errors to the Planck data reported in PIPXXX for the LR33 region. 
Figure 3 shows the EE and BB power spectra for our fiducial values of α_{M}, f_{M}, and N. The points represent the mean value computed over 1000 realizations. The errors are derived from the sample variance of the power in each ℓ bin. The fit from the analysis of PXXX and its 1σ error are overplotted. The simulations are able to reproduce the EE and BB dust power spectra. The asymmetry parameter f has a value smaller than unity. The factor f^{2} = 0.55 is close to the value of R_{BB} = 0.48 (Table 1). Within this model, unlike for that of Ghosh et al. (2017), the TE correlation accounts for only a small part of the EB asymmetry.
Figure 4 shows ratios between the different power spectra of the simulated maps. Each point represents the sample mean of the 1000 ratios , and of each bin and the error bars represent the sample standard deviation. For comparison, we plot the input values and uncertainties of the R_{TT}, R_{TE}, and R_{BB} ratios. The ratios computed on the simulated maps are consistent with the input values, as expected because the maps were constructed in such a way that their power spectra respect that covariance structure.
Fig. 4 Three ratios R_{TT}, R_{TE}, and R_{BB} computed on the simulated maps for the LR33 sky region. The blue diamonds are the mean ratios for each ℓ bin computed over 1000 realizations. For the R_{TT} and R_{TE} ratios, the dashed line and the light and dark gray regions represent the input value and the 1 and 2σ errors on the input value, respectively. For the R_{BB} ratio, the dashed line, the light and dark gray regions represent, respectively, the ratio between the fits of the BB and EE data spectra from PXXX and their 1 and 2σ uncertainties. 
5.2. Statistics of the power spectrum amplitudes
Our simulations allow us to compute the dispersion of the dust BB power within a given ℓ bin. Although the dust maps are computed from Gaussian realizations of the turbulent field, the various processes involved in the computation might make them nonGaussian. For example, we do not expect the distribution of the power at multipole ℓ to tend to a Gaussian distribution for ℓ → ∞ as quickly as it would for a Gaussian random field. For the same reason, the variance of the distribution of the power for a given ℓ is not necessarily the cosmic variance.
Figure 5 shows the distribution of the power within one multipole bin around ℓ = 110 with a bin width of Δℓ = 20. The power spectra were computed for the LR33 region for which the covered sky is roughly equally distributed around the north and south Galactic poles. The figure also presents a Gaussian fit to the histogram and the expected cosmic variance for the same bin if the maps were drawn from a Gaussian random field on the sphere. The actual dispersion is a few times larger than the cosmic standard deviation. This effect might be due to the nonstationarity of the intensity map. The LR33 region includes some bright structures in dust total intensity. These localized structures are likely to be the explanation for the enhanced dispersion in the simulations. If this is the right interpretation, the enhancement must apply to the true sky because we are using the Planck total dust intensity map in our model.
Fig. 5 Distribution of the power per bin computed on simulated maps is significantly broader than the cosmic variance. The solid red line represents a Gaussian fit to the distribution in the multipole bin ℓ = 110 with a width Δℓ = 20 (black histogram). The dashed line represents the distribution expected for a Gaussian random field in that same bin. 
In addition to the spread, we looked at the shape of the PDF of the power per bin. We made 10000 of our dust simulations and 10000 Gaussian random simulations. The power spectrum used to produce the Gaussian realizations is the sample mean power spectrum of the 10000 dust simulations. For the two cases, we computed the power spectra, binned them with a width of Δℓ = 20, and fitted a Gaussian function to the sample distribution. In both cases, the dust simulations and the Gaussian realizations, we see the same difference between the PDF of the power per bin and the Gaussian fit. We concluded that the shape of distribution of the power per bin of our simulations is very similar to that of a Gaussian random field.
5.3. Power variations over the sky
We now show that the simulations reproduce the Planck power spectra for the high latitude sky in general, not just for the specific sky region LR33 used as input. First, we compute the spectra of the simulated maps for the five other LRxx sky regions from PXXX. Second, as in PXXX, we compute the spectra for smaller sky patches at high Galactic latitude with f_{sky} = 1%. We compare the amplitudes of the simulations spectra with the Planck results.
The analysis of the simulations on the six regions provides six sample mean power spectra and their sample variances. The power spectra on each region are computed and are fitted in the same way as described in Sect. 5.1. In Table 3, we gather the results of the fits together with the corresponding Planck values collected from Table 1 of PXXX for comparison. Error bars on the data measurements are smaller than those of these noiseless simulations because the error bars on the simulation spectra contain the variance from multiple random realizations of the GMF that does not affect the data.
While the simulations are constructed such that they match the data on one particular sky region (LR33), Table 3 shows that they also agree with the data on the other five regions within a small difference, which we comment on below. The spectra amplitudes at ℓ = 80 increase with the sky fraction faster than what PXXX reported for the Planck data. This slight difference may arise from the fact that we assumed a fixed value of N independent of the dust total intensity and Galactic latitude. In models of stellar polarization data at low Galactic latitudes and in molecular clouds, Jones et al. (1992) and Myers & Goodman (1991) assumed that N scales linearly with the dust column density. While their model hypothesis would not work for the diffuse ISM, we could consider variations in N. Alternatively, the slight difference in scaling could come from another simplifying assumption of the method, as we ignore the variation of the mean GMF orientation with distance from the Sun. It will be possible to modify our model to test these two ideas but this is beyond the scope of the present paper.
Results of powerlaw fits to the power spectra computed on simulated dust maps for the six Galactic regions from PXXX.
For the analysis on the 1% sky patches, we perform the same procedure as in PXXX to derive the empirical law between the amplitude at ℓ = 80 of the power spectra and the total intensity. Figure 6 shows the amplitudes of the EE and BB spectra as a function of the mean intensity of each patch. We realized 100 simulations and each vertical black line represents the sample mean and sample dispersion amplitude of one 400deg^{2} patch. The empirical law derived from a linear fit in the log (I_{353})−log (A^{XX}) space is overplotted. From this fit, we find a slope value of 2.15 ± 0.03 for the EE spectrum and of 2.09 ± 0.03 for the BB spectrum. The values of the slopes are slightly larger than 1.9 ± 0.02, which is the value that was measured on the Planck data. This difference for the patches is similar to that observed for the large sky regions, where the amplitude of the power spectra increases with f_{sky} slightly faster in the simulation than in the data (see Table 3).
In PXXX, the authors found that the cosmic variance and their measurement uncertainties were not large enough to account for the dispersion around the fit of Fig. 6. For our simulations, the spread of the distribution of the power in a given ℓ bin shown in Fig. 5 can explain the scatter observed around the fit of Fig. 6, which is comparable to that seen in the data. As detailed in Sect. 5.2, the scatter in the model comes mainly from the turbulent component of the magnetic field. In particular, we checked that the spread around the line fit is correlated with the mean polarization fraction, which depends on the mean orientation of the magnetic field over a given sky patch.
Fig. 6 Amplitudes of the power spectra (A^{EE} and A^{BB}) plotted vs. the mean total dust intensity at 353GHz computed on each 1% sky region for 100 realizations. Each vertical black line represents the sample mean and sample standard deviation of one of the 400deg^{2} patch. The blue dashed and dotted lines represent the powerlaw fit with the 1σ uncertainty of the simulation results. For comparison the red line is the same fit to the Planck data for the same set of sky patches. 
6. Multifrequency simulations
So far we have discussed ways to simulate structures on the sky at a single reference frequency. Component separation methods for CMB experiments rely on multifrequency data. A common approach to multifrequency simulations is to simulate the sky structure and the SED separately. The SED can be simulated using templates or analytical forms relying on a set of parameters, such as a modified blackbody law. The simulated sky map at a given frequency is then extrapolated to other frequencies. This method could also be applied to our simulations, but it does not permit us to control the decorrelation between maps at different frequencies in harmonic space, which is a characteristic crucial for component separation as discussed in Planck Collaboration Int. L (2017). Indeed, the decorrelation has an impact on the relative weights between the principal foreground modes. Here we present a method for multifrequency simulations constrained to match a given set of auto and crosspower spectra.
6.1. Method
We follow a procedure close to that commonly used to compute pseudorandom Gaussian vectors with a desired covariance from vectors with unit covariance. We realize as many simulated dust polarization maps as the desired number of frequencies and rearrange them to form a new set of maps such that the covariance structure of the latter is exactly as wanted.
To build a set of N_{f} maps at frequencies , we proceed as follows:
 1.
Simulate N_{f} singlefrequency maps obtained as described in Sect. 3, whose polarization spherical harmonic coefficients are gathered in a 2N_{f}dimension (E and B for N_{f} maps) vector x_{ℓm} for each pair (ℓ,m).
 2.
Compute the auto and crosspower spectra of the maps and gather them in a matrix Σ_{ℓ}, which is 2N_{f} × 2N_{f} at each multipole ℓ.
 3.
Specify a covariance structure of the maps over the range of N_{f} frequencies in the form of a 2N_{f} × 2N_{f} matrix for each multipole ℓ.
 4.
For each multipole, ℓ, compute the Cholesky decomposition of Σ_{ℓ} and , i.e., where the superscript † denotes the transposition.
For each pair (ℓ,m), construct the 2N_{f}dimension vector (13)
6.2. Results
We applied the procedure to produce a multifrequency set of maps (I(ν,j),Q(ν,j),U(ν,j)) where ν = 70, 100, 143, 217, 353 GHz and j = 1...N_{p} is the pixel index. For both EE and BB, the diagonal of the imposed covariance is , where ν_{0} = 353 GHz, T_{0} = 19.6K, β = 1.6 (Planck Collaboration Int. XXII 2015) and is the power spectrum of simulations at frequency ν_{0}. The SEDindependent correlation ratio between two frequencies ν_{1} and ν_{2} is set to 1 below ℓ = 30 and set by the following equation above ℓ = 30: (14)This dependence applies if the variations of the SED can be parametrized with a spatially varying spectral index (Appendix B). The parameter σ is set in such a way that the correlation between the 353 and 217GHz channels is 0.9 within the range of values measured on Planck data (Planck Collaboration Int. L 2017). We then construct the SED map from (15)and compute the mean SED α_{ν} from (16)In Fig. 7, we plot the distribution of for each ν = 70,100,143, and 217 GHz. As expected, the distribution widens with the separation between ν and ν_{0} because the correlation coefficient R_{ℓ} decreases. The correlation between the normalized SED of the same four frequencies is given by
This matrix gives an estimation of the coherence of the normalized SED through frequencies. We do not control the way the SED of a given sky pixel varies with respect to the mean SED because we model the decorrelation in harmonic space statistically.
Fig. 7 Distribution of the relative difference to the mean SED, normalized to 1 at 353 GHz in the 217 GHz (red, the narrower), 143 GHz (yellow), 100 GHz (green), and 70 GHz (blue, the wider) maps. 
7. Astrophysical perspective
Our paper has so far focused on our contribution to component separation for CMB data analysis. We presented a phenomenological model that can be used to simulate dust polarization maps, which statistically match Planck observations and are noisefree. In this section, we discuss what we learn about the GMF in the local interstellar medium from the modeling of the dust polarization power spectra. We examine our model results from this astrophysical perspective. We also compare our results with those of a companion paper Ghosh et al. (2017), which uses Hi data to account for the multiphase structure of the diffuse ISM. In Sect. 7.1, we briefly review Planck power spectra of dust polarization and our model fit. In Sects. 7.2 and 7.3, we discuss the power spectrum of the GMF and its correlation with matter.
Fig. 8 High resolution power spectra on the LR63 region; simulation vs. data. From top to bottom: TE, EE, and BB power spectra of the Planck353GHz, CMBcorrected maps (black), and one high resolution (N_{side} = 2048, FWHM = 10′) realization of the model (red). 
7.1. Model fit of the dust polarization spectra
We computed one simulation at an angular resolution of 10′ (ℓ ≃ 1000) to illustrate the model fit of the Planck data over a wider range of multipoles than in Sect. 5. The TE, EE, and BB spectra are presented in Fig. 8. The data spectra are crossspectra computed over the LR63 region using the two halfmission maps at 353GHz of Planck (Planck Collaboration I 2016; Planck Collaboration VIII 2016) after subtraction of the corresponding halfmission SMICA CMB maps (Planck Collaboration IX 2016). The simulation is built for our fiducial parameters of the turbulence. The values of the four parameters (t,p_{0},ρ,f) were determined for this sky region and this specific realization to be (1.01,0.22,0.20,0.74).
The spectra in Fig. 8 are consistent with a single spectral exponent over multipoles 40 ≤ ℓ ≤ 1000. At ℓ> 1000, the Planck spectra are dominated by the noise variance. At ℓ< 40, the spectra we computed with the publicly available maps are not reliable due to uncorrected systematics. The EE and BBPlanck 353GHz spectra computed down to ℓ = 2 after systematics corrections are presented in Planck Collaboration Int. XLVI (2016). These spectra shown in their Fig. 20 indicate a flattening at ℓ< 20, which is more pronounced for EE than for BB; the E to B power ratio goes from about 2 to 1 toward low multipoles.
An effective distance to the emitting dust is necessary to convert multipoles into physical scales. Over the high Galactic latitude region LR63, we estimate the distance of the emitting dust to be in the range 100–200pc. This estimate is constrained by the distance to the edge of the local bubble (Lallement et al. 2014) and the scale height of the dust emission, 200pc at the solar Galactocentric radius from the model of Drimmel & Spergel (2001). For the upper value of this distance range, the multipole range 40–1000 corresponds to linear scales from 0.5 to 15pc.
7.2. Galactic magnetic field power spectrum
Three of the model parameters we use – f_{M}, N and p_{0} – were constrained in PXLIV. Within these constraints, we find that our model fits the dust polarization power spectra for a spectral exponent of the power spectrum α_{M} = −2.5 ± 0.1 (Sect. 5.1). Within the quoted uncertainty, this value matches the spectral exponent of − 2.42 ± 0.02 of the Planck dust that is measured over the same range of multipoles on the EE and BB353GHz Planck spectra. Thus, a main conclusion of our modeling is that the exponent of the dust polarization spectra is that of the spectrum. The same conclusion is reached by Ghosh et al. (2017) for a distinct modeling of the polarization layers. This conclusion holds within the common framework of these two models and the corresponding assumptions.
The spectral exponent α_{M} we derive from the data fit is significantly larger than the Kolmogorov value of − 11/3 that is the common reference in interstellar turbulence (Brandenburg & Lazarian 2013), which is observed to apply to the electron density over a huge range of physical scales (Armstrong et al. 1995; Chepurnov & Lazarian 2010). A similar difference has been reported for the GMF spectrum derived over a similar range of scales from the analysis of synchrotron emission (e.g., Iacobelli et al. 2013) and of Faraday rotation measures (Oppermann et al. 2012). As discussed theoretically for synchrotron emission by Chepurnov (1998) and Cho & Lazarian (2002), a shallower slope is expected for ℓ multipoles approaching , where L_{max} is the length of the emitting layer along the line of sight and L_{out} the outer scale of turbulence. Two given lines of sight cross independent turbulent cells when their separation angle approaches the angle ~. It is only for smaller separation angles that the power spectrum of the emission reflects that of the magnetic field. This explanation put forward for synchrotron emission and Faraday rotation in earlier studies could apply to our analysis of dust polarization too. The flattening observed at ℓ< 20 in the spectra presented by Planck Collaboration Int. XLVI (2016) supports this interpretation, but the Planck data do not have the sensitivity to fully test it by checking whether the dust polarization spectra steepen at ℓ> 1000. Alternatively, the exponent of the GMF spectra might follow from the correlation of the magnetic field with interstellar matter. Indeed, Ghosh et al. (2017) find an exponent of − 2.4 for the E map they computed assuming a perfect alignment between the magnetic field and filamentary structure of their cold neutral medium Hi map.
7.3. Correlation between matter and the GMF
In this section we relate the structure of the GMF to that of the gas density in the diffuse ISM. The two are expected to be correlated to the extent that the magnetic field is frozen in matter. We note that this assumption might not hold everywhere (Eyink et al. 2013). The dust total intensity at 353GHz is a tracer of interstellar matter within some limitations characterized in a number of studies (e.g., Planck Collaboration Int. XVII 2014; Planck Collaboration XI 2014), which are not a main concern for this discussion. The spectrum of the GMF we find is close to that measured for the dust total intensity. Over the same ℓ range, Planck Collaboration Int. XLVIII (2016) report an exponent of − 2.7 for the TT spectrum of their 353GHz map corrected for CIB anisotropies, and Ghosh et al. (2017) report a value of − 2.6 for their total dust intensity map built from Hi data.
Dust polarization data have been used to quantify the alignment of the magnetic field orientation with the filamentary structure of the diffuse ISM (Clark et al. 2014; Planck Collaboration Int. XXXII 2016). This is a striking facet of the correlation between matter and the GMF, which creates TE correlation and thereby EB power asymmetry (Clark et al. 2015; Planck Collaboration Int. XXXVIII 2016). Ghosh et al. (2017) presented a model of dust polarization where this correlation between matter and the GMF applies to one single polarization layer that is associated with the cold neutral medium as traced by narrow Hi spectral lines. In their model that layer accounts for both the TE correlation and the EB asymmetry measured over the sky region with the lowest dust column density in the southern sky they analyzed. In our model, the TE correlation is introduced by adding one dust emission layer, where polarization is only in Emodes and is fully correlated to the T map. This corresponds to the additive term proportional to the ρ parameter in the second equation in Eqs. (5). The dust filamentary structures are present in all layers and the polarization results from the addition of the signals. We checked on the simulated images that this process introduces a preferred alignment between the filamentary structure of the T map and the magnetic field orientation inferred from the polarization angle, but this alignment is not as tight as that reported by Planck Collaboration Int. XXXVIII (2016) from their analysis of the most conspicuous filaments at high galactic latitudes in the Planck data. This difference comes from the fact that we use the same intensity map for each layer. Planck Collaboration Int. XXXVIII (2016) shows that the filamentary structure of the cold neutral medium has a main contribution to the EB asymmetry but it does not exclude a significant contribution related to the generic anisotropy of MHD turbulence, as suggested by Caldwell et al. (2017). We stress here that our modeling of the EB power asymmetry is mathematical. It does not constrain its physical origin. In this respect, our model is a framework that we are using to match the data statistically, but without a predictive power for astrophysics.
8. Conclusion
We introduced a process to simulate dust polarization maps, which may be used to statistically assess component separation methods in CMB data analysis. We detailed the simulation of dust polarization maps at one frequency before we introduced a mathematical means to produce maps at several frequencies and matched a given set of auto and crossspectra. Our method and the main results obtained by analyzing the simulated maps are summarized here.
Our approach builds on earlier studies, i.e., the analysis of Planck dust polarization data and the model framework from PXLIV, which relate the dust polarization sky to the structure of the GMF and interstellar matter. The structure of interstellar matter is the dust total intensity map from Planck. The GMF is modeled as a superposition of a mean uniform field and a Gaussian random (turbulent) component with a powerlaw power spectrum of exponent α_{M}. The integration along the line of sight performed to compute the Stokes maps is approximated by a sum over a small number of emitting layers with different realizations of the random GMF component. The mean field orientation, the amplitude of the random GMF component with respect to the mean component, the spectral exponent α_{M}, and the number of polarization layers are parameters common to the model from PXLIV. To match the power spectra of dust polarization measured with the Planck data, we add two main parameters (ρ and f) that introduce mathematically the TE correlation and EB power asymmetry. They are determined by fitting the Planck 353GHz power spectra for ℓ> 40 on one sky region at high Galactic latitude, LR33 from PXXX.
The model allows us to compute multiple realizations of the Stokes Q and U maps for different realizations of the random component of the magnetic field and to quantify the dispersion of dust polarization spectra for any given sky area away from the Galactic plane. The simulations reproduce the scaling laws between the dust polarization power and the mean total dust intensity from Planck, including the observed dispersion around the mean relation.
This paper discusses what we learn about the GMF in the local interstellar medium from the modeling of the dust polarization power spectra. We find that the slopes of the EE and BB power spectra of dust polarization measured by Planck are matched for α_{M} = −2.5 ± 0.1. As in Ghosh et al. (2017), we find that, for our model, the exponent of the spectrum of is very close to that of the dust polarization spectra. This exponent is larger than the Kolmogorov value of − 11/3 but close to that measured for matter (− 2.7), over the same region and range of multipoles (ℓ = 40−1000), using the Planck dust total intensity at 353GHz as a tracer. Our model does not allow us to comment on the origin of the TE correlation and EB asymmetry.
It would be possible to extend the model we presented in several ways, which might lead to fruitful explorations. To fit dust polarization spectra down to the very low multipoles relevant for measuring E and Bmode CMB polarization associated with the Universe reionization, we might need to account for the injection scale of turbulence. Phenomenologically, this could be carried out by introducing a lowℓ cutoff in the power spectrum of the magnetic field in Eq. (4).
Further model changes could also provide a better match to the data, in particular toward low Galactic latitudes. We have used a constant orientation for the mean GMF. A 3D model of the density structure of the Galactic ISM can be used to assign distances to the shells, and, thereby, to take into account the 3D structure of the largescale magnetic field, as in, for example, Fauvet et al. (2011) and Planck Collaboration Int. XLII (2016). In this case the intensity maps will differ for each layer and the effective number of layers could be allowed to vary with, for example, Galactic latitude or dust column density. In such a model, it would be possible to introduce, for each layer, the correlation between matter and the GMF and distinct dust SEDs. This method of introducing the decorrelation of dust polarization maps with frequency might in essence better represent the lineofsight averaging of polarization data (Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017) than the mathematical means proposed here. Finally, our paper focuses on dust polarization but a similar approach could be applied to produce maps of synchrotron polarization that match the observed correlation with dust polarization (Planck Collaboration Int. XXII 2015; Choi & Page 2015).
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
Xpol is an algorithm for power spectrum estimation that is an extension to polarization of the Xspect method (Tristram et al. 2005).
These masks are available on the Planck Legacy Archive as HFI_Mask_GalPlaneapo2_2048_R2.00.fits and described in the Planck Explanatory Supplement 2015 accessible at the web page https://wiki.cosmos.esa.int/planckpla2015/index.php/Frequency_Maps#Galactic_plane_masks
Acknowledgments
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013)/ERC grant agreement No. 267934.
References
 Abazajian, K. N., Arnold, K., Austermann, J., et al. 2015, Astropart. Phys., 63, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Andersson, B.G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209 [NASA ADS] [CrossRef] [Google Scholar]
 BICEP2/Keck Array and Planck Collaborations 2015, Phys. Rev. Lett., 114, 101301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Brandenburg, A., & Lazarian, A. 2013, Space Sci. Rev., 178, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Caldwell, R. R., Hirata, C., & Kamionkowski, M. 2017, ApJ, 839, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Chepurnov, A. V. 1998, Astronomical and Astrophysical Transactions, 17, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Chepurnov, A., & Lazarian, A. 2010, ApJ, 710, 853 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, J., & Lazarian, A. 2002, ApJ, 575, L63 [NASA ADS] [CrossRef] [Google Scholar]
 Choi, S. K., & Page, L. A. 2015, JCAP, 12, 020 [NASA ADS] [CrossRef] [Google Scholar]
 Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Clark, S. E., Peek, J. E. G., & Putman, M. E. 2014, ApJ, 789, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Clark, S. E., Hill, J. C., Peek, J. E. G., Putman, M. E., & Babler, B. L. 2015, Phys. Rev. Lett., 115, 241302 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dolginov, A. Z., & Mitrofanov, I. G. 1976, Ap&SS, 43, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Drimmel, R., & Spergel, D. N. 2001, ApJ, 556, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Dunkley, J., Amblard, A., Baccigalupi, C., et al. 2009, in AIP Conf. Ser., 1141, eds. S. Dodelson, D. Baumann A. Cooray, et al., 222 [Google Scholar]
 Eyink, G., Vishniac, E., Lalescu, C., et al. 2013, Nature, 497, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Fauvet, L., MacíasPérez, J. F., Aumont, J., et al. 2011, A&A, 526, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghosh, T., Boulanger, F., Martin, P. G., et al. 2017, A&A, 601, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Guth, A. H. 1981, Phys. Rev. D, 23, 347 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hoang, T., & Lazarian, A. 2016, ApJ, 831, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Houde, M., Vaillancourt, J. E., Hildebrand, R. H., Chitsazzadeh, S., & Kirby, L. 2009, ApJ, 706, 1504 [NASA ADS] [CrossRef] [Google Scholar]
 Iacobelli, M., Haverkorn, M., Orrú, E., et al. 2013, A&A, 558, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jansson, R., & Farrar, G. R. 2012, ApJ, 757, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, T. J., Klebe, D., & Dickey, J. M. 1992, ApJ, 389, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Kalberla, P. M. W., Kerp, J., Haud, U., et al. 2016, ApJ, 821, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., & Kovetz, E. D. 2016, ARA&A, 54, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Lallement, R., Vergely, J.L., Valette, B., et al. 2014, A&A, 561, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lazarian, A., & Pogosyan, D. 2012, ApJ, 747, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, H. M., & Draine, B. T. 1985, ApJ, 290, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A. D. 1982, Phys. Lett. B, 108, 389 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Myers, P. C., & Goodman, A. A. 1991, ApJ, 373, 509 [NASA ADS] [CrossRef] [Google Scholar]
 O’Dea, D. T., Clark, C. N., Contaldi, C. R., & MacTavish, C. J. 2012, MNRAS, 419, 1795 [NASA ADS] [CrossRef] [Google Scholar]
 Oppermann, N., Junklewitz, H., Robbers, G., et al. 2012, A&A, 542, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XVII. 2014, A&A, 566, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XX. 2015, A&A, 576, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXI. 2015, A&A, 576, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXII. 2015, A&A, 576, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXII. 2016, A&A, 586, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLII. 2016, A&A, 596, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLIV. 2016, A&A, 596, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVI. 2016, A&A, 596, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. L. 2017, A&A, 599, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seljak, U., & Zaldarriaga, M. 1996, ApJ, 469, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Starobinskiǐ, A. A. 1979, Soviet J. Exp. Theor. Phys. Lett., 30, 682 [Google Scholar]
 Stein, W. 1966, ApJ, 144, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90 [NASA ADS] [CrossRef] [Google Scholar]
 Tristram, M., MaciasPerez, J. F., Renault, C., & Santos, D. 2005, MNRAS, 358, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Waelkens, A., Jaffe, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A, 495, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yang, M. 2008, Applied Economics Letters, 15, 737 [CrossRef] [Google Scholar]
 Zaldarriaga, M. 2001, Phys. Rev. D, 64, 103001 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Implementation of the method
In this appendix, we explain how we compute the dust polarization maps used in this paper. Sect. A.1 presents the procedure we use to derive the stage a maps using the framework in Sect. 2. Section A.2 describes how we produce the stage b maps that match the dust TE correlation and EB asymmetry measured by Planck, using the method described in Sect. 3.
Appendix A.1: Stage a maps
We explain how we produce the (I_{a},Q_{a},U_{a}) maps at a reference frequency ν_{0}, which we choose to be 353GHz, the bestsuited Planck channel to study dust polarization. These maps have no TE correlation and no EB asymmetry at ℓ> 40.
The intensity map I_{a} is not computed from Eqs. (2) but derived from observations. We use I_{a} = D_{353}, where D_{353} is the dust total intensity map at 353GHz of Planck Collaboration Int. XLVIII (2016) after separation from the CIB and CMB anisotropies. To compute Q_{a} and U_{a}, we need the set of angle maps γ_{i} and ψ_{i}, which determine the orientation of the magnetic field in the N layers. For each layer, we draw an independent Gaussian realization for each of the three components of in Eq. (3). The angle maps are computed for the total magnetic field B including the mean magnetic field B_{0}. With the set of angles maps γ_{i}, using the Stokes I equation in Eqs. (2), we compute the map S_{i}(ν) at the frequency ν_{0}, (A.1)where S_{i}(ν_{0}) has been assumed to be independent of the index i and p_{0} is set to a fiducial value of 0.25. Next, we combine S_{i}(ν_{0}) and the angle maps γ_{i} and ψ_{i} in the Stokes Q and U equations in Eqs. (2) to compute the ratio maps Q_{a}/ (p_{0} × I_{a}) and U_{a}/ (p_{0} × I_{a}) at the frequency ν_{0}. These ratio maps are independent of I_{a} and depend on p_{0} only through S_{i}(ν_{0}). They are computed at pixel resolution defined by the N_{side} = 256 HEALPix parameter. After multiplication by D_{353}, we obtain the two maps Q_{a}/p_{0} and U_{a}/p_{0}, which have an illdefined beam transfer function. The D_{353} map has a resolution that varies across the sky. We overcome this issue by smoothing (I_{a},Q_{a},U_{a}) to a resolution lower than the lowest resolution of the D_{353} map. The model maps used in the paper have N_{side} = 256 and a symmetric Gaussian beam with a full width at half maximum of 30′.
Appendix A.2: Stage b maps
From the harmonic coefficients of Eq. (5), we compute the power spectra of stage b maps at a given multipole ℓ, as functions of t,p_{0},ρ,f, and x, where (A.2), , and are the power spectra of stage a maps and ℰ[·] is an averaging over multipoles between ℓ = 60 and ℓ = 200. When the slope of the TT and polarization spectra are close to one another, the ratio x is close to being independent of multipole ℓ. Since this simplification approximately applies for dust emission (PXXX), we consider the ratio x to be constant over the relevant multipole range.
Assuming for X ≠ Y, the ratios of Eq. (6) can be expressed as follows: (A.3)where and . When the ratios R_{XY} are chosen, then the system A.3 becomes a system of equations in { f,y,z }. Although the system is not linear, it can be inverted, as long as , i.e., . When choosing values for the ratios R_{XY}, this condition has to be satisfied because the power spectra form a covariance, which must be positive definite. Restricting the set of solutions to positive reals, there is a unique solution, i.e., (A.4)From the solution of Eq. (A.4) and the normalization factor N_{B} = (p_{0}f)^{2}, we can compute the parameters ρ,f,t, and p_{0} as follows: (A.5)We note that ρ and the correlation coefficient between the T and E parts of stage b maps, noted , are related as follows: (A.6)We choose the BB normalization factor N_{B} such that the power spectrum of stage a divided by p_{0} map is adjusted to the fit of the power spectrum measured over the region LR33 in PXXX (noted ). Following the notation of PXXX, we have , where the values of the parameters and A^{BB,data} are taken from Table 1. In the case where N_{B} is ℓindependent, N_{B} is the solution of the minimization of the following chisquared (A.7)with the variance of , estimated from Monte Carlo simulations and (ℓ_{1},ℓ_{2}) = (60,200) as for x. The fit also provides the standard deviation on the normalization factor N_{B}.
Appendix A.3: Summary of the procedure
The following points sketch the procedure to produce our simulations:

1.
Draw Stokes maps Q_{a} and U_{a} divided by p_{0} as described in Sect. A.1.

2.
Mask the Galactic plane and compute the harmonic coefficients a_{ℓm}.

3.
Given a mask, compute the full sky power spectra .

4.
Evaluate x as defined in Eq. (A.2) and the BB normalization N_{B}.

5.
Choose values for the ratios R_{XY} of Eq. (6).

6.
Compute the corresponding solutions {f,y,z} of Eqs. (A.4).

Compute the parameters ρ,f,t, and p_{0} of Eqs. (A.5).

8.
Construct the harmonic coefficients b_{ℓm} according to their definition of Eqs. (5).

9.
Transform the b_{ℓm}’s to (I_{b},Q_{b},U_{b}).
The stage b maps thus constructed feature the desired twopoint statistics on the desired region of the sky. The procedure can be applied on separate multipole bins, which then gives scaledependent parameters.
Appendix B: Decorrelation due to a variable spectral index
This appendix shows how to compute the decorrelation in harmonic space between two frequency maps, when spectral differences about a mean SED may be parametrized with a spatially varying spectral index. This appendix restricts the proof to the simple case where the map that is scaled through frequencies and the spectral index map are correlated Gaussian white noise maps.
Let f(n) and δβ(n) be two Gaussian random fields on the sphere such that From f(n) and δβ(n) we construct a set of maps at frequencies ν_{i},(B.5)where ν_{0} is a reference frequency and K_{i} possibly contains the mean SED and unit conversion factors.
The aim is to compute the crossspectrum ( ) between the different f_{i}(n), i.e., where the Y_{ℓm}(n) represent the spherical harmonics; we assumed that two directions of the maps are uncorrelated and that the maps are statistically isotropic. We can rewrite the product f_{i}(n)f_{j}(n) as follows: (B.8)where , g(n) = f(n) /σ_{f} and δγ_{ij}(n) = σ_{ij}δβ(n) /σ_{β} with . It can be easily verified that (B.9)so that the expression in brackets of Eq. (B.8) has a normal lognormal mixture distribution as parametrized in, for example, Yang (2008). Thus, (B.10)and (B.11)We note that the correlation does not depend on the mean SED.
All Tables
Results of powerlaw fits to the power spectra computed on simulated dust maps for the six Galactic regions from PXXX.
All Figures
Fig. 1 Probability distribution functions of p^{2} and ψ (top and bottom plots) for stage a and stage b maps (red and blue histograms). The maps were computed using the fiducial values of α_{M}, f_{M}, and N and the corresponding parameters t, p_{0}, ρ, and f introduced in Sect. 4.2. The distributions are computed on the southern Galactic polar cap (b ≤ −60°) as in PXLIV. The very close match between the corresponding histograms shows that the inclusion of the TE correlation and the EB asymmetry does not alter the onepoint statistics of the simulated maps. 

In the text 
Fig. 2 Parameter p_{0} (top) and slopes of the EE (middle) and BB (bottom) power spectra of stage b maps vs. the slope α_{M} of the power spectrum of the turbulent component of the magnetic field (left) and vs. the relative strength of the turbulence f_{M} (right). In the left plots f_{M} = 0.8 and in the right plots α_{M} = −2.5. Red stars (blue squares) represent results for N = 4 (N = 7). The abscissae of the two sets of points are slightly shifted from their original values for a better visibility. The values observed in the data are represented by a gray shaded region for α_{EE} and α_{BB} (a dashed line for the mean, dark, and light gray for the 1 and 2σ uncertainties) and by a hatched regions for p_{0} (a dashed horizontal line for the mean, and a red 45° (resp. blue −45°) hatched region for 1σ uncertainty for N = 4 (resp. 7)). 

In the text 
Fig. 3 EE (top curve, diamond symbols) and BB (bottom curve, star symbols) power spectra of the simulated maps and their fits for the LR33 sky region. The diamonds and the stars represent the mean value computed over 1000 realizations. The 1σ error bars are derived from the sample standard deviation of the power in each ℓ bin. The blue dashed lines represent the fits to the mean spectra and the blue dotted lines the 1σ error on the fits. The red shade areas represent the powerlaw fit and the 1σ errors to the Planck data reported in PIPXXX for the LR33 region. 

In the text 
Fig. 4 Three ratios R_{TT}, R_{TE}, and R_{BB} computed on the simulated maps for the LR33 sky region. The blue diamonds are the mean ratios for each ℓ bin computed over 1000 realizations. For the R_{TT} and R_{TE} ratios, the dashed line and the light and dark gray regions represent the input value and the 1 and 2σ errors on the input value, respectively. For the R_{BB} ratio, the dashed line, the light and dark gray regions represent, respectively, the ratio between the fits of the BB and EE data spectra from PXXX and their 1 and 2σ uncertainties. 

In the text 
Fig. 5 Distribution of the power per bin computed on simulated maps is significantly broader than the cosmic variance. The solid red line represents a Gaussian fit to the distribution in the multipole bin ℓ = 110 with a width Δℓ = 20 (black histogram). The dashed line represents the distribution expected for a Gaussian random field in that same bin. 

In the text 
Fig. 6 Amplitudes of the power spectra (A^{EE} and A^{BB}) plotted vs. the mean total dust intensity at 353GHz computed on each 1% sky region for 100 realizations. Each vertical black line represents the sample mean and sample standard deviation of one of the 400deg^{2} patch. The blue dashed and dotted lines represent the powerlaw fit with the 1σ uncertainty of the simulation results. For comparison the red line is the same fit to the Planck data for the same set of sky patches. 

In the text 
Fig. 7 Distribution of the relative difference to the mean SED, normalized to 1 at 353 GHz in the 217 GHz (red, the narrower), 143 GHz (yellow), 100 GHz (green), and 70 GHz (blue, the wider) maps. 

In the text 
Fig. 8 High resolution power spectra on the LR63 region; simulation vs. data. From top to bottom: TE, EE, and BB power spectra of the Planck353GHz, CMBcorrected maps (black), and one high resolution (N_{side} = 2048, FWHM = 10′) realization of the model (red). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.