{'fraction_non_alphanumeric': 0.035076682422552456, 'fraction_numerical': 0.011969939130139082, 'mean_word_length': 4.323106163453706, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmic microwave background bispectrum on small angular scales

This article investigates the non-linear evolution of cosmological perturbations on sub-Hubble scales in order to evaluate the unavoidable deviations from Gaussianity that arise from the non-linear dynamics. It shows that the dominant contribution to modes coupling in the cosmic microwave background temperature anisotropies on small angular scales is driven by the sub-Hubble non-linear evolution of the dark matter component. The perturbation equations, involving in particular the first moments of the Boltzmann equation for the photons, are integrated up to second order in perturbations. An analytical analysis of the solutions gives a physical understanding of the result as well as an estimation of its order of magnitude. This allows to quantify the expected deviation from Gaussianity of the cosmic microwave background temperature anisotropy and, in particular, to compute its bispectrum on small angular scales. Restricting to equilateral configurations, we show that the non-linear evolution accounts for a contribution that would be equivalent to a constant primordial non-Gaussianity of order fNL~25 on scales ranging approximately from l~1000 to l~3000.


I. INTRODUCTION
The cosmic microwave background (CMB) offers a unique window on the physics of the early Universe, and in particular on inflationary models. The angular power spectrum of the CMB anisotropies has been extensively used to set constraints on the shape of the inflationary potentials; see e.g. Ref. [1]. The statistical properties of the temperature anisotropies and polarisation depend both on the inflationary period during which they were created and on the physics at play after Hubble-radius crossing and during the recombination. At linear order in metric perturbations, those latter physical processes amount to affect the metric perturbations by a multiplicative transfer function. The characteristic features observed in the temperature anisotropy spectrum originate from the development of acoustic oscillations that this transfer function encodes. The overall amplitude of the metric perturbation and its scale dependence are however determined by the inflationary phase.
At linear order, the calculation of the transfer function -and hence the detailed shape of the temperature power spectra -for generic inflationary models requires the identification of the relevant degrees of freedom during inflation (see e.g. Refs. [2,3,4]), as well as a full resolution of the dynamics up to recombination time. All these aspects are now fully understood (see e.g. Refs. [5,6] and references therein).
At this level of description, the metric perturbations are linearised so that the non-linear couplings that are inherently present in the Einstein equations are ignored. Therefore models that predict Gaussian initial metric fluctuations are expected to induce cosmic fields with Gaussian statistical properties. This is a priori the case for generic models of inflation. It has to be contrasted to models with active topological defects -such as cosmic strings -that have soon been recognised as a source of large non-Gaussianities [7,8,9,10]. The current data however clearly favor only mild non-Gaussianities although those might be larger than those induced by pure gravity couplings. This is not the case for single field slow-roll inflation for which it has been unambiguously shown in Ref. [11] that it can produce only very weak non-Gaussian signals, that are bound to be overridden by the gravity induced couplings. It has however been realised that some models of inflation might produce significant deviation from Gaussianity in the context of multiple-field inflation [12,13,14,15,16,17,18,19,20] or with non-standard kinetic terms [21]. The question of the observation of primordial non-Gaussianities is largely open.
In general however, primordial deviations from Gaussianity are in competition with the couplings induced during the non-linear evolution of the cosmic fields. It has triggered general studies aiming at characterising the bispectrum to be expected in the observation of the cosmic microwave background temperature anisotropies and polarisations whether it arises from inflation or from subsequent effects.
This task is multi-fold. It requires a proper identification of the mode couplings (at the quantum level) during the inflationary phase -so taking into account the usual gauge freedom -as well as a second order treatment of the post-inflationary evolution. While the former has been set on firm ground [11,22], the latter issue is still largely unexplored. This article proposes both numerical and analytical insights into it.
Hereafter, we assume that on super-Hubble scales, the only significant scalar perturbations are adiabatic and that they obey a nearly Gaussian statistics [47]. To be more precise, they are described in Fourier space by a single variable ζ 0 (k), k being a comoving wave-number, that satisfies and where "sym." stands for the two other terms obtained by permutation of the wave-numbers. This defines the primordial power spectrum P ζ (k) and the primordial mode coupling amplitude [48] f ζ NL . Considering an observable quantity θ related to the perturbation variables, the effect of evolution can generically be recapped [49] as where T (1) θ (k) is the linear transfer function and T (2) θ (k 1 , k 2 ) is the second order transfer function. θ can be thought as being e.g. the observed temperature anisotropies, but it could also stand for the CMB polarisation, or even cosmic shear surveys. When computing the bispectrum of θ there will be a contribution from the mode couplings induced by the second order transfer function T (2) and the possible initial non-Gaussianities, θ (k 1 )T (1) θ (k 2 )T (1) θ (k 3 ) P ζ (k 1 )P ζ (k 2 ) +sym. (4) where f θ NL is related to the second order transfer function by The full derivation of the details of T (2) θ is a fantastic task. It requires an understanding of the metric fluctuations behaviour at second order, from radiation dominated super-Hubble scales to matter dominated era at sub-Hubble scale, as well as a comprehension of the physics of recombination -through the Boltzmann equation -at a similar order. Such a task has been undertaken by several authors [50] and the multitude of effects at play needs to be sorted out. So, the goal of this article is not to provide an end to end calculation of T (2) θ , but to show that on small scales one can extract the dominant terms in order to get an insight into this physics at second order.
Modes coupling due to gravitational clustering is, by far, not a novel subject. It can be traced back to early works by Peebles [23] where the function T (2) for the non-linear sub-Hubble evolution of cold dark matter field (CDM) during a matter dominated era was derived. General modes coupling effects, within the same regime, has been extensively studied in the eighties and nineties where a whole corpus of results has been obtained (see e.g. Ref. [24] for an exhaustive review). On sub-Hubble scales, the second order mode coupling function for the gravitational potential reads in the particular case of an Einstein-de Sitter universe (here a is the scale factor and H the Hubble parameter). This well-established result proved to be useful for observational cosmology. The angular modulation it exhibits has indeed been observed in actual data sets; see e.g. Ref. [25]. The fact that on sub-Hubble scales, that is k 2 ≫ H 2 a 2 , the non-Gaussianity is driven by the nonlinearities of the CDM sector, and that they can start developing even before equality is one of the leading ideas of the present study. Indeed temperature anisotropies on small angular scales -i.e. beyond the first acoustic peak -mostly trace the gravitational potential [51] long after it has entered a sub-Hubble evolution and already during the matter dominated era. It is then natural to expect that the temperature anisotropies should be substantially determined by a form close to that of Eq. (6).
The goal of this paper is to evaluate how close we are from the behaviour (6) depending on scales, to which extent the temperature anisotropies trace this form and finally to estimate the amplitude of the temperature bispectrum on small angular scales. In this work two approaches will be compared: a full numerical integration of the second-order equations presented in § II, where the main approximation lies in the modelisation of the Compton scattering collision term at second order, see Eq. (30), and an approximate analytical resolution discussed in § III.
The bottom line of our analysis is that on small angular scales, the density perturbation of the cold dark matter starts to dominate the Poisson equation so at the time of decoupling we can assume that the system is split in (1) the evolution of CDM and (2) the evolution of the photons-baryons plasma which develops acoustic oscillations in the gravitational potential determined by the CDM component. As we shall demonstrate, at second order the dominant term of the temperature fluctuations is driven by the second order gravitational potential. Our approximation requires to consider a regime in which the Silk damping is efficient, that is wave-modes larger than the damping scales, hence corresponding to multipoles roughly larger that 2000. This picture will be shown to be in agreement with the numerical estimation (see § III D). We then proceed in § IV by a computation of the bispectrum in which we show that, for equilateral configurations, the non-linear dynamics has an amplitude equivalent to that of a primordial non-Gaussianity with constant f NL of order 25. A back-of-the-envelop argument allows us to understand the magnitude of this number.

II. PERTURBATION THEORY
This section is devoted to the presentation of the perturbation equations, up to second order, and of the initial conditions used in our study. We set the main notation and describe the background dynamics in § II A, and we define the perturbation variables in § II B. The perturbation equations and initial conditions are then presented in § II C and II D respectively.

A. The background dynamics
The background space-time is described by a Friedmann-Lemaître metric with scale factor a and cosmic time t. It is convenient to rescale the scale factor such that where ρ m and ρ r are the background matter and radiation energy densities respectively. The matter energy density can be decomposed as the sum of a cold dark matter component, that does not interact with normal matter, and a baryonic component, that can be coupled to radiation by Compton scattering prior to decoupling. We thus set ρ m = ρ c + ρ b where ρ c and ρ b refers to CDM and baryons respectively. It follows that when we neglect the contributions of the spatial curvature and of the cosmological constant, which are negligible for the whole history of the Universe but very recently. H ≡ a ′ /a is the conformal Hubble parameter and a prime refers to a derivative with respect to the conformal time η defined by dt = adη. H eq is the value of the H at equality, that is when y = 1.
The equation of state of the background fluid, composed of a mixture of non-relativistic matter and radiation, is w = 1/[3(1 + y)] and the density parameters of matter and radiation are Ω m = y 1 + y , Ω r = 1 1 + y (8) and indeed Ω c = Ω m (1 − f b ). Equality takes place at y = 1, from which we deduce that where Θ 2.7 ≡ T 0 /2.7 K is the temperature of the CMB today, and h the value of the Hubble constant in units of 100 km/s/Mpc. Equation (7) evaluated today implies that H eq ∼ H 0 √ 2y 0 so that The last scattering surface corresponds to a redshift [1] 1 + z LSS = 1090 ± 1 = y 0 /y LSS , and is mildly dependent of Ω c0 and Ω b0 . This implies that y LSS ∼ 3.3.
In Fourier space, a mode is super-Hubble when kη ≪ 1 and sub-Hubble otherwise. The mode becoming sub-Hubble at equality corresponds to a comoving wavelength of if we choose units such that a 0 = 1. We also introduce the parameter which will be useful to describe the physics of the baryons-photons plasma.

B. Perturbation variables
We focus on the dynamics of scalar perturbations (see e.g. Refs [26,27,28] for analysis of vector and tensor modes generation at second order). In Newtonian gauge, we can expand the metric as where Φ and Ψ are the two Bardeen potentials.
The various fluids contained in the universe will be described at the perturbation level by their density contrast δ and their velocity field. For the latter, we decompose the time-like tangent vector to the fluid worldlines according to where the first term accounts for the background Hubble flow. The perturbation v µ is further decomposed as v µ = (v 0 , v i ) with v i = ∂ i v and v 0 is constrained by the normalisation u µ u µ = −1.
When dealing with perturbations beyond first order, we assume that the perturbation variables are expanded according to where X (1) satisfies the first order field equations while the second order equations will involve purely second order terms, e.g. X (2) as well as terms quadratic in the first-order variables, e.g. [X (1) ] 2 . Thus, there shall never be any ambiguity about the order of perturbation variables involved as long as we know the order of the equation considered, and consequently we will usually omit the superscprit (1) or (2) which specifies the order of the perturbation. For a general discussions on second order perturbations and gauge issues, we forward to Refs. [29,30,31] C. Perturbation equations In this article we shall focus on the CDM-radiationbaryons system. Each component has a constant equation of state, that is P = wρ with w ′ = 0 so that c 2 s = w. Non-relativistic matter is described by a pressureless fluid with w m = w b = w c = 0 and radiation satisfies w r = 1 3 . In full generality, the evolution of each component can be obtained from the Boltzmann equation satisfied by the distribution function f a (x µ , p ν ) for this matter component. The stress-energy tensor can then be defined by integrating over momentum as where π + (p) is the volume element on the tangent space in x α such that p µ is non-spacelike and future directed (see e.g. Refs [29,32]) The first moment of the Boltzmann equation then gives a conservation equation of the form [32] where F ν a describes the force acting on the fluid labelled by a and satisfy F ν a u ν = 0 and a F ν a = 0, which is nothing but the action-reaction law (equivalently obtained from the Bianchi identity). Projecting along and perpendicular to u µ , we can extract respectively the continuity and Euler equations.

Linear order
Linear order calculations are used in particular to set the source terms of the second order equations. We closely follow the standard calculations and the main ingredients are recalled here. At linear order, the continuity equation for a fluid labelled by a takes the form [3,5,6] δ ′ a + 3H(c 2 s,a − w a )δ a + (1 + w a ) (∆v a − 3Ψ ′ ) = 0 , (17) while the Euler equation where π a is the contribution of the anisotropic pressure. In deriving Eq. (18), we have decomposed the force term vanishes at the background level, F a is gauge invariant. w a and c s,a are respectively the equation of state and sound speed of the component a.
In our analysis, we consider three components. Dark matter (label c) is described by a perfect fluid with w c = 0 and π c = 0 interacting only through gravity (F c = 0). Baryons and photons are coupled through Compton scattering so that F b and F r do not vanish. The actionreaction law (or equivalently the conservation of the total 1 3 Fr πr kinetic (8 moments) stress-energy tensor of matter) implies that F r = −F b , from which we deduce that At linear order, it is easily shown that where with n e being the free electrons number density and σ T the Thomson scattering cross-section. It follows that baryons will be described by a fluid (w b = 0 = π b = 0) interacting with radiation. In general radiation enjoys a non-vanishing anisotropic pressure (π r = 0) and should actually be describe by the full Boltzmann hierarchy. For the linear order calculations we choose to extend the fluid description by including the eight first moments of this hierarchy, including polarisation (see Appendix A for these equations). Our choices for the modelisation of the matter sector are summarised in table I.
The Einstein equations reduce to the set When the anisotropic pressure can be neglected, and in particular in the tight coupling regime discussed below, Eq. (24) implies that Φ = Ψ. Note that in this analysis we actually ignore the neutrinos the effect of which is thought to be marginal on the qualitative results we will obtain.

Second order
At second order, any first order equation, schematically written as D[X (1) ] = 0, of the first order perturbation variables X (1) will take the general form D[X (2) ] = S where S is a source term quadratic in the first order variables.
For the continuity and Euler equations respectively read As long as CDM is concerned, these source terms and Eqs. (17)(18) gives the full second order evolution of the fluid.
As already seen at first order, the fluid equations for the baryons and photons must include interaction terms, that is F (2) r and F (2) b , that derive from the Compton scattering collision term entering the Boltzmann equation for the radiation.
In the baryon rest-frame, this collision term includes only two types of contributions [33,34,35]. First, there is a term involving first order perturbation quantities and whose form is where C (1) is the first order collision term. This contribution involves the fluctuation of the visibility function, that is of the electron density n e and of the ionisation fraction x e . It accounts for the fact that a hotter or denser region decouples later. Its typical magnitude is of order 4C (1) δn e /n e . Second, there is a term involving second order perturbation variables and whose form is The forces derived from these two terms will satisfy by construction the action-reaction law (19), and this holds in any reference frame and at any order. This explains why the computation is easily carried out in the baryons rest-frame [34]. Then, when changing frame from the baryons rest-frame to the cosmological frame, where the computations are actually carried out, a second series of terms appears. They are of the form (1) ], where f (0) and f (1) are the background and first order distribution functions as well as similar terms for the polarisation (see Ref. [34] for the exact form of these terms.). Now, the contribution (28) is proportional to the collision term at first order. This implies that it will thus be negligible as long as tight coupling between baryons and photons is maintained at first order, i.e. as long as τ ′ /k ≫ 1. We are thus left only with the contribution (29). This term thus enforces the tight coupling regime at second order. Then, as long as tight coupling is effective, it is obvious that the second series of terms arising from the change of frames should compensate each other to give a vanishing contribution. We shall thus model the interaction term entering the Euler equations by that is by assuming that it keeps the same functional form as at first order. In conclusion, the continuity and Euler equations for baryons and photons with the interaction term (30) are the exact fluid limit of the full Boltzmann equation at second order as long as tight coupling is effective. This implies that at second order, and similarly as at first order, the two tightly coupled fluids are equivalent to a single perfect fluid; see section III A. Again, this stems from the action-reaction law which implies that Eq. (19) has to hold at any order and in any reference frame and that there cannot appear any external force acting on the resulting effective fluid since it is only coupled to other matter components through gravitation. When τ ′ /k becomes of order unity, tight coupling stops being effective. But thanks to Silk damping, the terms of (1) ], will still be negligible compared to the one we kept to obtain Eq. (30).
Let us now turn to the Einstein equations (22)(23)(24)(25). Their source terms read respectively This provides all the source terms appearing at second order.

