SIMULATIONS OF THE MICROWAVE SKY

We follow the method of Sehgal et al. (2007), constructing a primary microwave background map with the observed large-scale structure and small-scale structure consistent with theoretical expectations. For large scales, we use the WMAP 5-year Internal Linear Combination (ILC) map (Bennett et al. 2003; Hinshaw et al. 2009). For l < 20, we take the harmonic coefficients (a_lm) from the ILC map with no modification. At smaller scales, where the ILC map is smoothed significantly, a Gaussian random realization is added, so that the ensemble average power equals the theoretical power spectrum which maximizes the WMAP 5-year likelihood. Harmonic coefficients are computed to l_max = 8192, and we synthesize a map with HEALPix parameter N_side = 4096 (0859 pixels). Then we convert to a map with HEALPix parameter N_side = 8192 (0429 pixels) using the alter_alm subroutine included in the HEALPix distribution (Górski et al. 2005). The lower resolution (N_side = 4096) map is used to generate the lensed microwave background map (see Section 2.3), as that is more efficient than working with the N_side = 8192 map directly and results in no significant loss of accuracy. The lensed map is then also converted to an N_side = 8192 map.

A simulation of the large-scale structure of the universe is performed using the same methods as previously described in Bode et al. (2007) and Sehgal et al. (2007), and we refer the reader to their papers for additional details. Within a periodic box of comoving side length L = 1000 h⁻¹Mpc, 1024³ dark matter particles are evolved using the Tree-Particle-Mesh (TPM) N-body code (Bode & Ostriker 2003). The particle mass is 6.82 × 10¹⁰ h⁻¹M_☉ and the gravitational spline softening length is set to = 16.276 h⁻¹kpc. This simulation has been previously used by Bode et al. (2009) in developing a model for the intracluster medium.

A putative observer at redshift zero is placed at one corner of the periodic box, and the mass distribution along a past light cone, covering one octant (box coordinates x, y, z > 0) of the sky, is saved as the simulation is running. To create full-sky maps this octant is replicated, as described in Section 3.

At each simulation time step, spanning the redshift range z₁ > z > z₂, particles are found in a thin spherical shell at the comoving distance d corresponding to the light travel time to this observer, spanning the range d(z₁)>d > d(z₂). For z < 3 all of these particles are saved to disk and for z < 10 pixelized versions of the shells are saved, as described below. There are 579 such shells in all; the range Δz = z₁ − z₂ of each shell depends on the time step size set by the TPM code. It decreases from Δz ≈ 0.08 at z ≈ 9 to Δz ≈ 0.03 at z ≈ 3, and to Δz < 0.005 for z < 0.8. For comoving distances larger than the box size, d > 1000 h⁻¹Mpc, there can be some duplication in structures as the periodic box is repeatedly tiled. This tiling is done without random rotations to preserve the periodicity and prevent discontinuities in the large-scale structure. However, any repetition is usually at a different redshift and evolutionary state. In cases where a halo or filament appears twice at similar redshifts, it is viewed at two different angles, making each projected image unique.

For all the shells (z < 10) particles are subdivided by angular coordinates using the HEALPix scheme for pixelation of a sphere. The projected mass and momentum in each pixel of the shell are computed and saved to disk (as described previously in Das & Bode 2008). These pixelized shells are used to model the gravitational lensing (Section 2.3) and high-redshift SZ distortion (Section 2.4) of the microwave background.

For lower redshifts z < 3, the positions and velocities of all particles in the light cone are also saved to disk for postprocessing. The approximately 55 billion particles thus saved occupy 1.3 TB of disk space and are used to model the SZ effect (Section 2.4) in detail. A friends-of-friends (FoF) halo finder, with a standard linking length b_FoF equal to 1/5 of the mean interparticle spacing, is used to identify dark matter halos and tag the particles belonging to them. This halo catalog is used to model infrared (Section 2.5) and radio (Section 2.6) point sources.

As shown in Figure 1, the halo mass functions measured from particles in the light cone agree well with the semi-analytic fitting formula from Jenkins et al. (2001) for a linking length b_FoF = 0.2. The number of halos above a minimum mass per unit redshift per square degree is plotted for four minimum masses: M_min = 6.82 × 10¹², 2 × 10¹³, 1 × 10¹⁴, and 3 × 10¹⁴ h⁻¹M_☉. The first limit, the equivalent of 100 particles, is the lowest mass of halos that are populated with point sources. The second value corresponds to massive halos with enough resolution to model the intrahalo gas in greater detail. For any cluster survey, the minimum detectable mass is likely to be a function of redshift, and the last two minimum mass values reflect the range expected for upcoming SZ surveys. The data points in Figure 1 are measured using all halos in the octant to minimize sample variance, but then are normalized per square degree.

Zoom In Zoom Out Reset image size

Gravitational lensing is a secondary anisotropy of the microwave background caused by the deflection of primary microwave background photons by intervening large-scale-structure potentials (see Lewis & Challinor 2006, for a recent review). Lensing produces non-Gaussianities, smooths out features in the power spectrum and moves power from large to small scales. High-resolution microwave background experiments, such as ACT, SPT, and Planck, will probe the scales corresponding to the multipole range 500 ≲ ℓ ≲ 10, 000 with unprecedented precision. It will be necessary to include the effect of lensing for performing precision cosmological tests with such data sets. By studying the lensing distortions in the microwave sky, we should also be able to reconstruct the projected matter distribution or the deflection field in the universe (Hu & Okamoto 2002; Okamoto & Hu 2003; Hirata & Seljak 2003; Yoo & Zaldarriaga 2008). Such reconstruction will help break parameter degeneracies in the microwave background, making it sensitive to dark energy dynamics and neutrino mass (Smith et al. 2008; de Putter et al. 2009). Also, the cross-correlation of microwave background lensing with various tracers of large-scale structure will lead to better understanding of galaxy environments, structure formation, geometry, and gravity on large scales (Acquaviva et al. 2008; Das & Spergel 2009; Vallinotto et al. 2009). A realistic simulation of microwave lensing internally consistent with (i.e., based on the same large-scale-structure simulation as) the other simulated secondary components, such as the thermal and kinetic SZ effects and infrared and radio point sources, is essential for training the power spectrum and parameter estimation pipeline. Such simulations are also necessary ingredients in developing realistic lens reconstruction algorithms (Amblard et al. 2004; Perotto et al. 2009). Here we present the salient aspects of a full-sky microwave background lensing simulation at arcminute resolution.

A simulation of the lensed microwave background essentially consists of remapping the points of the primary microwave background sky according to deflection angles derived from a simulated convergence field. Several approaches have been proposed for implementing the remapping as well as for generating the deflection field (Lewis 2005; Das & Bode 2008; Carbone et al. 2007, 2008; Basak et al. 2008). In this work, we follow Das & Bode (2008) in all respects, except that we create a full-sky convergence map (instead of the polar cap geometry adopted in that paper) and lens a full-sky primary microwave background map. The details of the algorithm can be found in Das & Bode (2008)—here we delineate the main steps:

• 1.  
  As described in Section 2.2, at each time step of the large-scale-structure simulation, the positions of the particles in a thin shell (of width corresponding to the time step) are recorded onto the octant of a HEALPix pixelated sphere. Knowing the particle mass and the number of particles falling into each pixel, we generate a surface mass density plane (expressed in mass per steradian) by dividing the total mass in a pixel by the solid angle subtended by it. We then subtract the cosmological mean of the surface density to generate a surface overdensity on each shell.
