text:
# Precision era of the kinetic Sunyaev-Zeldovich effect: simulations, analytical models and observations and the power to constrain reionization

## 1 Introduction

Reionization is a critical phase in the evolution of the intergalactic medium (IGM). Recent constraints on the reionization epoch based on Lyman-$`\alpha `$ absorption and CMB appear to be in conflict. The WMAP measurement of the CMB temperature and polarization fluctuation constrains the Thomson optical depth to the last scattering surface $`\tau 0.17\pm 0.04`$. This large optical depth suggests an early reionization at $`z17`$ if the universe is re-ionized suddenly and remains completely ionized after that (Kogut et al., 2003; Spergel et al., 2003). This inference of the reionization redshift and optical depth are very robust: Thomson scattering is polarization dependent, and converts a temperature quadrupole into a linear polarization. The observed $`2\mu K`$ temperature polarization cross correlation are just the product of the optical depths and the local CMB quadrupole, $`Q_{\mathrm{rms},\mathrm{ps}}=20\mu `$K, which is known to high precision. Given the optical depth, and the known baryon content of the universe, one obtains a minimum redshift to accumulate the sufficient total column density of electrons.

But, detection of a Gunn-Peterson absorption trough in a SDSS quasar at redshift $`z=6.28`$ sets a constraint on the reionization redshift at $`z6`$ (Becker et al., 2001). The quasar radiation is seen to be fully absorbed at a Lyman-$`\alpha `$ redshift larger than 6.1. This corresponds to a lower bound of neutral hydrogen column density, corresponding to a neutral fraction greater than a percent. At face value, it would appear to be possible to have a small neutral fraction while the majority of the electrons are ionized. But if ionization occurs through the UV radiation from discrete sources, this radiation must be able to penetrate the IGM to ionize it. It is generally not possible for the radiation to pass through an opaque medium to ionize material behind it. One would generically expect the universe to be either mostly neutral, or fully transparent with neutral fractions less than $`10^3`$. Furthermore, in order to keep the IGM temperature $`T2\times 10^4`$ K at $`z2`$-$`4`$ required by the Lyman-$`\alpha `$ forest observation, a late time reionization must occur after $`z=10`$ (Hui & Haiman, 2003).

Exotic possibilities exist to reconcile a sudden reionization with both WMAP polarization and SDSS quasar data. The ionization could arise from decaying dark matter, which would be uniformly distributed. Or it could have ionized multiple times, and recombined in-between. All these arguments lead to more complicated reionization scenarios.

To make observational progress on this conundrum, more information about the fraction and dynamics of free electrons at $`z10`$ is needed. The primary CMB temperature measurement depends on the reionization history only through $`\tau `$ and thus has an intrinsic degeneracy. A precision CMB E-mode polarization measurement is in principle able to break this degeneracy (Hu & Holder, 2003). But the weak polarization signal makes such measurement unfeasible in the near future. Due to the huge Lyman-$`\alpha `$ optical depth which must be present in any multi-reionization scenarios, the Lyman-$`\alpha `$ forest can only detect the last reionization and makes it an indirect probe of the reionization history. The $`21`$ cm emission and absorption backgrounds in the dark age (Tozzi et al., 2000; Iliev et al., 2002, 2003; Chen & Miralda-Escude, 2003) would constrain the upper limit of the first reionization redshift. The contamination by free-free emission of ionizing halos (Oh & Mack, 2003) may be distinguished by their spectral features. Such a measurement may be possible with a square kilometer array, but its construction is still far in the future.

Free electrons resulting from reionization also scatter off CMB photons through Compton scattering and result in the kinetic Sunyaev Zel’dovich (SZ) effect. Compton scattering keeps information from high redshifts because it does not depend on redshift and is not affected by distance or the expansion of the universe. Since free electrons at high redshifts have higher number density, their contribution to the kinetic SZ effect is larger. So the kinetic SZ effect is a sensitive measure of the reionization process at high redshifts. Next generation CMB experiments such as the Atacama Cosmology Telescope (2003) are making its precision measurement feasible in the foreseeable future. Once we have an accurate model of this effect, the reionization history will be precisely constrained. It is the goal of this paper to build such a model and to quantify the power of the kinetic SZ effect to constrain the reionization history.

Unlike most radiative processes, the physics of the kinetic SZ effect is very simple. In contrast to models of the Lyman-$`\alpha `$ forest, the ionized gas peculiar momentum does not require an accurate understanding of the gas state such as metalicity, temperature and ionization equilibrium. In the linear regime, the exact model is known as the Vishniac effect (Vishniac, 1987). The non-linear regime is intrinsically more difficult and theoretical works in the literature have yet to converge quantitatively (Hu, 2000; da Silva et al., 2001a, b; Gnedin &Jaffe, 2001; Springel et al., 2001; Valageas et al., 2001; Ma & Fry, 2002; Zhang, Pen & Wang, 2002).

Hydrodynamical simulations are best equipped to capture the nonlinear IGM physics. Large box sizes are required to probe large-scale power while high resolution is needed to resolve small-scale, nonlinear structures. Most simulations to date have sacrificed one for the other. But the kinetic SZ effect requires the capturing of both large and small scale structures and puts a strong requirement on the simulation power. Limited computational resources introduce numerical artifacts which must be quantified. In contrast, analytical models can be constructed to span all relevant scales continuously, but they are often ad hoc procedures and must be tested and calibrated against simulations. In this paper, we run Eulerian hydrodynamical simulations with box sizes of 50 and 100 Mpc/$`h`$ to model the evolution of the IGM and measure the gas momentum for various redshifts. Numerical limitations such as finite simulation box size, finite resolution, and non-linearity are studied by running different resolution simulations, ranging from $`256^3`$ up to $`1024^3`$ fluid elements. Self-similar simulations are also used to quantify the numerical limitations. We then build an analytical model covering both linear and nonlinear regimes, motivated by the Vishniac effect and simulations. Our model is free of numerical artifacts and captures the main effect of the non-linearity at the same time.

In $`\mathrm{\S }`$2 we review the formalism for the KSZ effect. In $`\mathrm{\S }`$3 we describe the numerical simulations, including their results and their limitations. Motivated by the gas hydrodynamics and simulations, we build an analytical model in $`\mathrm{\S }`$4. This model is tested against simulations in §5. In $`\mathrm{\S }`$6 we calculate the KSZ power spectrum using our analytical model and in $`\mathrm{\S }`$7 we discuss the potential of the KSZ effect to constrain the reionization history.

## 2 The Kinetic SZ Effect: formalism

The kinetic SZ effect directly probes the ionized electron peculiar momentum $`𝐩(1+\delta _e)𝐯`$ through Compton scattering. The resulting CMB temperature fluctuation is given by Sunyaev & Zel’dovich (1980):