Note on our conventions
Since the first and second order equations are conveniently solved in Fourier space, the quadratic terms in the source terms can be written as a convolution on the wave-numbers k 1 and k 2 such that k 1 + k 2 = k. For simplicity of notations we will write the integral factor of the convolution as We also define µ = k 1 · k 2 /(k 1 k 2 ). Unless explicitly specified, we choose the convention of the Fourier transform in which the factors of (2π) n/2 are symmetric for the Fourier transform and its inverse, n being the dimension of space.

D. Initial conditions
To integrate this system of equations, we need to set the initial conditions both for the first and second order variables deep in the radiation era for super-Hubble modes at the initial time η init , that is modes such that kη init ≪ 1.

First order
At first order, we rely on the comoving curvature perturbation which is constant on super-Hubble scales. It is well-known [5,6] that for a perfect fluid with a timedependent equation of state (c 2 s = w), the comoving curvature perturbation, defined by is conserved on super-Hubble scales for adiabatic perturbations. Inflationary models predict the initial power spectrum of R (1) on those super-Hubble scales from which one can deduce the power spectrum of the gravitational potential. If we choose η init such that the decaying mode is negligible and Ψ (1) is constant on super-Hubble scales then, still neglecting the anisotropic pressure, Deep in the radiation era, this implies that for modes such that kη init ≪ 1. Since the density contrast of the total fluid δ ≃ δ r deep in the radiation era. From Eq. (22) we deduce that Now, assuming adiabatic initial perturbations, we must have . Since baryons and photons are tightly coupled deep in the radiation era, we deduce that This completely fixes the initial conditions for the set of first order perturbation equations.