• 2.  
  We then generate the effective convergence map (κ) up to the maximum redshift in the large-scale-structure simulation (z = 10) by combining the surface overdensities in the octant shells with proper geometrical weighting (see, Equations (15) and (20) in Das & Bode 2008). To accurately lens the microwave background, we need to include the contribution to the convergence from redshifts beyond the farthest simulated shell (z > 10). We do this by adding in a Gaussian random realization for this extra convergence based on a theoretical convergence power spectrum (calculating the contribution to the convergence power from z > 10 structure). At the end of this process, we have a realization of the effective convergence $\kappa (\hat{n})$ on the octant.
• 3.  
  We replicate the octant to create a full-sky convergence map. The detailed mapping of the octant to the full sky is described in Section 3.
• 4.  
  We invert the relation, $\kappa = \frac{1}{2} \nabla ^2 \phi$, between the convergence and the effective lensing potential ϕ in harmonic space to obtain the harmonic components of ϕ, and generate the deflection field, $\mbox{\boldmath $\alpha $}(\hat{n})$, using the relation α = −∇ϕ.
• 5.  
  We use the deflection field to find the source-plane positions, θ_s, corresponding to the observed positions, θ, of the rays (centers of pixels) according to: θ_s = θ + α.
• 6.  
  Finally, we sample the unlensed microwave background at the source-plane positions using an accurate interpolation technique to produce the lensed microwave background map.

For the results presented here, we have performed the above steps at the HEALPix resolution N_side = 4096. Then we convert to a map with HEALPix parameter N_side = 8192 using the alter_alm subroutine included in the HEALPix distribution (Górski et al. 2005).

Figure 2 compares the power spectrum of the simulated effective κ map with two theoretical models: one with and the other without nonlinear corrections (following Smith et al. 2003) to the matter power spectrum. The theoretical spectra have been computed with the publicly available CAMB¹¹ code. As expected, the simulated power spectrum agrees better with the theory curve with nonlinear corrections. There are two things worth noticing in the plot: the large scatter at the low multipole end of the power spectrum and the excess power in the simulation beyond ℓ ∼ 2000. The scatter is expected because we generated the full-sky convergence map by replicating an octant. The apparent high multipole excess is a consequence partly of the nonlinear corrections to the theory slightly underpredicting the nonlinear power, and partly of the shot noise due to the discrete nature of the particles.

Zoom In Zoom Out Reset image size

We computed the power spectrum of the lensed microwave background map and compared it with theoretical expectations. Figure 3 shows the excellent agreement between the theoretical power spectrum from CAMB and the power spectrum generated from the simulated lensed map. For completeness, we also show the power spectrum of the unlensed input map and the theoretical power spectrum used to generate it. The simulation confirms the theoretical prediction that lensing smoothes out acoustic features in the microwave background and moves power from large to small scales. Parenthetically, we would like to remark that the theoretical power spectra in Figures 2 and 3 were generated with CAMB using the parameter k_eta_max_scalar set to 20 times l_max_scalar. For k_eta_max_scalar = 2× l_max_scalar (which is normally used) the κ power spectrum, as well as the lensed microwave background power spectrum, become inaccurate at high multipoles; the κ power spectrum is suppressed beyond ℓ ∼ 1000, and the lensed microwave background power spectrum beyond ℓ ∼ 4000.

Zoom In Zoom Out Reset image size

A particularly interesting quantity is the difference between the lensed and the unlensed microwave background maps. Figure 4 shows this quantity on a 7° × 9° patch projected onto a cylindrical-equal-area grid from the HEALPix sphere. The figure shows the striking large-scale non-Gaussian features that lensing induces on the microwave sky. Lensing reconstruction techniques use these large-angular-scale correlations of small-scale features in the microwave background to recover the lensing potential.

Zoom In Zoom Out Reset image size

The SZ effect is a secondary distortion of the microwave background that results when cosmic microwave radiation scatters against free electrons. It is commonly considered to have two main components (Sunyaev & Zeldovich 1970, 1972). The thermal SZ (TSZ) term arises from inverse Compton scattering of the cosmic microwave background (CMB) with hot electrons, predominantly associated with shock-heated gas in galaxy clusters and groups. The kinetic SZ (KSZ) term is a Doppler term coming from scattering with electrons having fast peculiar motions. Another common distinction is that the KSZ effect has a nonlinear component associated with the small-scale intracluster medium and a more linear component coming from the large-scale IGM. The signal arising from the latter was first calculated for the linear regime by Ostriker & Vishniac (1986) and Vishniac (1987) and is often referred to as the OV effect. Below, we first write down the formalism for the nonrelativistic limit and then the more general relativistic case.

Modeling the SZ effect requires knowing the number density n_e, temperature T_e, and velocity v_e of the electron distribution. In the nonrelativistic limit, the change in the CMB temperature at frequency ν in the direction $\hat{n}$ on the sky is given by

$$\begin{equation} \frac{\Delta T}{T_{\rm CMB}}(\hat{n}) = \left(\frac{\Delta T}{T_{\rm CMB}}\right)_{\rm tsz} + \left(\frac{\Delta T}{T_{\rm CMB}}\right)_{\rm ksz} = f_\nu y - b, \end{equation} \tag{ 1 }$$

where the dimensionless Compton y and Doppler b parameters,

$$\begin{eqnarray} y &\equiv& \frac{k_B \sigma _T}{m_e c^2}\int n_e T_e {\rm d}l = \int \theta _{e}{\rm d}\tau, \\ b &\equiv& \frac{\sigma _{\rm T}}{c} \int n_e v_{\rm los}{\rm d}l = \int \beta _{\rm los}{\rm d}\tau, \end{eqnarray} \tag{ 2 }$$

are proportional to integrals of the electron pressure and momentum along the los, respectively. The dimensionless temperature, los peculiar velocity, and the optical depth through a path length dl are given by

$$\begin{eqnarray} \theta _e &\equiv& \frac{k_{\rm B}T_e}{m_e c^{2}} = 1.96\times 10^{-3} \left(\frac{T_e}{{\rm keV}}\right), \\ \beta _{\rm los} &\equiv& \frac{v_{\rm los}}{c} = 3.34\times 10^{-4} \left(\frac{v_{\rm los}}{100\ {\rm km\ s}^{-1}}\right), \\ {\rm d}\tau \equiv \sigma _{\rm T}n_e {\rm d}l &=& 2.05\times 10^{-3} \left(\frac{n_e}{10^{-3}\ {\rm cm}^3}\right)\left(\frac{{\rm d}l}{{\rm Mpc}}\right), \end{eqnarray} \tag{ 4 }$$

respectively. For a typical cluster, the Compton y parameter is expected to be approximately an order of magnitude larger than the Doppler b parameter.

The nonrelativistic TSZ component has a frequency dependence as specified by the function,

$$\begin{equation} f_\nu \equiv x\coth (x/2) - 4, \quad x\equiv h\nu /(k_{\rm B}T_{\rm CMB}), \end{equation} \tag{ 7 }$$

which has a null at ν ≈ 218 GHz. The distortion appears as a temperature decrement at lower frequencies and as an increment at higher frequencies relative to the null. The nonrelativistic KSZ component is independent of frequency, but the sign of the distortion depends on the sign of the los velocity. We chose the convention where v_los > 0 if the electrons are moving away from the observer.

In the general, relativistic case, the change in the CMB temperature at frequency ν in the direction $\hat{n}$ is given by

$$\begin{eqnarray} \frac{\Delta T}{T_{\rm CMB}}(\hat{n}) = \int \biggl [ & \theta _{e} (Y_0 + \theta _e Y_1 + \theta _e^2 Y_2 + \theta _e^3 Y_3 + \theta _e^4 Y_4) \biggr. \nonumber \\ \biggl. & + \beta ^2 \left[\frac{1}{3} Y_0 + \theta _e \left(\frac{5}{6} Y_0 + \frac{2}{3} Y_1\right) \right] \biggr. \nonumber \\ \biggl. & - \beta _{\rm los} (1 + \theta _e C_{1} + \theta _e^2 C_2) \biggr ]{\rm d}\tau, \end{eqnarray} \tag{ 8 }$$

