text:
\section{Introduction}
The three year data set of the Wilkinson Anisotropy Probe (WMAP) on Cosmic Microwave Background (CMB) temperature anisotropies
has offered a wealth of information about the evolution of the universe of unprecedented accuracy~\cite{kwmap3} and prepared the way for even more
ambitious future missions such as those of the {\it Planck} satellite~\cite{Planck} and CMB polarization observations. Such precise measurements
have been accompanied by increasing  theoretical efforts in predicting the signatures on the CMB from various cosmological scenarios. In particular
a lot of attention in the last years has been dedicated to the statistical properties of the CMB beyond the power spectrum,
in search for possible non-Gaussian signatures~\cite{review}.
Non-Gaussianity of the CMB owes its importance to the possibility of unveiling crucial aspects of the
physics of both the early and the late universe, which would be unreachable using only the CMB power spectra information. Different mechanisms for the
generation of the cosmological perturbations predict different amplitudes and shapes of primordial non-Gaussianity, thus a positive detection of NG or
an upper limit on its amplitude is powerful in discriminating among the various competing scenarios which would be indistinguishable otherwise~\cite{review}.
An illuminating example is given by the standard  models of inflation which typically predict a very low content of non-Gaussianity~\cite{ABMR,maldacena}
and as such they could be completely ruled out if a positive detection will be achieved.
CMB non-Gaussianity can also have a non-primordial origin due to secondary anisotropies
which arise when the CMB photons leave the last scattering surface and cross the large-scale structure of the Universe.
These include both secondary scatterings
(such as the thermal and kinetic Sunyaev-Zel' dovich effects, produced by the thermal and bulk motions of electrons in clusters,
and the Ostriker-Vishniac effect, due to bulk motions modulated by linear density perturbations) and gravitational secondaries such as the gravitational
lensing and the non-linear Integrated Sachs-Wolfe effect or, on smaller scales, the Rees-Sciama effect due to the non-linear evolution of the potentials.
Such non-linearities can give important and new information about the dark matter and the dark energy content of the universe~\cite{MGLM,Bern,
GS1,GS2,VerdeSDE,GBP,ks,fullT}. Statistics like the bispectrum (the Fourier transform of the three-point correlation function) and the trispectrum
(four-point correlation) can then be used to assess the level of primordial non-Gaussianity on various cosmological scales
and to discriminate it from the ones induced by the secondary anisotropies and systematic effects.

Up to now most of the attention has focused on the three-point statistics both for the primordial non-Gaussianity and for the analysis of the angular
bispectrum of the CMB temperature and polarization anisotropies~\cite{review,ks,kwmap,babichpol,Liguorietal,Liguorietal2}.
On the other hand, the four-point correlation function could display interesting features. From the observational point of view it has been argued that it
could be even more sensible to primordial non-Gaussianity than the bispectrum for very small angular scales in the next CMB
experiments~\cite{OkamotoHu,KogoKomatsu}.
>From the theoretical point of view the four-point connected correlation function for the large scale CMB anisotropies has been
computed in Ref.~\cite{Bartolo:2005fp} giving the theoretical predictions for the quadratic and cubic non-linearities
(which include generic non-Gaussian initial conditions), while Refs.~\cite{SL1,SL2} have computed
the trispectrum from single-field and uncoupled multiple fields in slow-roll inflation showing that it is of the order of the slow-roll parameters.
Interestingly in some cases, such as in some configurations of the curvaton scenarios~\cite{SasakiVW,BSW},
the main source of a non-Gaussian signal can come from the four-point correlations functions,
the three-point correlation function being suppressed (see also the phenomenological model discussed in Ref.~\cite{kwmap3} and Ref.~\cite{HuangShiu}).

However most of these computations deal with the trispectrum of the uniform density curvature perturbation $\zeta$
on large scales within a given inflationary model, but this is
not the physical quantity that is observed. The same is true of course for the bispectrum of the curvature perturbation.
But, if for the case of the curvature perturbation bispectrum one is sure
that the information about the primordial non-Gaussianity in the final observable quantity (the CMB anisotropies) is just the one obtained by evolving linearly
the curvature bispectrum, the same is not true any more for the trispectrum. Let us consider the standard way to characterize the level of non-Gaussianity
in the
gravitational potential by expanding it as
\begin{equation}
\label{fg}
\Phi=\Phi_L+f_{\rm NL}(\Phi^2_L-\langle \Phi^2_L\rangle)+g_{\rm NL} \Phi_L^3\, ,
\end{equation}
where $\Phi_L$ is the linear first-order Gaussian part,
$f_{\rm NL}$ and $g_{\rm NL}$ are the parameters which measure the quadratic and cubic non-linearities. In fact the non-linearity parameters might have
in general
a non trivial scale dependence. At linear order
the curvature perturbation is $\zeta^{(1)}=5 \Phi_L/3$ during the matter dominated epoch, and, for example, the Sachs-Wolfe effect tells us
that $\Delta T/T=-\Phi_L/3$.
The standard approach followed in the literature is to extend both of these two relations also at higher-order, so that for example for quadratic
non-Gaussianity
one writes $\zeta=5 \Phi/3$ using the expansion~(\ref{fg}) up to second-order. Evolving linearly the perturbations in such a way just accounts
for the primordial content of non-Gaussianity and is justified only when one assumes that the primordial level of non-Gaussianity
is much larger of the second-order corrections (in the form of first-order squared perturbations) which arise both in the relation between the curvature
perturbation
and the gravitational potential and when computing the CMB anisotropies. One is guaranteed that such primordial
contribution in the final CMB anisotropies will appear as computed in this way because it is already an intrinsically second-order non-linearity and
all the transfer
functions will therefore be the same as at linear order. For the CMB bispectrum this is also all what is necessary to account for the primordial
non-Gaussianity.
For the trispectrum the situation is different, however. Evolving linearly the primordial cubic
non-Gaussianity in Eq.~(\ref{fg}) would just lead to $\zeta^{(3)}/6=5g_{\rm NL}/3 \Phi_L^{3}$
(splitting $\zeta=\zeta^{(1)}+\zeta^{(2)}/2+...$ into a first- and higher order parts). Since the curvature perturbation remains constant on superhorizon
scales
(for adiabatic perturbations) it is a useful quantity to keep track of the primordial non-Gaussianity. Therefore, suppose to parametrize the primordial
non-Gaussianity as
\begin{equation}
\zeta=\zeta^{(1)}+ (a_{\rm NL}-1) \left( \1{\zeta} \right)^2 +(b_{\rm NL}-1) \left( \1{\zeta} \right)^3\, ,
\end{equation}
where the two non-linearity parameters depend on the physics of a given scenario for the generation of the perturbations (for example, for standard
single-field models of inflation $a_{\rm NL}=1$ and $b_{\rm NL}=1$, plus tiny corrections proportional to the slow-roll parameters). Thus the relation
$\zeta^{(3)}/6=5g_{\rm NL}/3 \Phi_L^{3}$ is equivalent to
\begin{equation}
\label{gr}
g_{\rm NL}=\frac{25}{9} (b_{\rm NL}-1)\, .
\end{equation}
However, this does not catch at all the whole information about the primordial non-Gaussianity that is contained in the CMB anisotropies up to third order.
This is due to two reasons: first, as a source for the evolution of the gravitational potentials at third-order now there are also the gravitational
potentials at second-order which inevitably contain primordial non-linearities proportional to the $f_{\rm NL}$ (or $a_{\rm NL}$) parameter;
and second, when computing the CMB anisotropies at third-order,
one still gets additional terms proportional to $f_{\rm NL} \Phi_L^3$.
This point is clear in the results we present in Sec.~\ref{3S} and~\ref{e3p}.
See for example the first line of Eq.~(\ref{T3scalars}), or Eqs.~(\ref{source}) and~(\ref{A}).
Another example can be found in Ref.~\cite{Bartolo:2005fp}, where a fully non-linear
expression for the Sachs-Wolfe effect has been obtained for generic non-Gaussian initial conditions. The corresponding parameter $g_{\rm NL}$ for CMB
anisotropies
has been computed (see Eq.~(67) of that Reference) by accounting for all the possible dependence on the primordial non-Gaussianity finding an expression
of the type
\begin{equation}
\label{gT}
g_{\rm NL}= \frac{25}{9} (b_{\rm NL}-1)+A({\bf k}_1,{\bf k}_2,{\bf k}_3;a_{\rm NL})\, ,
\end{equation}
which, contrary to Eq.~(\ref{gr}), depends also on the quadratic non-linearity parameter
(here ${\bf k}_i$ are the perturbation wavenumbers in Fourier space).
A similar expression to Eq.~(\ref{gT}) is given for the gravitational potentials at third-order in Eq.~(\ref{Phiba}). These examples show that even in the case
of a ``local'' model for quadratic non-Gaussianity (i.e. a constant $f_{\rm NL}$ parameter), at higher orders such non-linearities are modulated by
first-order perturbations generating a scale-dependent non-Gaussianity.

The best limits up to date on the $f_{\rm NL}$ parameter come from measurements of the CMB bispectrum on the WMAP data
giving $-36 \leq f_{\rm NL} \leq  100$~\cite{kwmap,kwmap3,Crem},
while there is at present no real bound on $g_{\rm NL}$. Given the increasing precision of future mesurements of CMB anisotropies,
it is clear that it is of fundamental importance to provide accurate predictions for all the cubic non-linearities that enter in the
evolution of the cosmological perturbations.
This is not only mandatory to be able to evaluate the trispectrum of CMB anisotropies,
but becomes crucial if one wants to account for the precise dependence of the trispectrum on the primordial non-Gaussianity.
Spurred by these considerations, in this paper we will focus on the CMB anisotropies up to
third-order from gravitational perturbations due to the redshift the photons suffer when they travel from the last scattering surface to the observer.
The computation will be performed following  Refs.~\cite{Pyne:1993np,Pyne:1995bs,Mollerach:1997up}
by perturbing at the desired order the photon geodesic
accounting for scalar, vector and tensor perturbations.
However we will also propose an alternative method to solve the photon geodesic equation (for scalar perturbations), which consists in
a non-perturbative computation. In the last part of the paper we focus on the large-scale approximation and on scalar perturbations for a Universe filled
with non-relativistic matter and a non-vanishing cosmological constant.
In fact one of the goals is to characterize the evolution of the gravitational potentials on large scales which is
responsible for the (late) Integrated Sachs-Wolfe effect, thus completing the analysis for the large scale CMB anisotropies at third-order already started in
Ref.~\cite{Bartolo:2005fp}.
Of course it is out of the goals of this paper to deal with the anisotropies generated from the dynamics at recombination
and due to the scattering terms~\footnote{A complete computation of the CMB anisotropies due to the non-linear dynamics
taking place at the last scattering epoch has been performed in Refs.~\cite{CMB2first,CMB2second},
where the full system of the Boltzmann equations up to second-order for photons, baryons and cold dark matter have been
presented togheter with analytical solutions in the tight coupling approximation.}.

\section{Temperature anisotropies}
\label{2S}
We are interested in the pattern of fluctuations of the CMB temperature as measured by an observer
in a perturbed flat Friedmann-Robertson-Walker spacetime.