Second order
At second order the previous procedure can be generalised [36,37]. It was shown that on super-Hubble scales and for adiabatic perturbations, the variable is a conserved quantity on super-Hubble scales for adiabatic perturbations. Once the decaying modes are negligible so that Φ (2) and Ψ (2) are constant, we can express them in terms of R (2) as In the case where the initial perturbations have been generated during a phase of one field inflation in slow-roll, it has been shown [11] that, for the variable defined in Eq. (38), R (2) ≃ −2R (1)2 + O(slow-roll parameters). Under this hypothesis, we can use the constancy of R (2) on super-Hubble scales and Eqs. (39)(40) to derive the initial conditions satisfied by the two Bardeen potentials at second order at an initial time η init deep in the radiation era. Up to small corrections of the order of the slow-roll parameters, they are given by We also need to determine the initial conditions for the energy density contrasts and the velocities of the different matter components. Again, we assume adiabatic initial conditions, which means that the total fluid behaves like a single fluid. While at linear order the pressure and density perturbations are simply related by the sound speed at second order we have This implies that the adiabaticity conditions at second order reads where use of the first order adiabaticity condition was made to get the last equality. It can also be shown from the perturbation equations that the condition (45) remains valid on super-Hubble scales, hence giving a conservation of the baryon-photon entropy on large scales, exactly as at linear order. Following the same procedure as at first order, we first use the Poisson equation (22) and Eq. (25) at second order to get that This completely fixes the initial conditions for the set of second order perturbation equations.

