GALACTIC FOREGROUNDS: SPATIAL FLUCTUATIONS AND A PROCEDURE FOR REMOVAL

Let us illustrate a possible procedure for the removal of Galactic foregrounds from the CMB signal. The cosmic microwave signals consist of the CMB signal I^CMB and foregrounds I^F. When we correlate the microwave signal at points "1" and "2," we get

$$\begin{equation} \big\langle \big(I^{\rm CMB}_1+I^{\rm F}_1\big)\big(I^{\rm CMB}_2+I^{\rm F}_2\big)\big\rangle = \big\langle I^{\rm CMB}_1I^{\rm CMB}_2\big\rangle +\big\langle I^{\rm F}_1 I^{\rm F}_2\big\rangle, \end{equation} \tag{ 1 }$$

where we assume I^CMB and I^F are uncorrelated. Therefore, the measured angular spectrum C^measured_l is just the sum of C^CMB_l and C^F_l, and we have

$$\begin{equation} C^{\rm CMB}_l=C^{\rm measured}_l-C^{\rm F}_l. \end{equation} \tag{ 2 }$$

This means that, if we know the angular spectrum of foregrounds C^F_l, we can obtain the CMB angular spectrum C^CMB_l.

How can one obtain C^F_l? A well-tested way of doing this is to use multi-frequency measurements of the CMB + foreground emission and separate the two components using the frequency templates of the foregrounds. This approach requires many measurements at different frequencies. In addition, for some foregrounds the frequency templates may be difficult to obtain. The so-called spinning dust foreground introduced in Draine & Lazarian (1998a, 1998b) presents an example of such a difficult-to-remove foreground. We also mention that the measurements of foregrounds at different frequencies may have different spatial resolutions and the use of the maps with different resolutions may also present a problem.

In this paper we address a somewhat different problem, which at its extreme⁶ can be formulated in the following way. Imagine that we separated the foreground and the CMB signals at low resolution l_low using the traditional multi-frequency approach. Is it possible to use this information to remove the foreground contribution from C^measured_l for l>l_low? For instance, the measurements of the Wilkinson Microwave Anisotropy Probe (WMAP) provide a high accuracy measure of C^F_l over a limited range of scales. If we know C^F_l as a function of l_low at scales smaller than those measured by WMAP, then one can extrapolate C^F_l to l>l_low. These values of C^F_l can be used to filter the microwave measurements from balloon-borne experiments using the procedure given by Equation (2). Note that the balloon-borne experiments usually have higher spatial resolutions, but not enough frequency coverage to remove foregrounds using the frequency templates.

The key question is to what extent we can predict C^F_l over a range of scales that is different from the range of scales at which C^F_l was measured. The answer to this question is trivial if C^F_l is a simple power law. While the actual spectra of foregrounds are more complex, in this paper we provide both theoretical arguments and the analysis of the foreground data, which support the notion that C^F_l can be successfully extended beyond the range of l that is measured.

We should stress that the filtering above is different from the accepted techniques of foreground removal using frequency-dependent templates. The outcome of the latter procedures are maps of foreground radiation and CMB radiation. The filtering described above is of statistical nature. Its result is C^CMB_l rather than emission maps.⁷

In view of the above, it is important to determine to what extent the spatial properties of C^F_l are predictable, in particular, explore to what extent the power-law approximation is applicable. This paper provides a study of the spatial statistical properties of fluctuations of synchrotron and dust emission and relates those to the properties of the underlying turbulence. It also provides an example of the filtration procedure that we advocate.

Another important issue is to determine the reasons for the change in the power-law behavior with latitude. This is what we study below for the synchrotron and dust foreground emissions. The ultimate goal of this research is to obtain models of Galactic foregrounds that provide a good fit for the foreground spatial spectrum at arbitrary scales and arbitrary latitudes. This paper is a step toward constructing such a model.

For the spatial removal procedure to work, we should understand the spatial spectra of foregrounds. Since the spatial fluctuations of foregrounds inherit the spectra of the underlying interstellar turbulence, we summarize the spectral behavior of interstellar turbulence.

It is generally accepted that the ISM is magnetized and turbulent (see reviews by Elmegreen & Scalo 2004 and McKee & Ostriker 2007). The so-called Big Power Law in the Sky corresponding to the Kolmogorov-type turbulence was reported in Armstrong et al. (1995). Recently, this law, based on the measurements of the radio scintillations arising from electron density inhomogeneities, has been extended to larger scales using the Wisconsin Hα Mapper H_α fluctuations (Chepurnov & Lazarian 2010).

Spectra of magnetic turbulence have been studied using Faraday rotation measurements (see Haverkorn et al. 2008) and starlight polarization (Hildebrand et al. 2009). The interpretation of the results is more challenging in those cases.

Molecular and atomic spectral lines present a very promising way of studying turbulence. The interstellar lines are known to be Doppler-broadened due to turbulent motions. Obtaining spectra from Doppler-broadened lines is not trivial, however. A significant portion of the research in this direction based on the use of the so-called velocity centroids, which were the main tool for studying velocity fluctuations, has been shown to produce erroneous results for supersonic turbulence (Lazarian & Esquivel 2003; Esquivel & Lazarian 2005; Esquivel et al. 2007). At the same time, new techniques based on the theoretical description of the position–position–velocity (PPV) data cubes, namely the Velocity Channel Analysis (VCA) and the Velocity Coordinate Spectrum (VCS), have been developed (Lazarian & Pogosyan 2000, 2004, 2006, 2008), tested (see Stanimirovic & Lazarian 2001; Lazarian et al. 2002; Esquivel et al. 2003; Chepurnov & Lazarian 2009; Padoan et al. 2006, 2009), and applied to the observational data to obtain the characteristics of the velocity turbulence (see Lazarian 2009 for a review).

The VCA and VCS techniques reveal that the velocities for interstellar turbulence may be somewhat steeper than the Kolmogorov one, while the density spectra of the fluctuations may be substantially shallower than the spectrum of Kolmogorov fluctuations (see Padoan et al. 2006, 2009; Chepurnov et al. 2010). This agrees well with the numerical studies of the magnetized supersonic turbulence (Beresnyak et al. 2005; Kowal et al. 2007; Kowal & Lazarian 2010). For the subsonic turbulence, the spectrum of magnetized media attains values close to the Kolmogorov index (see Goldreich & Sridhar 1995; Lazarian & Vishniac 1999; Cho & Vishniac 2000; Müller & Biskamp 2000; Maron & Goldreich 2001; Lithwick & Goldreich 2001; Cho et al. 2002; Cho & Lazarian 2002b, 2003; Boldyrev 2006; Beresnyak & Lazarian 2006, 2009). In view of that, we believe that the Kolmogorov spectrum can be used as a proxy for the underlying spectra of the velocity and magnetic field, while a more cautious approach should be adopted for the density in highly compressible environments, e.g., molecular clouds.

Below we shall show that the spectra of turbulence derived from the analysis of the foreground are consistent with both the results of dedicated observations of turbulence as well as the theoretical expectation for MHD turbulence.

We use the 408 MHz Haslam all-sky map (Haslam et al. 1982) and a model 94 GHz dust emission map that are available on NASA's LAMBDA Web site.⁸ Both maps were reprocessed for HEALPix (Górski et al. 2005) with nside = 512 (7' resolution).

The original Haslam data were produced by merging several different data sets. "The original data were processed in both the Fourier and spatial domains to mitigate baseline striping and strong point sources" (see the Web site for details). The angular resolution of the original Haslam map is ∼1°. Galactic diffuse synchrotron emission is the dominant source of emission at 408 MHz.

The 94 GHz dust emission map is based on the work of Schlegel et al. (1998) and Finkbeiner et al. (1999). Schlegel et al. (1998) combined 100 μm maps of the Infrared Astronomical Satellite (IRAS) and the Diffuse Infrared Background Experiment (DIRBE) on board the Cosmic Background Explorer (COBE) satellite and removed the zodiacal foreground and point sources to construct a full-sky map. Finkbeiner et al. (1999) extrapolated the 100 μm emission map and 100/240 μm flux ratio maps to submillimeter and microwave wavelengths. The 94 GHz dust map that we used is identical to the two-component model 8 of Finkbeiner et al. (1999). The angular resolution of the 100 μm map is ∼6' and that of the temperature correction derived from the 100/240 μm ratio map is ∼1° (Finkbeiner et al. 1999), which corresponds to l ∼ 180°/θ° ∼ 180.

In this section, we separately study the synchrotron emission from the Galactic halo (i.e., b ≳ 30°) and the Galactic disk (i.e., |b| ⩽ 2°). Our main goal is to see if the statistics of the synchrotron emission from the halo is consistent with turbulence models. When it comes to synchrotron emission from the Galactic disk, it is not easy to separate diffuse emission and emission from discrete sources. Therefore, we do not try to study turbulence in the Galactic disk. Instead, we will try to estimate which kind of emission is dominant in the Galactic disk.

