Non-Gaussianity and the Cosmic Microwave Background Anisotropies

We review in a pedagogical way the present status of the impact of non-Gaussianity (NG) on the Cosmic Microwave Background (CMB) anisotropies. We first show how to set the initial conditions at second-order for the (gauge invariant) CMB anisotropies when some primordial NG is present. However, there are many sources of NG in CMB anisotropies, beyond the primordial one, which can contaminate the primordial signal. We mainly focus on the NG generated from the post-inflationary evolution of the CMB anisotropies at second-order in perturbation theory at large and small angular scales, such as the ones generated at the recombination epoch. We show how to derive the equations to study the second-order CMB anisotropies and provide analytical computations to evaluate their contamination to primordial NG (complemented with numerical examples). We also offer a brief summary of other secondary effects. This review requires basic knowledge of the theory of cosmological perturbations at the linear level.


Introduction
Cosmic Microwave Background (CMB) anisotropies play a special role in cosmology, as they allow an accurate determination of cosmological parameters and may provide a unique probe of the physics of the early universe and in particular of the processes that gave origin to the primordial perturbations.
Cosmological inflation [1] is nowadays considered the dominant paradigm for the generation of the initial seeds for structure formation. In the inflationary picture, the primordial cosmological perturbations are created from quantum fluctuations "redshifted" out of the horizon during an early period of accelerated expansion of the universe, where they remain "frozen". They are observable through CMB temperature anisotropies (and polarization) and the large-scale clustering properties of the matter distribution in the Universe.
This picture has recently received further spectacular confirmations from the results of the Wilkinson Microwave Anisotropy Probe (WMAP) five year set of data [2]. Since the observed cosmological perturbations are of the order of 10 −5 , one might think that first-order perturbation theory will be adequate for all comparisons with observations. This might not be the case, though. Present [2,3] and future experiments [4] may be sensitive to the non-linearities of the cosmological perturbations at the level of secondor higher-order perturbation theory. The detection of these non-linearities through the non-Gaussianity (NG) in the CMB [5] has become one of the primary experimental targets.
There is one fundamental reason why a positive detection of NG is so relevant: it might help in discriminating among the various mechanisms for the generation of the cosmological perturbations. Indeed, various models of inflation, firmly rooted in modern particle physics theory, predict a significant amount of primordial NG generated either during or immediately after inflation when the comoving curvature perturbation becomes constant on super-horizon scales [5]. While standard single-field models of slowroll inflation [6] and -in general -two (multi)-field [7] models of inflation predict a tiny level of NG, "curvaton"-type models, in which a significant contribution to the curvature perturbation is generated after the end of slow-roll inflation by the perturbation in a field which has a negligible effect on inflation, may predict a high level of NG [8]. Alternatives to the curvaton model are models where a curvature perturbation mode is generated by an inhomogeneity in the decay rate [9,10], the mass [11] or the interaction rate [12] of the particles responsible for the reheating after inflation. Other opportunities for generating the curvature perturbations occur at the end of inflation [13], during preheating [14], and at a phase-transition producing cosmic strings [15]. Also, within single-field models of inflation, a high level of NG can be generated breaking the standard conditions of canonical kinetic terms and initially vaccum states: e.g. this is the case of Dirac-Born-Infeld (DBI) models of inflation [16], and initially excited states respectively [17]. For every scenario there exists a well defined prediction for the strength of NG and its shape [18,19] as a function of the parameters.
Statistics like the bispectrum and the trispectrum of the CMB can then be used to assess the level of primordial NG (and possibly its shape) on various cosmological scales and to discriminate it from the one induced by secondary anisotropies and systematic effects [5,20,21,22]. A positive detection of a primordial NG in the CMB at some level might therefore confirm and/or rule out a whole class of mechanisms by which the cosmological perturbations have been generated.
Despite the importance of evaluating the impact of primordial NG in a crucial observable like the CMB anisotropy, the vast majority of the literature has been devoted to the computation of the bispectrum of either the comovig curvature perturbation or the gravitational potential on large scales within given inflationary models. These, however, are not the physical quantities which are observed. One should instead provide a full prediction for the second-order radiation transfer function. A preliminary step towards this goal has been taken in Ref. [23] (see also [24,25,26,27,28,29]) where the full second-order radiation transfer function for the CMB anisotropies on large angular scales in a flat universe filled with matter and cosmological constant was computed, including the second-order generalization of the Sachs-Wolfe effect, both the early and late Integrated Sachs-Wolfe (ISW) effects and the contribution of the second-order tensor modes. 1 .
There are many sources of NG in CMB anisotropies, beyond the primordial one. The most relevant sources are the so-called secondary anisotropies, which arise after the last scattering epoch. These anisotropies can be divided into two categories: scattering secondaries, when the CMB photons scatter with electrons along the line of sight, and gravitational secondaries when effects are mediated by gravity [31]. Among the scattering secondaries we may list the thermal Sunyaev-Zeldovich effect, where hot electrons in clusters transfer energy to the CMB photons, the kinetic Sunyaev-Zeldovich effect produced by the bulk motion of the electrons in clusters, the Ostriker-Vishniac effect, produced by bulk motions modulated by linear density perturbations, and effects due to reionization processes. The scattering secondaries are most significant on small angular scales as density inhomogeneities, bulk and thermal motions grow and become sizeable on small length-scales when structure formation proceeds.
Gravitational secondaries arise from the change in energy of photons when the gravitational potential is time-dependent, the ISW effect, and gravitational lensing. At late times, when the Universe becomes dominated by the dark energy, the gravitational potential on linear scales starts to decay, causing the ISW effect mainly on large angular scales. Other secondaries that result from a time dependent potential are the Rees-Sciama effect, produced during the matter-dominated epoch by the time evolution of the potential on non-linear scales.
The fact that the potential never grows appreciably means that most second order effects created by gravitational secondaries are generically small compared to those created by scattering ones. However, when a photon propagates from the last scattering to us, its path may be deflected because of the gravitational lensing. This effect does not create anisotropies, but only modifies existing ones. Since photons with large wavenumbers k are lensed over many regions (∼ k/H, where H is the Hubble rate) along the line of sight, the corresponding second-order effect may be sizeable. The three-point function arising from the correlation of the gravitational lensing and ISW effects generated by the matter distribution along the line of sight [32,33] and the Sunyaev-Zeldovich effect [34] are large and detectable by Planck [35,36]. A crucial issue is the level of contamination to the extraction of the primordial NG the seconndary effects can produce. In Sec. 8 we briefly summarize some recent results about the level of CMB NG generated by some of these secondary effects.
Another relevant source of NG comes from the physics operating at the recombination. A naive esti- 1 A similar computation of the CMB anisotropies up to third-order from gravitational perturbations has been performed in Ref. [30], which is particularly relevant to provide a complete theoretical predition for cubic non-linearities characterizing the level of NG in the CMB through the connected four-point correlation function (trispectrum) [20]. mate would tell that these non-linearities are tiny being suppressed by an extra power of the gravitational potential. However, the dynamics at recombination is quite involved because all the non-linearities in the evolution of the baryon-photon fluid at recombination and the ones coming from general relativity should be accounted for. This complicated dynamics might lead to unexpected suppressions or enhancements of the NG at recombination. A step towards the evaluation of the three-point correlation function has been taken in Ref. [37] where some effects were taken into account in the so-called squeezed triangle limit, corresponding to the case when one wavenumber is much smaller than the other two and was outside the horizon at recombination. Refs. [38,39] (see also [40]), present the computation of the full system of Boltzmann equations, describing the evolution of the photon, baryon and Cold Dark Matter (CDM) fluids, at second order and neglecting polarization, These equations allow to follow the time evolution of the CMB anisotropies at second order on all angular scales from the early epochs, when the cosmological perturbations were generated, to the present time, through the recombination era. These calculations set the stage for the computation of the full second-order radiation transfer function at all scales and for a a generic set of initial conditions specifying the level of primordial non-Gaussianity. Of course, for specific effects on small angular scales like Sunyaev-Zel'dovich, gravitational lensing, etc., fully nonlinear calculations would provide a more accurate estimate of the resulting CMB anisotropy, however, as long as the leading contribution to second-order statistics like the bispectrum is concerned, second-order perturbation theory suffices.
The goal of this review is to summarize in a pedagogical form the present status of the evaluation of the impact of NG on the CMB anisotropies. This implies first of all determining how to set the initial conditions at second-order for the (gauge-invariant) CMB anisotropy when some source of primordial NG is present. The second step will be determining how primordial NG flows on small angular scales. In this review we will focus ourselves on the study of the second-order effects appearing at the recombination era when the CMB anisotropy is left imprinted. We will show how to derive the equations to evaluate CMB anisotropies, by computing the Boltzmann equations describing the evolution of the baryon-photon fluid up to second order. This permits to follow the time evolution of CMB anisotropies (up to second order) on all scales, from the early epoch, when the cosmological perturbations were generated, to the present time, through the recombination era. We will also provide the reader with some simplified analytical computation to evaluate the contamination of the recombination secondary effects onto the detection of primordial NG. The formalism for a more refined numerical analysis is also displayed and results for some worked examples will be also reported. The review is mainly based on a series of papers written by the authors along the past years on the subject (with various updates) and, as such, follows both a logic and a chronological order. It requires knowledge of the theory of cosmological perturbation at the linear level (which however we summarize in the Appendices). We have tried to write the different sections in a self-contained way. Nevertheless, we alert the reader that the level of complexity increases with the number of the sections.
The review is organized as follows. In Section 2 we provide a simple, but illuminating example to show why we do expect some NG present in the CMB anisotropy regardless if there is or not some primordial NG. In Section 3 we provide the reader with the necessary tools to study the dynamics at second-order in the gravity sector. In Section 4 we show how to set the initial conditions for the primordial NG, while in Section 5 we provide a gauge-invariant way to define the CMB temperature anisotropy at second-order on large scales. In Section 6 we go to small scales and present the full procedure to compute the Boltzmann equations necessary to follow the evolution of the non-linearities from the recombination epoch down to the present epoch. Section 7 presents some analytical solutions of the Boltzmann equations in the tight coupling limit, along the same lines of what is done at the linear level. The issue of contamination is addressed in Section 8, while in Section 9 we offer the reader with an analytical estimate of such a contamination. A more refined numerical work is presented in Section 10. Finally, in Section 11 some conclusions are given. This review has also some hopefully useful appendices: in Appendix A the reader can find the energy-momentm tensors at second-order, Appendix B gives the solutions of Eistein equations for the perturbed fluids up to second-order, while Appendix C offers the analytical solutions of the linearized Boltzmann equations in the tight coupling limit.
2 Why do we expect NG in the cosmological perturbations?
Before tackling the problem of interest -the computation of the cosmological perturbations at secondorder after the inflationary era-we first provide a simple, but insightful computation, derived from Ref. [28], which illustrates why we expect that the cosmological perturbations develop some NG even if the latter is not present at some primordial epoch. This example will help the reader to understand why the cosmological perturbations are inevitably affected by nonlinearities, beyond those arising at some primordial epoch. The reason is clear: gravity is nonlinear and it feeds this property into the cosmological perturbations during the post-inflationary evolution of the universe. As gravity talks to all fluids, this transmission is inevitable. To be specific, we focus on the CMB anisotropies. We will adopt the Poisson gauge which eliminates one scalar degree of freedom from the g 0i component of the metric and one scalar and two vector degrees of freedom from g ij . We will use a metric of the form where a(t) is the scale factor as a function of the cosmic time t, and ω i and χ ij are the vector and tensor peturbation modes respectively. Each metric perturbation can be expanded into a linear (first-order) and a second-order part, as for example, the gravitational potential Φ = Φ (1) + Φ (2) /2. However in the metric (1) the choice of the exponentials greatly helps in computing the relevant expressions, and thus we will always keep them where it is convenient. From Eq. (1) one recovers at linear order the well-known longitudinal gauge while at second order, one finds Φ (2) = φ (2) − 2(φ (1) ) 2 and Ψ (2) = ψ (2) + 2(ψ (1) ) 2 where φ (1) , ψ (1) and φ (2) , ψ (2) (with φ (1) = Φ (1) and ψ (1) = Ψ (1) ) are the first and second-order gravitational potentials in the longitudinal (Poisson) gauge adopted in Refs. [41,5] as far as scalar perturbations are concerned. We now consider the long wavelength modes of the CMB anisotropies, i.e. we focus on scales larger than the horizon at last-scattering. We can therefore neglect vector and tensor perturbation modes in the metric. For the vector perturbations the reason is that we are they contain gradient terms being produced as non-linear combination of scalar-modes and thus they will be more important on small scales (linear vector modes are not generated in standard mechanisms for cosmological perturbations, as inflation). The tensor contribution can be neglected for two reasons. First, the tensor perturbations produced from inflation on large scales give a negligible contribution to the higher-order statistics of the Sachs-Wolfe effect being of the order of (powers of) the slow-roll parameters during inflation (this holds for linear tensor modes as well as for tensor modes generated by the non-linear evolution of scalar perturbations during inflation).
Since we are interested in the cosmological perturbations on large scales, that is in perturbations whose wavelength is larger than the Hubble radius at last scattering, a local observer would see them in the form of a classical -possibly time-dependent -(nearly zero-momentum) homogeneous and isotropic background. Therefore, it should be possible to perform a change of coordinates in such a way as to absorb the super-Hubble modes and work with a metric of an homogeneous and isotropic Universe (plus, of course, cosmological perturbations on scale smaller than the horizon). We split the gravitational potential Φ as where Φ ℓ stands for the part of the gravitational potential receiving contributions only from the super-Hubble modes; Φ s receives contributions only from the sub-horizon modes where H is the Hubble rate computed with respect to the cosmic time, H =ȧ/a, and θ(x) is the step function. Analogous definitions hold for the other gravitational potential Ψ. By construction Φ ℓ and Ψ ℓ are a collection of Fourier modes whose wavelengths are larger than the horizon length and we may safely neglect their spatial gradients. Therefore Φ ℓ and Ψ ℓ are only functions of time. This amounts to saying that we can absorb the large-scale perturbations in the metric (1) by the following redefinitions a = a e −Ψ ℓ .
The new metric describes a homogeneous and isotropic Universe where for simplicity we have not included the sub-horizon modes. On super-horizon scales one can regard the Universe as a collection of regions of size of the Hubble radius evolving like unperturbed patches with metric (6) [43].
Let us now go back to the quantity we are interested in, namely the anisotropies of the CMB as measured today by an observer O. If she/he is interested in the CMB anisotropies at large scales, the effect of super-Hubble modes is encoded in the metric (6). During their travel from the last scattering surface -to be considered as the emitter point E -to the observer, the CMB photons suffer a redshift determined by the ratio of the emitted frequency ω E to the observed one ω O where T O and T E are the temperatures at the observer point and at the last scattering surface, respectively.
What is then the temperature anisotropy measured by the observer? The expression (7) shows that the measured large-scale anisotropies are made of two contributions: the intrinsic inhomogeneities in the temperature at the last scattering surface and the inhomogeneities in the scaling factor provided by the ratio of the frequencies of the photons at the departure and arrival points. Let us first consider the second contribution. As the frequency of the photon is the inverse of a time period, we get immediately the fully non-linear relation As for the temperature anisotropies coming from the intrinsic temperature fluctuation at the emission point, it maybe worth to recall how to obtain this quantity in the longitudinal gauge at first order. By expanding the photon energy density ρ γ ∝ T 4 γ , the intrinsic temperature anisotropies at last scattering are given by δ 1 T E /T E = (1/4)δ 1 ρ γ /ρ γ . One relates the photon energy density fluctuation to the gravitational perturbation first by implementing the adiabaticity condition δ 1 ρ γ /ρ γ = (4/3)δ 1 ρ m /ρ m , where δ 1 ρ m /ρ m is the relative fluctuation in the matter component, and then using the energy constraint of Einstein equations φ (1) = −(1/2)δ 1 ρ m /ρ m . The result is δ 1 T E /T E = −2Φ 1E /3. Summing this contribution to the anisotropies coming from the redshift factor (8) expanded at first order provides the standard (linear) Sachs-Wolfe effect δ 1 T O /T O = Φ 1E /3. Following the same steps, we may easily obtain its full non-linear generalization.
Let us first relate the photon energy density ρ γ to the energy density of the non-relativistic matter ρ m by using the adiabaticity conditon. Again here a bar indicates that we are considering quantities in the locally homogeneous Universe described by the metric (6). Using the energy continuity equation on large scales ∂ρ/∂t = −3H(ρ + P ), where H = d ln a/dt and P is the pressure of the fluid, one can easily show that there exists a conserved quantity in time at any order in perturbation theory [44] F ≡ ln a + 1 3 The perturbation δF is a gauge-invariant quantity representing the non-linear extension of the curvature perturbation ζ on uniform energy density hypersurfaces on superhorizon scales for adiabatic fluids [44]. Indeed, expanding it at first and second order one gets the corresponding definition ζ 1 = −ψ 1 − δ 1 ρ/ρ and the quantity ζ 2 introduced in Ref. [45]. We will come back to these definitions later. At first order the adiabaticity condition corresponds to set ζ 1γ = ζ 1m for the curvature perturbations relative to each component. At the non-linear level the adiabaticity condition generalizes to 1 3 or ln ρ m = ln ρ 3/4 γ .
Next we need to relate the photon energy density to the gravitational potentials at the non-linear level. The energy constraint inferred from the (0-0) component of Einstein equations in the matter-dominated era with the "barred" metric (6) is Using Eqs. (4) and (5) the Hubble parameter H reads where H = d ln a/dt is the Hubble parameter in the "unbarred" metric. Eq. (13) thus yields an expression for the energy density of the non-relativistic matter which is fully nonlinear, being expressed in terms of the gravitational potential Φ ℓ where we have droppedΨ ℓ which is negligible on large scales. By perturbing the expression (15) we are able to recover in a straightforward way the solutions of the (0-0) component of Einstein equations for a matter-dominated Universe in the large-scale limit obtained at second-order in perturbation theory. Indeed, recalling that Φ is perturbatively related to the quantity φ = φ (1) + φ (2) /2 used in Ref. [5] by The expression for the intrinsic temperature of the photons at the last scattering surface T E ∝ ρ 1/4 γ follows from Eqs. (11) and (15) Plugging Eqs. (8) and (17) into the expression (7) we are finally able to provide the expression for the CMB temperature which is fully nonlinear and takes into account both the gravitational redshift of the photons due to the metric perturbations at last scattering and the intrinsic temperature anisotropies From Eq. (18) we read the non-perturbative anisotropy corresponding to the Sachs-Wolfe effect Eq. (19) is one of the main results of this paper and represents at any order in perturbation theory the extension of the linear Sachs-Wolfe effect. At first order one gets and at second order 1 2 This result shows that the CMB anisotropies is nonlinear on large scales and that a source of NG is inevitably sourced by gravity.

Perturbing gravity
In this Section we provide the necessary tools to deal with perturbed gravity, giving the expressions for the Einstein tensor perturbed up to second-order around a flat Friedmann-Robertson-Walker background, and the relevant Einstein equations. In the following we will adopt the Poisson gauge which eliminates one scalar degree of freedom from the g 0i component of the metric and one scalar and two vector degrees of freedom from g ij . We rewrite the metric (1) as where a(η) is the scale factor as a function of the conformal time η. As we previously mentioned, for the vector and tensor perturbations, we will neglect linear vector modes since they are not produced in standard mechanisms for the generation of cosmological perturbations (as inflation), and we also neglect tensor modes at linear order, since they give a negligible contribution to second-order perturbations. Therefore we take ω i and χ ij to be second-order vector and tensor perturbations of the metric. Let us now give our definitions for the connection coefficients and their expressions for the metric (1). The number of spatial dimensions is n = 3. Greek indices (α, β, ..., µ, ν, ....) run from 0 to 3, while latin indices (a, b, ..., i, j, k, ....m, n...) run from 1 to 3. The total spacetime metric g µν has signature (−, +, +, +). The connection coefficients are defined as The Riemann tensor is defined as The Ricci tensor is a contraction of the Riemann tensor and in terms of the connection coefficient it is given by The Ricci scalar is given by contracting the Ricci tensor The Einstein tensor is defined as

