{'fraction_non_alphanumeric': 0.03619197482297404, 'fraction_numerical': 0.015718157181571817, 'mean_word_length': 4.30034287832453, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Bispectrum and the Trispectrum of the Ostriker and Vishniac Effect

We present analytical expressions for the Fourier analog of the CMB three-point and four-point correlation functions, the spatial bispectrum and trispectrum, of the Ostriker and Vishniac effect in the linear and mildly nonlinear regime. Through this systematic study, we illustrate a technique to tackle the calculation of such statistics making use of the effects of its small-angle and vector-like properties through the Limber approximation. Finally we discuss its configuration dependence and detectability in the context of Gaussian theories for the currently favored flat Lambda CDM cosmology.


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 [1, 2, 3, 4], 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 [5,6] and constrains possible non-Gaussian primordial contributions to the bispectrum and the trispectrum on degree and sub-degree angular scales using actual data [7,8,9,10]. 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 [11] was found to be the dominant linear secondary contribution to the CMB anisotropies below the Silk-damping scales at the arcminute level [12]. 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 [12,13,14] 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 A. In Sec. II B we review the theory of the OV effect and in Sec. II C we discuss the basic statistical properties of a general field through its n-point functions. In Sec. II D 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 E, we introduce the standard formalism of the signal-to-noise calculation by means of the Fisher matrix. For general consistency, in Sec. II F 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.

A. Cosmological model
We work in the context of the adiabatic cold dark matter (CDM) family of models. In units of the critical density, Ω 0 is the contribution to the non-relativistic-matter density, Ω b is the contribution to the baryonic matter density, Ω Λ is the contribution of the cosmological constant and H 0 = 100 h 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 theṅ where the over-dot denotes a derivative with respect to time. The scale factor is chosen such that a 0 H 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) = 1 2 z 0 dz ′ E(z ′ ) and at a time t(z) = 1 H0 ∞ z dz ′ (1+z ′ )E(z ′ ) . The conformal time is then obtained from dη = dt/a, such that the comoving distance to the horizon is the conformal time today, cη 0 = w(∞).
If the CDM density contrast at comoving position w at time t is δ( w, t), then the power spectrum P (k, t) is defined by the expectation value over all realizations <δ( k, t)δ * ( k ′ , t) >= (2π) 3 δ 3 D ( k − k ′ )P (k, t) where δ 3 D is the Dirac delta function. In linear theory, δ( w, t) = δ 0 ( w)D(t)/D(t 0 ), where t 0 is the age of the Universe, δ 0 ( w) ≡ δ( w, t 0 ), and the growth factor, as a function of redshift, is [15] The power spectrum is given by where T (q) is the CDM transfer function, k p = ka 0 = kH 0 /2c is the physical wave number with our conventions and Γ, the shape parameter, is defined as [16] Γ ≃ Ω 0 (h/0.5) exp(−Ω b − Ω b /Ω 0 ). As we chose a 0 H 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. [17] fitting formulae for CDM models instead of the improved version of Eisenstein and Hu [18], to facilitate comparison with previous work.
For δ H , we take the fits to the Cosmic Background Explorer given in [19]. In linear theory, the continuity equation relates the Fourier components of the velocity field and the density field A useful relation [15] isḊ Although we maintain generality in the derivations, we illustrate our results with a flat model, the ΛCDM model. The parameters for this model are Ω 0 = 0.35, Ω b = 0.05, Ω Λ = 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 ∆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 η that appear in the equations where S(x) = sinh x in an open universe and S(x) = sin x in a closed universe.