There exist several models for diffuse Galactic radio emission. Beuermann et al. (1985) showed that a two-component model, a thin disk embedded in a thick disk, can explain the observed synchrotron latitude profile. They claimed that the equivalent width of the thick disk is about several kiloparsecs and the thin disk has approximately the same equivalent width as the gas disk. They assumed that, in the direction perpendicular to the Galactic plane, the emissivity of each component follows

$$\begin{equation} \epsilon (z) = \epsilon (0) \mathrm{sech}^b(z/z_0), \end{equation} \tag{ 3 }$$

where z is the distance from the Galactic plane and (0), b, and z₀ are constants. The half-equivalent width of the disk, which is proportional to z₀, at the location of the Sun is ∼2 kpc.

Recently, several Galactic synchrotron emission models have been proposed in an effort to separate Galactic components from the WMAP polarization data (see, for example, Page et al. 2007; Sun et al. 2008; Miville-Deschênes et al. 2008; Waelkens et al. 2009). All the models mentioned above assume the existence of a thick disk component with a scale height equal to 1 kpc. Sun et al. (2008) considered an additional local spherical component motivated by the local excess of the synchrotron emission, which might be related to the "local bubble" (see, for example, Fuchs et al. 2009).

Detailed modeling of the Galactic synchrotron emission is beyond the scope of our paper. We will simply assume that there is a thick component with a scale height of ∼1 kpc. We will also assume that there could be an additional local spherical component. Then, in the following subsections, we will consider the relation between the spectrum of three-dimensional (3D) turbulence and the observed angular spectrum of the synchrotron emission.

In Appendix A, we discussed the relation between the 3D spatial MHD turbulence spectrum and the observed two-dimensional (2D) angular spectrum of the synchrotron emission (see Equation (A11) for a quick summary). But, in the appendix we assumed that the emission arises from a spherical region filled with homogeneous turbulence. In this section, we will show that the modulation of the synchrotron intensity of the emitting volume can also affect the observed angular spectrum of the synchrotron emission.

Earlier studies showed that the angular spectrum of the 408 MHz Haslam map has a slope close to −3: C_l ∝ l⁻³ (Tegmark & Efstathiou 1996; Bouchet et al. 1996). Recently, La Porta et al. (2008) performed a comprehensive angular power spectrum analysis of all-sky total intensity maps at 408 MHz and 1420 MHz. They found that the slope is close to −3 for high Galactic latitude regions. Other results also show slopes close to −3. For example, using Rhodes/HartRAO data at 2326 MHz (Jonas et al. 1998), Giardino et al. (2001b) obtained a slope of ∼2.92 for high Galactic latitude regions with |b|>20°. Giardino et al. (2001a) obtained a slope of ∼3.15 for high Galactic latitude regions with |b|>20° from the Reich & Reich (1986) survey at 1420 MHz. Bouchet & Gispert (1999) also obtained a slope of ∼l⁻³ spectrum from the 1420 MHz map.

In general, synchrotron emission from the Galactic disk makes it difficult to measure the angular spectrum of the synchrotron emission from the Galactic halo. In order to avoid the contamination by the Galactic disk, one may mask out the low Galactic latitude regions. This can be done, for example, by setting all synchrotron intensity to zero for regions with |b| < b_cut. However, the angular spectrum obtained with the Galactic mask exhibits spurious oscillations. Moreover, the spectrum obtained with a mask may not be the true one because it is contaminated by the mask. In principle, one may correct such oscillations and estimate the true spectrum using the convolution theorem: the Fourier coefficients (or, in this case, the spherical harmonic coefficients) of the masked data are a convolution of those of the true data and those of the mask. However, the practical implementation of the method is not simple.

Another approach is to first estimate a two-point correlation function and to then extract the angular power spectrum from it (Szapudi et al. 2001; see also Equation (10)). This method is free of artifacts caused by the mask. But the angular spectrum C_l obtained in this way is, in general, noisy and requires a large number of calculations to accurately measure the slope of the spectrum.

We are interested in the slope of the angular spectrum on small angular scales, which is the same as that of the underlying 3D spatial turbulence spectrum (see Appendix A). Therefore, in this paper, we use yet another approach. We first calculate the second-order angular structure function:

$$\begin{equation} D_2 (\theta) = \langle | I({\bf e}_1) -I({\bf e}_2) |^2 \rangle, \end{equation} \tag{ 4 }$$

where I(e) is the intensity of the synchrotron emission, e₁ and e₂ are unit vectors along the lines of sight, θ is the angle between e₁ and e₂, and the angle brackets denote an average taken over the observed region. Then, we extract the slope of the angular spectrum using the relation

$$\begin{equation} D_2(\theta) \propto \theta ^{m-2} \Leftrightarrow C_l \propto l^{-m} \Leftrightarrow E_{\rm 3D}\propto k^{-m} \end{equation} \tag{ 5 }$$

for small angular scales (see Appendix A.2).

We note that excessive care is required in the presence of white noise. In this case, the second-order structure function will be D₂(θ) = 〈|I(e₁) + δ₁ − I(e₂) − δ₂|²〉 = 〈|I(e₁) − I(e₂)|²〉 + 〈|δ₁ − δ₂|²〉, where δ₁ and δ₂ represent the noise (A. Chepurnov 2008, private communication). If the second term on the right (〈|δ₁ − δ₂|²〉) is sufficiently smaller than the first term on the right (〈|I(e₁) − I(e₂)|²〉), we can ignore the noise. Otherwise, the noise can interfere in the accurate measurement of the slope. In the Haslam map, the level of the white noise seems to be negligible. The reason is as follows. Our measurements show that D₂(0015) ∼ 0.05 in the Haslam map. This means that the 〈|δ₁ − δ₂|²〉 term is no larger than 0.05, which is sufficiently smaller than the values of D₂(θ) shown in Figure 1.

Zoom In Zoom Out Reset image size

In the left panel of Figure 1, we show the second-order structure function for the Galactic halo (i.e., |b|>30°). The slope of the second-order structure function lies between those of two straight lines. The steeper line has a slope of 4/3 and the other one has a slope of 1. The actual measured slope is ∼1.2. This result implies that the 3D turbulence spectrum is E_3D(k) ∝ k^−3.2, which is shallower than the Kolmogorov spectrum E_3D(k) ∝ k^−11/3.

Now the question arises: why is the slope shallower than that of the Kolmogorov spectrum? Chepurnov (1999) provided a discussion on the effects of density stratification on the slope. He used a Gaussian disk model and semi-analytically showed that the slope of the angular spectrum can be shallower than that of the Kolmogorov spectrum. In the following subsection, we present our understanding of the effect.

Below we take into account that the emission does depend on the distance from the Galactic plane. Synchrotron emission models (see the previous subsection) assume either an exponential ($e^{-z/z_0}$) or a square of hyperbolic secant (sech^b[z/z₀]) law for synchrotron emissivity, where z is the distance from the Galactic plane and b is a constant. For simplicity, we assume the observer is at the center of a spherical halo. That is, the geometry is not plane-parallel, but spherical. In what follows, we use r, instead of z, to denote the distance to a point.

To illustrate the effects of this inhomogeneity, we test the following three models.

• 1.  
  Homogeneous halo. Turbulence in the halo, and thus the emissivity, is homogeneous. Turbulence has a sharp boundary at d_max = 1 kpc. The outer scale of turbulence is 100 pc. Basically, this model is the same as the one we consider in Appendix A.
• 2.  
  Exponentially stratified halo. Emissivity shows an exponential decrease, $\epsilon (r) \propto e^{-r/r_0}$. We assume r₀ = 1 kpc and the halo truncates at r = 8 kpc. The outer scale of turbulence is 100 pc (see Figure 2).
• 3.  
  Two-component halo. Emissivity decreases as $\epsilon (r) \propto \epsilon _1 e^{-r/r_1}+\epsilon _2 e^{-r/r_2}$, where ₂ = 10₁, r₁ = 1 kpc, and r₂ = 100 pc. The halo truncates at r = 8 kpc. The outer scale of turbulence is 100 pc. The second component mimics the local enhancement of the synchrotron emissivity (see Figure 2).

Zoom In Zoom Out Reset image size

We numerically calculate the angular correlation function w(θ) and the second-order structure D₂(θ) from

$$\begin{eqnarray} &&w(\theta)=\int dr_1 \int dr_2 {\cal K}(|{\bf r}_1-{\bf r}_2|) \epsilon (r_1) \epsilon (r_2), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&D_2(\theta) \propto \overline{T} - w(\theta), \end{eqnarray} \tag{ 7 }$$