The line element can be written as
\begin{equation}
 \label{metric}
 ds^2 = a^2(\eta) g_{\mu \nu} dx^\mu dx^\nu = a^2(\eta) \left[ \0{g_{\mu \nu}} + \1{g_{\mu \nu}} + \2{g_{\mu \nu}} + \3{g_{\mu \nu}} + \dots \right]
 	dx^\mu dx^\nu \, ,
\end{equation}
where $a(\eta)$ is the scale factor in conformal time $d\eta=dt/a(t)$, $\0{g_{\mu \nu}}$ is the background Minkowski metric,
and $g^{(r)}_{\mu \nu}$ is the $r$-th order metric perturbation.

Photons travel along null geodesics $x^\mu(\lambda)$,
where $\lambda$ is an affine parameter in the conformal metric $g_{\mu \nu}$.
The photon path connects the point of observation, with coordinates $x^\mu_\O = (\eta_\O, \mathbf{x}_\O)$,
to the hypersurface of emission, defined as the spacelike hypersurface of constant conformal time $\eta_{\mathcal E}$.
The actual last scattering surface is the intersection of the observer's past light cone with this hypersurface.
We assume that on the constant $\eta_{\mathcal E}$ hypersurface every point $\mathbf{p}$ emits thermal radiation,
characterized by a temperature $T(\mathbf{p},\d)$ which depends on the point of emission $\mathbf{p}$
and on the direction of emission, described by the vector $\d$ normalized to unity in the conformal background metric.
The different photon paths are specified by the direction from which they arrive at $\O$,
described by a vector $\mathbf{\hat{e}}$, normalized to unity in the conformal background metric;
this unit vector can be thought as the direction toward which the observer is pointing an antenna.
The initial conditions $x^\mu_\O$, $\mathbf{\hat{e}}$ determine the point $\mathbf{p}$ and direction $\d$ of emission.

During their travel from the last scattering surface to the observer the CMB photons suffer a redshift determined by the ratio of the emitted
frequency $\omega_{{\mathcal E}}$ and the observed one $\omega_{\mathcal O}$. For a blackbody spectrum, the ratio $\omega / T$ is constant along the photon path,
and the temperature measured by an observer is given by
\begin{equation}
 \label{blackbody}
 T_\O (\mathbf{x}_\O, \mathbf{\hat{e}}) = \frac{\omega_\O}{\omega_{\mathcal E}} T_{\mathcal E} (\mathbf{p}, \d) \, .
\end{equation}
The expression for the frequency is
\begin{equation}
 \label{frequency}
 \omega = - g_{\mu \nu} u^\mu k^\nu \, ,
\end{equation}
where $u^\mu$ is the four-velocity of the observer or emitter, normalized to $a^2 g_{\mu \nu} u^\mu u^\nu = -1$,
and $k^\nu(\lambda) = \frac{d x^\nu}{d \lambda}$ is the photon wavevector, tangent to the null geodesic $x^\nu(\lambda)$.

Given the initial conditions $x^\mu_\O, \mathbf{\hat{e}}, \omega_\O$, we need to compute the quantities $\mathbf{p}, \d, \omega_{\mathcal E}$,
propagating the photons back from the observation to the emission surface.
These will depend on the photon path and the associated wavevector,
which we expand in series of the metric perturbations $g_{\mu \nu}^{(r)}$ and their derivatives:
\begin{align}
 x^\mu(\lambda) =& x^{(0)\mu}(\lambda) + x^{(1)\mu}(\lambda) + x^{(2)\mu}(\lambda) + x^{(3)\mu}(\lambda) + \dots \\
 k^\mu(\lambda) =& k^{(0)\mu}(\lambda) + k^{(1)\mu}(\lambda) + k^{(2)\mu}(\lambda) + k^{(3)\mu}(\lambda) + \dots
\end{align}

We find it useful to write the perturbed (conformal) metric as
\begin{equation}
\label{metricdef}
 g_{\mu \nu} dx^\mu dx^\nu = - e^{2 \Phi} d\eta^2 + 2 \omega_i d\eta dx^i
 	+ \left( e^{-2 \Psi} \delta_{ij} + \chi_{ij} \right) dx^i dx^j\, ,
\end{equation}
which is valid at any order in perturbation theory. The two gravitational potentials $\Phi$ and $\Psi$ correspond to scalar metric perturbations,
$\omega_i$ includes a scalar and a vector perturbation, and $\chi_{ij}$ contains another scalar, as well as vector and tensor perturbations.
To make contact with the usual perturbative calculations one expands every quantity as
$\Phi = \1{\Phi} + \frac{1}{2} \2{\Phi} + \frac{1}{6} \3{\Phi} + \dots$ and analogously for the others.
Notice that, unlike $\omega^{(r)}_i$ and $\chi^{(r)}_{ij}$, the quantities $\Phi^{(r)}$ and $\Psi^{(r)}$
are not the usual $r$-th order scalars $\phi^{(r)}$, $\psi^{(r)}$ appearing in Refs.~\cite{MMB,review}.
However, it is immediate to find the relation between these quantities. Up to third order
$\1{\Phi} = \1{\phi}$, $\1{\Psi} = \1{\psi}$, $\2{\Phi} = \2{\phi} - 2 \phi_{(1)}^2$,
$\2{\Psi} = \2{\psi} + 2 \psi_{(1)}^2$, $\3{\Phi} = \3{\phi} - 6 \1{\phi} \2{\phi} + 8 \phi_{(1)}^3$,
$\3{\Psi} = \3{\psi} + 6 \1{\psi} \2{\psi} + 8 \psi_{(1)}^3$. The form of the metric~(\ref{metricdef}) greatly helps in the intermediate computations,
however some results will be expressed in the variables $\phi^{(r)}$ and $\psi^{(r)}$ because they appear in a more compact form. Notice that in this
section we will not choose a particular gauge so all the following expressions for the CMB anisotropies are valid in any gauge.

The four-velocity will be expanded as
\begin{equation}
 u^\mu = \frac{1}{a} \left( \delta_0^\mu + v^{(1)\mu} + \frac{1}{2} v^{(2)\mu} + \frac{1}{6} v^{(3)\mu} + \dots \right)\, .
\end{equation}

The zero component of the velocity is fixed from the normalization condition; we find
\begin{align}
 v^0_{(1)} =& - \1{\phi} \, ; \\
 v^0_{(2)} =& - \2{\phi} + 3 \left( \1{\phi} \right)^2
  		+ 2 \1{\omega_i} v^i_{(1)} + v^i_{(1)} \1{v_i} \, ; \\
 v^0_{(3)} =& - \3{\phi} + 3 \2{\omega_i} v^i_{(1)} + 3 \1{\omega_i} v^i_{(2)} + 9 \2{\phi} \1{\phi}
  	+ 3 \left( \1{\chi_{ij}} - 2 \1{\psi} \delta_{ij} \right) v^i_{(1)} v^j_{(1)} \nonumber \\
  	&- 12 \1{\phi} \1{\omega_i} v^i_{(1)} - 15 \left( \1{\phi} \right)^3 - 3 \1{\phi} v^i_{(1)} \1{v_i}
  	+ 3 v^i_{(2)} \1{v_i} \, .
\end{align}

In order to compute the observed temperature up to third order in perturbations,
we expand the frequency as
\begin{equation}
 \omega = \0{\omega} \left( 1 + \1{\tilde{\omega}} + \2{\tilde{\omega}} + \3{\tilde{\omega}} + \dots \right) \, ,
\end{equation}
and the temperature at emission as
\begin{equation}
 T_{\mathcal E}(\mathbf{p}, \d) = \0{T_{\mathcal E}} \left( 1 + \tau(\mathbf{p}, \d) \right) \, ,
\end{equation}
where $\tau = \1{\tau} + \2{\tau} + \3{\tau} + \dots$ is the intrinsic temperature fluctuation at emission.
A calculation of this quantity up to third-order is beyond the goal of this paper, since
we are interested in the additional effects of gravity along the photon path.
However on large scales (bigger than the horizon at recombination) it has been computed in a non-pertubative way in Ref.~\cite{Bartolo:2005fp}.
Its fully non-linear expression is very simple: $\tau=e^{-2\Phi/3}-1$. Moreover, on smaller scales, a complete treatment of the CMB anisotropies
including the acoustic oscillations on the surface of last scattering up to second-order have been studied analytically
in Ref.~\cite{CMB2first,CMB2second}, and a full numerical evaluation
can be performed Ref.~\cite{numerico} using the set of Boltzmann equations provided in Ref.~\cite{CMB2first}.

Now, we have to take into account that we need to expand the point $\mathbf{p}$ and direction $\d$ at emission as
$\mathbf{p} = \0{\mathbf{p}} + \1{\mathbf{p}} + \2{\mathbf{p}} + \dots$ and $\d = \0{\d} + \1{\d} + \2{\d} + \dots$;
then we Taylor expand $\tau(\mathbf{p}, \d)$ around the background values, obtaining
\begin{equation}
\begin{split}
 \tau(\mathbf{p}, \d) =& \1{\tau} +
 	\left[ \2{\tau} + p_{(1)}^i \frac{\partial \1{\tau}}{\partial x^i} + d_{(1)}^i \frac{\partial \1{\tau}}{\partial d^i} \right]\\
 	&+ \left[ \3{\tau} + p_{(1)}^i \frac{\partial \2{\tau}}{\partial x^i} + d_{(1)}^i \frac{\partial \2{\tau}}{\partial d^i}
 		+ \frac{1}{2} p_{(1)}^i p_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial x^i \, \partial x^j}
 		+ p_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial x^i \, \partial d^j} \right. \\
 		& \left. \quad + \frac{1}{2} d_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial d^i \, \partial d^j}
 			+ p_{(2)}^i \frac{\partial \1{\tau}}{\partial x^i} + d_{(2)}^i \frac{\partial \1{\tau}}{\partial d^i} \right]
\end{split}
\end{equation}
where it is understood that all quantities are evaluated at $\0{\mathbf{p}}$, $\0{\d}$.
We find $p^{(1)i} = x^{(1)i} - k^{(0)i} x^{(1)0}$,
$p^{(2)i} = x^{(2)i} - k^{(0)i} x^{(2)0} + k^{(0)i} k^{(1)0} x^{(1)0} - k^{(1)i} x^{(1)0}$,
$d^{(1)i} = \frac{k^{(0)i} + k^{(1)i}}{\left|k^{(0)i} + k^{(1)i} \right|}
- \frac{k^{(0)i}}{\left|k^{(0)i}\right|}$,
$d^{(2)i} = \frac{k^{(0)i} + k^{(1)i} + k^{(2)i}}{\left|k^{(0)i} + k^{(1)i} + k^{(2)i} \right|}
- \frac{k^{(0)i} + k^{(1)i}}{\left|k^{(0)i} + k^{(1)i} \right|}$.
Notice that therefore the quantities $\hat{\bf d}^{(i)}$ are just the
difference between the photon normalized wavevector up to the $i$-th order and that at the $(i-1)$-th order.

Performing all the expansions in Eq.\eqref{blackbody}
and dividing by the background temperature $\frac{a_{\mathcal E}}{a_\O} \0{T_{\mathcal E}}$, we find
\begin{equation}
 \frac{\delta \1{T}}{T} = \1{\tilde{\omega}}_\O - \1{\tilde{\omega}}_{\mathcal E} + \1{\tau} \, ;