E. Integrating the evolution equations
Working in Fourier space, we first integrate the first order equations, that is • Eqs. (17)(18) for CDM and baryons assuming a source term of the form (20) for baryons; • the first eight moments of the Boltzmann equation for the radiation including the contribution of the polarisation. These equations are detailed in Appendix A; • the Einstein equations (22)(23)(24)(25).
The initial conditions are detailed in § II D 1. Technically, we have recast all these equations in order to use y as time variable and we remind that the time of decoupling is of the order of y = 3. We also remind that the modes of interest, that is k > k eq. , becomes sub-Hubble at y < 1. This first integration thus allows us determine all the first order perturbations as a function of time and wavenumber. An example of the results of the first order integration is presented on Fig.7 (Appendix C). Now, at second order, we integrate the same system of equations but supplemented by the source terms which are determined by the solutions of the previous integration. We thus solve • Eqs. (17)(18) for CDM and baryon with the source terms (26) and (27) respectively. We recall our main hypothesis which states that the coupling of baryons to radiation can be described by the interaction term (30); • radiation is described by the first four moments of the Boltzmann equation with the hypothesis (30) for the collision term; • the Einstein equations (22)(23)(24)(25) with the source terms (31)(32)(33)(34).
The initial conditions are detailed in § II D 2 so that we are finally able to compute the evolution of the perturbation variables with y for any wave-number. Fig. 1 shows the evolution of the two second order gravitational potentials. In particular, it shows that the solution is driven, as expected, toward Φ (2) = Ψ (2) . It is to be noted that the convergence takes place before equivalence, at y ≡ 1, as stressed in the following. On the other hand Fig. 2 depicts the evolution of the velocities and density contrasts and shows that the photonsbaryons plasma can safely be described as a single fluid, almost until the decoupling (shaded area on the figure).

III. ANALYTICAL INSIGHT
Before we proceed to describing the outcome of our numerical integrations, and in order to gain some insight into the physics of this intricate system, we present some analytic descriptions of its solutions.

A. Heuristic argument and hypothesis
Let us first assume that f b ≪ 1 so that the universe is mainly dominated by non-interacting cold dark matter and radiation components. When, in the radiation era, the CDM component is completely negligible the gravitational potential is determined by the density contrast of radiation. The latter however develops oscillations after Hubble-radius crossing while those in the CDM fluid increases. It follows that, while still formally in the radiation era (ρ r > ρ c ), the cold dark matter component is actually driving the gravitational potential. It then acts as an external driving term in the evolution equation of radiation. In such a scenario we then expect non-linearities that develop in the CDM sector to be transferred first to the gravitational potential and then to the radiation density fluctations.
To make this heuristic argument more quantitative we will thus assume that • we can first study the CDM-plasma system to determine the gravitational potential at second order where for simplicity the plasma is assumed to be radiation dominated; • then study the acoustic oscillations of the baryonphoton plasma driven by the gravitational potential derived this way, both at first and second order in the perturbations.
For the sake of simplicity we will work in the tight coupling approximation. We recall that the time of decoupling is both the time at which this tight coupling regime ceases to be valid and the time at which the radiation temperature is observed. The tight coupling approximation amounts to saying that the coupling terms F r and F b are so large that radiation and baryons behave as single fluid. It ensures that the two fluids have the same peculiar velocity (v r = v b ) and implies that the anisotropic pressure of radiation vanishes (π r = 0). From Eq. (17), it implies that At linear order, eliminating F a in Eq. (18) for radiation and baryons leads for the photons-baryons plasma to the continuity equation and the Euler equation where we have introduced the density contrast of the plasma δ pl = δρ pl /ρ pl with In the particular case at hand, it reduces to and the velocity perturbation are given by The equation of state and sound speed of the plasma are easily obtained from the fact that P pl = P r . They are explicitly given by and are time-dependent quantities (simply because the relative contribution of the two components changes with time). At second order, the density contrast and velocity perturbation of the plasma are given by b .
It can be checked that the plasma follows the equation (48) and (49) but supplemented with the source terms where S c,a=pl and S e,a=pl stand for S c,a and S e,a in which we take the fluid to be the plasma, that is δ a = δ pl , w a = v pl , etc.
The validity of the tight coupling approximation at first and second order can be checked from our numerical integration (which indeed does not make this assumption) respectively on Fig. 2 and Fig. 8 for the second and first orders.
B. CDM-radiation system

