format: text

text:
# The Bispectrum and the trispectrum of the Ostriker-Vishniac effect

## I Introduction

In recent years, with the prospect of the increase in the sensitivity and angular resolution of the forthcoming Cosmic Microwave Background (CMB) satellite and interferometry experiments MAP/PLANCK ; MINT ; CBI ; ACT , efforts have been driven to the theoretical study of the secondary anisotropies contributions to the temperature fluctuations on arcminute scales and below. While the primordial anisotropies from recombination are thought to be well understood secondary anisotropies from reionization are not.

As is well known, the current favored inflationary model of structure formation predicts a nearly Gaussian probability distribution for the primordial density fluctuations. In this case, the CMB is completely described by its two-point correlation function or power spectrum in Fourier space. All higher-order correlations can be expressed in terms of it. Primordial nonlinearities and secondary effects introduce deviations from Gaussianity, producing a detectable signal in both the power spectrum and higher-order statistics. Recent work provides the theoretical background for the calculation of estimators of these higher-order statistics Hu01 ; Komatsu:Spergel00  and constrains possible non-Gaussian primordial contributions to the bispectrum and the trispectrum on degree and sub-degree angular scales using actual data Kunz01 ; Santos01 ; Komatsu01 ; Troia03 . The interest is now in forecasting the expected signals on smaller scales due to secondary anisotropies, checking whether they are detectable and understanding how they can be separated from each other and from the primary anisotropies in light of the future data.

The Ostriker and Vishniac (OV) effect Ost:Vish86  was found to be the dominant linear secondary contribution to the CMB anisotropies below the Silk-damping scales at the arcminute level Hu94 . It is caused by Thomson scattering off of the CMB photons by moving electrons during the initial phase of reionization. It has the advantage of taking place during the linear regime of structure evolution and of being a small-angle effect enabling one to obtain analytical expressions for its higher-order correlation functions in the small-angle limit. Because of the highly predictive power of linear theory, any measurement of such statistics would be a sensitive probe of the reionization history of the universe, difficult to disentangle in measurements from nonlinear contributions. Several derivations for its power spectrum have been carried out Hu94 ; Dod:Jubas93 ; Vishniac87  but no one has fully addressed the calculation of its bispectrum or trispectrum.

It is then timely to obtain these expressions and to qualify and quantify them. We therefore extend and detail previous techniques used for the calculation of the OV power spectrum to the calculation of its higher-order statistics. As it will be shown, for the particular case of fields which are vector-like in nature, such as the OV effect, even moments will dominate over odd moments, making the trispectrum a much more sensitive statistics than the bispectrum.

Given the low redshifts of formation of structure, it is interesting to consider whether nonlinear effects can further enhance these statistics. So we will extend our study to the weakly nonlinear regime, allowing us to probe the most natural extension of the OV effect to nonlinear scales, the so-called kinetic Sunyaev-Zel’dovich (kSZ) effect. On small scales, both arise from the density modulation of the Doppler effect from large-scale bulk flows.

We review the relevant properties and parameters of the adiabatic cold dark matter (CDM) cosmology for structure formation in Sec. II.1. In Sec. II.2 we review the theory of the OV effect and in Sec. II.3 we discuss the basic statistical properties of a general field through its n-point functions. In Sec. II.4 we show how the homogeneous theory of turbulence combined with the Limber approximation enables one to infer the dominant contribution among n-point statistics of a vector field effect like the OV effect. In Sec. II.5, we introduce the standard formalism of the signal-to-noise calculation by means of the Fisher matrix. For general consistency, in Sec. II.6 we apply our method in detail to the calculation of the OV power spectrum as well as its nonlinear extension. In Sec. III and in Sec. IV we present the steps of the calculation of the OV bispectrum and of the trispectrum respectively and its nonlinear counterparts. We also present the results. Finally in Sec. V we conclude. In the Appendix A we generalize the Limber approximation to the 3-point and 4-point correlation functions. This may be useful for other cosmological studies.

## II General Considerations

### II.1 Cosmological model

We work in the context of the adiabatic cold dark matter (CDM) family of models. In units of the critical density, $`\mathrm{\Omega }_0`$ is the contribution to the non-relativistic-matter density, $`\mathrm{\Omega }_b`$ is the contribution to the baryonic matter density, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is the contribution of the cosmological constant and $`H_0=100h`$ km sec-1 Mpc-1 is the Hubble constant today. The Friedmann equations for the evolution of the scale factor of the Universe, $`a(t)`$, are then

$$\frac{\dot{a}}{a}=H_0E(z)H_0\sqrt{\mathrm{\Omega }_0(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }+(1\mathrm{\Omega }_0\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2},$$
(1)
$$\frac{\ddot{a}}{a}=H_0^2[\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Omega }_0(1+z)^3/2],$$
(2)
where the over-dot denotes a derivative with respect to time. The scale factor is chosen such that $`a_0H_0=2c`$.

Useful measures of distance (and time) are the conformal distance (and conformal time). If an observer is at the origin $`z=0`$ then an object at redshift $`z`$ is at a comoving distance, $`w(z)=\frac{1}{2}_0^z\frac{dz^{}}{E(z^{})}`$ and at a time $`t(z)=\frac{1}{H_0}_z^{\mathrm{}}\frac{dz^{}}{(1+z^{})E(z^{})}`$. The conformal time is then obtained from $`d\eta =dt/a`$, such that the comoving distance to the horizon is the conformal time today, $`c\eta _0=w(\mathrm{})`$.

If the CDM density contrast at comoving position $`\stackrel{}{w}`$ at time $`t`$ is $`\delta (\stackrel{}{w},t)`$, then the power spectrum $`P(k,t)`$ is defined by the expectation value over all realizations $`<\stackrel{~}{\delta }(\stackrel{}{k},t)\stackrel{~}{\delta }^{}(\stackrel{}{k}^{},t)>=(2\pi )^3\delta _D^3(\stackrel{}{k}\stackrel{}{k}^{})P(k,t)`$ where $`\delta _D^3`$ is the Dirac delta function. In linear theory, $`\delta (\stackrel{}{w},t)=\delta _0(\stackrel{}{w})D(t)/D(t_0)`$, where $`t_0`$ is the age of the Universe, $`\delta _0(\stackrel{}{w})\delta (\stackrel{}{w},t_0)`$, and the growth factor, as a function of redshift, is Peebles

$$D(z)=\frac{5\mathrm{\Omega }_0E(z)}{2}_z^{\mathrm{}}\frac{1+z^{}}{[E(z^{})]^3}𝑑z^{}.$$
(3)
The power spectrum is given by $`P(k,t)=P(k)(D/D_0)^2`$, where $`D_0D(t_0)`$. For $`P(k)P(k,t_0)`$ we use