where the Ys and Cs are known frequency-dependent coefficients (Nozawa et al. 1998). Note that Equation (1) contains only the first-order terms in Equation (8), and that f_ν = Y₀. We use Equation (8) to calculate the full SZ signal in these simulations.

Sky maps of the SZ effect are made by tracing through the simulated electron distribution and projecting the accumulated temperature fluctuations onto a HEALPix grid with N_side = 8192. Using the saved simulation data described in Section 2.2, SZ maps are constructed with contributions from three components: massive halos at z < 3, low-mass halos and the IGM at z < 3, and the high-redshift universe for 3 ⩽ z < 10. See Figure 5 for the resulting SZ power spectrum and Figure 6 for a section of the SZ map.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.4.1. Gas Model for Massive Halos at z < 3

The hot gas distribution associated with massive dark matter halos is calculated using a hydrostatic equilibrium model (Ostriker et al. 2005; Bode et al. 2007, 2009). Here we summarize the general method and refer the reader to Bode et al. (2009) for additional details, in particular on calibrating star formation and feedback against recent optical and X-ray observations. In the large-scale-structure simulation, massive dark matter halos with M_FoF > 2 × 10¹³ h⁻¹M_☉ and z < 3 are postprocessed to include intrahalo gas and stars. These halos are expected to be the dominant contribution to the TSZ terms in Equations (1) and (8) for the scales of interest (1000 ≲ l ≲ 10, 000). Lower mass halos are accounted for differently, as we explain later.

For each massive halo, the collection of N-body particles tagged by the FoF finder is enclosed in a nonperiodic, Cartesian mesh that is centered on the center of mass, has one axis parallel to the los, and a cell spacing that is twice the N-body spline softening length. The particles are assigned using a cloud-in-cell scheme to define the matter density ρ_m on the mesh. The gravitational potential ϕ is then calculated using a nonperiodic Fast Fourier Transform (FFT), and the location of the potential minimum is chosen as the halo center. This procedure preserves the concentration, substructure, and triaxiality of the dark matter halo.

Assuming hydrostatic equilibrium and a polytropic equation of state, the corresponding gas density ρ_g and pressure P_g on the same mesh are given by

$$\begin{eqnarray} \rho _g &=& \rho _{g,0}\psi ^{1/(\Gamma -1)}, \\ P_g &=& P_{g,0}\psi ^{\Gamma /(\Gamma -1)}, \end{eqnarray} \tag{ 9 }$$

respectively, where the dimensionless temperature variable,

$$\begin{equation} \psi \equiv 1-\left(\frac{\Gamma -1}{\Gamma }\right)\left(\frac{\rho _{g,0}}{P_{g,0}}\right)(\phi -\phi _0), \end{equation} \tag{ 11 }$$

is related to the gravitational potential. Bode et al. (2009) adopted a polytropic index Γ = 1.2 (consistent with the value of Γ adopted by other workers and the value found in hydrodynamic simulations) and determined the other two free parameters for the polytropic model, namely the central gas density ρ_g,0 and pressure P_g,0, by matching the gas mass fraction and temperature to recent X-ray observations (Vikhlinin et al. 2006; Sun et al. 2009). The electrons are assumed to have the same distribution of temperatures and velocities as the fully ionized gas.

For calculating the KSZ terms in Equations (1) and (8), the gas velocity on the mesh can be modeled in two ways. The velocity field can be taken from the dark matter particles; this approach allows for velocity substructure, although the dispersion would be inconsistent with the hydrostatic equilibrium assumption. Alternatively, the entire gas distribution can be assigned a single, mass-weighted average peculiar velocity; this is the approach we adopt for self-consistency. Comparing the two options, we underestimate the temperature variance 〈ΔT²_ksz〉 from massive halos alone by only ∼4% when ignoring velocity substructure. Correspondingly, the KSZ angular power spectrum from massive halos alone is also underestimated, but only by <1% at l = 1000 and ∼6% at l = 10, 000.

When a given halo is added to the SZ maps, mass and momentum are conserved by requiring that the baryonic mass be equal to the cosmic baryon fraction Ω_b/Ω_m times the halo mass M_FoF. This is accomplished by first sorting the mesh cells by their gravitational potential, from most negative to least negative. The mesh is then traversed in order and cells are tagged until the accumulated baryonic mass, in gas and stars, reaches the mass limit. The outer radius of the gas is usually beyond the cluster virial radius. Remaining cells do not contribute any SZ signal. This procedure results in a triaxial halo gas distribution with a physically motivated cutoff.

We find that the gas associated with these massive halos makes up approximately half of the TSZ flux, but it is the dominant contribution to the TSZ angular power spectrum. It accounts for ∼90% of the total TSZ power for the scales of interest (1000 ≲ l ≲ 10, 000). On the contrary, this gas accounts for only ∼15% of the total KSZ power because of the relative importance of the IGM to this signal, as discussed in Sections 2.4.2 and 2.4.3.

2.4.2. Low-mass Halos and Intergalactic Medium at z < 3

Low-mass halos and the IGM at z < 3 are modeled using all other saved N-body particles not associated with the massive dark matter halos discussed in Section 2.4.1. While we are able to identify halos down to M_FoF = 6.82 × 10¹² h⁻¹M_☉ (100 particles), there is insufficient resolution to accurately apply the hydrostatic equilibrium model. Instead, an alternative approach motivated by the virial theorem is adopted, one that can also be extended to particles not associated with identified halos. The IGM is expected to contribute only a small fraction of the TSZ signal, but it is very important to include it for calculating the KSZ signal.

According to the virial theorem, the internal energy of a virialized halo is twice its gravitational potential energy. During the formation of a halo, the dark matter undergoes violent relaxation and the gas gets shock-heated. The gas temperature can be estimated directly from the halo mass and used to model the TSZ signal (e.g., Cohn & White 2009), but we chose to calculate the gas temperature from the velocity dispersion of the dark matter. This approach is then applicable to all particles rather than just ones in identified halos.

For each particle, the local density ρ_m and velocity dispersion σ_v are calculated using 32 neighboring particles. For calculating the SZ terms, the effective pressure and momentum of the gas are taken as

$$\begin{eqnarray} P_g &=& \frac{\Omega _{\rm b}}{\Omega _{\rm m}} \left[ \frac{(\gamma -1)\rho _m\sigma _v^2}{2\mu m_{\rm H}k_{\rm B}} \right], \\ (\rho v)_g &=& \frac{\Omega _{\rm b}}{\Omega _{\rm m}}(\rho v)_m, \end{eqnarray} \tag{ 12 }$$

where γ = 5/3 is the ideal gas adiabatic index and μ = 4/(3 + 5X_H) = 0.588 is the mean atomic weight for a fully ionized gas. Helium reionization is assumed to occur instantaneously at z = 3. A temperature floor of 10⁴ K is imposed for the photoionized IGM.

When tracing through a given particle and projecting onto a HEALPix grid, the effective solid angle subtended by the particle is taken into account. A nearby, low-density particle can span a number of map pixels and simply assigning it to a single pixel will result in a large and incorrect temperature fluctuation. Mass and momentum are also conserved. Since the entire gas distribution at z < 3 is accounted for, the maps have a smooth transition from the concentrated halos to the diffuse IGM.