where |r₁ − r₂| = r²₁ + r²₂ − 2r₁r₂cos θ, (r) is the synchrotron emissivity, $\overline{T}=\lim _{\theta \rightarrow 0} w(\theta)$, and we use the spatial correlation function ${\cal K}(r)$ obtained from the relation

$$\begin{equation} {\cal K}(r) \propto \int _{0}^{\infty } 4\pi k^2 E_{\rm 3D}(k) \frac{ \sin kr }{ kr } dk, \end{equation} \tag{ 8 }$$

where the spatial spectrum of emissivity E_3D has the form

$$\begin{equation} E_{\rm 3D}(k) \propto \left\lbrace \begin{array}{@{}ll} \mbox{constant} & \mbox{if $k\le k_0$} \\ (k/k_0)^{-11/3} & \mbox{if $k\ge k_0$,} \end{array} \right. \end{equation} \tag{ 9 }$$

which is the same as the Kolmogorov spectrum for k ⩾ k₀ (∼1/L). The reason we use a constant spectrum for k ⩽ k₀ is explained in Appendix B (see also Chepurnov 1999). We obtain the angular spectrum from the relation

$$\begin{equation} C_l \propto \int P_l(\cos {\theta }) K(\cos \theta)\ d(\cos \theta), \end{equation} \tag{ 10 }$$

where P_l is the Legendre polynomial.

In Figure 3, we plot the calculation results. The angular correlation function w(θ) does not change much when θ is small, and follows ∼(π − θ)/sin θ ∼ 1/θ when θ is large. The critical angle is a few degrees for the homogeneous model (thick solid curve) and the single-component exponential model (dotted curve). As we discussed earlier, the critical angle for homogeneous turbulence is ∼(L/d_max)^rad ∼ 6°, where d_max (=1 kpc in our model) is the distance to the farthest eddy. In Figure 3 (left panel), we clearly see that the slope of w(θ) changes near θ ∼ 6°. The second-order structure function D₂(θ) also shows a change in the slope near the same critical angle (θ ∼ 6°). In the single-component exponential model (dotted curve), the value of d_max is not important. Instead, the scale height z₀ is a more important quantity, which is 1 kpc in our model. In the left and middle panels of Figure 3, we observe that the single-component exponential model also shows a change of slope near θ∼ a few degrees. Therefore, we can interpret that the critical angle for stratified turbulence is ∼L/z₀, instead of ∼L/d_max.

Zoom In Zoom Out Reset image size

Now, it is time to answer our earlier question of why the observed slope is shallower than that of Kolmogorov. Let us take a look at the right panel of Figure 3. All three models show that the slope of C_l is almost the Kolmogorov one for l ≳ 200>l_cr ∼ πd_max/L ∼ 30. However, if we measure the average slope of C_l between l = 10 and 200, we obtain slopes shallower than the Kolmogorov one. The single-component model and the homogeneous model give similar average slopes around −3. However, the homogeneous model gives a more abrupt change of slope near l ∼ 30. In fact, the right panel of Figure 3 shows a noticeable break near l ∼ 50. The average slope of the two-component model gives a more or less gradual change of the slope and the observed slope is very close to −3 for a broad range of multipoles l. It is difficult to tell which model is better because the models are highly simplified. Nevertheless, the two-component model looks the most promising, which is not so surprising because the two-component model has more degrees of freedom.

Note that, compared with the two-component model, the single-component model shows a more or less sudden change of slope near l ∼ l_cr ∼ πz₀/L ∼ 30. Therefore, if the single-component model is correct, the scale height z₀ cannot be much larger than ∼10 times the outer scale of turbulence L. If z₀ is much larger than ∼10L, l_cr becomes smaller and we will have an almost Kolmogorov slope for l ≳ 10. We also note that it is possible that the 3D spatial turbulence spectrum itself can be shallower than the Kolmogorov one. That is, it is possible that the spectrum of B(r), hence that of B²(r), can be shallower than the Kolmogorov one. For example, some recent studies show that strong MHD turbulence can have a ∼k^−3.5 spectrum, rather than k^−11/3 (Maron & Goldreich 2001; Boldyrev 2006; Beresnyak & Lazarian 2006).⁹ If this is the case, the observed angular spectrum can be slightly shallower than the Kolmogorov for l>l_cr.

To summarize this subsection, our simple model calculations imply that C_l will be very close to the underlying spectrum of magnetic turbulence for large values of l (l≳ a few times 100). The corresponding spectral slope is expected to be close to the Kolmogorov one (see Section 3). For intermediate values of l (e.g., 10 < l < 200), the average slope is shallower than the Kolmogorov one. Thus, our modeling shows the consistency of the observational spectra with the expectations. Studies of the observed spectra at higher l may be useful for better testing of our predictions.

In the right panel of Figure 1, we show how the second-order structure function changes with the Galactic latitude. The lower curve is the second-order angular structure function obtained from pixels in the range of 28° ⩽ b ⩽ 32°. The middle and upper curves are the second-order angular structure functions obtained from pixels in the range of 8° ⩽ b ⩽ 12° and −2° ⩽ b ⩽ 2°, respectively. The middle and upper curves clearly show the break in the slopes near θ ∼ 3° and ∼15, respectively. When the angular separation is larger than the angle of the break, the structure function becomes almost flat. As we move toward the Galactic plane, the sudden changes in the slopes happen at smaller angles.

There are at least two possible causes for the break in the slope. First, a geometric effect can cause it. As we discuss in Appendix A, the change of slope occurs near θ_c ∼ L/d_max. As we move toward the Galactic plane, the distance to the farthest eddy, d_max, will increase. As a result, the critical angle θ_c ∼ L/d_max will decrease. Therefore, we will have smaller θ_c toward the Galactic plane. This may be what we are observing in the right panel of Figure 1. Second, discrete synchrotron sources can cause flattening of the structure function on angular scales larger than their sizes. Although the map we use was reprocessed to remove strong point sources, there might be unremoved discrete sources. When filamentary discrete sources dominate synchrotron emission, the second-order structure function will be flat on scales larger than the typical width of the sources. In reality, both effects may work together.

In view of the variations of the spatial spectral slope of the synchrotron emission at low Galactic latitudes, the use of the foreground removal procedure discussed in Section 2 is more challenging at those latitudes. At the same time, it should be possible to reliably remove the high-l fluctuations corresponding to high latitudes with the procedure in Section 2 due to the observed regular power-law behavior.

Roughly speaking, the shape of the angular spectrum of polarized synchrotron emission will be similar to that of the total intensity at millimeter wavelengths. However, it is expected that at longer wavelengths, Faraday rotation and depolarization effects should cause flattening of the angular spectrum, which has been actually reported (see de Oliveira-Costa et al. 2003 and references therein). On the other hand, La Porta et al. (2006) analyzed the new Dominion Radio Astrophysical Observatory (DRAO) 1.4 GHz polarization survey and obtained angular power spectra with power-law slopes in the range [ − 3.0, − 2.5]. More observations on the polarized synchrotron foreground emission can be found in Ponthieu et al. (2005), Giardino et al. (2002), Tucci et al. (2002), and Baccigalupi et al. (2001). In this paper, we do not discuss the properties of polarized synchrotron emission. Readers may refer to recent models of polarized synchrotron emission (Page et al. 2007; Sun et al. 2008; Miville-Deschênes et al. 2008; Waelkens et al. 2009).

Let us begin with a model dust emission map created by Finkbeiner et al. (1999), which is available at the NASA LAMBDA Web site⁸. As we will explain later in this section, we can derive a polarized intensity map from this kind of total intensity map.

We note that the difference between the dust emission and synchrotron emission (discussed in the previous section) is expected. The origin of the synchrotron emission is related to cosmic-ray electrons, which are distributed within an extended magnetic halo (Ginzburg & Ptuskin 1976). At the same time, dust is expected to be localized mostly within the Galactic plane. In Figure 4, we present statistical properties of the map. The map shows a rough constancy of the emission for high Galactic latitude region when multiplied by sin b (left panel of Figure 4). The sin b factor also appears in the probability density function (PDF): it makes the PDF more symmetric (right panel of Figure 4). Therefore, it is natural to conclude that a disk component dominates the dust map.

Zoom In Zoom Out Reset image size

As in the Haslam map, we do not try to obtain the angular spectrum of the dust emission map directly. Instead, we use the second-order structure function to reveal the angular spectrum on small angular scales. Indeed, the second-order structure function of the dust map (see Appendix D) shows a slope of ∼0.6, which corresponds to an angular spectrum of ∼l^−2.6.