$`P(k)`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{8}}\delta _H^2(k/2)^nT^2(k_p\mathrm{Mpc}/h\mathrm{\Gamma }),`$ (4)
where $`T(q)`$ is the CDM transfer function, $`k_p=ka_0=kH_0/2c`$ is the physical wave number with our conventions and $`\mathrm{\Gamma }`$, the shape parameter, is defined as Hu:Sugiyama96  $`\mathrm{\Gamma }\mathrm{\Omega }_0(h/0.5)\mathrm{exp}(\mathrm{\Omega }_b\mathrm{\Omega }_b/\mathrm{\Omega }_0)`$. As we chose $`a_0H_0=2c`$, there is an extra factor of 8 in the denominator in Eq. (4). For the transfer function, we use the Bardeen *et al.* Bardeen86  fitting formulae for CDM models instead of the improved version of Eisenstein and Hu Eisenstein:Hu99 , to facilitate comparison with previous work. For $`\delta _H`$, we take the fits to the Cosmic Background Explorer given in Bunn:White97 .

In linear theory, the continuity equation relates the Fourier components of the velocity field and the density field

$$\stackrel{~}{\stackrel{}{v}}(\stackrel{}{k},t)=\frac{ia(t)}{k^2}\frac{\dot{D}}{D}\stackrel{}{k}\stackrel{~}{\delta }(\stackrel{}{k},t)=\frac{ia(t)}{k^2}\frac{\dot{D}}{D_0}\stackrel{}{k}\stackrel{~}{\delta }_0(\stackrel{}{k}).$$
(5)
A useful relation Peebles  is

$$\frac{\dot{D}}{D}=\frac{\ddot{a}}{\dot{a}}\frac{\dot{a}}{a}+\frac{5\mathrm{\Omega }_0}{2}\frac{\dot{a}}{a}\frac{(1+z)^2}{[E(z)]^2D(z)}.$$
(6)
Although we maintain generality in the derivations, we illustrate our results with a flat model, the $`\mathrm{\Lambda }CDM`$ model. The parameters for this model are $`\mathrm{\Omega }_0=0.35`$, $`\mathrm{\Omega }_b=0.05`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, $`h=0.65`$ and spectral index $`n=1`$. Concerning the reionization contribution we consider two reionization histories, both assuming steep reionization with ionization fraction $`x_e=1`$ and $`\mathrm{\Delta }z_r/(1+z_r)=0.1`$. In the first one, reionization takes place at $`z_r=8`$. The second one assumes $`z_r=17`$ and relies on the latest results from the Wilkinson Microwave Anisotropy Probe (WMAP) experiment (see below). Note that in an open or closed universe one replaces in the factors of $`\eta `$ that appear in the equations

$$\eta \frac{S(a_0H_0\eta \sqrt{1\mathrm{\Omega }_0\mathrm{\Omega }_\mathrm{\Lambda }})}{a_0H_0\sqrt{1\mathrm{\Omega }_0\mathrm{\Omega }_\mathrm{\Lambda }}}$$
(7)
where $`S(x)=\mathrm{sinh}x`$ in an open universe and $`S(x)=\mathrm{sin}x`$ in a closed universe.

### II.2 The Ostriker-Vishniac effect

The reionization of the Universe is one of the most important physical processes that took place in the early universe (see Barkana:Loeb01 ; Haiman:Knox99 ). The most accepted sources for reionization, which requires a source of ultra-violet photons, are an early generation of massive stars formed in dwarf galaxies or an early generation of quasars/AGNs in galaxies. In the currently favored adiabatic CDM class of models for structure formation, reionization is expected to occur in the range $`8z_r30`$. Measurements of the CMB anisotropies on sub-degree scales Bernardis00 ; Hanany00 ; Netterfield01  have been used to put an upper bound on the reionization redsfhit of $`z_r30`$ Griffiths99 ; Tegmark00 . Very recently, using polarization and temperature anisotropies of the CMB, WMAP has placed a fairly model-independent constrain Hu:Holder03  on the optical depth to electron scattering of $`\tau =0.17\pm 0.04`$ at $`68`$% CL which translates into $`z_r=17\pm 3`$ for instant reionization Kogutetal03 ; Spergeletal03 . Interestingly, the measurement of an increase of the neutral fractions with redshift in high-$`z`$ quasar spectra Songaila:Cowie02 ; Djorgovski01  and a first detection of the ’Gunn-Peterson trough’ Gunn:Peterson65  in a quasar spectra at $`z=6.28`$ by SLOAN Becker01  point to a reionization redshift of $`z_r6`$ Gnedin01 , in disagreement with WMAP results. However, even a fraction of neutral hydrogen as small as $`0.1`$% in the IGM could explain the result due to the large cross-section to Ly-$`\alpha `$ photons. Together with the results from WMAP, high-$`z`$ quasar measurements indicate that the reionization history is more complex than previously thought and attempts are being made to fully understand it Haiman&Holder03 ; Hui&Haiman03 ; Wyithe&Loeb03 .

Though extensive analytical (for a complete derivation see Hu94 ; Dod:Jubas93 ) and numerical studies (see references below) have been done to try to quantitatively understand many of the effects originated by reionization on the CMB, accuracy is difficult to reach and uncertainties still remain. Reionization will leave multiple distinctive imprints on the CMB anisotropies by bringing the CMB photons and the free moving electrons into scattering contact again. Through that window, the low-$`z`$ period of the universe evolution can be probed experimentally in more detail from the appearance of the first sources of ionization to the formation of the observable present large-scale structure. Studies have been done on the calculation of the contributions to the power spectrum of the CMB due to ionization induced effects like the Doppler effect on large angular scales Sunyaev78 ; Kaiser84 , the thermal SZ effect and its kinetic analog Bruscoli99 ; Ma:Fry01 ; Refregier:Teyssier00 ; Scaramella93 ; daSilva99 ; Springel00 ; Sunyaev:Zel80 ; Valageas00 , the Inhomogeneous Reionization Aghanim95 ; Benson00 ; Gnedin00 ; Gruzinov:Hu98 ; Knox98  and the OV effect on smaller scales Vishniac87 ; Hu:White96 ; Jaffe:Kam98 . Enlightening comparative studies between different effects can also be found Gnedin:Jaffe00 .

As ionization effects introduce non-Gaussianities in the anisotropies, further studies were done on the calculation of their possible contributions to the bispectrum. Many authors Goldberg:Spergel99 ; Cooray:Hu00 ; Cooray01a ; Cooray01b ; Komatsu02  investigated contributions to mixed bispectra due to couplings between lensing effects, Integrated Sachs-Wolfe effect, thermal SZ and Doppler effects, such as the OV effect. The trispectrum of ionization secondaries Cooray01b  hasn’t been explored very much. No one has addressed the calculation of the pure bispectrum and trispectrum of the OV effect until now.

In the linear regime, for the power spectrum, the dominant small-angular scale contribution from reionization was found to be the OV effect Ost:Vish86 ; Vishniac87 . It arises from the second order modulation of the Doppler effect by density fluctuations which affect the probability of scattering. Because of its density weighting, it peaks at small angular scales, typically arcminute scales in $`\mathrm{\Lambda }`$CDM models, and should produce $`\mu K`$ anisotropies. Its contributions to the temperature fluctuations along the line of sight can be written in the manner of Jaffe and Kamionkowski (JK) Jaffe:Kam98

$$\frac{\mathrm{\Delta }T}{\overline{T}}(\stackrel{}{\theta })=_0^{\eta _0}𝑑\eta g(\eta )\widehat{\theta }.𝐩(\widehat{\theta }\eta ,\eta )$$
(8)
where $`𝐩(\widehat{\theta }\eta ,\eta )𝐯(\widehat{\theta }\eta ,\eta )\delta (\widehat{\theta }\eta ,\eta )`$ and $`g`$ is the visibility function given by

$$g(\eta )=\frac{a(\eta )\sigma _T\overline{n_e}(\eta )}{c}e^{\tau (\eta )}=\frac{0.138\mathrm{\Omega }_bh}{c}(1+z(\eta ))^2x_e(\eta )e^{\tau (\eta )}$$
(9)
which gives the probability of scattering. The prefactor $`0.138`$ is obtained assuming that all the baryons are in the form of protons (if we use the fact that the mass fraction of helium is $`25`$% then one should multiply it by $`7/8`$). The visibility function is normalized such that $`_0^{\eta _0}g(\eta )𝑑\eta =1e^{\tau _r}`$ where $`\tau _r`$ is the optical depth to the surface of last scattering at recombination. Note that $`g`$ is only dependent on time, and not on position for the OV effect. The optical depth is given by $`\tau (\eta )=_0^\eta c𝑑\eta ^{}\sigma _Tn_e(\eta ^{})`$. Also $`\delta (\widehat{\theta }\eta ,\eta )`$ and $`𝐯(\widehat{\theta }\eta ,\eta )`$ are the density contrast and bulk velocity along the line of sight, $`\overline{n_e}(\eta )`$ is the mean electron density given by $`\overline{n_e}(\eta )=\mathrm{\Omega }_b\rho _cx_e(\eta )(1+z)^3/m_p`$, $`\sigma _T`$ is the Thomson cross-section, $`x_e`$ is the ionization fraction and $`m_p`$ the proton mass. We assume that the visibility function is approximated as a Gaussian in conformal time

$$g(\eta )=\frac{1e^{\tau _r}}{\sqrt{2\pi (\mathrm{\Delta }\eta _r)^2}}e^{\frac{1}{2}\frac{(\eta \eta _r)^2}{(\mathrm{\Delta }\eta _r)^2}}.$$
(10)
Following JK Jaffe:Kam98 , we choose a coordinate system such that $`\widehat{\theta }`$ represents a three-dimensional unit vector along the line of sight and $`\stackrel{}{\theta }`$ refers to a two-dimensional unit vector in the plane perpendicular to the line of sight. So we will have $`\stackrel{}{\theta }=(\theta _1,\theta _2,0)`$ and $`\widehat{\theta }=(\theta _1,\theta _2,\sqrt{1\theta _1^2\theta _2^2})(\theta _1,\theta _2,1)`$ where this approximation arises from the small-scale nature of the effect. Bold letters represent three-dimensional vectors.

The OV is a small-angle effect so we can work under the flat-sky approximation and expand the temperature perturbations in Fourier space

$`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K})`$ $`=`$ $`{\displaystyle _0^{\eta _0}}𝑑\eta g(\eta ){\displaystyle d^2\theta \frac{d^3q}{(2\pi )^3}\widehat{\theta }}.\stackrel{~}{𝐩}(𝐪,\eta )e^{i(\stackrel{}{K}.\stackrel{}{\theta }\eta 𝐪.\widehat{\theta })}`$ (11)
where $`𝐪(q_x,q_y,q_z)`$, $`\stackrel{}{K}(K_x,K_y,0)`$ and

$`\stackrel{~}{𝐩}(𝐪,\eta )`$ $`=`$ $`{\displaystyle \frac{ia(\eta )\text{}D}{2D_0^2}}{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\stackrel{~}{\delta }_0(𝐤)\stackrel{~}{\delta }_0(𝐪𝐤)\left(\frac{𝐪𝐤}{|𝐪𝐤|^2}+\frac{𝐤}{|𝐤|^2}\right)}`$ (12)
is the Fourier transform of $`𝐩(\widehat{\theta }\eta ,\eta )`$ (see JK Jaffe:Kam98 ). $`D`$ and  depend on $`\eta `$. We made use of the continuity equation in Fourier space (5) and of the linear evolution of the density field.

### II.3 Statistical properties of a general field

The statistical properties of a field can be characterized by the n-point correlation functions in real space or by the n-point spectra in Fourier space. If the field is Gaussian in nature, like the primordial density fluctuations field in the current favored inflationary cosmology, the connected part of the n-point functions disappears for $`n>2`$. The non-zero (even-n)-point correlation functions can be expressed with the 2-point correlation function. As a result, a Gaussian distribution is completely described by the two-point correlation function, or power spectrum, and any non-Gaussian field will be detectable by measuring the connected part of its n-point correlation function.

If we consider a general statistically homogeneous and isotropic 2-dimensional field $`\stackrel{}{p}`$ with zero mean, its power spectrum $`P`$, bispectrum $`B`$ and trispectrum $`T`$ are defined by the following equations in the appropriate Fourier convention:

$`<\stackrel{}{p}(\stackrel{}{k}_1)\stackrel{}{p}(\stackrel{}{k}_2)>`$ $`=`$ $`(2\pi )^2P(k_1)\delta _D^2(\stackrel{}{k}_1+\stackrel{}{k}_2)`$
$`<\stackrel{}{p}(\stackrel{}{k}_1)\stackrel{}{p}(\stackrel{}{k}_2)\stackrel{}{p}(\stackrel{}{k}_3)>_c`$ $`=`$ $`(2\pi )^2B(k_1,k_2,k_3)\delta _D^2(\stackrel{}{k}_1+\stackrel{}{k}_2+\stackrel{}{k}_3)`$ (13)
$`<\stackrel{}{p}(\stackrel{}{k}_1)\stackrel{}{p}(\stackrel{}{k}_2)\stackrel{}{p}(\stackrel{}{k}_3)\stackrel{}{p}(\stackrel{}{k}_4)>_c`$ $`=`$ $`(2\pi )^2T(k_1,k_2,k_3,k_4)\delta _D^2(\stackrel{}{k}_1+\stackrel{}{k}_2+\stackrel{}{k}_3+\stackrel{}{k}_4)`$
where the subscript $`c`$ stands for connected. The OV effect being a secondary effect will introduce non-Gaussianities in the original primordial Gaussian distributed temperature fluctuations. As a consequence, contributions to its bispectrum and trispectrum are expected.

### II.4 Dominant contributions among the statistics of the OV effect

Combining the homogeneous theory of turbulence with the Limber approximation enables one to infer the dominant contribution among n-point statistics of an isotropic and homogeneous vector-like field effect whose statistical properties vary slowly in time. In particular, we can apply this to the OV effect. In short, the theory of homogeneous turbulence shows how to build invariant spectral tensors of arbitrary order, corresponding to expectation values of arbitrary products of statistically homogeneous vector fields. It is based on techniques proposed in the area of homogeneous turbulence in the 1940s by Robertson Robertson40 . Relying on this theory, all expectation values of an odd product of an isotropic $`3`$-dimensional vector field $`𝐩(𝐪)`$ with $`𝐪(q_x,q_y,q_z)`$ must be proportional to at least one of the $`𝐪`$ vectors, contrary to the expectation values of even products.

Because of the Limber approximation, extended to higher-statistics in Appendix A, which states that the only contributions to the projected correlation function on the sky come from the Fourier modes perpendicular to the line of sight of the angular correlation function, all the $`q_{iz}`$ terms tend to be suppressed. There will be different levels of suppression depending on the order of the $`q_{iz}`$ dependence of our statistics.

Combining these two results, we can conclude that even correlation functions of the OV effect dominate over odd correlation functions making the trispectrum a much more sensitive statistics than the bispectrum. Also, we expect the correlation functions to obey the homogeneous and isotropy theory fully and thus to be the result of contributions of different orders in the $`q_{iz}`$ terms. We have developed a method which permits to calculate the dominant contribution, under the Limber approximation.

This is a characteristic of all effects physically described by an isotropic vector field and can thus be useful for other studies. As noted previously by Scannapieco Scannapieco00 , the alternation of dominant/subdominant/dominant higher order correlation functions provides a unique signal distinguishing the OV effect from other non-vector like secondary anisotropies and consequently can enable one to disentangle it from other contributions at similar angular scales.

### II.5 Signal-to-noise

A fundamental issue is to know how well we can separate the OV signal, which is non-Gaussian, from the Gaussian signal, noise and foregrounds which are always present in a measurement from an experiment. A way of quantifying this detection is to calculate the $`\chi ^2`$ statistics (as in Goldberg:Spergel99 ; Spergel:Goldberg99 ). To do so, one needs to calculate the Fisher information matrix $`F_{ij}`$ (for a good review see Tegmark96 ). If we think of the data $`𝐱`$ as a random variable with a likelihood function $`L(𝐱;\theta )`$ where $`\theta `$ is a vector of model parameters, the Fisher information matrix is defined as

$$F_{ij}\frac{^2lnL(𝐱;\theta )}{\theta _i\theta _j}.$$
(14)
By a very powerful theorem, called the Cramer-Rao inequality, it was shown Kenney:Keeping51 ; Kendall:Stuart69  that the variance of any unbiased estimator of a certain parameter in a model, can not be less than $`(F^1)_{ii}`$. As the signal calculated is expected to be rather small, we are interested in estimating its overall detectability as in Goldberg:Spergel99 ; Cooray:Hu00 ; Spergel:Goldberg99 . Therefore, we assume that the form of our model $`\theta `$ (in our case the bispectrum and the trispectrum) is correct and that the only interesting parameter is its amplitude $`A`$, where the true value of $`A=1`$. Then the Cramer-Rao inequality tells us that the variance of the measurement of $`A`$ is no less than $`\sigma ^2(A)=(F^1)_{AA}`$ and we define the $`\chi ^2`$ statistics as

$$\chi ^2\left(\frac{S}{N}\right)^2=\frac{1}{\sigma ^2(A)}=(F)_{AA}.$$
(15)
The calculation of the Fisher matrix $`(F)_{AA}`$ of the statistics of the OV effect involves the calculation of the contribution of the noise to the power spectrum $`C_{\mathrm{}}^{noise}`$, as it will be shown. The noise depends on the experiment characteristics.

We consider a hypothetical experiment which maps a fraction of the sky $`f_{sky}`$ with a Gaussian beam with full width at half maximum $`\theta _{fwhm}`$ and pixel noise $`\sigma _p=s/\sqrt{t_{pix}}`$, where $`s`$ is the detector sensitivity and $`t_{pix}`$ is the time spent observing each pixel. We use the inverse weight per solid angle, $`w^1(\sigma _p\theta _{fwhm}/T_0)^2`$, in order to have a measure which is independent of the pixel-size Jungman96 ; Knox95 . $`T_0=2.73K`$ is the CMB thermodynamic temperature. If only a fraction $`f_{sky}`$ of the sky is mapped, treating the pixel noise as Gaussian and ignoring any correlations between pixels a good estimate of the $`C_{\mathrm{}}^{noise}`$ Cooray:Hu00 ; Jungman96 ; Knox95  is

$$C_{\mathrm{}}^{noise}=f_{sky}w^1e^{\sigma _b^2\mathrm{}(\mathrm{}+1)}$$
(16)
where $`\sigma _b`$, in radians, is the width of the beam if we assume it has a Gaussian profile. It is related to $`\theta _{fwhm}`$, in arcminutes, by $`\sigma _b=\sqrt{8ln2}\theta _{fwhm}\times \pi /10800`$. Note that if an experiment maps the full sky and then a fraction $`1f_{sky}`$ is subtracted, one should not multiply $`w^1`$ by $`f_{sky}`$ (case of MAP and Planck). Hence we can estimate $`C_{\mathrm{}}^{noise}`$ for any experiment with characteristic $`f_{sky}`$, $`\theta _{fwhm}`$, $`s`$ and $`t_{pix}`$.

For the precise cases of MAP (renamed WMAP recently) and Planck, for which we used the specifications in table 1, we need to take into account their multi-frequency coverage with different characteristics. The $`C_{\mathrm{}}^{noise}`$ is then defined as Cooray:Hu00

$$\frac{1}{C_{\mathrm{}}^{noise}}=\underset{\nu }{}\frac{1}{C_{\mathrm{}}^{noise}(\nu )}$$
(17)
where the sum runs over all channels of the experiment and $`\nu `$ is the frequency of the channel.

### II.6 Power spectrum of the OV effect

#### II.6.1 Linear power spectrum

In Fourier space, the flat-sky power spectrum of the OV effect is related (Eq. (13)) to the following two-point expectation value of the OV temperature field perturbation $`\mathrm{\Delta }T/\overline{T}`$

$`<{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_1}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_2})>`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^{\eta _0}}𝑑\eta _1g(\eta _1){\displaystyle _0^{\eta _0}}𝑑\eta _2g(\eta _2){\displaystyle d^2\theta _1d^2\theta _2}`$ (18)
$`{\displaystyle \frac{d^3q_1}{(2\pi )^3}\frac{d^3q_2}{(2\pi )^3}\widehat{\theta }_{1i}\widehat{\theta }_{2j}(<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)>+<\stackrel{~}{p}_i(𝐪_2,\eta _2)\stackrel{~}{p}_j(𝐪_1,\eta _1)>)}`$
$`e^{i(\stackrel{}{K_1}.\stackrel{}{\theta }_1\eta _1𝐪_1\widehat{\theta }_1)}e^{i(\stackrel{}{K_2}.\stackrel{}{\theta }_2\eta _2𝐪_2.\widehat{\theta }_2)}`$
where $`\stackrel{~}{𝐩}`$ is defined as in Eq. (12). Many authors have derived the expression for this statistics Hu94 ; Dod:Jubas93 ; Vishniac87 ; Jaffe:Kam98 ; Scannapieco00 . We use the JK formalism but a different technique which will be useful in what follows.