We find that the low-mass halos and the IGM make up approximately half of the TSZ flux, but contribute to the angular power spectrum only by ∼10% at l = 1000 and ∼15% at l = 10, 000. This is in agreement with Hallman et al. (2007) who used a high-resolution hydrodynamic simulation and found that while roughly one-third of the TSZ flux comes from low-mass halos (M < 5 × 10¹³M_☉) and the IGM, these components add ≲10% to the TSZ angular power spectrum (for σ₈ = 0.9). Hallman et al. (2009) also found that the angular power spectrum of the filamentary warm-hot intergalactic medium (WHIM; Cen & Ostriker 1999; Davé et al. 2001) is only ≲1% of the total. Furthermore, Hernández-Monteagudo et al. (2006) found that while lower density gas (δ ≲ 10) makes up approximately half of the baryon budget, it gives rise to only ∼5% of the TSZ flux in the local universe.

For low-mass halos in the mass range 6.82 × 10¹² < M_FoF(h⁻¹M_☉) < 2 × 10¹³, there are adequate particle numbers (100–282) to estimate the halo velocity dispersion. The TSZ signal from this component is sufficiently well approximated with the above method, considering that these halos contribute ≲10% to the TSZ angular power spectrum. However, for filaments with typical overdensities of 10, there are only 10 particles per comoving (h⁻¹Mpc)³. While it is difficult to calculate the velocity dispersion accurately, an approximate solution is acceptable since the filaments contribute only ≲1% to the angular power spectrum.

The KSZ signal from the low-mass halos and the IGM are relatively more important. We find that at l = 1000, they account for ∼15% of the total power, comparable to the contribution from massive halos, and at l = 10, 000, they account for ∼40%. The majority of the KSZ signal is from the high-redshift universe (3 < z < 10), which we discuss in Section 2.4.3. The KSZ signal is relatively less sensitive to resolution since the bulk peculiar momentum is well captured in the large-scale-structure simulation. While velocity substructure is not well resolved, we have already demonstrated in Section 2.4.1 that ignoring it only results in percent level differences. However, given the modest resolution of the large-scale-structure simulation, we caution against using the maps for detailed studies of the SZ signal from the WHIM and IGM.

2.4.3. High-redshift Universe at 3 ⩽ z < 10

The high-redshift universe at 3 ⩽ z < 10 is modeled using the redshift shells containing the surface density fields for mass Σ_m and los momentum Σ_mv. At these higher redshifts, it is very reasonable to assume that the gas traces the matter distribution because the scales of interest are mostly still in the linear regime. Furthermore, since both the nonlinear mass and the collapsed mass fraction above a fixed mass rapidly decrease toward higher redshifts, the typical temperature of the gas is predominantly set by photoionization rather than shock-heating.

For each redshift shell, the electron temperature, los peculiar velocity, and column density are taken to be

$$\begin{eqnarray} T_e &=& 10^4\ {\rm K}, \\ v_{\rm los} &=& \frac{\Sigma _{mv}}{\Sigma _m}, \\ N_e &=& \left(\frac{\Omega _{\rm b}}{\Omega _{\rm m}}\right)\left(\frac{X_{\rm H}}{m_{\rm H}}+\frac{Y_{\rm He}}{m_{\rm He}}\right)\Sigma _m. \end{eqnarray} \tag{ 14 }$$

Hydrogen reionization is assumed to occur instantaneously at z = 10 and helium is only singly ionized. The Thomson optical depth for electron scattering is τ_T ≈ 0.08, consistent with the WMAP 5-year results (Dunkley et al. 2009). Since the electron distribution is warm rather than hot, the SZ temperature fluctuations are nonrelativistic and Equation (8) reduces to Equation (1). There is only a very small contribution to the TSZ signal from the warm electron distribution. The average Compton y parameter is only ∼10⁻⁷ over this redshift range. However, there is a substantial contribution to the KSZ signal because of the large, coherent velocities in the IGM.

Since the light cone is constructed by replicating and stacking from the periodic simulation box without random rotations, this procedure introduces excess power in the KSZ signal at l ≲ 500. The excess power shows up in the IGM component, but not in the halo component previously discussed in Section 2.4.1. Therefore, the IGM component of the KSZ, from the full redshift range 0 < z < 10, is filtered before adding it to the final SZ maps. A simple filter w(l) = 1 − exp [ − (l/500)²] is chosen to suppress the large-scale excess, but preserve power at l > 1000. The remaining low-l power is now more consistent with the level expected for the large-scale OV signal.

Since the simple filtering modifies the signal at l < 1000, the maps should not be used to predict the KSZ signal at these scales. Furthermore, we have neglected to model the contribution from the inhomogeneous reionization epoch. This component is expected to be larger than the contribution from the reionized universe (z ≲ 6) at l ≲ 1000 (e.g., Santos et al. 2003; McQuinn et al. 2005; Zahn et al. 2005). The KSZ signal on these larger scales is subdominant compared to the lensed microwave background and to the TSZ signal at frequencies away from the null.

Star formation in the universe leads to the production of dust grains which absorb the ultraviolet radiation in the environment and re-radiate it in the infrared range. In the context of microwave background observations, the redshifted infrared emission of these sources constitutes an important contaminant on small angular scales at frequencies higher than ∼150 GHz. In this section, we attempt to model the angular anisotropies that infrared (IR) galaxies introduce in the microwave sky by combining the cosmological simulation presented above with a model for the infrared galaxy population, which is partially based upon Righi et al. (2008).

2.5.1. The Basics of the Infrared Model

We first approximate the SED of the infrared galaxies by a modified black body of the form

$$\begin{equation} L_{\nu }\propto \phi (\nu, T_{\rm dust}) \equiv \nu ^{\beta }\left[B_{\nu }(T_{\mathrm{dust}})-B_{\nu }(T_{\mathrm{CMB}})\right], \end{equation} \tag{ 17 }$$

where β is a spectral index, T_CMB is the primary microwave background monopole temperature, and T_dust is the dust temperature evolving in redshift as in Blain (1999):

$$\begin{equation} T_{\mathrm{dust}}(z)= \big[T_{\mathrm{0}}^{4+\beta }+T_{\mathrm{CMB}}^{4+\beta }\big((1+z)^{4+\beta } -1\big)\big]^{1/(4+\beta)}. \end{equation} \tag{ 18 }$$

Here T₀ is the dust temperature if the galaxy was at z = 0, and the symbol B_ν(T) denotes the black body spectral intensity at temperature T. As we shall see below, the free model parameters β and T₀ will determine the spectral index for the effective power law describing the change of the infrared intensity versus frequency, I_ν ∝ ν^α (i.e., α = α(β, T₀)).

The star formation rate (SFR) in a given galaxy seems to be linearly correlated with the dust luminosity (Kennicutt relation, Kennicutt 1998). At the same time, there is observational evidence that the SFR reached a maximum at epochs corresponding to z ∼ 2 − 4 (Hopkins & Beacom 2006). These two factors combined suggest introducing a redshift window function for the infrared luminosities in galaxies close to the measured SFR behavior. In our model, we adopt the dependence

$$\begin{equation} W_{{\rm SFR}}(z) \propto \exp {(-1/z)}\; \exp {(-\left(z/3.5\right)^2)}. \end{equation} \tag{ 19 }$$

This window however gives more weight to the redsfhift range z ∈ [1, 2] if compared to the data of Hopkins & Beacom (2006), due mainly to the absence of halos in our simulation at z > 3.

In practice, given the simulation halo catalog, our approach consists of populating each halo with a number of galaxies with a given infrared luminosity. We adopt a mass range limited by M₁, M₂ (M₁ < M₂), in such a way that halos below M₁ will host no infrared galaxies. Halos with masses M ∈ [M₁, M₂] will typically host one infrared galaxy of total infrared luminosity L_IR(M, z) = W_SFR(z) · W_cool(M) · L_⋆ · (M/M₂). For more massive halos, the average number of infrared galaxies will be given by M/M₂, each of them with a total infrared luminosity given by L_IR(M, z) = W_SFR(z) · W_cool(M) · L_⋆. The function W_cool(M) = exp(−M/M_cool) accounts for the fact that bigger halos need longer times to cool down before being able to form stars. L_⋆ defines a luminosity parameter. The parameters M₁, M₂, L_⋆, and M_cool, combined with the spectral parameters β and T_dust introduced above, will define our infrared model. In the simulations, the actual number of infrared galaxies present in a given halo is driven by a Poissonian realization of its average number. We also introduced an effective scatter of 15% in the mean adopted value for the dust spectral index and in the infrared luminosity of each galaxy, partially motivated by the scatter in the effective spectral index of the intensity found in the analyses by Knox et al. (2004).