\end{equation}
\begin{equation}
 \frac{\delta \2{T}}{T} = \2{\tilde{\omega}}_\O - \2{\tilde{\omega}}_{\mathcal E}
 	+ \2{\tau} + p_{(1)}^i \frac{\partial \1{\tau}}{\partial x^i} + d_{(1)}^i \frac{\partial \1{\tau}}{\partial d^i}
 	+ \left( \1{\tau} - \1{\tilde{\omega}}_{\mathcal E} \right) \left( \1{\tilde{\omega}}_\O - \1{\tilde{\omega}}_{\mathcal E} \right) \, ;
\end{equation}
\begin{equation}
\begin{split}
 \frac{\delta \3{T}}{T} =& \, \3{\tilde{\omega}}_\O - \3{\tilde{\omega}}_{\mathcal E}
 	+ \3{\tau} + p_{(1)}^i \frac{\partial \2{\tau}}{\partial x^i} + d_{(1)}^i \frac{\partial \2{\tau}}{\partial d^i}
 	+ \frac{1}{2} p_{(1)}^i p_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial x^i \, \partial x^j} \\
 	& + p_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial x^i \, \partial d^j}
 	+ \frac{1}{2} d_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial d^i \, \partial d^j}
 		+ p_{(2)}^i \frac{\partial \1{\tau}}{\partial x^i} + d_{(2)}^i \frac{\partial \1{\tau}}{\partial d^i} \\
 	& + \left( \2{\tau} + p_{(1)}^i \frac{\partial \1{\tau}}{\partial x^i}
 		+ d_{(1)}^i \frac{\partial \1{\tau}}{\partial d^i} - \2{\tilde{\omega}}_{\mathcal E} \right)
 		\left( \1{\tilde{\omega}}_\O - \1{\tilde{\omega}}_{\mathcal E} \right) \\
 	& + \left( \1{\tau} - \1{\tilde{\omega}}_{\mathcal E} \right)
 		\left( \2{\tilde{\omega}}_\O - \2{\tilde{\omega}}_{\mathcal E} + (\1{\tilde{\omega}}_{\mathcal E})^2
 		- \1{\tilde{\omega}}_{\mathcal E} \1{\tilde{\omega}}_\O \right) \, .
\end{split}
\end{equation}

The expansion of the frequency using Eq.(\ref{frequency}) yields:
\begin{equation}
 \1{\tilde{\omega}_{\mathcal E}} = k^0_{1} + \1{\phi} + \1{\omega_i} e^i + \1{v_i} e^i \, ;
\end{equation}

\begin{equation}
\begin{split}
 \2{\tilde{\omega}_{\mathcal E}} =& k^0_{(2)} + \frac{1}{2} \2{\phi} + \frac{1}{2} \2{\omega_i} e^i
  	 + \frac{1}{2} v^i_{(1)} \1{v_i} - \frac{1}{2} (\1{\phi})^2
  	 + \1{\phi} k^0_{(1)} - \1{\omega_i} k^i_{(1)} - \1{v_i} k^i_{(1)} \\
  	 &- \1{\phi} \1{\omega_i} e^i - 2 \1{\psi} \1{v_i} e^i
  	 + \1{\chi_{ij}} v^i_{(1)} e^j + \frac12 \2{v_i} e^i \\
  	 &+ \1{\Delta \lambda} \frac{d k^0_{(1)}}{d \lambda}
  	 + p^i_{(1)} \left( \partial_i \1{\phi} + \partial_i \1{\omega_j} e^j + \partial_i \1{v_j} e^j \right) \, ;
\end{split}
\end{equation}

\begin{equation}
\begin{split}
 \3{\tilde{\omega}_{\mathcal E}} =& k^0_{(3)} + \1{\Delta \lambda} \frac{d k^0_{(2)}}{d \lambda}
  		+ \2{\Delta \lambda} \frac{d k^0_{(1)}}{d \lambda}
  		+ \frac{1}{2} \left( \1{\Delta \lambda} \right)^2 \frac{d^2 k^0_{(1)}}{d \lambda^2} \\
  &+ \frac16 \3{\phi} + \frac12 p^i_{(1)} \partial_i \2{\phi}
  		+ p^i_{(2)} \partial_i \1{\phi} + \frac12 p^i_{(1)} p^j_{(1)} \partial_i \partial_j \1{\phi} \\
  &+ e^i \left( \frac16 \3{v_i} + \frac12 p^j_{(1)} \partial_j \2{v_i}
  		+ p^j_{(2)} \partial_j \1{v_i}
  		+ \frac{1}{2} p^j_{(1)} p^k_{(1)} \partial_j \partial_k \1{v_i} \right) \\
  &- \frac12 \2{\phi} \1{\phi}
  	+ \frac12 \left( \1{\chi_{ij}} - 2 \1{\psi} \delta_{ij} \right) v^i_{(1)} v^j_{(1)} \\
  &+ \frac12 \left( \1{\phi} \right)^3 + \frac12 v^i_{(1)} \2{v_i}
  	+ \frac12 \1{\phi} \1{v_i} v^i_{(1)}
  	- \1{\phi} p^i_{(1)} \partial_i \1{\phi} \\
  &+ \frac12 k^0_{(1)} \left[ \2{\phi} + 2 p^i_{(1)} \partial_i \1{\phi}
  	- \left( \1{\phi} \right)^2 + v^i_{(1)} \1{v_i} \right]
  	+ \1{\phi} \left( k^0_{(2)} + \1{\Delta \lambda} \frac{d k^0_{(1)}}{d \lambda} \right) \\
  & + \frac16 \3{\omega_i} e^i + \frac12 p^i_{(1)} \partial_i \2{\omega_j} e^j
  		+ p^i_{(2)} \partial_i \1{\omega_j} e^j
  		+ \frac{1}{2} p^i_{(1)} p^j_{(1)} \partial_i \partial_j \1{\omega_k} e^k \\
  &- \left( k^i_{(1)} + e^i \1{\phi} \right) \left( \frac12 \2{\omega_i} + p^j_{(1)} \partial_j \1{\omega_i} \right)
  	- \1{\omega_i} \left( k^i_{(2)} + \1{\Delta \lambda} \frac{d k^i_{(1)}}{d \lambda} \right) \\
  &+ \1{\omega_i} \Bigg\{ \1{\phi} k^i_{(1)} + \frac12 e^i \Big[ - \2{\phi} - 2 p^j_{(1)} \partial_j \1{\phi}
  	+ 3 \left( \1{\phi} \right)^2 + 2 \1{\omega_j} v^j_{(1)} + v^j_{(1)} \1{v_j} \Big] \Bigg\} \\
  & - \delta_{ij} k^i_{(1)} \left( \frac12 v^j_{(2)} + p^k_{(1)} \partial_k v^j_{(1)} \right)
  	 	- \delta_{ij} v^i_{(1)} \left( k^j_{(2)} + \1{\Delta \lambda} \frac{d k^j_{(1)}}{d \lambda} \right) \\
  &+ \left( \1{\chi_{ij}} - 2 \1{\psi} \delta_{ij} \right)
  		\left( \frac12 v^i_{(2)} e^j + p^k_{(1)} \partial_k v^i_{(1)} e^j - v^i_{(1)} k^j_{(1)} \right) \\
  &+ \left[ \left( \frac12 \2{\chi_{ij}} - \2{\psi} \delta_{ij} \right)
  	 	+ p^k_{(1)} \partial_k \left( \1{\chi_{ij}} - 2 \1{\psi} \delta_{ij} \right) \right] v^i_{(1)} e^j \, .
\end{split}
\end{equation}
In these expressions, $\1{\Delta \lambda}$ and $\1{\Delta \lambda}+ \2{\Delta \lambda}$ are the differences
in affine parameter between the points where the first order (respectively, the second order)
and the background geodesics intersect the $\eta = \eta_{\mathcal E}$ hypersurface.
They are given by $\1{\Delta \lambda} = - x^{(1)0}$, $\2{\Delta \lambda} = - x^{(2)0} + k^{(1)0} x^{(1)0}$.

The next step is to obtain the null geodesics up to third order,
using the formalism set up in \cite{Pyne:1995bs,Pyne:1993np}. This will allow us to compute at the desired order
the photon path $x^{\mu}(\lambda)$ and the associated wavevector $k^{\mu}(\lambda)$
appearing in the previous expressions for the perturbations of the  photon frequency. The $a$-th order geodesic equation can be recasted
as the forced Jacobi equation~\cite{Pyne:1995bs,Pyne:1993np}
\begin{equation}
 \frac{d x^{(a)\mu}}{d \lambda^2} + 2 \Gamma^{(0)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(a)\beta}
 + \partial_\sigma \Gamma^{(0)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(0)\beta} x^{(a)\sigma} = f^{(a)\mu} \, .
\end{equation}

In a flat background, the solutions are
\begin{align}
\label{xk}
 x^{(a) \mu} =& \left(\lambda - \lambda_\O \right) k^{(a) \mu} (\lambda_\O)
 	+ \int_{\lambda_\O}^\lambda (\lambda - \tilde{\lambda}) f^{(a)\mu} (\tilde{\lambda}) d\tilde{\lambda} \, , \\
 k^{(a) \mu} =& k^{(a) \mu} (\lambda_\O) + \int_{\lambda_\O}^\lambda f^{(a)\mu} (\tilde{\lambda}) d\tilde{\lambda} \, .
\end{align}

The forcing terms, in a flat background, are given by
\begin{align}
\label{fterms}
 f^\mu_{(1)} =& - \Gamma^{(1)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(0)\beta} \, , \\
 f^\mu_{(2)} =& - 2 \Gamma^{(1)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(1)\beta}
 -  x^{(1)\sigma} \partial_\sigma \Gamma^{(1)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(0)\beta}
 - \Gamma^{(2)\mu}_{\alpha \beta} k^{(0)\alpha} k^{(0)\beta} \, , \\
\begin{split}
\label{fterms3}
 f^\mu_{(3)} =& - \Gamma^{(1) \mu}_{\alpha \beta} \left( 2 k^{\alpha}_{(2)} k^{\beta}_{(0)} + k^{\alpha}_{(1)} k^{\beta}_{(1)} \right)
 	- 2 k^{\alpha}_{(1)} k^{\beta}_{(0)} \left( \Gamma^{(2) \mu}_{\alpha \beta}
 	+ \partial_\sigma \Gamma^{(1) \mu}_{\alpha \beta} x^\sigma_{(1)} \right) \\
 	&- k^\alpha_{(0)} k^\beta_{(0)} \left( \Gamma^{(3) \mu}_{\alpha \beta}
 	+ \partial_\sigma \Gamma^{(1) \mu}_{\alpha \beta} x^\sigma_{(2)}
 	+ \partial_\sigma \Gamma^{(2) \mu}_{\alpha \beta} x^\sigma_{(1)}
 	+ \frac12 \partial_\sigma \partial_\tau \Gamma^{(1) \mu}_{\alpha \beta} x^\sigma_{(1)} x^\tau_{(1)} \right) \, .
\end{split}
\end{align}
 The Christoffel symbols are calculated in the Appendix,