As we see, this expression involves a double integration in time, angle and wavenumber, being numerically long to evaluate. It is useful to note that as the statistical properties of the field $`𝐩`$ vary slowly in time and as the OV effect is a small-angle effect we can employ the Limber approximation (see A1) to considerably simplify our derivations. We stress here that we are allowed to use this approximation as the OV effect takes place at sufficiently high $`l`$, where the difference between the approximation and the integral is very small.

As the two permutations $`<\stackrel{~}{p}_i\stackrel{~}{p}_j>`$ are symmetric due to statistical homogeneity, we only consider the first one and multiply the result by $`2`$. Using Eq. (12) for $`\stackrel{~}{p}`$ and the Wick theorem for the Gaussian 3-dimensional density field correlation function which states

$`<\delta (𝐤_1)\delta (𝐪_1𝐤_1)\delta (𝐤_2)\delta (𝐪_2𝐤_2)>=`$ (19)
$`(2\pi )^6P(k_1)P(𝐪_1𝐤_1)(\delta _D^3(𝐤_1+𝐪_2𝐤_2)\delta _D^3(𝐪_1+𝐤_2𝐤_1)+\delta _D^3(𝐤_1+𝐤_2)\delta _D^3(𝐪_1𝐤_1+𝐪_2𝐤_2))`$
where $`P(k)`$ is the power spectrum of density perturbations, we obtain two non-zero terms for $`<\stackrel{~}{p}_i\stackrel{~}{p}_j>`$ which can be written as

$$<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)>=2G(\eta _1)G(\eta _2)F_{ij}(𝐪_1)\delta _D^3(𝐪_1+𝐪_2)$$
(20)
with the time dependence functions gathered in $`G(\eta )=\left(\frac{iaD\dot{D}}{2D_0^2}\right)`$ and the general tensorial functions $`F_{\alpha \beta }`$ given by

$`F_{\alpha \beta }(𝐪_i)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)\left(\frac{a_\alpha }{a^2}+\frac{b_\alpha }{b^2}\right)\left(\frac{a_\beta }{a^2}+\frac{b_\beta }{b^2}\right)}`$ (21)
where $`𝐚=𝐊^{}`$ and $`𝐛=𝐪_i𝐊^{}`$. We could now replace the two previous expressions directly in Eq. (18) but, for clarity purposes as it will become obvious soon, we refrain from doing so and instead we keep working with $`F_{\alpha \beta }`$.

Indeed, in the small-sky approximation, for which the unit vector $`\widehat{\theta }(0,0,1)`$, we can contract the $`\widehat{\theta }`$’s of expression (18) with the vectors $`𝐚`$ and $`𝐛`$ of the last expression (21) such that we are left with the line of sight components of $`𝐚`$ and $`𝐛`$. We can thus define a new scalar function $`F`$ such that

$`F(𝐪_i)=\widehat{\theta }_{1\alpha }\widehat{\theta }_{2\beta }F_{\alpha \beta }`$ $``$ $`{\displaystyle d^3K^{}P(a)P(b)\left(\frac{a_z}{a^2}+\frac{b_z}{b^2}\right)^2}={\displaystyle d^3K^{}P(a)P(b)\left(\frac{a_z^2}{a^4}+\frac{2a_zb_z}{a^2b^2}+\frac{b_z^2}{b^4}\right)}.`$ (22)
The interesting step that follows is to expand this function in $`a_z=K_z^{}`$ and $`b_z=q_{iz}K_z^{}`$

$`F(𝐪_i)={\displaystyle d^3K^{}P(a)P(b)\left(K_z^2\left(\frac{1}{a^4}+\frac{1}{b^4}\frac{2}{a^2b^2}\right)+K_z^{}\left(\frac{2q_{iz}}{a^2b^2}\frac{2q_{iz}}{b^4}\right)+\left(\frac{q_{iz}^2}{b^4}\right)\right)}.`$ (23)
Expressions of this type will occur in the following and they illustrate the previous discussion in Sec. II.4. By looking at the integral, we see that the $`K_z^{}`$ integration for odd products of $`K_z^{}`$ is zero as the dependence on $`K_z^{}`$ of the terms with $`a`$ and $`b`$ is even. Therefore we are left with non-zero contributions from the first term (with $`K_z^2`$) and from the last term (with $`K_z^0`$). By applying the Limber approximation we can also infer that the dominant contribution has to come from the $`K_z^2`$ term as it has no explicit dependency on $`q_{iz}`$. As we know, in the Limber approximation, Fourier modes parallel to the line of sight (terms on $`q_{iz}`$) tend to be suppressed. We are then left with $`2`$ terms which we will label dominant and sub-dominant terms depending on the order of their cancellation. We are interested in the dominant one. This cancellation was forecasted in Sec. II.4 in view of the homogeneous turbulence theory. Indeed, the power spectrum of the OV effect is expected to have contributions only from terms with no dependency on $`q_{iz}`$ and with a $`q_{iz}q_{jz}`$ type dependency.

After some straightforward algebra, in the Limber approximation framework, we find for the dominant contribution

$`F(𝐪_i)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)\left(K_z^2\left(\frac{1}{a^4}+\frac{1}{b^4}\frac{2}{a^2b^2}\right)\right)}`$ (24)
$`=`$ $`2\pi q_i{\displaystyle _0^{\mathrm{}}}𝑑y_1{\displaystyle _1^1}𝑑\mu P(q_iy_1)P(q_iy_2){\displaystyle \frac{(1\mu ^2)(12\mu y_1)^2}{y_2^4}}`$
where we have performed a spherical coordinate transform such that $`\mu =\widehat{q}_i.\widehat{K^{}}`$, $`a=y_1q_i`$ and $`b=y_2q_i`$ with $`y_2=\sqrt{1+y_1^22y_1\mu }`$. To obtain the components of $`K_z^{}`$ in the chosen coordinate system, we used the Limber approximation to assume $`q_{iz}0`$, such that $`K_z^{}=K^{}\sqrt{1\mu ^2}`$. This assumption preserves our dominant term but suppresses any sub-dominant term that could naturally arise when calculating the integral. As a consequence, in Eq. (23), the dominant term (in $`q_{iz}^0`$) may contain hidden contributions to the subdominant term (in $`q_{iz}^2`$). That this indeed is the case can be understood by a very simple reasoning. Consider Eq. (22). It is easy to show that the terms in $`a_z^2`$ and $`b_z^2`$ give identical contributions to the integral, such that if we calculate twice the integral in $`a_z^2`$ we should obtain the same result at the end, i.e., various terms depending on different orders in $`q_{iz}`$. By doing this, our term in $`q_{iz}^2`$ present in Eq. (23) simply disappears. We might then expect it to show up in the integral calculation. But, most interestingly, when performing the calculations as previously using the Limber approximation, we obtain the same result, i.e., Eq. (24). As a consequence, by imposing $`q_{iz}0`$ in this example, we are in fact suppressing a sub-dominant term we know should be present. In conclusion, there are more terms contributing to the sub-dominant power spectrum than the one present in Eq. (23) and these aren’t so easily calculated. Hence we neglect all the sub-dominant terms in the Limber approximation. In the following, for the bispectrum and the trispectrum calculations, similar problems are present but they will not affect the lowest order terms in $`q_{iz}`$ and even order in $`K_z^{}`$, which are the dominant terms of interest to us.

We can finally combine Eqs. (18), (20) and (24). Applying Kaiser’s method and the Limber approximation as described in the Appendix, we obtain the well-known expression for the dominant contribution to the linear OV power spectrum

$$P_{dom}^{OV}(K)=\frac{1}{8\pi ^2}_0^{\eta _0}\frac{g^2(\eta )}{\eta ^2}[a(\eta )]^2\left(\frac{\text{}D}{D_0}\right)^2S(K/\eta )𝑑\eta $$
(25)
where

$$S(k)=k_0^{\mathrm{}}𝑑y_1_1^1𝑑\mu P(ky_1)P(k\sqrt{1+y_1^22y_1\mu })\frac{(1\mu ^2)(12\mu y_1)^2}{(1+y_1^22\mu y_1)^2}.$$
(26)
The Limber approximation reduced the dimension of the integral from $`6`$ to $`3`$ and helped to find an easier numerical and analytical expression. We note here that we obtain a difference of a factor $`1/2`$ compared to JK, a discrepancy pointed out by them when comparing to previous work.

For illustration we show plots of the linear dominant contribution of the OV power spectrum in figure (1) for the fiducial $`\mathrm{\Lambda }`$CDM model assuming $`z_r=8`$ and $`z_r=17`$. The correspondence between full-sky multipole moments $`C_{\mathrm{}}`$ and the flat-sky Fourier space $`P(K)`$ is straightforward $`C_{\mathrm{}}=P(K=\mathrm{})`$. As expected, the $`z_r=17`$ scenario (keeping $`x_e=1`$) increases the amplitude of the power spectrum, due to the rise of the optical depth, and shifts its peak towards smaller angular scales. Numerically, we find for the amplitude of the power spectrum the approximate scaling dependence $`C_{\mathrm{}\mathrm{\hspace{0.17em}500}}7.5\times 10^{18}x_e^2\mathrm{log}^{0.4}(1+z_r)`$. We advise the reader to consult JK paper Jaffe:Kam98  for the impact of changing the cosmological parameters (or the reionization history) on the power and peak of the effect as well as for the important detectability issues. The cosmological dependence applies as well for the higher-order statistics.

#### II.6.2 Nonlinear extension: the kSZ effect power spectrum

This extension was calculated previously Hu99  and we present it for the dominant power spectrum for the sake of consistency. The kinetic Sunyaev-Zel’dovich (kSZ) effect results from a Doppler effect suffered by the CMB photons as they travel through large scale structures emerged in a bulk flow. As pointed out by Hu Hu99  among others Ma:Fry01 ; Gnedin:Jaffe00 , in adiabatic CDM cosmologies, nonlinearities only affect the density field below the coherence scale of the bulk velocities and so the nonlinear density field is uncorrelated with the large-scale velocity field, which remains linear. Hence, the mild nonlinear extension of the density contribution in the OV effect expressions naturally becomes the kSZ effect from large scale structures.

As the result of this, we can use the previous calculations of the OV effect power spectrum to obtain a similar expression for the kSZ by introducing a nonlinear correction in the density field contribution. Following Hu Hu99  approach, we replace the linear density power spectrum with its nonlinear analogue but leave the contribution from the velocity power spectrum the same

$$P_{dom}^{NL}(K)=\frac{1}{8\pi ^2}_0^{\eta _0}\frac{g^2(\eta )}{\eta ^2}[a(\eta )]^2\left(\frac{\text{}D}{D_0}\right)^2\frac{P^{NL}(K/\eta )}{P(K/\eta )}S(K/\eta )𝑑\eta $$
(27)
where $`P^{NL}`$ stands for the nonlinear CDM power spectrum and where the mode coupling integral $`S`$ was calculated previously (Eq. (26)) under linear theory. This expression includes both the OV and the kSZ effects.

To calculate the nonlinear power spectrum, we assume that the baryonic gas traces the dark matter Hu99 . We follow Hamilton *et. al.* Hamilton91  who presented a scaling relation for the correlation function in the nonlinear regime that was generalized to its Fourier analog by Peacock and Dodds Pea:Dod94 . The basic hypothesis is that nonlinear fluctuations on a scale $`k`$ arise from linear fluctuations on a larger scale

$$k_{\mathrm{lin}}=[1+\mathrm{\Delta }_{\delta _c}^2(k)]^{1/3}k$$
(28)
so that there is a function relating the nonlinear and linear power spectra at these two scales

$$\mathrm{\Delta }_{\delta _c}^2(k)=f_{\mathrm{NL}}[\mathrm{\Delta }_{\delta _c}^{2(\mathrm{lin})}(k_{\mathrm{lin}})]$$
(29)
which can be fit to simulations. We use the Pea:Dod96  proposed expression for $`f_{\mathrm{NL}}`$. For stronger nonlinearities other corrections are necessary. See for example Zaldarriaga99 . One also needs the relation $`\mathrm{\Delta }_\delta ^2(k)=k^3/(2\pi ^2)P(k)`$.