2.5.2. The Angular Power Spectrum

Below, we calculate the angular power spectrum of our infrared model theoretically. This allows us to approximate the power spectrum efficiently, after adjustments to the model parameters, without generating HEALPix maps. For each IR galaxy, the spectral luminosity L_ν can be computed from the bolometric luminosity and the adopted form of the SED:

$$\begin{equation} L_{\nu } = \frac{L_{\rm IR}(M,z)}{\int d\nu \;\phi (\nu, T_{{\rm dust}})} \phi (\nu, T_{{\rm dust}}), \end{equation} \tag{ 20 }$$

where T_dust is a function of redshift (Equation (18)). For a flat universe, the measured spectral flux is therefore given by

$$\begin{equation} F_{\nu } = \frac{L_{\nu ^{\prime }}}{4\pi \;r^2(1+z)}, \end{equation} \tag{ 21 }$$

where ν = ν'/(1 + z) and r is the comoving distance to the IR galaxy placed at redshift z(r). The average spectral intensity generated by the IR galaxy population can then be written as the sum of the contribution of all galaxies per unit solid angle:

$$\begin{eqnarray} \hspace{-10pt}{\bar{I}}_{\nu } & = \int dr\; r^2\; \int dM\; \frac{dn}{dM}(r)\;N_g(M) \frac{a(r)L_{\nu ^{\prime }}(M,r)}{4\pi \;r^2}. \end{eqnarray} \tag{ 22 }$$

The comoving halo mass function is given by dn/dM, and N_g(M) provides the average number of infrared galaxies in a halo of mass M. The cosmological scale factor is expressed by a(r), and the typical spectral luminosity for each of the infrared galaxies in the halo is given by $L_{\nu ^{\prime }}(M,r)$, with ν' being the frequency observed at redshift z(r). $L_{\nu ^{\prime }}(M,r)$ is determined by the model described in Section 2.5.1 and Equation (20).

However, we are interested in the angular fluctuations of the infrared spectral intensity, and this is linked to the angular variation of the number of infrared sources. The number of infrared galaxies is driven a priori by Poisson statistics. However, halos tend to cluster in overdense regions of the large-scale matter density field, and this introduces another modulation in the infrared galaxy number density that is known as the correlation term (Haiman & Knox 2000; Song et al. 2003; Righi et al. 2008). In terms of the spectral intensity I_ν, it is easy to show that the Poisson term for the angular anisotropies generated by the infrared galaxy population can be written as (Righi et al. 2008)

$$\begin{equation} C_l^P = \int dr\; r^2 \int dM\; \frac{dn}{dM}(r)\;N_g(M) \left(\frac{a(r)L_{\nu ^{\prime }}(M,r)}{4\pi r^2}\right)^2, \end{equation} \tag{ 23 }$$

whereas the correlation or clustering term depends upon the halo power spectrum:

$$\begin{eqnarray} C_l^C &=& \int \frac{dr}{r^2} P_m(k,r) \biggl [ \int dM \frac{dn}{dM}(r) b(M,r) N_g(M)\nonumber\\ &&\times \frac{a(r)L_{\nu ^{\prime }}(M,r)}{4\pi } \biggr ]^2, \end{eqnarray} \tag{ 24 }$$

where k = (l + 1/2)/r. The bias factor b(M,r) expresses how halos are distributed compared to dark matter, in such a way that b²(M, r) is defined as the ratio of the halo power spectrum and the matter power spectrum (P_m(k, r)). Here we have adopted the model by Sheth & Tormen (1999), which provides a good fit (few percent accuracy) to the simulated halo power spectrum. The last two equations are obtained after using the Limber approximation, which should be accurate to better than 2% for l > 20 (e.g., Verde et al. 2000).

2.5.3. Observational Constraints

2.5.4. Final Model

We modeled two different source populations: the basic population essentially defines the properties of the whole IR model, whereas the very high flux population matches the high flux IR source counts detected by BLAST (BLAST team 2009, private communication). Each of these two models requires a different choice for the parameter set, as shown in Table 1.

Table 1. Parameters for the Two Populations Considered for the Infrared Point Source Model

Model	β	T0	M1	M2	Mcool	L⋆
		(K)	(M☉)	(M☉)	(M☉)	(L☉)
Basic	1.4	25	2.5 × 1011	3 × 1013	5 × 1014	1.25 × 1012
Very high flux	1.3	40	5 × 1014	1 × 1015	8 × 1014	6 × 1014

Download table as:  ASCIITypeset image

In order to avoid the limitation imposed by the minimum mass of the halo resolved in our simulation, we added a mock catalog of halos of masses between 10⁹ M_☉ and the minimum resolved halo mass, M_min = 6.82 × 10¹² h⁻¹M_☉ = 9.61 × 10¹² M_☉. In order to do this, we assumed that the mass function in this mass range was dictated by the fit of Jenkins et al. (2001). We used the population of low-mass halos resolved in the simulation as tracers of the mock halos appended in this low-mass range. We divided the mass range [10⁹ M_☉, M_min] in logarithmically spaced bins, and, in each of these bins, computed the average number of halos per resolved halo that is found in a thin mass interval centered at M_c = 1.10 × 10¹³ M_☉. That is, for each resolved halo of mass close to M_c, we computed the average number of halos present in each of those low-mass bins, and subsequently made a Poissonian realization of this average. The masses of these new mock halos are randomly distributed uniformly within the corresponding mass bin width. These mock halos were randomly distributed (with a Gaussian profile) in spheres of roughly 100 times the parent (resolved) halo virial radius, which should correspond to ∼10 Mpc. Therefore, for each real halo of mass close to M_c, we populated a cloud of lower mass halos in a sphere of roughly 10 Mpc radius and centered on the parent halo itself. This procedure provides our simulation with a low-mass halo catalog that is not Poisson distributed, but rather is correlated to the population of small but resolved halos close to M = M_c. In practice, only mock halos of mass above M₁ were generated, since they are the ones eventually hosting IR galaxies.

The very high flux population consists of few (∼300 in the whole octant) infrared galaxies above 200 mJy at 353 GHz. This very high flux population models the tail observed by BLAST at 250, 350, and 500 μm (BLAST team 2009, private communication). The high flux tail should be at or above the 10 mJy level in the 145–220 GHz channels, contributing significantly to the power spectrum. Some of this very high flux population includes sources believed to be lensed and magnified by the intervening structure (Lima et al. 2009).