The connection coefficients
For the connection coefficients we find

Einstein equations
The Einstein equations are written as G µν = κ 2 T µν , so that κ 2 = 8πG N , where G N is the usual Newtonian gravitational constant. They read Taking the traceless part of Eq. (32), we find where Q is defined by ∇ 2 Q = −P + 3N , with P ≡ P i i , and From Eq. (31), we may deduce an equation Here T µ ν is the energy momentum tensor accounting for different components, photons, baryons, dark matter. We will give the expressions later for each component in terms of the distribution functions.

Setting the initial conditions from the primordial NG
In this review we are concerned with the second-order evolution of the cosmological perturbations. This requires that we well define the initial conditions of the cosmological perturbations at second order. These boundary conditions may or may not contain already some level of NG. It they do, we say that there exist some primordial NG. The latter is usually defined in the epoch in which the comoving curvature perturbation remains constant on large superhorizon scales. In the standard single-field inflationary model, the first seeds of density fluctuations are generated on super-horizon scales from the fluctuations of a scalar field, the inflaton [1].

13
The crucial point is that the gauge-invariant curvature perturbation ζ remains constant on super-horizon scales after it has been generated during a primordial epoch and possible isocurvature perturbations are no longer present. Therefore, we may set the initial conditions at the time when ζ becomes constant. In particular, ζ (2) provides the necessary information about the "primordial" level of non-Gaussianity generated either during inflation, as in the standard scenario, or immediately after it, as in the curvaton scenario. Different scenarios are characterized by different values of ζ (2) . For example, in the standard single-field inflationary model ζ (2) = 2 ζ (1) 2 + O (ǫ, η) [6], where ǫ and η are the standard slow-roll parameters [1]. In general, we may parametrize the primordial non-Gaussianity level in terms of the conserved curvature perturbation as in Ref. [26] where the parameter a NL depends on the physics of a given scenario. For example in the standard scenario a NL ≃ 1, while in the curvaton case a NL = (3/4r) − r/2, where r ≈ (ρ σ /ρ) D is the relative curvaton contribution to the total energy density at curvaton decay [46,5]. Alternatives to the curvaton model are those models characterized by the curvature perturbation being generated by an inhomogeneity in the decay rate [9,10] or the mass [11] or of the particles responsible for the reheating after inflation. Other opportunities for generating the curvature perturbation occur at the end of inflation [13] and during preheating [14]. All these models generate a level of NG which is local as the NG part of the primordial curvature perturbation is a local function of the Gaussian part, being generated on superhorizon scales. In momentum space, the three point function, or bispectrum, arising from the local NG is dominated by the so-called "squeezed" configuration, where one of the momenta is much smaller than the other two.
Other models, such as DBI inflation [16] and ghost inflation [47], predict a different kind of primordial NG, called "equilateral", because the three-point function for this kind of NG is peaked on equilateral configurations, in which the lenghts of the three wavevectors forming a triangle in Fourier space are equal [18]. One of the best tools to detect or constrain the primordial large-scale non-Gaussianity is through the analysis of the CMB anisotropies, for example by studying the bispectrum [5]. In that case the standard procedure is to introduce the primordial non-linearity parameter f NL characterizing the primordial non-Gaussianity via the curvature perturbation [5] where the coefficient 3/5 arises from the first-order relation connecting the comoving curvature perturbation and the gravitational potential, ζ (1) = −5φ (1) /3, and the ⋆-product reminds us that one has to perform a convolution product in momentum space and that f NL is indeed momentum-dependent. To give the feeling of the resulting size of f NL when |a NL | ≫ 1, f NL ≃ 5a NL /3 (see Refs. [5,26]). Present limits on NG are summarized by −4 < f loc NL < 80 [48] and −125 < f equil NL < 435 [49] at 95% CL, where f loc NL and f equil NL stand for the non-linear parameter in the case in which the squeezed and the equilateral configurations dominate, respectively.

CMB anisotropies at second-order on large scales
In this section, we provide the exact expression for large-scale CMB temperature fluctuations at second order in perturbation theory. What this Section contains, therefore, should be considered as a more technical elaboration of what the reader can find in Section 2. The final expression we will find has various virtues. First, it is gauge-invariant. Second, from it one can unambiguously extract the exact definition of the nonlinearity parameter f NL which is used by the experimental collaborations to pin down the level of NG in the temperature fluctuations. Third, it contains a "primordial" term encoding all the information about the NG generated in primordial epochs, namely during or immediately after inflation, and depends upon the various fluctuation generation mechanisms. As such, the expression neatly disentangles the primordial contribution to the NG from that arising after inflation. Finally, the expression applies to all scenarios for the generation of cosmological perturbations.
In order to obtain our gauge-independent formula for the temperature anisotropies, we again perturb the spatially flat Robertson-Walker background. We expand the metric perturbations (22) in a first and a second-order part as Again, the functions φ (r) , ω ij , where (r) = (1, 2), stand for the rth-order perturbations of the metric. We can split ω The metric perturbations will transform according to an infinitesimal change of coordinates. From now on we limit ourselves to a second-order time shift where a prime denotes differentiation w.r.t. conformal time. In general a gauge corresponds to a choice of coordinates defining a slicing of spacetime into hypersurfaces (at fixed time η) and a threading into lines (corresponding to fixed spatial coordinates x), but in this Section only the former is relevant so that gauge-invariant can be taken to mean independent of the slicing. For example, under the time shift, the first-order spatial curvature perturbation ψ (1) transforms as ψ (1) (1) , and the traceless part of the spatial metric χ (1) ij turns out to be gauge-invariant. At second order in the perturbations we just give some useful examples like the transformation of the energy density and the curvature perturbation [41]. and ,i − Hα (2) .
We now construct in a gauge-invariant way temperature anisotropies at second order. Temperature anisotropies beyond the linear regime have been calculated in Refs. [24,41], following the photons path from last-scattering to the observer in terms of perturbed geodesics. The linear temperature anisotropies read [24,41] ∆T (1) where ij e i e j , the subscript E indicates that quantities are evaluated at last-scattering, e i is a spatial unit vector specifying the direction of observation and the integral is evaluated along the line-of-sight parametrized by the affine parameter λ. Eq. (45) includes the intrinsic fractional temperature fluctuation at emission τ E , the Doppler effect due to emitter's velocity v (1)i E and the gravitational redshift of photons, including the Integrated Sachs-Wolfe (ISW) effect. We omitted monopoles due to the observer O (e.g. the gravitational potential ψ (1) O evaluated at the event of observation), which, being independent of the angular coordinate, can be always recast into the definition of temperature anisotropies. Notice however that the physical meaning of each contribution in Eq. (45) is not gauge-invariant, as the different terms are gauge-dependent. However, it is easy to show that the whole expression (45) is gauge-invariant. Since the temperature T is a scalar, the intrinsic temperature fluctuation transforms as τ (1) , having used the fact that the temperature scales as T ∝ a −1 . Notice, instead, that the velocity v (1)i E does not change. Therefore, using the transformations of metric perturbations we find where we have used the fact that the integral is evaluated along the line-of-sight which can be parametrized by the background geodesics x (0)µ = λ, (λ O − λ E )e i (with dλ/dη = 1), and the decomposition for the total derivative along the path for a generic function f (λ, Eq. (46) shows that the expression (45) for first-order temperature anisotropies is indeed gauge-invariant (up to monopole terms related to the observer O). Temperature anisotropies can be easily written in terms of particular combinations of perturbations which are manifestly gauge-invariant. For the gravitational potentials we consider the gauge-invariant definitions ψ (1) GI = ψ (1) − Hω (1) and φ (1) GI = φ (1) + Hω (1) + ω (1) ′ . For the (0 − i) component of the metric and the traceless part of the spatial metric we define ω ij . For the matter variables we use a gauge-invariant intrinsic temperature fluctuation τ (1) Following the same steps leading to Eq. (46) one gets the linear temperature anisotropies in Eq. (45) in terms of these gauge-invariant quantities where A (1) e i e j and we omitted the subscript E. For the primordial fluctuations we are interested in the large-scale modes set by the curvature perturbation ζ (1) . Defining a gauge-invariant density perturbation δ (1) ρ GI = δ (1) ρ + ρ ′ ω (1) , we write the curvature perturbation as ζ (1) ρ GI /ρ ′ ). Since for adiabatic perturbations in the radiation (γ) and matter (m) eras (1/4)(δ (1) ρ γ /ρ γ ) = (1/3)(δ (1) ρ m /ρ m ), one can write the intrinsic temperature fluctuation as (1) ρ GI /ρ ′ ). In the large-scale limit, from Einstein equations, in the matter era φ (1) GI . Thus we obtain the large-scale limit of temperature anisotropies (47) GI /3, i.e. the usual Sachs-Wolfe effect.
At second order, the procedure is similar to the one described so long, though more lengthy and cumbersome. We only provide the reader with the main steps to get the final expression. The secondorder temperature fluctuations in terms of metric perturbations read [24,41] case, by choosing the shifts α (r) such that ω (r) = 0. For example, we consider the gauge-invariant gravitational potential [25] φ (2) where (1) ]. Expressing the second-order temperature anisotropies (48) in terms of our gauge-invariant quantities and taking the large-scale limit we find (having dropped the subscript E), and the gauge-invariant intrinsic temperature fluctuation at emission is We have dropped those terms which represent integrated contributions and other second-order small-scale effects that can be distinguished from the large-scale part through their peculiar scale dependence. At this point we make use of Einstein's equations. We take the expression for ζ (2) in Eqs. (37) and (38), and we use the (0 − 0) component and the traceless part of the (i − j) Einstein's equations (30) and (32) after having appropriately expanded the exponentials. Thus, on large scales we find that the temperature anisotropies are given by where we have defined a kernel Eq. (55) is the main result of this Section. It clearly shows that there are two contributions to the final nonlinearity in the large-scale temperature anisotropies. The contribution, [ζ GI ) 2 ], comes from the "primordial" conditions set during or after inflation. They are encoded in the curvature perturbation ζ which remains constant once it has been generated. The remaining part of Eq. (55) describes the postinflation processing of the primordial non-Gaussian signal due to the nonlinear gravitational dynamics, including also second-order corrections at last scattering to the Sachs-Wolfe effect [24,41]. Thus, the expression in Eq. (55) allows to neatly disentangle the primordial contribution to NG from that coming from that arising after inflation and from it we can extract the nonlinearity parameter f NL , see the expression (40) which is usually adopted to phenomenologically parametrize the NG level of cosmological perturbations and has become the standard quantity to be observationally The definition of f NL adopted in the analyses performed in Refs. [35] goes through the conventional Sachs-Wolfe formula ∆T /T = −Φ/3 where Φ is Bardeen's potential, which is conventionally expanded as (up to a constant offset, which only affects the temperature monopole) Φ = Φ L + f NL (Φ L ) ⋆ (Φ L ). Therefore, using ζ (1) GI during matter domination, from Eq. (55) we read the nonlinearity parameter in momentum space where with k 3 + k 1 + k 2 = 0 and k = |k 3 |. To obtain Eq. (57) we have made use of the expression (39) to set the initial conditions. In particular within the standard scenario where cosmological perturbations are due to the inflaton the primordial contribution to NG is given by a NL = 1 − 1 4 (n ζ − 1) [6], where the spectral index is expressed in terms of the usual slow-roll parameters as n ζ − 1 = −6ǫ + 2η [1]. The nonlinearity parameter from inflation now reads Therefore the main NG contribution comes from the post-inflation evolution of the second-order perturbations which give rise to order-one coefficients, while the primordial contribution is proportional to |n ζ − 1| ≪ 1. This is true even in the "squeezed" limit first discussed by Maldacena in [6], where one of the wavenumbers is much smaller than the other two, e.g. k 1 ≪ k 2,3 and K → 0.

CMB anisotropies at second-order at all scales
As we already mentioned in the Introduction, despite the importance of NG in CMB anisotropies, little effort has ben made so far to provide accurate theoretical predictions of it. On the contrary, the vast majority of the literature has been devoted to the computation of the bispectrum of either the comovig curvature perturbation or the gravitational potential on large scales within given inflationary models. These, however, are not the physical quantities which are observed. One should instead provide a full prediction for the second-order radiation transfer function. A preliminary step towards this goal has been taken in Refs. [23,24,28] (see also [29]) where the full second-order radiation transfer function for the CMB anisotropies on large angular scales in a flat universe filled with matter and cosmological constant was computed, including the second-order generalization of the Sachs-Wolfe effect, both the early and late Integrated Sachs-Wolfe (ISW) effects and the contribution of the second-order tensor modes. 3 We have partly reported about these works in the previous Sections.
In this Section we wish to offer a summary of some of the second-order effects in the CMB anisotropies on small scales. There are many sources of NG in CMB anisotropies, beyond the primordial one. The most relevant sources are the so-called secondary anisotropies, which arise after the last scattering epoch. These anisotropies can be divided into two categories: scattering secondaries, when the CMB photons scatter with electrons along the line of sight, and gravitational secondaries when effects are mediated by gravity [31]. Among the scattering secondaries we may list the thermal Sunyaev-Zeldovich effect, where hot electrons in clusters transfer energy to the CMB photons, the kinetic Sunyaev-Zeldovich effect produced by the bulk motion of the electrons in clusters, the Ostriker-Vishniac effect, produced by bulk motions modulated by linear density perturbations, and effects due to reionization processes. The scattering secondaries are most significant on small angular scales as density inhomogeneities, bulk and thermal motions grow and become sizeable on small length-scales when structure formation proceeds.
Gravitational secondaries arise from the change in energy of photons when the gravitational potential is time-dependent, the ISW effect, and gravitational lensing. At late times, when the Universe becomes dominated by the dark energy, the gravitational potential on linear scales starts to decay, causing the ISW effect mainly on large angular scales. Other secondaries that result from a time dependent potential are the Rees-Sciama effect, produced during the matter-dominated epoch by the time evolution of the potential on non-linear scales.
The fact that the potential never grows appreciably means that most second order effects created by gravitational secondaries are generically small compared to those created by scattering ones. However, when a photon propagates from the last scattering to us, its path may be deflected because of the gravitational lensing. This effect does not create anisotropies, but only modifies existing ones. Since photons with large wavenumbers k are lensed over many regions (∼ k/H, where H is the Hubble rate) along the line of sight, the corresponding second-order effect may be sizeable. The three-point function arising from the correlation of the gravitational lensing and ISW effects generated by the matter distribution along the line of sight [32,33] and the Sunyaev-Zeldovich effect [34] are large and detectable by Planck [35].
Another relevant source of NG comes from the physics operating at the recombination. A naive estimate would tell that these non-linearities are tiny being suppressed by an extra power of the gravitational potential. However, the dynamics at recombination is quite involved because all the non-linearities in the evolution of the baryon-photon fluid at recombination and the ones coming from general relativity should be accounted for. This complicated dynamics might lead to unexpected suppressions or enhancements of the NG at recombination. A step towards the evaluation of the three-point correlation function has been taken in Ref. [37] where some effects were taken into account in the in so-called squeezed triangle limit, corresponding to the case when one wavenumber is much smaller than the other two and was outside the horizon at recombination. This Section, which is based on Refs. [38,39], present the computation of the full system of Boltzmann equations, describing the evolution of the photon, baryon and Cold Dark Matter (CDM) fluids, at second order and neglecting polarization, These equations allow to follow the time evolution of the CMB anisotropies at second order on all angular scales from the early epochs, when the cosmological perturbations were generated, to the present time, through the recombination era. These calculations set the stage for the computation of the full second-order radiation transfer function at all scales and for a a generic set of initial conditions specifying the level of primordial non-Gaussianity. Of course on small angular scales, fully non-linear calculations of specific effects like Sunyaev-Zel'dovich, gravitational lensing, etc. would provide a more accurate estimate of the resulting CMB anisotropy, however, as long as the leading contribution to second-order statistics like the bispectrum is concerned, second-order perturbation theory suffices.

The collisionless Boltzmann equation for photons
We are interested in the anisotropies in the cosmic distribution of photons and inhomogeneities in the matter. Photons are affected by gravity and by Compton scattering with free electrons. The latter are tightly coupled to protons. Both are, of course, affected by gravity. The metric which determines the gravitational forces is influenced by all these components plus CDM (and neutrinos). Our plan is to write down Boltzmann equations for the phase-space distributions of each species in the Universe.
The phase-space distribution of particles g(x i , P µ , η) is a function of spatial coordinates x i , conformal time η, and momentum of the particle P µ = dx µ /dλ where λ parametrizes the particle path. Through the constraint P 2 ≡ g µν P µ P ν = −m 2 , where m is the mass of the particle one can eliminate P 0 and g(x i , P j , η) gives the number of particles in the differential phase-space volume dx 1 dx 2 dx 3 dP 1 dP 2 dP 3 . In the following we will denote the distribution function for photons with f . The photons' distribution evolves according to the Boltzmann equation where the collision term is due to the scattering of photons off free electrons. In the following we will derive the left-hand side of Eq. (60) while in the next section we will compute the collision term. For photons we can impose P 2 ≡ g µν P µ P ν = 0 and using the metric (1) in the conformal time η we find where we define From the constraint (61) Notice that we immediately recover the usual zero and first-order relations P 0 = p/a and P 0 = p(1 − Φ (1) )/a. The components P i are proportional to pn i , where n i is a unit vector with n i n i = δ ij n i n j = 1. We can write P i = Cn i , where C is determined by so that P i = p a n i e −2Ψ + χ km n k n m −1/2 = p a n i e Ψ 1 − 1 2 where the last equality holds up to second order in the perturbations. Again we recover the zero and first-order relations P i = pn i /a and P i = pn i (1 + Ψ (1) )/a respectively. Thus up to second order we can write Eq. (65) and (66) allow us to replace P 0 and P i in terms of the variables p and n i . Therefore, as it is standard in the literature, from now on we will consider the phase-space distribution f as a function of the momentum p = pn i with magnitude p and angular direction n i , f ≡ f (x i , p, n i , η). Thus, in terms of these variables, the total time derivative of the distribution function reads