As discussed in W.Hu Hu99 , this estimative should be seen as an upper limit of the kSZ because on very small scales the gas pressure, unaccounted for, smooths the gas density as compared to the the dark matter density. The assumption that the baryonic gas traces the dark matter was shown to break down at multipoles $`\mathrm{}10^4`$ Hu99 ; Gnedin:Hui98 . We show a plot of the power spectrum corrected for mild nonlinearities in figure (1) for illustration of the nonlinear enhancement of the power spectrum at small scales for $`z_r=8`$ and $`z_r=17`$.

## III The bispectrum of the OV effect

#### III.0.1 Linear bispectrum

In analogy with the power spectrum, the flat-sky bispectrum of the OV effect is connected by Eq. (13) to the following expectation value

$`<{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_1}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_2}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_3})>`$ $`=`$ $`{\displaystyle \frac{1}{6}}{\displaystyle _0^{\eta _0}}𝑑\eta _1g(\eta _1){\displaystyle _0^{\eta _0}}𝑑\eta _2g(\eta _2){\displaystyle _0^{\eta _0}}𝑑\eta _3g(\eta _3){\displaystyle d^2\theta _1d^2\theta _2d^2\theta _3}`$ (30)
$`{\displaystyle }{\displaystyle \frac{d^3q_1}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3q_2}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3q_3}{(2\pi )^3}}\widehat{\theta }_{1i}\widehat{\theta }_{2j}\widehat{\theta }_{3l}(<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)\stackrel{~}{p}_l(𝐪_3,\eta _3)>+perm.)`$
$`e^{i(\stackrel{}{K_1}.\stackrel{}{\theta }_1\eta _1𝐪_1.\widehat{\theta }_1)}e^{i(\stackrel{}{K_2}.\stackrel{}{\theta }_2\eta _2𝐪_2.\widehat{\theta }_2)}e^{i(\stackrel{}{K_3}.\stackrel{}{\theta }_3\eta _3𝐪_3.\widehat{\theta }_3)}`$
where the five permutations are with respect to the ordering $`(𝐪_1,𝐪_2,𝐪_3)`$. Proceeding in a similar way we did for the derivation of the power spectrum, this numerically heavy expression to integrate can be considerably simplified by using a generalization of the Limber equation to higher-order statistics (see Appendix).

We start by calculating a simplified expression for the first permutation $`<\stackrel{~}{p}_i\stackrel{~}{p}_j\stackrel{~}{p}_l>`$. At the end we generalize the results to include the total six permutations. Using expression (12) for $`\stackrel{~}{p}`$ and the Wick theorem for the Gaussian 3-dimensional density field $`6`$-point correlation function $`<\delta (𝐤_1)\delta (𝐪_1𝐤_1)\delta (𝐤_2)\delta (𝐪_2𝐤_2)\delta (𝐤_3)\delta (𝐪_3𝐤_3)>`$ we obtain $`C_2^6.C_2^4/3!=15`$ terms for $`<\stackrel{~}{p}_i\stackrel{~}{p}_j\stackrel{~}{p}_l>`$ of which $`8`$ are non-zero. After some simple calculations, these $`8`$ terms can be condensed into

$$<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)\stackrel{~}{p}_l(𝐪_3,\eta _3)>=4G(\eta _1)G(\eta _2)G(\eta _3)[F_{ijl}(𝐪_1,𝐪_2)+F_{ijl}(𝐪_1,𝐪_3)]\delta _D^3(𝐪_1+𝐪_2+𝐪_3)$$
(31)
with the time dependence function $`G(\eta )=\left(\frac{iaD\dot{D}}{2D_0^2}\right)`$ and the general tensorial functions $`F_{\alpha \beta \gamma }`$ given by

$`F_{\alpha \beta \gamma }(𝐪_i,𝐪_j)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)\left(\frac{a_\alpha }{a^2}\frac{b_\alpha }{b^2}\right)\left(\frac{c_\beta }{c^2}\frac{b_\beta }{b^2}\right)\left(\frac{a_\gamma }{a^2}\frac{c_\gamma }{c^2}\right)}`$ (32)
where $`𝐚=𝐊^{}`$, $`𝐛=𝐊^{}𝐪_i`$ and $`𝐜=𝐊^{}+𝐪_j`$. We concentrate on $`F`$ and combine all the terms in Eq. (31) when appropriate.

Again using the small-angle approximation for which $`\widehat{\theta }(0,0,1)`$, we can contract the $`\widehat{\theta }`$’s with the vectors $`𝐚`$, $`𝐛`$ and $`𝐜`$ such that we are left with the line of sight components of $`𝐚`$, $`𝐛`$ and $`𝐜`$. We can thus define a new scalar function $`F`$ such that