First order
Deep in the radiation era, the gravitational potential is mainly determined by the radiation density contrast and decays on sub-Hubble scales. The contribution of matter is negligible in the Poisson equation and it follows that the potential is given by where j 1 is a spherical Bessel function of order 1 and x = kη. The density contrast of radiation is given by δ r = −2Φ on super-Hubble scales (see § II D 1).
Let us now turn to the evolution of the CDM fluid during the radiation era. In terms of the variable y the continuity and Euler equations lead tö with a driving force determined by the gravitational potential where a dot stands for a derivative with respect to y. The general solution of Eq. (54) is of the form δ c (k, η) = A+B ln x+δ part where δ part is a particular solution given by where S Φ is given by Eq. (53) as long as δρ c ≪ δρ r . For y ≪ 1, the contribution of the particular solution is negligible so that A ≃ − 3 2 Φ(η init ) and B ≃ 0. The solution (53) shows that Φ vary mainly when η ∼ k −1 so that δ part ∼ A + B ln x. A and B can be obtained semi-analytically and are well approximated by A ≃ 6 and B ≃ −9 so that δ c ≃ Φ(k, η init )(−4.5 + 9 ln x).
This solution is valid as long as δρ c ≪ δρ r in the Poisson equation. However, on sub-Hubble scales δ r remains constant while, as we just saw, δ c grows logarithmically. Their contribution in the Poisson equation then become to be of the same order when y ∼ y ⋆ (k), where y ⋆ is solution of y ⋆ [−4.5 + 9 ln( √ 2y ⋆ k/k eq )] ∼ 6/(1 − f b ), where use has been made of x = √ 2ky/k eq as long as y ≪ 1,. The solution of this equation is depicted on Fig. 3. For most of the scales of interest, i.e. for k ≫ k eq. , the contribution of the CDM in the Poisson equation is dominant before equality, i.e. y ⋆ < 1).
For these modes, which became sub-Hubble during the radiation era, we shall consider that CDM dominates in the Poisson equation and neglect the contribution of the radiation density perturbation, so that the Poisson equation takes the form Neglecting the contribution of baryons, since their density contrast cannot grow because they are tightly coupled to the radiation, Eq. (54) then takes the form of the Mészáros equation [38] δ c + 2 + 3y 2y(1 + y)δ c − 3 2y(y + 1) δ c = 0 .
Its two solutions are a growing mode D + (y) = y + 2/3 (58) and the decaying mode

Second order
At second order, as at first order, while deep in the radiation era the gravitational potential generated by the density contrast of radiation decays when a mode becomes sub-Hubble. At this order the gravitational potential satisfies with where the source terms are given by Eqs. (31)(32)(33)(34). Up to a fast decaying solution, the general solution is with the Green function On sub-Hubble scales, the leading terms in S r are those quadratic in the first order velocity, of the form ∝ H 2 ∂ i v∂ i v ∼ η −2 , which behave as η −2 . All other terms in S r behave, at best, as k −2 η −4 . Using the first order solution, the second order gravitational potential asymptotically behaves as and it can be checked that this term is indeed regular in µ 2 = 1. Taking the homogeneous solution into account, Ψ (2) decays as (kη) −2 on sub-Hubble scales. Now, the evolution of the density contrast of CDM follows, using y as the time variable, the evolution equation where the source term is given by As in the previous section for the first order perturbations, CDM density perturbations grow faster than those of radiation so that, for any mode k that became sub-Hubble before η eq , there exists a time of order y ⋆ [k] such that for y > y ⋆ [k] the gravitational potential at second order is determined by the cold dark matter. Whereafter, even though we are still in the radiation era, we can neglect the contribution of the density perturbation of radiation so that the Poisson equation becomes In this regime, the evolution of the density contrast of CDM at second order can be derived from a second order Mészáros-like equation, in a similar way as at first order. Using Eq. (65) and the fact that the main contributions in S δc in this regime come from Eq. (64) takes the form with Now we shall neglect the effect of baryons, that is f b . This equation could then be called the second order Mészáros equation and describes a growth of CDM density perturbation in a regime where radiation dominates the dynamics of the background while its density perturbations are negligible in the Poisson equation. The Green function associated to this equation is obtained to be so that the general solution of Eq. (67) is In the limit where y ≫ 1, y ′ ≫ 1, the Green function behaves as G(y, y ′ ) ≃ 2 5 y 1 − y ′ y

5/2
, and the source term as with the Kernel (70) In the limit y ≫ 1, the particular solution dominates and our solution converges toward that is toward the standard result (6) describing the collapse of cold dark matter in a matter dominated era. It is compared to its expected late time behaviour (72). It is to be noted that the convergence toward this solution is extremely rapid and takes place as soon as equality is reached, e.g. y = 1. The results correspond to k1 = 6keq and k2 = 12keq. The difference in the amplitude of the function is due to the fact that the baryons component has been neglected in the derivation of Eq. (72) .
The second order gravitational potential is then obtained from the Poisson equation up to terms of order O(f b ). The convergence towards the solution (72) is explicitly depicted on Fig. 4 where the behaviour of the exact (numerically integrated) second order potential as a function of time (and angle) is compared to its expected late time behaviour (72). As detailed above, this solution is a better approximation for larger wave-numbers and at large y since we converge to this solution for y > y ⋆ (k). On this figure it can be observed that the convergence is extremely rapid and that the full kernel structure, including its angular dependence in (72), is indeed to be observed in Φ (2) .

C. Baryons-radiation system
We now want to understand the behaviour of the baryons-photons plasma, and in particular of its acoustic oscillation, in the regime in which the gravitational potential is determined by the solutions of the previous section.
We restrict our analysis to the tight coupling regime. And since it occurs for y < y LSS , our solution will gain in accuracy when the period between CDM domination, y = y ⋆ (k), and the last scattering surface, y = y LSS , is large, that is on the smallest scales.

First order
The computation at first order is well known [39] and we review its main steps to compare to the more unexplored second order case.
In the fluid limit, the baryons and photons both obey a continuity and conservation equations (17)(18) with a source term and in which the anisotropic stress of radiation can be neglected because of the tight-coupling approximation. As discussed in § III A, in the tight coupled regime, that is for k/τ ′ ≪ 1, it leads to the wave equation where R is defined in Eq. (13), the sound speed is given by Eq. (51). This is a wave equation with a forcing term on the r.h.s. which describes the oscillations of the plasma. For small wavelength modes, the variation of R and Φ is small compared to the period of the wave so that we can construct an adiabatic solution by resorting on a WKB approximation; see e.g. Ref. [39] for details. Defining The velocity field, v r = v b can then be determined from the Euler equation (49). This solution neglects the Silk damping effect that can be described by adding terms in k 2 /τ ′ in Eq. (73) so that the solution is exponentially suppressed by a factor where the damping scale is of the order of k D ∼ 15k eq . Since Φ decreases as (k/k eq ) −2 , we conclude that for large wavenumber, Θ SW (k, η) goes rapidly to zero.