21
In the following we will compute dx i /dη, dp/dη and dn i /dη. a) dx i /dη: From and from Eq. (65) and (66) b) dp/dη: For dp/dη we make use of the time component of the geodesic equation Using the metric (1) we find On the other hand the expression (66) of P 0 in terms of p and n i gives Thus Eq. (70) allows us express dp/dη as where in Eq. (71) we have replaced P 0 and P i by Eqs. (66) and (65). Notice that in order to obtain Eq. (73) we have used the following expressions for the total time derivatives of the metric perturbations and where we have taken into account that ω i is already a second-order perturbation so that we can neglect dn i /dη which is at least a first order quantity, and we can take the zero-order expression in Eq. (69), dx i /dη = n i . In fact there is also an alternative expression for dp/dη which turns out to be useful later and which can be obtained by applying once more Eq. (74) c) dn i /dη: We can proceed in a similar way to compute dn i /dη. Notice that since in Eq. (67) it multiplies ∂f /∂n i which is first order, we need only the first order perturbation of dn i /dη. We use the spatial components of the geodesic equations dP i /dλ = −Γ i αβ P α P β written as For the right-hand side we find, up to second order, while the expression (65) of P i in terms of our variables p and n i in the left-hand side of Eq. (77) brings Thus, using the expression (65) for P i and (63) for P 0 in Eq. (78), together with the previous result (73), the geodesic equation (77) gives the following expression dn i /dη (valid up to first order) To proceed further we now expand the distribution function for photons around the zero-order value f (0) which is that of a Bose-Einstein distribution where T (η) is the average (zero-order) temperature and the factor 2 comes from the spin degrees of photons. The perturbed distribution of photons will depend also on x i and on the propagation direction n i so as to account for inhomogeneities and anisotropies where we split the perturbation of the distribution function into a first and a second-order part. The Boltzmann equation up to second order can be written in a straightforward way by recalling that the total time derivative of a given i-th perturbation, as e.g. df (i) /dη is at least a quantity of the i-th order.
Thus it is easy to realize, looking at Eq. (67), that the left-hand side of Boltzmann equation can be written up to second order as where for simplicity in Eq. (83) we have already used the background Boltzmann equation (df /dη)| (0) = 0. In Eq. (83) there are all the terms which will give rise to the integrated Sachs-Wolfe effects (corresponding to the terms which explicitly depend on the gravitational perturbations), while other effects, such as the gravitational lensing, are still hidden in the (second-order part) of the first term. In fact in order to obtain Eq. (83) we just need for the time being to know the expression for dp/dη, Eq. (76).

Collision term 6.2.1 The Collision Integral
In this section we focus on the collision term due to Compton scattering where we have indicated the momentum of the photons and electrons involved in the collisions. The collision term will be important for small scale anisotropies and spectral distortions. The important point to compute the collision term is that for the epoch of interest very little energy is transferred. Therefore one can proceed by expanding the right hand side of Eq. (60) both in the small perturbation, Eq. (82), and in the small energy transfer.
The collision term is given (up to second order) by where a is the scale factor and 4 where E(q) = (q 2 +m 2 e ) 1/2 , M is the amplitude of the scattering process, ensures the energy-momentum conservation and g is the distribution function for electrons. The Pauli suppression factors (1 − g) have been dropped since for the epoch of interest the density of electrons n e is low. The electrons are kept in thermal equilibrium by Coulomb interactions with protons and they are non-relativistic, thus we can take a Maxwell-Boltzmann distribution around some bulk velocity v By using the three dimensional delta function the energy transfer is given by E(q) − E(q + p − p ′ ) and it turns out to be small compared to the typical thermal energies In Eq. (88) we have used E(q) = m e + q 2 /2m e and the fact that, since the scattering is almost elastic (p ≃ p ′ ), (p − p ′ ) is of order p ∼ T , with q much bigger than (p − p ′ ). In general, the electron momentum has two contributions, the bulk velocity (q = m e v) and the thermal motion (q ∼ (m e T ) 1/2 ) and thus the parameter expansion q/m e includes the small bulk velocity v and the ratio (T /m e ) 1/2 which is small because the electrons are non-relativistic. The expansion of all the quantities entering the collision term in the energy transfer parameter and the integration over the momenta q and q ′ is described in details in Ref. [53]. It is easy to realize that we just need the scattering amplitude up to first order since at zero order g(q ′ ) = g(q + p − p ′ ) = g(q) and δ(E(q) + p − E(q ′ ) − p ′ ) = δ(p − p ′ ) so that all the zero-order quantities multiplying |M | 2 vanish. To first order where cosθ = n · n ′ is the scattering angle and σ T the Thompson cross-section. The resulting collision term up to second order is given by [53] where we arrange the different contributions following Ref. [53]. The first order term reads while the second-order terms have been separated into four parts. There is the so-called anisotropy suppression term a term which depends on the second-order velocity perturbation defined by the expansion of the bulk flow as a set of terms coupling the photon perturbation to the velocity and a set of source terms quadratic in the velocity The last contribution are the Kompaneets terms describing spectral distortions to the CMB Let us make a couple of comments about the various contributions to the collision term. First, notice the term c (2) v (p, p ′ ) due to second-order perturbations in the velocity of electrons which is absent in Ref. [53]. In standard cosmological scenarios (like inflation) vector perturbations are not generated at linear order, so that linear velocities are irrotational v (1)i = ∂ i v (1) . However at second order vector perturbations are generated after horizon crossing as non-linear combinations of primordial scalar modes. Thus we must take into account also a transverse (divergence-free) component, v (2) (2)i T = 0. As we will see such vector perturbations will break azimuthal symmetry of the collision term with respect to a given mode k, which instead usually holds at linear order. Secondly, notice that the number density of electrons appearing in Eq. (90) must be expanded as n e =n e (1 + δ e ) and then gives rise to second-order contributions in addition to the list above, where we split δ e = δ (1) e /2 into a first-and second-order part. In particular the combination with the term proportional to v in c (1) (p, p ′ ) gives rise to the so-called Vishniac effect, as discussed in Ref. [53].

Computation of different contributions to the collision term
In the integral (90) over the momentum p ′ the first-order term gives the usual collision term where one uses the decomposition in Legendre polynomials where ϑ is the polar angle of n, cosϑ = n ·v.
In the following we compute the second-order collision term separately for the different contributions, using the notation C(p) = C (1) (p) + C (2) (p)/2. We have not reported the details of the calculation of the first-order term because for its second-order analog, c ∆ (p, p ′ ) + c (2) v (p, p ′ ), the procedure is the same. The important difference is that the second-order velocity term includes a vector part, and this leads to a generic angular decomposition of the distribution function (for simplicity drop the time dependence) such that Such a decomposition holds also in Fourier space. The notation at this stage is a bit confusing, so let us restate it: superscripts denote the order of the perturbation; the subscripts refer to the moments of the distribution. Indeed at first order one can drop the dependence on m setting m = 0 using the fact that the distribution function does not depend on the azimuthal angle φ. In this case the relation with f To perform the angular integral we write the angular dependence on the scattering angle cosθ = n · n ′ in terms of the Legendre polynomials where in the last step we used the addition theorem for spherical harmonics Using the decomposition (100) and the orthonormality of the spherical harmonics we find It is easy to recover the result for the corresponding first-order contribution in Eq. (226) by using Eq. (202). b) c (2) v (p, p ′ ): Let us fix for simplicity our coordinates such that the polar angle of n ′ is defined by µ ′ =v (2) · n ′ with φ ′ the corresponding azimuthal angle. The contribution of c (2) v (p, p ′ ) to the collision term is then We can use Eq. (103) which in our coordinate system reads so that dφ ′ 2π P 2 (n · n ′ ) = P 2 (n ·v (2) )P 2 (n ′ ·v (2) ) = P 2 (µ)P 2 (µ ′ ) .
By using the orthonormality of the Legendre polynomials and integrating by parts over p ′ we find As it is clear by the presence of the scalar product v (2) ·p the final result is independent of the coordinates chosen. c) c (2) ∆v (p, p ′ ): Let us consider the contribution from the first term where the velocity has to be considered at first order. In the integral (90) it brings The procedure to do the integral is the same as above. We use the same relations as in Eqs. (107) and (108) where now the angles are those taken with respect to the first-order velocity. This eliminates the integral over φ ′ , and integrating by parts over p ′ yields We now use the decomposition (98) and the orthonormality of the Legendre polynomials to find where we have used

29
In whose contribution to the collision term is This integration proceeds through the same steps as for C (2) ∆v(I) (p). In particular by noting that cosθ(1 − cosθ) = −1/3 + P 1 (cosθ) − 2P 3 (cosθ)/3, Eqs. (107) and (108) allows to compute and using the decomposition (98) we arrive at We then obtain As far as the remaining terms, these have already been computed in Ref. [53] (see also Ref. [54]) and here we just report them The term proportional to the velocity squared yield a contribution to the collision term The terms responsible for the spectral distortions give Finally, we write also the part of the collision term coming from Eq. (96) (120)

Final expression for the collision term
Summing all the terms we find the final expression for the collision term (90) up to second order and

31
Notice that there is an internal hierarchy, with terms which do not depend on the baryon velocity v, terms proportional to v · n and then to (v · n) 2 , v and v 2 (apart from the Kompaneets terms). In particular notice the term proportional to δ (1) e v · n is the one corresponding to the Vishniac effect. We point out that we have kept all the terms up to second order in the collision term. In Refs. [53,54] many terms coming from c (2) ∆v have been dropped mainly because these terms are proportional to the photon distribution function f (1) which on very small scales (those of interest for reionization) is suppressed by the diffusion damping. Here we want to be completely general and we have to keep them.
At first order it is useful to characterize the perturbations to the Bose-Einstein distribution function (81) in terms of a perturbation to the temperature as Thus it turns out that where we have used the fact that ∂f /∂Θ| Θ=0 = −p∂f (0) /∂p. In terms of this variable Θ (1) the linear collision term (122) will now become proportional to −p∂f (0) /∂p which contains the only explicit dependence on p, and the same happens for the left-hand side, Eq. (124). This is telling us that at first order Θ (1) does not depend on p but only on . This is well known and the physical reason is that at linear order there is no energy transfer in Compton collisions between photons and electrons. Therefore, the Boltzmann equation for Θ (1) reads where we made us of f ℓ , according to the decomposition of Eq. (98), and we have taken the zero-order expressions for dx i /dη, dropping the contribution from dn i /dη in Eq. (67) since it is already first-order.
Notice that, since Θ (1) is independent of p, it is equivalent to consider the quantity being ∆ (1) = 4Θ (1) at this order. The physical meaning of ∆ (1) is that of a fractional energy perturbation (in a given direction). From Eq. (83) another way to write an equation for ∆ (1) -the so-called brightness equation -is

Second order
The previous results show that at linear order the photon distribution function has a Planck spectrum with the temperature that at any point depends on the photon direction. At second order one could characterize the perturbed photon distribution function in a similar way as in Eq. (125) where by expanding Θ = Θ (1) + Θ (2) /2 + ... as usual one recovers the first-order expression. For example, in terms of Θ, the perturbation of f (1) is given by Eq. (126), while at second order However, as discussed in details in Refs. [54,53], now the second-order perturbation Θ (2) will not be momentum independent because the collision term in the equation for Θ (2) does depend explicitly on p (defining the combination −(p∂f (0) /∂p) −1 f (2) does not lead to a second-order momentum independent equation as above). Such dependence is evident, for example, in the terms of C (2) (p), Eq. (123), proportional to v or v 2 , and in the Kompaneets terms. The physical reason is that at the non-linear level photons and electrons do exchange energy during Compton collisions. As a consequence spectral distorsions are generated. For example, in the isotropic limit, only the Kompaneets terms survive giving rise to the Sunyaev-Zeldovich distorsions. As discussed in Ref. [54], the Sunyaev-Zeldovich distorsions can also be obtained with the correct coefficients by replacing the average over the direction electron v 2 with the mean squared thermal velocity v 2 th = 3T e /m e in Eq. (123). This is due simply to the fact that the distinction between thermal and bulk velocity of the electrons is just for convenience. This fact also shows that spectral distorsions due to the bulk flow (kinetic Sunyaev-Zeldovich) has the same form as the thermal effect. Thus spectral distorsions can be in general described by a global Compton y-parameter (see Ref. [54] for a full discussion of spectral distorsions). However in the following we will not be interested in the frequency dependence but only in the anisotropies of the radiation distribution. Therefore we can integrate over the momentum p and define [54,53] as in Eq. (128).
Integration over p of Eqs. (83)-(123) is straightforward using the following relations Here N = dpp 3 f (0) is the normalization factor (it is just proportional the background energy density of photonsρ γ ). At first order one recovers Eq. (129). At second order we find where we have expanded the angular dependence of ∆ as in Eq. (99) with where we recall that the superscript stands by the order of the perturbation. When going to Fourier space some convolutions will appear and the coefficients ∆ It is understood that on the left-hand side of Eq. (134) one has to pick up for the total time derivatives only those terms which contribute to second order. Thus we have to take 1 2 ,j + Ψ where we used Eqs. (69) and (80).

Hierarchy equations for multipole moments
Let us now move to Fourier space. In the following, for a given wave-vector k we will choose the coordinate system such that e 3 =k and the polar angle of the photon momentum is ϑ, with µ = cosϑ =k · n. Then Eq. (134) can be written as where S(k, n, η) can be easily read off Eq. (134). We now expand the temperature anisotropy in the multipole moments ∆ (2) ℓm in order to obtain a system of coupled differential equations. By applying the angular integral of Eq. (136) to Eq. (139) we find where the first term on the right-hand side of Eq. (140) has been obtained by using the relation is the angular mode for the decomposition (135) and κ ℓm = √ l 2 − m 2 . This relation has been discussed in Refs. [56,57] and corresponds to the term n i ∂∆ (2) /∂x i in Eq. (134). In Eq. (140) S lm are the multipole moments of the source term according to the decomposition (136). We do not show its complete expression here, since it is very long, and we prefer to make some general comments (for some specific examples about the terms S lm see Sec. 10). As it should be already clear from the calculations involving second-order perturbations performed so far, Eq. (140) has the same functional form of the Boltzmann equation at linear order (just replace the order of the linear perturbations (1) with that of intrinsically second-order terms (2)) with the exception of the source terms S lm which now contain both intrinsic second-order perturbations and also products of first-order perturbations. Therefore as expected, at second order we recover some intrinsic effects which are characteristic of the linear regime. In Eq. (140) the relation (141) represents the free streaming effect: when the radiation undergoes free-streaming, the inhomogeneities of the photon distribution are seen by the obsever as angular anisotropies. At first order it is responsible for the hierarchy of Boltzmann equations coupling the different ℓ modes, and it represents a projection effect of fluctuations on a scale k onto the angular scale ℓ ∼ kη. The term τ ′ ∆ (2) ℓm causes an exponential suppression of anisotropies in the absence of the source term S ℓm . The source term contains additional scattering processes and gravitational effects. The intrinsically second-order part of the source term just reproduces the expression of the first order case. Of course the dynamics of the second-order metric and baryon-velocity perturbations which appear will be different and governed by the second-order Einstein equations and continuity equations. The remaining terms in the source are second-order effects generated as non-linear combinations of the primordial (first-order) perturbations. Notice in particular that they involve the first-order anisotropies ∆ (1) ℓ and as a consequence such terms contribute to generate the hierarchy of equations (apart from the free-streaming effect). On large scales (above the horizon at recombination) we can say that the main effects are due to gravity, and they include the Sachs-Wolfe and the (late and early) Sachs-Wolfe effect due to the redshift photons suffer when travelling through the second-order gravitational potentials. These, toghether with the contribution due to the second-order tensor modes, have been already studied in details in Ref. [23] (see also [29]). Another important gravitational effect is that of lensing of photons as they travel from the last scattering surface to us. A contribution of this type is given by the last term of Eq. (138). On the other hand examples of second-order scattering effects are the terms proportional to the square of the baryon-velocity fluid (v 2 ), giving rise to the quadratic Doppler effect (like those in the last line of Eq. (134)), or those coupling the photon fluctuations to the baryon velocity (second-line from the bottom of Eq. (134)). The Vishniac effect corresponds to the terms proportional to δ (1) e . Finally notice that in the Boltzmann equation (134)) the second-order baryon velocity appears. At linear order the baryon velocity is irrotational, meaning that it is the gradient of a potential, and thus in Fourier space it is parallel tok, and following the same notation of Ref. [55], we write The second-order velocity perturbation will contain a transverse (divergence-free) part whose components are orthogonal tok = e 3 , and we can write where e i form an orhtonormal basis withk. The second-order perturbation ω i is decomposed in a similar way, with ω ±1 the corresponding components (in this case in the Poisson gauge there is no scalar component). Similarly for the tensor perturbation χ ij we have indicated its amplitudes as χ ±2 in the decomposition [57] In computing the source term S lm one has to take into account that in the gravitational part of the Boltzmann equation and in the collsion term there are some terms, like δ (1) e v, which still can be decomposed in the scalar and transverse parts in Fourier space as in Eq. (143). For a generic quantity f (x)v one can indicate the corresponding scalar and vortical components with (f v) m and their explicit expression is easily found by projecting the Fourier modes of f (x)v along thek = e 3 and (e 2 ∓ ie 1 ) directions Similarly for a term like f (x)∇g(x) one can use the notation

Integral solution of the second-order Boltzmann equation
As in linear theory, one can derive an integral solution of the Boltzmann equation (134) in terms of the source term S. Following the standard procedure (see e.g. Ref. [5]) for linear perturbations, we write the left-hand side as where η 0 stands by the present time. The expression of the photon moments ∆ (2) ℓm can be obtained as usual from Eq. (136). In the previous section we have already found the coefficients for the decomposition of source term S In Eq. (221) there is an additional angular dependence in the exponential. It is easy to take it into account by recalling that Thus the angular integral (136) is computed by using the decomposition of the source term (148) and Eq. (149) where the Wigner 3 − j symbols appear because of the Gaunt integrals Since the second of the Wigner 3-j symbols in Eq. (150) is nonzero only if m = m 2 , our solution can be rewritten to recover the corresponding expression found for linear anisotropies in Refs.
where j are the so called radial functions. Of course the main information at second order is included in the source term containing different effects due to the non-linearity of the perturbations.
In the total angular momentum method of Ref. [57] Eq. (152) is interpreted just as the intergration over the radial coordinate (χ = η 0 − η) of the projected source term. Another important comment is that, as in linear theory, the integral solution (150) is in fact just a formal solution, since the source term S contains itself the second-order photon moments up to l = 2. This means that one has anyway to resort to the hierarchy equations for photons, Eq. (140), to solve for these moments. Nevertheless, as in linear theory, one expects to need just a few moments beyond ℓ = 2 in the hierarchy equations, and once the moments entering in the source function are computed the higher moments are obtained from the integral solution. Thus the integral solution should in fact be more advantageous than solving the system of coupled equations (140).

The Boltzmann equation for baryons and cold dark matter
In this section we will derive the Boltzmann equation for massive particles, which is the case of interest for baryons and dark matter. These equations are necessary to find the time evolution of number densities and velocities of the baryon fluid which appear in the brightness equation, thus allowing to close the system of equations. Let us start from the baryon component. Electrons are tightly coupled to protons via Coulomb interactions. This forces the relative energy density contrasts and the velocities to a common value, δ e = δ p ≡ δ b and v e = v p ≡ v, so that we can identify electrons and protons collectively as "baryonic" matter.
To derive the Boltzmann equation for baryons let us first focus on the collisionless equation and compute therefore dg/dη, where g is the distribution function for a massive species with mass m. One of the differences with respect to photons is just that baryons are non-relativistic for the epochs of interest. Thus the first step is to generalize the formulae in Section 4 up to Eq. (80) to the case of a massive particle. In this case one enforces the constraint Q 2 = g µν Q µ Q ν = −m 2 and it also useful to use the particle energy E = q 2 + m 2 , where q is defined as in Eq. (62). Moreover in this case it is very convenient to take the distribution function as a function of the variables q i = qn i , the position x i and time η, without using the explicit splitting into the magnitude of the momentum q (or the energy E) and its direction n i . Thus the total time derivative of the distribution functions reads We will not give the details of the calculation since we just need to replicate the same computation we did for the photons. For the four-momentum of the particle notice that Q i has the same form as Eq. (65), while for Q 0 we find In the following we give the expressions for dx i /dη and dq i /dη. a) As in Eq. (69) dx i /dη = Q i /Q 0 and it turns out to be b) For dq i /dη we need the expression of Q i which is the same as that of Eq. (65) The spatial component of the geodesic equation, up to second order, reads Proceeding as in the massless case we now take the total time derivative of Eq. (156) and using Eq. (157) we find We can now write the total time derivative of the distribution function as This equation is completely general since we have just solved for the kinematics of massive particles. As far as the collision terms are concerned, for the system of electrons and protons we consider the Coulomb scattering processes between the electrons and protons and the Compton scatterings between photons and electrons where we have adopted the same formalism of Ref. [42] with p and p ′ the initial and final momenta of the photons, q and q ′ the corresponding quantities for the electrons and for protons Q and Q ′ . The integral over different momenta is indicated by and thus one can read c eγ as the unintegrated part of Eq. (85), and similarly for c ep (with the appropriate amplitude |M | 2 ). In Eq. (160) Compton scatterings between protons and photons can be safely neglected because the amplitude of this process has a much smaller amplitude than Compton scatterings with electrons being weighted by the inverse squared mass of the particles. At this point for the photons we considered the perturbations around the zero-order Bose-Einstein distribution function (which are the unknown quantities). For the electrons (and protons) we can take the thermal distribution described by Eq. (87). Moreover we will take the moments of Eqs. (160)-(161) in order to find the energy-momentum continuity equations.