The slope of the angular power spectrum of the model dust emission map is very similar to that of the original FIR data. Schlegel et al. (1998) found a slope of −2.5 for the original FIR data. On the other hand, other researchers found slopes close to −3 from other observations (see Tegmark et al. 2000 and references therein; see also Masi et al. 2001).

Dust density fluctuations mostly arise from cold dense phases of the ISM (see Draine & Lazarian 1998a for a list of idealized ISM phases). There the turbulence is known to be supersonic, which is vividly revealed, for instance, by Doppler broadening of observed molecular lines from giant molecular clouds (GMCs; see McKee & Ostriker 2007). The shallow spectrum of density fluctuations is consistent with the numerical simulations of supersonic MHD turbulence in Beresnyak et al. (2005). The shallow spectrum was later also reported in supersonic hydro turbulence in Kim & Ryu (2005), which indicates that the effect is not radically changed by the magnetic field. The latter makes the conclusion about the shallow spectrum independent of the degree of ISM magnetization and sub-Alfvénic versus super-Alfvénic character of the turbulence there. Thus, we claim that the observed spectra should be associated with the shallow spectra of underlying density fluctuations in the denser part of the ISM. We predict that the shallow power law arising from density fluctuations extends from the scales of the turbulence energy injection to the dissipation scales. As a result, power-law extension of the observed data and the corresponding filtering using Equation (2) is possible.

In general, it is advantageous to use all possible sources of information about foregrounds in order to improve their removal. In this section, we discuss how the starlight polarization maps can be used to construct the maps of dust polarized emission. In principle, we can construct a polarized dust emission map at millimeter wavelengths ($I_{{\rm {pol}},{\rm mm}}(l,b)$) from a total dust emission map (I_mm(l, b)) and a degree-of-polarization map (P_em,mm(l, b)) at millimeter wavelengths:

$$\begin{equation} I_{{\rm {pol}},{\rm mm}}(l,b)= P_{{\rm em},{\rm mm}}(l,b) I_{{\rm mm}}(l,b), \end{equation} \tag{ 11 }$$

where (l, b) denotes the Galactic coordinate. However, neither $I_{{\rm {pol}},{\rm mm}}(l,b)$ nor I_mm(l, b) is directly available. Therefore, we need indirect methods to get $I_{{\rm {pol}},{\rm mm}}(l,b)$ and I_mm(l, b).

Obtaining a total dust emission map (I_mm(l, b)) is relatively easy because total dust emission maps at FIR wavelengths are already available from the IRAS and COBE/DIRBE observations. Using the relation

$$\begin{equation} I_{{\rm mm}}(l,b) = I_{100\,\mathrm{\mu m}}(l,b)(1{\rm \,mm}/100\,\mathrm{\mu m})^{-\beta }, \end{equation} \tag{ 12 }$$

where 1 ≲ β ≲ 2, one can easily obtain an emission map at millimeter wavelengths (I_mm) from the maps at 100 μm or 240 μm. However, more sophisticated model dust emission maps at millimeter wavelengths already exist. For example, Finkbeiner et al. (1999) presented predicted full-sky maps of microwave emission from the diffuse interstellar dust using the FIR emission maps generated by Schlegel et al. (1998). In fact, the model dust emission map we analyzed in the previous section (Section 4.1) is one of the maps presented in Finkbeiner et al. (1999). Therefore, we can assume that the thermal dust emission map (I_mm(l, b)) is already available.

Obtaining a degree-of-polarization map at millimeter wavelengths (P_em,mm(l, b)) is relatively more complicated. We can use the measurements of starlight polarization at optical wavelengths to get P_em,mm(l, b). The basic idea is that the degree of polarization by emission at the millimeter (P_em,mm) is related to that at optical wavelengths ($P_{{\rm em},\rm optical}$), which in turn is related to the degree of polarization by absorption at optical wavelengths ($P_{\rm abs,\rm optical}$):

$$\begin{equation} P_{\rm abs,\rm optical} \rightarrow P_{{\rm em},{\rm optical}} \rightarrow P_{{\rm em},{\rm mm}}. \end{equation} \tag{ 13 }$$

We describe the relations in detail below.

When the optical depth is small, we have the following relation (see, for example, Hildebrand et al. 2000):

$$\begin{equation} P_{{\rm em},\rm opt} \approx -P_{\rm abs,\rm opt}/\tau, \end{equation} \tag{ 14 }$$

where $P_{{\rm em},\rm opt}$ is the degree of polarization by emission and τ is the optical depth (at optical wavelengths). We obtain polarization by emission at millimeter wavelengths (P_em,mm) using the relation

$$\begin{equation} P_{{\rm em},{\rm mm}}= P_{{\rm em},\rm opt} \left[ \frac{ C_{\rm max}-C_{\rm min} }{ C_{\rm max}+C_{\rm min} } \right]_{{\rm mm}}\Big / \left[ \frac{ C_{\rm max}-C_{\rm min} }{ C_{\rm max}+C_{\rm min} } \right]_{\rm opt}, \end{equation} \tag{ 15 }$$

where Cs are cross sections (of grains as projected on the sky) that depend on the geometrical shape (see, for example, the discussion in Hildebrand et al. 1999; see also Draine & Lee 1984) and dielectric function = ₁ + i₂ (see Draine 1985) of the grains.

For millimeter wavelengths, it is easy to calculate the ratio in Equation (15) because the wavelength λ is much greater than the grain size a (i.e., λ ≫ 2πa). In this case, if grains are oblate spheroids with a₁ < a₂ = a₃ and the short axes (a₁) of grains are perfectly aligned in the plane of the sky, we have

$$\begin{equation} C_j=\frac{ 2\pi V }{ \lambda } \frac{ \epsilon _2(\lambda) }{ (L_j [ \epsilon _1(\lambda)-1 ] +1)^2 + [ L_j \epsilon _2(\lambda) ]^2 }, \end{equation} \tag{ 16 }$$

where L values are defined by

$$\begin{eqnarray} L_1 &=& [(1+f^2)/f^2][1-(1/f)\arctan f], \nonumber \\ L_2 &=& L_3=(1-L_1)/2, \nonumber \\ f^2 &=& (a_2/a_1)^2-1 \end{eqnarray} \tag{ 17 }$$

(see, for example, Hildebrand et al. 1999).

However, for optical wavelengths, the condition λ ≫ 2πa is not always valid and, therefore, the expression in Equation (16) returns only approximate values. For accurate evaluation of the cross sections, one should use numerical methods. Fortunately, several numerical codes are publicly available for such calculations (for example, DDSCAT package by Draine & Flatau 1994, Draine & Flatau 2008; ampld.lp.f by Mishchenko 2000). We use ampld.lp.f to calculate the ratio in Equation (15). We assume that the grains are oblate spheroids (grain size is 0.1 μm, λ_optical = 0.5 μm, and λ_mm = 1000 μm). The left panel of Figure 5 shows that the ratio is around 1.5 when the magnetic field is perpendicular to the line of sight. It also shows that the ratio of $P_{{\rm em},{\rm mm}}/P_{{\rm em},\rm optical}$ is almost independent of the grain axis ratio.

Zoom In Zoom Out Reset image size

In this subsection, we describe a simple way to obtain a polarized map at millimeter wavelengths. However, actual implementation of the method can be more complicated due to the following reasons. First, we use an assumption that the grains that produce optical absorption also produce microwave emission. But, this is not true in general (see Whittet et al. 2008). Second, the expressions in Equations (15) and (16) are valid when the magnetic field direction is fixed and perpendicular to the line of sight and all the grains are perfectly aligned with the magnetic field. If this is not the case, Equation (15) will become $P_{{\rm em},{\rm mm}} \propto P_{{\rm em},\rm optical}$ with a constant of proportionality that depends on the magnetic field structure and the degree of grain alignment. The effect of partial alignment is expected to be less important.¹⁰

The effect of a non-perpendicular magnetic field can be potentially important. We perform a numerical calculation using ampld.lp.f to evaluate this effect. We assume that the grains are oblate spheroids (grain size is 0.1 μm, λ_optical = 0.5 μm, and λ_mm = 1000 μm). The right panel of Figure 5 shows that the polarization ratio drops from ∼1.5 to ∼1.1 when the angle (between the magnetic field and the plane of the sky) changes from 90° to ∼5°. Therefore, the effect is not very strong and can be potentially corrected for.¹¹

After we have constructed a map of the polarized emission from thermal dust, we can obtain the angular spectrum. However, if we are interested in only the shape of the angular spectrum, we do not need to construct the polarized thermal dust emission map. We can get the shape of the angular spectrum directly from the starlight polarization map $P_{\rm abs,\rm optical}(l,b)$.

Equation (11) tells us that $I_{{\rm {pol}},{\rm mm}}$ is given by P_em,mm times I_mm. From Equations (14) and (15), we have