$$\mathrm{\Theta }(\widehat{n})\frac{\mathrm{\Delta }T}{T_{\mathrm{CMB}}}=\chi _e\overline{n}_e\sigma _T\frac{𝐩\widehat{n}}{c}\mathrm{exp}[\tau (z)]adx=\frac{𝐩\widehat{n}}{c}\mathrm{exp}[\tau (z)]d\tau .$$
(1)
Here, $`\overline{n}_e`$ is the mean electron number density and $`\chi _e`$ is the ionization fraction. $`\sigma _T`$ is the Thomson cross section. $`\widehat{n}`$ is the direction on the sky and $`x`$ is the comoving distance. $`\tau (z)=_0^z\chi _e\overline{n}_e\sigma _Tadx`$ is the Thompson optical depth. $`𝐩`$ can be decomposed into a curl-free or gradient part $`𝐩_E`$ satisfying $`\times 𝐩_E=0`$ and a divergence-free or curl part $`𝐩_B`$ satisfying $`𝐩_B=0`$. The gradient term cancels out when integrating along the line of sight and has no contribution to the kinetic SZ effect. Thus the only part which contributes to the integral (1) is $`𝐩_B`$. We define the correlation functions $`\xi _{B,E}(r)𝐩_{B,E}(𝐱)𝐩_{B,E}(𝐱+𝐫)`$ and the corresponding power spectra $`P_{B,E}^2(k)|𝐩_{B,E}(𝐤)|^2`$. Throughout this paper, we use $`\mathrm{\Delta }_B^2(k,z)k^3P_B^2(k)/(2\pi ^2)`$ as the gas momentum curl part power spectrum instead of $`p_B^2(k)`$. $`\mathrm{\Delta }_B^2(k,z)`$ is a detailed description of the IGM state, whose 2D projection is the kinetic SZ power spectrum $`C_l`$ (Eq. 3).

Adopting Limber’s approximation and the fact that $`|\text{p}_B(\text{k})\widehat{n}|^2=\frac{1}{2}P_B^2(k)`$, we obtain the kinetic SZ angular correlation function:

$`w(\theta )`$ $``$ $`\left({\displaystyle \frac{\sigma _T}{c}}\right)^2\mathrm{cos}\theta {\displaystyle _0^{x_{\mathrm{re}}}}𝑑x[a\overline{n}_e\overline{\chi }_e]^2\mathrm{exp}[2\tau (z)]{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{2}}\xi _B(\sqrt{x^2\theta ^2+y^2})𝑑y`$
$``$ $`\left({\displaystyle \frac{\sigma _T}{c}}\right)^2{\displaystyle _0^{x_{\mathrm{re}}}}𝑑x[a\overline{n}_e\overline{\chi }_e]^2\mathrm{exp}[2\tau (z)]{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{2}}\xi _B(\sqrt{x^2\theta ^2+y^2})𝑑y.`$
Here, $`x_{\mathrm{re}}`$ is the comoving distance to the reionization epoch. We have set $`\chi _e=\overline{\chi }_e`$, namely, omitted the patchy reionization effect222We will discuss the patchy reionization in §7. following the suggestion of both simulations (Gnedin &Jaffe, 2001) and theory (Valageas et al., 2001). The last approximation introduces an error less than $`0.2\%`$ for $`\theta 1^0`$. Thus, the usual expression of the Limber’s equation (Peacock, 1999) still holds:

$$C_l=\frac{16\pi ^2}{(2l+1)^3}[\frac{\overline{n}_e(0)\sigma _T}{c}]^2\times _o^{z_{\mathrm{re}}}(1+z)^4\overline{\chi }_e^2\frac{1}{2}\mathrm{\Delta }_B^2(k,z)|_{k=l/x}\mathrm{exp}[2\tau (z)]x(z)\frac{dx(z)}{dz}dz.$$
(3)
$`𝐩_B`$ is zero in first order linear perturbation theory. Thus, the kinetic SZ effect is an intrinsically non-linear process which must be investigated by either hydro simulations or non-linear theories such as high order perturbation theory. After the pioneering work of the Vishniac effect (Vishniac, 1987), which is a robust prediction of the kinetic SZ effect in the linear regime, various authors (Hu, 2000; da Silva et al., 2001a, b; Gnedin &Jaffe, 2001; Springel et al., 2001; Valageas et al., 2001; Ma & Fry, 2002; Zhang, Pen & Wang, 2002) have studied the kinetic SZ effect in the non-linear regime with analytical models and simulations, but there are still some unclear key issues. Three key issues we will address in this paper are (1) How well do simulations capture $`\mathrm{\Delta }_B^2`$ and thus the large scale power spectrum of the kinetic SZ effect at large and small scales? $`\mathrm{\Delta }_B^2`$ has contributions from both large and small scales. Simulations have cutoffs at both scales due to the finite simulation box size and resolution, respectively. These cutoffs affect simulation results and need to be quantified. Then, how to quantify and overcome these effects to produce reliable results? (2) What is the effect of the non-linearity to $`\mathrm{\Delta }_B^2`$? How to model this effect? (3) How observable is the effect, and to what accuracy can it be used to model measured reionization by upcoming experiments? We will address these issues using simulations (§3) and our analytical model (§4).

## 3 Hydrodynamic simulations

We ran cosmological hydrodynamical simulations using a new Eulerian cosmological hydro code (Trac & Pen, 2003a, b). This Eulerian code (hereafter TP) is based on the finite-volume, flux-conservative total variation diminishing (TVD) scheme that provides high-order accuracy and high-resolution capturing of shocks. The hydrodynamics of the gas is simulated by solving the Euler system of conservation equations for mass, momentum, and energy on a fixed Cartesian grid. The gravitational evolution of the dark matter is simulated using a cloud-in-cell particle-mesh (PM) scheme (Hockney & Eastwood, 1988).

The robustness of the TP code has been tested by comparing the evolution of the dark matter and gas density power spectra from the simulations with the fitting formula of Smith et al. (2002). We also performed a code comparison by running the same initial conditions using the MMH code (Pen, 1998), which combines the shock capturing abilities of Eulerian schemes with the high dynamic range in density achieved by Lagrangian schemes. Power spectra are computed using FFTs. We find good agreement at all relevant scales and redshifts for both comparisons.

The KSZ effect is a non-local problem and has contributions from both large-scale and small-scale power. This makes it a challenging problem to model numerically because high resolution is needed at all scales. Eulerian schemes provide high mass resolution at all scales, unlike Lagrangian schemes based on SPH which achieve high dynamic range in high-density regions but do poorly in low-density regions. Eulerian schemes are ideal for simulating the evolution of the IGM to model the thermal and kinetic SZ effects and the Lyman alpha forest, because of their high speed, superior mass resolution, shock-capturing abilities. Furthermore, Eulerian algorithms are computationally very fast and memory friendly, allowing one to optimally use available computational resources.