$`F(𝐪_i,𝐪_j)=\widehat{\theta }_{1\alpha }\widehat{\theta }_{2\beta }\widehat{\theta }_{3\gamma }F_{\alpha \beta \gamma }`$ $``$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)\left(\frac{a_z}{a^2}\frac{b_z}{b^2}\right)\left(\frac{c_z}{c^2}\frac{b_z}{b^2}\right)\left(\frac{a_z}{a^2}\frac{c_z}{c^2}\right)}.`$ (33)
Expanding this expression in different orders in $`K_z^{}`$ will give us the different levels of contributions under the Limber approximation to the bispectrum

$`F(𝐪_i,𝐪_j)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)}`$ (34)
$`[K_z^3({\displaystyle \frac{1}{a^2b^4}}{\displaystyle \frac{1}{a^4b^2}}{\displaystyle \frac{1}{a^2c^4}}+{\displaystyle \frac{1}{b^2c^4}}+{\displaystyle \frac{1}{a^4c^2}}{\displaystyle \frac{1}{b^4c^2}})+`$
$`K_z^2\left(q_{iz}\left({\displaystyle \frac{2}{a^2b^4}}+{\displaystyle \frac{1}{a^4b^2}}{\displaystyle \frac{1}{b^2c^4}}+{\displaystyle \frac{2}{b^4c^2}}\right)+q_{jz}\left({\displaystyle \frac{2}{a^2c^4}}+{\displaystyle \frac{2}{b^2c^4}}+{\displaystyle \frac{1}{a^4c^2}}{\displaystyle \frac{1}{b^4c^2}}\right)\right)+`$
$`K_z^{}\left(q_{iz}^2\left({\displaystyle \frac{1}{a^2b^4}}{\displaystyle \frac{1}{b^4c^2}}\right)+q_{jz}^2\left({\displaystyle \frac{1}{b^2c^4}}{\displaystyle \frac{1}{a^2c^4}}\right)+q_{iz}q_{jz}\left({\displaystyle \frac{2}{b^4c^2}}{\displaystyle \frac{2}{c^4b^2}}\right)\right)+`$
$`(q_{iz}^2q_{jz}({\displaystyle \frac{1}{b^4c^2}})+q_{jz}^2q_{iz}({\displaystyle \frac{1}{b^2c^4}}))].`$
Again, and as detailed before in Sec. II.6, the integral in $`K_z^{}`$ for odd products of $`K_z^{}`$ is zero, and we are left with non-zero contributions from the second term (with $`K_z^2`$) and the last term (with $`K_z^0`$). The dominant contribution will come from the $`K_z^2`$ though it has a dependency on $`q_z`$ which, by the Limber approximation, tends to be suppressed compared to terms which don’t have such dependency. This result again confirms the discussion in Sec. II.4, which predicts that the bispectrum of the OV effect has contributions only from terms with a dependency on $`q_{iz}`$ or on products of $`3`$ $`q_z`$.

We are interested on the dominant $`K_z^2`$ term which can be written as

$`F(𝐪_i,𝐪_j)`$ $`=`$ $`q_{iz}f_1(𝐪_i,𝐪_j)+q_{jz}f_2(𝐪_i,𝐪_j)`$ (35)
where

$`f_1(𝐪_i,𝐪_j)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)K_z^2\left(\frac{2}{a^2b^4}+\frac{1}{a^4b^2}\frac{1}{b^2c^4}+\frac{2}{b^4c^2}\right)}`$ (36)
and

$`f_2(𝐪_i,𝐪_j)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)K_z^2\left(\frac{2}{a^2c^4}+\frac{2}{b^2c^4}+\frac{1}{a^4c^2}\frac{1}{b^4c^2}\right)}.`$ (37)
If we were to replace directly these last three expressions in (30) and then apply the Limber approximation we would have a direct cancellation. We wish to obtain an expression non-null which enhances the dominant character, as compared to the other sub-dominant terms, of this second-order contribution. To do so, we use again the small-angle approximation and consider $`\widehat{\theta }_i(0,0,1)`$ to write

$`F(𝐪_i,𝐪_j)`$ $`=`$ $`\widehat{\theta }_i.q_if_1(𝐪_i,𝐪_j)+\widehat{\theta }_j.q_jf_2(𝐪_i,𝐪_j).`$ (38)
Now we replace this particular equation in (31) and then in (30) and finally integrate once by parts the time dependence. This is the step that allows us to keep this term, which was expected to be suppressed under the Limber approximation. Then we follow Kaiser’s method and apply the Limber approximation, and using (13) we obtain the expression for one of the six possible permutation terms of the bispectrum

$$B_{dom}^{OV}(K_1,K_2,K_3)=\frac{i}{16\pi ^3}_0^{\eta _0}𝑑\eta \frac{h^2}{\eta ^4}\frac{h}{\eta }\left[f(\stackrel{}{q}_1,\stackrel{}{q}_2)+f(\stackrel{}{q}_1,\stackrel{}{q}_3)\right]$$
(39)
where $`h(\eta )=\left(g\frac{aD\dot{D}}{D_0^2}\right)`$ and $`f(\stackrel{}{q}_i,\stackrel{}{q}_j)=f_1(\stackrel{}{q}_i,\stackrel{}{q}_j)+f_2(\stackrel{}{q}_i,\stackrel{}{q}_j)`$.

In order to simplify the functions $`f_1`$ and $`f_2`$, we need to find an explicit relation between $`𝐚`$, $`𝐛`$, $`𝐜`$ and $`𝐪_i`$, $`𝐪_j`$ and $`𝐊^{}`$. To do so, we express both $`𝐊^{}`$ and $`𝐚`$, $`𝐛`$, $`𝐜`$ in the basis $`(\widehat{e}_z,\widehat{q}_i,\widehat{q}_j)`$ where we remind the reader that $`\widehat{q}_\alpha (\widehat{q}_\alpha ,0)`$ under the Limber approximation. Again, as for the power spectrum (see discussion in Sec. II.6), this suppresses any sub-dominant contribution that could arise in the process but leaves our dominant contribution intact. We obtain

$`f(𝐪_i,𝐪_j)`$ $`=`$ $`f(𝐪_i,𝐪_j,\widehat{q}_i.\widehat{q}_j)`$ (40)
$`=`$ $`{\displaystyle 𝑑K_z^{}𝑑u𝑑v\sqrt{1\mu ^2}P(a)P(b)P(c)K_z^2\left(\frac{1}{a^4b^2}+\frac{1}{a^4c^2}\frac{2}{a^2b^4}\frac{2}{a^2c^4}+\frac{1}{b^2c^4}+\frac{1}{b^4c^2}\right)}`$
where $`a=(K_z^2+u^2+v^2+2uv\mu )^{1/2}`$, $`b=(K_z^2+(uq_i)^2+v^2+2(uq_i)v\mu )^{1/2}`$, $`c=(K_z^2+(v+q_j)^2+u^2+2(v+q_j)u\mu )^{1/2}`$ and $`\mu =\widehat{q}_i.\widehat{q}_j`$. As this expression corresponds to only one of the permutation terms of expression (30), we need to include the other $`5`$ to obtain the final expression for the dominant flat-sky bispectrum

$$B_{dom}^{OV}(K_1,K_2,K_3)=\frac{i}{8\pi ^3}_0^{\eta _0}𝑑\eta \frac{h^2}{\eta ^4}\frac{h}{\eta }\underset{i,j=1;ij}{\overset{n=3}{}}f(\frac{\stackrel{}{K}_i}{\eta },\frac{\stackrel{}{K}_j}{\eta })$$
(41)
where $`f`$ is given by (40). We thus reduced our initial expression with twelve integrations to a four-dimensional integral, which can be numerically calculated for a chosen configuration of the wavenumbers $`K`$.

Note that we obtain the first order time derivative of $`h`$. It involves one single derivation, which indicates the $`q_{iz}`$ cancellation. As $`h`$ is smooth and slowly varying, the contribution from this term should be very small. Hence, the Limber cancellation at small scales reveals itself in the derivatives of the time dependent functions.

#### III.0.2 Nonlinear extension: the kSZ effect bispectrum

We follow the same approach as for the power spectrum in Sec. II.6 and apply it to the dominant contribution. As we have three CDM power spectrum showing up in the expression for the bispectrum, which we need to divide among linear/nonlinear contributions from density/velocity contributions we will assume the following bispectrum effect

$$B_{dom}^{0VNL}(K_1,K_2,K_3)=\frac{i}{8\pi ^3}_0^{\eta _0}𝑑\eta \frac{h^2}{\eta ^4}\frac{h}{\eta }\underset{i,j,k=1;ij}{\overset{n=3}{}}\left(\frac{P^{NL}(K_k/\eta )}{P(K_k/\eta )}\right)^{3/2}f(\frac{\stackrel{}{K}_i}{\eta },\frac{\stackrel{}{K}_j}{\eta }).$$
(42)
#### III.0.3 The Signal-to-noise

To define the $`\chi ^2`$ statistics of Eq. (15), we need to calculate the likelihood $`L`$ of observing the bispectrum elements $`B_\beta B_{\mathrm{}_1\mathrm{}_2\mathrm{}_3}`$ given the true parameters $`𝐩`$, and calculate the Fisher matrix as defined in Eq. (14)

$$F_{ij}\frac{^2lnL(𝐁;𝐩)}{p_ip_j}.$$
(43)
Assuming that the likelihood is Gaussian, we follow Cooray and Hu Cooray:Hu00  approach to calculate the $`\chi ^2`$ statistics

$$\chi ^2\left(\frac{S}{N}\right)^2=f_{sky}\underset{\mathrm{}_3\mathrm{}_2\mathrm{}_1}{}\frac{B_{\mathrm{}_1\mathrm{}_2\mathrm{}_3}^2}{\sigma _{\mathrm{}_1\mathrm{}_2\mathrm{}_3}^2},$$
(44)
where $`B_{\mathrm{}_1\mathrm{}_2\mathrm{}_3}`$ is the angular averaged bispectrum defined on the sphere, $`f_{sky}`$ represents the reduction in signal-to-noise due to incomplete sky coverage and $`\sigma _{\mathrm{}_1\mathrm{}_2\mathrm{}_3}^2`$ comes from the covariance matrix of the angular averaged bispectrum assuming a nearly Gaussian bispectrum and full-sky coverage Goldberg:Spergel99 ; Gangui:Martin00

$$\sigma _{\mathrm{}_1\mathrm{}_2\mathrm{}_3}^2=C_\mathrm{}_1^{tot}C_\mathrm{}_2^{tot}C_\mathrm{}_3^{tot}[1+\delta _D^3(\mathrm{}_1+\mathrm{}_2)+\delta _D^3(\mathrm{}_2+\mathrm{}_3)+\delta _D^3(\mathrm{}_3+\mathrm{}_1)+2\delta _D^3(\mathrm{}_1+\mathrm{}_2)\delta _D^3(\mathrm{}_1+\mathrm{}_3)]$$
(45)
where $`C_{\mathrm{}}^{tot}`$ stands for the sum of the power spectra of the primary cosmic signal, the thermal SZ (thSZ) effect which contributes significantly at the scales of interest, the linear OV effect, the detector noise and the foregrounds respectively

$$C_{\mathrm{}}^{tot}=C_{\mathrm{}}+C_{\mathrm{}}^{thSZ}+C_{\mathrm{}}^{OVL}+C_{\mathrm{}}^{noise}+C_{\mathrm{}}^{foregrounds}.$$
(46)
The thermal SZ effect was taken from Aghanim02  and was calculated semi-analytically at $`30`$ GHZ for a normalization factor $`\sigma _8=0.9`$. We include the linear OV effect contribution (see Sec. II.6) although its amplitude is small as compared with the primary and the thermal SZ signals at the scales considered in this work. The primary cosmic signal was computed with CMBFAST Seljak:Zaldarriaga96  and we will not consider any foreground. Indeed, for the case of MAP and Planck, studies indicate that the total $`C_{\mathrm{}}`$ should increase by $`10\%`$ maximum Tegmark99 . However, caution is required as this result assumes that foregrounds have a Gaussian distribution. The foregrounds contribution to the higher-order statistics could in fact be our main obstacle in measuring non-Gaussian effects and very little is known about them. For the noise power spectrum we use Eqs. (1617) applied to MAP and Planck.

We will study the contributions to the $`\chi ^2`$ per log interval in $`\mathrm{}`$. It gives us more sensibility on the angular scales of stronger detectability, depending of the effects considered. This also enables us to directly compare our results with Cooray and Hu Cooray:Hu00 . We calculate the total $`S/N`$ for each experiment by integrating over the multipole $`\mathrm{}`$.

A useful relation between the flat-sky bispectrum and the spherical harmonic angular averaged bispectrum Komatsu:Spergel00  is

$$B_{\mathrm{}_1\mathrm{}_2\mathrm{}_3}=\sqrt{\frac{(2\mathrm{}_1+1)(2\mathrm{}_2+1)(2\mathrm{}_3+1)}{4\pi }}\left(\begin{array}{ccc}\mathrm{}_1& \mathrm{}_2& \mathrm{}_3\\ 0& 0& 0\end{array}\right)B(K_1=\mathrm{}_1,K_2=\mathrm{}_2,K_3=\mathrm{}_3).$$
(47)
Since the Wigner-3j vanishes if $`\mathrm{}_1+\mathrm{}_2+\mathrm{}_3=odd`$, the full-sky bispectrum can only be estimated for even terms.

#### III.0.4 Results

We are concerned with the overall detectability of the dominant term, previously calculated. Hence we choose the simplest of the possible configurations, highly localized in Fourier space, $`K_1=K_2=K_3=l`$, for which the flat-sky dominant bispectrum becomes

$$B_{dom}^{OV}(K=\mathrm{})=\frac{3i}{4\pi ^3}_0^{\eta _0}𝑑\eta \frac{h^2}{\eta ^4}\frac{h}{\eta }f\left(\frac{\mathrm{}}{\eta }\right)$$
(48)
where $`f`$ is defined by Eq. (40). The nonlinear analog follows from the previous equation.

With this configuration, depending only on the multipole $`\mathrm{}`$, there is a simple way of calculating an estimate of the order of magnitude of the $`(S/N)^2`$ per bin of $`\mathrm{}`$ were we to include all the remaining $`\mathrm{}`$ modes Komatsu:Spergel00  of the full-sky bispectrum. Indeed in the Eq. (44), we see that the number of modes contributing to the $`(S/N)^2`$ per bin of $`\mathrm{}`$ increases as $`\mathrm{}^2`$ and in Eq. (47) $`\mathrm{}^3\left(\begin{array}{ccc}\mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0\end{array}\right)^2`$ increases as $`0.36\mathrm{}`$ so

$$\frac{d\chi ^2}{d\mathrm{}}f_{sky}\mathrm{}^2\frac{B_{\mathrm{}}^2}{\sigma _{\mathrm{}}^2}f_{sky}\mathrm{}^2\mathrm{}^3\left(\begin{array}{ccc}\mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0\end{array}\right)^2\frac{B^2(\mathrm{})}{\sigma _{\mathrm{}}^2}0.36f_{sky}\mathrm{}^3\frac{B^2(K=\mathrm{})}{\sigma _{\mathrm{}}^2}.$$
(49)
The $`\sigma _{\mathrm{}}^2`$ is calculated using Eq. (45) for all the $`\mathrm{}`$s equal and $`B(\mathrm{})`$ using Eq. (48).

We show the plots for the linear and the nonlinear flat-sky OV bispectrum in figure (2) for $`z_r=8`$ and $`z_r=17`$ in the context of our fiducial cosmological model. For $`z_r=8`$ the peak of the effect occurs around multipole $`l500`$ with an amplitude of $`l^3B(\mathrm{})/(2\pi )1.2\times 10^{23}`$ whereas for $`z_r=17`$ the peak takes place at a higher multipole of $`\mathrm{}700`$ with a stronger amplitude of $`\mathrm{}^3B(\mathrm{})/(2\pi )3\times 10^{23}`$, as expected. We find numerically for the amplitude of the bispectrum the approximate scaling dependence with reionization history $`B(\mathrm{}\mathrm{\hspace{0.17em}50})1.7\times 10^{29}x_e^3\mathrm{log}(1+z_r)`$. The dependence with $`z_r`$ is neither in agreement with the one obtained for the power spectrum (see Sec. II.6) nor with the one obtained for the trispectrum (see Sec. IV). This can be explained by looking at the bispectrum Eq. (48) where we have the first order time derivative of $`h`$, contrary to the corresponding expressions for the power spectrum (Eq. (25)) and the trispectrum (Eq. (70)). This derivative of $`h`$ introduces a stronger scaling relation with $`z_r`$ for the amplitude of the bispectrum as compared to the one for the other two statistics. For both reionization scenarios, the most interesting feature is the rapid drop of power after the peak which is not observed for the power spectrum (see figure (1)), which is a direct consequence of the Limber cancellation at small angular scales. The fact that the bispectrum is not considerable enhanced by nonlinearities, which take place at small angular scales, is also the result of the effect peaking at intermediate multipoles. The second, but expected, result is the low overall amplitude of the effect. The $`z_r=8`$ case (which we can easily compare with previous results in the literature) is lower by more than $`10`$ orders of magnitude than the lensing couplings presented by Cooray and Hu Cooray:Hu00 , and by a few orders of magnitude than most of the OV couplings involving the SZ effect though it should be comparable to the OV couplings involving the Doppler or ISW/SW effects. Despite the fact that such OV couplings suffer the Limber cancellation as well (they correspond to expectation values of an odd product of a vector field), this cancellation can be counterbalanced by the coupling to higher amplitude effects, like the SZ effect, and by the matching in redshift between the density and velocity fields of the OV effect and the secondaries to which is couples. This is the case of the hybrid coupling ISW-SZ-OV presented by Cooray:Hu00 , that has the largest signal of all secondaries that couple with OV.

A more revealing quantity than the power of the effect is the signal-to-noise ratio. We plot in figure (2), the estimated contributions to $`\chi ^2`$ per log interval in $`\mathrm{}`$ of the linear full-sky bispectrum had we included all the modes for MAP, Planck MAP/PLANCK  and no instrumental noise. We show the results for $`z_r=8`$ and $`z_r=17`$. The experimental specifications can be found in table 1. We found no necessity of plotting the corresponding $`d\chi ^2/d\mathrm{}`$ for the nonlinear bispectra due to its similarity to the linear bispectra. The structure in the $`d\chi ^2/d\mathrm{}`$ arises mainly from the structure of the CMB primary power spectrum at $`\mathrm{}>200`$. We point out that the continuous rise at $`\mathrm{}3000`$ for a perfect experiment is due to the fact that the signal decreases slower than the contributions to the $`C_{\mathrm{}}`$ from primary, linear OV and thermal SZ anisotropies up to these scales. Though it is not explicit in the figure, the $`d\chi ^2/d\mathrm{}`$ is zero when $`3\mathrm{}`$ is odd. As we can see, considering thermal SZ contributions to the noise evaluated at $`30`$ GHz, the signal-to-noise of the OV bispectrum is very small, even for a perfect experiment, for which the total $`S/N4.4\times 10^3`$ for $`z_r=8`$ and $`S/N1.4\times 10^2`$ for $`z_r=17`$. These values can be compared to the value $`S/N1.7`$ obtained by Cooray and Hu Cooray:Hu00  for the bispectrum generated by the coupling ISW-SZNL-OV. Hence, even for early reionization, no detection is to be expected from future experiments unable to remove the thermal SZ effect. For the ideal case of a perfect multi-frequency experiment observing in the millimeter and sub-millimeter capable of subtracting all of the thermal SZ effect, the total signal-to-noise increases to $`S/N0.1`$ for $`z_r=17`$, which means that not even for this case are we to expect a detection. Our results for the $`S/N`$ rely on various convenient assumptions required to considerably simplify the analytical expression for the estimate of the $`S/N`$ (see sections Sec. II.5 and Sec. III.0.3). This means that the values obtained for the $`S/N`$ are most probably upper limits on the real values expected and should thus be considered accordingly.

## IV The trispectrum of the OV effect

#### IV.0.1 Linear trispectrum

We follow the same procedure used for the power spectrum and the bispectrum. The trispectrum contribution to the temperature fluctuations is connected by Eq. (13) to

$`<{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_1}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_2}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_3}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }T}}{\overline{T}}}(\stackrel{}{K_4})>=`$ (50)
$`{\displaystyle \frac{1}{24}}{\displaystyle _0^{\eta _0}}𝑑\eta _1g(\eta _1){\displaystyle _0^{\eta _0}}𝑑\eta _2g(\eta _2){\displaystyle _0^{\eta _0}}𝑑\eta _3g(\eta _3){\displaystyle _0^{\eta _0}}𝑑\eta _4g(\eta _4){\displaystyle d^2\theta _1d^2\theta _2d^2\theta _3d^2\theta _4}`$
$`{\displaystyle }{\displaystyle \frac{d^3q_1}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3q_2}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3q_3}{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d^3q_4}{(2\pi )^3}}\widehat{\theta }_{1i}\widehat{\theta }_{2j}\widehat{\theta }_{3l}\widehat{\theta }_{4m}(<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)\stackrel{~}{p}_l(𝐪_3,\eta _3)\stackrel{~}{p}_m(𝐪_4,\eta _4)>+perm.)`$
$`e^{i(\stackrel{}{K_1}.\stackrel{}{\theta }_1\eta _1𝐪_1.\widehat{\theta }_1)}e^{i(\stackrel{}{K_2}.\stackrel{}{\theta }_2\eta _2𝐪_2.\widehat{\theta }_2)}e^{i(\stackrel{}{K_3}.\stackrel{}{\theta }_3\eta _3𝐪_3.\widehat{\theta }_3)}e^{i(\stackrel{}{K_4}.\stackrel{}{\theta }_4\eta _4𝐪_4.\widehat{\theta }_4)}.`$
The 24 total permutations arise from symmetries under permutation invariance. Again we concentrate on the first of the permutations and then generalize at the end. We obtain this time for the first permutation term of the previous Eq. $`<\stackrel{~}{p}_i\stackrel{~}{p}_j\stackrel{~}{p}_l\stackrel{~}{p}_m>`$ by the Wick theorem applied to the Gaussian 3-dimensional density field $`8`$-point correlation function $`C_2^8.C_2^6.C_2^4/4!=105`$ terms of which $`12`$ non-zero terms for the Gaussian contribution and $`48`$ non-zero terms for the connected part which interests us here. Performing some exhaustive and systematic calculations these $`48`$ terms can be written in the following condensed form

$`<\stackrel{~}{p}_i(𝐪_1,\eta _1)\stackrel{~}{p}_j(𝐪_2,\eta _2)\stackrel{~}{p}_l(𝐪_3,\eta _3)\stackrel{~}{p}_m(𝐪_4,\eta _4)>`$ $`=`$ $`4G(\eta _1)G(\eta _2)G(\eta _3)G(\eta _4)`$ (51)
$`[F_{imlj}(𝐪_1,𝐪_2,𝐪_3,𝐪_4)+F_{ijlm}(𝐪_1,𝐪_4,𝐪_2,𝐪_3)+F_{lijm}(𝐪_3,𝐪_4,𝐪_1,𝐪_2)+`$
$`F_{ljmi}(𝐪_3,𝐪_1,𝐪_2,𝐪_4)+F_{jmli}(𝐪_2,𝐪_3,𝐪_1,𝐪_4)+F_{jilm}(𝐪_2,𝐪_4,𝐪_1,𝐪_3)]`$
$`\delta ^3(𝐪_1+𝐪_2+𝐪_3+𝐪_4).`$
where $`G`$ is defined in the previous section and the general tensorial function $`F_{\alpha \beta \gamma \delta }`$ is defined as

$`F_{\alpha \beta \gamma \delta }(𝐪_i,𝐪_j,𝐪_l,𝐪_m)`$ $`=`$ $`f_{\alpha \beta \gamma \delta }(𝐪_i,𝐪_j,𝐪_l)+f_{\alpha \gamma \beta \delta }(𝐪_i,𝐪_j,𝐪_m)`$ (52)
where a general $`f_{xyzw}`$ is given by

$`f_{xyzw}(𝐪_i,𝐪_j,𝐪_l)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)P(d)\left(\frac{a_x}{a^2}+\frac{b_x}{b^2}\right)\left(\frac{d_y}{d^2}\frac{c_y}{c^2}\right)\left(\frac{b_z}{b^2}+\frac{d_z}{d^2}\right)\left(\frac{a_w}{a^2}\frac{c_w}{c^2}\right)}`$ (53)
where $`𝐚=𝐊^{}`$, $`𝐛=𝐊^{}𝐪_i`$, $`𝐜=𝐊^{}+𝐪_j`$ and $`𝐝=𝐊^{}+𝐪_j+𝐪_l`$. Again we will work with a single $`f`$ and at the end generalize the result.

Applying the small-angle approximation again we can remove the tensorial dependence of our expression and contract the $`\theta _i`$ with our vectors. We obtain the scalar $`f`$

$`f(𝐪_i,𝐪_j,𝐪_l)`$ $`=`$ $`\widehat{\theta }_{1\alpha }\widehat{\theta }_{2\beta }\widehat{\theta }_{3\gamma }\widehat{\theta }_{4\delta }f_{\alpha \beta \gamma \delta }(𝐪_i,𝐪_j,𝐪_l)`$ (54)
$``$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)P(d)\left(\frac{a_z}{a^2}+\frac{b_z}{b^2}\right)\left(\frac{d_z}{d^2}\frac{c_z}{c^2}\right)\left(\frac{b_z}{b^2}+\frac{d_z}{d^2}\right)\left(\frac{a_z}{a^2}\frac{c_z}{c^2}\right)}`$
and the scalar $`F`$