Energy continuity equations
We now integrate Eq. (159) over d 3 q/(2π) 3 . Let us recall that in terms of the distribution function the number density n e and the bulk velocity v are given by and where one can set E ≃ m e since we are considering non-relativistic particles. We will also make use of the following relations when integrating over the solid angle dΩ Finally notice that dE/dq = q/E and ∂g/∂q = (q/E)∂g/∂E. Thus the first two integrals just brings n ′ e and (n e v i ) ,i . Notice that all the terms proportional to the second-order vector and tensor perturbations of the metric give a vanishing contribution at second order since in this case we can take the zero-order distribution functions which depends only on η and E, integrate over the direction and use the fact that δ ij χ ij = 0. The trick to solve the remaining integrals is an integration by parts over q i . We have an integral like (the one multiplying (Ψ ′ − H)) after an integration by parts over q i . The remaining integrals can be solved still by integrating by parts over q i . The integral proportional to Φ ,i in Eq. (159) gives where we have used the fact that dE/dq i = q i /E. For the integral the integration by parts brings two pieces, one from the derivation of q i q k and one from the derivation of the energy E

40
The last integral in Eq. (169) can indeed be neglected. To check this one makes use of the explicit expression (87) for the distribution function g to derive and Thus it is easy to compute which is negligible taking into account that T e /m e is of the order of the thermal velocity squared.
With these results we are now able to compute the left-hand side of the Boltzmann equation (160) integrated over d 3 q/(2π) 3 . The same operation must be done for the collision terms on the right hand side. For example for the first of the equations in (160) this brings to the integrals c ep QQ ′ qq ′ + c eγ pp ′ qq ′ . However looking at Eq. (86) one realizes that c eγ pp ′ qq ′ vanishes because the integrand is antisymmetric under the change q ↔ q ′ and p ↔ p ′ . In fact this is simply a consequence of the fact that the electron number is conserved for this process. The same argument holds for the other term c ep QQ ′ qq ′ . Therefore the right-hand side of Eq. (160) integrated over d 3 q/(2π) 3 vanishes and we can give the evolution equation for n e . Collecting the results of Eq. (166) to (172) we find Similarly, for CDM particles, we find

