© 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, curl-like, signature in the polarization of the cosmic microwave background (CMB), referred to as primordial B-mode 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 B-mode 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¹353GHz polarization maps, we know that the primordial B-mode polarization of the CMB cannot be measured without subtracting the foreground emission, even in the faintest dust-emitting regions at high Galactic latitude (Planck Collaboration Int. XXX 2016, hereafter PXXX). The observed correlation between the B-mode 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 ground-based, balloon-borne, 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; Miville-Deschê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 signal-to-noise 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 Miville-Deschê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 (gradient-like) and B (curl-like) 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 B-mode 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 large-scale 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₀ 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₀, and turbulent (random), B_t, components, (3)where and are unit vectors in the directions of B₀ 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₀ has a fixed orientation in all layers. We ignore the structure of the GMF on galaxy-wide 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 power-law 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₀ and b₀ defining the orientation of , the number of layers N, the spectral exponent α_M, and the effective polarization fraction of the dust emission p₀. 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₀ = 70° ± 5° and b₀ = 24° ± 5° by fitting the large-scale pattern observed in the Stokes Q and U maps, and f_M = 0.9 ± 0.1, N = 7 ± 2 and p₀ = 26 ± 3% by fitting the distribution function (one-point statistics) of p², 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 E-B asymmetry, matching the Planck dust polarization power spectra in PXXX.

3. Introducing TE correlation and E-B 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 E-B asymmetry. The parameter t is a scaling factor for the intensity part and p₀ 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₀ and thus that the two ratios and in Eqs. (5) are independent of p₀.

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₀ 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₀ 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₀, ρ, 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 B-modes. 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₀ to the desired amplitude. The four parameters t, p₀, ρ, 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 two-point 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 ℓ₀ = 80, A^EE,data, were measured from power-law 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 power-law fit to the EE spectrum , minimizing the following chi-squared: (8)where is the standard deviation error on output from the Xpol power spectrum estimator².

The values we use as input for the simulations are gathered in Table 1.

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₀,b₀) = (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³ (Górski et al. 2005) grid with resolution parameter N_side = 256. Although the parameter p₀ 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.25. We used this value to compute stage 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 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₀, ρ and f are listed in Table 2 for our fiducial set of values for N, f_M and α_M. The value of p₀, 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⁴. 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 E-B asymmetry down to very low multipoles changes the one-point statistics of fraction and angle of polarization. To prevent this artefact, we introduce the E-B 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 E-B 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 (one-point statistics) of p and ψ computed 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 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 multipole-to-multipole 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(ℓ/ℓ₀^{)α_XX + 2} (X = E or B, ℓ₀ = 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₀ 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₀, α_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₀ with one-point statistics of the p² 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₀, α_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₀ 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₀ may also be used to constrain α_M. If the p₀ 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₀ 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₀ when α_M decreases.

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² = 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 E-B 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.

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 non-Gaussian. 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 non-stationarity 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.

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.

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² patch. The empirical law derived from a linear fit in the log (I₃₅₃)−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.

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 cross-power spectra.

6.1. Method

We follow a procedure close to that commonly used to compute pseudo-random 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 single-frequency 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 cross-power 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)

It can be easily verified that the set of maps whose spherical harmonics coefficients are gathered in y_ℓm has exactly the expected auto- and cross-spectra.

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 ν₀ = 353 GHz, T₀ = 19.6K, β = 1.6 (Planck Collaboration Int. XXII 2015) and is the power spectrum of simulations at frequency ν₀. The SED-independent correlation ratio between two frequencies ν₁ and ν₂ 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 ν₀ 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.

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 noise-free. 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.

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 cross-spectra computed over the LR63 region using the two half-mission maps at 353GHz of Planck (Planck Collaboration I 2016; Planck Collaboration VIII 2016) after subtraction of the corresponding half-mission 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₀,ρ,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 Galacto-centric 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₀ – 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 E-B 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 E-B 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 E-modes 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 E-B 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 E-B 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 cross-spectra. 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 power-law 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 E-B 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 E-B 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 B-mode 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 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 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 line-of-sight 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).

---

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 267934.

References

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 E-B 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 ν₀, which we choose to be 353GHz, the best-suited Planck channel to study dust polarization. These maps have no TE correlation and no E-B asymmetry at ℓ> 40.

The intensity map I_a is not computed from Eqs. (2) but derived from observations. We use I_a = D₃₅₃, where D₃₅₃ 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₀. With the set of angles maps γ_i, using the Stokes I equation in Eqs. (2), we compute the map S_i(ν) at the frequency ν₀, (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 in Eqs. (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 at pixel 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 overcome this issue by smoothing (I_a,Q_a,U_a) to a resolution lower than the lowest resolution of the D₃₅₃ 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₀,ρ,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₀f)², we can compute the parameters ρ,f,t, and p₀ 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₀ 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 chi-squared (A.7)with the variance of , estimated from Monte Carlo simulations and (ℓ₁,ℓ₂) = (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₀ 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₀ 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 two-point statistics on the desired region of the sky. The procedure can be applied on separate multipole bins, which then gives scale-dependent 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 ν₀ is a reference frequency and K_i possibly contains the mean SED and unit conversion factors.

The aim is to compute the cross-spectrum ( ) 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

All Figures