$`F(𝐪_i,𝐪_j,𝐪_l,𝐪_m)`$ $`=`$ $`f(𝐪_i,𝐪_j,𝐪_l)+f(𝐪_i,𝐪_j,𝐪_m).`$ (55)
Expanding $`f`$ in order of $`K_z^{}`$ and keeping the non-nulls terms (see previous section), i.e. the terms which are even in $`K_z^{}`$, we are left with

$`f(𝐪_i,𝐪_j,𝐪_l)`$ $`=`$ $`{\displaystyle d^3KP(a)P(b)P(c)P(d)}`$ (56)
$`[K_z^4({\displaystyle \frac{1}{a^2}}{\displaystyle \frac{1}{b^2}})({\displaystyle \frac{1}{c^2}}{\displaystyle \frac{1}{d^2}})({\displaystyle \frac{1}{d^2}}{\displaystyle \frac{1}{b^2}})({\displaystyle \frac{1}{c^2}}{\displaystyle \frac{1}{a^2}})`$
$`+K_z^2(q_{iz}^2({\displaystyle \frac{1}{b^4c^4}}{\displaystyle \frac{1}{a^2b^4c^2}}+{\displaystyle \frac{1}{a^4d^4}}{\displaystyle \frac{1}{a^2b^2d^4}}{\displaystyle \frac{1}{a^2c^2d^4}}+{\displaystyle \frac{1}{b^2c^2d^4}}{\displaystyle \frac{1}{a^2b^4d^2}}`$
$`{\displaystyle \frac{1}{a^4b^2d^2}}{\displaystyle \frac{1}{b^2c^4d^2}}+{\displaystyle \frac{1}{b^4c^2d^2}}+{\displaystyle \frac{2}{a^2b^2c^2d^2}})`$
$`+q_{jz}^2\left({\displaystyle \frac{1}{b^4c^4}}+{\displaystyle \frac{1}{a^2b^2c^4}}+{\displaystyle \frac{1}{a^2c^4d^2}}+{\displaystyle \frac{1}{b^2c^4d^2}}\right)+q_{lz}^2\left({\displaystyle \frac{1}{a^4d^4}}+{\displaystyle \frac{1}{a^2b^2d^4}}{\displaystyle \frac{1}{b^2c^2d^4}}{\displaystyle \frac{1}{b^2c^4d^2}}\right)`$
$`+q_{iz}q_{jz}\left({\displaystyle \frac{4}{b^4c^4}}{\displaystyle \frac{2}{a^2b^2c^4}}+{\displaystyle \frac{2}{a^2b^4c^2}}+{\displaystyle \frac{1}{a^4b^2c^2}}{\displaystyle \frac{2}{a^2c^2c^4}}{\displaystyle \frac{1}{b^2c^2d^4}}+{\displaystyle \frac{2}{a^2c^4d^2}}{\displaystyle \frac{1}{a^4c^2d^2}}\right)`$
$`+q_{iz}q_{lz}\left({\displaystyle \frac{2}{a^4d^4}}{\displaystyle \frac{2}{a^2c^2d^4}}{\displaystyle \frac{2}{a^2b^4d^2}}{\displaystyle \frac{1}{a^4b^2d^2}}{\displaystyle \frac{1}{b^2c^4d^2}}+{\displaystyle \frac{2}{b^4c^2d^2}}+{\displaystyle \frac{2}{a^2b^2c^2d^2}}\right)`$
$`+q_{jz}q_{lz}({\displaystyle \frac{2}{a^2c^2d^4}}{\displaystyle \frac{2}{b^2c^2d^4}}+{\displaystyle \frac{2}{a^2c^4d^2}}+{\displaystyle \frac{2}{b^2c^4d^2}}{\displaystyle \frac{1}{a^4c^2d^2}}{\displaystyle \frac{1}{b^4c^2d^2}}{\displaystyle \frac{2}{a^2b^2c^2d^2}}))+`$
$`(q_{iz}^3q_{jz}({\displaystyle \frac{1}{b^2c^2d^2}}{\displaystyle \frac{1}{b^4c^2d^2}})+q_{iz}^2q_{iz}^2({\displaystyle \frac{1}{b^4c^2}}{\displaystyle \frac{1}{b^2c^4d^2}})`$
$`+q_{iz}^2q_{jz}q_{lz}({\displaystyle \frac{2}{b^2c^2d^4}}{\displaystyle \frac{1}{b^4c^2d^2}})+q_{iz}q_{jz}^2q_{lz}({\displaystyle \frac{1}{b^2c^4d^2}})+q_{iz}q_{jz}q_{lz}^2\left({\displaystyle \frac{1}{b^2c^2d^4}}\right))].`$
We concentrate on the dominant term which can be simplified using the same method applied to the calculation of the dominant term of the power spectrum. It can be written as

$`f(𝐪_i,𝐪_j,𝐪_l)`$ $`=`$ $`{\displaystyle d^3K^{}P(a)P(b)P(c)P(d)K_z^4\left(\frac{1}{a^2}\frac{1}{b^2}\right)\left(\frac{1}{c^2}\frac{1}{d^2}\right)\left(\frac{1}{d^2}\frac{1}{b^2}\right)\left(\frac{1}{c^2}\frac{1}{a^2}\right)}.`$ (57)
To find an explicit relation between $`𝐚`$, $`𝐛`$, $`𝐜`$, $`𝐝`$ and $`𝐪_i`$, $`𝐪_j`$, $`𝐪_l`$ and $`𝐊^{}`$ we express both $`𝐊^{}`$ and $`𝐚`$, $`𝐛`$, $`𝐜`$, $`𝐝`$ in the basis $`(\widehat{e}_z,\widehat{q}_i,\widehat{q}_j)`$. We should say here that using the Limber approximation will allow only two of our initial vectors among the set of four ($`𝐪_1`$,$`𝐪_2`$,$`𝐪_3`$,$`𝐪_4`$) to be independent. Indeed all parallel components to the line of sight will be negligible and the four vectors will be inside the same plane, the one perpendicular to the line of sight. This justifies the use of the basis chosen. We stress one time more that this eliminates any sub-dominant term that could show up (see discussion in Sec. II.6). So in that basis our function $`f`$ can be expressed in the following simplified way

$`f(𝐪_i,𝐪_j,𝐪_l)`$ $`=`$ $`f(𝐪_i,𝐪_j,𝐪_l,\widehat{q}_i.\widehat{q}_j,\widehat{q}_i.\widehat{q}_l,\widehat{q}_j.\widehat{q}_l)`$ (58)
$`=`$ $`{\displaystyle 𝑑K_z^{}𝑑u𝑑v\sqrt{1\alpha ^2}P(a)P(b)P(c)P(d)K_z^4\left(\frac{1}{a^2}\frac{1}{b^2}\right)\left(\frac{1}{c^2}\frac{1}{d^2}\right)\left(\frac{1}{d^2}\frac{1}{b^2}\right)\left(\frac{1}{c^2}\frac{1}{a^2}\right)}`$
where $`\alpha =\widehat{q}_i.\widehat{q}_j`$, $`\beta =\widehat{q}_i.\widehat{q}_l`$ and $`\gamma =\widehat{q}_j.\widehat{q}_l`$. Also $`a=(K_z^2+u^2+v^2+2uv\alpha )^{1/2}`$, $`b=(K_z^2+(uq_i)^2+v^2+2(uq_i)v\alpha )^{1/2}`$, $`c=(K_z^2+(v+q_j)^2+u^2+2(v+q_j)u\alpha )^{1/2}`$, $`d=(K_z^2+(v+q_j+y)^2+(u+x)^2+2(v+q_j+y)(u+x)\alpha )^{1/2}`$ with $`y=q_l\frac{\gamma \beta \alpha }{(1\alpha ^2)}`$ and $`x=q_l\frac{\beta \gamma \alpha }{(1\alpha ^2)}`$. Our $`f`$ only depends on the norms of its arguments and on the angles between their directions. We can combine these expressions to simplify the expression (50). Following Kaiser’s method and proceeding with the Limber approximation, we obtain the expression for one of the $`24`$ possible permutation terms of the trispectrum

$`T_{dom}^{OV}(K_1,K_2,K_3,K_4)={\displaystyle \frac{1}{32\pi ^3}}{\displaystyle _0^{\eta _0}}𝑑\eta \left({\displaystyle \frac{aD\dot{D}}{D_0^2}}\right)^4{\displaystyle \frac{g(\eta )^4}{\eta ^6}}`$ (59)
$`[F(\stackrel{}{K}_1/\eta ,\stackrel{}{K}_2/\eta ,\stackrel{}{K}_3,/\eta \stackrel{}{K}_4/\eta )+F(\stackrel{}{K}_1/\eta ,\stackrel{}{K}_4/\eta ,\stackrel{}{K}_2/\eta ,\stackrel{}{K}_3/\eta )+F(\stackrel{}{K}_3/\eta ,\stackrel{}{K}_4/\eta ,\stackrel{}{K}_1/\eta ,\stackrel{}{K}_2/\eta )+`$
$`F(\stackrel{}{K}_3/\eta ,\stackrel{}{K}_1/\eta ,\stackrel{}{K}_2/\eta ,\stackrel{}{K}_4/\eta )+F(\stackrel{}{K}_2/\eta ,\stackrel{}{K}_3/\eta ,\stackrel{}{K}_1/\eta ,\stackrel{}{K}_4/\eta )+F(\stackrel{}{K}_2/\eta ,\stackrel{}{K}_4/\eta ,\stackrel{}{K}_1/\eta ,\stackrel{}{K}_3/\eta )]`$
where $`F`$ is defined by Eq. (55) and the $`f`$ by Eq. (58). Finally, by including all the permutation terms in Eq. (50) we obtain