where one can find their expressions for the conformal metric \eqref{metricdef}. Notice
that we have  recasted them in a compact form without specifying the order of the perturbations, even tough they can be used just up
to third-order (at least for the terms coming from the vector and tensor perturbations). The method that uses Eqs.~(\ref{fterms}) is an iterative method:
once one knows the geodesic equations of $r-th$ oder, he or she is able to determine those at the next-order through Eqs.~(\ref{xk}).
Thus, in order to have all the means to compute the photon geodesics up to third order in the following,
we report the expressions that we have computed for the wavevectors $k^{(\mu)}$ up to second-order :
\begin{equation}
\label{k10}
 k^0_{(1)} = - 2 \1{\Phi} - \1{\omega_i} e^i + \int_{\lambda_\O}^\lambda d\tilde{\lambda} \1{A}'(\tilde{\lambda})
    = - 2 \1{\Phi} - \1{\omega_i} e^i + I_1(\lambda)\, ,
\end{equation}
\begin{equation}
\label{k1i}
 k^i_{(1)} = - 2 \1{\Psi} e^i - \omega_{(1)}^i + \chi_{(1)}^{ij} e_j
    - \int_{\lambda_\O}^\lambda d\tilde{\lambda} \partial^i \1{A}(\tilde{\lambda})
    = - 2 \1{\Psi} e^i - \omega_{(1)}^i + \chi_{(1)}^{ij} e_j - I_1^i(\lambda) \, ,
\end{equation}
for the first-order wavevectors, and
\begin{equation}
\begin{split}
\label{k20}
 k^0_{(2)} =& - \2{\Phi} - 2 \1{\Phi}^2 - \frac12 \2{\omega_i} e^i - 2 k^0_{(1)} \1{\Phi}
    - 2 x_{(1)}^\mu \partial_\mu \1{\Phi} + \1{\omega_i} k^i_{(1)} - x_{(1)}^\mu \partial_\mu \1{\omega_i} e^i \\
    &+ \int_{\lambda_\O}^\lambda d\tilde{\lambda} \Bigg[ \frac12 \2{A}'
        + 2 \left( \1{\Phi} \1{\Phi}' - \1{\Psi} \1{\Psi}'\right)
        + \left( \1{\chi}'_{ij} e^j - \1{\omega_i}' \right) \left( e^i k^0_{(1)} + k^i_{(1)} \right) \\
    &\qquad \qquad + 2 k_{(1)}^0 \1{A}' + 2 \1{\Psi}' \1{A} + x_{(1)}^\mu \partial_\mu \1{A}' \Bigg]\, ,
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\label{k2i}
 k^i_{(2)} =& - e^i \2{\Psi} + 2 e^i \1{\Psi}^2 - 2 e^i \1{x}^\mu \partial_\mu \1{\Psi} + 2 \1{k}^i \1{\Psi} \\
    &+ \frac12 \chi^i_{(2)j} e^j + x^\mu \partial_\mu \chi^i_{(1)j} e^j - \chi^i_{(1)j} \1{k}^j
    - \frac12 \omega_{(2)}^i - x_{(1)}^\mu \partial_\mu \omega_{(1)}^i - k_{(1)}^0 \omega_{(1)}^i \\
    &+ \int_{\lambda_\O}^\lambda d\tilde{\lambda} \Bigg[- \frac12 \partial_i \2{A} - x^\mu_{(1)} \partial_\mu \partial^i \1{A}
    - 2 \1{k}^0 \partial^i \1{\Phi} + \partial^i \left( \1{\Psi}^2 - \1{\Phi}^2 \right) \\
    &\qquad \qquad + 2 e_j \1{k}^j \partial^i \1{\Psi} - e^j \1{k}^k \partial^i \1{\chi}_{jk}
    + k_{(1)}^j \partial^i \omega^{(1)}_j - k_{(1)}^0 \partial^i \1{\omega_j} e^j \Bigg] \, ,
\end{split}
\end{equation}
at second-order, where
\begin{equation}
 A^{(r)} \equiv \Phi^{(r)} + \Psi^{(r)} + \omega^{(r)}_i e^i - \frac12 \chi^{(r)}_{ij} e^i e^j \, .
\end{equation}

As a non-trivial test for the correctness of the above expressions, we have explicitly verified that the first
and second order wavevectors do satisfy the null vector condition $k^\mu k_\mu = 0$.
When expanded perturbatively, this reads at linear order
\begin{equation}
  k_{(1)}^0 + e_i k_{(1)}^i = - A^{(1)} \, ,
\end{equation}
and at second order
\begin{equation}
\begin{split}
 &k_{(2)}^0 + e_i k_{(2)}^i + \frac{1}{2} (k_{(1)}^0)^2 - \frac{1}{2} k_{(1)}^i k^{(1)}_i + 2 \Phi_{(1)} k_{(1)}^0 - \omega^{(1)}_i k_{(1)}^i
 	+ \omega^{(1)}_i e^i k_{(1)}^0 - 2 \Psi_{(1)} e_i k_{(1)}^i + \chi^{(1)}_{ij} e^i k_{(1)}^j \\
 	&+ \frac{1}{2} A_{(2)} + \Phi_{(1)}^2 - \Psi_{(1)}^2 + x_{(1)}^\mu \partial_\mu A_{(1)}= 0 \, .
\end{split}
\end{equation}
Both equations are satisfied when one substitutes the expressions~\eqref{k10}-\eqref{k2i}.

\section{CMB anisotropies on large scales from scalar perturbations}
\label{3S}

In this Section, we focus on scalar perturbations and we give the expressions for the CMB anisotropies on large scales.
Two main results are achieved. First we are able to give, using a non-perturbative method, an integral solution for the photon
geodesic equations which holds at any order in perturbation theory. Second,
we solve for the first time the evolution of the third-order gravitational potentials in a universe filled with non-relativistic
matter and a cosmological constant. The presence of a non-vanishing cosmological constant makes the gravitational potentials vary at late
times producing on large scales an Integrated-Sachs-Wolfe effect which is the counterpart of what happens
at linear order. Of course, due to the non linear evolution of the perturbations, other integrated terms will appear in the form of cubic
corrections. The ISW constitutes the main effect for large scale CMB anisotropies together with the Sachs-Wolfe effect originating at the last
scattering surface. In fact both of them keep memory of the initial primordial non-Gaussianity of the cosmological perturbations
(see details in Sec.~\ref{IC}), and as such  our expressions can be of particular relevance when trying to pin down the primordial
non-Gaussian content from  the higher-order statistics of the CMB anisotropies as the bispectrum or the trispectrum.
Let us see in details these two results.

 In the limit of large scales, it is convenient to adopt a non-perturbative formalism
to get the temperature anisotropies and the Einstein equations.
We take the metric to be $ds^2 = a^2(\eta)[- e^{2 \Phi} d\eta^2 + e^{-2 \Psi} \delta_{ij} dx^i dx^j]$
and for simplicity we take the points of emission and observation to be comoving
(in such a way we lose the Doppler effects, however they are important on scales smaller than those we are interested here). In writing this
metric we have neglected vector and tensor perturbation modes. For the vector perturbations the reason is that we are interested in
long-wavelength perturbations, while vector modes will 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 mechanism for comsological
perturbations, as inflation). For example the results of Ref.~\cite{MHMpol} show clearly this for second-order perturbations.
In order to study the CMB anisotropies from scalar perturbations in the large scale limit,
the tensor contribution can be negleted, since on large scales it has been proven to remain constant and to give a negligible
effect being of
the order of (powers of) the slow-roll parameters during inflation~\cite{SalopekBond,maldacena}.

We start from eq. \eqref{blackbody} and we write the frequency as
\begin{equation}
 \omega = - g_{\mu \nu} u^\mu k^\nu = e^{2 \Phi} u^0 k^0 = e^{\Phi} k^0 \, ,
\end{equation}
where $u^0=e^{-\Phi}$ is given by the normalization condition. The observed temperature is given by
\begin{equation}
 T_\O (\mathbf{x},\mathbf{\hat{e}}) = \frac{a_{\mathcal E}}{a_\O} e^{\Phi_\O-\Phi_{\mathcal E}} \frac{k^0_\O}{k^0_{\mathcal E}} T_{\mathcal E} (\mathbf{p},\d) \, .
\end{equation}

We make another simplification: namely, we reabsorb the terms computed at the point of observation
into a redefinition of $T_\O$.
This amounts to lose the monopole term, which however is unobservable.

At this point, we can perform the double expansion in perturbation orders and around the background geodesic
(i.e., the line of sight). We find
\begin{equation}
 \frac{\1{\delta} T}{T} = \1{\tau} - \1{\Phi_{\mathcal E}} - k^{(1)0}_{\mathcal E} \, ;
\end{equation}

\begin{equation}
\label{T2scalars}
\begin{split}
 \frac{\2{\delta} T}{T} =& \2{\tau} + p_{(1)}^i \frac{\partial \1{\tau}}{\partial x^i} + d_{(1)}^i \frac{\partial \1{\tau}}{\partial d^i}
 		- \frac12 \2{\Phi_{\mathcal E}} - p_{(1)}^i \frac{\partial \1{\Phi_{\mathcal E}}}{\partial x^i} + \frac{1}{2} (\1{\Phi_{\mathcal E}})^2 \\
 	& - k^{(2)0}_{\mathcal E} - \1{\Delta \lambda} \frac{d k^{(1)0}_{\mathcal E}}{d \lambda} + (k^{(1)0}_{\mathcal E})^2
 		- \1{\tau} \1{\Phi_{\mathcal E}} - \1{\tau} k^{(1)0}_{\mathcal E} + \1{\Phi_{\mathcal E}} k^{(1)0}_{\mathcal E} \, ;
\end{split}
\end{equation}