$$\begin{eqnarray} I_{{\rm {pol}},{\rm mm}}& = & P_{{\rm em},{\rm mm}} I_{{\rm mm}} \propto P_{{\rm em},\rm opt} I_{{\rm mm}} \nonumber \\ & \approx & (P_{\rm abs,\rm opt}/\tau) I_{{\rm mm}} \nonumber \\ & \propto & P_{\rm abs,\rm opt}. \end{eqnarray} \tag{ 18 }$$

Here we use the fact that τ ∝ I_mm. Note that the constant of proportionality does not affect the shape of the angular spectrum if grain properties do not vary much in the halo. Therefore, as to the power spectrum C_l of $I_{{\rm {pol}},{\rm mm}}$, we can use that of $P_{\rm abs,\rm optical}$:

$$\begin{equation} C_l \mbox{ of $I_{{\rm {pol}},{\rm mm}}$} \propto C_l \mbox{ of $P_{\rm abs,\rm opt}$}. \end{equation} \tag{ 19 }$$

Once we know the angular spectrum of $P_{\rm abs,\rm optical}$, we can estimate the angular spectrum of $I_{{\rm {pol}},{\rm mm}}$.

The angular spectrum of starlight polarization, P_abs,opt(l, b), is already available. Fosalba et al. (2002) obtained C_l ∼ l^−1.5 for starlight polarization. The stars used for the calculation are at different distances from the observer and most of the stars are nearby stars. The sampled stars are mostly in the Galactic disk. CL02 reproduced the observed angular spectrum numerically using a mixture of stars with a realistic distance distribution. CL02 also showed that the slope becomes shallower when only stars with a large fixed distance are used for the calculation. Therefore, it is clear that distance, or dust column density, to the stars is an important factor that determines the slope. We expect that, if we consider only the nearby stars with a fixed distance, the slope will be steeper. This means that, if we consider stars in the Galactic halo, the slope will be steeper.

The method described above requires measurements of polarization from many distant stars in the Galactic halo. Unfortunately, the number of stars outside the Galactic disk that can be used for this purpose are no more than a few thousands (Heiles 2000; see also discussions in Page et al. 2007 and Dunkley et al. 2009). When more observations are available, accurate estimation of $I_{{\rm {pol}},{\rm mm}}(l,b)$ (and C_l of $P_{\rm abs,\rm optical}$) will be possible.

We expect the fluctuations of the starlight polarization to arise primarily from the fluctuations of the magnetic fields.¹² The latter are expected to have a spectral index close to the Kolmogorov one (see Section 3).

As we discussed in the previous subsection, the angular spectrum, C_l, of $P_{\rm abs,\rm optical}$ for the Galactic halo will be different from the observed l^−1.5 spectrum for a mixture of stars with different distances in the Galactic disk. However, it is not clear exactly how the former is different from the latter. To deal with this problem we use numerical simulations again. We first generate two sets of magnetic fields on a 2D plane (8192 × 8192 grid points) using Kolmogorov 3D spectra.¹³ Since we need $P_{\rm abs,\rm optical}$ for stars well above the Galactic disk, we assume that the distance to the stars is fixed in each model. We consider the following three models.

• 1.  
  Case 1: Nearby stars in a homogeneous turbulent medium. We generate three (i.e., x, y, and z) components of the magnetic field on a 2D plane (8192 × 8192 grid points representing 400 pc × 400 pc), using the following Kolmogorov 3D spectrum: E_3D(k) ∝ k^−11/3 if k>k₀, where k₀ ∼ 1/100 pc. (The outer scale of turbulence is 100 pc.) We assume that the volume density of the dust is homogeneous. All stars are at a fixed distance of 100 pc from the observer.
• 2.  
  Case 2: Distant stars in a homogeneous turbulent medium. We generate three (i.e., x, y, and z) components of the magnetic field on a 2D plane (8192 ×  8192 grid points representing 4 kpc × 4 kpc), using the following Kolmogorov 3D spectrum: E_3D(k) ∝ k^−11/3 if k>k₀, where k₀ ∼ 1/100 pc. Other setups are the same as those of Case 1, but the distance to the stars is 2 kpc.
• 3.  
  Case 3: Stars in a stratified medium. We use the magnetic field generated in Case 1. The volume density of dust shows a sech²(z) decrease: ρ(r) = 4ρ₀/[exp(r/r₀) + exp(−r/r₀)]². We assume a spherical geometry and r₀ = 100 pc. The stars are at r = 200 pc from the observer. The outer scale of turbulence is 100 pc.

We assume that the dust grains are oblate spheroids. In the presence of a magnetic field, some grains (especially large grains) are aligned with the magnetic field (see Lazarian 2007 for a review). Therefore, cross sections parallel to and perpendicular to the magnetic field are different. We assume that the parallel cross section is ∼30% smaller than the perpendicular one. We use the following equations to follow changes of the Stokes parameters along the path (see Martin 1974 for original equations; see also Dolginov et al. 1996):

$$\begin{eqnarray} I^{-1} dI/ds & = & {-}\delta + \Delta \sigma Q/I, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} d(Q/I)/ds & = & \Delta \sigma - \Delta \sigma (Q/I)^2, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} d(U/I)/ds & = & {-}\Delta \sigma (Q/I)(U/I), \end{eqnarray} \tag{ 22 }$$

where δ = (σ₁ + σ₂), Δσ = (σ₁ − σ₂), and

$$\begin{eqnarray} 2\sigma _{1} &=& \sigma _{\perp }, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} 2\sigma _{2} &=& \sigma _{\perp } -(\sigma _{\perp }-\sigma _{\Vert }) \cos \gamma \end{eqnarray} \tag{ 24 }$$

(Lee & Draine 1985). Here, σ_⊥ and σ_|| are the extinction coefficients and γ is the angle between the magnetic field and the plane of the sky. After we get the final values of the Stokes parameters, we calculate the degree of polarization ($\sqrt{Q^2+U^2}/I$) and, then, the second-order angular structure function of the degree of polarization.

We show the result in Figure 6. When all the stars are at a distance of 100 pc (Case 1), the spectrum is consistent with the Kolmogorov one for small θ. The result for the stratified medium (Case 3) also shows a spectrum compatible with the Kolmogorov one for small θ. When the stars are far away (Case 2), the qualitative behavior is similar. However, if we measure average slope between θ = 02 and 20°, the result is different: the slope for Case 2 is substantially shallower. Note that θ = 02 and 20° correspond to l = 1000 and 10, respectively. This means that, when we only have distant stars, the angular spectrum will be shallower than the Kolmogorov one. When we have a mixture of distant and nearby stars, we will have an angular spectrum that is steeper than the case of the distant stars but shallower than the case of the nearby stars, which implies that the spectrum is shallower than the Kolmogorov one. Therefore, it is not surprising that Fosalba et al. (2002) obtained a shallow spectrum of ∼C_l ∝ l^−1.5 for a mixture of nearby and distant stars, mostly in the Galactic disk. Flattening of the spectrum (i.e., C_l ∝ l^−α with α ≈ 1.3–1.4) for polarized FIR dust thermal emission is also observed in Prunet et al. (1998; see also Prunet & Lazarian 1999).

Zoom In Zoom Out Reset image size

Note that the spectrum of the emission polarization for very large values of multipole l is expected to be steeper than the spectrum of starlight polarization measured at small values of l. For instance, from our model calculations, we predict that we should see the Kolmogorov spectrum of polarized emission, rather than ∼C_l ∝ l^−1.5. This difference must be taken into account if starlight polarization is used to filter the polarized microwave emission arising from dust. In fact, for filtering one should either make a model of the magnetic field distribution based on the extended samples of stars throughout the Galactic volume or use only distant stars to obtain the spectrum of polarization similar to that expected in the microwave range.

There could be systematic errors in transforming from the starlight polarization spectrum to that of emission. One such error may arise from the variation of the mean magnetic field direction and the line of sight. We expect that, while the starlight polarization spectrum is relatively insensitive to the mean magnetic field direction, the emission spectrum shows a stronger dependence on it. Our preliminary calculations show that such a systematic error is small. We will pursue this possibility in the future.

In the right panel of Figure 6, we plot the angular power spectrum of starlight polarization. As we mentioned earlier, the angular spectrum of the degree of starlight polarization should be similar to that of polarized thermal dust emission (see Equations (18) and (19)).

To obtain angular spectra, we use a Gauss–Legendre quadrature integration method as described in Szapudi et al. (2001). To be specific, we first generate magnetic fields from the three models we considered in the previous subsection. Then, we calculate angular correlation functions, K(cos θ). Finally, we obtain the angular spectra using Equation (10). Since C_l obtained in this way is very noisy, we plot C_l averaged over the multipole range (l/1.09, 1.09l). We do not show C_l for l>1000 because it is too noisy even with the averaging process. The second-order structure function in the left panel of Figure 6 implies that l(l + 1)C_l ∝ l^−5/3 for l>1000. The straight dashed line for l>1000 reflects this implication.