Second order
The former approach can be generalised at second order, but the behaviour of Θ SW will change mainly because the second order version of Eq (73) has a r.h.s. which is steadily growing in the range of interest, i.e. after η eq and on large k.
At second order, Eq. (73) will also contain terms coming from the second order Liouville equation [29] of the form This source term involves terms which are quadratic in the fluid perturbation variables (δ, v) and the potentials (Φ, Ψ). The former are exponentially suppressed due to Silk damping and the latter decrease as (k/k eq ) −2 . We can thus neglect this source term as long as we focus on small scales. Defining, the solution for Θ (2) SW will be similar to Eq. (75). When Silk damping is taken into account, Θ (2) SW +RΦ (2) is exponentially suppressed, exactly as at first order. The main difference with first order arises from the fact that the second order gravitational potentials are driven toward Ψ (2) ≃ Φ (2) ∼ (kη) 2 after equivalence (y > y eq ) so that we expect that Now, the velocity of radiation is given by Since Θ (2) SW +RΦ (2) is exponentially suppressed due to the Silk damping, we expect that the Doppler contribution is also negligible.
The equations (72) and (77) are the central results of the analytic insight of the non-linear regimes of the photons-baryons-CDM system at second order. They give the behaviour of the second order gravitational potentials at the time of decoupling together with the response of the photon-electron plasma. We also conclude that we expect Θ (2) SW = −RΦ (2) to dominate the CMB temperature anisotropies on small angular scales (see Appendix B for a discussion of the integrated Sachs-Wolfe contribution).
The bottom line of our analytic estimates is that on small angular scales, i.e. k/k eq ≫ 1, the density perturbation of CDM starts to dominate the Poisson equation from y ⋆ (k) so that from this time to decoupling we can assume that the system is split in (1) the evolution of CDM and (2) the evolution of the photons-baryons plasma which develops acoustic oscillations in the gravitational potential determined by the CDM component. Because of Silk damping, Θ SW + RΦ dies out on small scales. At first order this implies that Θ (1) SW ∼ −RΦ (1) which is suppressed by a factor k 2 eq /k 2 due to its evolution in the radiation era prior to y ⋆ (k). At second order however, we still have for the same reason that Θ (2) SW ∼ −RΦ (2) but now this term roughly grows as (kη) 2 . Note that since k D is of order 15k eq and that k eq roughly corresponds to a multipole ℓ ∼ 160, we expect our analysis to give a good description of the system for ℓ 2400.

D. Comparison to numerics
We now turn to the description of the numerical solutions of the system described in § II E. Its solutions will be described in the light of the analytic description we just developed. Fig. 5 shows the result of numerical integrations for second order quantities of interest. They are compared to our approximate formula that, we recall, is expected to be valid in the tight coupling regime. It shows indeed that modes with k > k eq relax temporarily toward the solution (77). The exact solution exhibits though large oscillations that are thought to be due to the acoustic oscillations that are present in the plasma at first order, but their average turns out to coincide with the proposed analytic formula (long dashed lines) as long as a strong coupling is ensured. The Silk damping effect is observed to play a key role to actually damping the oscillations. This effect is all the more important that k is large. Note that the impact of the oscillations on the observational quantities will also be damped by the finite width of the last scattering surface. The wave-number corresponding to this width is of order of 10k eq which is smaller than the damping scale.
When the coupling becomes loose, log 10 (y) approaching 0.6, the numerical solution departs from the expected solution (and it converges toward 0) as the full Boltzmann hierarchy is now at play. We also depict the velocity term which can be checked to be negligible, as expected.
From these set of results we can then argue that the approximate analytic solution described in (77) captures the physics of the dominant terms of the CMB anisotropies on small angular scales. Here we have explicitly checked that this form is consistent with the physics of recombination when the collision effects are taken into account. We limit though the collision effects to their first order expression. We expect nonetheless that an exact calculation, up to second order, would not significantly alter our conclusions, the collision physics playing a role only during a limited period of time. Although this is certainly desirable to do such a calculation, this is beyond the scope of this paper whose goal is to estimate the order of magnitude of the non-Gaussianity on these

scales.
In the following we explore the observational consequence of such a finding on the temperature bispectrum at small scale.

A. Flat sky approximation
Since our approximations hold on small angular scales, it is amply sufficient to treat the sky as flat to compute the properties of the CMB anisotropies. We thus decompose the CMB temperature anisotropies in 2D-Fourier space as so that On the other hand Θ(n) can be expanded in Fourier modes as where D LSS is the angular distance of the last scattering surface given by D LSS = y0 yLSS 1/(H[u]u)du. k ⊥ is the projection of k on the sky, ı.e. k = k r e + k ⊥ where e is the direction of the (flat) sky.
In Eq. (81) we refer toΘ(k, η LSS ) as considering the observed CMB anisotropies as a superposition of spheres of temperature anisotropy weighted by the visibility function v(η, η LSS ) which peaks at η LSS . The angular power spectrum is nothing but the two dimensional power spectrum of Θ(l), We now need to determine Θ(k, η) in terms of the perturbation variables. The CMB temperature anisotropies are usually split as an intrinsic Sachs-Wolfe effect, a Doppler effect and an integrated Sachs-Wolfe contribution. As discussed in Appendix B the integrated Sachs-Wolfe contribution is expected to be negligible on the scales of interest in our study.
At first order, the Fourier component of temperature anisotropy for a mode k, in a direction e emitted at a comoving distance η 0 −η is dominated by the Sachs-Wolfe and Doppler terms, The Sachs-Wolfe term is related to the perturbation variables by while the Doppler term is given Θ Dop ≡ −vk.e so that Θ (1) with k r = k · e. Using Eq.(81), we deduce thus that withĝ (1) and where k ⊥ = l/D LSS . Using the definition of the initial power spectrum Φ(k, 0)Φ(k ′ , 0) = δ (3) (k + k ′ )P (k), we finally get that C l is given by This reproduces the main features of the CMB angular power spectrum, as checked on Fig. 9. At second order, and for the scales of interest (k ≫ k eq and k > k D ), we stress again that the knowledge of the exact expression of the terms quadratic in the first order variables in the second order Sachs-Wolfe effect are not needed since they are suppressed because of the Silk damping or because of the decaying of the potential during the radiation era. The analysis of § III C 2 shows that the Doppler term is much smaller than the intrinsic Sachs Wolfe term of Eq. (77). Thus, the main contribution to the second order temperature anisotropy is well approximated by (2) . (90) We deduce that withĝ (2) It can be rewritten in terms of the initial first order gravitational potential aŝ hence defining f (Θ) NL .