\begin{equation}
\label{T3scalars}
\begin{split}
 &\frac{\3{\delta} T}{T} = \3{\tau} - k^{(3)0}_{\mathcal E} - \frac{1}{6} \3{\Phi_{\mathcal E}}
 	+ \frac{1}{2} \1{\Phi_{\mathcal E}} \2{\Phi_{\mathcal E}} - \frac{1}{6} \left( \1{\Phi_{\mathcal E}} \right)^3
 	+ 2 k^{(1)0}_{\mathcal E} k^{(2)0}_{\mathcal E} - (k^{(1)0}_{\mathcal E})^3 - \2{\tau} \1{\Phi} \\
 	&+ \left[ k^{(2)0}_{\mathcal E} - (k^{(1)0}_{\mathcal E})^2 \right] \left(\1{\Phi}- \1{\tau} \right)
 	+ \frac{1}{2} \1{\tau} \left[ - \2{\Phi_{\mathcal E}} + \left( \1{\Phi_{\mathcal E}} \right)^2 \right] \\
 	&+ k^{(1)0} \left[ - \2{\tau} + \frac{1}{2} \2{\Phi_{\mathcal E}} - \frac{1}{2} \left( \1{\Phi_{\mathcal E}} \right)^2
 		+ \1{\tau} \1{\Phi} \right] \\
 	&+ p_{(1)}^i \left[ \partial_i \2{\tau} - \frac{1}{2} \partial_i \2{\Phi} - \1{\tau} \partial_i \1{\Phi}
 		- \1{\Phi} \partial_i \1{\tau}
 		+ k^{(1)0} \left( \partial_i \1{\Phi}- \partial_i \1{\tau} \right)+ \1{\Phi} \partial_i \1{\Phi} \right] \\
 	&+ d_{(1)}^i \left[ \frac{\partial \2{\tau}}{\partial d^i}
 		- \left( k^{(1)0} + \1{\Phi} \right) \frac{\partial \1{\tau}}{\partial d^i} \right]
 	+ \frac{1}{2} p_{(1)}^i p_{(1)}^j \left( \partial_i \partial_j \1{\tau} - \partial_i \partial_j \1{\Phi} \right) \\
 	&+ \frac{1}{2} d_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial d^i \, \partial d^j}
 	+ p_{(2)}^i \left( \partial_i \1{\tau} - \partial_i \1{\Phi} \right)
 	+ d_{(2)}^i \frac{\partial \1{\tau}}{\partial d^i} + p_{(1)}^i d_{(1)}^j \frac{\partial^2 \tau_{(1)}}{\partial x^i \partial d^j} \\
 	&+ \1{\Delta \lambda} \left[ - \frac{d k^{(2)0}_{\mathcal E}}{d \lambda}
 		- \frac{1}{2} \1{\Delta \lambda} \frac{d^2 k^{(1)0}_{\mathcal E}}{d \lambda^2}
 		+ \left( 2 k^{(1)0}_{\mathcal E} + \1{\Phi} - \1{\tau} \right) \frac{d k^{(1)0}_{\mathcal E}}{d \lambda}\right]
 	- \2{\Delta \lambda} \frac{d k^{(1)0}_{\mathcal E}}{d \lambda} \, .
\end{split}
\end{equation}

In order to compute the photon wavevectors and the null geodesics, we will proceed in a different way than the usual standard perturbative method
introduced in \cite{Pyne:1995bs,Pyne:1993np} which makes use of the Eqs.~(\ref{fterms})-(\ref{fterms3}).
Let us instead start from a fully non-linear geodesic equation (time component) which can be written as
\begin{equation}
\label{geonp}
 \frac{d k^0}{d \lambda} = - 2 \frac{d \Phi}{d \lambda} k^0 + \left( \Phi' + \Psi' \right) \left(k^0\right)^2
\end{equation}
where we have used the null vector condition $g_{\mu \nu} k^\mu k^\nu = 0$. We can write the formal solution (formal since the potentials depend on the true geodesic)
as
\begin{equation}
\label{formalsol}
 k^0(\lambda) = e^{-2\Phi(\lambda)} \left[1
 	- \int_{\lambda_\O}^\lambda e^{-2\Phi(x(\t{\lambda}))}
 	\left( \Phi'(x(\t{\lambda})) + \Psi'(x(\t{\lambda})) \right) d\t{\lambda} \right]^{-1} \, .
\end{equation}

The advantage of Eq.~\eqref{formalsol} is that it can be straightforwardly expanded to find
the perturbations to the background wavevector (and, by integration, to the background geodesic),
up to any desired order. To make contact with the standard perturbative results, in addition to the expansion in
perturbative orders we need to perform a Taylor expansion of the true geodesic around the background geodesic
for the potentials, in order to express all the quantities on the background geodesic.
For example, if we have the non-linear quantity $\alpha(x)$ computed on the true
geodesic, its expansion up to second order will be
$\alpha(x) = \alpha^{(1)} (x^{(0)}) + \partial_\mu \alpha^{(1)} (x^{(0)})\, x^{(1)\mu} +
\alpha^{(2)} (x^{(0)})$, and similarly at higher-orders.

We find
\begin{equation}
\label{k01}
 k^{(1)0}(\lambda) = - 2 \1{\Phi} + I_1(\lambda) \, ,
\end{equation}

\begin{equation}
\label{k02}
 k^{(2)0}(\lambda) = - \2{\Phi} + 2 \left( \1{\Phi} \right)^2 - 2 \1{x}^\mu \partial_\mu \1{\Phi}
 + I_2(\lambda) - 2 \1{\Phi} I_1(\lambda) + I_1(\lambda)^2 \, ,
\end{equation}

\begin{equation}
\begin{split}
\label{k03}
 k^{(3)0}(\lambda) =& - \frac{1}{3} \3{\Phi} - \1{x}^\mu \partial_\mu \2{\Phi} - \1{x}^\mu \1{x}^\nu \partial_\mu \partial_\nu \1{\Phi}
 	- 2 \2{x}^\mu \partial_\mu \1{\Phi} + \frac{8}{3} \left( \1{\Phi} \right)^3 \\
 &+ \left( 2 \1{\Phi} - I_1(\lambda) \right)
 	\left[ \2{\Phi} - 2 \left( \1{\Phi} \right)^2 + 2 \1{x}^\mu \partial_\mu \1{\Phi} \right] \\
 & + I_3(\lambda) + 2 I_1(\lambda) I_2(\lambda) + I_1(\lambda)^3 - 2 \1{\Phi} \left(I_2(\lambda) + I_1(\lambda)^2 \right) \, ,
\end{split}
\end{equation}
where
\begin{align}
\label{I1}
 I_1(\lambda) =& \int_{\lambda_\O}^\lambda d\t{\lambda} \1{A}'(\t{\lambda}) \, , \\
\label{I2}
 I_2(\lambda) =& \int_{\lambda_\O}^\lambda d\t{\lambda} \left[ \frac{1}{2} \2{A}'(\t{\lambda})
 + \1{x}^\mu \partial_\mu \1{A}'(\t{\lambda}) - 2 \1{\Phi} \1{A}'(\t{\lambda}) \right] \, , \\
 \begin{split}