Momentum continuity equations
Let us now multiply Eq. (159) by (q i /E)/(2π) 3 and integrate over d 3 q. In this way we will find the continuity equation for the momentum of baryons. The first term just gives ( where we have used Eq. (171) and E = m e . The third term proportional to (H − Ψ ′ ) is where we have integrated by parts over q i . Notice that the last term in Eq. (176) is negligible being the same integral we discussewd above in Eq. (172). By the same arguments that lead to neglect the term of Eq. (172) it is easy to check that all the remaining integrals proportional to the gravitational potentials are negligible except for The integrals proportional to the second-order vector and tensor perturbations vanish as vector and tensor perturbations are traceless and divergence-free. The only one which survives is the term proportional to ω i′ + Hω i in Eq. (159). Therefore for the integral over d 3 qq i /E of the left-hand side of the Boltzmann equation (159) for a massive particle with mass m e (m p ) and distribution function (87) we find Now, in order to derive the momentum conservation equation for baryons, we take the first moment of both Eq. (160) and (161) multiplying them by q and Q respectively and integrating over the momenta. Since previously we integrated the left-hand side of these equations over d 3 qq i /E, we just need to multiply the previous integrals by m e for the electrons and for m p for the protons. Therefore if we sum the first moment of Eqs. (160) and (161) the dominant contribution on the left-hand side will be that of the protons Notice that the integral of the Coulomb collision term c ep (q i + Q i ) over all momenta vanishes simply because of momentum conservation (due to the Dirac function δ 4 (q + Q − q ′ − Q ′ )). As far as the Compton scattering is concerned we have that, following Ref. [42], still because of the total momentum conservation. Therefore what we can compute now is the integral over all momenta of c eγ p i . Notice however that this is equivalent just to multiply the Compton collision term C(p) of Eq. (86) by p i and integrate over d 3 p/(2π 3 ) where C(p) has been already computed in Eqs. (122) and (123). We will do the integral (181) in the following. First let us introduce the definition of the velocity of photons in terms of the distribution function where p γ = ρ γ /3 is the photon pressure and ρ γ the energy density. At first order we get where ∆ is the photon distribution anisotropies defined in Eq. (132). At second order we instead find Therefore the terms in Eqs. (122) and (123) proportional to f (1) (p) and f (2) (p) will give rise to terms containing the velocity of the photons. On the other hand the terms proportional to f where we are using the relations (165). Similarly all the terms proportional to v, (v·n) 2 and v 2 do not give any contribution to Eq. (181) and, in the second-order collision term, one can check that dΩY 2 (n)n i = 0.
0 /∂p for which we can use the rules (133) when integrating over p while the integration over the momentum direction is Finally from the second line of Eq. (123) we get three integrals. One is whereρ γ is the background energy density of the photons. The second comes from where we have used the rules (133), Eq. (165) and (dΩ/4π) n i n j n k n l = (δ ij δ kl + δ ik δ lj + δ il δ jk )/15. In fact the third integral exactly cancels the previous one. Summing the various integrals we find dp (2π) 3 C(p)p = n e σ Tργ Eq. (190) can be further simplified. Recalling that δ where the photon quadrupole Π ij γ is defined as Thus, our final expression for the integrated collision term (181) reads We are now able to give the momentum continuity equation for baryons by combining m p dg p /dη from Eq. (178) with the collision term (181) where ρ b is the baryon energy density and, as we previously explained, we took into account that to a good approximation the electrons do not contribute to the mass of baryons. In the following we will expand explicitly at first and second-order Eq. (194).

Second-order momentum continuity equation for baryons
At first order we find At second order there are various simplifications. In particular notice that the term on the right-hand side of Eq. (194) which is proportional to δ b vanishes when matched to expansion of the left-hand side by virtue of the first-order equation (195). Thus, at the end we find a very simple equation

Second-order momentum continuity equation for CDM
Since CDM particles are collisionless, at first order we find At second order we find 7 CMB anisotropies at second-order at all scales: analytical approach 7.1 Towards a second-order CMB radiation transfer function As pointed out in Sec. 6 various non-primordial sources of non-Gaussianity for the CMB anisotropies can arise from the non-linear evolution of the cosmological perturbations, including the Sunyaev-Z'eldovich, ISW (Rees-Sciama) effects, and the gravitational lensing, with possible correlations between these contributions. Here we will focus on another relevant source of non-Gaussianity: the non-linear effects operating at the recombination epoch. The dynamics at recombination is quite involved because all the non-linearities in the evolution of the baryon-photon fluid at recombination and the ones coming from general relativity should be accounted for.
The following section can be considered as an application of the Boltzmann equations found previously. Despite they depend on a high numbers of terms, because of the apperance within second-order perturbation theory of products of fist-oder perturbations, it is remarkable that an analytical study of the complicated recombination dynamics is possible. This allows to account for those effects that at the last scattering surface produce a non-Gaussian contribution to the CMB anisotropies that add to the primordial one. Such a contribution is so relevant because it represents a major part of the second-order radiation transfer function which must be determined in order to have a complete control of both the primordial and non-primordial part of NG in the CMB anisotropies and to gain from the theoretical side the same level of precision that could be reached experimentally in the near future [5,4,58].
In order to achieve this goal, we consider the Boltzmann equations derived in Sec. 6 at second-order describing the evolution of the photon, baryon and CDM fluids, and we manipulate them further under the assumption of tight coupling beteween photons and baryons. This leads to the generalization at secondorder of the equations for the photon energy density and velocity perturbations which govern the acoustic oscillations of the photon-baryon fluid for modes that are inside the horizon at recombination. The evolution is that of a damped harmonic oscillator, with a source term which is given by the gravitational potentials generated by the different species. An interesting result is that, unlike the linear case, at second-order the quadrupole moment of the photons is not suppressed in the tight coupling limit and it must be taken into account.
The analytical solutions for the acoustic oscillations of the photon-baryon fluid at second-order are derived adopting some simplifications which are also standard for an analytical treatment of the linear CMB anisotropies, and which nonetheless allow to catch most of the physics at recombination. One of these simplifications is to study separately two limiting regimes: intermediate scales which enter the horizon in between the equality epoch (η eq ) and the recombination epoch (η r ), with η −1 r ≪ k ≪ η −1 eq , and shortwave perturbations, with k ≫ η −1 eq , which enter the horizon before the equality epoch. Here our main concern is to provide a simple estimate of the quantitative behaviour of the non-linear evolution taking place at recombination, offering at the same time all the tools for a more accurate computation. We find that the second-order CMB anisotropies generated on the last scattering surface do not reduce only to the energy density and velocity perturbations of the photons evaluated at recombination, but a number of second-order corrections at last scattering arise from the Boltzmann equations of Ref. [38] in the form of products of first-oder perturbations.
We will see that the dynamics at recombination is indeed dominated on small scales by the nonlinear evolution of the second-order gravitational potentials feeded by the Cold Dark Matter density perturbations. The gravitational potentials determine the energy density fluctuations of the photons at recombination and their effects show up in the CMB anisotropies as which is the usual term due to the intrinsic fractional temperature fluctuation ∆ (2) 00 /4 on the last scattering surface (∆ (2) 00 is the monople of the photon distribution function) and the gravitational redshift due to the gravitational potential. However the analysis of the remaining contributions that come in the form of products of first-oder perturbations is equally important. The reason is that one of the central quantities we are interested in is the contamination to the primordial non-Gaussianity that is produced by the secondary effects. In that respect the reasonable question to ask is which kind of primordial non-Gaussianity a given secondary effect can contaminate most, wether it is of the so called "local type" or of the "equilateral type", for example. Therefore it might well be the case that, even if some secondary effects appear to be the dominant ones, they might give a high contamination to a given type of primordial non-Gaussianity, but a low contamination to a different kind of primordial non-Gaussianity, for which "subdominant " terms, on the other hand, represent the strongest contaminant. This section and the following two provide a clear example: we anticipate here that while the term in Eq. (199) will mainly mimic an equilateral type of primordial non-Gaussianity, second-order effects that come as products of first-order times first-order perturbations actually mainly contaminate local primordial non-Gaussianity.
Notice that the case k ≫ η −1 eq has been treated in two steps. First we just assume a radiation dominated universe, and then we provide a much more complete analysis by solving the evolution of the perturbations from the equality epoch onwards taking into account that the dark matter perturbations around the equality epoch tend to dominate the second-order gravitational potentials. As a byproduct, this last step provides the Meszaros effect at second-order. In deriving the analytical solutions we have accurately accounted for the initial conditions set on superhorizon scales by the primordial non-Gaussianity. In fact the primordial contribution is always transferred linearly, while the real new contribution to the radiation transfer function is given by all the additional terms provided in the source functions of the equations. Let us stress here that the analysis of the CMB bispectrum performed so far, as for example in Ref. [59], adopt just the linear radiation transfer function (unless the bispectrum originated by specific secondary effects, such as Rees-Sciama or Sunyaev-Zel'dovich effects, are considered, see on this Sec. 8). The formalism to provide in a systematic way the second-order CMB radiation transfer function will be reviewd in Sec. 10.

The Boltzmann equations in the tightly coupled limit
We now derive the moments of the Boltzmann equations for photons in the limit when the photons are tightly coupled to the the baryons (the electron-proton system) due to Compton scattering. This leads to the governing equations for the acoustic oscillations of the photon-baryon fluid. The well-known computation at linear order is briefly reviewed in C.1 under some standard simplifying assumptions. Here we will focus on the derivation of the equations at second-order in the perturbations, pointing out some interesting differences with respect to the linear case. In particular note that ,while we already know that the L.H.S. of the equations at second-order will have the same form as for the linear case, the source term on the R.H.S. of the moments of the Boltzmann equations will also consist of first-order squared terms.

Energy continuity equation
Let us start by integrating Eq. (134) over dΩ n /4π to get the evolution equation for the second-order photon energy density perturbations ∆ (2) 00 ∆ (2)′ 00 ,j + Ψ ,j )n i n j − (Φ ,i + Ψ ,i ) where we have used the explicit definition for the second-order velocity of the photons [38] 4 We can now make use of the tight coupling expansion to simplify Eq. (200). In the L.H.S. we use by using Eq. (C.4) and the evolution equation for the photon veleocity (C.5). Here we introduce the baryon-photon ratio We thus arrive at the following equation where the source term is given by Notice that the integral in the second line of Eq. (200) has been computed by expanding the linear anisotropies as ∆ (1) = ℓ (2ℓ+1)∆ ℓ P ℓ (v·n) and expressing explicitly the dependence on n i as P 1 (v·n) = v i n i , P 2 (v · n) = (3(v · n) 2 − 1)/2 and so on. This allows to perform the derivative ∂∆ (1) /∂n i . It turns out that the term proportional to the first-order quadrupole vanishes, and higher order terms can be neglected because of the tight coupling approximation. Therefore the only term that is non negligible is the dipole term and using ∆ (1) 1 = 4v/3 in the tight coupling limit, one obatins that the integral is equal to (16/

Velocity continuity equation
We now derive the second moment of the Boltzmann equation (134) and then we take its tight coupling limit. The integration of Eq. (134) over dΩ n n i /4π yields the continuity equation for the photon velocity ,j + Ψ ,j )n i n j − (Φ ,i + Ψ ,i ) The difference between the second-order baryon and photon velocities (v (2)i −v (2)i γ ) appearing in Eq. (206) is obtained from the baryon continuity equation which can be written as (see Ref. [38]) We want now to reduce Eq. (206) in the tightly coupled limit. We first insert the expression (207) in Eq. (206). Notice that the last three terms in Eq. (207) will cancel out. On the other hand in the tight coupling limit expansion one can set v (1) where in the tightly coupled limit we are neglecting the first-order quadrupole and (higher-order moments) of the photon distribution since it is suppressed by 1/τ with respect to the other terms. Moreover the integral in the fourth line of Eq. (206) has been computed following the same steps used for the analogous integral in Eq. (200). In this case it turns out that the contribution from the first-order dipole term vanishes, while the one from thequadrupole is non zero, but this and higher-order terms can be neglected because of the tight coupling approximation, so that the integral is in fact negligible. Next for the term like τ ′ δ (1) 00 /4 and we use the first order tight coupling equations (C.1) and (C.5) in order to further simplify Eq. (208). We finally obtain v (2)i′ where We have spent some time in giving the details of the computation for the photon Boltzmann equations at second-order in the perturbations. As a summary of the results obtained so far we refer the reader to Eqs. (204) and (209) as our master equations which we will solve in the next sections. In particular Eq. (209) is the second-order counterpart of Eq. (C.5) for the photon velocity in the tight coupling regime. Notice that there are two important differences with respect to the linear case. One is that, in Eq. (209), there will be a contribution not only from scalar perturbations but also from vector modes which, at second-order, are inevitably generated as non-linear combinations of first-order scalar perturbations. In particular we have included the vector metric perturbations ω i in the source term. Second, and most important, we have also kept in the source term the second-order quadrupole of the photon distribution Π (2)ij γ . At linear order we can neglect it together with higher order moments of the photons since they turn out to be suppressed with respect to the first two moments in the tight coupling limit by increasing powers of 1/τ . However in the next section we will show that at second order this does not hold anymore, as the photon quadrupole is no longer suppressed.
Finally following the same steps leading to Eq. (C.6) at linear order we can derive a similar equation for the second-order photon energy density perturbation ∆ (2) 00 which now will be characterized by the source terms S ∆ and S i where we have introduced the photon-baryon fluid sound speed c s = 1/ 3(1 + R).

Second-order quadrupole moment of the photons in the tight coupling limit
Let us now consider the quadrupole moment of the photon distribution defined in Eq. (192) and show that at second-order it cannot be neglected in the tightly coupled limit, unlike for the linear case. We first integrate the R.H.S. of Eq. (134) over dΩ n (n i n j − δ ij /3)/4π and then we set it to be vanishing in the limit of tight coupling. The integration involves various pieces to compute. For clarity we will consider each of them separately. The term ∆ so that the integral just brings Π (2)ij γ /10, since the only contribution in Eq. (212) comes from ∆ (2) /10 with all the other terms vanishing. The following nontrivial integral is where the baryon velocity appearing in P 2 (v · n) is first order and we make use of the relations (165) together with dΩ 4π n i n j n k n l = 1 15 (δ ij δ kl + δ ik δ lj + δ il δ jk ) .
The integrals of δ e (v · n) and v (2) · n vanish and where in the last step we take ∆ 1 = 4v/3 in the tight coupling limit. Similarly the integral of 14(v · n) 2 brings The integral of 2(v · n)∆ (1) can be performed by expanding the linear anisotropies as ∆ (1) where we have used Eq. (165) and O ℓ>2 indicates all the integrals coming from the multipoles ℓ > 2 in the expansion (for ℓ = 0 and ℓ = 2 they vanish.) In fact we have dropped the O(ℓ > 2) since they are proportional to firts-order photon moments ℓ > 2 which turn out to be suprressed in the tight coupling limit. Finally the term proportional to (v · n)∆ 2 (1 − P 2 (v · n)/5 gives a vanishing contribution. Collecting all the various pieces we find that the third moment of the R.H.S. of Eq. (134) is given by Therefore in the limit of tight coupling, when the interaction rate is very high, the second-order quadrupole moment turns out to be by setting Eq. (218) to be vanishing (the term multiplying δ e goes to zero in the tight coupling limit since it just comes from the first-order collision term). At linear order one would simply get the term 9τ ′ Π (1)ij γ /10 implying that, in the limit of a high scattering rate τ ′ , Π (1)ij γ goes to zero. However at second-order the quadrupole is not suppressed in the tight coupling limit becasue it turns out to be sourced by the linear velocity squared.

Second-order CMB anisotropies generated at recombination
The previous equations allow us to follow the evolution of the monopole and dipole of CMB photons at recombination. As at linear order, they will appear in the expression for the CMB anisotropies today ∆ (2) (k, n, η 0 ) together with various integrated effects. Our focus now will be to obtain an expression for the second-order CMB anisotropies today ∆ (2) (k, n, η 0 ) from which we can extract all those contributions generated specifically at recombination due to the non-linear dynamics of the photon-baryon fluid. This expression will not only relate the moments ∆ (2) ℓm today to the second-order monopole and dipole at recombination as it happens at linear order, but one has to properly account also for additional firstorder squared contributions. Let us see how to achieve this goal in some details.
As we have seen previously, it is possible to write down an integral solution of the photon Boltzmann equation (134) in Fourier space. We wrote in order to derive a solution of the form Here µ = cosϑ =k · n is the polar angle of the photon momentum in a coordinate system such that e 3 =k. At second-order the source term has been computed in Ref. [38] and can be read off Eq. (134) and Eq. (138) to be The key point here is to isolate all those terms that multiply the differential optical depth τ ′ . The reason is that in this case in the integral (221) one recognizes the visibility function g(η) = −e −τ τ ′ which is sharply peaked at the time of recombination and whose integral over time is normalized to unity. Thus for these terms the integral just reduces to the remaining integrand (apart from the visibility function) evaluated at recombination. The standard example that one encounters also at linear order is given by the first term appearing in the source S, Eq. (222), that is −τ ′ ∆ (2) 00 . The contribution of this term to the integral (221) just reduces to where η * is the epoch of recombination and, in the multipole decomposition (136), Eq. (223) brings the standard result ∆ having used the Legendre expansion e ik·x = ℓ (i) ℓ (2ℓ + 1)j ℓ (kx)P ℓ (k ·x). In Eq. (224) the monopole at recombination is found by solving the Boltzmann equations Eqs. (204)-(209) derived in tight coupling limit.
Looking at Eq. (222) we recognize immediately some terms which multiply explicitly τ ′ (the first one discussed in the example above and the last two lines of Eq. (222)). However it is easy to realize from the standard procedure adopted at the linear-order that such terms are not the only ones. This is clear by focusing, as an example, on the term −4n i Φ (2) ,i in the source S which appears in the same form also at linear order. In Fourier space one can replace the angle µ with a time derivative and thus this term gives rise to − 4ik where, in the last step, we have integrated by parts. In Eq. (225) the time derivative of the gravitational potential contributes to the Integrated Sachs-Wolfe effect, but also also a τ ′ results implying that we have also to evaluate Φ (2) at recombination. Thus, in the following, we look for those terms in the source (222) which give rise to a τ ′ factor in the same way as for −4n i Φ ,i . In particular let us consider the combination in Eq. (222) where for simplicity we are setting Φ (1) ≃ Ψ (1) . We already recognize terms of the form n i ∂ i (·). Moreover we can use the Boltzmann equation (129) to replace n i ∆ ,i in Eq. (226). This brings In fact we will not be interested for our purposes in the first three terms of Eq. (227), since they will not contribute to the anisotropies generated at recombination.
Therefore, as a result of Eqs. (222), (225) and (227), we can rewrite the source term (222) as and In Eq. (228) S * contains the contribution to the second-order CMB anisotropies created on the last scattering surface at recombination, while S ′ includes all those effects which are integrated in time from the last scattering surface up to now, including the second-order Integrated Sachs-Wolfe effect and the second-order lensing effect. Since the main concern of this Section is the CMB anisotropies generated at last scattering, from now on we will focus only on the contribution from the last scattering surface S * . In particular notice that, following the same steps that lead to Eq. (223), the first two terms of Eq. (229) give rise to the CMB temperature anisotropies as written in Eq. (199).

Tightly coupled solutions for the second-order perturbations
In this section we will solve the tightlty coupled limit of the Boltzmann equations (204) and (206) at second-order in perturbation theory. We will proceed as for the linear case, focusing on the two limiting cases of perturbation modes entering the horizon respectively much before and much after the time of equality. The solution of Eq. (211) can be written as where for simplicity we use the shorthand notation sin k [η, η ′ ] ≡ sin[k(r s (η) − r s (η ′ ))], and r s (η) is the sound horizon The source terms are given in Eq. (205) and (210). Notice that we can write S ′ ∆ + HR 1+R S ∆ = (S ∆ (1 + R)) ′ /1 + R so that we can perform an integration by parts in Eq. (231) leading to In order to give an analytical solution that catches most of the physics underlying Eq. (231) and which remains at the same time very simple to treat, we will make some simplifications which are also usually adopted in linear theory, see e.g. Ref. [87]. In particular we will treat R = R * as a constant evaluated at the time of recombination. In this way, with the presence of R in the varying cosines and sines, we can keep track of a damping of the photon velocity amplitude with respect to the case R = 0 which prevents the acoustic peaks in the power-spectrum to disappear. Treating R as a constant is justified by the fact that for modes within the horizon the time scale of the ocillations is much shorter than the time scale on which R varies. If R is a constant the sound speed is just a constant and the sound horizon is simply r s (η) = c s η.
There is also another good reason not to forget about R, even tough R < 1, which is particularly relevant when estimating the second-order contributions to the CMB anisotropies in (199). As we will see in Sec. 9, that combination on small scales comes actually as proportional to R, so that neglecting R would lead to miss one of the most significant CMB contribution at second order.
Another simplification consists in solving the evolutions of the perturbations in two well distinguished limiting regimes. One regime is for those perturbations which enter the Hubble radius when matter is the dominant component, that is at times much bigger than the equality epoch, with k ≪ k eq ∼ η −1 eq , where k eq is the wavenumber of the Hubble radius at the equality epoch. The other regime is for those perturbations with much smaller wavelenghts which enter the Hubble radius when the universe is still radiation dominated, that is perturbations with wavenumbers k ≫ k eq ∼ η −1 eq . In fact we are interested in perturbation modes which are within the horizon by the time of recombination η * . Therefore we will further suppose that η * ≫ η eq in order to study such modes in the first regime. Even though η * ≫ η eq is not the real case, it allows to obtain some analytical expressions.

Setting the initial conditions: primordial non-Gaussianity
The integration constants A and B are fixed according to the initial conditons for the second-order cosmological perturbations. These refer to the values of the perturbations on superhorizon scales deep in the radiation dominated period. We will consider the case of initial adiabatic perturbations, for which there exist some useful conserved quantities on large scales which as such carry directly the information about the initial conditons.
As explained in Sec. 4 such a conserved quantity is given by the curvature perturbation ζ and its conserved value allows to set the initial conditions for the metric and matter perturbations accounting for the primordial contributions. At linear order during the radiation-dominated epoch and on large scales ζ (1) = −2Ψ (1) /3. On the other hand, after some calculations, one can easily compute ∆ζ (2) for a radiation dominated epoch where in Eq. (38) one uses that on large scales δ (1) ρ γ /ρ γ = −2Ψ (1) and the energy continuity equation where we are evaluating the quantities in the large scale limit for η → 0. Using the parametrization (39) at the initial times the quantity ∆ 00 − 4Ψ (2) is given by Since for adiabatic perturbations such a quantity is conserved on superhorizon scales, it follows that the constant B = 0 and A = 2(9a NL − 7)Ψ (1)2 (0). Eqs. (231) and (233) are analytical expressions describing the acoustic oscillations of the photonbaryon fluid induced at second-order for perturbation modes within the horizon at recombination. In the following we will adopt similar simplifications already used for the linear case in order to provide some analytical solutions. In particular, if in Eq. (233) we treat R as a constant we can write, using the initial conditions determined above,

Perturbation modes with k ≫ k eq
In order to study the contribution to the second-order CMB anisotropies coming from perturbation modes that enter the horizon during the radiation dominated epoch, we will assume that the second-order gravitational potentials are the ones of a pure radiation dominated universe throughout the evolution.
Though not strictly correct, this approximation will give us the basic picture of the acoustic oscillations for the baryon-photon fluid occuring for these modes. Also for the second-order case, in Section 7.7 we will provide the appropriate corrections accounting for the transition from radiation to matter domination which is indeed (almost) achieved by the recombination epoch. Before moving into the details a note of caution is in order here. At second order in the perturbations all the relevant quantities are expressed as convolutions of linear perturbations, bringing to a mode-mode mixing. In some cases in our treatment for a given regime under analysis (k ≫ k eq or k ≪ k eq ) we use for the first-order perturbations the solutions corresponding to that particular regime, while the mode-mode mixing would require to consider in the convolutions (where one is integrating over all the wavenumbers) a more general expression for the first-order perturbations (which analytically does not exist anyway). For the computation of the CMB bispectrum this would be equivalent to consider just some specific scales, i.e. all the three scales involved in the bispectrum should correspond approximately to wavenumbers k ≫ k eq or k ≪ k eq , and not a combination of the two regimes (a step towards the evaluation of the three-point correlation function has been taken on Ref. [37] where it was computed in the in so-called squeezed triangle limit, when one mode has a wavelength much larger than the other two and is outside the horizon).
In Appendix C we show how to solve the Boltzmann equations at the linear level in the various regimes. In fact for k ≫ k eq the acoustic oscillations can be solved in an alternative way (see Ref. [39]) for this procedure at linear order). One can start directly from Eq. (204) where we can neglect the gravitational potential term Ψ (2) ′ . The reason is that, as it happens at linear order, the second-order gravitational potentials decay at late times as η −2 , while the second-order velocity v (2)i γ oscillates in time.
Let us now see that in some details.
The evolution equation for the gravitational potential Ψ (2) is given by Eq. (B.13) and is characterized by the source term S γ , Eq. (B.14). In particular the source term contains the second-order quadrupole moment of the photons Π (2)ij γ . We saw in Section 7.2.3 that at second-order the quadrupole moment is not suppressed in the tight coupling limit, being fed by the non-linear combination of the first-order velocities, Eq. (219). For the pertubation modes we are considering here the velocity at late times is oscillating being given in Appendix C in Fourier space. Since the linear gravitational potential decays in time and for a radiation dominated period H = 1/η, it is easy to check that the dominant contribution at late times to the source term S γ simply reduces to where and the sound speed is c s = 1/ 3(1 + R). Before proceeding further let us explain the notation that we are using. The equivalence symbol will be used to indicate that we are evaluating the expression in Fourier space. At second-order in perturbation theory most of the Fourier transforms reduce to some convolutions. We will not indicate these convolutions explicitly but just through their kernel. For example in Eq. (239) by F (k 1 , k 2 , k) we actually indicate the convolution operator In the specific case of Eq. (239) the kernel is given by The choice of these conventions is due not only for simplicity and to keep our expressions shorter, but also because at the end we will be interested to the bispectrum of the CMB anisotropies generated at recombination, and the relevant expressions entering in the bispectrum are just the kernels of the convolution integrals.
Having determined the leading contribution to the source term at late times, we can now solve the evolution equation (B.13). Since the source term scales like η −2 it is useful to introduce the rescaled variable χ = η 2 Ψ (2) . Eq. (B.13) then reads For perturbation modes which are subhorizon with kη ≫ 1 the solution of the homogeneous equation is given by from which we can build the general solution where W = −kc s is the Wronskian, and χ + = cos(kc s η), χ − = sin(c s kη). Using Eq. (239) the integrals involve products of sines and cosines which can be performed giving Thus the gravitational potential Ψ (2) at late times is given by where we have set the integration constant B = 0 and A = −3Ψ (2) (0)/(kc s ) 2 in order to match the homogenoeus solution at late times. Here Ψ (2) (0) is the intial condition for the gravitational potential taken on large scales deep in the radiation dominated era which will be determined in Section 7.5.2.
Eq. (247) shows the result that we anticipated: also at second order the gravitational potential varies in time oscillating with an amplitude that decays as η −2 . Let us then take the divergence of the (i − 0)

Einstein equation (31) expanded at second-order
which, using the first-order (i − 0) Einstein equation and Φ (1) ≃ Ψ (1) , reduces to Since Ψ (1) during a radiation dominated period is given by Eq. (B.11) and at late times it decays oscillating, it is easy to see that (Ψ (1)′ ∂ i Ψ (1) ) will be oscillating and decaying as η −4 and thus can be neglected with respect to Ψ (2)′ , which oscillates with an amplitude decaying as η −2 . Also HΦ (2) turns out to be subdominant. Recall that Φ (2) = Ψ (2) − Q (2) and Q (2) is dominated by the second-order quadrupole of the photons in Eq. (B.14), so that Φ (2) scales like Ψ (2) but there is the additional damping factor of the Hubble rate H = 1/η. Thus the dominant terms give Eq. (250) allows to proceed further in a similar way as for the linear case by using the results found so far, Eqs. (247) which, using the first-order equation (C.1), further simplifies to where we have kept only the dominant terms at late times. The gravitational potential Ψ (2) is given in Eq. (247), so the integration of Eq. (252) gives Needless to say, modes for k ≫ k D , where k −1 D indicates the usual damping length, are supposed to be multiplied by an exponential e −(k/k D ) 2 (see, e.g. [42]).

Vector perturbations
So far we have discussed only scalar perturbations. However at second-order in perturbation theory an unavoidable prediction is that also vector (and tensor) perturbation modes are produced dynamically as non-linear combination of first-order scalar perturbations. In particular notice that the second-order velocity appearing in Eq. (229), giving rise to a second-order Doppler effect at last scattering, will contain a scalar and a vector (divergence free) part. Eq. (250) provides the scalar component of the second-order velocity. We now derive an expression for the velocity that includes also the vector contribution.
The (second-order) vector metric perturbation ω i when radiation dominates can be obtained from Eq. (B.16) where we have dropped the gravitational potentials Ψ (1) ≃ Φ (1) which are subdominant at late times.
On the other hand from the velocity continuity equation (209) neglecting the term proportional to R and the decaying gravitational potentials. Using the tight coupling equations at first-order and integrating over time one finds We can thus plug Eq. (256) into Eq. (254) to find that at late times (for kη ≫ 1) We will come later to the explicit expression for the term on the R.H.S. of Eq. (257). Here it is enough to notice that the second-order quadrupole oscillate in time and thus ω i will decay in time as H 2 = 1/η 2 . This shows that ω i in Eq. (256) can be in fact neglected with respect to the other terms giving It can be useful to compute the combination on the R.H.S. of Eq. (257) (δ i j − ∂ i ∂ j /∇ 2 )∂ k Π (2)kj γ . The second-order quadrupole moment of the photons in the tightly coupled limit is given by Eq. (219), and where in the last step we have used that the linear velocity is the gradient of a scalar perturbation. We thus find Notice that if we split the quadrupole moment into a scalar, vector (divergence-free) and tensor (divergencefree and traceless) parts as then it turns out that where Π (2)i γ is the vector part of the quadrupole moment. Therefore one can rewrite Eq. (257) as

Initial conditions for the second-order gravitational potentials
In order to complete the study of the CMB anisotropies at second-order for modes k ≫ k eq we have to specify the initial conditions Ψ The conserved value of ζ (2) is parametrized by ζ (2) = 2a NL ζ (1)2 , where, as explained in Section 7.4 the parameter a NL specifies the level of primordial non-Gaussianity depending on the particular scenario for the generation of the cosmological perturbations. On the other hand , during radiation-domination where where we are evaluating Eq. (B.14) in the limit kη ≪ 1. The contribution from the second-order quadrupole moment in this limit reads where F and C are defined in Eqs. (241) and (240). Therefore we find that in Fourier space and from Eq. (264) we read off the intial condition as (convolution products are understood)

Perturbation modes with k ≪ k eq
Let us consider the photon perturbations which enter the horizon between the equality epoch and the recombination epoch, with wavelenghts η −1 * < k < η −1 eq . In fact, in order to find some analytical solutions, we will assume that by the time of recombination the universe is matter dominated η eq ≪ η * . In this case the gravitational potentials are sourced by the dark matter component and their evolution is given in Sec .B.1. At linear order the gravitational potentials remain constant in time, while at second-order they are given by Eq. (B.6). In turn the gravitational potentials act as an external force on the CMB photons as in the equation (211) describing the CMB energy density evolution in the tightly coupled regime.
For the regime of interest it proves convenient to use the solution of Eq. (211) found in (238). The source functions S ∆ and S i V are given by Eqs. (205) and (210), respectively. In particular S ∆ at early times -S ∆ (0) appearing in Eq. (238)-vanishes. For a matter dominated period where we have used the linear evolution equations (C.1) and (C.5) with Φ (1) = Ψ (1) , and As at linear order we are evaluating these expressions in the limit In S i V we have used the expression (219) for the second-order quadrupole moment Π (2)ij γ of the photons in the tight coupling limit, with the velocity v (1) = v (1) γ . Notice that, for the modes crossing the horizon at η > η eq , we have expressed the gravitational potential during the matter dominated period in terms of the initial value on superhorizon scales deep in the radiation dominated epoch as Ψ (1) = 9Ψ (1) (0)/10.
As for the second-order gravitational potentials we have to compute the combination Φ (2) + Ψ (2) appearing in Eq. (238). The gravitational potential Ψ (2) is given by Eq. (B.6), while Φ (2) is given by according to the relation (33), where for a matter-dominated period We thus find which in Fourier space reads where the kernels of the convolutions are given by Eq. (242) and In Eq. (276) Ψ (2) m (0) is the initial condition for the gravitational potential fixed at some time η i > η eq . For the regime of interest it corresponds to the value of the gravitational potential on superhorizon scales during the matter-dominated epoch.
Notice a property that will be useful later on. By looking at Eqs. (275) and the explicit solution for Ψ (2) (B.6) which grows as η 2 , it is easy to realize that on very small scales, for kη ≫ 1, the two gravitational potentials are equal with We are now able to compute the integrals entering in the solution (238). The one involving the second-order gravitional potentials is straightforward to compute where for brevity F and G stand for the convolutions (242) and (278). Notice also the we have isolated the terms proportional to R as it will be useful later on. For the two remaining integrals, in the following we will show only the terms that in the final expression for ∆ (2) 00 and the second-order velocity v (2)i γ give the dominant contributions for kη ≫ 1, even though we have perfomed a fully computation. The integral over the source function S ∆ yields a sum of oscillating functions (cosines) which turn out to be subdominant, so we do not display the full result. For the last integral we find ((1 ↔ 2) stands by an exchange of indices) where the terms that have been dropped vary in time as sines and cosines. We have written the contribution in Eq. (281) because, upon integration over time, it will give a non-negligible contribution to the velocity v (2)i γ . From the general solution (238) and the expression (B.6) for the second-order gravitational potential Ψ (2) we thus obtain ∆ (2) 00 We warn the reader that in writing Eq. (282) we have kept all those terms that contain the primordial non-Gaussianity parametrized by a NL , and the terms which dominate at late times for kη ≫ 1.

Initial conditions for the second-order gravitational potentials
The initial condition Ψ (2) m (0) for the modes that cross the horizon after the equality epoch is fixed by the value of the gravitational potential on superhorizon scales during the matter dominated epoch. To compute this value we use the conservation on superhorizon scales of the curvature perturbation ζ (2) defined in Eq. (37). For a matter-dominated period the curvature perturbation on large-scales turns out ot be where we used the energy continuity equation (1) in the superhorizon limit.
From the (0 − 0) Einstein equation on large scales δ (2) The conserved value of ζ (2) is parametrized as in Eq.
we have expressed the gravitational potential during the matter dominated period Ψ (1) in terms of the initial value on superhorizon scales after the equality epoch as Ψ where F is the kernel defined in Eq. (242). We can use the explicit expression for Ψ (2) m (0) in Eq. (282), still keeping only the terms that contain the primordial non-Gaussianity parametrized by a NL , and the terms which dominate at late times for kη ≫ 1 to find ∆ (2) 00

Second-order photon velocity perturbation
The second-order velocity of the photons can be obtained from Eq. (209) Notice that for simplicity in writing Eq. (288) we have dropped off the dependence on R. In fact the main conclusions of this Subsection remains unchanged. The second-order gravitional potential in matter-dominated universe can be obtained from Eqs. (274)-(275) and Eq. (B.6) as In Fourier space this becomes where the kernels of the convolutions are given by Eqs. (242) and (278). The integral over Φ (2) in Eq. (288) is then easily computed where as usual the equivalence symbol means that we are evaluating a given expression in Fourier space. For the integral over the source function S i V we use its expression in Fourier space, Eq. (273), and the dominant terms for kη ≫ 1 are Notice that, in order to compute this integral, we must know the second-order vector metric petturbation ω i . This is easily obtained for a matter-dominated universe from Eq. (B.7). Using Eqs. (B.4) and (B.2) one finds giving rise to the third term in Eq. (292). Finally for the integral over ∆ (2) 00 some caution is nedeed. Since in the final expression for v (2)i γ the dominant terms at late times turn out to be proportional η, one has to use an expression for ∆ (2) 00 that keep track of all those contributions that, upon integration, scale like η. Thus we must use the expression written in Eq. (282), plus Eq. (281), and some terms of Eq. (280) that have been previously neglected in Eq. (282). Then we find for kη ≫ 1 Using Eqs. (291), (292) and (294) To obtain Eq. (295) we have also used the explicit expression (286) for Ψ (2) m (0) and we have kept the terms depending on a NL parametrizing the primordial non-Gaussianity and the terms that dominate at late times for kη ≫ 1.

Perturbation modes with k ≫ k eq : Cold Dark Matter perturbations at second order and improved analytical solutions for CMB anisotropies
In Sec. 7.5 we have computed the perturbations of the CBM photons at last scattering for the modes that cross the horizon at η < η eq under the approximation that the universe is radiation-dominated. However around the equality epoch, through recombination, the dark matter component will start to dominate. In this section we will account for its contribution to the gravitational potential and for the resulting perturbations of the photons from the equality epoch onwards. This leads to a more realistic and accurate abalytical solutions for the acoustic oscillations of the photon-baryon fluid for the modes of interest. The starting point is to consider the density perturbation in the dark matter component for subhorizon modes during the radiation dominated epoch. Its value at the equality epoch will fix the magnitude of the gravitational potential at η eq and hence the initial conditions for the subsequent evolution of the photons fluctuations during the matter dominated period. At linear order the procedure is standard (see, e.g [42]), and we will use a similar one at second-order in the perturbations. Thus this section serves also as a guide through the evolution of the CDM density perturbations at second-order accounting for those modes that enter the horizon during the radiation dominated epoch. This allows to determine the second-order transfer function for the density perturbations with a generalization at second-orde of the Meszaros effect. 7.7.1 Subhorizon evolution of CDM perturbations for η < η eq From the energy and velocity continuity equations for CDM it is possible to isolate an evolution equation for the density perturbation δ d = δρ d /ρ d , where the subscript d stands for cold dark matter. In Ref. [38] we have obtained the Boltzmann equations up to second-order for CDM. The number density of CDM evolves according to [38] At linear order n d =n d + δ (1) n d and one recovers the usual energy continuity equation with δ The CDM velocity at the same order of perturbation obeys [38] v (1) Perturbing n CDM up to second-order we find The R.H.S. of this equation can be further manipulated by using the linear equation (297) to replace v where we use Φ (1) = Ψ (1) . In Ref. [38] the evolution equation for the second-order CDM velocity perturbation has been already obtained At linear order we can take the divergence of Eq. (298) and, using Eq. (297) to replace the velocity perturbation, we obtain a differental equation for the CDM density contrast which can be rewritten as where When the radiation is dominating the gravitational potential is mainly due to the perturbations in the photons, and a(η) ∝ η. For subhorizon scales Eq. (303) can be solved following the procedure introduced in Ref. [?]. Using the Green method the general solution to Eq. (303) (in Fourier space) is given by where the first two terms correspond to the solution of the homogeneous equation. At early times the density contrast is constant with having used the adiabaticity condition, and thus we fix the integration constant as and C 2 = 0. The gravitational potential during the radiation-dominated epoch starts to decay as a given mode enters the horizon. Therefore the source term S (1) behaves in a similar manner and this implies that the integrals over η ′ reach asymptotically a constant value. Once the mode has crossed the horizon we can thus write the solution as where the constants A (1) and B (1) are defined as and The upper limit of the integrals can be taken to infinity because the main contribution comes form when kη ∼ 1 and once the mode has entered the horizon the result will change by a very small quantity. Performing the integrals in Eq. (309) and (310) one finds that A (1) = −9.6 and B (1) ≃ 0.44. Before moving to the second-order case, a useful quantity to compute is the CDM velocity in a radiation dominated epoch. From Eq. (298) it is given by where the last equality holds in Fourier space and we have used Eq. (B.11) (and the fact that a(η) ∝ η when radiation dominates).
Combining Eq. (300) and (301) we get the analogue of Eq. (303) at second-order in perturbation theory δ where the source function is In fact we write Eq. (312) in a more convenient way as where for simplicity we have introduced the two functions and In this way we get an equation of the same form as (303) in the variable [δ (2) with source s 2 on the R.H.S.. Its solution in Fourier space therefore is just as Eq. (305) As we will see, Eq. (317) provides the generalization of the Meszaros effect at second-order in perturbation theory.

Initial conditions
In the next two sections we will compute explicitly the expression (317) for the second-order CDM density contrast on subhorizon scales during the radiation dominate era. First let us fix the constants C where in the last step we have used Eq. (306). Therefore we can fix C 2 = 0 and Eq. (237) gives ∆ (2) 00 (0) − 4Ψ (2) (0) in terms of the primordial non-Gaussianity parametrized by a NL , and the expression for Ψ (2) (0) have been already computed in Eq. (269). Thus we find (in Fourier space) ∆ (2) 00 = 2(3a NL − 1) + 44 and from Eqs. (318) we derive the initial density contrast for CDM at second-order Eq. (321) togheter with Eq. (269) allows to compute the constant C 1 as

Meszaros effect at second order
We now compute the integrals over the functions s 1 and s 2 appearing in Eq. (317). Let us first focus on the integral η 0 dη ′ s 1 (η ′ ). Notice that, using the linear equations (297) and (298) for the CDM density and velocity perturbations, the function s 1 (η ′ ) can be written in a more convenient way as and then In Eq. (324) all the quantities are known being first-order perturbations: the linear gravitional potential Ψ (1) for a radiation dominated era is given in Eq. (B.11), the CDM velocity perturbation corresponds to Eq. (311) and the CDM density contrast is given by Eq. (308). Thus the integral in Eq. (324) reads (in Fourier space) Let us recall that we are interested in the evolution of the CDM second-order density contrast on subhorizon scales during the radiation dominated epoch. Therefore once we compute the integrals we are interested in the limit of their expression for late times (kη ≫ 1). For the first contribution to Eq. (325) we find that at late times it is well approximated by the expression We have computed also the remaining integral in Eq. (325), but it turns out to be negligible compared to Eq. (326). The reason is that the integrand oscillates with an amplitude decaying in time as η −3 , and thus it leads just to a constant (the argument is the same we used at linear order to compute the integrals in Eq. (305)). Thus we can write We now compute the integrals over the function s 2 (η) given in Eq. (316). Since at late times φ (1)2 ∼ 1/η 4 and (Ψ (1) ′ v (1)i d ) ,i ∼ 1/η 3 the main contribution to the integral will come from the two remaining terms, Φ (2) and v (1)2 d , whose amplitudes scale at late times as 1/η 2 Two are the integrals that we have to compute and the one multiplying ln(kη) Let us first consider the contributions from ∇ 2 v (1)2 d . The second integral is easily computed using the expression (311) for the linear CDM velocity. We find that at late times the dominant term is (for kη ≫ 1) The first integral can be computed by making the following approximation: we split the integral between 0 < kη < 1 and kη > 1 and for 0 < kη < 1 we use the asymptotic expression while for kη > 1 we use the limit The the integral for 0 < kη < 1 just gives a constant, while the integral for kη > 1 brings the dominant contribution at late times being proportional to [ln(kη)] 2 so that we can write (kη ≫ 1) As far as the contribution to the integrals (329) and (330) due to ∇ 2 Φ (2) is concerned we have just to keep track of the initial condition provided by the primordial non-Gaussianity. We have verified that all the other terms give a negligible contribution. This is easy to understand: the integrand function on large scale is a constant while at late times it oscillates with decreasing amplitudes as η −2 , and thus the integrals will tend asymptotically to a constant. We find that where γ = 0.577... is the Euler constant, and Φ (2) (0) is given by Eq. (265). Therefore, from Eqs. (334), (331), and (335)-(336) we find that for kη ≫ 1 Let us collect the results of Eqs. (322), (327) and (337) into Eq. (317). We find that for kη ≫ 1 Notice that in Eq. (317) we have neglected Ψ (2) , which decays on subhorizon scales during the radiation dominated epoch (see Eq. (247), and we have used Eqs. (306) and (308). Eq. (338) represents the secondorder Meszaros effect: the CDM density contrast on small scales (inside the horizon) slowly grows starting from the initial conditions that, at second-order, are set by the primordial non-Gaussianity parameter a NL . As one could have guessed the primordial non-Gaussianity is just transferred linearly. The other terms scale in time as a logarithm squared. We stress that the computation of these terms allows one to derive the full transfer function for the matter perturbations at second order accounting for the dominant second-order corrections. In the next section we will use (338) to fix the initial conditions for the evolution on subhorizon scales of the photons density fluctuations ∆ (2) 00 after the equality epoch.

Second-order CMB anisotropies for modes crossing the horizon during the radiaton epoch
In this section we derive the energy density perturbations ∆ (2) 00 of the photons during the matter dominated epoch, for the modes that cross the horizon before equality. In Sec. 7.6 we have already solved the problem assuming matter domination for modes crossing the horizon after equality. Thus it is sufficient to take the result (282) and replace the initial conditions ∆ (2) 00 where we have restored the generic integration constants A and B, Ψ (1) is the linear gravitational potential (which is constant for the matter era) and Ψ (2) m (0) represents the initial condition for the second-order gravitational potential fixed at some time η i > η eq . Eq. (338) allows to fix the proper initial conditions for the gravitational potentials on subhorizon scales (accounting for the fact that around the equality epoch they are mainly determined by the CDM density perturbations). At linear order this is achieved by solving the equation for δ On small scales one neglects the time derivatives of the gravitational potential in Eqs. (302) and (340) to obtain δ where we have also dropped the contribution to the gravitational potential from the radiation component. The solution of this equation is matched to the value that δ (1) d has during the radiation dominated epoch on subhorizon scales, Eq. (308), and one finds that for η ≫ η eq on subhorizon scales the gravitational potential remains constant with We skip the details of the derivation of Eq. (342) since it is a standard computation that the reader can find, for example, in Refs. [42]. Since around η eq the dark matter begins to dominate, an approximation to the result (342) can be simply achieved by requiring that during matter domination the gravitational potential remains constant to a value determined by the density contrast (308) at the equality epoch from Eq. (340) on small scales, leading to where we used a(η) ∝ η 2 during matter domination and Eq. (308) with A 1 = −9.6 and B 1 = 0.44. At second-order we follow a similar approximation. The general solution for the evolution of the the second-order gravitational potential Ψ (2) for η > η eq is given by Eq. (B.6). We have to determine the initial conditions for those modes that cross the horizon during the radiation epoch.

Secondary effects and contamination to primordial NG
Given that a detection of a sizable primordial non-Gaussianity (and its shape) would represent a real breakthrough into the understanding of the dynamics of the universe during its very first stages, it is crucial that all sources of contamination to the primordial signal are well understood and kept under control. In fact any non-linearities can make initially Gaussian perturbations non-Gaussian. Such nonprimordial effects can thus complicate the extraction of the primordial non-Gaussianity: we have to be sure we are not ascribing a primordial origin to a signal that is extracted from the CMB (or LSS) data using estimators of non-Gaussianity when that signal has a different origin. Moreover, as stressed in Section. 7.1, we must always specify of which primordial non-Gaussianity we study the contamination from the non-primordial sources (e.g. primordial non-Gaussianity of "local" , "equilateral" or "folded" shape). Broadly speaking, non-primordial sources of non-Gaussianity can be classified into four categories: instrumental systematic effects; residual foregrounds and unresolved point sources; some well known secondary CMB anisotropies, such as Sunyaev-Zel'dovich (SZ) effect, gravitational lensing, Rees-Sciama effect, and finally previously unknown effects coming from non-linearities in the Boltzmann equations, which are related to the non-linear nature of General Relativity and to the non-linear dynamics of the photon-baryon system. This Review focuses on the last category, but in this Section we offer also a summary of some recent results of the other types of possible contaminations. Before doing that, however, it is important to be more precise about how secondary effects may impact the extraction of a primordial non-Gaussian signal. Mainly they act in two ways. They might "mimic" a three-point correlation function similar in shape to the primordial one. This would produce a bias or a contamination to the estimator of primordial non-Gaussianity. On the other hand secondary effects might increase the variance of the estimator without contributing to the signal-to-noise ratio. In other words, in this case they degrade the signal-to-noise ratio. Let us define in a quantitative way how to characterize these effects.

Signal-to-noise ratio, shapes, and contamination to primordial non-Gaussianity
The angular bispectrum, B ℓ 1 ℓ 2 ℓ 3 , the harmonic transform of the angular three-point function, of the CMB anisotropies is often used to measure non-Gaussianity (see for example [5]). From the usual spherical harmonic coefficients of the temperature anisotropies the angle-averaged bispectrum is given by (for more details see, e.g., Ref. [35]) where the brackets indicate ensemble average. We shall quantify the degree to which the primordial and the secondary bispectra are correlated, as well as their expected signal-to-noise ratio following a method that as now become standard [35,33,65]. Namely, in the limit of weak non-Gaussianity, one can introduce the Fisher matrix for the amplitudes of the bispectra, F ij , given by where the variance of the bispectrum is and ∆ ℓ 1 ℓ 2 ℓ 3 takes values 1, 2, and 6 when all ℓ's are different, two of them are equal and all are the same, respectively. The power spectrum, C ℓ , is the sum of the theoretical CMB and the detector noise.
The signal-to-noise ratio is given by and we define the cross-correlation coefficient between different shapes i and j, r ij , as which does not depend on the amplitudes of the different bispectra, but on the shapes, thus measuring how similar are two bispectra i and j. One can also define a degradation parameter d i = F ii F −1 ii for a degraded (S/N ) due to r ij .
If a secondary bispectra has some correlation with the primordial one, one expects that the signal-tonoise ratio will be degraded: if one does not account for any secondary bispectra (S/N ) 0 = F prim,prim while, marginalizing over a (single) secondary bispectrum, (S/N ) prim gets modified from its zero-order value to (S/N ) prim = (S/N ) 0 (1−r 2 ij ) 1/2 . This means that the minimum detectable value of the primordial non-linearity parameter f NL which is obtained by imposing (S/N ) prim = 1 (or the 1-σ uncertainty on f NL given by δf NL = (F ) −1 prim,prim | f NL =1 ) gets shifted by a quantity ∆f NL /(f NL ) 0 = (1 − r 2 ij ) −1/2 − 1. How much does a given secondary bispectra contaminate the extraction of the primordial bispectrum? If, for example, the predicted shape of the secondary bispectra is sufficiently different from that of the primordial bispectrum, then one would hope that the contamination would be minimal. We can quantify the contamination of the primordial bispectrum as follows: we fit the primordial bispectrum template to the second-order bispectrum, and find the best-fitting f con N L ("con" stands for contamination) by minimizing the χ 2 given by with respect to f NL . Here, B prim ℓ 1 ℓ 2 ℓ 3 is the primordial bispectrum with amplitude f NL = 1 [35]. We obtain This is the value of f NL one would find, if one did not know that the primordial bispectrum did not exist but there was only the secondary bispectrum. In other words the effective non-linearity parameter f con NL gives the value that a bispectrum estimator designed for constraining primordial non-Gaussianity would measure due to the presence of the secondary bispectrum B 2nd ℓ 1 ℓ 2 ℓ 3 . As one can see f con NL turns out to be proportional to the mixed entry of the Fisher matrix, since this governs the correlation between different types of bispectra on the basis of their shape dependence.

Some "well-known" secondary sources of non-Gaussianity
Instrumental systematic effects are one of the principal issues in the search for primordial non-Gaussianity. Since they depend on the specific instrument, here we do not intend to go into any detail, but we refer the reader to Ref. [60] as an interesting example of this kind of effects for the currently flying mission of the satellite Planck. Residual foregrounds and unresolved point sources represent other important possible contaminants and particular attention have been devoted to them in the analysis of the WMAP data [59,36,48]. For an experiment like WMAP secondary effects such as Sunyaev-Zel'dovich effect, gravitational lensing, Integrated-Sachs-Wolfe (or Rees-Sciama) effect should be not so relevant for ℓ < 500 with a bias of just 1.5 to 2 [70] for local primordial non-Gaussianity (see also [35] and [36]). 5 However for higher resolution and more sensitive experiments like Planck, such effects must be taken into account with care. Here we report the results of some recent analyses about the contamination to primordial non-Gaussianity from some secondary effects. Such a brief summary is by no means exhaustive. An interesting case is given by the cross-correlation of the ISW and the lensing effects, see Ref. [33]. In the detailed analysis of [71] it is shown that the ISW-lensing bispectrum produces a bias of f cont NL ≃ 9 to local primordial non-Gaussianity (see also [70]), while the bias to equilateral primordial non-Gaussianity is negligible (see also [36,70]). The impact is mainly on the local type of primordial non-Gaussianity because the ISW-lensing correlation correlates the large-scale gravitational potential fluctuations sourcing the ISW effect with the small scale lensing effects of the CMB, thus producing a bispectrum which peaks in the squeezed configurations as the local shape. As far as the increase of the variance of the estimator of primordial non-Gaussianity induced by CMB lensing is concerned it turns out to be of the order of 20% for an experiment like Planck (see Refs. [71,80]). Notice that, interestingly enough, gravitational lensing can also have another peculiar effect on the extraction of primordial non-Gaussianity: it can modify its shape by smoothing its acoustic features. However, according to the analysis of Ref. [71] (and contrary to the results of [61]) such a "distorsion" of the shape should be small, of the order of 10% for l < 2000.
A similar analysis has been performed looking also at scales where non-linear effects can become relevant. In this case a bispectrum from the correlation of the Rees-Sciama and lensing effects has been studied in Ref. [63], and it has been found that in this case the bias to a local primordial non-Gaussianity corresponds to f cont NL ≃ 10, agreeing with the results of [71] where the two analyses can be compared. Another secondary source of contamination that has been studied quite recently in Ref. [62] deals with bispectra generated by correlations of number density and lensing magnification of radio and SZ point sources with the ISW effect. Also in this case the large-scale modulation of small-scale number density due to fluctuations which source the ISW generates a bispectrum which peaks on squeezed configurations and it has been found that it corresponds to a contamination of local primordial non-Gaussianity of f cont NL ≃ 1.5, which should be not so relevant for an experiment like Planck. The novelty of this effect is that it is different from the usual bispectrum of point sources when treated with a Poisson distribution (see, e.g., [35]).
Other secondary effects include, for example, the SZ-SZ-SZ bispectrum [64] and non-Gaussianities from the kinetic SZ and Ostriker-Vishniac effect [66,65]. However the SZ-SZ-SZ bispectrum can be considered as an extra Poisson contribution to the unresolved point sources for l < 1500 [70]. Moreover recently an analysis of the bispectrum generated from inhomogeneous reionization has been performed in [67], where it has been shown that it can give rise to a very small contamination to the local primordial non-Gaussianity of f cont NL ≃ −0.1. Finally let us also recall some recent studies of a specific effect arising at second-order in perturbation theory, when looking at the fluctuations in the Boltzmann equations for the photon-baryon fluid. In Refs. [72,68,73,74] the perturbations to the phase of recombination between electrons and baryons has been considered ("inhomogeneous recombination") . It turns out that the electron density perturbation can be of the order of 5 times bigger than the one in the density of the baryons. This gives rise to the prospect of a non-Gaussianity which corresponds to |f cont NL | ≃ 5. In fact the detailed study of this effect shows that the contamination is much smaller. The bispectrum from the perturbed recombination peaks in the squeezed configuration with a contamination to the primordial non-Gaussianity of the local type of |f cont NL | ≃ 0.7. Let us conclude this section by mentioning that the non-linearities emerging from secondary effects can be very interesting by themeselves, since they carry a lot of information about the evolution of the universe and the growth of structures after recombination till very low redshift. Detailed examples and reviews about this aspect can be found, e.g., in Refs. [75,76,77]. 9 How to analytically estimate the signal-to-noise ratio from the NG at recombination and its contamination to the primordial NG As we have seen in Sec. 7, the dynamics at recombination is quite involved because all the non-linearities in the evolution of the baryon-photon fluid at recombination and the ones coming from general relativity should be accounted for. Such a contribution is so relevant because it represents a major part of the second-order radiation transfer function which must be determined in order to have a complete control of both the primordial and non-primordial part of NG in the CMB anisotropies. The NG generated at the surface of last scattering comprises various effects. In this Section, we devote our attention to one particular relevant contribution, the one coming from the non-linear evolution of the second-order gravitational potential which grows in time on small scales. The analysis of this contribution offers the opportunity to understand how some secondary effects can contaminate mostly a given type of primordial NG, while leaving almost untouched other different forms of primordial NG. In this case this effect is a causal one, developing on small scales, so we expect that the NG it generates will be mainly of the equilateral type, rather than of local type. Therefore, a reasonable question is to which extent the NG from recombination alters the possibile detection of the primordial NG of the equilateral type. The goal of this Section is therefore to illustrate how to estimate in a semi-analytical way the contribution to NG from recombination.

Signal-to-Noise ratio for the primordial equilateral bispectrum
In this subsection we wish to recover the estimate for the signal-to-noise ratio (S/N ) given in Ref. [36] for the primordial bispectra of "equilateral" type [78] by adopting a simple model. In other words, we test the goodness of the semi-analytical model we will be using in the next Section to estimate the bispectrum from the recombination era. Our starting point is the primordial equilateral bispectrum [78] Φ where and the permutations act only on the last term in parentheses. The parameter f equil NL quantifies the level of NG while A = 17.46 × 10 −9 is the amplitude of the primordial gravitional potential power spectrum computed at first-order with P (k) = A/k 3 . Since the signal-to-noise ratio (S/N ) will be some function of the maximum multipole a given experiment can reach, ℓ max ≫ 1, we can use the flat-sky approximation [79,80] and write for the bispectrum where where k ′ means k evaluated such that k = ℓ/(η 0 − η r ) and is the radiation transfer function defined by the CMB source function S(k, η). In this notation, η 0 and η r represent the present-day and the recombination conformal time, respectively and k z and k are the momentum components perpendicular and parallel respectively to the plane orthogonal to the line-ofsight. The (S/N ) ratio in the flat-sky formalism is [79,80] S N where f sky stands for the portion of the observed sky. In order to compute the bispectrum B equil (ℓ 1 , ℓ 2 , ℓ 3 ) and the power spectrum C(ℓ) we adopt the following model where we mimic the effects of the transfer function on small scales as i.e. a simple exponential and a normalization coefficent a to be determined to match the amplitude of the angular power spectrum at the characteristic scale ℓ ≃ ℓ * = k * (η 0 − η r ). 6 It is important to make clear what are the reasons underlying the choice of such a model. When computing the (S/N ), Eq. (365) with ℓ * = k * (η 0 − η r ) ≃ 750 and a ≃ 3 is able to account for the combined effects of "radiation driving", which occours at ℓ > ℓ eq ≃ 160 and boosts the angular power spectrum with respect to the Sachs-Wolfe plateau, and the effects of Silk damping which tend to suppress the CMB anisotropies for scales ℓ > ℓ D ≃ 1300.
The combination of these effects produces a decrease in the angular power spectrum from a scale ℓ * ≃ 750. 7 The power spectrum in the flat-sky approximation is given by a( l 1 )a( l 2 ) = (2π) 2 δ (2) ( l 12 )C(ℓ 1 ) with The exponential of the transfer function for Eq. (365) allows to cut off the integral for k ≃ k * and one finds (see also Ref. [80]) where the last equality holds for ℓ ≫ ℓ * . To compute the bispectrum we proceed in a similar way. One first uses the Dirac deltas, δ (1) (k z 123 ) and δ (2) ( ℓ 123 ). Then it proves to be useful the change of variable k z 1 = In this way the transfer functions become∆ T (ℓ i , k z i ) ∝ e −1/2(|x i |ℓ i /ℓ * ) 1.2 which allows to cut the integrals over x i (i = 1, 2) at ℓ * /ℓ i . Now, as a good approximation to see the effects of the transfer functions, we can take ℓ ≫ ℓ * and thus the integral over x i can be easily computed by just evaluating the integrand in x i = 0 times 4(ℓ * /ℓ 1 )(ℓ * /ℓ 2 ). With this approximation the integral in k z i is easily obtained and we get for the bispectrum where The coefficient f 1 ≃ 1/1.4 = 0.7 is a fudge factor that improves the matching between our approximation for the bispectrum and numerical results that have been consistenly checked. Notice that, according to our approximation, the equilateral structure of Eq. (359) is preserved in ℓ space. 8 ℓ 1 = ℓ 2 , in correspondence of the maximum following expression In computing the signal-to-noise ratio, consistency with our approximation (369) requires that we integrate over ℓ 1 , ℓ 2 starting from a minimum ℓ min > ℓ * up to ℓ max and paying attention to the fact that even ℓ 3 in Eq. (370) must be larger than ℓ min . The scaling with ℓ max with respect to the case of a local type bispectrum turns out to be much milder. While for the local type (S/N ) 2 ∝ ℓ 2 max [80], for the equilateral bispectrum (358) we find 9 (S/N ) 2 ∝ ℓ max and, setting 7 The choice of the exponent 1.2 derives from the study of the diffusion damping envelope in Ref. [81]. 8 The expression (369) can be also written as , θ being the angle between ℓ 1 and ℓ 2 .

Non-Gaussianity from recombination
Comforted by the goodness of our model, in this Section we wish to estimate the level of NG generated at the recombination era. One can check that on small scales the second-order anisotropies are dominated by the second-order gravitational potential Φ (2) which grows as η 2 , as first pointed out in Ref. [82]. Let us first briefly show how this result can be obtained using the calculations of Sec. 7. Consider the tight coupling regime, and small scales, well inside the diffusion Silk-damping scale λ D . Then look at Eq. (229). One realizes that, apart from the combination ∆ (2) 00 + 4Φ (2) , all the other terms scale at most as η or they are suppressed either because in the tight coupling regime or because of the damping diffusion on small scales. On the other hand one finds that the main contribution to the bispectrum generated at recombination comes from which is the usual term appearing in the CMB anisotropies due to the intrinsic photon energy density fluctuations ∆ (2) 00 and the gravitational redshift due to the potential. In fact such a term on large scales reduces to the Sachs-Wolfe effect, while on small scales at recombination, using Eq. (279) and (282), we find 10 We thank M. Liguori for discussions about the minimum value of f equil NL detectable by Planck and for its scaling with ℓ max .
where we have evaluated the expression at the recombination time η r and R = 3ρ b /4ρ γ is the baryon-tophoton energy density ratio. In writing Eq. (374) we have dropped off the contribution from primordial non-Gaussianity (which, anyway, always propagates linerarly) since here we are focusing on the relevant secondary effects that can constitute a source of contamintation to it (see Eq.(199) and the discussion below). Notice also that this expression has been obtained assuming all the momenta much larger than k eq . Eq. (374) has been confirmed numerically in Ref. [82]. There is also another way to obtain the analytical expression in the first line of Eq. (374), as discussed in details in Ref. [82], which exploits some well-known results at linear order extending them at second-order. At second-order in the perturbations the solution for the acoustic oscillations will have a form very similar to the linear solution written in Eq.(C.18), except, as usual, of some source terms S made by products of fistr-order perturbations. ∆ (2) 00

= ∆
(2) 00 Now on sufficiently small scales the products of first-order terms are indeed suppressed, due either to Silk damping, or because gravitational potentials are suppressed as (k/k equ ) 2 . Recall that only acoustic oscillations are damped by Silk damping diffusion, so that the cosine in Eq.(375) is multiplied by e −(k/k D ) 2 k ≫ k D , where k −1 D = λ D indicates damping length, but the baryon induced shift (1 + R)Φ (1) is left untouched. Therefore on small scales the combination of the damping effects and the growth of the secondorder gravitational potential Φ (2) as η 2 single out the dominant effect in Eq.(375) to be Θ

Bispectrum from recombination due to the non-linear growth of the gravitational potential
Let us start to go into some details of the form of the bispectrum induced by the contribution (374). In Eq. (374) the kernel G is given by From the form of the kernel we see that the NG at recombination is dominated by an equilateral configuration, as expected from the fact that its origin is gravitational. Here and in the following we are implicitly assuming that a convolution is acting on the kernel as The reader should remember that, at first-order in perturbation theory, the combination Θ is exponentially suppressed by the Silk damping, but still greater than the term RΦ (1) (which does not suffer the damping) for the maximum multipole of interest, ℓ max ∼ 2000. This is meanly due to the fact that the first-order gravitational potential rapidly decays on small scales. On the contrary, at secondorder in perturbation theory, the gravitational potential grows like the scale factor on small scales and it turns out that the RΦ (2) dominates on small scales (see Ref. [82]).
The gravitational potential at linear order can be expressed as usual in terms of the transfer function T (k) where the last step is an approximation valid on scales smaller the the equivalence scale, k ≫ k eq . In the following we will account for the logarithmic growth just with a coeffcient T 0 (k) = 12 ln[k/8k eq ] ≈ 11 for the scales of interest.
In the flat-sky approximation one arrives at an expression similar to (362), where now one of the linear transfer functions must replaced by a transfer function at second-order. Specifically one finds By using our model (365) and∆ for the second-order radiation transfer function, we find At this point we proceed further by employing the same approximation described after Eq. (368). We use the Dirac delta to replace the variable k 3z , and the exponential allow us to evaluate the integral for k 1z = k 2z = 0, for scales ℓ i ≫ ℓ * . This leads to A 2 a 2 T 2 0 (k eq η r ) 2 ℓ 2 eq ℓ 2 * e −1/2(ℓ 1 /ℓ * ) 1.2 e −1/2(ℓ 2 /ℓ * ) 1.2 Again here f 2 is a coefficient to better calibrate our approximations with numerical results that we have performed in order to test the validity of our approach. Not surprisingly, it turns out that f 2 ≃ f 1 ≃ 1/1.4.

Contamination to primordial non-Gaussianity from recombination: Fisher matrices
Our goal now is to quantify the level of NG coming from the recombination era and to estimate the level of degradation it causes on the possible measurement of the equilateral primordial bispectrum. The reader should keep in mind that, given the form of the kernel function (376), the NG from recombination is expected to be of the equilateral type. Following Sec. 8.1, a rigorous procedure is to define the Fisher matrix (in flat-sky approximation) as where i (or j)= (rec, equil), and to define the signal-to-noise ratio for a component i, (S/N ) i = 1/ F −1 ii , and the degradation parameter d i = F ii F −1 ii due to the correlation bewteen the different components The first entry F equil,equil of the Fisher matrix corresponds to the (S/N ) 2 ratio computed in Eq. (371) which does not account for any kind of cross-correlation. Due to the equilateral form of the NG generated at recombination we expect that the minimum value detectable for f equil NL will be higher that the one reported in Eq. (372). For the mixed entry we find where ℓ 3 is given by Eq. (370). The factor 3 in front of this expression comes from cyclic permutations. The integral can be performed numerically and, integrating from a minimum ℓ min ≃ 1200 up to ℓ max = 2000, and by taking R ≃ 0.3 when evaluated at recombination, a ≃ 3, T 0 ≃ 11, (k eq η r ) 2 ≃ 26, ℓ eq = 150, ℓ * = 750, we find F rec,equil ≃ 9.4 × 10 −4 . Finally for the entry F rec,rec we get and we find a value F rec,rec ≃ 0.014. We are now able to compute the entries of inverse of the Fisher matrix, F −1 ij . In the following we report our results for the signal-to-noise ratios and the degradation parameters As a confirmation of our expectations, we find that the NG of the type given by Eq. (376) has a quite high correlation with an equilateral primordial bispectrum. This translates into a degradation in the mimimum detectable value for f equil NL with respect to the value given in (372). In fact from the signal-tonoise ratio (385) we find a minimum detectable value of imposing that (S/N ) equil = 1 with an an increase of ∆f equil NL = O(10). This corresponds to an increase of the 1-σ uncertainty on f equil NL of 20% (see Sec. 8.1). 11 Following Sec. 8.1, in order to measure the contamination to the primordial bispectra one can define that effective non-linearity parameter f con NL which minimizes the χ 2 defined as to find and an analogous expression to compute the contamination to the local primordial bispectrum. In this case we find an effective non-linearity parameter f cont NL ≃ 5. 12 Similarly we can compute the Fisher matrix accounting for the NG generated at recombination and the primordial NG of the local type where (394) 11 Due to a non-vanishing correlation, r ij , (S/N ) gets modified from its zero-order value to (S/N ) = (S/N ) 0 (1 − r 2 ij ) 1/2 , so that the minimum detectable value of f equil NL gets shifted by a quantity ∆f equil NL /(f equil NL ) 0 = (1−r 2 ij ) −1/2 −1. 12 Notice that sometimes in the literature one can find also an effective non-linearity parameter f eff.
NL defined in such a way that the equilateral (or the local) bispectrum has the same Fisher matrix errors as the recombination bispectrum However this is not the proper quantity to define the contamination level to primordial NG; here we are just comparing signal-to-noise ratios, while the contamination f con NL defined in Eq. (400) contains a somewhat richer information, since we are asking what is the value of equilateral (local) f NL which best mimics the bispectrum from recombination.
The bispectrum and the signal-to-noise ratio as defined in Eq. (364) have already been computed in the flat-sky approximation in Ref. [80]. The result is that (S/N ) 2 loc = 4π −2 f sky (ℓ * /ℓ min )(f loc NL ) 2 A ℓ 2 max , corresponding to a minimum detectable value of f loc NL = O(7) for ℓ max = 2000 (when other possible sources of NG are ignored). We can compute the off-diagonal entry of the Fisher matrix in a similar way to what we have described in this section, and we get F rec,loc ≃ 8 × 10 −3 f loc NL . Finally the entry F rec,rec ≃ 0.014 has already been computed above. From inverting the Fisher matrix, we get the following signal-to-noise ratios and the degradation parameters In particular, from Eq. (395) we see that now the minimum detectable value of f loc NL remains basically unchanged in the presence of the recombination signal. Similarly the effective f con NL reads which is much smaller than the effective non-linearity parameter (400) for the equilateral case. We have also checked the cross-correlation between the primordial local and equilateral bispectra finding a value of r loc,equil ≃ 0.23, which is in agreement with the value reported in Ref. [36]. This reflects the fact that the primordial local and equilateral signals are not fully uncorrelated. The reason is due to the fact that the equilateral and local bispectrum (394) and (359) approach the same shape in the equilateral configuration. This is also the reason why the cross-correlation between the primordial local and recombination bispectra is not so small.
10 How to perform a numerical analysis of the CMB bispectrum produced by second-order perturbations In the previous section, through a specific example, we have learned some basic principles to determine a reasonable and quite fast analytic estimate for the contamination to the primordial NG. In this section we provide all the tools necessary to a full numerical implementation of the second-order Boltzmann equations which allow to obtain a systematic numerical evaluation of the CMB angular bispectrum produced by second-order cosmological perturbations. This section is mainly based on the results of Ref. [84].
In particular we will apply this formalism to numerically evaluate the contamination to primordial NG from the second-order fluctuations of the Boltzmann equations that come as products of the first-order perturbations, and ignore the intrinsically second-order terms, or the effects of the perturbed recombination [72,73,74]. 13 Here therefore we come to the discussion of the second example mentioned in Sec. 7.1. As anticipated there we will see that, unlike the intrisically second-order contribution considered in the previous section, the products of the first-order perturbations give a CMB bispectrum that peaks in the squeezed configuration. Therefore we shall study the contamination of the primordial NG of the local type.

.1.1 Definitions
Again, we expand the temperature fluctuation into the linear (first-order) part and the second-order part as ∆T (n) The spherical harmonic coefficients of temperature anisotropy, a ℓm = T −1 d 2n Y * ℓm (n)∆T (n), are therefore expanded as Recall that to expand the Boltzmann equation up to the second order in perturbations, we had to expand the distribution function, up to the second order in perturbations: Θ = Θ (1) + Θ (2) /2 + . . . , and accordingly f = f (0) + f (1) + f (2) /2 + . . .. We then had to compute the fractional perturbation in photon's energy density at the i-th order in perturbations, ∆ (i) , by multiplying f (i) by p, and integrating over p 2 dp: At the linear order, we recovered the usual relation between the linear fractional temperature fluctuation, Θ (1) = ∆T (1) /T , and the linear fractional energy density perturbation, ∆ (1) = δρ (1) γ /ρ γ , i.e., ∆ (1) = 4Θ (1) . At the second order we have which is related to the second-order temperature fluctuation as ∆T T (2) where we have subtracted the average of the temperature fluctuation so that the average of ∆T (2) /T vanishes. We compute a (2) ℓm from ∆T (2) /T using where we defineã Here the matrix is the Wigner 3j symbol. The CMB angular-averaged bispectrum, B ℓ 1 ℓ 2 ℓ 3 , is related to the ensemble average of a ℓ 1 m 1 a ℓ 2 m 2 a ℓ 3 m 3 as This definition guarantees rotational invariance for the bispectrum, and the Wigner 3j symbol ensures that the bispectrum must satisfy triangle conditions: ℓ i − ℓ j | ≤ ℓ k ≤ ℓ i + ℓ j for all permutations of indices, and selection rules: m 1 + m 2 + m 3 = 0. The ensemble average is given by where cyclic means that we have to sum the cyclic permutations of Eq. (410) for indices (1, 2, 3) → (3, 1, 2) → (2, 3, 1). As we assume that a ℓm 's in Eq. (410) is given by the sum of products of all possible pairs. Each pair gives the angular power spectrum, C ℓ : We obtain a (1) Substituting the right hand side of equation (412) for the second term of equation (410), and using ℓ 1 + ℓ 2 + ℓ 3 = even, we obtain the angular averaged bispectrum, where we have defined the quantities, andB

The second-order CMB radiation transfer function
In this section we show how it is possible to define in a rigorous way the radiation transfer function for CMB anisotropies at second-order in the perturbations. It is a generalization of the well known quantity used at linear-order and as such it allows us to develop a systematic numerical analysis of the angular bispectrum produced by second-order perturbations in the very same way as at linear-order various numerical codes are nowadays available for the computation of the CMB power spectrum, such as CMBFAST or CAMB. The starting point is the Boltzmann equation that governs the evolution of ∆ (1) (k, µ, η) and ∆ (2) (k,n, η), where µ =k ·n and n is the direction of propagation of photons. Note that for the linear perturbation there is azimuthal symmetry such that ∆ (1) depends only on the angle between k and n; however, for the second-order perturbation there is no such symmetry. We write again the Boltzmann equations in Fourier space where the primes denote derivatives with respect to the conformal time ∂/∂η, S (1) and S (2) are the source functions at the first and the second orders, respectively, and τ ′ is the differential optical depth which is defined by using the mean electron number density,n e , the Thomson scattering cross-section, σ T , and the scale factor, a, as We expand again the angular dependence of ∆ (i) as and that of the source terms as where i = 1, 2.
The source functions relate the observed a ℓm 's to the primordial curvature perturbations in comoving gauge, ζ(k). The relations contain the linear radiation transfer function, g ℓ (k), and the second-order radiation transfer function, F ℓ ′ m ′ ℓm (k), and are given by The linear transfer function is given by where u ≡ k(η 0 − η) and S (1) ℓm is the standard linear source function where Φ (1) (k, η) and Ψ (1) (k, η) are the metric perturbations at the linear order in the longitudinal gauge and ∆ 1 (k, η), and ∆ 2 (k, η) are the coefficients of the expansion in Legendre polynomials of ∆ (1) (k, µ, η). These coefficients ∆ The first-order velocity perturbation, v 0 (k, η), is the irrotational part of the baryon velocity defined by v(k) = −iv 0 (k)k.
The new piece, the second-order radiation transfer function, is the line-of-sight integral of the secondorder source terms in the Boltzmann equation: Here, we have introduced a new function, S ℓm (k ′ , k ′′ , k, η), which is defined by the following equation: Basically, S ℓm (k ′ , k ′′ , k, η) is the second-order source function divided by ζ(k ′ )ζ(k ′′ ). The explicit expression of the second-order source function can be be read from Eq. (134).
Using equation (421) and (422), we calculate the first term in Eq. (413),B ℓ 1 ℓ 2 ℓ 3 , as follows: where, as in the previous section we will sometimes use the notation k 123 = k 1 + k 2 + k 3 . P ζ (k) is the power spectrum of ζ given by the usual definition: In order to perform the integral over angles,k, we expand the three-dimensional δ-function using Rayleigh's formula, and also expand the angular dependence of S ℓm (k 1 , k 2 , k 3 , η) by introducing the transformed source function, S µ 1 µ 2 µ 3 λ 1 λ 2 λ 3 (k 1 , k 2 , k 3 , η), as Finally, by adding the remaining term in the full bispectrum, Eq. (413), we obtain The remaining task is to calculate the angular-averaged source term, S λ 1 λ 2 λ 3 (k 1 , k 2 , k 3 , r) according to Eq. (438). Up to now our scope has been to offer the reader all the general tools that allow a systematic treatment of the CMB angular bispectrum from second-order perturbations. In the following we will give just some examples of how such tools can be implemented in order to obtain numerical results about signal-to-noise ratios and contamination to primordial non-Gaussianity.

Source Term: an example
The explicit expression for S (2) ℓm (k, η) can be obtained from Eq. (134). Here we do not want to report the complete expression S (2) ℓm (k, η) of the source function at second-order, rather, for the goal of this Review, we think it is more instructive to show explicitly how the calculation of S λ 1 λ 2 λ 3 proceeds focusing on just one simple example. Therefore consider, for example, the term of the second-order source term from Eq. (134) given by First we compute the multipole coefficients of the source term S (2) ℓm (k) from this contribution, as defined in Eq. (420). They are given by the convolution Now for ∆ (1) ℓm (k) we use so that We now compute the corresponding "angular-averaged source function" coefficients S λ 1 λ 2 λ 3 (k 1 , k 2 , r) defined by Eq. (438). From Eq. (445) you read the kernel defined in Eq. (429) From here S (2) λ 1 λ 2 λ 3 (k 1 , k 2 , r) = i λ 1 +λ 2 2λ 1 + 1 4π where we have used Y * ℓm = (−1) −m Y ℓ−m , and the orthonormality of the spherical harmonics. Using the property of the Wigner 3-j symbols we find S This term therefore corresponds to coefficients S 0λ 2 λ 2 (k 1 , k 2 , r).
In general the perturbation variables of the source term can be split into two parts (see, for example, Eq.(134)). A part containing perturbations that are intrinsically second-order (these perturbations have superscripts (2), and ω m and χ m are also intrinsically second-order). Solving for these terms requires solving the full second-order Boltzmann equations coupled with the Einstein equations.
Another part contains terms that are products of two linear variables, as the example that we have just considered. Evaluation of these terms is much easier than that of the intrinsically second-order terms, as the first-order variables have already been calculated using the standard linearized Boltzmann code such as CMBFAST.
10.2 Second-order bispectrum from products of the first-order terms 10

.2.1 A worked example
We shall now focus only on the products of the first-order perturbations and we choose to analyze just some of these contributions to offer the reader an example of the analysis one can perform. We warn the reader that the full analysis is under progress and that the full results will be presented in Ref. [85], including the intrinsically second-order perturbations which are equally important, and the contribution from perturbing the recombination history [72,73,74].
For the products of the first-order perturbations, from now on we will consider the following non-zero four cases for the source terms, S λ 1 λ 2 λ 3 , which have been analyzed in Ref. [84] (for notational simplicity we shall omit the superscripts (1)): 14 In particular it is recognizable the contribution to S 000 from the example just discussed. From these results we find that S λ 1 λ 2 λ 3 does not depend on k 3 , i.e., S λ 1 λ 2 λ 3 = S λ 1 λ 2 λ 3 (k 1 , k 2 , r). Note also that S 011 (k 1 , k 2 , r) = S 101 (k 2 , k 1 , r).

Bispectrum from products of the first-order terms
For the four non-vanishing combinations of λ 1 , λ 2 , and λ 3 in Eq. (450), we rewrite the expression for the bispectrum, Eq. (441), as (1,0,1) where we have used B (0,1,1) (1,0,1) and To proceed further, it turns out to be useful to simplify the expression by introducing the following notation for the integral over k that appears many times: This function corresponds to the existing functions in the literature in the appropriate limits. For example, for x(k, r) = π/2, this function is the same as β (n) ℓℓ ′ (r) introduced in [86]. In fact, we find that an order-of-magnitude estimate of [x] ℓℓ ′ (r)/π × x(k = ℓ ′ /r, r) for a smooth function of x(k, r). As β (n) ℓℓ ′ (r) is a sharply peaked function at the decoupling epoch, r = r * , we find that [x] (n) ℓℓ ′ (r) is also sharply peaked at r = r * . With these tools in hand, one can calculate B (0,0,0) (1,0,1) ℓ 1 ℓ 2 ℓ 3 , and B (1,1,2) ℓ 1 ℓ 2 ℓ 3 . We do not bother here the reader with the details of this computation whose details can be found in Ref. [84]. Rather here we prefer to go straight to the results of the analysis of these bispectra which can be particularly instructive as far as their shape and their contamination to primordial non-Gaussianity are concerned. 10.3 Shape and signal-to-noise of the second-order bispectrum from products of the first-order terms One of the motivations for calculating the second-order bispectrum is to see how much the second-order effects in gravity and the photon-baryon fluid contaminate the extraction of the primordial bispectrum. If, for example, the predicted shape of the second-order bispectrum is sufficiently different from that of the primordial bispectrum, then one would hope that the contamination would be minimal. To investigate this, we shall compare the numerical results of the second-order bispectrum with the so-called "local" model of the primordial bispectrum. We extract the first-order perturbations from the CMBFAST code. We use the following cosmological parameters: Ω Λ = 0.72, Ω m = 0.23, Ω b = 0.046, h = 0.70, and assume a power law spectrum, P ζ ∝ k n−4 , with n = 1. We determine the decoupling time, η * , from the peak of the visibility function. In this model we have cη 0 = 14.9 Gpc and cη * = 288 Mpc. While the most of the signal is generated in the region of the decoupling epoch, in the low-ℓ regime we must also take into account the late time contribution due to the late integrated Sachs-Wolfe effect; thus, we integrate over the line-of-sight, r, in the following regions: c(η 0 − 5η * ) < r < c(η 0 − 0.7η * ) for ℓ > 100, and 0 < r < c(η 0 − 0.7η * ) for l ≤ 100. The step size is ∆r = 0.1η * around the decoupling epoch, and we use the same time steps used by CMBFAST after the decoupling epoch.
The local primordial bispectrum is given by [35] B ℓ 1 ℓ 2 ℓ 3 = 2I ℓ 1 ℓ 2 ℓ 3 ∞ 0 r 2 drb L ℓ 1 (r)b L ℓ 2 (r)b N L ℓ 3 (r) + cyclic , Note that our linear transfer function, g ℓ (k), is related to that of [35], g KS T ℓ (k), by g ℓ (k) = 3 5 g KS T ℓ (k). Figure 1 shows the shape of the bispectrum generated by the products of the first-order terms f selected in Eq. (450), and compares it to the primordial local bispectrum, for ℓ 3 = 200. Both shapes (second-order and primordial) have the largest signals in the squeezed triangles, ℓ 1 ≪ ℓ 2 ≈ ℓ 3 . This is an expected result: both the local primordial bispectrum and the second-order bispectrum that we have computed here arises from the products of the first-order terms, also products in position space. However, these two shapes are slightly different when ℓ 1 /ℓ 3 is not so small (ℓ 1 /ℓ 3 = O(0.1)): the ways in which the radiation transfer function (which gives the acoustic oscillations) enters into the bispectrum are different for the products of the first-order terms and the primordial bispectrum. The primordial bispectrum contains j ℓ (kr * )g ℓ (k), whereas the second-order bispectrum contains j ℓ (kr * )g ℓ (k)x(k, r * ) where x = ∆ 0 , v 0 , etc., also has the oscillations. Therefore, the second-order bispectrum has more interferences between multiple radiation transfer functions. Moreover, the second-order effects contain derivatives that the local primordial effects do not have, which also makes the details of the two shapes different.
Notice, in particular, that most of these gradients in the source term, Eq. (134), are contracted with the direction vector,n. There is only one term that has a scalar product of two wave-vectors, k 1 · k 2 , which vanishes in the squeezed limit. The resulting bispectrum, Eq. (451), resembles that of a local form, except for the extra powers of k coming from the derivatives. These extra powers of k will affect the scale-dependence of the bispectrum, i.e., the second-order bispectrum is no longer scale-invariant. Nevertheless, the largest signal of the bispectrum still comes from the squeezed configurations, as the number of extra powers of k from the derivatives in the source term is not large enough to change the fact that we have the largest contribution when one of k 1 , k 2 , and k 3 is very small. In other words, schematically the bispectrum looks like B(k 1 , k 2 , k 3 ) ∼ (k m 1 1 k m 1 2 )/(k 3 1 k 3 2 ) + cyclic, where m 1 and m 2 are the extra powers of k from the derivatives. Therefore, the largest contribution is in the squeezed configurations as long as m i < 3. Figure 2 shows the same for ℓ 3 = 1000. The results are similar to those for ℓ 3 = 200, but the acoustic oscillations are more clearly visible.
We can quantify the degree to which the second-order and the primordial bispectra are correlated, as well as the expected signal-to-noise ratio of the second-order bispectrum, following the general definitions summarized in Sec. 8. Notice that we are ignoring the noise contribution. In other words, we shall only consider ideal cosmic-variance limited experiments with full sky coverage, which however, at least for the multipole maximum multipole of l max ∼ 2000 we will consider, is a good reference for an experiment like Planck (see, e.g., [35]). products of the first-order primordial (for f NL =1) Figure 3: Signal-to-noise ratios for the local primordial bispectrum for f NL = 1 (dashed), and the second-order bispectrum from the products of the first-order terms (solid), for an ideal full-sky and cosmic-variance-limited (noiseless) experiment.
The signal-to-noise ratio is given by Eq. (354). In Fig. 3 we show the cumulative signal-to-noise ratio, summed up to a maximum multipole of ℓ max , of the primordial bispectrum, assuming f NL = 1 and ignoring the second-order bispectrum, i.e., (S/N ) prim = (F prim,prim ) 1/2 , as well as that of the secondorder bispectrum, ignoring the primordial bispectrum, i.e., (S/N ) 2nd = (F 2nd,2nd ) 1/2 . In both cases S/N increases roughly as S/N ∝ ℓ max (or ∝ N pix where N pix is the number of independent pixels in the map). A larger contribution to the second-order bispectrum at ℓ ≤ 50 comes from the terms involving the Integrated Sachs-Wolfe effect. The signal-to-noise ratio of the second-order bispectrum reaches ∼ 0.4 at ℓ max = 2000; thus, this signal is undetectable. While our calculation includes the temperature anisotropy only, including polarization would increase the signal-to-noise by a factor of two at most, which would not be enough to push the signal-to-noise above unity.
While the total signal-to-noise does not exceed unity, it may still be instructive to show which terms of B (λ 1 ,λ 2 ,λ 3 ) ℓ 1 ℓ 2 ℓ 3 and B C ℓ ℓ 1 ℓ 2 ℓ 3 are more important than the others. To do this we show the following quantity: where a, b = 1, 2, 3, 4, and 0 correspond to (0, 0, 0), (1, 1, 0), (1, 0, 1), (1, 1, 2), and C ℓ , respectively. The results are shown in Fig. 4. We find that (S/N ) 2nd is dominated by B (λ 1 ,λ 2 ,λ 3 ) ℓ 1 ℓ 2 ℓ 3 for ℓ ≤ 100, whereas it is dominated by B C ℓ ℓ 1 ℓ 2 ℓ 3 for ℓ > 100 (see the top panel). Among B (λ 1 ,λ 2 ,λ 3 ) ℓ 1 ℓ 2 ℓ 3 , the most dominant term is (1, 0, 1) (the bispectrum from the second-order dipole created by the first-order dipole and monopole). The second most dominant is (0, 0, 0) (from the second-order monopole created by the  l max Figure 5: The cross-correlation coefficient between the second-order bispectrum from the products of the first-order terms and the local primordial bispectrum. first-order monopole) for l ≤ 400 and (1, 1, 0) (from the second-order monopole created by the first-order dipole) for l > 400. The cross terms (middle and bottom panels) are sub-dominant compared to the auto terms (top panel) at all multipoles. How similar are the second-order and the primordial bispectra? In Fig. 5 we show the cross-correlation coefficient, between the local bispectrum and second-order bispectrum from the products of the first-order terms given in Eq. (450). The cross-correlation coefficient (see Eq. (355), reaches ∼ 0.5 for ℓ max = 200, and the shapes for ℓ 3 = 200 are shown in Fig. 1. After ℓ max = 200 the correlation weakens, and reaches ∼ 0.35 at ℓ max = 1000, and the shapes for ℓ 3 = 1000 are shown in Fig. 2. These results show that the second-order bispectrum from the products of the first-order perturbations and the local primordial bispectrum are fairly similar, with a sizable correlation coefficient.
How large is the contamination of the primordial bispectrum? The level of contamination is measured by the effective non-linearity parameter given by Eq. (357).
In Fig. 6 we show f con NL as a function of the maximum multipoles, ℓ max . We find that f con NL reaches the maximum value, ∼ 0.9, when the correlation coefficient reaches the maximum at ℓ max ∼ 200, but then decreases to ∼ 0.5 at ℓ max ∼ 2000. Therefore, we conclude that the contamination of the primordial bispectrum due to the second-order bispectrum from the terms in Eq. (450) is negligible for CMB experiments.
Finally, one can also calculate the 1-σ uncertainty of f NL , δf NL , with the second-order bispectrum marginalized over. This is given by δf NL = (F ) −1 prim,prim | f NL =1 . Fig. 7 shows that an increase in the uncertainty of f NL due to marginalization is totally negligible. l max Figure 6: Contamination of the local primordial bispectrum as measured by f con NL (Eq (357)).

Conclusions
In these review we have addressed a basic question in cosmology: how a primordial NG propagates into an observable like the CMB anisotropy. Answering this question is fundamental as it will help us in getting some knowledge about the way the primordial cosmological perturbation was generated at the very early stages of the evolution of the universe. In the first sections we have shown how to set the initial conditions at second-order for the (gauge-invariant) CMB anisotropy when some source of primordial NG is present. This was more or less straightforward because on large angular scales is basically gravity which dictates the non-linear dynamics. On small angular scales the computation of the second-order effects in the CMB anisotropy is far more difficult, as so many sources of non-linearities appear. In this review, we have focussed ourselves on the study of the second-order effects appearing at the recombination era when the CMB anisotropy is left imprinted. We have shown how to derive the equations which allow to evaluate CMB anisotropies, by computing the Boltzmann equations describing the evolution of the baryon-photon fluid up to second order. This allows to follow the time evolution of CMB anisotropies (up to second order) on all scales, from the early epoch, when the cosmological perturbations were generated, to the present time, through the recombination era. Through some analytical and simplified example we have also shown how estimate the contamination of the recombination secondary effects onto the detection of primordial NG. More refined numerical results confirm our estimates and are also reported here. It goes without saying that this line of research should be pursued until the level of accuracy in the theoretical prediction is reached and is comparable to the one provided by the current or future satellite experiments. where δ (1) k (0) is the initial condition in the matter-dominated period. At second-order, using Eqs. (35) and (33) and the fact that the first-order gravitational potential is constant, we find and equation for the gravitational potential at second-order Ψ (2) Ψ (2) ′′ + 3HΨ (2) ′ = S m , (B.5) whose solution is with W (η) = W 0 /a 3 (a 0 = 1) the Wronskian obtained from the corresponding homogeneous equation. In Eq. (B.6) Ψ (2) m (0) represents the initial condition (taken conventionally at η → 0) deep in the matterdominated phase.
In order to give an analytical we will use some simplifications following Ref. [89,42]. First, for simplicity, we are going to neglect the ratio R wherever it appears, except in the arguments of the varying cosines and sines, where we will treat R = R * as a constant evaluated at the time of recombination. In this way we keep track of a damping of the photon velocity amplitude with respect to the case R = 0 which prevents the acoustic peaks in the power-spectrum to disappear. Treating R as a constant is justified by the fact that for modes within the horizon the time scale of the ocillations is much shorter than the time scale on which R varies. If R is a constant the sound speed is just a constant c s = 1/ 3(1 + R * ), and the sound horizon is simply r s (η) = c s η.
Second, we are going to solve for the evolutions of the perturbations in two well distinguished limiting regimes. One regime is for those perturbations which enter the Hubble radius when matter is the dominant component, that is at times much bigger than the equality epoch, with k ≪ k eq ∼ η −1 eq , where k eq is the wavenumber of the Hubble radius at the equality epoch. The other regime is for those perturbations with much smaller wavelenghts which enter the Hubble radius when the universe is still radiation dominated, that is perturbations with wavenumbers k ≫ k eq ∼ η −1 eq . In fact we are interested in perturbation modes which are within the horizon by the time of recombination η * . Therefore we will further suppose that η * ≫ η eq in order to study such modes in the first regime. Even though η * ≫ η eq is not the real case, it allows to obtain some analytical expressions.
Before solving for these two regimes let us fix our initial conditions, which are taken on large scales deep in the radiation dominated era (for η → 0). During this epoch, for adiabatic perturbations, the gravitational potentials remain constant on large scales (we are neglecting anisotropic stresses so that Φ (1) ≃ Ψ (1) ) and from the (0 − 0)-component of Einstein equations where ω 0 = kc s .

C.2 Perturbation modes with k ≪ k eq
This regime corresponds to perturbation modes which enter the Hubble radius when the universe is matter dominated at times η ≫ η eq . During matter domination the gravitational potential remains constant (both on super-horizon and sub-horizon scales), as one can see for example from Eq. (B.1), and its value is fixed to Φ (1) (k, η) = 9 10 Φ (1) (0), where Φ (1) (0) corresponds to the gravitational potential on large scales during the radiation dominated epoch. Since we are interested in the photon anisotropies around the time of recombination, when matter is dominating, we can perform the integral appearing in Eq. (C.10) by taking the gravitational potential equal to its value during matter domination so that it is easily computed This regime corresponds to perturbation modes which enter the Hubble radius when the universe is still radiation dominated at times η ≪ η eq . In this case an approximate analytical solution for the evolution of the perturbations can be obtained by considering the gravitational potential for a pure radiation Notice that the solutions (C.20)-(C.21) are actually correct only when radiation dominates. Indeed, between the epoch of equality and recombination, matter starts to dominate. Full account of such a period is given e.g. in Section 7.3 of Ref. [42], while its consequences for the CMB anisotropy evolution can be found e.g. in Ref. [90]. We refer to [39], where an alternative way to solve for the acoustic oscillations in this regime is displayed which turns out to be useful for the corresponding computation at second order in Sec. 7.5.

Symbol
Definition Equation