C. Numerical computation
In order to perform the previous integrals, we need to specify the initial power spectrum. We assume that the power spectrum is scale invariant and we normalize it using the results of WMAP, that is We then compute the bispectrum of an equilateral configuration for which all momentums are equals; l 1 = l 2 = l 3 . The only free parameter for such configuration is the norm l of the three vectors. The result is depicted on Fig. 6 and is compared to the bispectrum one would obtain from a initial constant f Φ NL assuming a linear transfer function. It appears that on scales that range from l = 1000 to l = 3000 the bispectrum resembles that of an effective constant primordial f Φ NL of order 25.
The order of magnitude of the amplitude of the bispectrum can be understood from the following rule of thumb for modes larger than k D . According to our analysis, considering the equilateral configuration where k 1 = k 2 = k, the second order temperature anisotropy on the last scattering surface is of order where we have assumed that, on average, the Kernel K(k 1 , k 2 ) is of order unity. Now, assuming a constant primordial f NL evolved with the linear transfer function, the second order temperature anisotropy would roughly be of order since, for these modes, the integral on the visibility function keeps only the average of the Sachs-Wolfe contribution. The ratio of the contribution of the non-linear dynamics compared to a primordial non-Gaussianity is 2 3 Since for large modes the gravitational potential has been decaying as (kη) −2 in the radiation dominated era, and growing logarithmically when the potential started to be determined by the cold dark matter component(see the analysis of section III B 1), we deduce that the first order transfer function is typically given by where A(k) is a steadily growing function. At the Silk damping scale we find numerically A(k D ) ≃ 10. We thus conclude that in the bispectrum, the evolution for l ≃ k D D LSS ≃ 2400 is equivalent to a primordial evolved linearly. This estimates of the order of magnitude is in complete agreement with Fig. 6 where it can be read that for multipole ranging from 2000 to 3000 the amplitude of the bispectrum is comparable with the one that would be obtained from a constant f NL ranging between 10 and 50.
There is no guarantee however that for arbitrary geometries the shape dependence of the temperature bispectrum would be that of a constant f NL . It is rather determined by the kernel shape of the form (77).

V. CONCLUSION
This article investigates the non-Gaussianity that arises in the CMB temperature anisotropies due to the post-inflationary non-linear dynamics during the radiation and matter dominated era. More specifically it aims at identifying the leading mechanisms and leading terms that determine the shape of the CMB bispectrum on small angular scales.
The driving idea that we have pursued throughout the paper is that at small angular scales the second order CMB anisotropies trace the second order gravitational potential as it is shaped by the CDM component during its sub-Hubble evolution. To give support to this picture, we have developed both analytical insights into the joint evolution of the density potentials and the temperature fluctuations and numerical tools.
We have thus solved numerically the joint evolution equations of the cosmic fluids up to second order. We have been able to check that at the time of decoupling, the second order potential indeed traces its expected shape. This conclusion is summarised and illustrated on Fig. 4. At this stage, and as long as one restricts these results to the tight coupling regime, no approximations have been made. The accuracy with which the k-dependence of the matter dominated mode coupling kernel, i.e. Eq. (6), is recovered is truly remarkable. We stress that this is due to the fact that for the physics at work at small scales, i.e. k/k eq ≫ 1, the density perturbations of the CDM component start to dominate the Poisson equation much before equality. This implies that the non-linearities developed by the CDM can be transferred very efficiently to the gravitational potential even before the beginning of the matter era.
Determining exactly how this mode coupling kernel is actually transferred to the source term of the CMB anisotropies relies on further numerical integrations through the recombination era. At this stage, the only approximation we make concerns the Compton scattering collision term entering the Boltzmann equation for radiation at second order. We did not use its full second order expression but we argue that it can be reduced to its formal first order form, namely to Eq. (30). Such an assumption is clearly valid in the tight coupling regime. Actually this is the only term appearing in the collision term in the baryons rest-frame as long as tight coupling at first order is efficient. We then argue that when the coupling drops, Silk damping effects effectively suppress all other contributions. This leads to the behaviour depicted on Fig. 5 for the main source term of the temperature anisotropies. In particular Θ (2) SW , the monopole of the second order source term, is found to be attracted toward a non-vanishing and non-oscillatory term, −RΦ (2) , where we recall that R is the baryon to radiation ratio, while the dipole contribution, and thus the Doppler effect, vanishes. As a result the main contribution to the CMB temperature anisotropies at second order is found to be directly proportional to the second order gravitational potential. It has to be noted that the efficiency with which the second order term converges to this form is considerably accelerated by the Silk damping effects which efficiently suppress the oscillatory parts of the solution. This results are clearly illustrated on Fig. 5. We argued from order of magnitude arguments that, because the damping scale is typically of order 15k eq , our description shall be valid for l 2400. This observation is the basis of the main result of this paper. Actually,as Fig. 6 tends to show, it seems that this description may be valid at lower multipoles.
Finally we explore the consequence on the CMB bispectrum. For obvious reasons we use the small angle approximation to perform the numerical integrations. The bispectrum for equilateral configurations is illustrated on Fig. 6. We show that for these configurations, its amplitude corresponds to what a primordial non-Gaussian potential of f Φ NL of order 25 would have given (also for an equilateral configuration). As shown in the text, this number can easily be recovered from back-ofthe-envelop calculations. The first lesson that can be drawn from this result is that it gives a signal larger than what a model with a primordial f NL of order unity would give! The second lesson is that the l-dependence of the bispectrum is expected to be different from the one induced by primordial mode couplings. It is expected to have a specific shape as encoded in the CDM kernel expression.
In conclusion, this work offers a breakthrough insight into the physics of CMB in the non-linear regime and on small angular scales. It identifies what is, as we argued, the main small scale contribution of the bispectrum, hence filling the gap with the standard results that have been obtained in the weakly non-linear regime of gravitational clustering of dark matter. We did not check this result against a (yet non-existing) full second order Boltzmann code, and this is probably desirable, but we argue that, given the amplitude of the effects, all other contributions will be subdominant. With such a large signal, detection of this bispectrum should be easily within reach of future CMB experiments! unique and one could have chosen the Bardeen potential.
With such a choice however the f ζ NL incorporates only the inflation dependent couplings -f ζ NL is proportional to the slow roll parameter in single field inflation for instance. The other coupling terms induce by the change of variable can be incorporated into T (2) ; see Eq. (5) below.
[48] This is the expression for the bispectrum obtained assuming ζ could be expanded as ζ = ζG + f ζ NL ζGζG where ζG is assumed to obey Gaussian statistics. This is not however a valid description when the bispectrum originates from multiple-field couplings or from quantum calculation. The formal expression (2) is always valid though; see Refs. [22,41].
[49] Things are actually slightly more complicated since usually observables cannot be decomposed into 3D Fourier modes. The functions T θ (k) and T (2) θ (k1, k2) should then be thought as projection operators. This is in particular the case for temperature anisotropies and polarisations. This does not affect however the general point we want to make in this introduction.
[50] Early derivations are to be found in Ref. [33,42,43]. A more rigorous and comprehensive calculation -including a proper derivation of the Boltzmann coupling terms and taking into account the polarisation effects -is to be found in Refs. [29,34].
[51] At least to some extent as we shall see in the course of this paper.