$$T_{dom}^{OV}(K_1,K_2,K_3,K_4)=\frac{3}{8\pi ^3}_0^{\eta _0}𝑑\eta \left(\frac{aD\dot{D}}{D_0^2}\right)^4\frac{g^4(\eta )}{\eta ^6}\underset{i,j,l=1;ijl}{\overset{n=4}{}}f(\frac{\stackrel{}{K}_i}{\eta },\frac{\stackrel{}{K}_j}{\eta },\frac{\stackrel{}{K}_l}{\eta })$$
(60)
where $`f`$ is given by (58). The power of the Limber approximation was to reduce an almost impossible integration to a $`4`$-dimensional integral, which can be numerically calculated for a chosen configuration of the wavenumbers $`K`$.

#### IV.0.2 Nonlinear extension: the kSZ effect trispectrum

We follow the same approach as for the power spectrum and bispectrum and again we calculate the nonlinear extension for the dominant contribution. As we have four CDM power spectrum showing up in the expression for the trispectrum, which we need to divide among linear/nonlinear contributions from density/velocity contributions we will assume the following trispectrum statistic

$$T_{dom}^{OVNL}(K_1,K_2,K_3,K_4)=\frac{3}{8\pi ^3}_0^{\eta _0}𝑑\eta \left(\frac{aD\dot{D}}{D_0^2}\right)^4\frac{g^4(\eta )}{\eta ^6}\underset{i,j,l,k=1;ijl}{\overset{n=4}{}}\left(\frac{P^{NL}(K_k/\eta )}{P(K_k/\eta )}\right)^2f(\frac{\stackrel{}{K}_i}{\eta },\frac{\stackrel{}{K}_j}{\eta },\frac{\stackrel{}{K}_l}{\eta }).$$
(61)
#### IV.0.3 Signal-to-noise

As pointed out by Zaldarriaga Zaldarriaga00  and later by Hu Hu99 , the maximal signal-to-noise can be proven to be

$$\chi ^2\left(\frac{S}{N}\right)^2=f_{sky}\underset{L}{}\underset{\mathrm{}_1\mathrm{}_2\mathrm{}_3\mathrm{}_4}{}\frac{1}{2L+1}\frac{T_{\mathrm{}_1\mathrm{}_2\mathrm{}_3\mathrm{}_4}(L)^2}{C_\mathrm{}_1^{tot}C_\mathrm{}_2^{tot}C_\mathrm{}_3^{tot}C_\mathrm{}_4^{tot}}$$
(62)
where $`T_{\mathrm{}_1\mathrm{}_2\mathrm{}_3\mathrm{}_4}(L)`$ is one of the possible configurations of the full-sky trispectrum as defined in what follows. The covariance matrix used to obtain the Fisher matrix is calculated assuming full-sky coverage and $`\mathrm{}_1\mathrm{}_2<\mathrm{}_3\mathrm{}_4`$ by Komatsu Komatsu02  in the weakly non-Gaussian limit. If $`\mathrm{}_1\mathrm{}_2<\mathrm{}_3\mathrm{}_4`$ is not respected, the covariance would be distributed across many $`L`$s and can lead to overestimates of the signal-to-noise Hu01 . By not respecting this last constrain in Eq. (62) we are calculating an upper limit of the $`S/N`$ estimate.

Concerning the equivalence between the full-sky and flat-sky formalisms, we follow Hu’s Appendix Hu01  where it is argued that we can find a relation between the two formalisms by breaking up the trispectrum in the three possible combinations in each configuration defined by $`(\mathrm{}_1,\mathrm{}_2,\mathrm{}_3,\mathrm{}_4)`$

$$T_{\mathrm{}_1\mathrm{}_2\mathrm{}_3\mathrm{}_4}=T_{(\mathrm{}_1\mathrm{}_2)(\mathrm{}_3\mathrm{}_4)}(L_{12})+T_{(\mathrm{}_1\mathrm{}_3)(\mathrm{}_2\mathrm{}_4)}(L_{13})+T_{(\mathrm{}_1\mathrm{}_4)(\mathrm{}_3\mathrm{}_2)}(L_{14})$$
(63)
where $`L_{ij}`$ correspond to the side of the triangle of sides $`(\mathrm{}_i,\mathrm{}_j,L_{ij})`$. Each of the $`T_{(\mathrm{}_i\mathrm{}_j)(\mathrm{}_k\mathrm{}_m)}(L_{ij})`$ is then related to the flat-sky equivalent by

$`T_{(\mathrm{}_i\mathrm{}_j)(\mathrm{}_k\mathrm{}_m)}(L_{ij})`$ $`=`$ $`{\displaystyle \frac{2L+1}{4\pi }}\sqrt{(2\mathrm{}_i+1)(2\mathrm{}_j+1)(2\mathrm{}_l+1)(2\mathrm{}_m+1)}\left(\begin{array}{ccc}\mathrm{}_i& \mathrm{}_j& L_{ij}\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}\mathrm{}_k& \mathrm{}_m& L_{ij}\\ 0& 0& 0\end{array}\right)`$ (69)
$`T((K_i=\mathrm{}_i,K_j=\mathrm{}_j),(K_k=\mathrm{}_k,K_m=\mathrm{}_m)(L_{ij})).`$
As for the bispectrum, the Wigner-3j vanishes if $`\mathrm{}_i+\mathrm{}_j+L_{ij}=odd`$. So the full-sky trispectrum can only be estimated for even terms.

#### IV.0.4 Results

We are interested in the numerical evaluation of the dominant contribution of the trispectrum as it may be of cosmological interest in the near future due to the next generation of experiments. We choose the trapezoidal configuration in Fourier space $`K_1=K_2=K_3=K_4=l`$ with an angle $`\theta `$ between two consecutive sides. This will give rise to the following flat-sky trispectrum

$`T_{dom}^{OV}(\mathrm{},\theta )`$ $`=`$ $`{\displaystyle \frac{3}{2\pi ^3}}{\displaystyle _0^{\eta _0}}𝑑\eta \left({\displaystyle \frac{aD\dot{D}}{D_0^2}}\right)^4{\displaystyle \frac{g^4(\eta )}{\eta ^6}}`$ (70)
$`(f({\displaystyle \frac{\mathrm{}}{\eta }},1,ϵ,ϵ))+f({\displaystyle \frac{\mathrm{}}{\eta }},1,ϵ,ϵ))+f({\displaystyle \frac{\mathrm{}}{\eta }},ϵ,1,ϵ))`$
$`+f({\displaystyle \frac{\mathrm{}}{\eta }},ϵ,1,ϵ))+f({\displaystyle \frac{\mathrm{}}{\eta }},ϵ,ϵ,1))+f({\displaystyle \frac{\mathrm{}}{\eta }},ϵ,ϵ,1)))`$
where $`ϵ=\mathrm{cos}(\theta )`$ and $`f`$ is defined by Eq. (58). The various $`f`$ correspond to the different specific configurations due to the ordering of the vectors $`\stackrel{}{K}`$s in the sum over the function $`f`$ in Eq. (60). Note that the analytical expression of $`T_{dom}^{OV}`$ is symmetric around $`\theta =\pi /2`$ reflecting the symmetry of the configuration. The nonlinear analog follows from the previous equation.

With this configuration which depends only on the multipole $`\mathrm{}`$ and the angle $`\theta `$ between two sides one can calculate an estimate of the order of magnitude of the $`(S/N)^2(\theta )`$ per bin of $`\mathrm{}`$ were we to include all the remaining $`\mathrm{}`$ modes of the full-sky trispectrum. Indeed in the Eq. (62), we see that the number of modes contributing to the $`(S/N)^2`$ per bin of $`\mathrm{}`$ increases as $`\mathrm{}^3`$ so we will have the following relation for $`\chi ^2`$ in a bin of l

$$\frac{d\chi ^2(\theta )}{d\mathrm{}}\mathrm{}^3f_{sky}\underset{i=1}{\overset{n=3}{}}\frac{T_{\mathrm{}}(L_i)^2}{(2L_i+1)(C_{\mathrm{}}^{tot})^4}$$
(71)
where $`T_{\mathrm{}}(L_i)`$ is one of the $`3`$ possible configurations of the full-sky trispectrum for a given multipole $`\mathrm{}`$ and angle $`\theta `$ which needs to be calculated from its flat-sky counterpart (Eq. (70)) using Eqs. (63) and (69). The $`C_{\mathrm{}}^{tot}`$ is calculated as for the bispectrum (see Eq. (46)). Again, $`f_{sky}`$ represents the reduction in signal-to-noise due to incomplete sky coverage.

We show the plots for the linear and the nonlinear OV flat-sky trispectrum for $`0.95<ϵ<0.00`$, assuming both $`z_r=8`$ and $`z_r=17`$, in figure (3). As the trispectrum is symmetric around $`ϵ=0.00`$ we only plotted the negative $`ϵ`$s, but this choice was arbitrary. For $`z_r=8`$ the linear trispectrum peaks around multipole $`\mathrm{}2.10^3`$ and has a maximum amplitude of $`\mathrm{}^4T(\mathrm{})/(2\pi )2.4\times 10^{27}`$, regardless of the configuration. Choosing an earlier reionization of $`z_r=17`$ shifts the peak to $`\mathrm{}3.10^3`$ and increases the overall amplitude of the trispectrum to $`l^4T(\mathrm{})/(2\pi )1.2\times 10^{26}`$. Numerically we find for the amplitude of the trispectrum the approximate scaling dependence with reionization history $`T(\mathrm{}\mathrm{\hspace{0.17em}170})2.6\times 10^{31}x_e^4\mathrm{log}^{0.8}(1+z_r)`$, in agreement with the relation found for the power spectrum. The effect of increasing the value of angle between two consecutive sides of the trapezoid is to increase the power at small scales ($`\mathrm{}>10^4`$), such that only small scales are sensitive enough to probe different configurations. At those small scales, the power is maximum for $`ϵ=0.00`$ and then decreases as we decrease the angle. So at small angular scales the square configuration is the one contributing the most. We could have expected this behavior as the OV effect has a quite symmetric morphological signature, such that a more filamentary structure probed by the collapsed configuration of the trispectrum is not as likely.

Here we do not observe the sharp drop in power at small scales due to the Limber cancellation which enables one to speak of the trispectrum of the kSZ effect. Indeed, concerning the nonlinear trispectrum, we observe an interesting feature. Contrary to the bispectrum and similarly to the power spectrum, the trispectrum is strongly affected by the weak nonlinear enhancement due to formation of structure at small scales. This enhancement has the power to broaden the shape, to shift the peak of the effect to smaller angular scales and to increase slightly the amplitude in a configuration dependent way. So (for $`z_r=8`$/$`z_r=17`$) we measure a higher and higher maximum amplitude ($`2.47.2\times 10^{27}`$/$`1.94.8\times 10^{26}`$) at a smaller and smaller angle ($`\mathrm{}14.10^4`$) when you go from the collapsed trapezoid to the square configuration. It is worth noticing that, whereas the amplitude of the nonlinear trispectra increases, the multipoles corresponding to the peak of the effect remain the same for both reionization scenarios. This is because the nonlinearities take place at the same instant in time for both reionization scenarios, leaving an imprint at the same characteristic angular scale.

But the most important quantity is the signal-to-noise. We show in figure (4) the estimated contributions to $`\chi ^2`$ per log interval in $`\mathrm{}`$ for the linear and the nonlinear OV full-sky trispectrum ($`z_r=8`$ and $`z_r=17`$) had we included all the modes for $`ϵ=0.95`$ with no instrumental noise, Planck noise and MAP noise included in the variance MAP/PLANCK . We chose $`ϵ=0.95`$, i.e. the collapsed configuration, because it generates the highest contributions to the full-sky trispectrum near the peak of the effect, contrary to the flat-sky trispectrum for which the square configuration was the one producing the highest amplitude. The square configuration continues to provide the strongest contribution at very small angular scales, but at larger angular scales the collapsed configuration dominates. This is due to the angular averaged factors relating flat- to full-sky trispectra (see Eq. (69)). Though it is not explicit in the figure, the $`d\chi ^2/d\mathrm{ln}\mathrm{}`$ is zero when $`2\mathrm{}+L`$ is odd. Again, as for the bispectrum case, the structure in the $`d\chi ^2/d\mathrm{}`$ arises mainly from the structure of the CMB primary power spectrum at $`\mathrm{}>200`$. The continuous rise at $`\mathrm{}3000`$ for a perfect experiment is due to the fact that the signal decreases slower than the contributions to the $`C_{\mathrm{}}`$ from primary, OV and thermal SZ anisotropies up to these scales. Contrary to the common and most naive expectation, for Planck, an eventual detection is possible at arcminute scales. Indeed, the total $`S/N4.7`$ for $`z_r=17`$. Most of the contributions come from multipoles between $`\mathrm{}=12.10^3`$. This probably is the most important conclusion of this work and illustrates our predictions. This result should be taken with caution as it corresponds to a very optimistic upper limit on the $`S/N`$ (see sections Sec. II.5 and Sec. IV.0.3). Firstly, the use of the Fisher matrix formalism gives the minimum variance for our statistic (see section Sec. II.5). Secondly, the $`\chi ^2`$ method assumes that the form of the model is correct, which may not be the case. Thirdly, when calculating the covariance matrix of the Fisher matrix, two simplifying assumptions were used: that the main contribution to the covariance was Gaussian in nature and that we observed with full-sky coverage. Finally, other physical mechanisms, like for example lensing effects or the thermal SZ effect Hu01 ; Cooray01b ; Zaldarriaga00 ; Cooray:Kesden02 , and unaccounted foregrounds Tegmark99  can contribute to the trispectrum at this level and so the separability problem needs to be addressed in due time. Of course, uncertainties of roughly an order of magnitude in the modeling of the thermal SZ signal are also a source of error in our estimates. The forecasted ability of future multi-frequency experiments to remove most of the thermal SZ contributions would minimize these uncertainties and would much favor a detection. Last but not least, further progress in the implementation of optimal unbiased trispectrum estimators to probe such small scales and power is required.