We ran a total of six simulations with the best fit WMAP-alone cosmology $`\mathrm{\Omega }_m=0.268`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.752`$, $`\mathrm{\Omega }_b=0.044`$, $`h=0.71`$, and $`\sigma _8=0.84`$ (Spergel et al., 2003). Two box sizes of 50 and 100 Mpc/$`h`$ are used and for each we perform the simulations on fixed grids with $`256^3`$, $`512^3`$, and $`1024^3`$ cells. The ratio of dark matter particles to fluid elements is 1:8 in all six $`\mathrm{\Lambda }`$CDM simulations. In the highest resolution simulation, we achieve a spatial resolution of $`\mathrm{\Delta }x50`$ kpc/$`h`$ and a dark matter particle mass resolution of $`\mathrm{\Delta }m7\times 10^7M_{}`$. The initial conditions are generated by sampling from an initial power spectrum computed using CMBFAST (Seljak & Zaldarriaga, 1996). The simulations are started at a redshift of $`z=100`$ and evolved down to $`z=1`$, with data outputs at $`z=`$ 20, 10, 6, 3, and 1. The KSZ effect, for angular scales $`l1000`$, has contributions from electron scattering primarily from the $`z1`$ IGM and therefore, we do not push the simulations below $`z=1`$ for the purposes of this paper. Furthermore, the final prediction of the kinetic SZ effect is calculated using our analytical model, instead of using these simulations directly, where the contribution from $`z<1`$ IGM is included. The highest resolution simulation takes approximately 700 time-steps to evolve from $`z=100`$ down to $`z=1`$. On a GS320 Compaq Alpha server with 32 cpus and total theoretical peak speed of 32 Gflops, the run takes approximately 2 days.

The peculiar gas momentum is decomposed into the rotational (B) and irrotational terms (E) in Fourier space. In Fig. 1 we plot the simulated dimensionless $`\mathrm{\Delta }_B(k,z)k/H(z)`$ at $`z=1`$ and $`10`$ and compare them to the linear Vishniac effect, which is a robust prediction in the linear regime (§4). $`H(z)`$ is the Hubble parameter. $`\mathrm{\Delta }_B(k,z)k/H(z)`$ is the velocity of a wave expressed in units of wave lengths per expansion time. It reaches $`1`$ when the overdensity reaches the mildly nonlinear value $`\delta 1`$. We note that our convention differs from the literature that we plot the root-mean-square standard deviation instead of a variance.

Three key conclusions can be drawn from our simulations. (1) Finite simulation box size causes $`\mathrm{\Delta }_B^2(k,z)`$ to lose power at large scales. At redshift $`z=10`$, where linear theory is a good approximation, we would expect the simulation results to agree well with the Vishniac prediction. But it turns out that they agree only at a very limited $`k`$ range. More counter-intuitively, the discrepancy diverges toward large scales. The simulated $`\mathrm{\Delta }_B^2(k,z)`$ loses power toward large scales. 50 Mpc/$`h`$ simulations lose more power at large scales than 100 Mpc/$`h`$ simulations at the same scales. This suggests the effect of the finite simulation box size. Simulations cut off density fluctuation modes larger than half simulation box size and thus affect $`\mathrm{\Delta }_B^2`$ through its non-local dependence on the density and velocity fields. (2) Finite simulation resolution causes $`\mathrm{\Delta }_B^2(k,z)`$ to lose power at small scales (less than $`10`$-$`20`$ grid separation). $`\mathrm{\Delta }_B^2(k,z)`$ truncates at small scales and the truncating scale $`k_{\mathrm{trunc}}`$ keeps increasing with respect to the simulation resolution. Gas pressure is expected to decrease power at very small scales. If this truncation is caused by the gas pressure, we would expect $`k_{\mathrm{trunc}}`$ to decrease with the simulation resolution since higher resolution results in higher gas pressure. Since we observe opposite behavior, we believe that this truncation is caused by the simulation resolution, which predicts that $`k_{\mathrm{trunc}}`$ to be roughly proportional to the resolution, as apparent in figure. 1. (3) The non-linearity increases $`\mathrm{\Delta }_B^2(k,z)`$ with respect to the linear prediction. This is visible in the top panel of figure. 1.

We further investigate these issues by self similar simulations. We ran two self similar $`512^3`$ simulations with $`n=1`$ and $`n=2`$, respectively. $`\mathrm{\Omega }_m=1.0`$. $`\mathrm{\Omega }_b=0.164`$ is adopted such that $`\mathrm{\Omega }_b/\mathrm{\Omega }_m`$ is the same as the WMAP cosmology . The effects of finite simulation box size and finite simulation resolution scale as expected. The non-linear effect turns out to be more subtle. For $`n=1`$, simulated $`\mathrm{\Delta }_B^2`$ is smaller than the Vishniac prediction in the non-linear regime. The non-linearity decreases $`\mathrm{\Delta }_B^2`$. For $`n=2`$, simulated $`\mathrm{\Delta }_B^2`$ is larger than the Vishniac prediction. Thus the non-linearity increases $`\mathrm{\Delta }_B^2`$. To quantify simulation artifacts and the non-linearity effect, we build an analytical model in §4. Therefore we postpone the discussion and explanation of the self similar simulations until in §5.

## 4 The analytical model

In this section, we will extend the Vishniac result to the grid space to model the simulation effect (§4.1) and to the non-linear regime to model the non-linearity effect (§4.2).

### 4.1 Linear regime

In the linear regime, the motion of IGM is determined by the gravitational potential $`\varphi `$. Thus $`𝐯\varphi `$ and is curl-free. Since $`^2\varphi \delta `$, $`𝐯(𝐤)\delta _k𝐤/k^2`$. From this relation, $`\mathrm{\Delta }_B^2`$ in the linear regime was originally calculated by Vishniac (1987). With the correction of the gas window function $`W_g`$, we obtain

$$\mathrm{\Delta }_B^2(k)=\frac{k^3}{2\pi ^2}\frac{1}{2}\frac{(a\dot{D}/D)^2}{(2\pi )^3}_0^{\mathrm{}}P_L(k_2,z)P_L(|\text{k}\text{k}_2|,z)\left[W_g(|\text{k}\text{k}_2|)\beta (\text{k,k}_2)+W_g(k_2)\beta (\text{k},\text{k}\text{k}_2)\right]^2\mathrm{d}^3k_2.$$
(4)
$`D`$ is the linear growth factor of the dark matter over-density, $`\dot{D}`$ is its time derivative and $`a`$ is the scale factor of the universe. $`P_L(k,z)`$ is the linear dark matter power spectrum. Gas pressure originated from various non-gravitational heating processes smooths out the gas density at small scales. This effect can be modeled by a window function $`W_g(k)`$ such that $`P_g(k)=W_g^2(k)P_{\mathrm{DM}}(k)`$, which can be choose as a Gaussian window function (see e.g. Hui & Gnedin (1997); Gnedin (2000); Gnedin et al. (2003)). In an adiabatically evolving universe, gas traces dark matter down to very small scales and $`W_g`$ would be effectively unity. One can model the simulation resolution by introducing an effective grid smoothing window $`W_g`$. This procedure works well but the modeling of $`W_g`$ is arbitrary. For simplicity, we omit this gas window effect. Since the gas peculiar velocity $`𝐯_𝐠`$ is mainly determined by the dark matter gravitational potential field, we have assumed that $`𝐯_𝐠`$ traces the dark matter peculiar velocity at all scales. The kernel $`\beta (𝐤,𝐤_\mathrm{𝟐})=(𝐤_\mathrm{𝟐}\frac{𝐤(𝐤𝐤_\mathrm{𝟐})}{𝐤^\mathrm{𝟐}})/𝐤_\mathrm{𝟐}^\mathrm{𝟐}`$. This form originates from the gradient operators in $`𝐯\varphi `$ and $`^2\varphi \delta `$. Eq. (4) explicitly preserves the symmetry between $`𝐤_1𝐤𝐤_2`$ and $`𝐤_2`$.

Eq. (4) is exact in the linear regime and an ideal simulation should be able to reproduce these results. But we notice from this equation that $`\mathrm{\Delta }_B^2`$ has contributions from both large and small scales. For a real simulation, the simulation box size and resolution cut off contributions from these scales and the simulation result may deviate from Eq. (4). To address this issue and to estimate the limitation of simulations, we calculate $`\mathrm{\Delta }_B^2`$ on discrete grids. For a simulation with box size $`L`$, resolution $`N`$ and the periodic boundary condition, the Fourier component of the gradient operator $``$ is $`\stackrel{~}{}=i\mathrm{sin}(𝐤\mathrm{\Delta }L)/\mathrm{\Delta }Li\{\mathrm{sin}(k_x\mathrm{\Delta }L)/\mathrm{\Delta }L,\mathrm{sin}(k_y\mathrm{\Delta }L)/\mathrm{\Delta }L,\mathrm{sin}(k_z\mathrm{\Delta }L)/\mathrm{\Delta }L\}`$. Here, $`\mathrm{\Delta }LL/N`$. Comparing to $`\stackrel{~}{}=i𝐤`$ in the continuous case, we obtain the grid version of $`\mathrm{\Delta }_B^2`$:

$$\mathrm{\Delta }_B^2(k)=\frac{k^3}{2\pi ^2}\frac{1}{2}(a\dot{D}/D)^2L^3\underset{\text{k}_2}{}P_L(k_2,z)P_L(|(\text{k}\text{k}_2)_p|,z)\left[W_g(|(\text{k}\text{k}_2)_p|)\beta _G(\text{k,k}_2)+W_g(k_2)\beta (\text{k},(\text{k}\text{k}_2)_p)\right]^2.$$
(5)
Here, $`\beta _G(𝐤,𝐤_\mathrm{𝟐})=\beta (𝐤\mathrm{sin}(𝐤\mathrm{\Delta }L)/\mathrm{\Delta }L,𝐤_\mathrm{𝟐}\mathrm{sin}(𝐤_\mathrm{𝟐}\mathrm{\Delta }L)/\mathrm{\Delta }L)`$, $`k_{x,y,z},k_{2x,y,z}=2\pi /L\times (N/2,N/21,\mathrm{},N/21)`$. $`(\text{k}\text{k}_2)_p`$ means $`\text{k}\text{k}_2`$ under the periodic boundary condition. Namely, $`\text{k}\text{k}_2`$ needs to be converted to the same range as $`\text{k}_2`$ by the periodic boundary condition. This grid version represents simulation discretization effect at the first order approximation.

### 4.2 Non-linear regime

In the non-linear regime, the perturbation theory does not hold and Eq. (4) & (5) break down. Hu (2000); Ma & Fry (2002) argue that the velocity field is less non-linear than the density field and one may substitute $`P_L(k,z)`$ in Eq. (4) introduced by the density field with the corresponding non-linear density power spectrum $`P_{NL}(k,z)`$ and keep the other $`P_L`$ introduced by the velocity field unchanged. But the generation of curl in the velocity field by shell crossing in the non-linear regime weakens this argument. In the presence of curl, $`\times 𝐩_B=\delta \times 𝐯_E+\delta \times 𝐯_B+(1+\delta )\times 𝐯_B`$. The first term in the right side of this equation is well described by the linear Vishniac effect and may be well described by the method of Hu (2000); Ma & Fry (2002) in the non-linear regime. But the remaining terms bring non-negligible contribution in the non-linear regimes. In our simulations, we find that in the highly non-linear regime, velocity field reaches equi-partition and $`\mathrm{\Delta }_{𝐯_B}^2=2\mathrm{\Delta }_{𝐯_E}^2`$. This suggests that the contribution from $`𝐯_B`$ may be dominant in the non-linear regime. Thus, the method of Hu (2000); Ma & Fry (2002) underestimates $`\mathrm{\Delta }_B^2`$. The generation of $`𝐯_B`$ through shell crossing involves complicated gas hydrodynamics and makes the analytical understanding of $`𝐯_B`$ intractable. Its statistics in principle can only be probed by simulations. But since the generation of $`𝐯_B`$ is caused by the non-linearity, to model its effect, one may phenomenologically substitute the linear density power spectrum in Eq. (4) introduced by the velocity field with the corresponding non-linear density power spectrum. This procedure at least qualitatively captures the effect of $`𝐯_B`$ generation. Thus we obtain the final continuous form of $`\mathrm{\Delta }_B^2(k,z)`$:

$$\mathrm{\Delta }_B^2(k,z)=\frac{k^3}{2\pi ^2}\frac{1}{2}\frac{(a\dot{D}/D)^2}{(2\pi )^3}_0^{\mathrm{}}P_{NL}(k_2,z)P_{NL}(|\text{k}\text{k}_2|,z)\left[W_g(|\text{k}\text{k}_2|)\beta (\text{k,k}_2)+W_g(k_2)\beta (\text{k},\text{k}\text{k}_2)\right]^2\mathrm{d}^3k_2$$
(6)
and the corresponding grid form

$$\mathrm{\Delta }_B^2(k)=\frac{k^3}{2\pi ^2}\frac{1}{2}(a\dot{D}/D)^2L^3\underset{\text{k}_2}{}P_{NL}(k_2,z)P_{NL}(|(\text{k}\text{k}_2)_p|,z)[W_g(|(\text{k}\text{k}_2)_p|)\beta _G(\text{k,k}_2)++W_g(k_2)\beta (\text{k},(\text{k}\text{k}_2)_p)]^2.$$
(7)
In our model, the nonlinearity effect is explicitly quadratic and larger than that of Hu (2000); Ma & Fry (2002). We will adopt the code of Smith et al. (2002) to calculate $`P_{NL}(k,z)`$.

In §5, we will test our model against simulations. Our logic is, if the grid version of our model, which captures key simulation artifacts, gives a good description of the simulation, then, its asymptotic case, namely the continuous version of our model would represent an ideal simulation, be free of these simulation artifacts and describe the real $`\mathrm{\Delta }_B^2`$.

## 5 Comparison between analytical model and simulations

We first compare the grid version prediction of our model Eq. (7) with simulations. Self similar simulations are ideal to test the effects of box size, resolution and non-linearity due to the scaling relation between redshifts and scales (See e.g. Peebles (1980)). The grid prediction of our model fits simulations very well at all reliable scales. The cut off due to the finite simulation box size can be seen from the large scale behavior of $`\mathrm{\Delta }_B^2`$ and the resolution effect can be seen from the small scale behavior of $`\mathrm{\Delta }_B^2`$ (Fig. 2). The dependence of the non-linearity effect on the initial density power index $`n`$ is naturally explained by our model. The effect of the non-linearity to the density power spectrum depends on $`n`$. For $`n=1`$, the non-linearity decreases the non-linear density power spectrum with respect to the linear density power spectrum. Our model then predicts a decreasing in $`\mathrm{\Delta }_B^2`$ with respect to the Vishniac prediction by the non-linearity. For $`n=2`$, it is the opposite case. For the WMAP cosmology, the effective power index in simulation scales is around $`2`$ and an amplification of $`\mathrm{\Delta }_B^2`$ by the non-linearity is expected, as shown in Fig. 1.

Since our ultimate goal is to provide an accurate description of $`\mathrm{\Delta }_B^2`$ for our real universe, we test our model in detail against the $`1024^3`$ simulations with WMAP cosmology. In order to do that, we need to estimate the variance of $`\mathrm{\Delta }_B^2(k)`$. In the linear regime $`\delta `$ is Gaussian distributed and $`p_B(\text{k})\delta \delta \frac{\text{k}}{k}`$, we obtain $`p_B^4(\text{k})=\frac{105}{9}p_B^2(\text{k})^2`$ and estimate the uncertainty of the simulated $`\mathrm{\Delta }_B^2`$ by this relation. Our result (Fig. 3) shows that the grid calculation agrees with simulated $`\mathrm{\Delta }_B^2`$ at all output redshifts and the full $`k`$ range that simulations can reliably probe. Our model passed the test of all self similar and WMAP simulations we have run and shows a wider range of applicability than all previous analytical models. Thus we believe that the grid version of our model is a suitable tool to describe simulated $`\mathrm{\Delta }_B^2`$. Its continuous version (Eq. 6) then should be able to describe an ideal simulation with infinite box size and resolution and produce the real $`\mathrm{\Delta }_B^2`$.

Comparing our analytical prediction and simulations, we concluded that the simulated $`\mathrm{\Delta }_B^2`$ loses power at large scales, gains power at small scales due to non-linearity until approaching the resolution limit. Thus it has a steeper slope. We address several simulation limitations here. (1) The finite simulation box size causes the simulated $`\mathrm{\Delta }_B^2`$ to lose power at large scales. From Eq. (7), $`\mathrm{\Delta }_B^2`$ has contributions from all scales. The finite box size cuts off the density fluctuation modes with size larger than half box size. It is these missing modes causing the simulated $`\mathrm{\Delta }_B^2`$ to lose power at large scales. The smaller the box size, the more $`\mathrm{\Delta }_B^2`$ power is lost, as can be seen from Fig. 3. (2) The finite resolution causes $`\mathrm{\Delta }_B^2`$ to lose power at less than about $`10`$-$`20`$ grid separation. (3) The periodic boundary condition causes the simulated $`\mathrm{\Delta }_B^2(k)`$ to gain power. This can be seen by the following argument. For the cosmology we adopt, $`P(k)`$ is a decreasing function for $`k>0.02`$ Mpc/h. Comparing Eq. (6) and Eq. (7), the periodic boundary condition produces a smaller $`|\text{k}\text{k}_2|`$ and thus a larger $`P(|\text{k}\text{k}_2|)`$ for those $`\text{k}\text{k}_2`$ which are beyond the simulation k range. The case is the same for the gas window function and is similar for the kernel function $`\beta (\text{k,k}_2)`$. But this effect is negligible comparing to the first two.

Our modeling of the non-linearity effect works surprisingly well. For WMAP cosmology, the non-linearity does increase $`\mathrm{\Delta }_B^2`$ significantly (fig. 3). As expected, the non-linearity effect is larger than what Hu (2000); Ma & Fry (2002) predicted since their models omit the contribution of $`𝐯_B`$. Our model passes the test of all of our simulations and suggests that it not only captures the contribution of $`𝐯_B`$ qualitatively but also quantitatively. Though it is hard to explain why it works well quantitatively, we believe that it gives a reasonable description of $`\mathrm{\Delta }_B^2`$ in the non-linear regime.

Our model could be understood from the equi-partition of $`𝐩`$ in the non-linear regime, which states $`\mathrm{\Delta }_B^2=2\mathrm{\Delta }_E^2`$, as found in our simulations, and the mass conservation equation, which states $`\dot{\delta }+𝐩_E/a=0`$. Defining the cross correlation coefficient $`r\mathrm{\Delta }_{\delta \dot{\delta }}^2/\sqrt{\mathrm{\Delta }^2\mathrm{\Delta }_{\dot{\delta }}^2}`$ between $`\delta `$ and $`\dot{\delta }`$, one then has

$$\mathrm{\Delta }_B^2=\frac{1}{2}\frac{a^2}{k^2}\frac{\dot{\mathrm{\Delta }^2(k)}^2}{r^2\mathrm{\Delta }^2(k)}.$$
(8)
$`r`$ in the linear theory and the stable clustering regime can be calculated analytically. But since possibly due to the ongoing merger process, the stable clustering may not hold even in highly non-linear regime (Smith et al., 2002), for most non-linear regimes, $`r`$ is hard to estimate, even in an interpolation way. In our simulations, we find that in highly non-linear regime, $`\delta `$ and $`\dot{\delta }`$ becomes weakly correlated . $`r^2`$ roughly scales as $`1/\mathrm{\Delta }^2`$ (Fig. 4). This scaling gives $`\mathrm{\Delta }_B^2\mathrm{\Delta }^4`$, which is consistent with our previous model (Eq. 6 and 7).

Simulation limitations put strong requirements on the simulation power. For example, from Eq. (7), we estimate that, in order to produce a reliable $`\mathrm{\Delta }_B^2`$ at $`0.1`$ h/Mpc $`<k<10`$ h/Mpc, a hydro simulation with a box size $`L400`$ Mpc/h and resolution $`N>4096`$ is required, which is not realistic at present. But the continuous version of our analytical model (Eq. 6) is free of these simulation limitations and is able to produce a reliable $`\mathrm{\Delta }_B^2`$ at all $`k`$ and $`z`$ range relevant to the kinetic SZ effect. Hereafter, we will adopt Eq. (6) to calculate the kinetic SZ power spectrum.

## 6 The kinetic Sunyaev Zel’dovich effect

The kinetic SZ power spectrum $`C_l`$ is very sensitive to the reionization history. To see this, we estimate its dependence on the reionization redshift. For simplicity, we assume that the universe is reionized at once at $`z_{\mathrm{reion}}`$ and remains completely ionized after that. Since $`\tau 1`$, we omit its effect temporarily for the estimation. First, we consider the linear case. Combining Eq. (3) & (4), we obtain

$`C_l(z_{\mathrm{reion}})`$ $``$ $`{\displaystyle ^{z_{\mathrm{reion}}}}(1+z)^4(a\dot{D}D)^2(\chi (z))^{1n_B}{\displaystyle \frac{\mathrm{d}\chi }{\mathrm{d}z}}𝑑z{\displaystyle ^{z_{\mathrm{reion}}}}({\displaystyle \frac{dD}{da}}D)^2\sqrt{(\mathrm{\Omega }_0(1+z)^3+\mathrm{\Omega }_v)}(\chi (z))^{1n_B}𝑑z`$
$``$ $`{\displaystyle ^{z_{\mathrm{reion}}}}(1+z)^{(1n_B)\alpha 0.5}𝑑z(1+z_{\mathrm{reion}})^{(1n_B)\alpha +0.5}\mathrm{for}z1`$
Here, $`n_B`$ is the power index of $`\mathrm{\Delta }_B^2(k,z)`$ at $`k=l/\chi (z)`$. For $`l10^3`$ and $`z>1`$, $`n_B1`$. We approximate $`\chi (z)(1+z)^\alpha `$. For $`z2`$, $`0<\alpha <0.5`$. For $`z10`$, $`0<\alpha <0.1`$. When $`z2`$, $`Da`$ and $`\frac{dD}{da}1`$. Then, for $`z>2`$, $`C_l`$ diverges toward high $`z_{\mathrm{reion}}`$ with a scaling relation $`C_l(1+z_{\mathrm{reion}})^{0.5}`$. The non-linearity increases the contribution from low redshifts and makes the $`C_l`$ redshift dependence weaker. The damping caused by the optical depth $`\tau `$ causes further suppression. We show $`C_l`$ calculated for various $`z_{\mathrm{reion}}`$ by Eq. (3) & (6) in Fig. 5. The $`C_l`$ behavior has three distinct regions. (1) $`l1000`$. $`C_l`$ in this region is mainly contributed by low redshift linear regions. Its dependence on $`z_{\mathrm{reion}}`$ is weak. (2) $`l10000`$. $`C_l`$ in this region is mainly contributed by low redshift highly non-linear regions. Its dependence on $`z_{\mathrm{reion}}`$ is also weak. (3) $`1000l10000`$. $`C_l`$ in this region has large contributions from the high redshift universe and thus has the highest sensitivity to $`z_{\mathrm{reion}}`$. As we will discuss in §7, this region is also observationally accessible in the near future. For $`6<z_{\mathrm{reion}}<20`$ and $`3000<l<9000`$, our result can be fitted by

$$l^2C_l/(2\pi )0.91\times 10^{12}[(1+z_{\mathrm{reion}})/10]^{0.34}(l/5000)^{0.230.015(z_{\mathrm{reion}}9)}.$$
(10)
The error introduced by this fitting formula is less than $`2\%`$. The corresponding temperature fluctuation $`\mathrm{\Delta }T`$ is then

$$\mathrm{\Delta }T_l2.60[(1+z_{\mathrm{reion}})/10]^{0.17}(l/5000)^{0.1150.0075(z_{\mathrm{reion}}9)}\mu \mathrm{K}.$$
(11)
For the WMAP favored $`z_{\mathrm{reion}}=16.5`$, we show $`C_l`$ predicted by our model and by the corresponding Vishniac effect in fig 6. The non-linearity increases the power spectrum at small angular scales. At $`l4000`$, this amplification reaches $`2`$. The non-linearity then produces a broad peak in the power spectrum, which extends from $`l3000`$ to $`l10^5`$ with an amplitude $`\mathrm{\Delta }T2.7\mu `$K. The large $`l`$ behavior will be certainly affected by the gas window function $`W_g(k)`$. Its effect can be easily incorporated, as discussed in §4.1. Since it mainly affects the $`l>10^4`$ region, where it is hard to observe, we will not discuss its effect here.

## 7 The power for the kinetic SZ effect to constrain the reionization history

For $`z_{\mathrm{reion}}7`$, as favored by the Lyman-$`\alpha `$ absorption observation versus $`z_{\mathrm{reion}}17`$, as inferred from WMAP, $`C_l`$ would differ by $`45\%`$. This sensitive dependence of $`C_l`$ on the reionization history opens a window for future CMB observations to constrain the reionization history and break the degeneracy met in the primary CMB argument. The power of the kinetic SZ effect to constrain the reionization history depends on the accuracy of the kinetic SZ effect measurement. In this section, we take the Atacama cosmology telescope (ACT) as our target to address this issue.

ACT will operate at three bands, $`145`$ Ghz, $`225`$ Ghz and $`265`$ Ghz, measure over $`100`$ square degree sky with noise error $`\sigma _N2\mu `$K per $`1.7^{^{}}\times 1.7^{^{}}`$ pixel. At the $`225`$ Ghz band, the thermal SZ signal is negligible and we take this band to estimate the sensitivity of ACT to measure the kinetic SZ effect. The noise power spectrum is given by $`C_N=4\pi f_{\mathrm{pix}}\sigma _N^2`$. Here, $`f_{\mathrm{pix}}`$ is the fractional sky coverage of each pixel. The error in the kinetic SZ power spectrum measurement comes from three sources: instrumental noise, primary CMB and cosmic variance. If the primary CMB power spectrum and the noise spectrum are known and can be deducted, we have then

$$\frac{\mathrm{\Delta }C_l}{C_l}=\sqrt{\frac{(\eta +2)+2(C_N/C_l/W_N^2)^2+2(C_{CMB}/C_l)^2}{(2l+1)\mathrm{\Delta }lf_{\mathrm{sky}}}}.$$
(12)
$`W_N(l)`$ is the Fourier transform of the natural beam function. For the $`225`$ Ghz band, the ACT resolution is $`1.1^{^{}}`$. We then approximate the beam function as a top-hat function with an effective radius $`\sqrt{1.1^{^{}}\times 1.1^{^{}}/\pi }`$. We normalize it to have $`W_N(l1)1`$. $`f_{\mathrm{sky}}`$ is the fractional sky coverage of the whole survey area, which we adopt as $`100`$ square degree. $`\eta `$ is a measure of the non-Gaussianity, which will be taken as $`105/93=26/3`$, as in the 3D case obtained from the second order perturbation theory (§4.1). Due to the projection, the 2D kinetic SZ non-Gaussianity would be weaker and the error obtained above should serve as an upper limit.

We consider three reionization scenarios: (A) the universe reionized once at $`z_{\mathrm{reion}}=7`$, (B) the universe reionized completely once at $`z_{\mathrm{reion}}=16.5`$, (C) the universe reionized completely at $`z=21`$ and became neutral at $`z=13`$ and reionized again at $`z=7`$. The last two has roughly the same optical depth $`\tau =0.174`$ to the last scattering surface and are degenerate in the primary CMB measurement. In Fig. 7, we forecast the accuracy for ACT to measure the kinetic SZ $`C_l`$ assuming the reionization scenario B. $`\mathrm{\Delta }l=l/10`$ is adopted. At $`l2500`$, primary CMB dominates and prohibits the extraction of the kinetic SZ signal. At $`l10^4`$, the instrumental noise dominates. At $`l6000`$, ACT will reach its highest sensitivity with a $`3\%`$ accuracy in the kinetic SZ $`C_l`$ measurement. Since the kinetic SZ power spectrum is almost flat at $`l>3000`$, one can bin the data points ($`10`$ points, fig. 7) at $`3000<l<8000`$ and do not lose any crucial information. Binning in this way, the ACT-measured $`C_l`$ accuracy will reach better than $`1\%`$.

The accuracy of inferred reionization history depends on both the precision of KSZ power spectrum measurement and theoretical predictions. Given perfect theoretical predictions, for the single reionization scenario, the uncertainty in the reionization redshift $`\mathrm{\Delta }z_{\mathrm{reion}}/z_{\mathrm{reion}}\mathrm{\Delta }C_l/C_l/0.343\%`$. The two CMB-degenerate reionization scenario B & C can be distinguished with $`>10\sigma `$ confidence333We assume different bins are uncorrelated. This assumption is exact for a whole sky survey. For a survey with limited sky coverage, this assumption only roughly holds in our case. For two successive data points of Fig. 7 at $`3000<l<8000`$, the corresponding spacial separation at $`z10`$ is $`1`$ Mpc/h, which is about the same as the density correlation length at $`z=10`$. So, one can safely neglect cross correlations in different bins.. But uncertainties in the KSZ signal extraction and theoretical prediction may degrade the power of ACT to measure the kinetic SZ effect.

One is the lensing effect. Since lensing changes CMB photon positions but not CMB temperature, it results in a smoothing of the CMB power spectrum. Because the primary CMB power spectrum drops quickly after $`l2000`$, lensing amplifies it by large factors at small angular scales (Fig. 7). In our estimation, we have assumed that this lensed CMB power spectrum can be predicted and subtracted. Since the lensing effect at these scales involves non-linear dark matter clustering, the precise prediction of its power spectrum could be difficult. If we take this possibility, the residual CMB (lensed subtracted unlensed) would prohibit the measurement of KSZ at $`l<4000`$ and increase the error of KSZ $`C_l`$ measurement at $`l[4000,6000]`$ to be larger than $`10\%`$. But its effect at $`l>6000`$ is small. Binning the $`4`$ data points at $`7000<l<9000`$ would still give a better than $`2\%`$ accuracy. Furthermore, the future lensing survey such as the CFHT legacy survey (2003) will measure the dark matter distribution to $`1\%`$ accuracy and allow a precise prediction of the lensed CMB power spectrum. So, the lensing effect of the primary CMB could be well handled. Lensing also changes the KSZ power spectrum. But since the KSZ power spectrum is almost flat, this lensing effect is negligible.

In our discussion, we have assumed the thermal SZ effect to be negligible at $`225`$ Ghz. This is true in the non relativistic limit. The relativistic correction (see e.g. Dolgov et al. (2001)) produces a temperature decrement at this band and its effect on ACT measurement of cluster SZ effect has been discussed (Aghanim et al., 2002). This relativistic correction is of the order $`k_BT_e/(m_ec^2)T_e/500`$keV. Though it is non-trivial for the measurement of high temperature cluster SZ effect, it is negligible for a blank sky SZ power spectrum measurement. The mean electron temperature at present due to gravitational heating is expected to be less than $`0.5`$ keV (Pen, 1999; Zhang & Pen, 2001). COBE/FIRAS puts an upper limit $`1.5\times 10^5`$ (Fixsen et al., 1996) on the mean ’y’ parameter and basically excludes dramatic feedback with $`T1`$ keV. One then expects the mean electron temperature at present to be $`1`$ keV. Since the blank sky thermal SZ power spectrum at $`l>4000`$ is mainly contributed by these low temperature IGM gas at $`z0.5`$ (Zhang & Pen, 2001), we expect that the relativistic correction is much less than $`1\%`$ and is thus negligible.

We have not considered the contamination of unresolved IR and radio sources, whose effects are estimated to be non-negligible (White & Majumdar, 2003). At small angular scales, these sources are present as Poisson noise in the KSZ measurement. From the large $`l`$ power spectrum behavior, the power spectrum of the sum of instrumental noise and this Poisson noise can in principle be measured and subtracted. Multi-frequency observations also help to clean them out.

We estimate the dependence of the KSZ power spectrum on cosmological parameters. One can easily work out that $`\mathrm{\Delta }_B^2(k,z)\beta ^2H^2(z)\sigma _8^{46}`$ ($`4`$ in the linear regime and $`6`$ in the highly nonlinear regime where stable clustering holds). Here, $`\beta adD/da/D\mathrm{\Omega }_m^{0.6}(z)`$. In our interested $`3000l10000`$, $`C_l`$ contributions mainly come from $`z>1`$, where $`\beta 1`$, $`H(z)\mathrm{\Omega }_0^{1/2}`$ and the comoving distance $`x`$ is roughly $`\mathrm{\Omega }_0^\gamma `$ with $`0<\gamma <1/2`$. Here, $`\mathrm{\Omega }_0\mathrm{\Omega }_m(z=0)`$. Applying the above relation to Eq. 3, we find that $`C_l(\mathrm{\Omega }_bh)^2\mathrm{\Omega }_0^{1/2+\gamma (n_B1)}\sigma _8^{46}(\mathrm{\Omega }_bh)^2\mathrm{\Omega }_0^{1/2}\sigma _8^{46}`$. As a reminder, $`n_B`$ is the power index of $`\mathrm{\Delta }_B^2(k,z)`$ at $`k=l/x`$ and $`0<n_B1`$ for $`l10^3`$ and $`z1`$. So there is a degeneracy between $`z_{\mathrm{reion}}`$, $`\mathrm{\Omega }_bh`$ and $`\sigma _8`$ and the recovered $`z_{\mathrm{reion}}(\mathrm{\Omega }_bh)^6\sigma _8^{15}`$. Thus, current uncertainties in cosmological parameters degrade the accuracy of the extracted $`z_{\mathrm{reion}}`$ significantly. But this degenerancy is not that bad as appears since it is the matter density power spectrum $`P_{NL}(k)`$ at Mpc/h scale, instead of $`\sigma _8`$, that directly determines the kinetic SZ effect at $`l`$several thousands through $`C_lP_{NL}^2`$. $`P_{NL}(k)`$ at Mpc/h scale will be measured by CFHT legacy survey with $`1\%`$ accuracy (see Pen et al. (2003) for the discussion of relevant length scales in the recovered 3D power spectrum). Future experiments such as Planck, SNAP and CFHT would significantly reduce the above uncertainties in cosmological parameters and density power spectrum (Tegmark & Zaldarriaga, 2002) and allow a precise recovery of the reionization history from the kinetic SZ effect.

At the presence of various noise sources in KSZ measurement discussed above, ACT would still be able to measure $`C_l`$ with $`2\%`$ accuracy. Given a sufficient understanding of cosmic baryonic matter density and the matter density power spectrum, which will be measured precisely by future experiments, the reionization redshift could be constrained with several percent accuracy, some multiple reionization scenarios degenerated in the primary CMB temperature and T-E polarization measurement could be distinguished by $`10\sigma `$ confidence.

By far we only discussed the uniformly reionized universe. Patchy reionization introduces extra power to the KSZ effect and may degrade the power of KSZ to constrain the reionization history. But future observations may be able to separate it by its characteristic $`C_l`$ behavior. Furthermore, by cross correlating the KSZ effect with the cosmic 21 cm background, one can further recover the redshift information of the KSZ effect at high redshifts and extract the contribution from the patchy reionization. Since in the literature, the significance of patchy reionization is under debate (e.g. Gnedin &Jaffe (2001); Valageas et al. (2001); Santos et al. (2003)), we postpone its further discussion.

## 8 Summary

The dependence of the kinetic SZ effect on the reionization history allows a detailed probe of the reionization process and the feasibility to break the degeneracy posed in the primary CMB measurement. This relies on the precision of the kinetic SZ effect measurement, which will be realized by future CMB experiments such as ACT and precision modeling of the kinetic SZ effect. We have built an analytical model for the kinetic SZ effect and tested against our high resolution hydro simulations through the 3D power spectrum $`\mathrm{\Delta }_B^2(k,z)`$ of the momentum curl part, which completely determines the kinetic SZ effect angular power spectrum. The grid version of our model, which captures various limitations of simulations such as the finite box size effect and thus reproduces simulation results at all available simulations scales and redshifts with better than several percent accuracy. The continuous version, which represents an ideal simulation with infinite box size and resolution and is free of simulation limitations, should describe the real $`\mathrm{\Delta }_B^2`$ to a high accuracy. Comparing between simulations and analytical predictions, we found that

* The simulated $`\mathrm{\Delta }_B^2`$ loses power at large scales due to the cutoff of the density fluctuation beyond the simulation box. This effect causes the simulated $`\mathrm{\Delta }_B^2`$ to be steeper than the real one at large $`k`$. The smaller the simulation box, the severe the problem is.
* The simulated $`\mathrm{\Delta }_B^2`$ loses power at small scales due to the finite resolution.
* Correspondingly, the simulated kinetic SZ power spectrum loses power at both large and small scales. In order for simulation to directly produce reliable $`C_l`$ at $`1000<l<10000`$ without reference to our models, a hydro simulation with box size $`>400`$ Mpc/h and resolution $`>4096^3`$ is required.
* The effect of the non-linearity on $`\mathrm{\Delta }_B^2`$ has the same tendency as its effect on the density power spectrum. Its effect can be modeled by substituting both linear density power spectra in the Vishniac integral (Eq. 4 & 5) with the corresponding non-linear ones (Eq. 6 & 7). For WMAP cosmology, the non-linearity increases the density power spectrum and $`\mathrm{\Delta }_B^2`$.
* Our analytical model in the grid version captures all main features of simulated $`\mathrm{\Delta }_B^2`$ and reproduces the simulation results with high accuracy. Its continuous version is an accurate description of the real $`\mathrm{\Delta }_B^2`$ and therefore the KSZ effect over a wider parameter range in power spectrum slope than the past models proposed in the literature.
* For the WMAP cosmology and the single reionization scenario, the kinetic SZ power spectrum can be approximated with better than several percent accuracy as $`l^2C_l/(2\pi )0.91\times 10^{12}[(1+z_{\mathrm{reion}})/10]^{0.34}(l/5000)^{0.230.015(z_{\mathrm{reion}}9)}`$ for $`6<z_{\mathrm{reion}}<20`$ and $`3000<l<9000`$. This strong reionization redshift dependence will allow a precision measurement of the reionization history. For the WMAP favored $`z_{\mathrm{reion}}16.5`$, the kinetic SZ power spectrum has a broad peak extending from $`l3000`$-$`10^5`$ with an amplitude $`2.7\mu `$k. Among cosmological parameters, $`\sigma _8`$ and $`\mathrm{\Omega }_b`$ are the two that $`C_l`$ is most sensitive to. $`C_l`$ scales as $`(\mathrm{\Omega }_bh)^2\sigma _8^{46}`$.
* ACT will measure the kinetic SZ effect to $`1\%`$ accuracy. Given precise measurements of cosmic baryonic density $`\mathrm{\Omega }_bh`$ and matter density power spectrum, for the single reionization scenario, the reionization redshift $`z_{\mathrm{reion}}`$ can be constrained with several percent accuracy. The kinetic SZ effect can further distinguish some more complicated reionization scenarios and breaks the degeneracy met in the primary CMB measurement with $`10\sigma `$ confidence.

## Acknowledgments

We would like to thank Uros Seljak and Steen Hansen for helpful discussions. We thank Lyman Page for providing information about ACT. We thank Asantha Cooray, Lloyd Knox and Mario Santos for the discussion of patchy reionization. P.J. Zhang thanks the department of Astronomy & Astrophysics and CITA at the University of Toronto where part of the work was done. P.J. Zhang is supported by the DOE and the NASA grant NAG 5-10842 at Fermilab. Computations were performed on the CITA Pscinet computers funded by the Canada Foundation for Innovation.