\label{I3}
 I_3(\lambda) =& \int_{\lambda_\O}^\lambda d\t{\lambda} \Big[ \frac{1}{6} \3{A}' + \frac{1}{2} \1{x}^\mu \partial_\mu \2{A}'
 		+ \frac{1}{2} \1{x}^\mu \1{x}^\nu \partial_\mu \partial_\nu \1{A}' + \2{x}^\mu \partial_\mu \1{A}' \\
 	&\qquad \quad - 2 \1{\Phi} \left( \frac{1}{2} \2{A}' + \1{x}^\mu \partial_\mu \1{A}' \right)
 		+ \1{A}' \left( 2 \1{\Phi}^2 - \2{\Phi} - 2 \1{x}^\mu \partial_\mu \1{\Phi} \right)\Big] \, ,
 \end{split} \\
\label{An}
 A^{(n)} =& \Phi^{(n)} + \Psi^{(n)} \, .
\end{align}
The above expressions are completely known once we find the geodesic $x^{\mu}(\lambda)$ (it is sufficient to compute them up to second-order,
as it is clear by looking at Eq.(\ref{I3})). For the time component the
integration of  $k^{(0)}(\lambda)$ yields
\begin{equation}
 x^{(1)0}(\lambda) = \int_{\lambda_\O}^{\lambda} d\t{\lambda} \left[- 2 \1{\Phi}
 	+ \left(\lambda - \t{\lambda}\right) \1{A}' \right] \, ,
\end{equation}

\begin{equation}
\begin{split}
 x^{(2)0}(\lambda) =& \int_{\lambda_\O}^{\lambda} d {\tilde \lambda}
 	\left[- \2{\Phi} + 2 \left( \1{\Phi} \right)^2 - 2 \1{x}^\mu \partial_\mu \1{\Phi} - 2 \1{\Phi} I_1(\t{\lambda}) + I_1(\t{\lambda})^2 \right. \\
 &\qquad \quad + \left. \left(\lambda - \t{\lambda}\right) \left( \frac{1}{2} \2{A}'(\t{\lambda}) + \1{x}^\mu \partial_\mu \1{A}'(\t{\lambda})
 	- 2 \1{\Phi} \1{A}'(\t{\lambda}) \right) \right] \, .
\end{split}
\end{equation}
Notice that from Eq.(\ref{geonp}) to Eq.(\ref{I3}) we have never used the large scale approximation, so these expressions are general and
we recover exactly the corresponding expressions for the second-order quantities already found in Refs.~\cite{Pyne:1993np,Mollerach:1997up}.
On the other hand to compute the spatial components $x^i(\lambda)$ we do make use of the large scale approximation in the following way.
We start from the non-perturbative spatial geodesic equation
\begin{equation}
\label{gsnp}
	\frac{d k^i}{d \lambda} = 2 \frac{d \Psi}{d \lambda} k^i - \left( \partial^i \Psi + \partial^i \Phi \right) \delta_{jk} k^j k^k \, .
\end{equation}
To find the solution of this equation we first decompose $k^i = k_{||} e^i + k^i_{\perp}$ where
$k_{||} = e_j k^j$ and $k^i_\perp = (\delta^i_j - e^i e_j) k^j$ are the parallel and transverse parts of the wavevector
with respect to the background geodesic. Then, in the geodesic equation we approximate
$\delta_{ij} k^j k^k = k_{||}^2 + \delta_{jk} k_\perp^j k_\perp^k \simeq k_{||}^2$,
since the second term will be negligible on large angular scales. Then, projecting Eq.(\ref{gsnp}), we can split it in the two equations
\begin{align}
\frac{d k_{||}}{d \lambda} =& 2 \frac{d \Psi}{d \lambda} k_{||} - e_i \left( \partial^i \Psi + \partial^i \Phi \right) k_{||}^2 \, , \\
	\frac{d k_\perp^i}{d \lambda} =& 2 \frac{d \Psi}{d \lambda} k_\perp^i
		- (\delta^i_j - e^i e_j) \left( \partial^j \Psi + \partial^j \Phi \right) k_{||}^2 \, ,
\end{align}
which can be solved. The solution for $k_{||}$ is analogous to the one for $k^0$
\begin{equation}
\label{solkp}
k_{||} = e^{2 \Psi} \left[ - 1 + \int_{\lambda_\O}^\lambda d \tilde{\lambda} e^{2 \Psi}
	e^i \partial_i \left( \Psi + \Phi \right) \right]^{-1} \, .
\end{equation}
from which the solution for $k^i_{\perp}$ is easily obtained
\begin{equation}
\label{solkt}
k_\perp^i = - e^{2 \Psi} \int_{\lambda_\O}^\lambda d \tilde{\lambda} e^{2 \Psi}
	(\delta^{ij} - e^i e^j) \partial_j \left( \Psi + \Phi \right)
	\left[ - 1 + \int_{\lambda_\O}^{\tilde{\lambda}} d \tilde{\lambda} e^i \partial_i \left( \Psi + \Phi \right) \right]^{-2} \, ,
\end{equation}
since as initial condition we have $k_\perp^i (\lambda_\O) = 0$. Perturbing Eqs.~(\ref{solkp}) and~(\ref{solkt}) up to second-order we find
\begin{equation}
\label{ki1}
k^i_{(1)} = - 2 \1{\Psi}e^i - \int_{\lambda_\O}^\lambda d \tilde{\lambda} \partial^i \1{A}\, ,
\end{equation}
\begin{equation}
\begin{split}
\label{ki2}
k^i_{(2)} =& - \2{\Psi} e^i - 2 x_{(1)}^\mu \partial_\mu \1{\Psi} e^i - 2 \1{\Psi}^2 e^i
	- 2 \1{\Psi} \int_{\lambda_\O}^\lambda d\tilde{\lambda} \, \partial^i \1{A} \\
	&+ \int_{\lambda_\O}^\lambda d\tilde{\lambda} \left[ - \frac12 \partial^i \2{A} - x^\mu_{(1)} \partial_\mu \partial^i \1{A}
	- 2 \1{\Psi} \partial^i \1{A} - 2 \partial^i \1{A} \int_{\lambda_\O}^{\tilde{\lambda}} d\lambda' e^j \partial_j \1{A} \right]\, .
\end{split}
\end{equation}
Notice that in fact the above expressions for $k^{i}$ are exact, as we have checked explicitly by using the pertubative method of Eqs.~(\ref{fterms})-(\ref{fterms3}).
The reason is that up to this order the large-scale approximation we used also corresponds to neglect higher-order terms in the perturbations.
Moreover, as a further consistency check, we have explicitly verified that the expressions~\eqref{k01}-\eqref{k03} and~\eqref{ki1}-\eqref{ki2}
are the same one gets with the perturbative method employed in Sec~\ref{2S}.

The integrals~(\ref{I1})-(\ref{I3}) along the line of sight are at the origin of the late (and early) Integrated Sachs-Wolfe effect, since
they express the redshift the photons suffer from the last scattering surface to the observer while traveling through a time varying gravitational
potential. In particular the first integrand terms, $A^{(2)'}$ and $A^{(3)'}$ in each of the integrals at second and third order,
$I_2(\lambda)$ and $I_3(\lambda)$, are just the straightforward extension of the well-known first order ISW effect. However one must consider
all the additional contributions which are integrated terms and which depend on the time variation of the gravitational potentials.
So terms of this type will appear also from the first two lines of Eq.(\ref{k03}). Notice in particular that terms containing spatial gradients
like  $x^{(r)j} \partial_j A^{(m)'}$ (and similar) coming from the combination $x^{(r)\mu} \partial_\mu A^{(m)'}=x^{(r)j} \partial_j A^{(m)'}+x^{(r)0}
\partial_0A^{(m)'}$   cannot be neglected a priori even in the large scale approximation because they are integrated along the line of sight.
\footnote{In fact one can check that for a term like, for example, $x^{(r)j} \partial_jA^{(m)'}$ only the parts of such gradients
transverse to the background geodesic can be neglected on large scales, while the longitudinal part combines with  $x^{(r)0} \partial_0 A^{(m)'}$
to give a term which does not contain spatial gradients.}

Finally plugging the expressions~(\ref{k01})-(\ref{k03}) for $k^{(0)}$ up to third-order into Eq.(\ref{T3scalars})
the expression for the third order temperature fluctuations on large scales is
\begin{equation}
\label{T3scalarsls}
\begin{split}
  &\frac{\3{\delta} T}{T} = \3{\tau} + \frac{1}{6} \3{\Phi_{\mathcal E}}+
\2{\tau} \1{\Phi_{\mathcal E}} + \1{\tau} \frac{1}{2} \2{\Phi_{\mathcal E}}
 	+ \1{\tau} \frac{1}{2} \1{\Phi_{\mathcal E}}^2+ \frac{1}{2} \1{\Phi_{\mathcal E}} \2{\Phi_{\mathcal E}}+ \frac{1}{6} \left( \1{\Phi_{\mathcal E}} \right)^3\\
& + x^{(1)0}  \2{\Phi_{\mathcal E}}' + \left( x^{(1)0} \right)^2 \1{\Phi_{\mathcal E}}''
 	+ 2  x^{(2)0} \1{\Phi_{\mathcal E}}'  + 2 x^{(1)0} \1{\Phi_{\mathcal E}}' \left(\1{\Phi_{\mathcal E}} + \1{\tau} \right) \\
        &- I_{3{\mathcal E}} - I_{2{\mathcal E}} \left[ \1{\tau} + \1{\Phi_{\mathcal E}} \right]
        - I_{1{\mathcal E}} \left[ \2{\tau} + \1{\tau} \1{\Phi_{\mathcal E}} + \frac{1}{2} \2{\Phi_{\mathcal E}} + \frac{1}{2} \1{\Phi_{\mathcal E}}^2 + 2 x^{(1)0} \1{\Phi_{\mathcal E}}' \right] \\
 	&+ d_{(1)}^i \left[ \frac{\partial \2{\tau}}{\partial d^i}
 		+ \left( \1{\Phi_{\mathcal E}} - I_{1{\mathcal E}} \right) \frac{\partial \1{\tau}}{\partial d^i} \right]
+ \frac{1}{2} d_{(1)}^i d_{(1)}^j \frac{\partial^2 \1{\tau}}{\partial d^i \, \partial d^j} + d_{(2)}^i \frac{\partial \1{\tau}}{\partial d^i} \\
 	&+ \1{x}^0 \Bigg[ \partial_0\left(-\2{\Phi_{\mathcal E}} + 2 \1{\Phi_{\mathcal E}}^2 - 2 x^{(1)0} \1{\Phi_{\mathcal E}}' \right)
 		+ \frac{1}{2} \2{A_{\mathcal E}}' + x^{(1)0}  \1{A_{\mathcal E}}'' - 4 \1{\Phi_{\mathcal E}} \1{A_{\mathcal E}}' \\
 		&\quad \qquad + I_{1{\mathcal E}} \1{A_{\mathcal E}}'
 		+ I_{1{\mathcal E}} \left( - 2 \1{\Phi_{\mathcal E}}'+ \1{A_{\mathcal E}}' \right)
 		- \frac{1}{2} \1{x}^0 \left( \1{A_{\mathcal E}}''- 2 \1{\Phi_{\mathcal E}}'' \right) \\
 		& \qquad \quad - \left( 2 I_{1{\mathcal E}} - 3 \1{\Phi_{\mathcal E}} - \1{\tau} \right)
 		\left( \1{A_{\mathcal E}}' - 2 \1{\Phi_{\mathcal E}}' \right) \Bigg]
 	+ 2 \2{\Delta \lambda} \1{\Phi_{\mathcal E}}' - \2{\Delta \lambda} A_{\mathcal E}^{(1)'}\, ,
\end{split}
\end{equation}
where $A^{(r)}$ is given by Eq.~(\ref{An}) and we recall that $\2{\Delta \lambda} = - x^{(2)0} + k^{(1)0} x^{(1)0}$. We recognize various contributions to the
temperature anistropies in Eq.~(\ref{T3scalarsls}).
The first line includes the Sachs-Wolfe effect, which combines the intrinsic temperature fluctuations
$\tau$ with the photon redshift on the last scattering surface due to the gravitational potential perturbations ${\Phi_{\mathcal E}}$.
The terms in the third line correspond to the Integrated  Sachs-Wolfe effect, while in the fourth line there are contributions due to a lensing effect
at the last scattering surface, being dependent on the vector ${\bf \hat{d}}$ which specifies the direction of emission of the photon.
Finally all the remaining terms are due to a possible time dependence of the
gravitational potentials at the last scattering epoch, due to the fact that by that time the universe is still not completely
matter dominated (similarly to the early integrated Sachs-Wolfe effect). At second-order, in Eq.~(\ref{T2scalars}), the lensing like
contributions appear as well. Notice that in order to make a comparison with the bispectrum computed using the different technique of Ref.~\cite{CZRec},
one has to take into account also these lensing terms along with the Sachs-Wolfe effect.

\section{Third-order scalar perturbations of a flat $\Lambda$CDM Universe}
\label{e3p}
\subsection{Evolution of the gravitational potentials on large scales}
We now consider a spatially flat Universe filled with a cosmological constant $\Lambda$
and a non-relativistic perfect pressureless fluid and we solve the third order Einstein equations for the scalar perturbations,
in the large scale limit. The energy momentum tensor of the matter component reads $T^\mu_\nu = \rho u^\mu u_\nu$,
with energy density $\rho$ and four-velocity $u^\mu$.
The derivation of the relevant equations is sketched in the Appendix.

The evolution equation for the third order gravitational potential $\3{\psi}$ reads
\begin{equation}
\label{psi}
    \psi''_{(3)} + 3 \H \psi'_{(3)} + a^2 \Lambda \psi_{(3)} = S(\eta, \mathbf{x}) \, ,
\end{equation}
where the source term is given by
\begin{equation}
\label{source}
 \begin{split}
  S(\eta, \mathbf{x}) =& \frac{4}{\H^2 \Omega_m}
  	\nabla^{-2} \Big\{ \left[ \left( \H^2 + a^2 \Lambda \right) E + \H E' \right]
  	\left[ 3 \nabla^{-2} \partial_i \partial^j \left( \varphi_0 \partial^i \varphi_0 \partial_j \varphi_0 \right)
  		- \varphi_0 (\partial^l \varphi_0) (\partial_l \varphi_0) \right] \\
  	&\qquad + \left[ \left( \H^2 + a^2 \Lambda \right) F + \H F' \right]
  		\left[ 3 \nabla^{-2} \partial_i \partial^j \left( \partial^i \alpha_0 \partial_j \varphi_0 + \partial^i \varphi_0 \partial_j \alpha_0 \right)
  		- 2 (\partial^l \alpha_0) (\partial_l \varphi_0) \right] \\
  	&\qquad - \left( g' + 2 \H g \right) \left[ 3 \H (g')^2 + a^2 \Lambda g^2 (2 f -1) \right] \varphi_0 \alpha_0 \Big\} \\
  	&- 6 (g')^2 \varphi_0^3 + 3 \left( g' A' \varphi_0^3 + g' B' \varphi_0 \alpha_0 \right) + \textrm{gradients} \, ,
 \end{split}
\end{equation}
where $\nabla^{-2}$ stands for the inverse of the Laplacian operator.
Here $g(\eta)$ is the linear growth function, defined by $\varphi(\mathbf{x},\eta) = g(\eta) \varphi_0(\mathbf{x})$,
where $\varphi_0(\mathbf{x})$ is the peculiar gravitational potential linearly extrapolated to the present time,
and $f(\eta) = 1 + \frac{g'(\eta)}{\H g(\eta)}$.
The growth-suppression factor is given by $g(\eta)=D_+(\eta)/a(\eta)$ where $D_+(\eta)$ is the linear growing-mode
of density fluctuations in the Newtonian limit.
The exact form of $g$ can be found in Refs.~\cite{lahav,Carroll,Eisenstein}.
A very good approximation for $g$ as a function of redshift $z$ is given in Refs.~\cite{lahav,Carroll}
\begin{equation}
g \propto \Omega_m\left[\Omega_m^{4/7} - \Omega_\Lambda +
\left(1+ \Omega_m/2\right)\left(1+ \Omega_\Lambda/70\right)\right]^{-1} \;,
\end{equation}
with $\Omega_m=\Omega_{0m}(1+z)^3/E^2(z)$,
$\Omega_\Lambda=\Omega_{0\Lambda}/E^2(z)$,
$E(z) \equiv (1+z) {\mathcal H}(z)/{\mathcal H}_0 = \left[\Omega_{0m}(1+z)^3 +
\Omega_{0\Lambda}\right]^{1/2}$ and
$\Omega_{0m}$, $\Omega_{0\Lambda}=1-\Omega_{0m}$, the present-day
density parameters of non-relativistic matter and cosmological constant,
respectively. We will normalize the growth-suppression factor so that
$g(z=0)=1$. The function $f(\eta)$ can be written as a function of $\Omega_m$ as $f(\Omega_m)
\approx \Omega_m(z)^{4/7}$ ~\cite{lahav,Carroll}. In the $\Lambda=0$ case $g=1$ and $f(\eta)=1$.

In Eq.(\ref{source}) $A(\eta)$, $B(\eta)$ enter in the solutions for the second order gravitational potentials
on large scales \cite{fullT}:
\begin{align}
 \2{\Psi}(\mathbf{x},\eta) =& A(\eta) \varphi_0^2(\mathbf{x}) + B(\eta) \alpha_0(\mathbf{x}) \, \\
 \2{\Phi}(\mathbf{x},\eta) =& A(\eta) \varphi_0^2(\mathbf{x}) + C(\eta) \alpha_0(\mathbf{x}) \,
\end{align}
and their explicit expressions are
\begin{align}
  \label{al0} \alpha_0(\mathbf{x}) =& \nabla^{-2} \left[ \partial^i \varphi_0 \partial_i \varphi_0
      - 3 \nabla^{-2} \partial_i \partial^j (\partial^i \varphi_0 \, \partial_j \varphi_0) \right] \, , \\
  \label{A} A(\eta) =& B_1(\eta) + 2 g^2(\eta) - 2 g_m g(\eta)
      - \frac{10}{3} \left( a_{NL} - 1 \right) g_m g(\eta) \, , \\
  \label{B} B(\eta) =& B_2(\eta) - \frac{4}{3} g_m g(\eta) \, , \\
  \label{C} C(\eta) =&  B(\eta) + \frac{4}{3} g_m g(\eta) \left( \frac{f^2(\eta)}{\Omega_m(\eta)} + \frac{3}{2} \right) \, ,
\end{align}
where the functions $B_1(\eta)$ and $B_2(\eta)$ are given in Ref.~\cite{fullT} while $g(\eta_m)=g_m$ is the value of the growth suppression
factor during matter domination, when the cosmological constant was still negligible. A good approximation is
$g_m \approx \frac{2}{5} \Omega_{0m}^{-1}(\Omega_{0m}^{4/7} + \frac{3}{2}\Omega_{0m})$.

Notice that in $A(\eta)$ appears the parameter $a_{\rm NL}$ specifying the
level of quadratic primordial non-Gaussianity which depends on the particular scenario for the generation of cosmological perturbations, as we will
discuss
later in more detail. For example, for standard single-field models of slow-roll inflation $a_{\rm NL}=1+{\cal O}(\epsilon, \eta)$
~\cite{ABMR,maldacena,BMR2,review}, where $\epsilon$ and $\eta$ are the standard slow-roll parameters~\cite{lrreview}. The reason it
appears in Eq.~(\ref{source}) is due to the fact that at third-order some of the source terms contain products of linear and
second-order gravitational potentials, as detailed in the Appendix.

We have also introduced the functions
\begin{align}
  \begin{split}
  \label{E}
  E(\eta) =& 2 A' g' + \left( 5 \H^2 - a^2 \Lambda \right) A g - 2 \H g' g^2
  + 2 \H (A' g + g' A) - 2 g (g')^2 \\
  &- 2 g (g' + \H g)^2 + 2 \frac{g f}{\Omega_m} \left( g' + \H g \right)^2 \, ,
  \end{split}\\
  \label{F}
  F(\eta) =& \frac{1}{2} B' g' + \frac{1}{4} \left( 5 \H^2 - a^2 \Lambda \right) C g
        + \frac{1}{2} \H (C g' + B' g) \, .
\end{align}
Notice that when $\Lambda=0$ then $B_1(\eta)$ and $B_2(\eta)\rightarrow 0$,$B(\eta)\rightarrow-4/3$ and
$A(\eta)\rightarrow 2 -(10/3) a_{\rm NL}$, and therefore $E(\eta)$ simplifies to $E(\eta)=-(50/3) a_{\rm NL} {\cal H}^2$ while
$F(\eta)=(5/2){\cal H}^2$. Thus, looking at Eq.(\ref{source}), it is simple to see that all the terms in the last two lines are specific only to
the case of a non-vanishing cosmological constant.

The solution of Eq.\eqref{psi} is given by the Green's formula
\begin{equation}
\label{green}
 \3{\psi}(\eta) = \frac{g}{g_m} \3{\psi_m}
 	+ \psi_+ (\eta) \int_{\eta_m}^\eta d \eta' \frac{\psi_- (\eta')}{W(\eta')} S(\eta')
 	+ \psi_- (\eta) \int_{\eta_m}^\eta d \eta' \frac{\psi_+ (\eta')}{W(\eta')} S(\eta') \, ,
\end{equation}
where
\begin{equation}
 \psi_+ (\eta) = g(\eta) \, , \qquad
 \psi_- (\eta) = \frac{\H(\eta)}{a^2(\eta)} \, , \qquad
 W(\eta) = \frac{\H_0^2}{a^3(\eta)} \left( f_0 + \frac{3}{2} \Omega_{0m} \right) \, ,
\end{equation}
are respectively the growing and decaying mode solutions and the Wronskian of the homogeneous equation.
Here the suffix `0' stands for the value of the corresponding quantities at the
present time, while $\psi^{(3)}_{m} \equiv \psi^{(3)}(\eta_m)$ represents the
initial condition taken deep in the matter dominated era on super-horizon
scales, $\eta_m$ being the epoch when full matter domination starts. It is such an initial value
that must be properly determined in order to account for the primordial cubic non-Gaussianity in the cosmological perturbations.

The solution for $\3{\phi}$ is then obtained from the relation between $\3{\psi}$ and $\3{\phi}$ obtained in the Appendix (Eq.~\ref{tracelessnp})
\begin{equation}
\label{relationphipsi}
 \3{\phi} = \3{\psi} + 12 g(\eta) A(\eta) \varphi_0^3 + 6 g(\eta) \left( B(\eta) + C(\eta)\right) \alpha_0 \varphi_0
 	- \frac{4}{\H^2 \Omega_m} \left[ E(\eta) \mu_0 + F(\eta) \nu_0 \right] \, ,
\end{equation}
where we have defined the following kernels:
\begin{align}
 \mu_0 (\mathbf{x}) =& 3 \nabla^{-4} \partial_i \partial^j \left(\varphi_0 \, \partial^i \varphi_0 \, \partial_j \varphi_0 \right)
        - \nabla^{-2} \left( \varphi_0 \, \partial^k \varphi_0 \, \partial_k \varphi_0 \right) \, ,\\
 \nu_0 (\mathbf{x}) =& 3 \nabla^{-4} \partial_i \partial^j \left( \partial^i \alpha_0 \, \partial_j \varphi_0 + \partial^i \varphi_0 \partial_j \alpha_0 \right)
       - 2 \nabla^{-2} (\partial^k \alpha_0 \, \partial_k \varphi_0) \, .
\end{align}

\subsection{Initial conditions from primordial non-Gaussianity}
\label{IC}
We now discuss the key issue of the initial conditions,
which are conveniently fixed at the time when the relevant modes of the perturbations are well
outside the Hubble radius.
In order to follow the super-horizon evolution of the density perturbations produced during inflation,
we use the curvature perturbations on uniform-density hypersurfaces $\zeta$,
which will be expanded as $\zeta = \1{\zeta} + \frac{1}{2} \2{\zeta} + \frac{1}{6} \3{\zeta} + \dots$, where at linear order
$\zeta^{(1)}=-\psi^{(1)}-{\cal H} \delta^{(1)} \rho/\bar{\rho}$ was first introduced in Ref.~\cite{BST}. For our purposes it is useful
to adopt a non-perturbative generalization of $\zeta$ given by (see Refs.~\cite{SalopekBond,KMNR,Bartolo:2005fp})
\begin{equation}
\label{znp}
 \zeta = - \Psi + \frac{1}{3(1+w)} \ln \frac{\rho}{\bar{\rho}} \, ,
\end{equation}
where just for simplicity in this expression we have assumed a constant equation of state parameter $w=\bar{p}/\bar{\rho}$,
$\bar\rho$ and $\bar{p}$ being the background energy density and pressure respectively.
This variable has two crucial features: it is \emph{gauge invariant} and, most important, it is \emph{constant} in time
for adiabatic, super-horizon perturbations.\footnote{
As a consistency test, we have used Eq.~(\ref{znp}) together with the Einstein equations contained in the Appendix,
to verify that $\zeta'$ turns out to be proportional to the L.H.S. of Eq.~(\ref{psi}) minus the very same source term
appearing in Eq.~(\ref{source}). This provides an independent check for the expression of the source term which turns out
to be consistent with $\zeta'=0$ on large scales.
}
This allows us to set the initial conditions at the time when $\zeta$ becomes constant,
and follow them until the perturbation mode reenters the horizon.
The primordial non-Gaussianity can be parametrized in terms of the curvature perturbation as~\cite{Bartolo:2005fp}
\begin{equation}
 \2{\zeta} = 2 (a_{\rm NL}-1) \left( \1{\zeta} \right)^2 \, , \qquad \3{\zeta} = 6 (b_{\rm NL}-1) \left( \1{\zeta} \right)^3 \, ,
\end{equation}
where the two non-linearity parameters depend on the physics of a given scenario of generation of the perturbations. For example for standard
single-field models of inflation $a_{\rm NL}=1$ and $b_{\rm NL}=1$ (plus tiny corrections proportional to the slow-roll parameters), while for other
scenarios they might well be non-negligible.

On the other hand, the physical observable quantity is given by the CMB anisotropies, which are one of the best tool to detect or constrain the
primordial non-Gaussianity generated on large scales~\cite{review}.
The standard procedure is to introduce the quadratic and cubic non-linearity parameters, $f_{\rm NL}$ and $g_{\rm NL}$
(which can be in fact also non-trivial kernels in Fourier and harmonic space) characterizing the non-Gaussianity in the large scale temperature
anisotropies~\cite{ks,OkamotoHu,kwmap,LLMR,kwmap3}. In the limit of large non-Gaussianity, $|a_{\rm NL}| \gg 1$ and $|b_{\rm NL}| \gg 1$
~\cite{review,prl,Bartolo:2005fp},
the resulting size of the non-Gaussianity in the CMB anisotropies can be estimated as $f_{\rm NL}\sim 5 a_{\rm NL}/3$ and $g_{\rm NL} \sim
25 b_{\rm NL}/9$ if accounting only for the contribution to the CMB anisotropies from the Sachs-Wolfe effect. In fact for the cubic
non-linearities entering in the Sachs-Wolfe effect such an estimate is too rough,
and the correct expression relating the observable quantity $g_{\rm NL}$ to the primordial
non-linearity parameters $a_{\rm NL}$ and $b_{\rm NL}$ is given by Eq.~(67) of Ref.~\cite{Bartolo:2005fp}.

Our results allow to take into account also
the contribution from the (late) Integrated Sachs-Wolfe effect. The memory of the initial non-Gaussianity is kept in the matter-dominated value of the
large-scale gravitational potential $\psi^{(3)}_m$, and also, as noticed before, in the source term~(\ref{source}). In order to determine
$\psi^{(3)}_m$ we use the expression~(\ref{znp}) for the curvature perturbation in the matter dominated era
\begin{equation}
 \zeta_m = - \Psi_m - \frac{2}{3} \Phi_m \, ,
\end{equation}
where we have employed the energy constraint~(\ref{00}) to replace for the energy density.

We can expand both sides: at first order we find the usual result
\begin{equation}
 \1{\zeta_m} = - \frac{5}{3} g_m \varphi_0 \, ;
\end{equation}
at second order
\begin{equation}
 \2{\phi_m} = - \frac{3}{5} \2{\zeta_m} + 2 g_m^2 \left( \varphi_0^2 + \alpha_0 \right)
	= 2 g_m^2 \left[  - \frac{5}{3} (a_{\rm NL}-1) \varphi_0^2 + \varphi_0^2 +\alpha_0 \right] \, ;
\end{equation}
and at third order
\begin{equation}
 \3{\phi_m} = - \frac{3}{5} \3{\zeta_m} + 6 \varphi_m \2{\phi_m} - 8 \varphi_m^3
 	+ \frac{12}{5 \H^2 \Omega_m} \left[ E(\eta) \mu_0 + F(\eta) \nu_0 \right] \, .
\end{equation}

With the use of Eq.~\eqref{relationphipsi} evaluated in the matter dominated period, we get as initial conditions
\begin{align}
\label{Phiba}
 \3{\phi_m} =& 6 g_m^3 \left[ \left(\frac{25}{9} (b_{\rm NL}-1) - \frac{10}{3} (a_{\rm NL}-1) + \frac{2}{3} \right) \varphi_0^3
 	+ 2 \varphi_0 \alpha_0 - \frac{20}{3} (a_{\rm NL}-1) \mu_0 + \nu_0 \right] \, , \\
 \3{\psi_m} =& 6 g_m^3 \left[ \left(\frac{25}{9}(b_{\rm NL}-1)  + \frac{10}{3} (a_{\rm NL}-1)  + \frac{2}{3} \right) \varphi_0^3
 	- \frac{4}{3} \varphi_0 \alpha_0 - \frac{160}{9}  (a_{\rm NL}-1) \mu_0 + 8 \nu_0 \right] \, .
\end{align}

Finally, from Eqs.~\eqref{green}, \eqref{relationphipsi}, the large-scale solutions for $\3{\psi}$, $\3{\phi}$ are given by
\begin{align}
\begin{split}
 &\3{\psi}(\eta,\mathbf{x}) = 6 g(\eta) g_m^2 \left[ \left(\frac{25}{9} (b_{\rm NL}-1)  + \frac{10}{3} (a_{\rm NL}-1)
		+ \frac{2}{3} \right) \varphi_0^3(\mathbf{x})
 	- \frac{4}{3} \varphi_0(\mathbf{x}) \alpha_0(\mathbf{x}) - \frac{160}{9} (a_{\rm NL}-1) \mu_0(\mathbf{x}) + 8 \nu_0(\mathbf{x}) \right] \\
	&+ \H_0^{-2} \left( f_0 + \frac{3}{2} \Omega_{0m} \right)^{-1}
	\left[ g(\eta) \int_{\eta_m}^\eta d\t{\eta} a(\t{\eta}) \H(\t{\eta}) S(\t{\eta},\mathbf{x})
	- \frac{\H(\eta)}{a^2(\eta)} \int_{\eta_m}^\eta d\t{\eta} a^3(\t{\eta}) g(\t{\eta}) S(\t{\eta},\mathbf{x})\right] \, ;
\end{split}\\
\begin{split}
 &\3{\phi}(\eta,\mathbf{x}) = \3{\psi}(\eta,\mathbf{x}) + 12 g(\eta) A(\eta) \varphi_0^3(\mathbf{x})
 	+ 6 g(\eta) \left(B(\eta) + C(\eta) \right) \alpha_0(\mathbf{x}) \varphi_0(\mathbf{x}) \\
 	&- \frac{4}{\H^2 \Omega_m} \left[E(\eta) \mu_0(\mathbf{x}) + F(\eta) \nu_0(\mathbf{x}) \right] \, .
\end{split}
\end{align}

\section{Conclusions}
In this paper we have computed the expression for the CMB anisotropies due to the
redshift the photons suffer when they travel from the last scattering surface to the observer up to third-order in the gravitational
perturbations.
We thus have completed the findings of Ref.~\cite{Bartolo:2005fp}, where a fully non-linear expresssion for the Sachs-Wolfe effect
has been obtained,
by including also the Integrated Sachs-Wolfe effect and lensing effects
at the last scattering surface. To achieve  this goal we have proposed an
alternative method to the standard perturbative one in order to solve for the geodesic photon equation which provides a fully
non-linear integral solution.
Moreover we have studied the evolution of the gravitational potentials on large scales, allowing for generic non-Gaussian
initial conditions. Our results, together with those
of Ref.~\cite{Bartolo:2005fp}, are of particularly relevance when facing the trispectrum of the CMB anisotropies, in that,
as pointed out in the
Introduction, they include all the relevant cubic non-linearities for a coherent prediction of such statistic from various
cosmological scenarios.
In particular they include those contributions related to primordial non-Gaussianity, some of which would be inevitably lost
if one would stick to the linear evolution
approximation that is often used in the literature.

\section*{Acknowledgments}
The authors wish to thank the anonymous Referee for the suggestions which helped to improve the paper.

\newpage

redpajama_set_name: RedPajamaArXiv


I remember hearing about Surefire flashlights on all the knife forums and thinking to myself, "how good can they be?" I soon found out after taking the plung on a previously owned satin grey Surefire E2. That light didn't leave my pocket from the first day I got it and now in addition to my knife, I always have a flashlight with me. The E2 has long since moved on, but there are many that have taken its place. You'd be surprised (I sure was) how useful a flashlight can be even in the daytime. You don't need to spend a hundred bucks on a light now as there are really great affordable lights out there, but there sure are a lot of expensive lights as well. Great if you are afraid of the dark or just want to "light up the room."
Custom Lights Custom lights are a fairly new product. As opposed to modified lights that take existing production lights and modify them to make them different/better, custom lights are completely custom made. Usually these are designed from the ground up from body to circuitry, but they also sometimes use existing parts (the bulbs are pretty much always existing parts).
Modified Lights We've all had things where we thought to ourselves, "if only it did this," or "if only I could change this." That's basically the theory behind the modified lights. Modders take an existing flashlight and try to improve on the design. There are quite a bit of talented modders out there and as a result, I have much better lights to go along with my much lighter wallet.
Production Lights There used to be only a choice between a Maglite, or a plastic bodied 2 dollar flashlight, but there are quite a few production flashlight companies who have entered the fray. Now consumers have much better options for affordable production lights and with all of the recent disasters the world has seen, it has become more and more apparent that a good light is essential.


Did you know that men who see female physicians are less likely to fill their prescriptions?
The U.S. wasted $317.4 billion last year treating unnecessary medical complications that could have been avoided if patients had taken their medications as prescribed. That's more money than the country spent treating diabetes, heart disease and cancer, combined. It's true – non-adherence is the costliest health condition we face.
You are currently browsing the archives for the Clinical Laboratory category.


Russia has called on Libya to stop using force against civilians and fully comply with UN resolutions.
Russian Foreign Minister Sergei Lavrov, who met two envoys of Libyan leader Muammar Gaddafi in Moscow, said the answer "cannot be called negative".
The envoys, he added, urged the rebels to show the same restraint and called on Nato to stop its bombing campaign.
Moscow has refused to accept the rebels as the legitimate power, and still has formal ties with Col Gaddafi.
After Tuesday's meeting with the Libyan envoys, Mr Lavrov said Tripoli should withdraw armed groups from cities, in line with UN Security Council resolutions.
"The main thing at the moment is to agree with terms and timeframe of a truce," he told reporters.
Earlier on Tuesday, explosions were heard near Colonel Gaddafi's residence in Tripoli.
Libyan officials later took reporters to the area and showed them the buildings, which were on fire after being hit in what appeared to be Nato air strikes.
A government spokesman said one of buildings was a security services office.
In a separate development, the International Criminal Court (ICC) chief prosecutor pressed on Monday for the arrest of Colonel Gaddafi.
Prosecutor Luis Moreno-Ocampo said Colonel Gaddafi, his son Saif al-Islam, and intelligence chief Abdullah al-Sanussi bore the greatest responsibility for "widespread and systematic attacks" on civilians.
ICC judges must still decide whether or not to issue warrants for their arrest. The Libyan government has already said it will ignore the announcement.
The uprising against Col Gaddafi's rule began in February, sparked by similar uprisings in neighbouring Tunisia and Egypt.
It developed into an armed conflict pitting rebels - who control Benghazi and other eastern cities - against the forces of Col Gaddafi, who retain Tripoli and the west.
In March the United Nations approved a resolution authorising "all necessary measures" to protect civilians.
Russia was one of five nations to abstain in the UN vote, and has argued that the Nato campaign goes beyond the UN mandate.


This guest post is by John Foley, former lead solo pilot of the US Navy's Blue Angels. He is the founder and CEO of John Foley Inc., which provides keynote speaking, coaching, and hands on training services for Fortune 500 companies, professional sports teams, healthcare service providers and beyond. The core of his teaching is The Diamond Performance Framework, which translates the Blue Angels' culture of excellence and peak execution for the World of Business.
The flying that I did with the Blue Angels was not dangerous, but it was inherently unforgiving, and nowhere was there less margin for error than "The Knife Edge Pass" maneuver. Imagine two jets headed straight for one another on a near collision course, 100 feet off the ground, closing in at 1000mph! As a pilot, it's your job to rotate your jet at the last split second, brushing past the opposing plane by a mere wingspan! There is no room for improvisation, and Blue Angels master this routine by leveraging a very simple tool called CenterPoint Alignment. By briefing on a CenterPoint, the two opposing pilots know exactly where to be, and exactly when to be there, and it's this precision alignment that transforms the improbable into the possible.
How many times in the Business world do we enter a meeting room with an investor, partner, or customer without being fully prepared and aligned with your team? In a Blue Angel environment, you would very quickly put yourself into a highly dangerous zone. Having the right set of information and sharing it between pilots is the first step to getting aligned on a CenterPoint. When the two solo pilots attend a Blue Angel briefing, they're given all of the necessary information, given the chance to visualize their actions, and mentally carry out the task at hand. Today, that's were Spoke comes in to play. Spoke.com is a great tool to help anyone collect the business information they need and get their team aligned around an information set, just like CenterPoint Alignment. Users scour the web to gather business insights on industries, trends, partners, customers and competitors, and create action-worthy intelligence to share with their team members. With these detailed industry maps, business users and teams understand where their CenterPoint is, and can fully align themselves along what needs to get done.
When I speak about the Blue Angels, I'm speaking about the best of the best, about the top .01% of the top 1% of all pilots. Keeping your team aligned on various CentertPoints can be the difference between performing in the top 1%, and soaring with the top .01% in your industry.
To learn more about John and his cutting-edge business philosophies, visit his website, or find him on Spoke.com.