## V Conclusions

Due to its strong predictive power, linear theory is a very sensitive probe of the early stages of the reionization history through the Ostriker and Vishniac effect. Analytical expressions for its correlation functions can be derived and their measurement would be of high value to our present knowledge of that still unclear epoch of the universe evolution. We have presented detailed calculations of the three Fourier statistics of interest of the OV effect, the power spectrum, the bispectrum and the trispectrum. For that purpose we have developed a new technique that allows one to obtain their dominant contributions under the Limber approximation framework. This method is applicable to the derivation of any statistics involving correlations among vector-like effects. It illustrates what was expected under statistical homogeneity and isotropy assumptions and the vector and small scale nature of the OV effect. We also evaluated numerically, as a function of scale and for a specific configuration (equilateral for the bispectrum and trapezoidal for the trispectrum), these statistics for a flat $`\mathrm{\Lambda }`$CDM cosmology and two reionization scenarios. The first one is based on our pre-WMAP knowledge ($`z_r=8`$) and the second one takes into account the high values for the electron optical depth measured by WMAP ($`z_r=17`$). We numerically obtained approximate scaling relations for the amplitude of the OV statistics on the reionization history considered and found that the dependence is stronger on the ionization fraction than on the redshift of reionization. We have also studied their detectability in view of future satellite experiments. The alternation of dominant/subdominant/dominant higher order correlation functions was numerically shown. While the bispectrum is undetectable even by a perfect multi-frequency experiment capable of subtracting the thermal SZ contributions, the trispectrum could be measured by Planck or by interferometer experiments targeting arcminute scales with high sensitivity and for a sufficiently long period of time. This provides a unique signal distinguishing the OV effect from other non vector-like secondary anisotropies and could be useful when trying to separate different physical mechanisms imprinting themselves on these measurable statistics. One should bear in mind that despite this useful characteristic signature, our results are quite optimistic, although encouraging, as they rely on various analytically helpful idealized assumptions, as described previously. Also, other contributions to the signal are to be expected at the arcminute level and thus further study of CMB small-scale secondary anisotropies and foregrounds contributions to the trispectrum is required and much justified. In order to obtain an upper limit on the possible kSZ contributions, we also extended our calculations to the mildly nonlinear regime. We found that, contrary to the bispectrum, there is a noticeable enhancement of the contributions of the trispectrum in a morphological dependent way and that this enhancement reflects itself on the calculations of the signal-to-noise. Hence nonlinearities are expected to enhance the even non-Gaussian signals produced by the OV effect and further complicate its disentanglement from inevitable model-dependent nonlinear effects arising from structure formation. Conversely, having a template of the OV effect can help in extracting the nonlinear contributions to reionization at those angular scales, providing another possible window on the complex physics associated with reionization.

###### Acknowledgements.

I am much indebted to Pedro G. Ferreira whose suggestions, comments and constant support have been invaluable during the preparation of this work. I thank also Asantha Cooray for very helpful discussions and comments and for the thermal SZ data used to facilitate direct comparison with his results, Evan Scannapieco for helping to improve the manuscript and to correct a fundamental error in the calculations, Nabila Aghanim for her support and the thermal SZ data used throughout the article, Andrew Jaffe for his comments and Martin Kunz for his encouragement. This work was supported by the Fundação para a Ciência e a Tecnologia under the reference PRAXIS XXI/BD/21249/99.

## Appendix A Limber Approximation

We review the fundamental steps of the limber approximation as used in the text. The Limber Eq. Limber54  describes the two-point statistics of a field which is the two-dimensional projection on the sky of a three-dimensional field whose statistical properties vary slowly along the line of sight. The Fourier space analog of this result was calculated by Kaiser Kaiser92 . Later it was further extended from the flat-sky approximation to an all-sky approach for spatially flat cosmologies by Hu and White Hu:White96  and to open cosmologies by Hu Hu99 . Buchalter *et. al.* Buchalter00  derived the Fourier space analog of the Limber’s equation for the bispectrum.

Here we review the bispectrum derivation of Buchalter *et. al.* and generalize it to the trispectrum. This applies as well to higher-order statistics. The error introduced by the Limber approximation is inferior to $`1\%`$ for effects with $`l>200`$ Cooray:Hu00 ; Pollo01 .

Suppose we observe the projection along the line of sight of a three-dimensional statistically homogeneous and isotropic random field $`f`$

$$p(\stackrel{}{\theta })=_0^{\mathrm{}}𝑑\eta q(\eta )f(\eta \widehat{\theta })$$
(72)
where $`\stackrel{}{\theta }=(\theta _1,\theta _2,0)`$ and $`\widehat{\theta }=(\theta _1,\theta _2,1)`$. We propose ourselves to find a relation between the spatial bispectrum $`B_p(K_1,K_2,K_3)`$ of $`p`$ and the spatial bispectrum $`B_f(K_1,K_2,K_3)`$ of $`f`$. Following Kaiser method, we consider $`p`$ to be the sum of the contributions from narrow shells with a width $`\mathrm{\Delta }\eta `$ much bigger than the relevant wavelength, that is $`\theta <<\mathrm{\Delta }\eta /\eta <<1`$. This choice allows us to look at fluctuations on scales much less than the characteristic scale over which $`q`$ varies and to assume that contributions from different shells are statistically independent. We then calculate the contributions to the bispectrum from each of the shells. At the end we can sum the power for all the shells relying on their statistical independence.

Assuming that $`q`$ varies very little along the shell and that the section of the shell is plane-symmetric the contribution from the shell of width $`\mathrm{\Delta }\eta `$ centered at $`\eta _0`$ is

$$\mathrm{\Delta }p(\stackrel{}{\theta })=q(\eta _0)_{\eta _0\mathrm{\Delta }\eta /2}^{\eta _0+\mathrm{\Delta }\eta /2}𝑑\eta f(\eta _0\theta _1,\eta _0\theta _2,\eta )$$
(73)
Decomposing the fields in Fourier space, the spectrum of $`\mathrm{\Delta }p`$ is

$$\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K})=q(\eta _0)\frac{d^3k}{(2\pi )^3}\stackrel{~}{f}(𝐤)d^2\theta e^{i(\eta _0\stackrel{}{k}\stackrel{}{K}).\stackrel{}{\theta }}_{\eta _0\mathrm{\Delta }\eta /2}^{\eta _0+\mathrm{\Delta }\eta /2}𝑑\eta e^{ik_z\eta }$$
(74)
where $`\stackrel{}{K}=(K_x,K_y,0)`$ and $`𝐤=(k_x,k_y,k_z)=(\stackrel{}{k},k_z)`$. We can perform the $`d^2\theta `$ integral which is $`(2\pi )^2/\eta _0^2\delta _D^2(\stackrel{}{k}\stackrel{}{K}/\eta _0)`$ and the time integral which is $`\mathrm{\Delta }\eta j_0(k_z\mathrm{\Delta }\eta /2)`$ to yield

$$\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K})=\frac{\mathrm{\Delta }\eta q(\eta _0)}{\eta _0^2}\frac{dk_z}{(2\pi )}\stackrel{~}{f}(\frac{K_x}{\eta _0},\frac{K_y}{\eta _0},k_z)j_0(k_z\mathrm{\Delta }\eta /2)$$
(75)
Using $`<\stackrel{~}{f}(𝐤_1)\stackrel{~}{f}(𝐤_2)\stackrel{~}{f}(𝐤_3)>=(2\pi )^3B_f(k_1,k_2,k_3)\delta _D^3(𝐤_1+𝐤_2+𝐤_3)`$ the three-point spectrum is

$`<\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_1)\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_2)\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_3)>`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }\eta ^3q^3(\eta _0)}{\eta _0^6}}{\displaystyle 𝑑k_{1z}𝑑k_{2z}𝑑k_{3z}B_f(\sqrt{\frac{K_1^2}{\eta _0^2}+k_{1z}^2},\sqrt{\frac{K_2^2}{\eta _0^2}+k_{2z}^2},\sqrt{\frac{K_3^2}{\eta _0^2}+k_{3z}^2})}`$ (76)
$`j_0(k_{1z}\mathrm{\Delta }\eta /2)j_0(k_{2z}\mathrm{\Delta }\eta /2)j_0(k_{3z}\mathrm{\Delta }\eta /2)\delta _D^2(\stackrel{}{K}_1/\eta _0+\stackrel{}{K}_2/\eta _0+\stackrel{}{K}_3/\eta _0)\delta _D(k_{1z}+k_{2z}+k_{3z})`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }\eta ^3q^3(\eta _0)}{\eta _0^4}}{\displaystyle 𝑑k_{1z}𝑑k_{2z}B_f(\sqrt{\frac{K_1^2}{\eta _0^2}+k_{1z}^2},\sqrt{\frac{K_2^2}{\eta _0^2}+k_{2z}^2},\sqrt{\frac{K_3^2}{\eta _0^2}+(k_{1z}+k_{2z})^2})}`$
$`j_0(k_{1z}\mathrm{\Delta }\eta /2)j_0(k_{2z}\mathrm{\Delta }\eta /2)j_0((k_{1z}k_{2z})\mathrm{\Delta }\eta /2)\delta _D^2(\stackrel{}{K}_1+\stackrel{}{K}_2+\stackrel{}{K}_3)`$
The important simplification comes from the fact that the major contribution from the first two Bessel functions comes from $`k_{1z}<1/\mathrm{\Delta }\eta `$ and $`k_{2z}<1/\mathrm{\Delta }\eta `$. But by assumption $`K_1/\eta _0>>1/\mathrm{\Delta }\eta `$, $`K_2/\eta _0>>1/\mathrm{\Delta }\eta `$ and $`K_3/\eta _0>>1/\mathrm{\Delta }\eta `$. Therefore, $`k_{1z}<<K_1/\eta _0`$, $`k_{2z}<<K_2/\eta _0`$ and $`k_{1z}+k_{2z}<<K_3/\eta _0`$. So we can neglect all Fourier modes parallel to the line of sight. To a very good approximation $`B_f(\sqrt{\frac{K_1^2}{\eta _0^2}+k_{1z}^2},\sqrt{\frac{K_2^2}{\eta _0^2}+k_{2z}^2},\sqrt{\frac{K_3^2}{\eta _0^2}+(k_{1z}+k_{2z})^2})B_f(\frac{K_1}{\eta _0},\frac{K_2}{\eta _0},\frac{K_3}{\eta _0})`$. The integration of the Bessel functions gives

$$𝑑u𝑑vj_0(u)j_0(v)j_0(u+v)=\pi ^2$$
(77)
We obtain finally

$`<\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_1)\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_2)\stackrel{~}{\mathrm{\Delta }p}(\stackrel{}{K}_3)>`$ $`=`$ $`4\pi ^2{\displaystyle \frac{\mathrm{\Delta }\eta q^3(\eta _0)}{\eta _0^4}}B_f({\displaystyle \frac{K_1}{\eta _0}},{\displaystyle \frac{K_2}{\eta _0}},{\displaystyle \frac{K_3}{\eta _0}})\delta _D^2(\stackrel{}{K}_1+\stackrel{}{K}_2+\stackrel{}{K}_3)`$ (78)
Summing over the shells, and using $`<\stackrel{~}{p}(\stackrel{}{K}_1)\stackrel{~}{p}(\stackrel{}{K}_2)\stackrel{~}{p}(\stackrel{}{K}_3)>=(2\pi )^2B_p(K_1,K_2,K_3)\delta _D^2(\stackrel{}{K}_1+\stackrel{}{K}_2+\stackrel{}{K}_3)`$

$`B_p(K_1,K_2,K_3)`$ $`=`$ $`{\displaystyle 𝑑\eta \frac{q^3(\eta )}{\eta ^4}B_f(\frac{K_1}{\eta },\frac{K_2}{\eta },\frac{K_3}{\eta })}`$ (79)
The exact same reasoning can be applied to the calculation of the trispectrum. This time we obtain for the same projection

$`T_p(K_1,K_2,K_3,K_4)`$ $`=`$ $`{\displaystyle 𝑑\eta \frac{q^4(\eta )}{\eta ^6}T_f(\frac{K_1}{\eta },\frac{K_2}{\eta },\frac{K_3}{\eta },\frac{K_4}{\eta })}`$ (80)