APPENDIX A: DESCRIPTION OF RADIATION
To describe the evolution of radiation, we use the first moments of the Boltzmann hierarchy (see e.g. Refs. [44]) including polarisation. The hierarchy reads where the first moments are related to the fluid variables by The Boltzmann hierarchy is infinite and we truncated it after the multipole ℓ = 8, when computing the first order and after ℓ = 3 when computing the second order. In order to cut the hierarchy without numeric reflection [45], we use the free-streaming solution of this hierarchy, and use it to express in the last equation the multipole ℓ + 1 in function of the multipoles for ℓ and ℓ − 1. Explicitly, the closure relation reads

APPENDIX B: THE INTEGRATED SACHS-WOLFE EFFECT CONTRIBUTION TO THE SMALL SCALE BISPECTRUM
At linear order, the contribution of the integrated Sachs-Wolfe effect on small scales is usually small because the time dependence of the potential vanishes in the matter dominated era. This is no more the case at second order. It is thus legitimate to investigate the impact of the time dependence of the second order gravitational potential on the amplitude of the bispectrum.
At linear order the expression of the temperature anisotropies is while at second order an extra source term should be included. It is formally given by Θ ISW (l) = ηLSS 0 dη d dη Ψ (2) [x(η), η] + Φ (2) [x(η), η] (B2) assuming instantaneous recombination at η = η LSS . To estimate the magnitude of this effect, we assume that Ψ (2) = Φ (2) and that the time dependence is the one otained in Eq. (72), where the time dependence has been explicited. Obviously, expression (B2) gives an extra term contributing to the bispectrum. Consistently with our former analysis, let us evaluate this contribution in the small angle approximation. This leads to that has to be compared with Eq.(95). The integral over η can then be performed to give If one examines the UV convergence properties of this expression (for the integrals over k r ), it appears that the integral over k r1 + k r2 converges at a scale given by the inverse of η LSS , e.g.
due to the oscillatory behaviour of w, whereas the integral over k r1 − k r2 converges because of the power spectrum shape and therefore at a scale which is of the order of l 1 /η LSS or l 2 /η LSS (whichever is smaller). Thus if the power spectrum is approximated by a power law, where the spectral index n s varies a priori from 1 (at very large scale) to −3 at very small ones -it is a priori of the order of say −2 at the scales of interest -then the integral over k r1 − k r2 leads to the factor, ∞ −∞ dk P 2 Φ ( k 2 r + l 2 /η 2 LSS ) = l 2(ns−3)−1 √ πη 7−2ns LSS Γ 7 2 − n s Γ(4 − n s ) (B9) that is 63π/256 (η LSS /l) 11 for n s = −2. This is to be compared with the amplitude of the intrinsic effects we have computed. The latter differs in Eq. (95) because of the absence of filtering function w. The amplitude of the bispectrum is then roughly given by, so that its amplitude is dominated by the square of which is equal to (3π/8) 2 (η LSS /l) 10 for n s = −2. The ratio of the two contributions scales then as 1/(Rl) in favor of the intrinsic effect.

APPENDIX C: CHECK OF THE NUMERICAL INTEGRATION
We report in this appendix the results of the first order numerical integration. We first report in Figs. 7 and 8 the evolution of the perturbed quantities where it can be seen that for y > y ⋆ (k) the gravitational potential tends to be determined by the cold dark matter density perturbation. We also report the angular power spectrum obtained from the flat sky approximation using the expression (89). The linear dynamics is then used to calculate the bispectrum arising from a constant primordial f NL evolved linearly. It can be checked that the form obtained on Fig. 9 is completely consistent with the literature (see e.g. Ref. [46]).
r . It can be seen that the adiabaticity condition holds until recombination, hence justifying the approximation of § III. FIG. 9: Left: in solid red line, the angular power spectrum using our code and the flat sky approximation (which does not includes the late ISW). The black dashed line represents the spectrum obtained using CAMB, which also takes into account the late ISW, but no reionisation. Middle: The precision on small scale depends highly on the computation of the visibility function. On small scales, our transfer function is approximately 30% smaller than the one predicted by CAMB and is thus a good approximation. Right: the bispectrum obtained from a primordial constant f NL = 100 evolved linearly.