We normalize the spectra using the condition ∑¹⁰₂(l + 1)C_l/2π = 3(μK²) for the case of the stars in the stratified medium. This normalization is based on the values given in Page et al. (2007, their Equation (25)).¹⁴ We assume that the observed band is W band (ν = 94GHz).

The plot shows that the slopes for nearby stars (thick solid line) and stars in a stratified medium are shallower than that of the Kolmogorov spectrum for l < 1000. This result is consistent with that obtained with the angular structure function. Note that the case of the distant stars has a much flatter spectrum for l < 1000. We believe that our toy model for the stratified medium (Case 3) better represents the actual situation for polarized emission from thermal dust in the Galactic halo. Therefore, we expect that the polarized thermal emission from thermal dust in high-latitude Galactic halo has a spectrum slightly shallower than the Kolmogorov spectrum for l < 1000. Our calculations do not tell us about the slopes for l>1000. However, judging from the behavior of the angular structure function for θ ≲ 01, we expect that C_l ∝ l^−11/3 for l>1000 (see the straight dashed line for l>1000 in the right panel of Figure 6; see also discussions in Section 6.4).

We also show the polarized CMB "EE" spectrum in the figure (data from CMBFAST online tool at http://lambda.gsfc.nasa.gov/). The figure shows that the EE spectrum dominates polarized thermal dust emission from high-latitude Galactic halo for l ≳ 100. The EE spectrum is expected to be sub-dominant when l>5000.

Removal of Galactic foregrounds has always been a major concern for CMB studies. The challenge is only going to substantially increase now, when CMB polarization studies are attempted.

The knowledge that foregrounds are not an arbitrary noise, but have well-defined statistical properties in terms of their spatial power spectra, is an important additional information that can be utilized to evaluate and eventually eliminate the foreground contribution.

Utilizing the information about underlying turbulence power spectrum is not straightforward, however. Our study shows that the observed power spectrum may depend on the geometry of the emitting volume. Therefore, detailed modeling of the foreground fluctuations should involve accounting for the geometry of the emitting volume.

The latter point stresses the synergy of Galactic foreground and CMB studies. Indeed, our fitting of the power spectra in Figure 3 shows that on the basis of their variations we may distinguish between different models of the emitting turbulent volume. As soon as this is achieved, one can predict, for instance, the level of fluctuations that are expected from the foreground at scales smaller than those studied.

While the previous statements are true in general terms, a number of special cases in which simpler analysis is applicable are available. For instance, a simplification that is expected at higher resolutions that are currently available is that at sufficiently small scales the statistics should become independent of the large-scale distributions of the emitting matter. Moreover, simple power laws are expected and observed (see Figure 1) for the foregrounds at high Galactic latitudes. Therefore, modeling of the distribution of the Galactic emission and the filtering of it using the approach in Section 3 are not necessarily interlocked problems.

We remind the reader that the approach in Section 3 requires first separating a particular component of a foreground over a range of scales,¹⁵ e.g., using the traditional technique of frequency templates, and then extending the spatial scaling of the foreground's C^F_l to higher l. Consider a few examples of utilizing this approach.

For instance, high-resolution measurements of the South Pole Telescope (SPT; see Lueker et al. 2010) and the Atacama Cosmology Telescope (ACT; see Fowler et al. 2010) will provide measurements at high l, but will have limited frequency coverage to remove the foregrounds. Therefore, the spatial extrapolation of C^F_l obtained with other low-resolution experiments, which, however, provide good frequency coverage required for the C^F_l identification, may be advantageous.

Prior to the release of the Planck data, the extension of C^F_l obtained, for instance, with WMAP to higher l may be useful for the foreground filtering for the suborbital missions. After the release of the Planck data, studies of the advocated approach can still be useful (consider, for instance, Figure 6). If Planck measures the spatial spectrum up to l = 2000, then for a higher resolution balloon mission one can evaluate the level of foreground contamination by extrapolating the expected foreground spectrum. Moreover, the study of the B-modes would require new experiments with higher sensitivity and the extrapolation procedure can be useful again.

In this subsection, we demonstrate how the filtering process works. For simplicity we use a Cartesian coordinate system. For demonstration purposes, we consider the CMB EE spectrum (see the thin solid line in Figure 6(b)) and Galactic dust foreground for the case of a stratified medium (see the dashed line in Figure 6(b)).

Suppose that low-resolution all-sky maps are already available for many frequency channels. Let us assume that the resolution of the maps is Δθ ∼ 21', which is similar to the WMAP resolution. Since there are many channels, one may use the usual filtering techniques to separate the CMB and the foreground signals. However, since we already know both the CMB and the foreground spectra in advance in this example, we do not follow the usual filtering techniques. Instead, we just assume that we already know the CMB and the foreground spectra for l ≲ 500.

The "all-sky" map in the upper left panel of Figure 7 represents this low-resolution map. The map is defined on a grid of 1024 × 1024 points. The angular resolution of the map is Δθ ∼ 21'. We generate the map from a much higher resolution all-sky map (see upper middle panel), which has a dimension of 8196 × 8196 in this example. When we obtain the low-resolution all-sky map from the much higher resolution all-sky map, we apply a circular top-hat beam pattern with radius ∼105 to mimic the actual observation. We calculate the spectra of the CMB and the foreground by direct Fourier transform of the low-resolution map data.

Zoom In Zoom Out Reset image size

Then, suppose that we perform a high-resolution balloon experiment that can cover only part of the sky. In such balloon experiments, frequency channels are usually limited and removing foregrounds is a challenging task.

Here we show that, if we know the foreground spectrum, we can easily obtain the CMB spectrum from the balloon data. The upper right panel of Figure 7 represents the high-resolution balloon data. The size of the data is 128 × 128. The angular resolution is Δθ ∼ 52 and the map covers 11° × 11° in the sky. We obtain the map by skipping every other point in each direction in the lower right corner of the original map with 8192 × 8192 points.

Our goal is to obtain the angular spectrum of the CMB signal using the balloon data. We first make the balloon data periodic by proper reflections and translations and obtain the periodic data on a grid of 4096 × 4096 points. Then we multiply the data by a Gaussian profile of width ∼11°. The center of the Gaussian profile should locate near the center of the newly constructed 4096 × 4096 pixel data and coincides with the center of the 128 × 128 pixel original data. Then, we perform Fourier transformation of the resulting data of size 4096 × 4096. In this way, we obtain an angular spectrum of the total (i.e., CMB+foreground) fluctuations.

Now, we derive an angular spectrum of the foreground. We first take the angular spectrum from the low-resolution map, which is already available (see the lower left panel of Figure 7). The dashed curve is the foreground spectrum. The foreground spectrum is not defined for l ≳ 500. In principle, we know the foreground spectrum when we know the geometry and turbulence spectrum. Here, we simply assume that the spectrum for l ≳ 500 is Kolmogorov. In this way, we can construct the foreground spectrum for all values of l.

The remaining task is just subtraction: when we subtract the foreground spectrum from that of the total fluctuations, we get the CMB spectrum. The solid curve in Figure 8 is the CMB spectrum obtained in this way. The solid curve shows very good agreement with the original CMB spectrum used for generating the original all-sky data on a grid of 8192 × 8192 points (see the middle panels). The rms relative percentage error (100ΔC_l/C^CMB_l) for 100 < l < 2000 is 24%, where ΔC_l is the difference between the estimated CMB spectrum (C^estimated_l; solid curve) and the true CMB spectrum (C^CMB_l; dotted curve).

Zoom In Zoom Out Reset image size

Note that the foreground spectrum depends on the geometry of the emitting regions and the underlying turbulence spectrum. In earlier sections, we discussed the geometry of light-emitting regions. The turbulence spectrum may have some uncertainties (see the Introduction). The Kolmogorov spectrum may be good for a number of cases. However, it is possible that the spectrum deviated from the Kolmogorov one in some regions, such as the molecular clouds. However, such deviation will not affect our filtering process much if the resolution of the balloon data is not far beyond that of the low-resolution all-sky map.

In the previous subsection, we demonstrated a reconstruction of the CMB spectrum. The result is very promising because the dust foreground model we used is not far from a realistic one. Note that we normalized the amplitude of the dust foreground spectrum using the value quoted in an earlier work (Page et al. 2007).

However, there are a couple of issues regarding the previous demonstration. First, in the previous demonstration, we assumed that we know the exact power-law index of the foreground spectrum. In general, we may not know the exact power-law index of the foregrounds. When this is the case, we need to extrapolate the foreground spectrum to larger values of multipoles and figure out to what extent we can safely extrapolate. Second, it is necessary to test how well our method works when the CMB and the foreground spectra are comparable.

In this subsection, we develop a method for deriving the most probable extrapolation of the foreground spectrum. We also test the performance of our technique for the case where the CMB and the foreground spectra are comparable. For this purpose, we scale up the amplitude of the foreground spectrum by a factor of 1000.¹⁶

As in the previous subsection, we assume that multi-channel observations are available for l ≲ 500. Therefore, we have foreground observation for l ≲ 500. However, unlike the previous subsection, we assume that observations are available only for selected values of the multipole l. We mark four such values of l in Figure 9. We use the same foreground spectrum as in the previous subsection to generate the four observed points. We also show the 1σ observation error bars for them, which are arbitrarily assigned in this example.

Zoom In Zoom Out Reset image size

Using the four observed points for foreground, we find the most probable extrapolation for l>500. It is tempting to use the linear least-squares fit in a log –log  plane, in which the x-axis is a logarithm of l and the y-axis is that of l(l + 1)C_l/2π. But, it is not easy to estimate uncertainties in this case.

Therefore, we adopt a slightly different method. We work on the log –log  plane. In this subsection, we use the following conventions: x ≡ log l, y ≡ log [l(l + 1)C_l/2π], and (x_i, y_i) (i = 1, ..., N, where N = 4 in this example) denotes an observed point. For a given x, we want to find the most probable value and the 1σ uncertainty of y. In order to find these, we follow the following procedure. First, we select an arbitrary y and consider all possible lines that pass through (x, y). Let the slope of such a line be a and the y-intercept be b. Then, we calculate the following probability for the chosen y:

$$\begin{equation} P(y)\equiv \exp \left[ -\sum _{i=1}^{N} \frac{(y_i - z_i)^2}{2\sigma _i^2} \right], \end{equation} \tag{ 25 }$$

where z_i = ax_i + b and σ_i is the 1σ observation error for i. Second, we repeat similar calculations for different values of y. Third, we find the value of y where P_tot(y) ≡ ∑_all linesP(y) has the maximum. Fourth, we find 1σ of the distribution P_tot(y). Fifth, we repeat similar calculations for different values of x.

We plot the results of this procedure in Figure 9(a). The solid line denotes the most probable values of "y" and the shaded region represents the 1σ uncertainty. The 1σ uncertainty gets larger as l gets larger when l>500. Therefore, it may not be a good idea to extrapolate for l's an order of magnitude larger than 500, which is the maximum l of the low-resolution multi-channel maps in our current example. However, it seems to be accurate to extrapolate for l's a few times larger than 500 in our current example. We can reduce the 1σ uncertainty by reducing the observation error bars, which may enable us to extrapolate further for larger l's.

As in the previous subsection, we subtract the most-likely foreground spectrum from the observed (i.e., CMB+foreground) spectrum. We plot the result in Figure 9(b). Figure 9(b) shows that our technique recovers the CMB spectrum quite well for l ≳ 550. Note that Figure 9(a) tells us that the CMB and the foreground signals are almost same at l ∼ 550. This example implies that our technique works when the CMB spectrum is slightly larger than the foreground spectrum. It is not likely that our technique in the current form works if the foreground spectrum is larger than the CMB one. However, even in this case, it is possible that our technique can help to reduce observational uncertainties when combined with conventional multi-channel methods.

The spatial filtering that we advocate here may be used as part of a more general filtration procedure that uses both spatial and frequency information. Indeed, it is well recognized that the studies of tensor B-modes provide an excessively severe challenge to the precision of the removal of foregrounds. In this situation, it is important to reduce the errors in the determination of the foreground signal. The customarily used frequency templates provide filtering that is limited by both systematic and measurement errors. In this situation, any additional information that can help in decreasing the errors is highly valuable. The spatial power spectrum of the fluctuations that we discussed in this paper does provide such information. This approach can be illustrated in Figure 9, which currently only has points for low l. However, it can be seen that the initial error bars of the points can be decreased if we require the points to correspond to a power-law spectrum. It is evident that if we had more points at higher l, the uncertainties could be further reduced. In other words, the constraint that the foreground fluctuations follow a power law allows us to partially remove foregrounds at scales at which no foreground templates based on multi-frequency measurements are available and, at the same time, increase the precision of the foreground removal if the templates are available.

The major claim in this paper is that the statistical regularity of the foreground fluctuations enables one to extend their spectra from the scales where observations are available to the scales with no observations. In other words, the spatial spectra of foregrounds are predictable. This, in turn, makes spatial filtering of foregrounds possible.

The predictability of Galactic foreground fluctuations stems from the fact that they are due to ubiquitous Galactic turbulence. MHD turbulence is known to have well-defined statistical properties both in compressible and incompressible limits (see Cho & Lazarian 2005 for a review). Thus, one expects to see a power-law behavior, which in the case of a magnetic field is expected to correspond to the 3D power spectrum close to the Kolmogorov one. In the case of density, the spectrum is Kolmogorov for the subsonic turbulence, but gets shallower than the Kolmogorov one for supersonic turbulence (see Section 4.1). Irrespective of the underlying 3D spectrum being Kolmogorov or not we claim that one can predict the entire spectrum of spatial 2D foreground fluctuations when the measurements over a limited range of scales are available. For this purpose, we need to know the geometric properties of the volume. Our expectations of the underlying spectrum to be, in some instances, e.g., for magnetic fields, close to Kolmogorov only help in increasing the accuracy of our prediction of the 2D foreground spectrum.

With the 2D predicted spectrum, one can then use the procedure of filtering the foreground using Equation (2) (see Section 5.3 for demonstration of the procedure). Alternatively, the good correspondence of the observed spatial spectrum of the foregrounds may serve as an additional proof of the accuracy of the foreground removal procedure. Needless to say, the information of Galactic turbulence and the geometry of the emitting region, which can be a by-product of CMB research, is of astrophysical significance.

A note of warning is due, however. The procedure of statistical filtering that we demonstrated in this paper, as any other procedure of foreground removal, is not ideal. The power-law approximation of C^F_l is definitely not exact. In this paper, we have analyzed the causes for such deviations and provided explanations for the most notable features characterizing the change. More detailed modeling of Galactic turbulence is required.

In this paper, we have discussed angular spectra of Galactic foregrounds. We have focused on synchrotron total intensity and polarized thermal dust emission. Our current study, as well as earlier studies (Chepurnov 1999; CL02), predict that C_l will reveal a true 3D turbulence spectrum on small angular scales.

Our model calculations that take into account stratification effects imply that

• 1.  
  θ< a few times 01 (or l> a few times 100) for synchrotron emission (see Figure 3) and
• 2.  
  θ ≲ 01 (or l ≳ 1000) for polarized emission from thermal dust (see Figure 6).

On larger angular scales, spectra are expected to be shallower.

In this paper, we have analyzed the Haslam 408 MHz map and a model dust emission map, and compared the results with model calculations. We found the following.

• 1.  
  The Haslam map for the high Galactic latitude (b>30°) can be explained by MHD turbulence in the Galactic halo. The measured second-order angular structure function is proportional to θ^1.2, which corresponds to an angular spectrum of l^−3.2. The high-order statistics for high Galactic latitude (b>30°) is consistent with that of incompressible MHD turbulence. Our model calculations show that a two-component model (see Section 3.3 and Figure 3) can naturally explain the observed angular spectrum. The one-component model can also explain the observed slope. But, the slope of the spectrum shows a more abrupt change near l ∼ 30.
• 2.  
  The model dust emission map may not have anything to do with turbulence on large angular scales. That is, we do not find signatures of turbulence in the map.
• 3.  
  Both maps show flat high-order structure functions for the Galactic plane. This kind of behavior is expected when discrete structures dominate the map.

We have described how we can obtain the angular spectrum of the polarized emission from thermal dust in high Galactic latitude regions. Our model calculations show that the starlight polarization arising from dust in high Galactic latitude regions will have a Kolmogorov spectrum, C_l ∝ l^−11/3, for l ≳ 1000 and a shallower spectrum for l ≲ 1000 (Figure 6). We expect that the polarized emission from the same dust also has a similar angular spectrum. That is, we expect that the angular spectrum of the polarized emission from thermal dust is close to a Kolmogorov one for l ≳ 1000.

We have described a new technique of filtering CMB foregrounds. When we have (1) low angular resolution full-sky measurements (such as WMAP data) and (2) high angular resolution partial-sky measurements with limited frequency channels, we can use this technique to derive the CMB angular spectrum for the high-resolution data. In Sections 5.3 and 5.4, we have demonstrated the technique.

The existing confusion in the literature includes naive identification of the 2D spectra of foregrounds with the spectra of the underlying fluctuations. Therefore, for instance, from the fact that the spectral slope of C^F_l differs from the Kolmogorov one, the conclusion about the nature of the fluctuations is made.

In this paper, we have shown that the 2D spectra may have a spectral slope different from the underlying spectral slope of the turbulence. We showed that it is essential to take into account the non-trivial geometry of observations with the observer sampling turbulence along the diverging lines of sight and within the volume where the density of emitters change. For the case of the synchrotron fluctuations, we showed that the observed non-Kolmogorov value of the spectral index of C^F_l can be reconciled with the Kolmogorov-type turbulence.

A notable difference between our study and CL02 is that in the latter we tried to necessarily associate the spectra of foregrounds with the underlying Kolmogorov or the Goldreich & Sridhar (1995) turbulence with a spectral slope of −11/3. The limitations of this approach are evident in view of the establishment of the shallow spectral index of density fluctuations in supersonic MHD turbulence (Beresnyak et al. 2005). These fluctuations, according to our present study, can explain the observed spectra of Galactic dust emission.

The procedure of spatial statistical filtering discussed in Section 2 is very simplified. It should provide satisfactory results for a simple power-law behavior of C^F_l (see Figure 2). In more complex cases, detailed modeling of the emitting volume and/or the use of the filtering as part of a more sophisticated foreground removal procedure is required. We partially addressed this issue in Section 5. Especially, in Section 5.4, we described a method to estimate uncertainties stemming from a power-law modeling of the foreground spectrum. Nevertheless, more improvement is needed for the new technique. At the same time, we believe that our approach can be applied to the removal of foregrounds not only from the CMB data, but, for instance, from the high-z hydrogen statistics studies.

For synchrotron radiation, emissivity at a point r is given by (r) ∝ n(e)|B_⊥|^γ, where n(e) is the electron number density, B_⊥ is the component of the magnetic field perpendicular to the line of sight. The index γ is approximately 2 for radio synchrotron frequencies (see Smoot 1999). If electrons are uniformly distributed over the scales of magnetic field inhomogeneities, the spectrum of synchrotron intensity reflects the statistics of the magnetic field. For small amplitude perturbations (δb/B ≪ 1; this is true for scales several times smaller than the outer scale of turbulence if we interpret B as the local mean magnetic field strength and δb as the random fluctuating field in the local region), if δb has a power-law behavior, the synchrotron emissivity will have the same power-law behavior (see Getmantsev 1959; Lazarian & Shutenkov 1990; Chepurnov 1999). Therefore, we expect that the angular spectrum of synchrotron intensity also reflects the spectrum of 3D MHD turbulence.

When an observer is located inside a turbulent medium, the angular correlation function, and hence the power spectrum, shows two asymptotic behaviors. When the angle is larger than a critical angle, we can show that the angular correlation shows a universal θ⁻¹ scaling. In contrast, when the angular separation is smaller than the critical angle, the angular correlation reflects statistics of turbulence. In this small-angle limit, we can show that the angular power spectrum is very similar to that of turbulence. The critical angle is determined by the geometry. Let the outer scale of turbulence be L and the distance to the farthest eddies be d_max. Then the critical angle is

$$\begin{equation} \theta \sim L/d_{\rm max}. \end{equation} \tag{ A1 }$$

When the angle between the lines of sight is small (i.e., θ < L/d_max), the angular spectrum C_l has the same slope as the 3D energy spectrum of turbulence. Lazarian & Shutenkov (1990) showed that if the 3D spatial spectrum of a variable follows a power law, E_3D(k) ∝ k^−m, then the 2D spectrum of the variable projected on the sky also follows the same power law,

$$\begin{equation} C_l \propto l^{-m}, \end{equation} \tag{ A2 }$$

in the small θ limit. For Kolmogorov turbulence (E_3D ∝ k^−11/3), we expect

$$\begin{equation} C_l \propto l^{-11/3}, \mbox{ if $\theta <L/d_{\rm max}$.} \end{equation} \tag{ A3 }$$

Note that l ∼ π/θ.

In some cases, when we have data with incomplete sky coverage, we need to infer C_l from the observation of the angular correlation function

$$\begin{equation} w(\theta)=\langle I({\bf e}_1) I({\bf e}_2) \rangle, \end{equation} \tag{ A4 }$$

where I(e) is the intensity of synchrotron emission, e₁ and e₂ are unit vectors along the lines of sight, θ is the angle between e₁ and e₂, and the angle brackets denote an average taken over the observed region. As we discuss in Appendix B, when the underlying 3D turbulence spectrum is ∝k^−m (e.g., m = 11/3 for Kolmogorov turbulence), the angular correlation function w(θ) is given by

$$\begin{equation} w(\theta) \propto \langle I^2 \rangle -\mbox{const}\; \theta ^{m-2}, \mbox{ if $\theta < L/d_{\rm max}$.} \end{equation} \tag{ A5 }$$

It is sometimes inconvenient to use the angular correlation function in practice to study turbulence statistics because of the constant 〈I²〉.

A better quantity in small-angle limit would be the second-order angular structure function:

$$\begin{equation} D_2 (\theta) = \langle | I({\bf e}_1) -I({\bf e}_2) |^2 \rangle \end{equation} \tag{ A6 }$$

$$\begin{equation} \hspace{1.3pc}= 2\langle I^2 \rangle - 2w(\theta). \end{equation} \tag{ A7 }$$

Thus, in homogeneous turbulence with 3D spatial spectrum of E(k) ∝ k^−m, we have

$$\begin{equation} D_2 (\theta) \propto \theta ^{m-2}. \end{equation} \tag{ A8 }$$

When we measure the slope of the angular structure function, we can infer the slope of the 3D spatial power spectrum of turbulence.

In this limit, the angular correlation function is more useful than the structure function. Following Lazarian & Shutenkov (1990), we can show that the angular correlation function for θ>L/d_max follows

$$\begin{eqnarray} w(\theta) & = & \int \int dr_1 dr_2 {\cal K}(|{\bf r}_1-{\bf r}_2|), \nonumber \\ & = & \frac{1}{\sin {\theta }} \int _0^{\infty }dr \, r {\cal K}(r) \int _{\theta /2}^{\pi -\theta /2}d\psi \propto \frac{\pi -\theta }{\sin \theta } \sim \frac{\rm const}{\theta },\nonumber \\ \end{eqnarray} \tag{ A9 }$$

where ${\cal K}(r)$ is the 3D spatial correlation function and we change variables: (r₁, r₂) → (r, ψ), which is clear from Figure 10. The Jacobian of the transformation is r/sin θ. We can qualitatively understand the 1/θ behavior as follows. When the angle is large, points along the lines of sight near the observer are still correlated. These points extend from the observer over a distance ∝1/sin(θ/2).

Zoom In Zoom Out Reset image size

If we assume L/d_max < θ ≪ 1, we can get the angular power spectrum C_l using Fourier transform:

$$\begin{eqnarray} C_l & \sim & \int \int w(\theta) e^{-i{\bf l}\cdot {\bf \theta }} d\theta _x d\theta _y \nonumber \\ & \sim & \int d\theta \, \theta J_0(l\theta) w(\theta) \propto l^{-1}, \end{eqnarray} \tag{ A10 }$$

where θ = (θ²_x + θ²_y)^1/2, J₀ is the Bessel function, and we use w(θ) ∝ θ⁻¹.

In summary, for homogeneous Kolmogorov turbulence (i.e., E(k) ∝ k^−11/3), we expect from Equations (A3) and (A10) that

$$\begin{equation} C_l \propto \left\lbrace \begin{array}{@{}ll} l^{-11/3} & \mbox{if $l>l_{{\rm cr}}$} \\ l^{-1} & \mbox{if $l<l_{{\rm cr}}$,} \end{array} \right. \end{equation} \tag{ A11 }$$

which means that the power index α of C_l is¹⁷−1 ⩽ α ⩽ −11/3. We expect the following scaling for the second-order angular structure function:

$$\begin{equation} D_2(\theta) \propto \left\lbrace \begin{array}{@{}ll} \theta ^{5/3} & \mbox{if $\theta <L/d_{\rm max}$ } \\ \mbox{constant} & \mbox{if $\theta >L/d_{\rm max}$.} \end{array} \right. \end{equation} \tag{ A12 }$$

The critical angle θ_cr ∼ L/d_max depends on the size of the large turbulent eddies and on the length of the line of sight. If we assume that turbulence is homogeneous along the lines of sight and has L ∼ 100 pc corresponding to the typical size of a supernova remnant, and that d_max ∼ 1 kpc for a synchrotron halo (see Smoot 1999), we get θ_cr ∼ 6°.