The total cumulative source counts above a given flux threshold for the sum of the basic and the very high flux populations are displayed in the left panel of Figure 7. We show fluxes above 1 mJy (which is about SCUBA's detection threshold), although our IR source population extends to dimmer fluxes. It is actually the numerous but dim IR source population that is responsible for most of the CIB. In our model, our IR source population generates between 50% and 80% of the total CIB detected by FIRAS (middle panel). This is consistent with BLAST observations (Marsden et al. 2009; Devlin et al. 2009), which indicate that most, and perhaps all, of the CIB is generated by infrared sources.

Zoom In Zoom Out Reset image size

In the right panel, we display the expected angular power spectrum of our total IR source population. At very high multipoles the Poisson term must dominate, whereas the clustering term takes over at l ∼ 1000. Again due to the construction of the light cone by replicating and stacking the periodic N-body simulation box without random rotations, excess power is introduced in the infrared point source signal at l ≲ 300, as was the case for the KSZ signal from the IGM (see Section 2.4.3). We do not filter this signal to suppress the low-l excess as we did for the KSZ, as this signal consists of discrete sources instead of a smooth component, and such a filter would disturb the real space point source distribution. For l > 300, the theoretical expectation of the angular power spectrum coincides with the actual power spectrum from the simulated infrared map.

The particular values for the model parameters for each of the two populations are displayed in Table 1. Note that the masses in this table refer to FoF halo masses. For the basic population, the values of the dust temperature and infrared luminosity are not far from those found by Dye et al. (2009) when comparing BLAST infrared sources to radio and mid-infrared counterparts in the Chandra Deep Field South. For the few galaxies in the very high flux population, characteristic masses and infrared luminosities must be increased by roughly an order of magnitude, whereas the typical dust temperature increases from 25 to 40 K.

Our model for radio point sources follows in spirit the treatment of the subject by Wilman et al. (2008, hereafter W08). Similar to W08,we separate the radio-loud AGNs into two populations, a high-power family that evolves strongly with redshift, and a low-power one with mild redshift evolution. These two populations can be roughly identified with the type II and I classification of radio sources by Fanaroff & Riley (1974), respectively (however, see discussions in Willott et al. 2001). We also follow W08 and Jackson & Wall (1999) to employ a relativistic beaming scheme to separate the contributions to the total flux from the compact core and the extended lobes of radio sources.

Unlike W08, we do not consider radio-quiet AGNs and star-forming galaxies in the radio source model. As our primary goal is to simulate the microwave sky at relatively bright flux limits (e.g., ≳0.1 mJy), ignoring these populations has negligible effects. Furthermore, the infrared source modeling presented in Section 2.5 is a self-contained model for the star-forming galaxies, and therefore the combined infrared and radio source models should adequately account for the majority of sources of interest to the current generation of microwave experiments.

Our model is a way to populate dark matter halos with radio sources in a Monte Carlo fashion, according to simple prescriptions on (1) the SED of the sources, (2) the number of radio sources as a function of halo mass and redshift, and (3) the (radio) LD of the sources. Below we outline the most important aspects of the model.

2.6.1. Compact Versus Extended Emission of Radio Sources

The extended lobes of radio sources are characterized by a steep falling spectrum. On the other hand, the emission from the core region of a radio source is usually flatter and can exhibit a more complex spectral shape (see, e.g., Lin et al. 2009b). In addition, when the jet axis is close to the observer's los, the fluxes from the core (and/or the base of the jets) will be Doppler-boosted.

Following W08, we assume the observed core-to-lobe flux ratio, R_obs, is related to the intrinsic value, R_int, via R_obs = R_intB(θ), where θ is the angle between the jet axis and the los, B(θ) = [(1 − βcos θ)⁻² + (1 + βcos θ)⁻²]/2, $\beta =\sqrt{1-\gamma ^{-2}}$, and γ is the Lorentz factor of the jet.

Given the intrinsic luminosity L_int = L_c,int + L_l,int of a source, where the subscripts c and l refer to the core and lobes, the beamed luminosity is then

$$\begin{eqnarray} L_{\rm beam} & =& L_{\rm c,beam}+L_{\rm l,int} \nonumber \\ & =& (1+R_{\rm obs}) L_{\rm l,int} \nonumber \\ & =& L_{\rm int} (1+R_{\rm obs})/(1+R_{\rm int}). \end{eqnarray} \tag{ 25 }$$

We use empirically derived values of R_int and γ from W08, and we assume a uniform distribution of cos θ in order to separate L_c,beam and L_l,int given L_beam.

Since the compact core emission is usually observed to exhibit a complex spectral shape (Lin et al. 2009b), we model it with a third-order polynomial, log f_ν = ∑³_ia_i(log ν)ⁱ, where ν is in GHz. We assume the extended sources obey f_ν ∝ ν^−0.8, which is also consistent with observations (Lin et al. 2009b).

2.6.2. Halo Occupation Number and the 151 MHz RLF

Our first task is to reproduce the local 151 MHz RLF, which is effectively a convolution of the halo mass function, the halo occupation number (HON), and the radio LD. The main reason for targeting the RLF at 151 MHz is because at such low frequencies, one can safely assume that the fluxes are mostly dominated by the steep spectrum, extended sources, and the resulting RLF is not very sensitive to biases due to orientation effects. In addition, this RLF is not greatly affected by synchrotron/inverse Compton losses (Blundell et al. 1999). For both FR I and FR II populations, we assume their HON is given by N(M) = N₀(M/M₀)^α. As for the LD, we assume it is independent of halo mass, and adopt two-segment power laws that roughly mimic the observed RLF:

$$\begin{equation} p(L) = \left\lbrace \begin{array}{@{}ll} (L/L_b)^m & \mbox{if $L>L_b$} \\ (L/L_b)^n & \mbox{otherwise}. \end{array} \right. \end{equation} \tag{ 26 }$$

We explore the parameter space spanned by the HON and LD for both radio source populations, until the combined RLF fits the observed one at z ∼ 0.1, which is based on the 3CRR survey (Laing et al. 1983). A comparison of the model prediction and the observed 151 MHz RLF is shown in Figure 8.

Zoom In Zoom Out Reset image size

2.6.3. Redshift Evolution and Source Number Counts

Once our model successfully reproduces the low-redshift 151 MHz RLF, the radio sources are divided into two populations, roughly corresponding to the type I and II classification of Fanaroff & Riley (1974). These populations are given different number density evolutions with redshift, so that both the high-z RLFs and source counts at 151 MHz can be fitted. For FR I sources, we consider an evolution of the form (1 + z)^δ up to z_p (and constant at z > z_p); for FR IIs, we find that a much stronger redshift evolution is needed and, for convenience, adopt an asymmetric Gaussian centered at z_p, with dispersions σ_l and σ_h at z ⩽ z_p and z > z_p, respectively. (Note that there are many more FR Is and FR IIs at high redshift; for example, the comoving density of model FR II sources is ∼200 times higher at z_p = 1.3 compared to at z ∼ 0.) In our model, z_p is a free parameter, which is adjusted to match observations.

Having matched the high-z 151 MHz RLF, we assume each source at 151 MHz is described by L_beam = L_c,beam + L_l,int, where the lobe component dominates. Given R_int and γ from W08, and cos θ picked from a uniform distribution, we determine L_c,beam. We next adjust the shape of the core SEDs (by varying the coefficients a_i) so that the source counts at ν>151 MHz are reproduced. These source counts at higher frequencies are dominated by the radio cores. In Figure 9, we show a comparison of our model prediction (red solid points) with observed source counts at 1.4, 33, 95, and 145 GHz (see caption for observation references). In general, the model provides a good fit to the data.¹² We note that for the 95 GHz panel, we are comparing our model to the predictions of the de Zotti et al. (2005) model.

Zoom In Zoom Out Reset image size

There are 24 free parameters in the model, which fit 24 observational constraints (12 RLFs at various frequencies and redshifts and 12 source count observations at 12 frequencies). (Note that a₀, a₁, a₂, a₃, and α are set to be the same for both FR Is and FR IIs.) Several of the parameters are degenerate; in Table 2, we present a parameter combination that can reproduce the results presented here.

Table 2. Model Parameters for Radio Sources

Type	Rint	γ	a0	a1a	a2a	a3a	N0	M0	α	Lb	m	n	δ	zp	σl	σh
								(h−171M☉)		(W/Hz)
FR I	10−2.6	6	0	−0.12, 0.07	−0.34, 0.99	−0.75, − 0.23	1	4 × 1013	0.1	1024	−1.55	0	3	0.8	⋅⋅⋅	⋅⋅⋅
FR II	10−2.8	8	0	−0.12, 0.07	−0.34, 0.99	−0.75, − 0.23	0.015	3 × 1015	0.1	1027.5	−1.6	−0.65	⋅⋅⋅	1.3	0.4	0.73

Note. ^aWe draw random uniform variates in this range when assigning the spectral shape for cores.

Download table as:  ASCIITypeset image

For a more complete description of the model and further comparison with available observational constraints see Y.-T. Lin et al. (2010, in preparation).

For the contribution of Galactic thermal dust emission, we use the "model 8" prediction from Finkbeiner et al. (1999), an extrapolation to microwave frequencies of dust maps from Schlegel et al. (1998). We use the software included with these maps to interpolate onto an oversampled HEALPix grid with N_side = 8192. This model is a two-component fit to IRAS, DIRBE, and FIRAS data, and is shown by Bennett et al. (2003) to be a reasonable template for dust emission in the WMAP maps. See Figure 10 for a section of the Galactic dust map.

Zoom In Zoom Out Reset image size

The Galactic synchrotron and free–free emission is comparable to the Galactic dust emission at ∼70 GHz and becomes dominate/subdominant at lower/higher frequencies (see, e.g., Smoot 1998). Thus these signals are not expected to be a significant contaminant for ACT or SPT (which surveys target regions away from the Galaxy), although they will be important to consider for the Planck low-frequency channels. We leave the inclusion of these signals to later work, and note work by, for example, Leach et al. (2008), which discusses the inclusion of these signals in microwave simulations for Planck.

Below we list some general properties of the maps themselves. The full-sky HEALPix maps for this simulation are in the Celestial coordinate system with RING ordering (Górski et al. 2005). The conversion between the HEALPix θ and ϕ coordinates and right ascension and declination is ϕ = R.A., 90 − θ = decl. We store all maps in units of flux per solid angle (as Jy sr⁻¹). To convert to temperature (as μK) we use the conversion:

$$\begin{equation} \frac{dT}{dI}=\frac{c^2}{2 k_{B} \nu ^{2}}\frac{(e^x - 1)^2}{x^2 e^x}, \end{equation} \tag{ 27 }$$

where x = hν/k_BT_CMB. To convert from SI to more convenient units requires μK sr Jy⁻¹ = 10²⁰ K sr m² Hz W⁻¹. When calculating the power spectrum with Polspice¹³ or anafast,¹⁴ a correction for the pixel window function is required.

The properties of the halos that are specified in the halo catalog are redshift, R.A., decl., position of the halo potential minimum, peculiar velocity, FOF mass, virial mass, virial gas mass, virial radius, integrated Compton y parameter, integrated kinetic SZ, integrated full SZ at the six frequencies (the latter three integrated within the virial radius, R₂₀₀, and R₅₀₀), the stellar mass, and the central gas density, temperature, pressure, and potential. In the radio and infrared galaxy catalogs, the properties given for each source are redshift, R.A., decl., and the flux in mJy at the six frequencies.

The large-scale structure octant is mirrored into full-sky maps as follows:

$$\begin{eqnarray} &&{\displaystyle \begin{array}{@{}l@{}} \theta _{N} = \theta _{0} \\ \phi _{N} = \phi _{0} + n\displaystyle\frac{\pi }{2} \\ \theta _{S} = \pi - \theta _{0} \\ \phi _{S} = \displaystyle\frac{\pi }{2} - \phi _{0} + n\displaystyle\frac{\pi }{2}, \end{array}} \end{eqnarray} \tag{ 28 }$$

where n = 0, 1, 2, and 3, and the subscripts N and S refer to the northern and southern hemispheres of the map. The subscript 0 refers to the coordinates in the original octant. For the power spectrum analysis, this method for repetition in the Southern hemisphere is preferable to reflecting the north. Although simple mirroring eliminates real space discontinuities across the equator, it makes the simulated sky an even function with respect to the polar axis. In the spherical harmonic transform, the associated Legendre functions are even or odd depending on l. Integrating mirrored hemispheres over them yields zero for odd-l harmonic coefficients. Because rotating the sky transforms harmonics with common l into each other, odd-l harmonics will be zero regardless of the orientation of the mirroring.

In Figure 11, we show the power spectra of the lensed microwave background, infrared galaxies, radio galaxies, and the full SZ signal from the simulated maps at 148 GHz. The SZ power peaks at around l ∼ 3000 and the clustering of the infrared galaxies is apparent for l < 1000, consistent with measurements by Spitzer (Lagache et al. 2007) and BLAST (Viero et al. 2009). No infrared or radio sources have been removed in creating the power spectra, which were made with Polspice on the full-sky maps. The power from the infrared and radio sources, without masking out bright sources, is comparable to the lensed microwave background between l ∼ 2000 and l ∼ 3000. In Figure 11, the magenta line through the blue radio source power spectrum gives the radio source power suggested by the Toffolatti et al. (1998) model scaled by 0.64, which was found by Wright et al. (2009) to be a good fit to the WMAP 5-year source counts. The two power spectra are in excellent agreement.

Zoom In Zoom Out Reset image size

We investigate the SZ flux versus cluster mass scaling relation for this simulation below. An accurate calibration of this relation off observations is important for assessing the selection function of a potential SZ selected cluster sample via simulations. The gas physics input into the simulations will determine the surface brightness of a given cluster based on its mass and redshift, which in turn determines how likely it is to be detected by a given observation strategy.

To characterize the SZ flux—cluster mass scaling relation in this simulation we compare the integrated Compton y parameter to cluster mass for the halos in the one unique octant. The integrated Compton y parameter is given by Y = (∫ydΩ)d²_A(z), where y is defined in Equation (2), dΩ is the solid angle of the cluster, and d²_A(z) is the cluster's angular diameter distance. For the self-similar model, we expect a virialized halo of mass M_vir to have a virial temperature of

$$\begin{eqnarray} &&T_{\mathrm{vir}}\propto [M_{\mathrm{vir}} E(z)]^{2/3}, \end{eqnarray} \tag{ 29 }$$

where for a flat ΛCDM cosmology

$$\begin{eqnarray} &&E(z)=[\Omega _{m}(1+z)^3 + \Omega _{\Lambda }]^{1/2}. \end{eqnarray} \tag{ 30 }$$

If galaxy clusters were isothermal, then the Compton y parameter would satisfy the relation

$$\begin{eqnarray} &&Y \propto f_{\mathrm{gas}} M_{\mathrm{halo}} T, \end{eqnarray} \tag{ 31 }$$

where f_gas is the cluster gas mass fraction. Combining the above equations, we would expect

$$\begin{eqnarray} &&Y\propto f_{\mathrm{gas}} M_{\mathrm{vir}}^{5/3} E(z)^{2/3}, \end{eqnarray} \tag{ 32 }$$

for the self-similar relation. Deviations from the self-similar relation are not unexpected as clusters are not isothermal, not always in virial equilibrium, and f_gas need not be constant with variations in cluster mass and redshift.

The gas physics in this simulation has been calibrated off of X-ray observations of cluster gas fractions as a function of temperature as discussed in Section 2.4 and Bode et al. (2009). For each halo with M_FoF > 2 × 10¹³ h⁻¹M_☉ and z < 3 (whose gas model is described in Section 2.4.1) we calculate Y₂₀₀ and Y₅₀₀, which are the projected Compton y parameter in a disk of radius R₂₀₀ and R₅₀₀, respectively. Here R₂₀₀ and R₅₀₀ are the radii within which the cluster mean mass density is 200 and 500 times the critical density at the cluster redshift. We then rank the clusters by mass and bin them in mass bins of Δlog(M₂₀₀) = 0.05 and Δlog(M₅₀₀) = 0.1 for M₂₀₀ and M₅₀₀, respectively. The mean of M_Δ as well as the mean and standard deviation on the mean of Y_Δ × E(z)^−2/3 are then calculated for each bin. These values, for M₂₀₀ and Y₂₀₀, are plotted in Figure 12, along with all the individual halos with M₂₀₀ > 2 × 10¹⁴M_☉ in the octant.

Zoom In Zoom Out Reset image size

The binned values are then fit to the power-law relation

$$\begin{eqnarray} &&\frac{Y_{\Delta }}{E(z)^{2/3}} = 10^{\beta } \Big (\frac{M_{\Delta }}{10^{14} M_{\odot }} \Big)^{\alpha }. \end{eqnarray} \tag{ 33 }$$

For M₂₀₀ > 2 × 10¹⁴M_☉, we find for the Y₂₀₀ versus M₂₀₀ relation α = 1.682 ± 0.004 and β = −5.648 ± 0.002, with a reduced χ² of 1.2 for ∼13, 000 clusters divided into 19 bins. This slope is slightly steeper than that expected for the self-similar case. It is also consistent with the Y–M relations derived from hydrodynamical cluster simulations (e.g., White et al. 2002; da Silva et al. 2004; Motl et al. 2005; Nagai 2006; Shaw et al. 2008). The Y–M relation calculated within a disk of radius R₂₅₀₀ is also consistent with the observed Y₂₅₀₀ − M_gas(<R₂₅₀₀) relation of Bonamente et al. (2008), as shown in Bode et al. (2009).

We also find the scatter in this Y–M relation about the best-fit relation. Specifically we calculate

$$\begin{eqnarray} &&\sigma _{YM} = \left[ \frac{\sum ^{N}_{i=1} ({\rm ln} \tilde{Y}_{\Delta } - {\rm ln} Y_{\Delta })_i^2 }{N-2} \right]^{1/2} \end{eqnarray} \tag{ 34 }$$

as discussed in, e.g., Shaw et al. (2008), where $\tilde{Y}_{\Delta }$ is the projected SZ flux measured within a disk of radius R_Δ, Y_Δ is the fitted SZ flux of for a cluster of mass M_Δ using the best-fit parameters for Equation (33), and N is the total number of clusters. We find a scatter for the Y₂₀₀–M₂₀₀ relation of about 14%, and a slightly larger scatter of about 17% for the Y₅₀₀–M₅₀₀ relation. We find a lower scatter of roughly 13% for the Y₅₀₀–M₂₀₀ relation, which has been suggested by Shaw et al. (2008) to have less scatter than the Y₂₀₀–M₂₀₀ relation. We list the best-fit parameters for Equation (33) and the scatter for the Y₂₀₀–M₂₀₀, Y₅₀₀–M₅₀₀, and Y₅₀₀–M₂₀₀ relations in Table 3.

Concerning the scatter given above, we assume hydrostatic equilibrium for the gas; however, we do not assume the same entropy level for all the halos. The kinetic energy of the gas is derived from that of the dark matter, so for a non-virialized, merging halo with a high kinetic energy relative to the cluster potential, the gas will also be given a high energy and thus temperature. Thus merging of halos is taken into account at some level, although a full treatment of the gas physics would give more scatter. Also, the repetition of halos as the periodic box is tiled to fill the octant, even though these halos are often at different redshifts and evolutionary states, may serve to slightly decrease the scatter.

Note that the Y values discussed above are intrinsic to the cluster and do not include any contamination from foreground or background SZ signals along the los. Inclusion of this projection contamination will increase the Y–M scatter somewhat, but for clusters with masses greater than roughly 2 × 10¹⁴M_☉, an optimistic mass limit of current SZ surveys, this projection contamination should be below ∼20% over a wide range of possible σ₈ values (0.7–1.0) and redshifts (z < 2) (Holder et al. 2007; Hallman et al. 2007).

An important question for both determinations of an SZ cluster sample's selection function and measurements of the SZ power spectrum is the correlation and contamination of the SZ cluster signal by infrared and radio galaxies. Observations suggest we should expect some preference for radio and infrared galaxies to reside in galaxy clusters (e.g., Coble et al. 2007; Lin & Mohr 2007; Daddi et al. 2009), and these sources could potentially fill in the SZ decrement at frequencies below the null. This point source contamination of the SZ flux could both affect cluster detection and bias a measured SZ flux–mass scaling relation. If these effects are significant and unaccounted for, then inaccurate cosmological conclusions would result from measured SZ cluster samples. Regarding measurements of the SZ power spectrum and derived values of σ₈, typically a significant number of high flux point sources are masked out of microwave maps before determining the power spectrum. If these sources reside preferentially in massive clusters, a significant portion of the SZ signal could be masked out and a biased determination of σ₈ would arise.

Here we investigate, via these simulations, the fraction of SZ clusters that may be substantially contaminated by radio or infrared galaxy populations. We first add the flux of all galaxies of a given population residing in a given halo. Then, using each galaxy cluster's integrated SZ flux within R₂₀₀ (in units of mJy), we calculate the fraction of halos whose combined flux from resident galaxies equals 20% or more of the halo's SZ flux decrement at 30, 90, and 148 GHz. Such halos will be harder to detect by SZ surveys, and their signal is more likely to be masked out in attempts to measure the SZ power spectrum. We choose a 20% contamination fraction as that is roughly where the uncertainty in the expected SZ flux for a cluster of a given mass becomes dominated by point source contamination as opposed to the intrinsic scatter in the Y–M relation or projection effects.

For the contamination due to radio sources, we find that for halos with M₂₀₀ > 1 × 10¹⁴M_☉, roughly 3% or less of clusters with z < 1.6 have their SZ flux decrement at 148 GHz filled in by more than 20% (see Figure 13). At 90 GHz, this fraction increases slightly, and less than 4% of clusters with z < 1.6 have a contamination of 20% or more. At 30 GHz, one may expect 20% to 30% of clusters with M₂₀₀ ∼ 1 × 10¹⁴M_☉ to be so contaminated, with a decrease in this fraction as the halo mass is increased. These contamination fractions exhibit a weak redshift dependence, with slightly larger fractions at higher redshifts.

Zoom In Zoom Out Reset image size

Infrared galaxies seem to be a significant source of contamination for SZ clusters at 90 and 148 GHz (see Figures 14 and 15). However, these simulations suggest that, at 90 GHz, the fraction of halos with M₂₀₀ > 2.5 × 10¹⁴M_☉ and z < 1.2, that are contaminated at a level of 20% or more, is less than 20%. In all mass bins, the general trend is for the contamination fraction to increase at higher redshifts where more infrared galaxies reside. There is no appreciable contamination from infrared galaxies at 30 GHz. At 148 GHz, less than 20% of halos with M₂₀₀ > 2.5 × 10¹⁴M_☉ and z < 0.8 have their SZ decrements filled in at a level of 50% or more (see Figure 15). It should be noted that the above estimates only include galaxies that reside in halos and does not include contamination from galaxies along a halo's los. In this regard, these fractions should be viewed as lower limits. The one exception to this is for the very high flux infrared galaxies that most likely represent a lensed population of galaxies. These galaxies are placed in massive halos and are the main source of contamination in the most massive bins. These sources should more properly be placed behind massive halos, however, the contamination effect in projection will still be similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Such contamination fractions as discussed above should be factored into sample completeness estimates and the scatter in the Y–M relation in order to use SZ detected clusters for precision cosmology. The models presented here should give a reasonable approximation to the expected levels of contamination from the different galaxy populations. However, these estimates should be cross-checked and made more precise with multi-wavelength observations of SZ clusters.