Statistical simulations of the dust foreground to CMB polarization
Flavien Vansyngel, François Boulanger, Tuhin Ghosh, Benjamin D. Wandelt, Jonathan Aumont, Andrea Bracco, François Levrier, Peter G. Martin, Ludovic Montier
AAstronomy & Astrophysics manuscript no. final c (cid:13)
ESO 2017April 28, 2017
Statistical simulations of the dust foreground to cosmic microwavebackground polarization
F. Vansyngel , F. Boulanger , T. Ghosh , , B. Wandelt , , J. Aumont , A. Bracco , , F. Levrier , P. G. Martin , L.Montier (A ffi liations can be found after the references) Preprint online version: April 28, 2017
ABSTRACT
The characterization of the dust polarization foreground to the cosmic microwave background (CMB) is a necessary step toward the detection ofthe B -mode signal associated with primordial gravitational waves. We present a method to simulate maps of polarized dust emission on the spherethat is similar to the approach used for CMB anisotropies. This method builds on the understanding of Galactic polarization stemming from theanalysis of Planck data. It relates the dust polarization sky to the structure of the Galactic magnetic field and its coupling with interstellar matterand turbulence. The Galactic magnetic field is modeled as a superposition of a mean uniform field and a Gaussian random (turbulent) componentwith a power-law power spectrum of exponent α M . The integration along the line of sight carried out to compute Stokes maps is approximated bya sum over a small number of emitting layers with di ff erent realizations of the random component of the magnetic field. The model parametersare constrained to fit the power spectra of dust polarization EE , BB , and T E measured using
Planck data. We find that the slopes of the E and B power spectra of dust polarization are matched for α M = − .
5, an exponent close to that measured for total dust intensity but larger than theKolmogorov exponent -11 /
3. The model allows us to compute multiple realizations of the Stokes Q and U maps for di ff erent realizations of therandom component of the magnetic field, and to quantify the variance of dust polarization spectra for any given sky area outside of the Galacticplane. The simulations reproduce the scaling relation between the dust polarization power and the mean total dust intensity including the observeddispersion around the mean relation. We also propose a method to carry out multifrequency simulations, including the decorrelation measuredrecently by Planck , using a given covariance matrix of the polarization maps. These simulations are well suited to optimize component separationmethods and to quantify the confidence with which the dust and CMB B -modes can be separated in present and future experiments. We alsoprovide an astrophysical perspective on our phenomenological modeling of the dust polarization spectra. Key words.
Polarization – ISM: general – cosmology: comic background radiation – Galaxy: ISM – submillimeter: ISM
1. Introduction
An era of exponential expansion of the universe, dubbed cos-mic inflation, has been proposed to explain why the universeis almost exactly Euclidean and nearly isotropic (Guth 1981;Linde 1982). One generic prediction of this theoretical paradigmis the existence of a background of gravitational waves, whichproduces a distinct, curl-like, signature in the polarization ofthe cosmic microwave background (CMB), referred to as pri-mordial B -mode polarization (Starobinskiˇi 1979). The detec-tion of this signal would have a deep impact on cosmology andfundamental physics, motivating a number of experiments de-signed to measure the sky polarization at microwave frequen-cies. Current projects have achieved the sensitivity required todetect the CMB B -mode signal predicted by the simplest mod-els of inflation (Abazajian et al. 2015; Kamionkowski & Kovetz2015). Yet, any detection has relied on the proper removal ofmuch brighter Galactic foregrounds.Thermal emission from aspherical dust grains aligned withrespect to the Galactic magnetic field (GMF) is the dominantpolarized foreground for frequencies higher than about 70 GHz(Dunkley et al. 2009; Planck Collaboration X 2016). From theanalysis of the Planck
353 GHz polarization maps, we know Planck ( ) is a project of theEuropean Space Agency (ESA) with instruments provided by two sci-entific consortia funded by ESA member states and led by PrincipalInvestigators from France and Italy, telescope reflectors provided that the primordial B -mode polarization of the CMB can-not be measured without subtracting the foreground emission,even in the faintest dust-emitting regions at high Galactic lat-itude (Planck Collaboration Int. XXX 2016, hereafter PXXX).The observed correlation between the B -mode signal detectedby 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 polar-ization signals, CMB experiments must rely on multifrequencyobservations. Component separation is a main challenge becausethe spatial structure of dust polarization is observed to vary withfrequency (Planck Collaboration L 2016). This introduces twoquestions that motivate our work. What design of CMB exper-iments and combination of ground-based, balloon-borne, andspace observations is best to achieve an optimal separation? Howcan confidence in the subtraction of foregrounds be quantified?To provide quantitative answers, we must be able to simulateobservations of the polarized sky combining Galactic and CMBpolarization. This paper presents a statistical model with a fewparameters to simulate maps of dust polarization in a way similarto what is available for CMB anisotropies (Seljak & Zaldarriaga1996). through a collaboration between ESA and a scientific consortium ledand funded by Denmark, and additional contributions from NASA(USA). 1 a r X i v : . [ a s t r o - ph . C O ] A p r . Vansyngel et al.: Simulations of dust foreground to CMB polarization The analysis of the Wilkinson Microwave Anisotropy Probe(
WMAP ) data and the preparation of the
Planck project moti-vated a series of models of the polarized synchrotron and ther-mal dust emission at microwave frequencies (Page et al. 2007;Miville-Deschˆenes et al. 2008; Fauvet et al. 2011; O’Dea et al.2012; Delabrouille et al. 2013). These early studies followedtwo distinct approaches. The first is to produce a sky that isas close as possible to the observed sky combining data tem-plates and a spectral model. Prior to
Planck , for dust polar-ization this was performed using stellar polarization data byPage et al. (2007), and
WMAP observations of synchrotron po-larization by Delabrouille et al. (2013). More recently, the sim-ulations presented in Planck Collaboration XII (2016) use the
Planck
353 GHz data to model dust polarization. This first ap-proach is limited by the signal-to-noise ratio of available data,which for dust polarization is low at high Galactic latitudeeven after smoothing to one degree angular resolution. The sec-ond approach is to simulate the polarization sky from a 3Dmodel of the GMF and of the density structure of the interstel-lar medium (ISM), both its regular and turbulent components,as carried out by Miville-Deschˆenes et al. (2008), Fauvet et al.(2011) and O’Dea et al. (2012). This method connects the mod-eling of the microwave polarized sky to broader e ff orts tomodel the GMF (Waelkens et al. 2009; Jansson & Farrar 2012;Planck Collaboration Int. XLII 2016). Planck polarization maps have been used to charac-terize the structure (Planck Collaboration Int. XIX 2015;Planck Collaboration Int. XX 2015) and the spectral en-ergy distribution (SED) of polarized thermal emissionfrom Galactic dust (Planck Collaboration Int. XXI 2015;Planck Collaboration Int. XXII 2015). Several studies have es-tablished the connection between the structure of the magneticfield and matter (Clark et al. 2014; Planck Collaboration Int. XX2015; Planck Collaboration Int. XXXII 2016; Martin et al.2015; Kalberla et al. 2016). The power spectra analysispresented in PXXX decomposes dust polarization into E (gradient-like) and B (curl-like) modes (Zaldarriaga 2001;Caldwell et al. 2016). This analysis led to two unexpectedresults: a positive T E correlation and a ratio of about 2 betweenthe E and B dust powers over the (cid:96) range 40 to 600. Clark et al.(2015) and Planck Collaboration Int. XXXVIII (2016) haveshowed that the observed T E correlation and asymmetrybetween E − and B -mode power amplitudes for dust polarizationcould be both accounted for by the preferred alignment betweenthe filamentary structure of the total intensity map and theorientation of the magnetic field inferred from the polarizationangle.The work presented here makes use of the model frameworkintroduced in Planck Collaboration Int. XLIV (2016) (hereafterPXLIV). By analyzing Planck dust polarization maps towardthe southern Galactic cap, the part of the sky used for CMBobservations from Antartica and Atacama, PXLIV related thelarge-scale patterns of the maps to the mean orientation of themagnetic field, and the scatter of the dust polarization angle andfraction ( ψ and p ) to the amplitude of its turbulent component.In this paper, we extend their work to produce Stokes maps thatfit dust polarization power spectra including the T E correlationand the
T T / EE and EE / BB power ratios at high and intermedi-ate Galactic latitudes. In a companion paper Ghosh et al. (2016),the dust polarization of the southern sky region with the lowestdust column density is modeled using Hi observations and astro-physical insight to constrain their model parameters. In essence,our approach is more mathematical but it allows us to model dustpolarization over a larger fraction of the sky. The two approaches are complementary and compared in this paper. We also presenta mathematical process to introduce spatial decorrelation acrossmicrowave frequencies via the auto and cross spectra of dust po-larization. By doing this, we obtain a model to compute indepen-dent realizations of dust polarization sky maps at one or multiplefrequencies with a few parameters adjusted to fit the statisticalproperties inferred from the analysis of the Planck data awayfrom the Galactic plane.The paper is organized as follows. Sects. 2 and 3 presentthe framework we use to model dust polarization in generalterms. Our method is illustrated by producing simulated maps at353 GHz presented in Sect. 4. We show that these maps success-fully match the statistical properties of dust polarization derivedfrom the analysis of
Planck data (Sect. 5). One method to com-pute dust polarization maps at multiple frequencies is presentedin Sect. 6. We discuss the astrophysical implications of our workin Sect. 7. The main results of the paper are summarized inSect. 8. Appendix A details how the simulated maps used inthis study are computed. Appendix B shows how to computethe cross correlation between two frequency maps, when spec-tral di ff erences about a mean SED may be parametrized with aspatially varying spectral index.
2. Astrophysical framework
To model dust polarization we used the framework introducedby PXLIV, which we briefly describe here. We refer to PXLIVfor a detailed presentation and discussion of the astrophysicalmotivation and the simplifying assumptions of our modeling ap-proach.The polarization of thermal dust emission results from thealignment of aspherical grains with respect to the GMF (Stein1966; Lee & Draine 1985; Planck Collaboration Int. XXI 2015).Within the hypothesis that grain polarization properties, includ-ing alignment, are homogeneous, the structure of the dust polar-ization sky reflects the structure of the magnetic field combinedwith that of matter. We assume that this hypothesis applies tothe di ff use ISM where radiative torques provide a viable mech-anism to align grains e ffi ciently (Dolginov & Mitrofanov 1976;Andersson et al. 2015; Hoang & Lazarian 2016).To compute the Stokes parameters I , Q , and U describ-ing the linearly polarized thermal dust emission, we startfrom the integral equations in Sect. 3.2 and Appendix B ofPlanck Collaboration Int. XX (2015) for optically thin emissionat frequency ν , i.e., I ( ν ) = (cid:90) S ( ν ) (cid:34) − p (cid:32) cos γ − (cid:33)(cid:35) d τ ν ; Q ( ν ) = (cid:90) p S ( ν ) cos (2 φ ) cos γ d τ ν ; (1) U ( ν ) = (cid:90) p S ( ν ) sin (2 φ ) cos γ d τ ν . where S ( ν ) is the source function, τ ν the optical depth, p a pa-rameter related to dust polarization properties (the grain crosssections and the degree of alignment with the magnetic field), γ the angle that the local magnetic field makes with the planeof the sky, and φ the local polarization angle (see Fig. 14 inPlanck Collaboration Int. XX (2015)).As in PXLIV, the integration along the line of sight is ap-proximated by a sum over a finite number N of layers. This sum
2. Vansyngel et al.: Simulations of dust foreground to CMB polarization is written as I ( ν ) = N (cid:88) i = S i ( ν ) (cid:34) − p (cid:32) cos γ i − (cid:33)(cid:35) ; Q ( ν ) = N (cid:88) i = p S i ( ν ) cos(2 φ i ) cos γ i ; (2) U ( ν ) = N (cid:88) i = p S i ( ν ) sin(2 φ i ) cos γ i ;where S i ( ν ) is the integral of the source function over layer i ,and γ i and φ i define the magnetic field orientation within eachlayer. As discussed in PXLIV, the layers are a phenomenologi-cal means to model the density structure of the interstellar matterand the correlation length of the GMF. This approach accountsfor both signatures of the turbulent magnetic field component inGalactic polarization maps: the depolarization resulting from theintegration along the line of sight of emission with varying po-larization orientations, and the scale invariant structure of the po-larization maps across the sky reflecting the power spectrum ofthe turbulent component of the magnetic field (Cho & Lazarian2002; Houde et al. 2009). It overcomes the di ffi culty of gener-ating realizations of the turbulent component of the magneticfield in three dimensions over the celestial sphere. Ghosh et al.(2016) uses Hi data to associate the layers with di ff erent phasesof the ISM, each of which provide a di ff erent intensity map. Onthe contrary, in the simulations presented in this paper, like inPXLIV, the term S i ( ν ) in Eqs. 2 is a sky map assumed to be thesame in each layer, i.e., it is independent of the index i . Thuswe do not address the question of the physical meaning of thelayers.Through the angles γ i and φ i , the model relates the dust po-larization to the structure of the GMF. The magnetic field B isexpressed as the sum of its mean (ordered), B , and turbulent(random), B t , components, B = B + B t = | B | ( ˆ B + f M ˆ B t ) , (3)where ˆ B and ˆ B t are unit vectors in the directions of B and B t ,and f M a model parameter that sets the relative strength of therandom component of the field. To simulate dust as a foregroundto the CMB we need a description of the GMF within the so-lar neighborhood. We follow PXLIV in assuming that B hasa fixed orientation in all layers. We ignore the structure of theGMF on galaxy-wide scales because the dust emission arisesmainly from a thin disk with a relatively small scale height andwe are interested in modeling dust polarization away from theGalactic plane. This scale height is not measured directly in thesolar neighborhood but modeling of the dust emission from theMilky Way indicates that it is ∼
200 pc at the solar distance fromthe Galactic center (Drimmel & Spergel 2001).Each component of the vector field ˆ B t in 3D, in each layer, isobtained from independent Gaussian realizations of a power-lawpower spectrum, which is written as C (cid:96) ∝ (cid:96) α M for (cid:96) ≥ . (4)Our modeling of ˆ B t is continuous over the celestial sphere anduncorrelated between layers. The coherence of the GMF orien-tation along the line of sight comes from the mean field and iscontrolled by the parameter f M .The model has six parameters: the Galactic longitude andlatitude l and b defining the orientation of ˆ B , the factor f M , the number of layers N , the spectral exponent α M , and the e ff ectivepolarization fraction of the dust emission p . The PXLIV authorsused the same model to analyze the dust polarization measuredby Planck at 353 GHz over the southern Galactic cap (Galacticlatitude b < − ◦ ). They determined l = ◦ ± ◦ and b = ◦ ± ◦ by fitting the large-scale pattern observed in the Stokes Q and U maps, and f M = . ± . N = ± p = ± p , thesquare of the dust polarization fraction p , and of the polarizationangle ψ , computed after removal of the regular pattern from theordered 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 modelmaps have equal EE and BB power, and no T E correlation at (cid:96) (cid:38)
30. This follows from the fact that our modeling does not in-clude the alignment observed between the filamentary structureof the di ff use ISM and the GMF orientation. Some T E corre-lation is present at low (cid:96) because the mean GMF orientation isclose to being within the Galactic disk, and we take into accountthe latitude dependence of the total dust intensity. In the nextsection, we explain how we modify the spherical harmonic de-composition of the stage a maps to introduce the T E correlationand the E - B asymmetry, matching the Planck dust polarizationpower spectra in PXXX.
3. Introducing TE correlation and E-B asymmetry
Our aim is to simulate maps that match given observablesbased on dust angular power spectra, namely the
T E correlation,
T T / EE and EE / BB ratios, and the BB spectrum without alteringthe statistics of p and ψ of the Stokes maps from stage a . Wedescribe 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 beapplied on a full, or a masked, sky.We start with the Stokes maps ( I a , Q a , U a ) obtained as de-scribed in Sect. 2. We compute the spherical harmonic coe ffi -cients of the stage b maps from those of the stage a maps asfollows: b T (cid:96) m = ta T (cid:96) m b E (cid:96) m = p ( a E (cid:96) m / p + ρ a T (cid:96) m ) b B (cid:96) m = p ( f a B (cid:96) m / p ) . (5)where a X (cid:96) m and b X (cid:96) m denote the coe ffi cients of the X = T , E , B har-monic decomposition of stage a and b maps, respectively. Theparameter ρ introduces the T E correlation and the factor f the E - B asymmetry. The parameter t is a scaling factor for the in-tensity part and p is the polarization parameter introduced inEqs. 1. These parameters control the amplitude of the T T , EE , BB , and T E power spectra of stage b maps. We note that Q a and U a scale linearly with p and thus that the two ratios a E (cid:96) m / p and a B (cid:96) m / p in Eqs. 5 are independent of p .At this stage b , our modeling of the random compo-nent of the magnetic field is anisotropic, which is a fun-damental 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 pa-rameter ρ , has a filamentary structure and thus is anisotropic.This amounts to adding an extra polarization layer that is per-fectly aligned with the filamentary structure of the matter and issimilar to what is carried out by Ghosh et al. (2016) for their coldneutral medium map. Second, the factor f breaks the symmetrybetween E and B , whereas the power is expected to be equally
3. Vansyngel et al.: Simulations of dust foreground to CMB polarization distributed between E and B modes in the case of isotropic turbu-lence (Caldwell et al. 2016). Through the parameter f , the ran-dom component of B is anisotropic in all layers and everywhereon the sky, unlike in Ghosh et al. (2016) where anisotropy is in-troduced in only one layer.In the simplest case, ρ, f , t , and p are constants over thewhole multipole range and in the most general case they arefunctions of (cid:96) and m . We find that the statistics of p and ψ foundusing the stage a maps are lost at stage b if f (cid:44) ρ (cid:44) t and p are constants but f and ρ depend on (cid:96) and tendtoward 1 and 0 for very low (cid:96) values, respectively.The power spectra of stage b maps are noted C XY (cid:96) with X , Y = T , E , B and use the quantity D XY (cid:96) ≡ (cid:96) ( (cid:96) + C XY (cid:96) / (2 π ). The t , p , ρ, and f coe ffi cients in Eqs. 5 are chosen such that the powerspectra of stage b maps match a given set of averaged ratios asfollows: R TT ≡ E (cid:104) D TT (cid:96) / D EE (cid:96) (cid:105) , R TE ≡ E (cid:104) D TE (cid:96) / D EE (cid:96) (cid:105) , (6) R BB ≡ E (cid:104) D BB (cid:96) / D EE (cid:96) (cid:105) , where E [ · ] is a given averaging process over multipoles. Theabsolute scaling is performed by matching the amplitude of onepower spectrum. For this purpose, we use the BB spectrum be-cause the main motivation of the simulations is to produce po-larized dust skies for component separation of B -modes. Thus,to Eqs. 6 we add the fourth constraint N B ˆ = ( p f ) , (7)where N B is an overall factor that scales the BB power spectrumof stage a maps divided by p to the desired amplitude. The fourparameters t , p , ρ, and f can be derived analytically from thefour input parameters R TT , R TE , R BB , and N B . One can chooseany values for R TT , R TE , R BB , and N B , as long as the normaliza-tion is positive and the ratios respect the condition R TT > R TE forced by the positive definiteness of the power spectra covari-ance.We construct the b T (cid:96) m , b E (cid:96) m , and b B (cid:96) m according to their defi-nitions in Eqs. 5. The final product is a triplet of Stokes maps( I b , Q b , U b ) that have the desired two-point statistics.
4. Simulated maps
To illustrate our method, we apply the formalism presented in theprevious sections and simulate dust polarization maps that fit the
Planck power spectra. The input values for R TT , R TE , and R BB inEqs. 6 are derived from Planck data (Sect. 4.1). We introduce thesimulated maps in Sect. 4.2. The method used to compute thesemaps is detailed in Appendix A.
The EE , BB , T E , and T B angular power spectra of dust polariza-tion were measured using the
Planck maps at 353 GHz on the sixlarge regions at high and intermediate Galactic latitude definedin PXXX. The e ff ective sky fraction f sky , after a 5 ◦ (FWHM)apodization, ranges from f sky =
24% to f sky = EE and BB spectra reported in PXXX are well fittedby power laws with exponents α dataEE , BB = − . ± .
02, with nosystematic dependence on the sky region. The amplitudes of the
Table 1.
Input values for the simulations. f sky R TT R TE R BB α dataBB A BB , data µ K
33% 44.2 ± ± ± ± ± Table 2.
Values of the parameters t , p , ρ, and f corresponding tothe ratio and normalization values of Table 1 and to our fiducialset of values for N , f M , and α M . t p ρ f ± ± ± ± spectra at a reference multipole (cid:96) = A EE , data , were mea-sured from power-law fits over the range 40 < (cid:96) <
600 with anindex fixed to its mean value of − .
42. These amplitudes are ob-served to increase with the mean total dust intensity in the mask, I dust , following the law A XX , data ∝ ( I dust ) . ± . ( X = E , B ). Wecombine the amplitude A EE , data and the EE to BB ratios listed inTable 1 of PXXX for their LR33 mask to compute the amplitude A BB , data of the D BB , data (cid:96) spectrum at (cid:96) = R TT and R TE ratios are not listed in Table 1of PXXX. To determine these values, we combine the fit to the EE spectrum from PXXX, the TE spectrum plotted in Fig. B.1 ofPXXX, and the TT spectrum we computed using the Planck dustmap at 353 GHz obtained by Planck Collaboration Int. XLVIII(2016) after separation from the cosmic infrared background(CIB) anisotropies. The C (cid:96) data points and error bars of the TEspectrum were provided to us by the contact author of PXXX.The spectra are binned between (cid:96) =
40 and (cid:96) =
600 with ∆ (cid:96) =
20 and the binned spectra are noted C XY , datab ( XY = T T , T E ). Wecompute the ratios R XY by comparing the measured power spec-tra C XY , datab with the power-law fit to the EE spectrum C EE , datab ,minimizing the following chi-squared: χ ( R ) = (cid:88) b (cid:16) C XY , datab − R C EE , datab (cid:17) / ( σ XY , datab ) , (8)where σ XY , datab is the standard deviation error on C XY , datab outputfrom the Xpol power spectrum estimator .The values we use as input for the simulations are gatheredin Table 1. Here and in Appendix A, we introduce the simulated maps anddescribe how we produce them.We have analyzed the simulated maps over a larger skyarea than in PXLIV. We have not, however, attempted to fit thePXLIV model of the mean field to the
Planck data over a largerregion. In particular, the adopted mean field direction is given bythe same Galactic coordinates ( l , b ) = (70 ◦ , ◦ ). Although thisspecific choice a ff ects the Q a and U a maps, it has no critical im-pact on the statistical results presented in the paper. Our fiducialset of values for N , f M , and α M is 4, 0.9, and − .
5, respectively.To quantify the impact of these parameters on the model powerspectra, we computed simulated maps for several combinationsaround the fiducial values within the constraints set by PXLIV.For N we considered two values 4 and 7, and for f M the range Xpol is an algorithm for power spectrum estimation that is an ex-tension to polarization of the Xspect method (Tristram et al. 2005)4. Vansyngel et al.: Simulations of dust foreground to CMB polarization α M from − . − . a and b mapsis described in Sects. A.1 and A.2. We produced our simulationsat an angular resolution of 30 (cid:48) on a HEALPix (G´orski et al.2005) grid with resolution parameter N side = p was computed at stage b , we needed an ini-tial guess in order to compute the total intensity map of stage a maps (see Eqs. 2 for I ( ν )). Based on PXLIV, we took p = . b maps from Q a / p and U a / p that do not depend much on p (Sect. A.1 ).The parameters t , p , ρ, and f used to construct stage b mapswere determined by the ratios R TT , R TE , and R BB and the ampli-tude of the BB spectrum (Sect. A.2). We used the BB amplitudeand the ratio values computed on the LR33 mask (Table 1). Thecorresponding values of the stage b parameters t , p , ρ and f arelisted in Table 2 for our fiducial set of values for N , f M and α M .The value of p , 0 . ± .
05 agrees with that derived by PXLIVfrom 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 spaceinduces leakage of power from the Galactic plane to high lat-itudes. In order to avoid this artifact, the brightest part of theGalactic plane must be masked before performing the transfor-mation. The Planck collaboration provides eight Galactic masksfor general purposes. They are derived from the 353 GHz in-tensity map by gradually thresholding the intensity after havingsubtracted the CMB. These masks are then apodized with a 2degree Gaussian kernel and cover respectively 15, 33, 51, 62,72, 81, 91, and 95% of the sky . The precise choice of the maskis not critical. We chose the mask corresponding to f sky = E - B asymmetry downto very low multipoles changes the one-point statistics of frac-tion and angle of polarization. To prevent this artefact, we intro-duce the E - B asymmetry and the T E correlation smoothly fromlow multipoles. In practice, the parameters ρ and f are functionsof (cid:96) as follows: (cid:40) ρ ( (cid:96) ) = ρ w ( (cid:96) ) f ( (cid:96) ) = − (1 − f ) w ( (cid:96) ) . (9)Here w ( (cid:96) ) is a window function going smoothly from 0 to 1around multipole (cid:96) c and is defined as follows: w ( (cid:96) ) = (cid:96) ≤ (cid:96) c − δ(cid:96)/ (cid:16) − sin (cid:16) (cid:96) c − (cid:96)δ(cid:96) π (cid:17)(cid:17) / (cid:96) c − δ(cid:96)/ < (cid:96) < (cid:96) c + δ(cid:96)/
21 if (cid:96) c + δ(cid:96)/ (cid:54) (cid:96) , (10)where we set (cid:96) c =
30 and δ(cid:96) =
30. After this modification, the E - B power ratio tends to 1 for (cid:96) < (cid:96) c in agreement with the EE and BB Planck
353 GHz power spectra presented in Fig. 20 http://healpix.sourceforge.net These masks are available on the Planck Legacy Archive as
HFI Mask GalPlane-apo2 2048 R2.00.fits and described in thePlanck Explanatory Supplement 2015 accessible at the web page https://wiki.cosmos.esa.int/planckpla2015/index.php/Frequency_Maps − −
50 0 50 1000.0000.0050.0100.0150.0200.025 − −
50 0 50 1000.0000.0050.0100.0150.0200.025 P D F P D F p ( ⇥ ) (deg) Fig. 1.
Probability distribution functions of p and ψ (top andbottom plots) for stage a and stage b maps (red and blue his-tograms). The maps were computed using the fiducial values of α M , f M , and N and the corresponding parameters t , p , ρ, and f introduced in Sect. 4.2. The distributions are computed on thesouthern Galactic polar cap ( b ≤ − ◦ ) as in PXLIV. The veryclose match between the corresponding histograms shows thatthe inclusion of the T E correlation and the E - B asymmetry doesnot alter the one-point statistics of the simulated maps.of Planck Collaboration Int. XLVI (2016) at (cid:96) <
30. Figure 1shows that the distributions (one-point statistics) of p and ψ com-puted around the southern Galactic pole of the stage a and b maps are very similar.
5. Model power spectra
In this section, we show that our simulated stage b maps repro-duce the Planck EE , BB , and T E dust spectra constraining theexponent α M of the magnetic field power spectrum (Sect. 5.1),provide the statistical variance of the dust polarization power ina given (cid:96) bin (Sect. 5.2), and match the observed scaling betweenthe spectra amplitude and the mean dust total intensity for bothlarge and small sky regions (Sect. 5.3). To compare our model results directly with the analysis of the
Planck data in PXXX, we compute power spectra of the sim-ulated maps over the LR33 mask. The power spectra are com-puted using the PolSpice estimator (Chon et al. 2004) that cor-rects for multipole-to-multipole coupling due to the masking.We checked that we obtain very similar results when the spectraare computed with the Xpol estimator.
5. Vansyngel et al.: Simulations of dust foreground to CMB polarization − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − f M ↵ M ↵ EE ↵ BB p Fig. 2.
Parameter p (top) and slopes of the EE (middle) and BB (bottom) power spectra of stage b maps vs. the slope α M of thepower 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 = . α M = − .
5. Red stars (blue squares) represent results for N = N = α 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 (a dashed horizontal line for the mean, and a red 45 ◦ (resp. blue -45 ◦ ) hatched regionfor 1 σ uncertainty for N = N = f M and α M as follows. First, we keep f M fixed to 0.9 and let α M vary from − . − . − . α M to − . f M vary from 0.7 to 1 insteps of 0.1. For each set of parameters, we compute a sampleof 1000 realizations with the procedure described in Sect. 4.2and Appendix A. The power spectra of stage b maps are binnedfrom (cid:96) =
60 to 200 with a bin width of ∆ (cid:96) =
20. We fit themodel A XX ( (cid:96)/(cid:96) ) α XX + ( X = E or B , (cid:96) =
80) to the sample meanspectrum D XX (cid:96) . The weights used in the fit are the entries of thesample covariance matrix. For each pair of ( f M , α M ) values, wecan derive the mean and covariance of ( A XX , α XX ) from the fit.Fig. 2 shows the changes in the parameter p and the spectralindices α EE and α BB when varying either f M or α M , for N = p , α EE , and α BB and the error bars represent the sample standarddeviation. The results are compared to the data values reportedin PXLIV and in PXXX. In PXLIV, the authors constrain thevalue of p with one-point statistics of the p and ψ around thesouth pole at a fixed number of layers N (see middle plot ofFig. 10 of PXLIV). Over the range of values we consider, p , α EE , and α BB are mostly sensitive to α M . The comparison of thepower spectra between simulations and data does not constraint f M nor N . The parameter f M a ff ects both the dispersion of ψ and p through depolarization along the line of sight (PXLIV). Thesetwo e ff ects modify the variance of the dust polarization in oppo-site directions. The fact that the parameter p is independent of f M (see top right panel of Fig. 2) suggests that they compensateeach other over the range of values we are considering.The measured values of α EE and α BB constrain α M to be − . B t spectra ( α M ≤ − . α EE and α BB are roughly constant with mean values lower than the ob-served values. In this regime, turbulence is not significant overthe (cid:96) range used in this analysis. The dust total intensity mapand the changing orientation with respect to the line of sight ofthe mean magnetic field dominate the variance of the polarizedmaps. For α M ≥ − . α EE and α BB are roughly equal to α M within a small positive o ff set of about 0.1. In other words, theexponents of the dust polarization spectra reproduce the expo-nent of the magnetic field power spectrum.The parameter p may also be used to constrain α M . If the p values from PXLIV for N =
6. Vansyngel et al.: Simulations of dust foreground to CMB polarization
Multipole ` D ` ( µ K )
70 100 20010100 70 100 20010100
Fig. 3. EE (top curve, diamond symbols) and BB (bottom curve,star symbols) power spectra of the simulated maps and their fitsfor the LR33 sky region. The diamonds and the stars representthe mean value computed over 1000 realizations. The 1 σ errorbars are derived from the sample standard deviation of the powerin each (cid:96) bin. The blue dashed lines represent the fits to the meanspectra and the blue dotted lines the 1 σ error on the fits. The redshade areas represent the power-law fit and the 1 σ errors to the Planck data reported in PIPXXX for the LR33 region.we find that the model fit constrains α M to be − . p with α M follows from dispersion of the B t ori-entation on angular scales corresponding to multipoles (cid:96) > f M , this dispersion decreases as the power spectrumof B t steepens (i.e., toward low values of α M ). Hence, the ob-served amplitude of the BB spectrum is matched for increasingvalues of p when α M decreases.Fig. 3 shows the EE and BB power spectra for our fiducialvalues of α M , f M , and N . The points represent the mean valuecomputed over 1000 realizations. The errors are derived fromthe sample variance of the power in each (cid:96) bin. The fit from theanalysis of PXXX and its 1 σ error are overplotted. The simula-tions are able to reproduce the EE and BB dust power spectra.The asymmetry parameter f has a value smaller than unity. Thefactor f = .
55 is close to the value of R BB = .
48 (Table 1).Within this model, unlike for that of Ghosh et al. (2016), the
T E correlation accounts for only a small part of the E - B asymmetry.Fig. 4 shows ratios between the di ff erent power spectra ofthe simulated maps. Each point represents the sample mean ofthe 1000 ratios D TT (cid:96) / D EE (cid:96) , D TE (cid:96) / D EE (cid:96) and D BB (cid:96) / D EE (cid:96) of each binand the error bars represent the sample standard deviation. Forcomparison, we plot the input values and uncertainties of the R TT , R TE , and R BB ratios. The ratios computed on the simulatedmaps are consistent with the input values, as expected becausethe maps were constructed in such a way that their power spectrarespect that covariance structure. Our simulations allow us to compute the dispersion of the dust BB power within a given (cid:96) bin. Although the dust maps are com-puted from Gaussian realizations of the turbulent field, the var-ious processes involved in the computation might make themnon-Gaussian. For example, we do not expect the distribution ofthe power at multipole (cid:96) to tend to a Gaussian distribution for (cid:96) → ∞ as quickly as it would for a Gaussian random field. For
70 100 200Multipole l R BB R T E R TT Multipole ` R TT R T E R BB Fig. 4.
Three ratios R TT , R TE , and R BB computed on the simu-lated maps for the LR33 sky region. The blue diamonds are themean ratios for each (cid:96) bin computed over 1000 realizations. Forthe R TT and R TE ratios, the dashed line and the light and darkgray regions represent the input value and the 1- and 2 σ errorson the input value, respectively. For the R BB ratio, the dashedline, the light and dark gray regions represent, respectively, theratio between the fits of the BB and EE data spectra from PXXXand their 1- and 2 σ uncertainties.the same reason, the variance of the distribution of the power fora given (cid:96) is not necessarily the cosmic variance.Fig. 5 shows the distribution of the power within one mul-tipole bin around (cid:96) =
110 with a bin width of ∆ (cid:96) =
20. Thepower spectra were computed for the LR33 region for which thecovered sky is roughly equally distributed around the north andsouth Galactic poles. The figure also presents a Gaussian fit tothe histogram and the expected cosmic variance for the samebin if the maps were drawn from a Gaussian random field onthe sphere. The actual dispersion is a few times larger than thecosmic standard deviation. This e ff ect might be due to the non-stationarity of the intensity map. The LR33 region includes somebright structures in dust total intensity. These localized structuresare likely to be the explanation for the enhanced dispersion in thesimulations. If this is the right interpretation, the enhancementmust apply to the true sky because we are using the Planck totaldust intensity map in our model.In addition to the spread, we looked at the shape of the PDFof the power per bin. We made 10000 of our dust simulations and10000 Gaussian random simulations. The power spectrum usedto produce the Gaussian realizations is the sample mean powerspectrum of the 10000 dust simulations. For the two cases, wecomputed the power spectra, binned them with a width of ∆ (cid:96) = , and fitted a Gaussian function to the sample distribution. In
7. Vansyngel et al.: Simulations of dust foreground to CMB polarization c o un t C BB ` Fig. 5.
Distribution of the power per bin computed on simulatedmaps is significantly broader than the cosmic variance. The solidred line represents a Gaussian fit to the distribution in the multi-pole bin (cid:96) =
110 with a width ∆ (cid:96) =
20 (black histogram). Thedashed line represents the distribution expected for a Gaussianrandom field in that same bin.both cases, the dust simulations and the Gaussian realizations,we see the same di ff erence between the PDF of the power per binand the Gaussian fit. We concluded that the shape of distributionof the power per bin of our simulations is very similar to that ofa Gaussian random field. We now show that the simulations reproduce the
Planck powerspectra for the high latitude sky in general, not just for the spe-cific sky region LR33 used as input. First, we compute the spec-tra of the simulated maps for the five other LRxx sky regionsfrom PXXX. Second, as in PXXX, we compute the spectra forsmaller sky patches at high Galactic latitude with f sky = Planck results.The analysis of the simulations on the six regions providessix sample mean power spectra and their sample variances. Thepower spectra on each region are computed and are fitted in thesame way as described in Sect. 5.1. In Table 3, we gather theresults of the fits together with the corresponding
Planck valuescollected from Table 1 of PXXX for comparison. Error bars onthe data measurements are smaller than those of these noiselesssimulations because the error bars on the simulation spectra con-tain the variance from multiple random realizations of the GMFthat does not a ff ect the data.While the simulations are constructed such that they matchthe data on one particular sky region (LR33), Table 3 shows thatthey also agree with the data on the other five regions within asmall di ff erence, which we comment on below. The spectra am-plitudes at (cid:96) =
80 increase with the sky fraction faster than whatPXXX reported for the
Planck data. This slight di ff erence mayarise from the fact that we assumed a fixed value of N indepen-dent of the dust total intensity and Galactic latitude. In models ofstellar polarization data at low Galactic latitudes and in molec-ular clouds, Jones et al. (1992); Myers & Goodman (1991) as-sumed that N scales linearly with the dust column density. Whiletheir model hypothesis would not work for the di ff use ISM, wecould consider variations in N . Alternatively, the slight di ff er-ence in scaling could come from another simplifying assumption
100 100010
100 100010 EE BB I ( µ K CMB ) A EE ( µ K C M B ) A BB ( µ K C M B ) Fig. 6.
Amplitudes of the power spectra ( A EE and A BB ) plottedvs. the mean total dust intensity at 353 GHz computed on each1% sky region for 100 realizations. Each vertical black line rep-resents the sample mean and sample standard deviation of oneof the 400 deg patch. The blue dashed and dotted lines repre-sent the power-law fit with the 1 σ uncertainty of the simulationresults. For comparison the red line is the same fit to the Planck data for the same set of sky patches.of the method, as we ignore the variation of the mean GMF ori-entation with distance from the Sun. It will be possible to modifyour model to test these two ideas but this is beyond the scope ofthe present paper.For the analysis on the 1% sky patches, we perform the sameprocedure as in PXXX to derive the empirical law between theamplitude at (cid:96) =
80 of the power spectra and the total intensity.Fig. 6 shows the amplitudes of the EE and BB spectra as a func-tion of the mean intensity of each patch. We realized 100 simula-tions and each vertical black line represents the sample mean andsample dispersion amplitude of one 400 deg patch. The empiri-cal law derived from a linear fit in the log( I ) − log( A XX ) spaceis overplotted. From this fit, we find a slope value of 2.15 ± EE spectrum and of 2.09 ± BB spectrum. Thevalues of the slopes are slightly larger than 1.9 ± Planck data. This di ff erencefor the patches is similar to that observed for the large sky re-gions, 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 andtheir measurement uncertainties were not large enough to ac-count for the dispersion around the fit of Fig. 6. For our sim-ulations, the spread of the distribution of the power in a given (cid:96) bin shown in Fig. 5 can explain the scatter observed aroundthe 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 mainlyfrom the turbulent component of the magnetic field. In particu-lar, we checked that the spread around the line fit is correlatedwith the mean polarization fraction, which depends on the meanorientation of the magnetic field over a given sky patch.
8. Vansyngel et al.: Simulations of dust foreground to CMB polarization f e ff sky A EE ( µ K ) 30.6 ± ± ± ± ±
12 419 ± α EE -2.37 ± ± ± ± ± ± r A α, E -0.85 -0.84 -0.83 -0.83 -0.82 -0.81 χ ( N dof =
5) 0.1 0.7 0.4 0.8 1.4 1.0 A BB ( µ K ) 14.3 ± ± ± ± ± ± α BB -2.35 ± ± ± ± ± ± r A α, B -0.82 -0.82 -0.81 -0.81 -0.81 -0.81 χ ( N dof =
5) 0.6 0.1 0.7 0.9 0.5 1.4 A EE , data ( µ K ) 37.5 ± ± ± ± ± ± α dataEE -2.40 ± ± ± ± ± ± A BB , data ( µ K ) 18.4 ± ± ± ± ± ± α dataBB -2.29 ± ± ± ± ± ± Table 3.
Results of power-law fits to the power spectra computed on simulated dust maps for the six Galactic regions from PXXX.The quantities r A α, X and χ are the correlation between A XX and α XX and the value of the χ at the fit values, respectively. Valuesfrom the Planck data taken from PXXX are given for comparison.
6. Multifrequency simulations
So far we have discussed ways to simulate structures on thesky at a single reference frequency. Component separation meth-ods for CMB experiments rely on multifrequency data. A com-mon approach to multifrequency simulations is to simulate thesky structure and the SED separately. The SED can be simu-lated using templates or analytical forms relying on a set ofparameters, such as a modified blackbody law. The simulatedsky map at a given frequency is then extrapolated to other fre-quencies. This method could also be applied to our simulations,but it does not permit us to control the decorrelation betweenmaps at di ff erent frequencies in harmonic space, which is acharacteristic crucial for component separation as discussed inPlanck Collaboration L (2016). Indeed, the decorrelation has animpact on the relative weights between the principal foregroundmodes. Here we present a method for multifrequency simula-tions constrained to match a given set of auto- and cross-powerspectra. We follow a procedure close to that commonly used to computepseudo-random Gaussian vectors with a desired covariance fromvectors with unit covariance. We realize as many simulated dustpolarization maps as the desired number of frequencies and re-arrange them to form a new set of maps such that the covariancestructure of the latter is exactly as wanted.To build a set of N f maps at frequencies (cid:110) ν i , i = . . . N f (cid:111) , weproceed as follows:1. Simulate N f single-frequency maps obtained as described inSect. 3, whose polarization spherical harmonic coe ffi cientsare gathered in a 2 N f -dimension (E and B for N f maps) vec-tor x (cid:96) m for each pair ( (cid:96), m ).2. Compute the auto- and cross-power spectra of the maps andgather them in a matrix Σ (cid:96) , which is 2 N f × N f at each mul-tipole (cid:96) .3. Specify a covariance structure of the maps over the range of N f frequencies in the form of a 2 N f × N f matrix C (cid:96) for eachmultipole (cid:96) .4. For each multipole, (cid:96) , compute the Cholesky decompositionof Σ (cid:96) and C (cid:96) , i.e., Σ (cid:96) = L (cid:96) L † (cid:96) , (11) C (cid:96) = M (cid:96) M † (cid:96) , (12) − − − c o un t ↵ j ⌫ / ↵ ⌫ Fig. 7.
Distribution of the relative di ff erence 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, thewider) maps.where the superscript † denotes the transposition.5. For each pair ( (cid:96) , m ), construct the 2 N f -dimension vector y (cid:96) m = M (cid:96) L − (cid:96) x (cid:96) m . (13)It can be easily verified that the set of maps whose spherical har-monics coe ffi cients are gathered in y (cid:96) m has exactly the expectedauto- and cross-spectra. We applied the procedure to produce a multifrequency set ofmaps ( I ( ν, j ) , Q ( ν, j ) , U ( ν, j )) where ν = , , , , j = . . . N p is the pixel index. For both EE and BB , the diagonal of the imposed covariance is C ν × ν(cid:96) = ( ν/ν ) β ( B ν ( T ) / B ν ( T )) C ν × ν (cid:96) , where ν =
353 GHz, T = . β = . C ν × ν (cid:96) is the power spectrum of simulations at frequency ν . TheSED-independent correlation ratio R (cid:96) = C ν × ν (cid:96) / (cid:113) C ν × ν (cid:96) C ν × ν (cid:96) between two frequencies ν and ν is set to 1 below (cid:96) =
30 andset by the following equation above (cid:96) = R (cid:96) = exp − σ (cid:34) log (cid:32) ν ν (cid:33)(cid:35) . (14)
9. Vansyngel et al.: Simulations of dust foreground to CMB polarization
This dependence applies if the variations of the SED can beparametrized with a spatially varying spectral index (AppendixB). The parameter σ is set in such a way that the correlation be-tween the 353 and 217 GHz channels is 0.9 within the range ofvalues measured on Planck data (Planck Collaboration L 2016).We then construct the SED map α j ν from α j ν = (cid:112) Q ( ν, j ) + U ( ν, j ) (cid:112) Q ( ν , j ) + U ( ν , j ) (15)and compute the mean SED α ν from α ν = (cid:89) j α j ν / N p . (16)In Fig. 7, we plot the distribution of α j ν /α ν − ν = , , , and 217 GHz. As expected, the distributionwidens with the separation between ν and ν because the corre-lation coe ffi cient R (cid:96) decreases. The correlation between the nor-malized SED of the same four frequencies is given by
217 143 100 70217 .
91 0 .
80 0 . .
91 1 0 .
94 0 . .
80 0 .
94 1 0 . .
74 0 .
86 0 .
96 1 . This matrix gives an estimation of the coherence of the nor-malized SED through frequencies. We do not control the way theSED of a given sky pixel varies with respect to the mean SEDbecause we model the decorrelation in harmonic space statisti-cally.
7. Astrophysical perspective
Our paper has so far focused on our contribution to compo-nent separation for CMB data analysis. We presented a phe-nomenological model that can be used to simulate dust polar-ization maps, which statistically match
Planck observations andare noise-free. In this section, we discuss what we learn about theGMF in the local interstellar medium from the modeling of thedust polarization power spectra. We examine our model resultsfrom this astrophysical perspective. We also compare our resultswith those of a companion paper Ghosh et al. (2016), which uses Hi data to account for the multiphase structure of the di ff useISM. In Sect. 7.1, we briefly review Planck power spectra of dustpolarization and our model fit. In Sects. 7.2 and 7.3, we discussthe power spectrum of the GMF and its correlation with matter.
We computed one simulation at an angular resolution of10 (cid:48) ( (cid:96) (cid:39) Planck data over a wider range of multipoles than in Sect. 5.The
T E , EE , and BB spectra are presented in Fig. 8. Thedata spectra are cross-spectra computed over the LR63 re-gion using the two half-mission maps at 353 GHz of Planck (Planck Collaboration I 2016; Planck Collaboration VIII 2016)after subtraction of the corresponding half-mission
SMICA
CMBmaps (Planck Collaboration IX 2016). The simulation is built forour fiducial parameters of the turbulence. The values of the fourparameters ( t , p , ρ, f ) were determined for this sky region andthis specific realization to be (1 . , . , . , . D l ( µ K )
100 1000multipole (l)10 Multipole ` D ` ( µ K ) T E EE BB Fig. 8.
High resolution power spectra on the LR63 region; sim-ulation vs. data. From top to bottom:
T E , EE , and BB powerspectra of the Planck
353 GHz, CMB-corrected maps (black),and one high resolution (N side = = (cid:48) ) realiza-tion of the model (red).The spectra in Fig. 8 are consistent with a single spectral ex-ponent over multipoles 40 ≤ (cid:96) ≤ (cid:96) > Planck spectra are dominated by the noise variance. At (cid:96) <
40, the spec-tra we computed with the publicly available maps are not reliabledue to uncorrected systematics. The EE and BB Planck
353 GHzspectra computed down to (cid:96) = (cid:96) < EE than for BB ; the E to B powerratio goes from about 2 to 1 toward low multipoles.An e ff ective distance to the emitting dust is necessary toconvert multipoles into physical scales. Over the high Galacticlatitude region LR63, we estimate the distance of the emittingdust to be in the range 100-200 pc. This estimate is constrainedby the distance to the edge of the local bubble (Lallement et al.2014) and the scale height of the dust emission, 200 pc at the so-lar Galacto-centric radius from the model of Drimmel & Spergel(2001). For the upper value of this distance range, the multipolerange 40-1000 corresponds to linear scales from 0.5 to 15 pc. Three of the model parameters we use – f M , N and p – wereconstrained in PXLIV. Within these constraints, we find that ourmodel fits the dust polarization power spectra for a spectral ex-
10. Vansyngel et al.: Simulations of dust foreground to CMB polarization ponent of the ˆ B t power spectrum α M = − . ± . − . ± .
02 of the
Planck dust that is measured overthe same range of multipoles on the EE and BB
353 GHz
Planck spectra. Thus, a main conclusion of our modeling is that the ex-ponent of the dust polarization spectra is that of the ˆ B t spectrum.The same conclusion is reached by Ghosh et al. (2016) for a dis-tinct modeling of the polarization layers. This conclusion holdswithin the common framework of these two models and the cor-responding assumptions.The spectral exponent α M we derive from the datafit is significantly larger than the Kolmogorov value of − / ff erence has been reported for the GMF spectrum de-rived over a similar range of scales from the analysis of syn-chrotron emission (e.g., Iacobelli et al. 2013) and of Faradayrotation measures (Oppermann et al. 2012). As discussed the-oretically for synchrotron emission by Chepurnov (1998) andCho & Lazarian (2002), a shallower slope is expected for (cid:96) mul-tipoles approaching π L max L out , where L max is the length of the emit-ting layer along the line of sight and L out the outer scale of tur-bulence. Two given lines of sight cross independent turbulentcells when their separation angle approaches the angle ∼ L out L max .It is only for smaller separation angles that the power spectrumof the emission reflects that of the magnetic field. This explana-tion put forward for synchrotron emission and Faraday rotationin earlier studies could apply to our analysis of dust polarizationtoo. The flattening observed at (cid:96) <
20 in the spectra presented byPlanck Collaboration Int. XLVI (2016) supports this interpreta-tion, but the
Planck data do not have the sensitivity to fully testit by checking whether the dust polarization spectra steepen at (cid:96) > − . E map they computed assuming a perfect alignment betweenthe magnetic field and filamentary structure of their cold neutralmedium Hi map. In this section we relate the structure of the GMF to that of thegas density in the di ff use ISM. The two are expected to be cor-related to the extent that the magnetic field is frozen in mat-ter. We note that this assumption might not hold everywhere(Eyink et al. 2013). The dust total intensity at 353 GHz is a tracerof interstellar matter within some limitations characterized ina number of studies (e.g., Planck Collaboration Int. XVII 2014;Planck Collaboration XI 2014), which are not a main concernfor this discussion. The spectrum of the GMF we find is close tothat measured for the dust total intensity. Over the same (cid:96) range,Planck Collaboration Int. XLVIII (2016) report an exponent of − . T T spectrum of their 353 GHz map corrected forCIB anisotropies, and Ghosh et al. (2016) report a value of -2.6for their total dust intensity map built from Hi data.Dust polarization data have been used to quantify thealignment of the magnetic field orientation with the fil-amentary structure of the di ff use ISM (Clark et al. 2014;Planck Collaboration Int. XXXII 2016). This is a striking facetof the correlation between matter and the GMF, which creates T E correlation and thereby E - B power asymmetry (Clark et al. 2015; Planck Collaboration Int. XXXVIII 2016). Ghosh et al.(2016) presented a model of dust polarization where this corre-lation between matter and the GMF applies to one single polar-ization layer that is associated with the cold neutral medium astraced by narrow Hi spectral lines. In their model that layer ac-counts for both the T E correlation and the E - B asymmetry mea-sured over the sky region with the lowest dust column density inthe southern sky they analyzed. In our model, the T E correlationis introduced by adding one dust emission layer, where polariza-tion is only in E -modes and is fully correlated to the T map. Thiscorresponds to the additive term proportional to the ρ parameterin the second equation in Eqs. 5. The dust filamentary structuresare present in all layers and the polarization results from the ad-dition of the signals. We checked on the simulated images thatthis process introduces a preferred alignment between the fila-mentary structure of the T map and the magnetic field orientationinferred from the polarization angle, but this alignment is notas tight as that reported by Planck Collaboration Int. XXXVIII(2016) from their analysis of the most conspicuous filaments athigh galactic latitudes in the Planck data. This di ff erence comesfrom the fact that we use the same intensity map for each layer.Planck Collaboration Int. XXXVIII (2016) shows that the fila-mentary structure of the cold neutral medium has a main contri-bution to the E - B asymmetry but it does not exclude a signifi-cant contribution related to the generic anisotropy of MHD tur-bulence, as suggested by Caldwell et al. (2016). We stress herethat our modeling of the E - B power asymmetry is mathemati-cal. It does not constrain its physical origin. In this respect, ourmodel is a framework that we are using to match the data statis-tically, 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 separationmethods in CMB data analysis. We detailed the simulation ofdust polarization maps at one frequency before we introduced amathematical means to produce maps at several frequencies andmatched a given set of auto- and cross-spectra. Our method andthe main results obtained by analyzing the simulated maps aresummarized here.Our approach builds on earlier studies, i.e., the analysis of
Planck dust polarization data and the model framework fromPXLIV, which relate the dust polarization sky to the structure ofthe GMF and interstellar matter. The structure of interstellar mat-ter is the dust total intensity map from
Planck . The GMF is mod-eled as a superposition of a mean uniform field and a Gaussianrandom (turbulent) component with a power-law power spec-trum of exponent α M . The integration along the line of sight per-formed to compute the Stokes maps is approximated by a sumover a small number of emitting layers with di ff erent realiza-tions of the random GMF component. The mean field orienta-tion, the amplitude of the random GMF component with respectto the mean component, the spectral exponent α M , and the num-ber of polarization layers are parameters common to the modelfrom PXLIV. To match the power spectra of dust polarizationmeasured with the Planck data, we add two main parameters ( ρ and f ) that introduce mathematically the T E correlation and E - B power asymmetry. They are determined by fitting the Planck353 GHz power spectra for (cid:96) >
40 on one sky region at highGalactic latitude, LR33 from PXXX.The model allows us to compute multiple realizations of theStokes Q and U maps for di ff erent realizations of the randomcomponent of the magnetic field and to quantify the dispersion
11. Vansyngel et al.: Simulations of dust foreground to CMB polarization of dust polarization spectra for any given sky area away fromthe Galactic plane. The simulations reproduce the scaling lawsbetween the dust polarization power and the mean total dust in-tensity from
Planck , including the observed dispersion aroundthe mean relation.This paper discusses what we learn about the GMF in thelocal interstellar medium from the modeling of the dust polar-ization power spectra. We find that the slopes of the EE and BB power spectra of dust polarization measured by Planck arematched for α M = − . ± .
1. As in Ghosh et al. (2016), wefind that, for our model, the exponent of the spectrum of ˆ B t isvery close to that of the dust polarization spectra. This exponentis larger than the Kolmogorov value of − / − . (cid:96) = − Planck dust total intensityat 353 GHz as a tracer. Our model does not allow us to commenton the origin of the
T E correlation and E - B asymmetry.It would be possible to extend the model we presented in sev-eral ways, which might lead to fruitful explorations. To fit dustpolarization spectra down to the very low multipoles relevant formeasuring E and B -mode CMB polarization associated with theuniverse reionization, we might need to account for the injectionscale of turbulence. Phenomenologically, this could be carriedout by introducing a low- (cid:96) cuto ff in the power spectrum of themagnetic field in Eq. 4.Further model changes could also provide a better match tothe data, in particular toward low Galactic latitudes. We haveused a constant orientation for the mean GMF. A 3D model ofthe density structure of the Galactic ISM can be used to assigndistances to the shells, and, thereby, to take into account the 3Dstructure of the large-scale magnetic field, as in, for example,Fauvet et al. (2011) and Planck Collaboration Int. XLII (2016).In this case the intensity maps will di ff er for each layer andthe e ff ective number of layers could be allowed to vary with,for example, Galactic latitude or dust column density. In sucha model, it would be possible to introduce, for each layer, thecorrelation between matter and the GMF and distinct dust SEDs.This method of introducing the decorrelation of dust polarizationmaps with frequency might in essence better represent the line-of-sight averaging of polarization data (Tassis & Pavlidou 2015;Planck Collaboration L 2016) than the mathematical means pro-posed here. Finally, our paper focuses on dust polarizationbut a similar approach could be applied to produce maps ofsynchrotron polarization that match the observed correlationwith dust polarization (Planck Collaboration Int. XXII 2015;Choi & Page 2015). Acknowledgements. The research leading to these results has received fund-ing from the European Research Council under the European Union’s SeventhFramework Programme (FP7 / / ERC grant agreement No. 267934.
References
Abazajian, K. N., Arnold, K., Austermann, J., et al. 2015, Astroparticle Physics,63, 55Andersson, B.-G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209BICEP2 / Keck Array and Planck Collaborations. 2015, Phys. Rev. Lett., 114,101301Brandenburg, A. & Lazarian, A. 2013, Space Sci. Rev., 178, 163Caldwell, R. R., Hirata, C., & Kamionkowski, M. 2016, ArXiv e-printsChepurnov, A. & Lazarian, A. 2010, ApJ, 710, 853Chepurnov, A. V. 1998, Astronomical and Astrophysical Transactions, 17, 281Cho, J. & Lazarian, A. 2002, ApJ, 575, L63Choi, S. K. & Page, L. A. 2015, ArXiv e-printsChon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350,914 Clark, S. E., Hill, J. C., Peek, J. E. G., Putman, M. E., & Babler, B. L. 2015,Physical Review Letters, 115, 241302Clark, S. E., Peek, J. E. G., & Putman, M. E. 2014, ApJ, 789, 82Delabrouille, J., Betoule, M., Melin, J.-B., et al. 2013, A&A, 553, A96Dolginov, A. Z. & Mitrofanov, I. G. 1976, Ap&SS, 43, 291Drimmel, R. & Spergel, D. N. 2001, ApJ, 556, 181Dunkley, J., Amblard, A., Baccigalupi, C., et al. 2009, in American Instituteof Physics Conference Series, Vol. 1141, American Institute of PhysicsConference Series, ed. S. Dodelson, D. Baumann, A. Cooray, J. Dunkley,A. Fraisse, M. G. Jackson, A. Kogut, L. Krauss, M. Zaldarriaga, & K. Smith,222–264Eyink, G., Vishniac, E., Lalescu, C., et al. 2013, Nature, 497, 466Fauvet, L., Mac´ıas-P´erez, J. F., Aumont, J., et al. 2011, A&A, 526, A145Ghosh, T., Boulanger, F., Martin, P. G., et al. 2016, Submitted to A&AG´orski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759Guth, A. H. 1981, Phys. Rev. D, 23, 347Hoang, T. & Lazarian, A. 2016, ArXiv e-printsHoude, M., Vaillancourt, J. E., Hildebrand, R. H., Chitsazzadeh, S., & Kirby, L.2009, ApJ, 706, 1504Iacobelli, M., Haverkorn, M., Orr´u, E., et al. 2013, A&A, 558, A72Jansson, R. & Farrar, G. R. 2012, ApJ, 757, 14Jones, T. J., Klebe, D., & Dickey, J. M. 1992, ApJ, 389, 602Kalberla, P. M. W., Kerp, J., Haud, U., et al. 2016, ArXiv e-printsKamionkowski, M. & Kovetz, E. D. 2015, ArXiv e-printsLallement, R., Vergely, J.-L., Valette, B., et al. 2014, A&A, 561, A91Lazarian, A. & Pogosyan, D. 2012, ApJ, 747, 5Lee, H. M. & Draine, B. T. 1985, ApJ, 290, 211Linde, A. D. 1982, Physics Letters B, 108, 389Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153Miville-Deschˆenes, M.-A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093Myers, P. C. & Goodman, A. A. 1991, ApJ, 373, 509O’Dea, D. T., Clark, C. N., Contaldi, C. R., & MacTavish, C. J. 2012, MNRAS,419, 1795Oppermann, N., Junklewitz, H., Robbers, G., et al. 2012, A&A, 542, A93Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335Planck Collaboration XI. 2014, A&A, 571, A11Planck Collaboration I. 2016, A&A, 594, A1Planck Collaboration VIII. 2016, A&A, 594, A8Planck Collaboration IX. 2016, A&A, 594, A9Planck Collaboration X. 2016, A&A, 594, A10Planck Collaboration XII. 2016, A&A, 594, A12Planck Collaboration Int. XVII. 2014, A&A, 566, A55Planck Collaboration Int. XIX. 2015, A&A, 576, A104Planck Collaboration Int. XX. 2015, A&A, 576, A105Planck Collaboration Int. XXI. 2015, A&A, 576, A106Planck Collaboration Int. XXII. 2015, A&A, submitted, 576, A107Planck Collaboration Int. XXX. 2016, A&A, 586, A133Planck Collaboration Int. XXXII. 2016, A&A, 586, A135Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141Planck Collaboration Int. XLII. 2016, A&A, submittedPlanck Collaboration Int. XLIV. 2016, A&A, submittedPlanck Collaboration Int. XLVI. 2016, A&A, submittedPlanck Collaboration Int. XLVIII. 2016, A&A, in pressPlanck Collaboration L. 2016, A&A, submittedSeljak, U. & Zaldarriaga, M. 1996, ApJ, 469, 437Starobinskiˇi, A. A. 1979, Soviet Journal of Experimental and Theoretical PhysicsLetters, 30, 682Stein, W. 1966, ApJ, 144, 318Tassis, K. & Pavlidou, V. 2015, MNRAS, 451, L90Tristram, M., Macias-Perez, J. F., Renault, C., & Santos, D. 2005, Mon. Not.Roy. Astron. Soc., 358, 833Waelkens, A., Ja ff e, T., Reinecke, M., Kitaura, F. S., & Enßlin, T. A. 2009, A&A,495, 697Yang, M. 2008, Applied Economics Letters, 15, 737Zaldarriaga, M. 2001, Phys. Rev. D, 64, 103001 Appendix A: Implementation of the method
In this appendix, we explain how we compute the dust polar-ization maps used in this paper. Sect. A.1 presents the proce-dure we use to derive the stage a maps using the framework inSect. 2. Sect. A.2 describes how we produce the stage b mapsthat match the dust T E correlation and E - B asymmetry mea-sured by Planck , using the method described in Sect. 3.
12. Vansyngel et al.: Simulations of dust foreground to CMB polarization
A.1. Stage a maps We explain how we produce the ( I a , Q a , U a ) maps at a referencefrequency ν , which we choose to be 353 GHz, the best-suited Planck channel to study dust polarization. These maps have no
T E correlation and no E - B asymmetry at (cid:96) > I a is not computed from Eqs. 2 but derivedfrom observations. We use I a = D , where D is the dust totalintensity map at 353 GHz of Planck Collaboration Int. XLVIII(2016) after separation from the CIB and CMB anisotropies. Tocompute 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 real-ization for each of the three components of ˆ B t in Eq. 3. Theangle maps are computed for the total magnetic field B includ-ing the mean magnetic field B . With the set of angles maps γ i ,using the Stokes I equation in Eqs. 2, we compute the map S i ( ν )at the frequency ν , S i ( ν ) = D / N (cid:88) i = (cid:34) − p (cid:32) cos γ i − (cid:33)(cid:35) , (A.1)where S i ( ν ) has been assumed to be independent of the index i and p is set to a fiducial value of 0.25. Next, we combine S i ( ν )and the angle maps γ i and ψ i in the Stokes Q and U equations inEqs. 2 to compute the ratio maps Q a / ( p × I a ) and U a / ( p × I a )at the frequency ν . These ratio maps are independent of I a and depend on p only through S i ( ν ). They are computed atpixel resolution defined by the N side =
256 HEALPix parameter.After multiplication by D , we obtain the two maps Q a / p and U a / p , which have an ill-defined beam transfer function. The D map has a resolution that varies across the sky. We over-come this issue by smoothing ( I a , Q a , U a ) to a resolution lowerthan the lowest resolution of the D map. The model maps usedin the paper have N side =
256 and a symmetric Gaussian beamwith a full width at half maximum of 30 (cid:48) . A.2. Stage b maps From the harmonic coe ffi cients of Eq. 5, we compute the powerspectra of stage b maps at a given multipole (cid:96) , as functions of t , p , ρ, f , and x , where x ˆ = E (cid:104) A EE (cid:96) / A TT (cid:96) (cid:105) = E (cid:104) A BB (cid:96) / A TT (cid:96) (cid:105) , (A.2) A TT (cid:96) , A EE (cid:96) , and A BB (cid:96) are the power spectra of stage a mapsand E [ · ] is an averaging over multipoles between (cid:96) =
60 and (cid:96) = T T and polarization spectra areclose to one another, the ratio x is close to being independent ofmultipole (cid:96) . Since this simplification approximately applies fordust emission (PXXX), we consider the ratio x to be constantover the relevant multipole range.Assuming A XY (cid:96) = X (cid:44) Y , the ratios of Eq. 6 can beexpressed as follows: R TT = z + y R TE = z + y R BB = f y + y , (A.3)where y ˆ = x / ( p ρ ) and z ˆ = t / ( p ρ ). When the ratios R XY are cho-sen, then the system A.3 becomes a system of equations in { f , y , z } . Although the system is not linear, it can be inverted, aslong as R TT > R TE , i.e., Det (cid:32) D TT (cid:96) D TE (cid:96) D TE (cid:96) D EE (cid:96) (cid:33) >
0. When choosingvalues for the ratios R XY , this condition has to be satisfied be-cause the power spectra form a covariance, which must be pos-itive definite. Restricting the set of solutions to positive reals,there is a unique solution, i.e., f = (cid:113) R BB R TT / (cid:16) R TT − R TE (cid:17) y = (cid:113)(cid:16) R TT − R TE (cid:17) / R TE z = R TT / R TE . (A.4)From the solution of Eq. A.4 and the normalization factor N B = ( p f ) , we can compute the parameters ρ, f , t , and p asfollows: f = f , p = (cid:112) N B / f , ρ = x / ( p y ) , t = zp ρ . (A.5)We note that ρ and the correlation coe ffi cient between the T and E parts of stage b maps, noted r TE = R TE / R . TT , are related asfollows: ρ = xp (cid:115) r TE − r TE . (A.6)We choose the BB normalization factor N B such that thepower spectrum A BB (cid:96) of stage a divided by p map is adjustedto the fit of the power spectrum measured over the region LR33in PXXX (noted C BB , data (cid:96) ). Following the notation of PXXX, wehave (cid:96) ( (cid:96) + C BB , data (cid:96) = π A BB , data ( (cid:96)/ α dataBB + , where the val-ues of the parameters α dataBB and A BB , data are taken from Table 1.In the case where N B is (cid:96) -independent, N B is the solution of theminimization of the following chi-squared χ ( u ) = (cid:96) (cid:88) (cid:96) = (cid:96) (cid:32) A BB (cid:96) − u C BB (cid:96) (cid:33) /σ (cid:96) , (A.7)with σ (cid:96) the variance of A BB (cid:96) , estimated from Monte Carlo simu-lations and ( (cid:96) , (cid:96) ) = (60 , x . The fit also provides thestandard deviation on the normalization factor N B . A.3. Summary of the procedure
The following points sketch the procedure to produce our simu-lations:1. Draw Stokes maps Q a and U a divided by p as described inSect. A.12. Mask the Galactic plane and compute the harmonic coe ffi -cients a (cid:96) m
3. Given a mask, compute the full sky power spectra A (cid:96)
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. 66. Compute the corresponding solutions { f , y , z } of Eqs. A.47. Compute the parameters ρ, f , t , and p of Eqs. A.58. Construct the harmonic coe ffi cients b (cid:96) m according to theirdefinition of Eqs. 59. Transform the b (cid:96) m ’s to ( I b , Q b , U b )The stage b maps thus constructed feature the desired two-point statistics on the desired region of the sky. The procedurecan be applied on separate multipole bins, which then givesscale-dependent parameters.
13. Vansyngel et al.: Simulations of dust foreground to CMB polarization
Appendix B: Decorrelation due to a variablespectral index
This appendix shows how to compute the decorrelation in har-monic space between two frequency maps, when spectral di ff er-ences about a mean SED may be parametrized with a spatiallyvarying spectral index. This appendix restricts the proof to thesimple case where the map that is scaled through frequenciesand the spectral index map are correlated Gaussian white noisemaps.Let f ( n ) and δβ ( n ) be two Gaussian random fields on thesphere such that (cid:104) f ( n ) (cid:105) = (cid:104) δβ ( n ) (cid:105) = , (B.1) (cid:10) f ( n ) f ( n (cid:48) ) (cid:11) = δ ( n − n (cid:48) ) σ f , (B.2) (cid:10) δβ ( n ) δβ ( n (cid:48) ) (cid:11) = δ ( n − n (cid:48) ) σ β , (B.3) (cid:10) f ( n ) δβ ( n (cid:48) ) (cid:11) = δ ( n − n (cid:48) ) r σ f σ β . (B.4)From f ( n ) and δβ ( n ) we construct a set of maps at frequencies ν i , f i ( n ) = K i f ( n ) (cid:32) ν i ν (cid:33) δβ ( n ) , (B.5)where ν is a reference frequency and K i possibly contains themean SED and unit conversion factors.The aim is to compute the cross-spectrum C ν i × ν j (cid:96) ( (cid:96) (cid:62)
1) be-tween the di ff erent f i ( n ), i.e., C ν i × ν j (cid:96) = (cid:90) d n d n (cid:48) (cid:68) f i ( n ) f j ( n (cid:48) ) (cid:69) Y ∗ (cid:96) m ( n ) Y (cid:96) m ( n (cid:48) ) (B.6) = (cid:68) f i ( n ) f j ( n ) (cid:69) , (B.7)where the Y (cid:96) m ( n ) represent the spherical harmonics; we as-sumed that two directions of the maps are uncorrelated and thatthe maps are statistically isotropic. We can rewrite the product f i ( n ) f j ( n ) as follows: f i ( n ) f j ( n ) = κ i j (cid:32) g ( n ) exp (cid:34) δγ i j ( n ) (cid:35)(cid:33) , (B.8)where κ i j = K i K j σ f , g ( n ) = f ( n ) /σ f and δγ i j ( n ) = σ i j δβ ( n ) /σ β with σ i j = log( ν i ν j /ν ) σ β . It can be easily verified that (cid:32) g ( n ) δγ i j ( n ) (cid:33) ∼ N (cid:32)(cid:34) (cid:35) , (cid:34) r σ i j r σ i j σ i j (cid:35)(cid:33) (B.9)so that the expression in brackets of Eq. B.8 has a normal lognor-mal mixture distribution as parametrized in, for example, Yang(2008). Thus, (cid:68) f i ( n ) f j ( n ) (cid:69) = κ i j exp σ i j (cid:16) + r σ i j (cid:17) (B.10)and C ν i × ν j (cid:96) (cid:113) C ν i × ν i (cid:96) C ν j × ν j (cid:96) = exp − σ β (cid:34) log (cid:32) ν i ν j (cid:33)(cid:35) × (cid:16) + r σ i j (cid:17)(cid:113)(cid:16) + r σ ii (cid:17) (cid:16) + r σ j j (cid:17) . (B.11) We note that the correlation does not depend on the mean SED. Institut d’Astrophysique Spatiale, CNRS, Univ. Paris-Sud, Universit´eParis-Saclay, Bˆat. 121, 91405 Orsay cedex, France California Institute of Technology, Pasadena, California, U.S.A. Sorbonne Universit´es, UPMC Univ Paris 6 et CNRS, UMR 7095,Institut dAstrophysique de Paris, 98 bis bd Arago, 75014 Paris,France Sorbonne Universit´es, Institut Lagrange de Paris (ILP), 98 bisBoulevard Arago, 75014 Paris, France Laboratoire AIM, IRFU / Service d’Astrophysique - CEA / DSM -CNRS - Universit´e Paris Diderot, Bˆat. 709, CEA-Saclay, F-91191Gif-sur-Yvette Cedex, France LERMA, Observatoire de Paris, PSL Research University, CNRS,Sorbonne Universits, UPMC Univ. Paris 06, Ecole normale su-prieure, F-75005, Paris, France CITA, University of Toronto, 60 St. George St., Toronto, ON M5S3H8, Canada8