B. 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 [20,21]). 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 8 ≤ z r ≤ 30. Measurements of the CMB anisotropies on sub-degree scales [22,23,24] have been used to put an upper bound on the reionization redsfhit of z r ≃ 30 [25,26]. Very recently, using polarization and temperature anisotropies of the CMB, WMAP has placed a fairly model-independent constrain [27] on the optical depth to electron scattering of τ = 0.17 ± 0.04 at 68% CL which translates into z r = 17 ± 3 for instant reionization [28,29]. Interestingly, the measurement of an increase of the neutral fractions with redshift in high-z quasar spectra [30,31] and a first detection of the 'Gunn-Peterson trough' [32] in a quasar spectra at z = 6.28 by SLOAN [33] point to a reionization redshift of z r ≃ 6 [34], 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-α 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 [35,36,37].
Though extensive analytical (for a complete derivation see [12,13]) 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 [38,39], the thermal SZ effect and its kinetic analog [40,41,42,43,44,45,46,47], the Inhomogeneous Reionization [48,49,50,51,52] and the OV effect on smaller scales [14,53,54]. Enlightening comparative studies between different effects can also be found [55].
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 [56,57,58,59,60] 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 [59] 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 [11,14]. 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 ΛCDM models, and should produce µK anisotropies. Its contributions to the temperature fluctuations along the line of sight can be written in the manner of Jaffe and Kamionkowski (JK) [54] ∆T where p(θη, η) ≡ v(θη, η)δ(θη, η) and g is the visibility function given by 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 0 g(η)dη = 1 − e −τr where τ 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 τ (η) = η 0 c dη ′ σ T n e (η ′ ). Also δ(θη, η) and v(θη, η) are the density contrast and bulk velocity along the line of sight,n e (η) is the mean electron density given byn e (η) = Ω b ρ c x e (η) (1 + z) 3 /m p , σ 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 Following JK [54], we choose a coordinate system such thatθ represents a three-dimensional unit vector along the line of sight and θ refers to a two-dimensional unit vector in the plane perpendicular to the line of sight. So we will have θ = (θ 1 , θ 2 , 0) andθ = (θ 1 , θ 2 , 1 − θ 2 1 − θ 2 2 ) ≃ (θ 1 , θ 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 spacẽ where q ≡ (q x , q y , q z ), K ≡ (K x , K y , 0) and is the Fourier transform of p(θη, η) (see JK [54]). D andḊ depend on η. We made use of the continuity equation in Fourier space (5) and of the linear evolution of the density field.

C. 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 p with zero mean, its power spectrum P , bispectrum B and trispectrum T are defined by the following equations in the appropriate Fourier convention: 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.

D. 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 [61]. Relying on this theory, all expectation values of an odd product of an isotropic 3-dimensional vector field p(q) with q ≡ (q x , q y , q z ) must be proportional to at least one of the q 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 [62], 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.

E. 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 χ 2 statistics (as in [56,63]). To do so, one needs to calculate the Fisher information matrix F ij (for a good review see [64]). If we think of the data x as a random variable with a likelihood function L(x; θ) where θ is a vector of model parameters, the Fisher information matrix is defined as By a very powerful theorem, called the Cramer-Rao inequality, it was shown [65,66] 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 [56,57,63]. Therefore, we assume that the form of our model θ (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 σ 2 (A) = (F −1 ) AA and we define the χ 2 statistics as 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 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 θ f whm and pixel noise σ p = s/ √ 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 ≡ (σ p θ f whm /T 0 ) 2 , in order to have a measure which is independent of the pixel-size [67,68]. 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 noise ℓ [57,67,68] is where σ b , in radians, is the width of the beam if we assume it has a Gaussian profile. It is related to θ f whm , in arcminutes, by σ b = √ 8ln2 θ f whm × π/10800. Note that if an experiment maps the full sky and then a fraction 1 − f sky is subtracted, one should not multiply w −1 by f sky (case of MAP and Planck). Hence we can estimate C noise ℓ for any experiment with characteristic f sky , θ f whm , s and t pix .
For the precise cases of MAP (renamed WMAP recently) and Planck, for which we used the specifications in table I, we need to take into account their multi-frequency coverage with different characteristics. The C noise ℓ is then defined as [57] 1 C noise where the sum runs over all channels of the experiment and ν is the frequency of the channel.  F. Power spectrum of the OV effect

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 ∆T /T wherep is defined as in Eq. (12). Many authors have derived the expression for this statistics [12,13,14,54,62]. 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 p 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 <p ipj > are symmetric due to statistical homogeneity, we only consider the first one and multiply the result by 2. Using Eq. (12) forp and the Wick theorem for the Gaussian 3-dimensional density field correlation function which states where P (k) is the power spectrum of density perturbations, we obtain two non-zero terms for <p ipj > which can be written as with the time dependence functions gathered in G(η) = iaDḊ 2D 2 0 and the general tensorial functions F αβ given by where a = K ′ and b = q i − K ′ . 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 αβ . Indeed, in the small-sky approximation, for which the unit vectorθ ≃ (0, 0, 1), we can contract theθ's of expression (18) with the vectors a and b of the last expression (21) such that we are left with the line of sight components of a and b. We can thus define a new scalar function F such that The interesting step that follows is to expand this function in Expressions of this type will occur in the following and they illustrate the previous discussion in Sec. II D. 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 ′2 z ) and from the last term (with K ′0 z ). By applying the Limber approximation we can also infer that the dominant contribution has to come from the K ′2 z 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 D 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 where we have performed a spherical coordinate transform such that µ =q i .K ′ , a = y 1 q i and b = y 2 q i with y 2 = 1 + y 2 1 − 2y 1 µ. 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 ′ 1 − µ 2 . This assumption preserves our dominant term but suppresses any subdominant term that could naturally arise when calculating the integral. As a consequence, in Eq. (23), the dominant term (in q 0 iz ) may contain hidden contributions to the subdominant term (in q 2 iz ). 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 2 z and b 2 z give identical contributions to the integral, such that if we calculate twice the integral in a 2 z 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 2 iz 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 where 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 ΛCDM model assuming z r = 8 and z r = 17. The correspondence between full-sky multipole moments C ℓ and the flat-sky Fourier space P (K) is straightforward C ℓ = P (K = ℓ). 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 ℓ≃ 500 ≃ 7.5 × 10 −18 x 2 e log 0.4 (1 + z r ). We advise the reader to consult JK paper [54] 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.

Nonlinear extension: the kSZ effect power spectrum
This extension was calculated previously [69] 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 [69] among others [41,55], 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 [69] approach, we replace the linear density power spectrum with its nonlinear analogue but leave the contribution from the velocity power spectrum the same where P N L 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 [69]. We follow Hamilton et. al. [70] who presented a scaling relation for the correlation function in the nonlinear regime that was generalized to its Fourier analog by Peacock and Dodds [71]. The basic hypothesis is that nonlinear fluctuations on a scale k arise from linear fluctuations on a larger scale so that there is a function relating the nonlinear and linear power spectra at these two scales which can be fit to simulations. We use the [72] proposed expression for f NL . For stronger nonlinearities other corrections are necessary. See for example [73]. One also needs the relation ∆ 2 δ (k) = k 3 /(2π 2 )P (k). As discussed in W.Hu [69], 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 ℓ ≥ 10 4 [69,74]. 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.

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 where the five permutations are with respect to the ordering (q 1 , q 2 , q 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 <p ipjpl >. At the end we generalize the results to include the total six permutations. Using expression (12) forp and the Wick theorem for the Gaussian

3-dimensional density field 6-point correlation function
2 .C 4 2 /3! = 15 terms for <p ipjpl > of which 8 are non-zero. After some simple calculations, these 8 terms can be condensed into with the time dependence function G(η) = iaDḊ 2D 2 0 and the general tensorial functions F αβγ given by where a = K ′ , b = K ′ − q i and c = K ′ + q j . We concentrate on F and combine all the terms in Eq. (31) when appropriate.
Again using the small-angle approximation for whichθ ≃ (0, 0, 1), we can contract theθ's with the vectors a, b and c such that we are left with the line of sight components of a, b and c. We can thus define a new scalar function F such that Expanding this expression in different orders in K ′ z will give us the different levels of contributions under the Limber approximation to the bispectrum Again, and as detailed before in Sec. II F, 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 ′2 z ) and the last term (with K ′0 z ). The dominant contribution will come from the K ′2 z 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 D, 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 ′2 z term which can be written as where and 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θ i ≃ (0, 0, 1) to write 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 where h(η) = g aDḊ D 2 0 and f ( q i , q j ) = f 1 ( q i , q j ) + f 2 ( q i , q j ). In order to simplify the functions f 1 and f 2 , we need to find an explicit relation between a, b, c and q i , q j and K ′ . To do so, we express both K ′ and a, b, c in the basis (ê z ,q i ,q j ) where we remind the reader thatq α ≃ (q α⊥ , 0) under the Limber approximation. Again, as for the power spectrum (see discussion in Sec. II F), this suppresses any sub-dominant contribution that could arise in the process but leaves our dominant contribution intact. We obtain where a = (K ′2 z +u 2 +v 2 +2uvµ) 1/2 , b = (K ′2 z +(u−q i ) 2 +v 2 +2(u−q i )vµ) 1/2 , c = (K ′2 z +(v+q j ) 2 +u 2 +2(v+q j )uµ) 1/2 and µ =q i .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 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.

Nonlinear extension: the kSZ effect bispectrum
We follow the same approach as for the power spectrum in Sec. II F 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

The Signal-to-noise
To define the χ 2 statistics of Eq. (15), we need to calculate the likelihood L of observing the bispectrum elements B β ≡ B ℓ1ℓ2ℓ3 given the true parameters p, and calculate the Fisher matrix as defined in Eq. (14) Assuming that the likelihood is Gaussian, we follow Cooray and Hu [57] approach to calculate the χ 2 statistics where B ℓ1ℓ2ℓ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 σ 2 ℓ1ℓ2ℓ3 comes from the covariance matrix of the angular averaged bispectrum assuming a nearly Gaussian bispectrum and full-sky coverage [56,76] where C 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 ( The thermal SZ effect was taken from [75] and was calculated semi-analytically at 30 GHZ for a normalization factor σ 8 = 0.9. We include the linear OV effect contribution (see Sec. II F) 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 [77] and we will not consider any foreground. Indeed, for the case of MAP and Planck, studies indicate that the total C ℓ should increase by 10% maximum [78]. 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. (16)(17) applied to MAP and Planck.
We will study the contributions to the χ 2 per log interval in ℓ. 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 [57]. We calculate the total S/N for each experiment by integrating over the multipole ℓ.

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 The σ 2 ℓ is calculated using Eq. (45) for all the ℓs equal and B(ℓ) 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 l ≃ 500 with an amplitude of l 3 B(ℓ)/(2π) ≃ 1.2 × 10 −23 whereas for z r = 17 the peak takes place at a higher multipole of ℓ ≃ 700 with a stronger amplitude of ℓ 3 B(ℓ)/(2π) ≃ 3 × 10 −23 , as expected. We find numerically for the amplitude of the bispectrum the approximate scaling dependence with reionization history B(ℓ ≃ 50) ≃ 1.7 × 10 −29 x 3 e log (1 + z r ). The dependence with z r is neither in agreement with the one obtained for the power spectrum (see Sec. II F) 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. (65)). 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 [57], 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 [57], 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 χ 2 per log interval in ℓ of the linear full-sky bispectrum had we included all the modes for MAP, Planck [1] and no instrumental noise. We show the results for z r = 8 and z r = 17. The experimental specifications can be found in table I. We found no necessity of plotting the corresponding dχ 2 /dℓ for the nonlinear bispectra due to its similarity to the linear bispectra. The structure in the dχ 2 /dℓ arises mainly from the structure of the CMB primary power spectrum at ℓ > 200. We point out that the continuous rise at ℓ ≃ 3000 for a perfect experiment is due to the fact that the signal decreases slower than the contributions to the C ℓ from primary, linear OV and thermal SZ anisotropies up to these scales. Though it is not explicit in the figure, the dχ 2 /dℓ is zero when 3ℓ 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/N ∼ 4.4 × 10 −3 for z r = 8 and S/N ∼ 1.4 × 10 −2 for z r = 17. These values can be compared to the value S/N ∼ 1.7 obtained by Cooray and Hu [57] for the bispectrum generated by the coupling ISW-SZ N L -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/N ∼ 0.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 E and Sec. III 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.

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 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. <p ipjplpm > by the Wick theorem applied to the Gaussian 3-dimensional density field 8-point correlation function C 8 2 .C 6 2 .C 4 2 /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 where G is defined in the previous section and the general tensorial function F αβγδ is defined as where a general f xyzw is given by where a = K ′ , b = K ′ − q i , c = K ′ + q j and d = K ′ + q j + q 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 θ i with our vectors. We obtain the scalar f and the scalar F Expanding f in order of K ′ z and keeping the non-nulls terms (see previous section), i.e. the terms which are even in 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 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.

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 3. Signal-to-noise As pointed out by Zaldarriaga [79] and later by Hu [69], the maximal signal-to-noise can be proven to be where T ℓ1ℓ2ℓ3ℓ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 ℓ 1 ≤ ℓ 2 < ℓ 3 ≤ ℓ 4 by Komatsu [60] in the weakly non-Gaussian limit. If ℓ 1 ≤ ℓ 2 < ℓ 3 ≤ ℓ 4 is not respected, the covariance would be distributed across many Ls and can lead to overestimates of the signal-to-noise [5]. 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 [5] 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 (ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 ) where L ij correspond to the side of the triangle of sides (ℓ i , ℓ j , L ij ). Each of the T (ℓiℓj )(ℓ k ℓm) (L ij ) is then related to the flat-sky equivalent by As for the bispectrum, the Wigner-3j vanishes if ℓ i + ℓ j + L ij = odd. So the full-sky trispectrum can only be estimated for even terms.

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 θ between two consecutive sides. This will give rise to the following flat-sky trispectrum where ǫ = cos(θ) and f is defined by Eq. (58). The various f correspond to the different specific configurations due to the ordering of the vectors Ks in the sum over the function f in Eq. (60). Note that the analytical expression of T OV dom is symmetric around θ = π/2 reflecting the symmetry of the configuration. The nonlinear analog follows from the previous equation.
With this configuration which depends only on the multipole ℓ and the angle θ between two sides one can calculate an estimate of the order of magnitude of the (S/N ) 2 (θ) per bin of ℓ were we to include all the remaining ℓ 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 ℓ increases as ℓ 3 so we will have the following relation for χ 2 in a bin of l where T ℓ (L i ) is one of the 3 possible configurations of the full-sky trispectrum for a given multipole ℓ and angle θ which needs to be calculated from its flat-sky counterpart (Eq. (65)) using Eqs. (63) and (64). The C 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 ℓ ≃ 2.10 3 and has a maximum amplitude of ℓ 4 T (ℓ)/(2π) ≃ 2.4 × 10 −27 , regardless of the configuration. Choosing an earlier reionization of z r = 17 shifts the peak to ℓ ≃ 3.10 3 and increases the overall amplitude of the trispectrum to l 4 T (ℓ)/(2π) ≃ 1.2 × 10 −26 . Numerically we find for the amplitude of the trispectrum the approximate scaling dependence with reionization history T (ℓ ≃ 170) ≃ 2.6 × 10 −31 x 4 e 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 (ℓ > 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.4 − 7.2 × 10 −27 /≃ 1.9 − 4.8 × 10 −26 ) at a smaller and smaller angle (ℓ ≃ 1 − 4.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 χ 2 per log interval in ℓ 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 [1]. 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. (64)). Though it is not explicit in the figure, the dχ 2 /d ln ℓ is zero when 2ℓ + L is odd. Again, as for the bispectrum case, the structure in the dχ 2 /dℓ arises mainly from the structure of the CMB primary power spectrum at ℓ > 200. The continuous rise at ℓ ≃ 3000 for a perfect experiment is due to the fact that the signal decreases slower than the contributions to the C ℓ 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/N ∼ 4.7 for z r = 17. Most of the contributions come from multipoles between ℓ = 1 − 2.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 E and Sec. IV 3). Firstly, the use of the Fisher matrix formalism gives the minimum variance for our statistic (see section Sec. II E). Secondly, the χ 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 [5,59,79,80], and unaccounted foregrounds [78] 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 Λ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. statistical properties vary slowly along the line of sight. The Fourier space analog of this result was calculated by Kaiser [82]. Later it was further extended from the flat-sky approximation to an all-sky approach for spatially flat cosmologies by Hu and White [53] and to open cosmologies by Hu [69]. Buchalter et. al. [83] 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 [57,84].
Suppose we observe the projection along the line of sight of a three-dimensional statistically homogeneous and isotropic random field f where θ = (θ 1 , θ 2 , 0) andθ = (θ 1 , θ 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 ∆η much bigger than the relevant wavelength, that is θ << ∆η/η << 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.