---
abstract: 'In the context of cosmological perturbation theory, we derive the second order Boltzmann equation describing the evolution of the distribution function of radiation without a specific gauge choice. The essential steps in deriving the Boltzmann equation are revisited and extended given this more general framework: i) the polarisation of light is incorporated in this formalism by using a tensor-valued distribution function; ii) the importance of a choice of the tetrad field to define the local inertial frame in the description of the distribution function is emphasized; iii) we perform a separation between temperature and spectral distortion, both for the intensity and for polarisation for the first time; iv) the gauge dependence of all perturbed quantities that enter the Boltzmann equation is derived, and this enables us to check the correctness of the perturbed Boltzmann equation by explicitly showing its gauge-invariance for both intensity and polarization. We finally discuss several implications of the gauge dependence for the observed temperature.'
author:
- 'Atsushi Naruko$^{1}$, Cyril Pitrou$^{2,3}$, Kazuya Koyama$^{4}$, Misao Sasaki$^{5}$'
date: 'May 8, 2013'
title: 'Second order Boltzmann equation : gauge dependence and gauge invariance'
---

Introduction
============

The non-Gaussianity in the Cosmic Microwave Background (CMB) has been one of the hottest topics in cosmology because it could open a new window for probing the primordial universe. Recent CMB observations, especially WMAP [@Komatsu:2010fb] and Planck [@Ade:2013ydc], have confirmed to a very high accuracy that the primordial curvature perturbations have a nearly scale invariant initial power spectrum and the associated statistics is nearly Gaussian. These observations are consistent with the predictions of an early inflationary era driven by a single slow-rolling scalar field.

The possibility of non-Gaussianity in the primordial curvature perturbations was discussed for the first time quantitatively by Komatsu and Spergel [@Komatsu:2001rj]. They parameterized the level of non-Gaussianity in the potential $\Phi(x)$ (the curvature potential in the Newton or Poisson gauge) by $$\Phi (x) = \Phi_{{\mathrm{L}}} (x) + f_{{\mathrm{NL}}}^{{\mathrm{local}}} \Bigl[ \Phi^2_{{\mathrm{L}}} (x)
 - \langle \Phi_{{\mathrm{L}}}^2 (x) \rangle \Bigr] \,.$$ Here, $\Phi_{{\mathrm{L}}} (x)$ denotes the Gaussian part of the perturbation, or in perturbation theory its linear part, and $\langle \cdots \rangle$ designates the statistical average. This type of non-Gaussianity leads to a non-vanishing three point correlation function, or equivalently in reciprocal space to a non-vanishing bispectrum. The prediction for the possible values of this parameter $f^{{\mathrm{local}}}_{{{\mathrm{NL}}}}$ from a phase of single-field slow-roll inflation was first performed by Maldacena [@Maldacena:2002vr], and it was shown that it is of order of the slow-roll parameters and thus highly suppressed. Hence, if a $f^{{\mathrm{local}}}_{{{\mathrm{NL}}}}$ of order unity or greater is detected, this simplest model of inflation will be ruled out.

Unfortunately, the interpretation of the measured non-Gaussianity is not so straightforward because we do not observe directly the primordial non-Gaussianity in the curvature perturbation but its effect on the CMB fluctuations. Therefore, to relate the primordial curvature perturbation to CMB, we first have to compute the evolution of perturbations after inflation. These effects can be split unambiguously in two parts: i) a linear transfer that cannot create a non-Gaussian signal if the initial conditions are purely Gaussian, and ii) a non-linear transfer that generates a non-Gaussian signal in the observables even if the initial conditions are purely Gaussian. The resulting non-Gaussian signal from i) is often called [*primordial non-Gaussianity*]{} and all the possible sources of non-linear evolutions which enter the category ii) are called [*secondary non-Gaussianity*]{}. Recently the Planck collaboration provided a constraint $f^{{\mathrm{local}}}_{{{\mathrm{NL}}}} = 2.7 \pm 5.8$ [@Ade:2013ydc]. This result was obtained by subtracting one of the secondary non-Gaussianities that arises from the correlation between the lensing and integrated Sachs-Wolfe effect. This clearly demonstrates the importance of subtracting all the secondary non-Gaussianity consistently in order to obtain an accurate constraint on the primordial non-Gaussianity.

The evolution of the perturbations on super-horizon scales is well understood, even fully non-linearly using either a covariant approach [@Langlois:2005ii; @Langlois:2005qp; @Enqvist:2006fs; @Pitrou:2007xy], or a separate universe approach with the so-called $\delta N$ formalism [@Comer:1994np; @Kolb:2004jg; @Lyth:2004gb], since it leads to a conservation law for the curvature perturbations in the case of adiabatic perturbations. On the other hand, the evolution for modes below the horizon scale is not so simple analytically, especially at the non-linear order, and the use of a kinetic description cannot be avoided on small scales since radiation starts to develop an anisotropic stress. In order to obtain numerical results for the non-linear evolution, we need to derive and solve without approximations the coupled system of non-linear a) Einstein equation for the metric, b) conservation and Euler equations for fluids and c) Boltzmann equation for radiation (photons and neutrinos). Note that the conservation and Euler equations can always be deduced from the lowest moments of the Boltzmann equation, and the full set of equations is often only referred to as [*Einstein-Boltzmann system*]{} of equations.

In order to follow this roadmap, the second order Boltzmann equation was written down in the Poisson gauge in Refs. [@Bartolo:2006cu; @Bartolo:2006fj; @Pitrou:2008hy; @Pitrou:2008ut; @Pitrou:2010sn; @Beneke:2010eg]. The gauge dependence of the distribution function was obtained at linear order [@Durrer:1993db] and then at second order [@Pitrou:2007jy] but leaving aside the problem of polarisation. It was then extended to include polarised light in Ref. [@Pitrou:2008hy]. The system of equations was then solved numerically in Fourier space in Poisson gauge in [@Pitrou:2010sn], and it was reported that the secondary effects around the last-scattering surface could mimic a primordial signal of $f_{{\mathrm{N L}}}^{{\mathrm{local}}} \sim 4$. Recently there have been a huge progress in improving the numerical calculations and clarifying the amplitude of various secondary non-Gaussianities at recombination [@Huang:2012ub; @Su:2012gt; @Pettinari:2013he], and a consensus emerged that when including all the non-linear effects around recombination *and* the integrated early effects after recombination, it could mimic a primordial signal of $f_{{\mathrm{N L}}}^{{\mathrm{local}}} \sim 0.8$, as expected from analytic approximations [@Creminelli:2011sq].

The description of the spectral dependence of the distribution function is also crucial at second order. Indeed, at first order there are no spectral distortions and the perturbation of the photon distribution function can be understood as a single, spectrum-independent temperature fluctuation. However, at second order, there appears a deviation from the Planck distribution, resulting in a continuum of spectral distortions, which in principle must superimpose the thermal Sunyaev-Zeldovich effect [@Ade:2013qta]. To describe this distortion, we use a direction and position dependent Compton $y$ parameter [@Stebbins:2007ve; @Pitrou:2009bc], and also introduce a similar tensor-valued variable to describe the distortion in the polarisation, thus extending the formalism introduced in Ref. [@Pitrou:2009bc]. In this paper, we derive the second order Boltzmann equation without restricting to a specific gauge, and including polarisation. Reflecting on the above, our motivation is two-fold.

First, since the structure of the second order Einstein and Boltzmann equations depends very much on a choice of the gauge, we have to find a gauge in which we can numerically solve this system accurately and quickly. Therefore, it is preferable not to specify the gauge from the beginning but to formulate the equations without specifying it. We can impose different gauge restrictions in their final form to explore the stability and efficiency of the numerical integration. Second, we would like to check the equations derived in Refs. [@Pitrou:2008hy; @Pitrou:2010sn; @Beneke:2010eg]. As a direct check, we recover them in the specific case of the Poisson gauge. Then, as an indirect check, we revisit the transformation properties of the distribution function and the metric perturbations and confirm that the perturbed Boltzmann equation is gauge-invariant up to second order in perturbations, thus increasing our confidence in the rather lengthy derivation.

In Ref. [@Pettinari:2013he], it was found that the inclusion or omission of certain line of sight terms can make a large impact on the estimation of the bias to the primordial non-Gaussianity due to the secondary non-Gaussianity. In Refs [@Huang:2012ub; @Pettinari:2013he] all physical effects were included except for lensing and time-delay. These time-integrated effects require a separate analysis because at later times small-scale multipoles get excited and numerically it is very difficult to evolve the equations. In this paper, we point out that the separation of these effects depends on a gauge. Given that the lensing-ISW cross correlation gives the largest bias to the primordial local type non-Gaussianity, one should bear this gauge dependence in mind when separating these time integrated effects in the calculations.

The choice of the gauge and the associated choice of the tetrad field for the distribution function is also crucial in the interpretation of the quantities as observables. These subtle details do not affect our interpretation of observables in the linear theory since it is only relevant for the monopole and the dipole. However it is no longer the case at second order in perturbations. We must understand the transformation properties of the distribution function under a gauge transformation or a change of the inertial frame and determine what is a gauge and a choice of inertial frame that is related to CMB experiments.

The structure of this paper is as follows. In section \[sec:def\], we give the definitions of the variables that we use for the metric, momentum and distribution function. Especially, to express the perturbation of the metric, we use a geometrical $(3+1)$ decomposition, or the ADM [@Arnowitt:1962hi] parametrisation of the metric. At first order, there is no particular advantage in using this formalism, but various expressions are simplified at second order for the choice of the inertial frame that we make. In section \[sec:Boltz\], we derive the second order Boltzmann equation with polarisation without restricting to a specific gauge. In section \[sec:gauge\], we discuss the gauge dependence of the variables. We carefully investigate the gauge transformation of the metric, momentum and the distribution function. We then check explicitly the gauge invariance of the perturbed Boltzmann equation up to second order as a consistency test. Finally, in Section \[sec:conc\], we summarize our results and we comment briefly on the relevance of our formalism for the observed CMB anisotropy. Useful technical details are gathered in the appendices.

Definitions {#sec:def}
===========

In this section we build all the tools which are used for the description of polarized radiation in cosmology. We first review briefly the parametrization of cosmological perturbations, and explain how a photon momentum can be uniquely described by its energy and direction once a suitable tetrad choice has been made. We then introduce the tensor-valued distribution function which is used to treat statistically a gas of polarized photons, and which is the key object in the Boltzmann equation, and we finally present how it can be decomposed into its main spectral components.

Spacetime coordinates and local inertial frame {#ssec:deftet}
----------------------------------------------

We shall use the ADM formalism to write down the expression of the perturbed metric where the metric can be decomposed as $$\begin{aligned}
\label{MetricADM}
 {{\rm d}}s^2 &= a^2(\eta) \Bigl[ - N^2 {{\rm d}}\eta^2
 + \gamma_{ij} ({{\rm d}}x^i + \beta^i {{\rm d}}\eta)
 ({{\rm d}}x^j + \beta^j {{\rm d}}\eta) \Bigr] \nonumber\\
 &= a^2(\eta) \Bigl[ - (N^2 - \gamma_{i j} \beta^i \beta^j) {{\rm d}}\eta^2
 + 2 \gamma_{ij} \beta^j {{\rm d}}x^i {{\rm d}}\eta
 + \gamma_{ij} {{\rm d}}x^i {{\rm d}}x^j \Bigr] \,,\end{aligned}$$ where $N$ is the lapse function, $\beta^i$ is the shift vector, $\gamma_{i j}$ is the spatial metric, and indices of the Latin type ($i,j,k\cdots$) run from $1$ to $3$. To describe the perturbations around the flat Friedmann-Lemaître-Robertson-Walker (FLRW) space-time, perturbation variables, $\alpha$ and $h_{i j}$, are introduced as $$\label{DefPertNNi}
 N \equiv 1 + \alpha \,, \qquad
 \gamma_{i j} \equiv \delta_{i j} + 2 h_{i j} \,.$$ For simplicity, we use the following definitions $$\beta_i\equiv \delta_{ij} \beta^j \,, \quad
 h^i{}_j \equiv \delta^{ik} h_{kj} \,, \quad
 h^{ij} \equiv \delta^{ik} \delta^{jl}h_{kl} \,,$$ where the spatial indices are raised and lowered with $\delta_{ij}$ and $\delta^{ij}$, rather than with $\gamma_{ij}$ and $\gamma^{ij}$. As is clearly seen below, the ADM form of the metric perturbation will simplify the expressions of the perturbed Boltzmann equation. Any perturbation $X$ will be expanded into its first and second order parts as $$X= X^{(1)} + \frac{1}{2} X^{(2)}\,.$$ The relations between the ADM variables and the usual definitions of cosmological perturbations are provided in Appendix \[app:ADM\].

The Boltzmann equation is better formulated by explicitly using a local inertial frame at every point of the space-time and this can be achieved by using a tetrad field. It is a set of four vector fields which satisfy $$\eta_{(a) (b)} = g_{\mu \nu} e_{(a)}{}^\mu e_{(b)}{}^\nu \,, \qquad
 g_{\mu \nu} = \eta_{(a) (b)} e^{(a)}{}_\mu e^{(b)}{}_\nu \,.$$ These conditions determine the choice of tetrad only up to rotations and boosts. Here the following particular tetrads are chosen up to second order accuracy $$\begin{aligned}
 e^{(0)}{}_\mu &= a (- N, 0, 0, 0) \,,\qquad
 e^{(i)}{}_\mu = a \left( \beta^i + h^i{}_j \beta^j,
 \delta^i{}_j + {h^i}_j - \frac{1}{2} h^{i k} h_{k j} \right) \,,
\label{tetrads}\end{aligned}$$ and the inverse tetrads are given by $$\begin{aligned}
 e_{(0)}{}^\mu = - e^{(0) \mu} &= - \frac{1}{a} \left( \frac{1}{N}, - \frac{\beta^i}{N} \right) \,, \qquad
 e_{(i)}{}^\mu = \frac{1}{a} \left[ 0, \delta_{(i) (k)} \left( \delta^{j k}
 - h^{j k} + \frac{3}{2} h^{j l} h_l{}^k \right) \right] \,.\end{aligned}$$ The time-like tetrad is chosen to be orthogonal to the constant time hypersurfaces since ${\bm e}^{(0)} \propto {{\rm d}}\eta$. As for the spatial tetrads, this choice corresponds to asking that there is no rotation between the background and the perturbed tetrads [@Pitrou:2007jy].

Momentum
--------

To facilitate the separation between the magnitude of the momentum and its direction in a covariant manner, let us consider the projection of the momentum of photon $p^\mu$ onto the set of tetrads, $$\begin{aligned}
 p^{(a)} = e^{(a)}{}_\mu p^\mu \,.
 \end{aligned}$$ We introduce the conformal momentum of photon rather than the physical momentum $p^{(a)}$ $$\label{defpP}
 q^{(a)} \equiv a p^{(a)} \,.$$ Since the momentum of photon satisfies the null condition $p^\mu p_\mu = 0$, or equivalently $q^{(a)} q_{(a)} = 0$, only three components among four are independent, that is $$\begin{aligned}
 q^{(a)} q_{(a)} = 0 \,, \quad \Leftrightarrow \quad
 (q^{(0)})^2 = \delta_{(i) (j)} q^{(i)} q^{(j)}  \,.
 \end{aligned}$$ Thus the three spatial components $q^{(i)}$ can be regarded as such independent variables. Furthermore, $q^{(i)}$ can be decomposed into its magnitude $q$ and direction $n^{(i)}$ as $$\begin{aligned}
 q \equiv \sqrt{\delta_{(i) (j)} q^{(i)} q^{(j)}} = |q^{(0)}| \,, \qquad
 n^{(i)} \equiv \frac{q^{(i)}}{q} \,.
 \end{aligned}$$ Physically the above $q$ can be understood as the conformal (re-scaled) energy, $q = a E_{{\mathrm{phys}}}$, seen by an observer orthogonal to time constant hypersurfaces.

From Eq. (\[defpP\]), the components of momentum $p^{\mu}$ are expressed as functions of $(q, n^{(i)})$ up to the second order as

$$\begin{aligned}
\label{p0pi}
 p^0 &= \frac{q}{a^2} (1 - \alpha + \alpha^2) \,, \\
 p^i &= \frac{q}{a^2} \left( n^{(i)} - \beta^i - {h^i}_j n^{(j)}
 + \alpha \beta^i + \frac{3}{2} h^{i k} h_{k j} n^{(j)} \right) \,.
\label{momentum}\end{aligned}$$

Conversely, ($q, n^{(i)}$) are given by the components of momentum as

$$\begin{aligned}
 q &= a^2 (1 + \alpha) p^0 \,,
\label{def:q} \\
 n^{(i)} &= \left[ (1 - \alpha + \alpha^2) \delta^i{}_j + (1 - \alpha) h^i{}_j
 - \frac{1}{2} h^{ik} h_{k j} \right] \frac{p^j}{p^0}
 + (1 - \alpha) \beta^i + \beta^j {h^i}_j \,.
\label{def:n^i}\end{aligned}$$

One can introduce a projection operator in terms of $e^{(0)}{}_\mu$ and $n_\mu$. The projection operator, often called the screen projector, is defined as $$\begin{aligned}
 S_{\mu \nu}
 &\equiv g_{\mu \nu} + e^{(0)}{}_\mu e^{(0)}{}_\nu - n_\mu n_\nu \,,\end{aligned}$$ where $n^\mu$, the direction vector of photon, is defined by $$\begin{aligned}
 n^\mu \equiv e_{(i)}{}^\mu n^{(i)} \,.
 \end{aligned}$$ Clearly $S_{\mu \nu}$ is a projection of the tangent space onto a two dimensional plane orthogonal to both $e^{(0)}{}_\mu$ and $n_\mu$ since $S^{\mu \nu} e^{(0)}{}_\mu$ and $S^{\mu \nu} n_\mu$ vanish. Its expression in tetrad components reduces necessarily to the identity of the two-dimensional subspace which is left invariant by the projector, that is $$\begin{aligned}
 S_{(i) (j)} = \delta_{(i) (j)} - n_{(i)} n_{(j)} \,,\qquad
 S_{(0) (0)}=S_{(0) (i)}=0 \,.\end{aligned}$$

Distribution function for photons and Stokes parameters {#ssec:f}
-------------------------------------------------------

In order to describe the polarisation of radiation, we introduce a tensor-valued distribution function $f_{\mu \nu}$, which is complex valued and Hermitian. The construction of this distribution function is discussed in Appendix \[app:Dist\]. It is independent of the choice of the electromagnetic gauge and contains only four physical degrees of freedom since it satisfies the conditions $$f_{\mu\nu} e_{(0)}{}^\mu = f_{\mu\nu} e_{(0)}{}^\nu
 = f_{\mu\nu} n^\mu = f_{\mu\nu} n^\nu = 0 \,.$$ Note that the distribution function depends on the observer’s velocity, $u^\mu \equiv e_{(0)}{}^{\mu}$, used in its definition. As long as no confusion arises from such dependence, we omit to specify it. In the case where this is needed, mainly when studying the transformation properties of such a quantity, we shall use the notation $f_{\mu\nu}^{\bm{e}_{(0)}}$ to stress that the tensor-valued distribution function is dependent on the observer’s velocity and thus on the choice of the tetrad field.

The four degrees of freedom can be extracted by decomposing $f_{\mu \nu}$ into a trace part, a symmetric traceless part and an antisymmetric part as $$f_{\mu \nu}  \equiv \frac{1}{2} I
 S_{\mu \nu} + P_{\mu \nu}+ \frac{{{\rm i}}}{2} \epsilon_{\rho \mu \nu \sigma}
 e_{(0)}{}^\rho n^\sigma V \,,
 \label{Pol}$$ where the antisymmetric tensor is defined by $$\begin{gathered}
 \epsilon_{\alpha \beta \gamma \delta} = \epsilon_{[ \alpha \beta \gamma \delta] } \,,
 \quad \epsilon_{0 1 2 3} = \sqrt{- g} \,, \qquad {\rm or} \qquad
 \epsilon_{(0)(1)(2)(3)}=-\epsilon^{(0)(1)(2)(3)} = 1 \,.\end{gathered}$$ [*I*]{} is the intensity and $V$ is the degree of circular polarisation. $P_{\mu \nu}$ encodes the two degrees of linear polarisation (so called $Q$ and $U$ Stokes parameters). All these functions, together with the original tensor-valued distribution function, are functions of the position on space-time $x^\mu=(\eta,x^i)$ and on the point in tangent space. This point in the tangent space can be chosen to be parametrized either by the components $p^\mu$ in the basis canonically associated with the coordinates system, or alternatively by their Cartesian counterparts $p^{(a)}$. In fact we will choose to parametrize the tangent space by the components of the conformal momentum in tetrad space, $q^{(i)} = a
p^{(i)}$, expressed in their spherical coordinates $q$ and $n^{(i)}$, as this leads to the most simple form for the Boltzmann equation as we shall see further.

Spectral distortion
-------------------

On the background space-time, the distribution function, which is characterized only by the intensity $I$, is given by a Planck distribution whose temperature $\bar T$ depends only on $\eta$ due to the symmetries of the FLRW universe. As we will check later, the background temperature scales as $\propto 1/a$. We thus have $$\bar I(\eta,q) = {{I_{\rm BB}}}\left[\frac{q}{a(\eta)\bar T(\eta)}\right] \,, \quad
 \text{with} \quad {{I_{\rm BB}}}(x) \equiv \frac{2}{( e^x - 1 )} \,.$$

At first order in perturbation, the fluctuation of intensity can be described as a fluctuation of temperature $\delta T$ which is independent of $q$. There are two reasons for this. First, as we shall discuss further, gravitational interactions do not induce spectral distortions in the sense that they shift all wavelengths by the same ratio. Second, the collisions at linear order in perturbation do not induce spectral distortions and the redistribution of the photon directions resulting from it can be described by a direction dependent temperature. A similar procedure can be followed for the description of polarisation at first order.

However at second order the situation becomes more complicated since the Compton scattering at this order of perturbation induces spectral distortions which cannot be reabsorbed in a simple direction dependent temperature. As a result, the photon distribution is not described by a Planck distribution function, but fortunately it is sufficient to use two direction dependent quantities. The first remains the temperature and the second describes the type of spectral distortion generated at second order. Actually in general, at the $n$-th order, $n$ directional dependent functions would be needed [@Stebbins:2007ve; @Pitrou:2009bc] to characterize fully the spectrum.

In order to parametrize this distortion, we introduce on top of the temperature $T$, the so-called Compton $y$ parameter. In this section we will omit the dependence of all quantities on the coordinates $x^\mu$ and we will focus on the dependence on the tangent space coordinates $(q,n^{(i)})$. The distribution function can be expanded around a Planck distribution in the so-called Fokker-Planck expansion as  [@Stebbins:2007ve] $$\begin{aligned}
 \label{defy}
 I \Bigl( q,n^{(i)} \Bigr)
 &\simeq {{I_{\rm BB}}}\left( \frac{q}{a T} \right)
 + y \bigl( n^{(i)} \bigr) q^{-3} {{\frac{\partial}{\partial \ln q}}}\left[
 q^3 {{\frac{\partial}{\partial \ln q}}}{{I_{\rm BB}}}\left( \frac{q}{a T} \right) \right] \notag\\
 &= {{I_{\rm BB}}}\left( \frac{q}{a T} \right) + y \bigl( n^{(i)} \bigr)
 {{\cal D}}_q^2{{I_{\rm BB}}}\left( \frac{q}{aT} \right) \,,\end{aligned}$$ where $${{\cal D}}_q^2 \equiv q^{-3} {{\frac{\partial}{\partial \ln q}}}\left( q^3 {{\frac{\partial}{\partial \ln q}}}\right)
 = {{\frac{\partial^2}{\partial \ln q^2}}}+3 {{\frac{\partial}{\partial \ln q}}}\,.$$ Because the number density of photon is given by $n \propto a^{-3}\int I q^2 {{\rm d}}q$, the $y$ term does not contribute to the photon number density and the temperature $T$ is the temperature of the black-body that would have the same number density (see Ref. [@Pitrou:2010sn] for a discussion on other possible definitions for the temperature) and we call it here [*number density temperature*]{}. It can be expanded around the background temperature as $$T \bigl( n^{(i)} \bigr)
 \equiv \bar T(\eta) \left[ 1 + \Theta \bigl( n^{(i)} \bigr) \right] \,.$$ Note that the expansion (\[defy\]) is not the same as Eq. ($11$) nor Eq. ($15$) of Ref. [@Stebbins:2007ve]. Indeed, the temperature of the Planck spectrum around which we expand is neither the physically motivated logarithmic averaged temperature of Ref. [@Stebbins:2007ve] nor a fiducial temperature, but another physically motivated temperature (the number density temperature) that suits better to describe the spectral distortion of the type that appears in CMB.

However, when performing perturbations in cosmology, we need to refer to the background space-time temperature $\bar T$, not to the local number density temperature. Thus it is convenient to expand the distribution function around a Planck distribution at $\bar T$ rather than $T$. Expanding Eq. (\[defy\]) in $\Theta$ up to the second order, we obtain the expansion as $$\begin{aligned}
 I &= {{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}- \bigl( \Theta+\Theta^2 \bigr) {{\frac{\partial}{\partial \ln q}}}{{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}+\left( y + \frac{1}{2} \Theta^2 \right) {{\cal D}}_q^2 {{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}\,,
\label{intensity}\end{aligned}$$ where we used the fact that $y$ is at least a second order quantity. Here, in order to simplify the notation, it is implied that $\Theta$ and $y$ depend on $x^\mu$ and $n^{(i)}$. For a given $I$, the spectral components $\Theta$ and $y$ can be extracted by performing different types of integrals on $q$ (see appendix \[AppExtraction\] for details). This expansion is similar to Eq. ($11$) of Ref. [@Stebbins:2007ve] when only second derivatives of the Planck distribution are kept.

Now we want to obtain a similar decomposition for polarisation. Indeed, when dealing with polarisation we also need to expand its spectral dependence in a way similar to what has been performed for the intensity in Eqs. (\[defy\]) and (\[intensity\]), that is we want to separate the polarisation tensor into a spectral distortion $Y_{\mu\nu}$ and non-distorted component ${{\cal P}}_{\mu \nu}$. However, this separation is slightly different given that there is no polarisation on the background and hence there is no term corresponding to the first term in Eq. (\[intensity\]). In the appendix of Ref. [@Stebbins:2007ve], it has been shown that the expansion should be $$\label{defyPbase}
 P_{\mu \nu} \Bigl( q, n^{(i)} \Bigr)
 \simeq - {{\cal P}}_{\mu\nu} \bigl( n^{(i)} \bigr) {{\frac{\partial}{\partial \ln q}}}{{I_{\rm BB}}}\left( \frac{q}{a T} \right)
 + Y_{\mu \nu} \bigl( n^{(i)} \bigr) {{\cal D}}_q^2 {{I_{\rm BB}}}\left( \frac{q}{a T} \right) \,,$$ which is just a consequence of the fact that there is no background polarisation. We will check that $Y_{\mu\nu}$ vanishes at first order as it is not generated by collisions at this order. Similarly to the expansion of the intensity part, we want to expand the distribution function around a Planck spectrum at the background temperature $\bar T$ rather than the local number density temperature $T$. Thus we expand Eq. (\[defyPbase\]) in $\Theta$ up to first order to get $$\begin{aligned}
\label{defyP}
 P_{\mu \nu} &= - (1 + 3 \Theta) {{\cal P}}_{\mu\nu} {{\frac{\partial}{\partial \ln q}}}{{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}+ (Y_{\mu \nu} + \Theta {{\cal P}}_{\mu\nu}) {{\cal D}}_q^2 {{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}\,.\end{aligned}$$ Again here, in order to simplify the notation, it is implied that $\Theta$, ${{\cal P}}_{\mu \nu}$ and $Y_{\mu\nu}$ depend on $x^\mu$ and $n^{(i)}$.

Boltzmann equation {#sec:Boltz}
==================

Now that we have all the tools at hand, we are ready to formulate the Boltzmann equation for polarized radiation in the cosmological context, and extract its spectral components. This section is entirely dedicated to this task. Given that the complete and detailed derivation can be rather lengthy, all details which are not necessary in a first reading are gathered in Appendix \[app:Boltz\_pol\]. We first present the general expression of the Boltzmann equation for a tensor-valued distribution function. Since the Boltzmann equation is the description of how this distribution function evolves along a photon geodesic, it is necessary to perturb the geodesic equation up to second order. We then show how the Boltzmann equation can be split into its main spectral components, that is into a temperature and a distortion. Finally we write the explicit forms of the free-streaming part and the collision part of the Boltzmann equation.

Boltzmann equation {#boltzmann-equation}
------------------

The evolution of the tensor-valued distribution function is dictated by the Boltzmann equation [@Tsagas:2007yx] $$\label{Evolfmunu}
 S_\mu{}^\rho S_\nu{}^\sigma \frac{{{\cal D}}f_{\rho \sigma}}{{{\cal D}}\lambda}
 = C_{\mu \nu} \,,$$ where ${{\cal D}}/ {{\cal D}}\lambda$ is the covariant derivative along a photon trajectory $x^{\mu}(\lambda)$ and the momentum and $C_{\mu \nu}$ is the associated collision term. The explicit form of ${{\cal D}}f_{\mu \nu}/ {{\cal D}}\lambda$ is $$\frac{{{\cal D}}f_{\mu \nu}}{{{\cal D}}\lambda}
 \equiv \nabla_\rho f_{\mu \nu} \frac{{{\rm d}}x^\rho}{{{\rm d}}\lambda}
 + \frac{{{\partial}}f_{\mu \nu}}{{{\partial}}q^{(i)}} \frac{{{\rm d}}q^{(i)}}{{{\rm d}}\lambda} \,,$$ where $\nabla_\mu$ indicates a covariant derivative associated with $g_{\mu \nu}$. Using spherical coordinates $(q,n^{(i)})$ instead of $q^{(i)}$ for the momentum space, the Liouville operator, that is the l.h.s of Eq. (\[Evolfmunu\]), reads $$\label{DDLoperatortensor}
 S_\mu{}^\rho S_\nu{}^\sigma \frac{{{\cal D}}f_{\rho \sigma}}{{{\cal D}}\lambda}
 = S_\mu{}^\rho S_\nu{}^\sigma \nabla_\tau f_{\rho \sigma}
 \frac{{{\rm d}}x^\tau}{{{\rm d}}\lambda}
 + \frac{{{\partial}}f_{\mu \nu}}{{{\partial}}\ln q} \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 + D_{(i)} f_{\mu \nu} \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} \,,$$ where $D_{(i)}$ is a covariant derivative for momentum. The detail of the construction of such derivative is discussed in Appendix \[sapp:covp\]. Thanks to the operation of the projection $S_\mu{}^\rho S_\nu{}^\sigma$ onto $ {{\cal D}}f_{\rho \sigma} / {{\cal D}}\lambda$, the left hand side of the Boltzmann equation can be decomposed into the $I$, $V$ and $P_{\mu \nu}$ parts similarly to Eq. (\[Pol\]); $$\begin{gathered}
 S_\mu{}^\rho S_\nu{}^\sigma \frac{{{\cal D}}f_{\rho \sigma}}{{{\cal D}}\lambda}
 = \frac{1}{2} L[I] S_{\mu \nu} + L[{\bf P}\,]_{\mu \nu}
 + \frac{{{\rm i}}}{2} L[V] \epsilon_{\rho \mu \nu \sigma} e_{(0)}{}^\rho n^\sigma \,,\end{gathered}$$ where the corresponding Liouville operators are defined as $$\label{defofLs}
 L[I] \equiv \frac{{{\cal D}}I}{{{\cal D}}\lambda} \,, \qquad
 L[{\bf P} \, ]_{\mu \nu} \equiv  S_\mu{}^\rho S_\nu{}^\sigma
 \frac{{{\cal D}}P_{\rho \sigma}}{{{\cal D}}\lambda} \,, \qquad
 L[V]\equiv \frac{{{\cal D}}V}{{{\cal D}}\lambda} \,.$$ Here, the operator ${{\cal D}}/ {{\cal D}}\lambda$ on a scalar distribution function $f$, takes the simpler form $$\frac{{{\cal D}}f}{{{\cal D}}\lambda} \Bigl( x^\mu, q, n^{(i)} \Bigr)
 \equiv {{\partial}}_\mu f \frac{{{\rm d}}x^\mu}{{{\rm d}}\lambda}
 + \frac{{{\partial}}f}{{{\partial}}\ln q} \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 + D_{(i)} f \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda}\,.$$ In a similar manner, one can also decompose the collision term, that is the r.h.s of Eq. (\[Evolfmunu\]) as $$\begin{gathered}
 C_{\mu \nu} \equiv \frac{1}{2} C^I S_{\mu \nu} + C^P_{\mu \nu}
 + \frac{{{\rm i}}}{2} C^V \epsilon_{\rho \mu \nu \sigma} e_{(0)}{}^\rho n^\sigma \,,
 \label{PolC}\end{gathered}$$ so that after extracting the trace, symmetric traceless and antisymmetric parts we get obviously $L[I]=C^I$, $L[{\bf P}]_{\mu\nu}
= C^P_{\mu\nu}$ and $L[V]=C^V$. Beware that this does not mean that the intensity, linear polarization and circular polarization evolve independently, since for instance $C^I$ is the “intensity part” of the collision term but it may involve in general all components $I$, $P_{\mu\nu}$ and $V$ of the tensor-valued distribution function. As a matter of fact, Compton collision does indeed intermix intensity and linear polarization, whereas circular polarization evolves independently.

Geodesic equation and momentum evolution
----------------------------------------

From Eq (\[p0pi\]) and the definition of momentum, ${{\rm d}}x^\mu/{{\rm d}}\lambda=p^\mu$ we obtain $$\begin{aligned}
 \frac{{{\rm d}}\eta}{ {{\rm d}}\lambda} &= \frac{q}{a^2} (1 - \alpha) \,,
\label{eta-evo} \\
 \frac{{{\rm d}}x^i}{ {{\rm d}}\lambda} &= \frac{q}{a^2} \left( n^{(i)} - \beta^i
 - {h^i}_j  n^{(j)} \right) \,,
\label{xi-evo}\end{aligned}$$ where only the first order terms are kept. In terms of $(q, n^{(i)})$, the geodesic equation leads to the evolution equation for the conformal energy $$\begin{aligned}
 \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 &= \frac{q}{a^2} \Bigl[ - \alpha_{,i} n^{(i)} + \beta_{i, j} n^{(i)} n^{(j)}
 - h_{ij}{}' n^{(i)} n^{(j)} + \alpha \Bigl( \alpha_{,i} n^{(i)}
 - \beta_{i, j} n^{(i)} n^{(j)} + h_{ij}{}'  n^{(i)} n^{(j)} \Bigr) \notag\\
 & \qquad \qquad
 + \alpha_{, j} h^j{}_i n^{(i)} + \beta^k h_{i j, k} n^{(i)} n^{(j)}
 + (\beta_{k, i} - \beta_{i ,k} + 2 h_{i k}{}') h^k{}_j n^{(i)} n^{(j)} \Bigr]
 \,.
\label{q-evo}\end{aligned}$$ As for the direction evolution, up to first order in perturbations, we obtain $$\frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda}
 = -\frac{q}{a^2} S^{(i) (j)} \Bigl[ \alpha_{,j} - (\beta_{k ,j}
 - h_{j k}{}') n^{(k)} + (h_{j l,k} - h_{k l,j}) n^{(k)}n^{(l)} \Bigr] \,.
\label{n-evo}$$ It is obvious that $n_{(i)}{{\rm d}}n^{(i)}/{{\rm d}}\lambda=1$ as it ought to be since $n^{(i)}$ is a unit vector. We need these expressions only at first order, except for the evolution of $q$ because the background distribution function is constant in space-time (see below).

Before closing this subsection, we mention that we are free to choose another affine parameter than $\lambda$, to label a point on a geodesic. A convenient choice is to take the conformal time $\eta$ at each point of space-time crossed by the geodesic. The advantage of such choice, is that for photons having the same direction (and which thus follow the same path), but not the same energy, the same conformal time $\eta$ would correspond to the same point of the geodesic. We can trade $\lambda$ for $\eta$ using Eq (\[eta-evo\]), that is with $$\begin{aligned}
 \frac{{{\rm d}}\eta}{ {{\rm d}}\lambda} = \frac{q}{a^2} (1 - \alpha) \quad
 & \Longrightarrow \quad
 \frac{{{\rm d}}\lambda}{ {{\rm d}}\eta} = \frac{a^2}{q} (1 + \alpha ) \,.\end{aligned}$$ The evolution of position, conformal energy, and direction, take then the form

$$\begin{aligned}
 \frac{{{\rm d}}x^i}{{{\rm d}}\eta}
 = \frac{{{\rm d}}x^i}{ {{\rm d}}\lambda} \frac{{{\rm d}}\lambda}{{{\rm d}}\eta}
 &= n^{(i)} + \alpha n^{(i)} - \beta^i - {h^i}_j n^{(j)} \,, \\
 \frac{{{\rm d}}\ln q}{{{\rm d}}\eta}
 = \frac{{{\rm d}}\ln q}{ {{\rm d}}\lambda} \frac{{{\rm d}}\lambda}{{{\rm d}}\eta}
 &= - \alpha_{,i} n^{(i)} + \beta_{i, j} n^{(i)} n^{(j)}
 - h_{ij}{}' n^{(i)} n^{(j)} \notag\\
 & \qquad
 + \alpha_{, j} h^j{}_i n^{(i)} + \beta^k h_{i j, k} n^{(i)} n^{(j)}
 + (\beta_{k, i} - \beta_{i ,k} + 2 h_{i k}{}') h^k{}_j n^{(i)} n^{(j)} \,, \\
 \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\eta}
 = \frac{{{\rm d}}n^{(i)}}{ {{\rm d}}\lambda} \frac{{{\rm d}}\lambda}{{{\rm d}}\eta}
 &= - S^{(i) (j)} \left[ \alpha_{,j} - (\beta_{k ,j} - h_{j k}{}') n^{(k)}
 + (h_{j l,k} - h_{k l,j}) n^{(k)}n^{(l)} \right] \,.\end{aligned}$$

Spectral decomposition of the Boltzmann equation
------------------------------------------------

We are now in position of writing down explicitly the Boltzmann equation, expanding the orders of perturbations, and separating the spectral components. Let us first look at the formal structure of the Boltzmann equation, especially focusing on the spectral decomposition. At the background level, the Boltzmann equation yields $$\label{EqbackgroundI}
 \frac{{{\cal D}}}{{{\cal D}}\lambda} {{I_{\rm BB} \left( \frac{q}{a \bar T} \right)}}= \frac{q}{a^2} \left. \frac{{{\partial}}{{I_{\rm BB}}}(x)}{{{\partial}}\eta}
 \right|_{q}
 = - \frac{{{\rm d}}\ln (a \bar T)}{{{\rm d}}\lambda}
 \left.\frac{{{\rm d}}{{I_{\rm BB}}}(x)}{{{\rm d}}\ln x} \right|_{x=q/(a\bar T)}
 = 0 \,.$$ This implies that $\bar I$ has no time dependence and $\bar T$ scales as $1/a$. One can conclude that the Planck distribution does not change in time if the initial distribution is given by the Planck one. This does not mean that the radiation is not losing energy as the universe expands. Indeed, since the physical energy of a photon is not the conformal energy $q$ but $q/a$, then $\bar \rho \propto \int \bar I (q/a)^3 {{\rm d}}q/a \propto a^{-4}$ as expected. This background result for the scaling of $\bar T$ is useful as it implies that only the partial derivative with respect to $q$ on ${{I_{\rm BB}}}[q/(a \bar T)]$ are relevant, and this motivates our use of the conformal energy.

Now the action of the Liouville operator on the intensity, Eq (\[intensity\]), is expanded up to the second order as $$\begin{aligned}
 L[I] &= \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}- L \Bigl[ \Theta + \Theta^2 \Bigr] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}- \Theta \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \frac{{{\partial}}^2 {{I_{\rm BB}}}}{{{\partial}}\ln q^2}
 + L \left[ y + \frac{1}{2} \Theta^2 \right] {{\cal D}}_q^2 {{I_{\rm BB}}}\notag\\
 &= - \left[ L [\Theta] - \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 - \Theta \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 + 2 \Theta \left( L [\Theta] - \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \right)
 \right] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \left[ L [y] + \Theta \left( L [\Theta]
 - \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \right) \right] {{\cal D}}^2_q {{I_{\rm BB}}}\,,
\label{original}\end{aligned}$$ where we used ${{\partial}}{{I_{\rm BB}}}/{{\partial}}\eta = 0$.

When we compare this formulation of the Liouville operator with the spectral decomposition (\[intensity\]), we are tempted to say that the expression inside the first square brackets contributes to the evolution of the temperature due to free-streaming, and that the expression inside the second square brackets is very closely related to the evolution of the distortion. It order to give a clear meaning to this assertion we decompose Eq. (\[original\]) according to $$\label{spectraldecI}
 L[I] \equiv \frac{q}{a^2} \left[ - \Bigl( {{\cal L}}^\Theta + 2 \Theta {{\cal L}}^\Theta \Bigr) {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \Bigl( {{\cal L}}^Y + \Theta {{\cal L}}^\Theta \Bigr) {{\cal D}}^2_q {{I_{\rm BB}}}\right] \,.$$ This spectral decomposition is motivated by the fact that i) in the case where the conformal energy is not affected by free-streaming, ${{\rm d}}\ln q / {{\rm d}}\lambda = 0$, then $(q/a^2) {{\cal L}}^\Theta$ simply reduces to $L[\Theta]$; and ii) the prefactor $q/a^2$ is introduced because the Liouville term is expected to be proportional to $q/a^2$, as it can be inferred from the explicit form (\[eta-evo\]) of ${{\rm d}}\eta/{{\rm d}}\lambda$. From a comparison of Eq (\[original\]) with this decomposition, we have $$\begin{aligned}
\label{SpectraldecLI}
\frac{q}{a^2}{\cal L}^\Theta&=L[\Theta] -(1+ \Theta )\frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}\,,\\
 \frac{q}{a^2}{\cal L}^Y&=L[y]\,.\end{aligned}$$ We must bear in mind that in these expressions, even though the operator $L[.]$ has been defined in Eq. (\[defofLs\]) for functions of $(\eta,x^i,q,n^{(i)})$, it is applied on the spectral components $\Theta$ and $y$ which do not depend on $q$.

Since the Liouville operator is equated to the collision term in the Boltzmann equation, it is convenient to decompose the collision term in the same manner as the Liouville term. That is, it is decomposed as $$C^I \equiv \frac{q}{a^2} \left[ - \Bigl( {{\cal C}}^\Theta + 2 \Theta {{\cal C}}^\Theta \Bigr) {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \Bigl( {{\cal C}}^Y + \Theta {{\cal C}}^\Theta \Bigr) {{\cal D}}^2_q {{I_{\rm BB}}}\right] \,,$$ such that the spectral components of the Boltzmann equation can be formally very simple and are given by $$\begin{aligned}
 {{\cal L}}^\Theta= {{\cal C}}^\Theta \,,\qquad {{\cal L}}^Y={{\cal C}}^Y.
 \end{aligned}$$ This decomposition means that once the spectral decomposition of the collision term is known (${\cal C}^\Theta$ and ${\cal C}^Y$), then we only need to obtain the spectral decomposition of the Liouville term from Eqs (\[SpectraldecLI\]).

We follow the same logic for polarization. First, the corresponding Liouville operator reads $$\begin{aligned}
 L [ {\bf P} ]_{\mu \nu}
 &= - L \Bigl[ (1 + 3 \Theta) {{\cal P}}_{\mu\nu} \Bigr] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}- {{\cal P}}_{\mu \nu} \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \frac{{{\partial}}^2 {{I_{\rm BB}}}}{{{\partial}}\ln q^2}
 + L \Bigl[ Y_{\mu \nu} + \Theta {{\cal P}}_{\mu \nu} \Bigr] {{\cal D}}_q^2 {{I_{\rm BB}}}\notag\\
 &= - \left[ (1 + 3 \Theta) L [ {\bm {{\cal P}}} ]_{\mu \nu}
 + 3 \left(L[\Theta]- \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}\right) {{\cal P}}_{\mu \nu} \right] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \left[ L [ {\bf Y} ]_{\mu \nu} + \left(L[\Theta]- \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}\right){{\cal P}}_{\mu \nu}
 + \Theta L [{\bm {{\cal P}}}]_{\mu \nu} \right] {{\cal D}}_q^2 {{I_{\rm BB}}}\,. \end{aligned}$$ For the same reasons as in the case of intensity, it appears natural to decompose this Liouville operator into spectral components according to $$\label{spdec:calP_lhs}
 L [ {\bf P} ]_{\mu \nu}\equiv\frac{q}{a^2} \left\{ - \Bigl[ (1 + 3 \Theta) {{\cal L}}_{\mu \nu}^P
 + 3 {{\cal L}}^\Theta {{\cal P}}_{\mu \nu} \Bigr] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \Bigl( {{\cal L}}^Y_{\mu \nu} + {{\cal L}}^\Theta {{\cal P}}_{\mu \nu}
 + \Theta {{\cal L}}_{\mu \nu}^P \Bigr) {{\cal D}}_q^2 {{I_{\rm BB}}}\right\} \,,$$ which implies that the spectral components are given by $$\label{SpectraldecLP}
\frac{q}{a^2} {{\cal L}}_{\mu \nu}^P \equiv L [ {\bm {{\cal P}}} ]_{\mu \nu},\qquad \frac{q}{a^2} {{\cal L}}_{\mu \nu}^Y \equiv L [ {\bf Y} ]_{\mu \nu} \,.$$ The collision term must then follow the same type of decomposition, that is $$\begin{aligned}
 C^P_{\mu \nu}
 &= \frac{q}{a^2} \left\{ - \Bigl[ (1 + 3 \Theta) {{\cal C}}_{\mu \nu}^P
 + 3 {{\cal C}}^\Theta {{\cal P}}_{\mu \nu} \Bigr] {{\frac{{{\partial}}I_{\rm BB}}{{{\partial}}\ln q}}}+ \Bigl( {{\cal C}}^Y_{\mu \nu} + {{\cal C}}^\Theta {{\cal P}}_{\mu \nu}
 + \Theta {{\cal C}}_{\mu \nu}^P \Bigr) {{\cal D}}_q^2 {{I_{\rm BB}}}\right\} \,,
\label{spdec:calP_rhs}\end{aligned}$$ so that, again, the spectral components of the polarized part of the Boltzmann equation take the formally simple form $$\label{eq:Boltz_pol_dis}
{{\cal L}}^P_{\mu \nu}={{\cal C}}^P_{\mu \nu} \,,\qquad {{\cal L}}^Y_{\mu \nu} ={{\cal C}}^Y_{\mu \nu} \,.$$ Again, this decomposition means that once the spectral decomposition of the collision term is known (${\cal C}^P_{\mu\nu}$ and ${\cal C}^Y_{\mu\nu}$), then we only need to obtain the spectral decomposition of the Liouville term from Eqs (\[SpectraldecLP\]) bearing in mind that the Liouville operator $L[.]$ applies on functions which do not depend on $q$, but only on $(\eta,x^i,n^{(i)})$.

Temperature and spectral distortion of Liouville operators
----------------------------------------------------------

Now that the spectral separation of the Boltzmann equation is performed, it is time to expand the equations obtained in orders of perturbations. In the next two sections, we present such expansion for the temperature and spectral distortion parts of the Boltzmann equation. The case of polarization is reported in Appendix \[sapp:Boltzeq\_pol\].

At first order in perturbation, with Eqs. (\[eta-evo\]) and (\[q-evo\]), the Boltzmann equation leads to $$\begin{aligned}
 {\cal L}^\Theta&= \Theta' + \Theta_{, i} n^{(i)}
 + \alpha_{, i} n^{(i)} - \beta_{i, j} n^{(i)} n^{(j)}
 + h_{i j}{}' n^{(i)} n^{(j)} \,.\end{aligned}$$ Note that there is absolutely no q-dependence, nor scale factor $a$ in this expression, meaning that our spectral decomposition performed in Eq. (\[spectraldecI\]) is adequate. Concerning the spectral distortion part, the Liouville part at first order is ${\cal L}^Y= y'+y_{,i}
n^{(i)}$, but since ${\cal C}^Y=0$ at first order (see section \[ssec:Collision\]), one can conclude that only the temperature part evolves at first order and no spectral distortion is induced.

Up to the second order, the Boltzmann equation for the temperature is given by $$\begin{aligned}
\label{LCThetanotOpened}
 {\cal L}^\Theta
 &= \frac{a^2}{q}\left[L[\Theta]
 - (1 + \Theta) \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}\right] \notag\\
 &= \Theta' + \Theta_{, i} n^{(i)}
 + \alpha_{, i} n^{(i)} - \beta_{i, j} n^{(i)} n^{(j)}
 + h_{i j}{}' n^{(i)} n^{(j)} \notag\\
 & \qquad 
 +\frac{a^2}{q}
 \left[ \left. \frac{{{\rm d}}\eta}{{{\rm d}}\lambda} \right|^{(1)} \Theta'
 + \left. \frac{{{\rm d}}x^i}{{{\rm d}}\lambda} \right|^{(1)} \Theta_{, i}
 + \left. \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} \right|^{(1)} D_i \Theta
 - \left. \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \right|^{(1) \times (1)}
 - \Theta \left. \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda} \right|^{(1)} \right] 
  = {{\cal C}}^\Theta \,.\end{aligned}$$ As for the spectral distortion, the Boltzmann equation is given, even at second order by $$\begin{aligned}
 {\cal L}^Y &= y' + y_{, i} n^{(i)} = {{\cal C}}^Y \,.\end{aligned}$$ This means simply that gravitational effects do not induce spectral distortions, and this result holds actually non-perturbatively.

Before ending this subsection, for the sake of completeness, we shall write down the most general form of the second order Boltzmann equation for the intensity. Using Eqs. (\[eta-evo\]), (\[q-evo\]), (\[xi-evo\]) and (\[n-evo\]), the detailed form of the evolution equation for temperature obtained in Eq. (\[LCThetanotOpened\]) is $$\begin{aligned}
 {\cal L}^\Theta=& \Theta' + \Theta_{, i} n^{(i)} + \alpha_{,i} n^{(i)}
 - \beta_{i ,j} n^{(i)} n^{(j)} + h_{i j}{}' n^{(i)} n^{(j)} \notag\\
 & \quad
 - \alpha \Theta' - (\beta^i + h^i{}_j n^{(j)}) \Theta_{, i}
 - \Bigl[ \alpha_{,i} - (\beta_{j ,i} - h_{i j}{}') n^{(j)}
 + (h_{k i, j} - h_{j k, i}) n^{(j)} n^{(k)} \Bigr] D^i \Theta \notag\\
 & \quad
 - \Bigl( \alpha \alpha_{,i} + \alpha_{, j} h^j{}_i \Bigr) n^{(i)}
 - \Bigl[ \alpha (- \beta_{i ,j} + h_{i j}{}') + \beta^k h_{i j, k}
 + (\beta_{k, i} - \beta_{i, k} + 2 h_{i k}{}') h^k{}_j \Bigr] n^{(i)} n^{(j)}
 \notag\\
 & \quad
 - \Theta n^{(i)} \Bigl( - \alpha_{,i} + \beta_{i ,j} n^{(j)} - h_{i j}{}' n^{(j)} \Bigr) 
 = {{\cal C}}^\Theta \,.\end{aligned}$$

Temperature and spectral distortion of collision terms {#ssec:Collision}
------------------------------------------------------

The expression of the collision term has been derived by taking into account only the intensity in [@Dodelson1993; @Bartolo:2006cu] and then it was extended to include the effect of polarisation in [@Pitrou:2008hy; @Pitrou:2008ut; @Beneke:2010eg]. Here we summarize the result obtained by Beneke et al. [@Beneke:2010eg] applying the decomposition of the distribution function into intensity and linear polarisation.

The complete expression of the collision term for intensity is given by $$\begin{aligned}
 {{\cal C}}^\Theta
 &= a \, \bar{n}_e \sigma_T \left( - \Theta + \langle \Theta \rangle
 - \frac{3}{4} S^{(i) (j)} \Bigl[ \langle \Theta m_{(i) (j)} \rangle
 - 2 \langle {{\cal P}}_{(i) (j)} \rangle \Bigr]
 + v^{(i)} n_{(i)} + {{\cal S}}^T + S^{(i) (j)} {{\cal Q}}^T_{(i) (j)}
  + \delta_e {{\cal C}}^\Theta \right) \,, \\
 {{\cal C}}^Y &= a \, \bar{n}_e \sigma_T \left( - y + \langle y \rangle
 - \frac{3}{4} S^{(i) (j)} \Bigl[ \langle y m_{(i) (j)} \rangle
 - 2 \langle Y_{(i) (j)} \rangle \Bigr]
 + {{\cal S}}^Y + S^{(i) (j)} {{\cal Q}}^Y_{(i) (j)} \right) \,,\end{aligned}$$ where ${{\cal S}}^T$, ${{\cal S}}^Y$, ${{\cal Q}}^T_{(i) (j)}$ and ${{\cal Q}}^Y_{(i) (j)}$ are quadratic contributions defined by $$\begin{aligned}
 {{\cal S}}^T
 &= \Theta^2 + \langle \Theta^2 \rangle 
  - 2 \Theta \langle \Theta \rangle - \Theta v^{(i)} n_{(i)}
  + 2 \langle \Theta \rangle v^{(i)} n_{(i)} 
  - 2  \langle \Theta n_{(i)} \rangle v^{(i)} - \frac{1}{5} v^{(i)} v_{(i)}
 + (v^{(i)} n_{(i)})^2 \,, \\
 {{\cal S}}^Y
 &= \frac{1}{2} \Theta^2 +  \frac{1}{2} \langle \Theta^2 \rangle
 - \Theta \langle \Theta \rangle - \Theta v^{(i)} n_{(i)} 
 + \langle \Theta \rangle v^{(i)} n_{(i)}
 - \langle \Theta n_{(i)} \rangle v^{(i)} 
 + \frac{1}{5} v^{(i)} v_{(i)} + \frac{1}{2} (v^{(i)} n_{(i)})^2 \,,
 \end{aligned}$$ $$\begin{aligned}
 {{\cal Q}}^T_{(i) (j)}
 &= - \frac{3}{4} \langle \Theta^2 m_{(i) (j)} \rangle
 + \frac{3}{2} \Theta \langle \Theta m_{(i) (j)} \rangle
 + \frac{3}{4} \Bigl[ \langle \Theta n_{(j)} \rangle v_{(i)}
 + \langle \Theta n_{(i)} \rangle v_{(j)} \Bigr] \notag\\
 & \qquad
 - \frac{3}{4} n^{(k)} \Bigl[ v_{(j)} \langle \Theta m_{(i) (k)} \rangle
 + v_{(i)} \langle \Theta m_{(j) (k)} \rangle \Bigr]
 - \frac{3}{4} v^{(k)} \Bigl[ 2 n_{(k)} \langle \Theta m_{(i) (j)} \rangle
 + \langle \Theta m_{(i) (j)} n_{(k)} \rangle \Bigr]
 - \frac{1}{5} v_{(i)} v_{(j)} \notag\\
 & \qquad
 + \frac{9}{2} \langle \Theta {{\cal P}}_{(i) (j)} \rangle
 - 3 \Theta \langle {{\cal P}}_{(i) (j)} \rangle
 + \frac{3}{2} n^{(k)} \Bigl[ v_{(j)} \langle {{\cal P}}_{(i) (k)} \rangle
 + v_{(i)} \langle {{\cal P}}_{(j) (k)} \rangle \Bigr] \notag\\
 & \qquad
 + \frac{3}{2} v^{(k)} \Bigl[ \langle {{\cal P}}_{(i) (k)} n_{(j)} \rangle
 + \langle {{\cal P}}_{(j) (k)} n_{(i)} \rangle \Bigr]
 - \frac{3}{2} v^{(k)} \Bigl[ \langle {{\cal P}}_{(i) (j)} n_{(k)} \rangle
 - 2 n_{(k)} \langle {{\cal P}}_{(i) (j)} \rangle \Bigr] \,, \\
 {{\cal Q}}^Y_{(i) (j)}
 &= - \frac{3}{8} \langle \Theta^2 m_{(i) (j)} \rangle
 + \frac{3}{4} \Theta \langle \Theta m_{(i) (j)} \rangle
 - \frac{3}{4} v^{(k)} \Bigl[ n_{(k)} \langle \Theta m_{(i) (j)} \rangle
 - \langle \Theta m_{(i) (j)} n_{(k)} \rangle \Bigr]
 - \frac{1}{20} v_{(i)} v_{(j)} \notag\\
 & \qquad
 + \frac{3}{2} \langle \Theta {{\cal P}}_{(i) (j)} \rangle
 - \frac{3}{2} \Theta \langle {{\cal P}}_{(i) (j)} \rangle
 + \frac{3}{2} v^{(k)} \Bigl[ n_{(k)} \langle {{\cal P}}_{(i) (j)} \rangle
 - \langle {{\cal P}}_{(i) (j)} n_{(k)} \rangle \Bigr] \,,\end{aligned}$$ and $m_{(i) (j)} \equiv n_{(i)} n_{(j)} - \delta_{(i) (j)}/3$. Note that we have introduced the notation $\langle Q \rangle = \int_{\Omega} Q = \int {{\rm d}}^2 n^{(i)} Q$, which corresponds to a multipole extraction that we do not perform explicitly here. Note also that $\delta_e=\delta n_e/\bar n_e$ is the fractional perturbation of the baryons number density, and $v^{(i)}$ are the tetrad components of the baryons spatial velocity. Again, we also defer the expression of the collision term for polarisation to Appendix \[sapp:Boltzeq\_pol\].

Gauge dependence of the distribution function {#sec:gauge}
=============================================

Now that we have established the Boltzmann equation, up to second order, with its spectral components separated, we investigate the gauge dependence of its constituents. Eventually the Boltzmann equation itself should be gauge-invariant, so if we are able to check explicitly that the Boltzmann equation is gauge-invariant, this means that it is very likely that i) the perturbative expansion of the equation is correct; and ii) the gauge transformation rules for all its constituents (metric and distribution function perturbations) are correctly understood. We thus consider this verification as a *consistency test*. This section is dedicated entirely to this task. We first review the gauge dependence for tensors, and deduce how it can be extended to a scalar distribution function. The case of a tensor-valued distribution function, even though it is the less trivial part, is treated in appendix \[AppGTtensor\]. We then infer what should be the transformation rule of the Liouville and collision operators, and in order to complete the consistency test, we check that the perturbed expressions of the Liouville and collision operators do indeed transform following these rules.

Coordinates on the manifold
---------------------------

We need to specify how the functional dependence of the quantities appearing in the Boltzmann equation (and in the Einstein equation) is obtained. If we consider a scalar function $f:{\cal M}\mapsto {\mathbb{R}}$ on the space-time manifold ${\cal M}$, then the choice of a coordinates system [^1] $c:{\mathbb{R}}^4 \mapsto {\cal M}$ does not affect the geometrical meaning of this function, but it affects its functional form $f \circ c:{\mathbb{R}}^4 \mapsto {\mathbb{R}}$, where $\circ$ designates the composition rule, in the sense that for another coordinates system $\tilde c:{\mathbb{R}}^4 \mapsto {\cal M}$, then $f\circ \tilde c \neq f \circ c$. The confusion only arises from the fact that we often refer to $f \circ c$ as $f$ only. A distribution function $f$ (that we take as a scalar-valued for simplicity here) is a function on the tangent bundle $T {\cal M}$ of the manifold and can be regarded as a function on manifold ${\cal M}$ describing the space-time and on the tangent space ${\mathbb{R}}^4$ (more precisely a restriction to the mass shell ${\mathbb{R}}^3$) of each point. It is thus a function $$f: (T{\cal M})\mapsto {\mathbb{R}}\,.$$ Again the particular choice of coordinates on the tangent space does not affect the geometrical meaning of the function but its functional form. Once a choice $c$ of coordinates on the manifold ${\cal M}$ has been made, there is a natural basis, called canonical basis, that is made of the partial derivatives with respect to the coordinates. This leads to a natural coordinates system $Tc$ for the tangent space at each point. Thus $(c,Tc) : {\mathbb{R}}^4 \times {\mathbb{R}}^3 \mapsto T{\cal M}$, is a coordinates system for the tangent bundle. In order to simplify the notation we will note $(c,Tc)$ as simply as $c$.

Furthermore, in this paper we use coordinates in the tangent space described by the tetrad components. More specifically we use the components of the conformal momentum in spherical coordinates in this tetrad basis, $(q,n^{(i)})$. There is a problem with such a choice since the tetrad basis is not unique. However, once a choice $c$ of coordinates on the manifold ${\cal M}$ has been made, the tetrad basis might be completely fixed from the metric through a prescription described in § \[ssec:deftet\]. Once a coordinates system $c$ has been chosen, and the tetrad is fixed thanks to this choice, we obtain the functional form of $f$ in the form $f \circ c:{\mathbb{R}}^4\times {\mathbb{R}}^3\mapsto {\mathbb{R}}$.

Geometrical interpretation of the gauge
---------------------------------------

When performing perturbations around a background FLRW space-time, we need to have a one-to-one correspondence between the background space-time $\overline{\cal M}$ and the physical (and perturbed) space-time ${\cal M}$. This can be completely defined geometrically [@Bruni:1996im] but we take a shorter approach. If we have two sets of coordinates [^2] $\bar c:{\mathbb{R}}^4 \mapsto \overline{\cal M}$ and $c:{\mathbb{R}}^4 \mapsto {\cal M}$ on the background and the perturbed space-time, then we identify points with the same coordinates, that is, we identify points with $c \circ \bar c^{-1}:\overline{\cal M}\mapsto {\cal M}$. Since we could have chosen different sets of coordinates, there is some freedom in this choice, which is known as the gauge freedom.

On the background, the symmetries can justify that we can find a preferred choice of coordinates. For instance for a flat FLRW space-time within a given background cosmology, it is enough to choose that the time coordinate is the proper time of observers with 4-velocity orthogonal to the homogeneous surfaces, that is the proper time of comoving observers. On a given homogeneous surface, there are preferred choices of Cartesian coordinates, since it is conformally related to ${\mathbb{R}}^3$, and all these Cartesian systems on the spatial homogeneous surfaces are related by global translation and rotation in ${\mathbb{R}}^3$ which are irrelevant given the homogeneity. And once a coordinate system has been chosen on a homogeneous surface it can be Lie dragged by the comoving observers to any homogeneous surface. So essentially there is a unique mapping $\bar c$ from ${\mathbb{R}}^4$ to the background manifold $\overline{\cal M}$. This point is illustrated in the left part of Fig. \[fig1\].

However on the physical space-time we could consider another coordinates system $\tilde c:{\mathbb{R}}^4 \mapsto {\cal M}$, and this leads to a different identification through $\tilde c \circ \bar c^{-1}$. The fact that we fix the system of coordinates on the background space-time but there is still some freedom on the physical space-time for the choice of coordinates leads to a freedom in the identification between points of these two space-times. With $c$, a given point $P \in {\cal M}$ would be labelled by the coordinates $x^\mu$, that is $c( {\bm x}) = P$, and with $\tilde c$ it would be labeled by the coordinates $\tilde x^\mu$, that is $\tilde c( \tilde{\bm x}) = P$ (see the upper part of Fig. \[fig1\]). For every point, there exist four numbers ${\xi}^\mu= ({ T},{ L}^i)$ such that $$\begin{gathered}
\label{Trulexmu}
 \tilde{x}^\mu( {\bm x}) = x^\mu + {\xi}^\mu( {\bm x}) \,.\end{gathered}$$ We should note that, although there is an index in the notation, ${\xi}^\mu$ is not a vector field on ${\cal M}$ nor on $\overline{\cal M}$, but it can be seen as a vector field on ${\mathbb{R}}^4$. In the literature, a vector field $\zeta^\nu$ is often used to generate a coordinate transformation  [@Bruni:1996im] $$\tilde x^\mu = \exp\left({\cal L}_{{{\boldsymbol{\zeta}}}}\right)x^\mu
 =x^\mu+\zeta^\mu+\frac{1}{2}\zeta^\mu{}_{,\nu} \zeta^\nu+\cdots\,.$$ In this paper we adopt the former definition Eq. (\[Trulexmu\]) and care must be taken when comparing our transformation rules with those in the literature. Note that in the rest of this section, we will use extensively the notation $$\tilde{\bm x} \equiv \tilde x^\mu\,,\qquad {\bm x}\equiv x^\mu\,,$$ even though $x^\mu$ and $\tilde x^\mu$ are not vectors but just coordinates.

![\[fig1\] In this figure, we represent all four dimensional spacetimes with only two dimensions. First, we have noted that there is essentially a unique way to map ${\mathbb{R}}^4$ to the background manifold, that is to relate the top-left to the bottom left. However, there are several ways to relate ${\mathbb{R}}^4$ to the perturbed manifold, and consequently to relate the background manifold to the perturbed manifold. We represented a coordinates system $c$ and another coordinates system $\tilde c$ which relate ${\mathbb{R}}^4$ (top-left) to the physical manifold (top-right). For each point $P$ of the physical manifold, there is a set of four numbers $x^\mu$ and another set of four numbers $\tilde
          x^\mu$ such that $c(x^\mu) = \tilde c(\tilde
          x^\mu)$. Furthermore, each coordinates system has different surfaces of constant time, and thus different set of tetrads, given that the null tetrad is always chosen to be orthogonal to the constant-time surfaces. When using the tetrad field to extract the components of a given momentum (that is a point in the tangent bundle), this will lead to different components, depending on the coordinates system chosen, but since this is the same momentum at the same point and the tetrad are normalized, we can relate the components by a Lorentz transformation. For a given point of the tangent bundle, that is, for a given point of the manifold and a given momentum, we always have $c(x^\mu,p^{(a)})=\tilde c(\tilde x^\mu, p^{\tilde{(a)}})$.](PlotGauge.pdf){width="\textwidth"}

Metric transformation
---------------------

In general the coordinate transformation rule of any tensorial quantity is given by $$\label{GTTensor}
 {\bm T} \circ \tilde c{(\tilde {\bm x})} = {\bm T}\circ c({\bm x})\,.$$ This means that it is invariant under a coordinates transformation since it is geometrically defined and independent of the coordinates used to parametrise the manifold. From this we can deduce the gauge transformation of its components, which is defined as the transformation when the tensors are compared at the same coordinate in two different coordinates systems. Let us consider in particular the metric. Since $g_{\mu \nu} \equiv {\bm g}({{\partial}}/{{\partial}}x^\mu,{{\partial}}/{{\partial}}x^\nu)$ and $g_{\tilde \mu \tilde \nu} \equiv
 {\bm g}({{\partial}}/{{\partial}}\tilde x^\mu,{{\partial}}/{{\partial}}\tilde x^\nu)$, we then deduce the usual coordinate transformation rule of a 2-form $$\label{Transfog1}
g_{\tilde \mu \tilde \nu} \circ \tilde c{(\tilde {\bm x})} = \frac{{{\partial}}x^\alpha}{{{\partial}}\tilde x^\mu} \frac{{{\partial}}x^\beta}{{{\partial}}\tilde x^\nu}  g_{\alpha \beta}\circ c({\bm x})\,.$$ Expanding the left hand side which is evaluated at $\tilde {\bm x}$ around ${\bm x}$ leads the gauge transformation rule up to the second order $$\begin{gathered}
\label{EqTmetric}
g_{\tilde \alpha \tilde \beta}\circ \tilde c{\left({\bm x}\right)} =
g_{\alpha \beta}\circ c({\bm x})
 - {{\cal L}}_{{{{\boldsymbol{\xi}}}}} g_{\alpha \beta}\circ c({\bm x}) + \frac{1}{2}
 {{\cal L}}^2_{{{{\boldsymbol{\xi}}}}} g_{\alpha \beta}\circ c({\bm x})
 + \frac{1}{2} {{\cal L}}_{{{{\boldsymbol{\xi}}}} {{{\boldsymbol{\xi}}}}} g_{\alpha \beta}\circ c({\bm x}) \,,\end{gathered}$$ where ${{{\boldsymbol{\xi}}}} {{{\boldsymbol{\xi}}}}$ designates the quantity ${\xi}^\mu{}_{,\nu}{\xi}^\nu$. We emphasize that the coordinate transformation  is the transformation of components at the same point of space-time, whereas the gauge transformation  is the transformation of the components at different points which have the same coordinates in two different coordinates systems. Since the notation can become rather cumbersome if we specify which coordinates system $c$ or $\tilde c$ is to be used, we shall indicate it only when it is the new coordinates system $\tilde c$.

Throughout this paper, we use the symmetrization and anti-symmetrization definitions $$X_{(\mu\nu)} \equiv \frac{1}{2}\left(X_{\mu\nu}+ X_{\nu\mu}\right) \,, \qquad
 X_{[\mu\nu]} \equiv \frac{1}{2}\left(X_{\mu\nu}-X_{\nu\mu}\right)\,.$$ Up to second order, from the transformation (\[EqTmetric\]) and the decomposition (\[MetricADM\]), we obtain the gauge transformation for the ADM variables as

\[TrulemetricADM\] $$\begin{aligned}
 \tilde{\alpha}
 &= \alpha - {{\cal H}}T - T' + \frac{1}{2} ({{\cal H}}^2 + {{\cal H}}') T^2
 + {{\cal H}}(2 T' T + T_{,i} L^i) - {{\cal H}}\alpha T \notag\\
 & \qquad - \alpha T' - \alpha' T - \alpha_{,i} L^i + \beta^i T_{,i}
 + T'' T + T'{}^2 + T'_{,i} L^i + \frac{1}{2} T^{, i} T_{, i} \,, \\
 \tilde{\beta}^i
 &= \beta^i + T^{,i} - L^i{}' + 2 \alpha T^{, i} - \beta^i T'
 - \beta^i{}' T + \beta^j L^i{}_{,j} - \beta^i{}_{,j} L^j - 2 h^{i j} T_{,j}
 \notag\\
 & \qquad - (2 T' T^{, i} + T T'{}^{, i}) + T L^i{}'' + T' L^i{}'
 + T_{,j} L^{i ,j} - T^{, i}{}_{,j} L^j + L^i{}_{,j}{}' L^j \,, \\
 2 \tilde{h}_{i j} &= 2 h_{i j} - 2 {{\cal H}}T \delta_{i j} - 2 L_{(i ,j)}
 - 4 {{\cal H}}T h_{i j} + (2 {{\cal H}}^2 + {{\cal H}}') T^2 \delta_{i j}
 + 2 {{\cal H}}( T T' + T_{,k} L^k) \delta_{i j} \notag\\
 & \qquad
 + 4 {{\cal H}}T L_{(i ,j)} - 2 \beta_{(i} T_{,j)} - 2 h_{i j}{}' T
 - 2 h_{i j, k} L^k - 4 h_{(i| k} L^k{}_{,|j)} \notag\\
 & \qquad
 - T_{,i} T_{,j} + 2 T L_{(i ,j)}{}'
 + 2 T_{(,i} L_{j)}{}' + 2 L_{(i ,j) k} L^k
 + 2 L_{(i|, k} L^k{}_{,|j)} + L_{k, i} L^k{}_{, j} \,.\end{aligned}$$

These relations must be understood as follows. $a^2 \alpha$ at first order is the $00$ component of the first order metric in the $c$ coordinates system, and taken at the point of coordinates ${\bm x}$, and is thus equal to $-\frac{1}{2}g^{(1)}_{00}\circ c({\bm x})$. Instead, $a^2 \tilde \alpha$ which means $a^2 \tilde \alpha({\bm x})$, when considered at first order is the $\tilde 0 \tilde 0$ component of the first order metric in the $\tilde c$ coordinates system, but also taken at the point of coordinates ${\bm x}$, that is $-\frac{1}{2}g^{(1)}_{\tilde 0 \tilde
  0}\circ \tilde c{({\bm x})}$.

Tangent space basis and tetrads
-------------------------------

As discussed in Section \[ssec:f\], the transformation of the basis on the tangent space is entirely linked to the coordinates change in the base manifold. Usually, the canonical basis ${{\partial}}/{{\partial}}x^\mu$ and the corresponding forms ${{\rm d}}x^\mu$ are used as a basis of the tangent space and they transform according to $$\label{TruleCanonicalBasis}
\frac{{{\partial}}}{{{\partial}}\tilde x^\mu}\circ \tilde c{(\tilde{\bm x})}
 = \frac{{{\partial}}x^\nu}{{{\partial}}\tilde x^\mu} \frac{{{\partial}}}{{{\partial}}x^\nu}{({\bm x})} \,,
 \qquad {{\rm d}}\tilde x^\mu \circ \tilde c {(\tilde{\bm x})}
 = \frac{{{\partial}}\tilde x^\mu}{{{\partial}}x^\nu} {{\rm d}}x^\nu {c({\bm x})} \,.$$ When we consider the components of the metric, these components refer to this canonical basis. However we will use the tetrad field as a basis for the tangent space. We thus need to relate $\tilde{\bm
  e}_{(a)}\circ \tilde c{(\tilde{\bm x})}$ with ${\bm  e}_{(a)}\circ c{({\bm x})}$ in a similar fashion. Since this is a relation between two orthonormal basis at the same point of space-time, there exists a Lorentz transformation $\Lambda^{(a)}_{~~ (b)}$ such that we can relate the two tetrad fields associated with the coordinates systems $c$ and $\tilde c$ (see the bottom right part of Fig. \[fig1\] for an illustration of this) as $$\tilde{\bm e}_{(a)}\circ \tilde c{(\tilde{\bm x})}
 = \Lambda_{(a)}^{~~ (b)} ({\bm x}) {\bm e}_{(b)} ({\bm x}) \,,\qquad
  \tilde{\bm e}^{(a)} \circ \tilde c{(\tilde{\bm x})}
 = \Lambda^{(a)}_{~~ (b)} ({\bm x}) {\bm e}^{(b)} ({\bm x})\,,\qquad
 \eta_{(c) (d)} \Lambda^{(c)}_{~~ (a)} \Lambda^{(d)}_{~~ (b)} = \eta_{(a) (b)} \,.$$ The components of this Lorentz transformation are given, up to first order, by $$\label{Reltetradssamephysicalpoint}
 \Lambda^{(0)}_{~~ (0)} = 1 \,,\qquad
 \Lambda^{(0)}_{~~ (i)} = \Lambda^{(i)}_{~~ (0)} = {{\partial}}_i T \,,\qquad
 \Lambda^{(i)}_{~~ (j)} = \delta^{(i)}_{~~ (j)} +L^{[i}_{~~ , j]}\,.$$ At second order, it proves more useful to relate the components of these quantities, that is, to relate $\widetilde{e}_{(a)}{}^{\widetilde{\mu}}\circ \tilde c{(\tilde{\bm x})}
 \equiv {{\rm d}}\tilde x^\mu\circ \tilde c{(\tilde{\bm    x})}
 [\widetilde{\bm e}_{(a)}]$ to $
 e_{(a)}{}^\mu({\bm x})\equiv {{\rm d}}x^\mu({\bm x})[{\bm e}_{(a)}]$. For the former components, we must use  with $a(\tilde{\eta})$, $\tilde \alpha(\tilde{\bm x})$, $\tilde \beta_i(\tilde{\bm x})$ and $\tilde h_{ij}(\tilde{\bm x})$ \[that is $g_{\tilde \mu \tilde \nu}\circ \tilde c{(\tilde{\bm x})}$\] which can be deduced from the rule  just by shifting the argument of the left hand side of the rules  from ${\bm x}$ to $\tilde{\bm x}$. For the latter we must use  with $a(\eta)$, $\alpha({\bm x})$, $\beta_i({\bm x})$ and $h_{ij}({\bm x})$ \[that is $g_{\mu \nu}{({\bm
     x})}$\]. We do not report the corresponding expression since we will work instead directly with the perturbation components of the metric in the next section.

Momentum, energy $q$ and direction $n^{(i)}$
--------------------------------------------

The components of the momentum of a particle in the canonical basis transform as $$\label{Transfopmu}
 p^{\widetilde{\mu}}\circ \tilde c(\tilde{\bm x})
 = \frac{{{\partial}}\tilde{x}^\mu}{{{\partial}}x^\nu} p^\nu({\bm x})\,, \qquad
 \text{with} \qquad p^{\widetilde{\mu}} \equiv {{\rm d}}\tilde
 x^\mu({\bm p})\,,\quad p^{\mu} \equiv {{\rm d}}x^\mu({\bm p})\,.$$ However, we are going to use the tetrad basis in the tangent space rather than the canonical basis. We thus want to express $p^{\widetilde{(a)}}\circ \tilde c{(\tilde{\bm x})}
 \equiv \widetilde{\bm e}^{{(a)}}\circ \tilde c{(\tilde{\bm x})}[{\bm p}]$, \[or $\tilde q\circ \tilde c{(\tilde{ \bm x})}$ [^3] and $n^{\widetilde{(\imath)}}\circ \tilde c{(\tilde{ \bm x})}$\] as a function of $p^{(a)}({\bm x}) \equiv {\bm e}^{{(a)}} ({\bm x})[{\bm p}]$ \[or $q({\bm x})$ and $n^{(i)}({\bm x})$\]. Eventually, we will prefer to use $q$ and $n^{(i)}$ rather than $p^{(i)}$. From the definition of $q$, Eq (\[def:q\]) we then obtain the desired transformation relation $$\begin{aligned}
 \tilde{q}\circ \tilde c{(\tilde{\bm x})}
 &\equiv q ({\bm x})+\delta q ({\bm x})
 \equiv q ({\bm x})[1+\delta \ln q ({\bm x})] \,,
 \end{aligned}$$ as $$\begin{aligned}
\label{EqTp}
 \tilde{q}\circ \tilde c{(\tilde{\bm x})}
 &\equiv a^2(\tilde \eta) \tilde{N} p^{\tilde{0}}
 \circ \tilde c{(\tilde{\bm x})} \notag\\
 &= q ({\bm x}) \left[ 1 + {{\cal H}}T + T_{,i}n^{(i)}
 + \frac{1}{2} ({{\cal H}}' + {{\cal H}}^2) T^2 + {{\cal H}}T T_{,i}n^{(i)}
 - T' T_{,i}n^{(i)} + \frac{1}{2}  T^{,i} T_{,i}
 + T_{,i} \left( \alpha n^{(i)} - h^i{}_j n^{(j)} \right)\right] \,,\end{aligned}$$ As for $n^{(i)}$ from Eq (\[def:n\^i\]), its transformation $$\begin{aligned}
 n^{\widetilde{(i)}}\circ \tilde c{(\tilde{\bm x})}
 &\equiv  n^{(i)}({\bm x}) + \delta n^{(i)}({\bm x}) \,,\end{aligned}$$ is given by $$\begin{aligned}
 n^{\widetilde{(i)}}\circ \tilde c{(\tilde{\bm x})}
 &\equiv \left[ \tilde{\beta}^i + \left( \frac{1}{\tilde{N}} \delta^i{}_j
 + \tilde{h}^i{}_j \right) \frac{p^{\tilde{j}}}{p^{\tilde{0}}} \right]
 \circ \tilde c{(\tilde{\bm x})} \notag\\
 &= n^{(i)}({\bm x}) + S^{(i)(j)} T_{,j} + L^{[i,j]}  n_{(j)} \,.
 \end{aligned}$$ Here when the argument is not specified, it is $({\bm x})$. Note that for future use, we have defined in these expressions the differences $\delta q$, $\delta \ln q$, and $\delta n^{(i)}$. The first order and second order perturbation of these can be read directly from the expressions above. We also define $\delta q^{(i)}\equiv q\Bigl[(\delta \ln q) n^{(i)}
 + \delta n^{(i)}\Bigr]$.

It is worth stressing that these differences are measured at the same point. For instance $\delta n^{(i)}({\bm x})=  n^{\widetilde{(i)}} \circ
 \tilde c{(\tilde{\bm x})} - n^{(i)}({\bm x}) $, and they do not vanish because in one case we use the tetrads ${\bm e}^{(i)}$ associated with the coordinate system $c$ to obtain the components, and in another case we use the tetrads $\tilde{\bm e}^{(i)}$ associated with the coordinate system $\tilde c$. It is the basis at a given point of space-time that changes when we change the coordinate system, not the momentum itself. This point is illustrated in the bottom right part of Fig. (\[fig1\])

Scalar distribution function
----------------------------

If we were using the natural basis associated with a coordinate system (the canonical basis) for the tangent space, then any scalar function on the tangent bundle $T{\cal M}$, that is a function of the space-time position and of the tangent space at each point, would transform as $$I \circ \tilde c{(\tilde {\bm x},p^{\widetilde{\mu}})} = I({\bm x},p^\mu) \,,$$ where $\tilde {\bm x}$ and ${\bm x}$ are related by  and $p^{\widetilde{\mu}}$ and $p^{\mu}$ are related by . Again this rule is a statement that the function is invariant under a change of coordinates because it is defined purely geometrically.

However, as mentioned earlier, we use the basis of the tetrad field, ${\bm e}_{(i)}$ and ${\bm e}^{(i)}$, to obtain the components of momentum not the canonical basis. The tetrads are also completely determined by the choice of coordinates due to our prescription . The tetrad field, though being of tensorial nature, is not invariant as in Eq. . Furthermore, as mentioned earlier, we also work with the conformal momentum ${\bm q}$ rather than the momentum itself ${\bm p}$. Given this choice for the basis of the tangent space, the scalar function transforms as $$\label{Trulef}
 I\circ \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 = I ({\bm x},q^{(\imath)})\,,$$ that is, it is unchanged when it is evaluated at the same point of the tangent bundle. On the other hand, the gauge transformation is a transformation rule at the same coordinate point and it is the relation between $I\circ \tilde c{({\bm x},q^{{(i)}})}$ and $I{( {\bm x},q^{{(i)}})}$. We then need the expressions of $q^{\widetilde{(\imath)}}\circ \tilde c{(\tilde{\bm x})}$ in terms of $q^{{(\imath)}}({\bm x})$ in spherical coordinates, which are derived in the previous section in . At first order, using the fact that the background distribution function cannot depend on the direction $n^{(i)}$, we obtain $$I \circ \tilde c{({\bm x},q,n^{(i)})}
 + \left(\xi^\mu{{\partial}}_\mu +\delta q \frac{{{\partial}}}{{{\partial}}q}\right) I \circ
 \tilde c{({\bm x},q,n^{(i)})}= I ({\bm x},q,n^{(i)})\,.$$ Given that the background distribution function also depends neither on time nor on space but only on $q$, we obtain $$I^{(1)}\circ \tilde c{({\bm x},q,n^{(i)})}
 = I^{(1)} ({\bm x},q,n^{(i)})-\delta \ln q \,{{\frac{\partial \bar I{(q)}}{\partial \ln q}}}\,.$$ At second order, we obtain $$I \circ \tilde c{({\bm x},q,n^{(i)})}
 + \left(\xi^\mu{{\partial}}_\mu +\delta q \frac{{{\partial}}}{{{\partial}}q} +\delta n^{(i)} D_{(i)}
 +\frac{1}{2}\delta q^2\frac{{{\partial}}^2 }{{{\partial}}q^2}\right) I \circ
 \tilde c{({\bm x},q,n^{(i)})}= I ({\bm x},q,n^{(i)}) \,.$$ Using the first order expressions, the gauge transformation rule reads $$\begin{aligned}
\label{GaugeTRuleI2}
 \frac{1}{2}I^{(2)}\circ \tilde c{({\bm x},q,n^{(i)})}
 &=& \frac{1}{2} I^{(2)} ({\bm x},q,n^{(i)})
 - \left(\xi^\mu{{\partial}}_\mu +\delta \ln q {{\frac{\partial}{\partial \ln q}}}+\delta n^{(i)} D_{(i)}\right)
 I^{(1)}{({\bm x},q,n^{(i)})}\nonumber\\*
 &&+ \left(\xi^\mu{{\partial}}_\mu +\delta n^{(i)}\frac{{{\partial}}}{{{\partial}}n^{(i)}}\right)
 (\delta \ln q) {{\frac{\partial}{\partial \ln q}}}\bar I{(q)}
 +\frac{1}{2}(\delta \ln q)^2 \left({{\cal D}}_q^2-2 {{\frac{\partial}{\partial \ln q}}}\right)\bar I(q)\,,\end{aligned}$$ where we used ${{\partial}}\delta \ln q / {{\partial}}\ln q = 0$.

This transformation rule can be applied to $I$ or $V$ since these are scalar valued distribution functions. However we are interested in the transformation rule for the spectral components of $I$. Using the decomposition Eq (\[intensity\]) for $I$ we obtain that the temperature is transforming under a gauge transformation as (noting for simplicity $\tilde \Theta \equiv \Theta \circ \tilde c$) $$\begin{aligned}
 \tilde{\Theta}^{(1)} &= \Theta^{(1)} + (\delta \ln q)^{(1)} \,,\\
 \frac{1}{2} \tilde{\Theta}^{(2)} &= \frac{1}{2} \Theta^{(2)}
 + \frac{1}{2} (\delta \ln q)^{(2)} +\Theta^{(1)} (\delta \ln q)^{(1)}
 - \left(\xi^\mu \frac{{{\partial}}}{{{\partial}}x^\mu}+\delta n^{(i)} D_{(i)} \right)
 \Bigl[ \Theta^{(1)}+(\delta \ln q)^{(1)} \Bigr]\,,\end{aligned}$$ where it is implied that all quantities are evaluated either at ${\bm x}$ or at $({\bm x},q,n^{(i)})$. The detailed form of the transformation rule can then be obtained just by considering the perturbations of $q$ and $n^i$, $\delta \ln q$ and $\delta n^{(i)}$, which have been obtained in Eqs. (\[EqTp\]). For completeness we report it here $$\begin{aligned}
 \tilde{\Theta}
 &= \Theta + {{\cal H}}T + T_{,i}n^{(i)} + \frac{1}{2} ({{\cal H}}^2 - {{\cal H}}') T^2
 + {{\cal H}}T (T_{, i} n^{(i)} - T') - (T' T_{,i} + T T_{, i}') n^{(i)} \notag\\
 &\qquad
 + T_{,i} \Bigl( \alpha n^{(i)} - h^i{}_j n^{(j)} \Bigr) - L^i ({{\cal H}}T_{, i}
 + T_{, i j} n^{(j)}) + \left( n^{(i)}n^{(j)} - \frac{1}{2} \delta^{i j}
 \right) T_{,j} T_{, i} - L^{[i}{}_{,j]} T_{, i} n^{(j)} \notag\\
 &\qquad
 + ({{\cal H}}T + T_{, i} n^{(i)}) \Theta - \xi^\mu \Theta_{, \mu}
 - \Bigl( S^{(i) (j)} T_{,j} + L^{[i}{}_{,j]} n^{(j)} \Bigr) D_{(i)} \Theta \,.\end{aligned}$$ Finally, the gauge transformation of $y$ is trivial. Since $y$ vanishes at first order, $y$ is gauge invariant at second order, and it can also be checked directly by extracting $y$ out of the transformation rule of $I$ at second order .

Baryons fluid description
-------------------------

In order to obtain the complete gauge transformation of the collision term, we need the gauge transformation rule of the baryons fluid velocity in the tetrad frame up to second order, since it appears in the collision term. We obtain $$\begin{aligned}
 v^{\widetilde{(\imath)}}\circ \tilde c{({\bm x})}
 &= v^{(i)}({\bm x}) + T^{, i} + \alpha T^{, i} - h^i{}_j T^{, j} - T' T^{, i}
 + L^{[i}{}_{, j]} (v^{(j)} + T^{, j})
 - T (v^{(i)}{}' + T^{, i}{}') - L^j (v^{(i)} + T^{, i})_{, j} \,.\end{aligned}$$ We also need the gauge transformation rule up to first order of the electrons density and it is easily obtained to be $$\delta_e\circ \tilde c{({\bm x})}
 = \delta_e({\bm  x}) + 3 {\cal H} T n_e \,.$$

Gauge dependence of the Boltzmann equation and Gauge invariant form {#ssec:lc}
-------------------------------------------------------------------

Having derived all the necessary gauge transformation rules, it is now possible to check the gauge dependence of the derived second order Boltzmann equation and collision term explicitly. More precisely we shall check that the Liouville and the collision terms of the Boltzmann equation transform as they should do. The Liouville term and the Collision term are distribution functions (scalar or tensor valued depending whether or not we are considering the intensity or polarisation). From the transformation rule of a scalar distribution function  and the spectral decomposition (\[spectraldecI\]) we deduce that ${\cal L}^Y$ and ${{\cal C}}^Y$ are gauge invariant and ${{\cal L}}^\Theta$ and ${{\cal C}}^\Theta$ should transforms up to the second order as (noting $\widetilde{{{\cal L}}^\Theta} \equiv {{\cal L}}^\Theta \circ \tilde c$ and $\widetilde{{{\cal C}}^\Theta} \equiv {{\cal C}}^\Theta \circ \tilde c$ ) $$\begin{aligned}
 \widetilde{{{\cal L}}}^\Theta
 &= {{\cal L}}^\Theta - \Bigl( \xi^\mu {{\partial}}_\mu + \delta n^{(i)} D_{(i)} \Bigr) {{\cal L}}^{\Theta\,(1)}
 + 2 {{\cal H}}T {{\cal L}}^{\Theta\,(1)} \,, \\
 \widetilde{{{\cal C}}^\Theta}
 &= {{\cal C}}^\Theta - \Bigl( \xi^\mu {{\partial}}_\mu + \delta n^{(i)} D_{(i)} \Bigr) {{\cal C}}^{\Theta (1)}
 + 2 {{\cal H}}T {{\cal C}}^{\Theta (1)} \,.
 \end{aligned}$$

After very long and tedious but straightforward calculations using all the transformation rules derived so far for the distribution function, the metric components and the baryons velocity and energy density, we have checked that the Liouville operator ${{\cal L}}^\Theta$ and the collision term ${{\cal C}}^\Theta$ actually transform as in the above equations. This completes the consistency test of the Boltzmann equation that we have derived as well as all the gauge transformation rules obtained for its constituents.

The same property is found of course for the circular polarisation since in that case the collision term vanishes. As for polarisation, we have also checked that the Liouville operator ${{\cal L}}^P_{(i)(j)}$ and the collision term ${{\cal C}}^P_{(i)(j)}$ transform like tensor valued quantities (see the details in Appendix \[AppGTLC\]). More importantly, we have checked that the Liouville and collision terms for the spectral distortions (${{\cal L}}^Y_{(i)(j)}$ and ${{\cal C}}^Y_{(i)(j)}$) are gauge invariant as it should be since they vanish on the background and first-order spacetimes.

The gauge invariance of the Boltzmann equation, as in the case of Einstein equation and in general for covariant equations, enables us to write it down in terms of gauge invariant variables. In practice, it is equivalent to completely fix the gauge and write down the equations in term of the perturbation in this gauge.

Summary and discussion {#sec:conc}
======================

In this paper we derived the second order Boltzmann equation in the most general manner incorporating polarisation and without fixing a gauge. In order to describe the polarisation of photon, we used a formalism based on a tensor-valued distribution function. We performed the separation between temperature and spectral distortion for the intensity and we also extended this separation to polarisation.

We then derived the gauge transformation rules for the metric, the momentum and the distribution function to see how those quantities are mixed under the gauge transformation. As an application, we checked the gauge dependence of the derived Boltzmann equation under a gauge transformation and obtained consistent transformation rules. This is a non-trivial check of the correctness of the derived equations as well as the gauge transformation rules.

We now discuss two issues related to the gauge dependence in the Boltzmann equation.

Gauge dependence of lensing term
--------------------------------

It is well known that the lensing term $$\begin{aligned}
 {{\cal L}}\supset \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} D_i \Theta
 ~ : ~ ({{\mathrm{lensing ~ term}}}) \,,\end{aligned}$$ which is written in terms of the conventional lensing potential in the Newtonian gauge significantly affects the bispectrum of CMB [@Hanson:2009kg; @Lewis:2011fk; @Ade:2013ydc]. Indeed, the correlation between the lensing and ISW effect is the dominant contribution to the bias for the local-type non-Gaussianity in Planck. However it is very hard to include this contribution in the line of sight integration and evaluate it until today. Usually, the lensing effect is added separately to the final result obtained in the Poisson gauge. However the effect from this lensing term depends on the gauge choice. Actually, as we have seen above, the lensing term is mixed with other terms under the gauge transformation. This means that some of lensing effects in a specific gauge are absorbed into other effects in another gauge. In principle, there exists a gauge where we can avoid the difficult computation of this lensing term to some extent by evaluating other more tractable terms. Since we have derived the Boltzmann equation without choosing any specific gauge, it should be possible to investigate this possibility further.

Observed temperature anisotropies
---------------------------------

Here we make a comment on the observed temperature anisotropies. In the main part of this paper, we have shown that the second order Boltzmann equation is gauge invariant and thus it can be written in terms of gauge invariant quantities. However, there is a subtlety in the meaning of “gauge invariance”. This originates from our choice of the local inertial frame.

As is clear from Eq (\[tetrads\]), we always choose the local inertial frame so that the three-velocity vanishes, $\hat{v}^i=0$ and there is no rotation of the spatial axis relative to the background spatial coordinate axis, let us call $\theta_i=0$. In order to achieve this, the local inertial frame has to be changed when we perform a gauge transformation. If we were to identify this local inertial frame as the one of an observer, we would be lead to consider different observers in different gauges. This is clear from the gauge transformation at the first order: $$\Theta \to \Theta + {\cal H} T + T_{,i} n^i\,.
\label{1sttempg}$$ The last term comes from a change of local inertial frame. By fixing the gauge we can promote $\Theta$ to the gauge invariant temperature fluctuations but these are temperature fluctuations observed by an observer with $\hat{v}^i=0$ and $\theta_i=0$ in this gauge. In order to evaluate temperature fluctuations observed by a different observer, we need to change the local inertial frame. Alternatively, we can perform a gauge transformation keeping the conditions $\hat{v}^i=0$ and $\theta_i=0$ for the local inertial frame so that this frame coincides with the one of the observer.

In all the literature, the second order temperature anisotropies are calculated in the Poisson gauge so far, with a specific choice of the local inertial frame. Strictly speaking this is not the temperature anisotropies that we observe as there is no reason for us to be comoving with the local inertial frame associated with such a gauge. One thus needs to change the local inertial frame or change a gauge. At first order, this was not an issue. As is clear from (\[1sttempg\]), the change of gauge and the local inertial frame only affect the monopole $\ell =0$ and dipole $\ell =1$ if we expand the temperature anisotropies into multipole components. Thus, the $\ell \geq 2$ modes are not affected by the change of observers. However this is no longer the case at the second order. In the second order gauge transformation, there are terms that are convolutions of the first order temperature anisotropies and the gauge transformation; $$\begin{aligned}
 \Theta \to \Theta - \xi^\mu \Theta_{, \mu} - \delta n^{(i)} D_{(i)} \Theta
 + \cdots \,.
 \end{aligned}$$ These terms affect the observed temperatures even for the $\ell \geq 2$ modes.

In order to define the "observed temperature anisotropies”, we should keep the conditions $\hat{v}^i=0$ and $\theta_i=0$ for the local inertial frame and specify the gauge so that this local inertial frame coincides with our local inertial frame where we perform experiments. Thus special care must be taken when we compare theoretical predictions to observations. Our formula for the gauge transformation will be useful to investigate this issue further.

A.N is grateful to Shuichiro Yokoyama and Ryo Saito for their continuous encouragement and fruitful discussion. This work was supported in part by Monbukagaku-sho Grant-in-Aid for the Global COE programs, "The Next Generation of Physics, Spun from Universality and Emergence” at Kyoto University, by JSPS Grant-in-Aid for Scientific Research (A) No. 21244033, and by the Long-term Workshop at Yukawa Institute on Gravity and Cosmology 2012, YITP-T-12-03. AN is partly supported by Grant-in-Aid for JSPS Fellows No. 21-1899 and JSPS Postdoctoral Fellowships for Research Abroad. CP was supported by the STFC (UK) grant ST/H002774/1 during the first part of this research, and was then supported by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. K.K. is supported by STFC grant ST/H002774/1, ST/K0090X/1, the European Research Council and the Leverhulme trust.

ADM variables and usual perturbation variables {#app:ADM}
==============================================

From the form of the metric in the ADM parametrisation (\[MetricADM\]), and the perturbation of the lapse function and the spatial metric given in the equations (\[DefPertNNi\]), the perturbed metric is expressed up to second order as $${{\rm d}}s^2 = a^2 (\eta) \Bigl[ - (1 + 2 \alpha + \alpha^2 - \beta_i \beta^i)
 {{\rm d}}\eta^2 + 2 (\delta_{i j} + 2 h_{i j}) \beta^j {{\rm d}}x^i {{\rm d}}\eta
 + (\delta_{i j} + 2 h_{i j}) {{\rm d}}x^i {{\rm d}}x^j \Bigr] \,.$$ This has to be compared with the usual parametrisation of the perturbations of the metric which is in the form $$\begin{gathered}
 {{\rm d}}s^2 = a^2 (\eta) \Bigl[ - (1 + 2 A) {{\rm d}}\eta^2 + 2 B_i {{\rm d}}x^i {{\rm d}}\eta
 + (\delta_{i j} + 2 C_{i j}) {{\rm d}}x^i {{\rm d}}x^j \Bigr] \,,
 \end{gathered}$$ where $C_{i j}$ can be further split into 2 scalar, 2 vector and 2 tensor degrees of freedom. By a direct comparison of these two parametrisations, the relation between the two metric parametrisations is

$$\begin{aligned}
 g_{0 0} &: 1 + 2 \alpha + \alpha^2 - \beta_i \beta^i = 1 + 2 A\,,\\
 g_{0 i} &: \beta_i + 2 h_{i j} \beta^j = B_i\,, \\
 g_{i j} &: \delta_{i j} + 2 h_{i j} = \delta_{i j} + 2 C_{i j}\,.\end{aligned}$$

At first order we obtain $$\begin{gathered}
 A^{(1)} = \alpha^{(1)} \,, \qquad B_i^{(1)} = \beta^{(1)}_i \,, \qquad
 C_{i j}^{(1)} = h_{i j}^{(1)} \,,\end{gathered}$$ and the two parametrisations are the same. However, at second order we get the relations $$\begin{gathered}
 A^{(2)}=  \alpha^{(2)} + \alpha^{(1) 2} - \beta^{(1)}_i \beta^{(1) i}  \,, \qquad
 B_i^{(2)} = \beta^{(2)}_i + 4 h_{i j}^{(1)} \beta^{(1) j}  \,,\qquad
 C_{i j}^{(2)} = h_{i j}^{(2)} \,.\end{gathered}$$

Construction of the distribution function for polarised light {#app:Dist}
=============================================================

We consider a two-dimensional polarisation plane defined by two unit complex vectors $\hat {{\boldsymbol{\epsilon}}}_{({{\mathrm{I}}})}$ and $\hat {{\boldsymbol{\epsilon}}}_{({{\mathrm{II}}})}$, which are mutually orthogonal, $\hat \epsilon_{({{\mathrm{A}}})}{}^{\star\mu} \hat \epsilon_{({{\mathrm{B}}})}{}^\nu g_{\mu\nu}
 =\delta_{({{\mathrm{A}}})({{\mathrm{B}}})}$. Any polarisation $\boldsymbol{\epsilon}$ can be represented by a superposition of $\hat {{\boldsymbol{\epsilon}}}_{(A)}$ in the form $${{\boldsymbol{\epsilon}}}= \sum_{A={{\mathrm{I}}},{{\mathrm{II}}}} \epsilon^A \hat {{\boldsymbol{\epsilon}}}_{(A)}
 \label{pol}\,.$$ The orthogonal vectors $\hat {{\boldsymbol{\epsilon}}}_{(A)}$ define a polarisation plane and we choose them to be orthogonal to the direction of the photon $n^{(i)}$ and to the observer velocity ${\bm e}_{(0)}$, $$\label{EqProjepsilons}
 \hat \epsilon_{(A)}{}^\mu n^{\nu} g_{\mu \nu}
 = \hat \epsilon_{(A)}{}^\mu {e}_{(0)}{}^\nu g_{\mu\nu}=0\,.$$ We can also associate canonically polarisation forms through $\hat \epsilon^{(A)}{}_\mu \equiv g_{\mu\nu} \hat \epsilon_{(A)}{}^\nu$, and they will be also complex unit forms and mutually orthogonal. The polarisation density matrix $f_{A B}$ is defined so that the expected number of photon in a phase-space element with a polarisation state ${{\boldsymbol{\epsilon}}}$ is given by $$\begin{gathered}
 f( {\bf x}, {\bf p}, {{\boldsymbol{\epsilon}}})
 \equiv f_{A B}({\bf x}, {\bf p}) \epsilon^{\star A} \epsilon^B \,.
 \end{gathered}$$ With such a parametrisation, all the electromagnetic gauge degrees of freedom have been fixed and we parametrise the physical degrees of freedom of this density matrix by the usual Stokes parameters as $$\begin{gathered}
 f_{A B} = \frac{1}{2}
\begin{pmatrix}
 I + Q & U - i V \\
 U + i V & I - Q
\end{pmatrix} \,,\end{gathered}$$ where it is implied that $f_{AB}$ and the Stokes parameters depend on the position $x^\mu$ and on the momentum $p^\mu$ \[or $(q,n^{(i)})$ in spherical coordinates\]. $f_{AB}$ is a Hermitian matrix since $f_{AB}=f_{BA}^\star$.

From the four-dimensional point of view, the polarisation density matrix is a tensor-valued distribution function. It is a 2-form defined by $$\begin{gathered}
 f_{\mu \nu} \equiv f_{A B} \epsilon^{\star(A)}{}_\mu \epsilon^{(B)}{}_\nu \,,\end{gathered}$$ and the expected number of photon in a phase-space element for a polarisation state ${{\boldsymbol{\epsilon}}}$ is given by $$\label{Eqdefffromfmunu}
 f ({\bf x}, {\bf p}, {{\boldsymbol{\epsilon}}})
 \equiv f_{\mu\nu} ({\bf x}, {\bf p}) \epsilon^{\star\mu} \epsilon^\nu \,.$$ This can be viewed as a multipolar expansion in the polarisation state. From (\[EqProjepsilons\]), the tensor-valued distribution function is a projected quantity such that $$f_{\mu \nu} =S_\mu{}^\alpha S_{\mu}{}^\beta f_{\alpha \beta}\,.$$ It is then straightforward to realize that it can be decomposed according to (\[Pol\]).

Boltzmann equation for the tensor-valued distribution function {#app:Boltz_pol}
==============================================================

From a scalar valued to a tensor-valued distribution function
-------------------------------------------------------------

In this section, we explain in detail how the Boltzmann equation for the tensor-valued distribution function can be obtained from the Boltzmann equation of a scalar distribution function. Since this scalar distribution function $f$ depends on $x^\mu,p^\mu$ but also on $\epsilon^\mu$, the action of the Liouville operator is given by $$\label{CrudeBoltzmann}
 \frac{{{\cal D}}}{{{\cal D}}\lambda} f (x^\mu, p^\mu, \epsilon^\mu)
 = \frac{{{\rm d}}x^\alpha}{{{\rm d}}\lambda} \frac{{{\partial}}f}{{{\partial}}x^\alpha}
 + \frac{{{\rm d}}p^\alpha}{{{\rm d}}\lambda} \frac{{{\partial}}f}{{{\partial}}p^\alpha}
 + \frac{{{\rm d}}\epsilon^\alpha}{{{\rm d}}\lambda} \frac{{{\partial}}f}{{{\partial}}\epsilon^\alpha}
 = C [f]\,.$$ In the geometric optics approximation, $p^\mu$ and $\epsilon^\mu$ are parallel transported and we obtain $$\begin{aligned}
 0 =& \frac{{{\cal D}}p^\alpha}{{{\cal D}}\lambda}
 = \frac{{{\rm d}}p^\alpha}{{{\rm d}}\lambda}
 + \Gamma^\alpha_{\beta \gamma} p^\beta p^\gamma \,, \\
 0 =& \frac{{{\cal D}}\epsilon^\alpha}{{{\cal D}}\lambda}
 = \frac{{{\rm d}}\epsilon^\alpha}{{{\rm d}}\lambda}
 + \Gamma^\alpha_{\beta \gamma} \epsilon^\beta p^\gamma \,.\end{aligned}$$ Using (\[Eqdefffromfmunu\]), the term involving the evolution of polarisation is obtained as $$\begin{gathered}
 \frac{{{\rm d}}\epsilon^\alpha}{{{\rm d}}\lambda} \frac{{{\partial}}f}{{{\partial}}\epsilon^\alpha}
 = f_{\alpha \beta} \frac{{{\rm d}}\epsilon^{\alpha}}{{{\rm d}}\lambda} \epsilon^{* \beta}
 + f_{\alpha \beta} \frac{{{\rm d}}\epsilon^{\star \beta}}{{{\rm d}}\lambda} \epsilon^{\alpha}
 = - \Gamma^{\alpha}_{\gamma \delta} \epsilon^\gamma p^\delta f_{\alpha \beta}
 \epsilon^{* \beta}
 - \Gamma^{\beta}_{\gamma \delta} \epsilon^{\star\gamma} p^\delta f_{\alpha \beta}
 \epsilon^{\alpha}\,.\end{gathered}$$ Combining this result with the space-time derivative term of the Liouville operator, we get $$\begin{gathered}
 p^\alpha \frac{{{\partial}}f}{{{\partial}}x^\alpha} + \frac{{{\rm d}}\epsilon^\alpha}{{{\rm d}}\lambda}
 \frac{{{\partial}}f}{{{\partial}}\epsilon^\alpha}
 = p^\gamma \epsilon^\alpha (\nabla_\gamma f_{\alpha \beta}) \epsilon^{* \beta} \,,\end{gathered}$$ and thus the Boltzmann equation (\[CrudeBoltzmann\]) can be rewritten as $$\label{EqBoltzmannfee}
 \frac{{{\cal D}}f}{{{\cal D}}\lambda}
 = \epsilon^\mu \left( p^\alpha \nabla_\alpha f_{\mu \nu}
 + \frac{{{\rm d}}p^\alpha}{{{\rm d}}\lambda} \frac{{{\partial}}f_{\mu \nu}}{{{\partial}}p^\alpha} \right)
 \epsilon^{* \nu}
 = C [f] \equiv \epsilon^\mu C_{\mu \nu} \epsilon^{* \nu} \,.$$ The last equality defines the tensor-valued collision term. If we do fix the electromagnetic gauge condition for the collision term in the same manner as what we did for $f_{\mu \nu}$, that is, if $C_{\mu \nu} =S_\mu{}^\alpha S_{\nu}{}^\beta C_{\alpha \beta}$, then the Boltzmann equation for $f_{\mu\nu}$ is given by $$\label{EqBaseBoltzmann}
 S_\mu{}^\alpha S_{\mu}{}^\beta \frac{{{\cal D}}f_{\alpha \beta}}{{{\cal D}}\lambda}
 \equiv S_\mu{}^\alpha S_{\mu}{}^\beta \left(p^\sigma \nabla_\sigma f_{\alpha \beta}
 + \frac{{{\rm d}}p^\sigma}{{{\rm d}}\lambda}
 \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}p^\sigma}\right)
 = C_{\mu\nu}\,.$$ Note that the use of the projectors is required because the components of the equation which are not in the polarisation plane are not fixed by (\[EqBoltzmannfee\]).

From the canonical basis to the tetrad basis {#sapp:covp}
--------------------------------------------

In the equation (\[EqBaseBoltzmann\]), the Greek indices refer to a given coordinate system and its canonical basis for the tangent space, and the distribution function is a function of $(x^\mu,p^\mu)$. If we want to use instead an orthonormal basis for the tangent space, that is, to use the components $p^{(i)}=e^{(i)}{}_\mu p^\mu$ or the conformal momentum components $q^{(i)} = a p^{(i)}$ in the tetrad basis, then the Boltzmann equation can be modified accordingly. In order to do so, we need to be explicit about the partial derivatives to emphasize which variables are to be kept constant when the partial derivatives are evaluated. The Boltzmann equation reads indeed $$\label{Boltzbaseopen}
S_\mu{}^\alpha S_{\nu}{}^\beta \frac{{{\cal D}}f_{\alpha \beta}}{{{\cal D}}\lambda}
 = S_\mu{}^\alpha S_{\nu}{}^\beta \left( p^\gamma
 \left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}x^\gamma}\right|_{p^\mu}
 - p^\gamma \Gamma^\delta_{\gamma \alpha} f_{\delta \beta}
 - p^\gamma \Gamma^\delta_{\gamma \beta} f_{\alpha \delta}
 + \left.\frac{{{\rm d}}p^\gamma}{{{\rm d}}\lambda}
 \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}p^\gamma}\right|_{x^\mu}\right) \,.$$ Using the properties $$\left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}x^\mu} \right|_{p^\mu}
 = \left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}x^\mu}\right|_{q^{(i)}}
 + \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}q^{(i)}}
 \frac{{{\partial}}(a e^{(i)}{}_\nu)}{{{\partial}}x^\mu}p^\nu \,, \qquad
 \left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}p^\mu}\right|_{x^\mu}
 = \left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}q^{(i)}}\right|_{x^\mu} ae^{(i)}{}_\mu \,,
 \qquad
 \frac{{{\rm d}}q^{(i)}}{{{\rm d}}\lambda}
 = a e^{(i)}{}_\mu \frac{{{\rm d}}p^\mu}{{{\rm d}}\lambda}
 + p^\mu \frac{{{\rm d}}(a e^{(i)}{}_\mu)}{{{\rm d}}\lambda} \,,$$ we then deduce that $$\begin{aligned}
\label{Boltzbasetetrad}
 S_\mu{}^\alpha S_{\nu}{}^\beta \frac{{{\cal D}}f_{\alpha \beta}}{{{\cal D}}\lambda}
 &=S_\mu{}^\alpha S_\nu{}^\beta  \left(p^\gamma
 \left.\frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}x^\gamma} \right|_{q^{(i)}}
 - p^\gamma \Gamma^\delta_{\gamma \alpha} f_{\delta \beta}
 - p^\gamma \Gamma^\delta_{\gamma \beta} f_{\alpha\delta}
 + \left.\frac{{{\rm d}}q^{(i)}}{{{\rm d}}\lambda}
 \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}q^{(i)}} \right|_{x^\mu}\right) \notag\\
 &= S_\mu{}^\alpha S_{\mu}{}^\beta
 \left(p^\gamma \nabla_\gamma f_{\alpha \beta} + \frac{{{\rm d}}q^{(i)}}{{{\rm d}}\lambda}
  \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}q^{(i)}}\right)\,.\end{aligned}$$ Comparing Eqs. (\[Boltzbaseopen\]) and (\[Boltzbasetetrad\]), we notice that the notation $p^\gamma \nabla_\gamma f_{\alpha \beta}$ could be ambiguous. Indeed if we use the canonical coordinate system for the tangent space, then ${{\partial}}/{{\partial}}x^\mu$ is to be taken at $p^\mu$ fixed, but if we take the tetrad basis (or another coordinate system for the tangent space), then ${{\partial}}/{{\partial}}x^\mu$ is to be taken with $q^{(i)}$ fixed.

Eq. (\[Boltzbasetetrad\]) is not exactly the desired form of the Boltzmann equation when the distribution function depends on $(x^\mu,q^{(i)})$. In fact, the use of the tetrad basis makes it natural to work with spherical coordinates in the tangent space. In order to introduce them, we first relate the Cartesian derivative in the tangent space, that is, the derivative with respect to $q^{(i)}$, to the covariant derivative on the unit sphere which is described by the possible directions $n^{(i)}$ of the momentum. We must stress that at any point of space-time, $x^\mu$, a distribution function (tensor-valued like $f_{\mu\nu}$ or scalar valued like its trace $I$) which depends on $(x^\mu,q^{(i)})$ can be considered as a field in the tangent space because the tangent space at a given point can be considered as a flat three-dimensional manifold whose points are labelled by $q^{(i)}$ and the natural covariant derivative in this manifold is ${{\partial}}/{{\partial}}q^{(i)}$. Using that $q^{(i)} = q n^{(i)}$, and the property $$\label{ExtrinsicK}
\frac{{{\partial}}n^{(i)}}{{{\partial}}q^{(j)}} = \frac{1}{q} S^{(i)}{}_{(j)} \,,$$ it is possible to show the following relations; $$\begin{aligned}
\label{Expandderpi}
 q \frac{{{\partial}}T_{\mu\nu}}{{{\partial}}q^{(i)}}
 &= \frac{{{\partial}}T_{\mu\nu}}{{{\partial}}\ln q} n_{(i)} + D_{(i)} T_{\mu\nu}
 - e_{(i)}{}^\rho (T_{\mu \rho} n_{\nu} + T_{\nu \rho}n_{\mu}) \,,
 &\quad\text{with}\quad
 & D_{(i)} T_{\mu\nu} \equiv q S_{(i)}{}^{(j)} S_\mu{}^\alpha S_\nu{}^\beta
 \frac{{{\partial}}T_{\alpha\beta}}{{{\partial}}q^{(j)}} \,, \\*
 q \frac{{{\partial}}S}{{{\partial}}q^{(i)}}
 &= \frac{{{\partial}}S}{{{\partial}}\ln q} n^{(i)} + D_{(i)} S \,,
 &\quad\text{with}\quad
 & D_{(i)} S \equiv q S_{(i)}{}^{(j)} \frac{{{\partial}}S}{{{\partial}}q^{(j)}} \,,\end{aligned}$$ where $T_{\mu\nu}(q^{(i)})$ is a projected tensor field in the tangent space such that $T_{\mu\nu}(q^{(i)}) n^{\nu}=T_{\mu\nu}(q^{(i)}) n^{\nu}=0$, and $S(q^{(i)})$ is a scalar field in the tangent space. Here $D_{(i)}$ is the covariant derivative on the two-sphere associated with the unit direction vector $n^{(i)}$, and it appears naturally as an induced derivative on the sphere, given that this is the surface orthogonal to $n^{(i)}$ (see Ref. [@Gourgoulhon:2007ue] for more details on induced derivatives). We can then deduce the useful property $$\label{KeyRuleForCovDSpherical}
 q S_\mu{}^\alpha S_\nu{}^\beta \frac{{{\partial}}f_{\alpha \beta}}{{{\partial}}q^{(i)}}
 = \frac{{{\partial}}f_{\mu\nu}}{{{\partial}}\ln q} n_{(i)} + D_{(i)} f_{\mu\nu} \,.$$ Given that $$\frac{{{\rm d}}q^{(i)}}{{{\rm d}}\lambda}
 = q \left(\frac{{{\rm d}}\ln q }{{{\rm d}}\lambda} n^{(i)}
 + \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} \right)\,,$$ from Eqs. (\[Expandderpi\]) and Eq. (\[Boltzbasetetrad\]), we then find that the Boltzmann equation takes the form $$S_\mu{}^\alpha S_\nu{}^\beta \frac{{{\cal D}}f_{\alpha \beta}}{{{\cal D}}\lambda}
 = S_\mu{}^\alpha S_\nu{}^\beta \nabla_\gamma f_{\alpha \beta}
 \frac{{{\rm d}}x^\gamma}{{{\rm d}}\lambda}
 + \frac{{{\partial}}f_{\mu \nu}}{{{\partial}}\ln q} \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 + D_{(i)} f_{\mu \nu} \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} = C_{\mu \nu} \,.$$

Decomposition of the Boltzmann equation
---------------------------------------

In order to obtain equations for the components of $f_{\mu\nu}$, $I, P_{\mu \nu}$ and $V$, we want to apply the same type of decomposition on the equation itself. Applying ${{\cal D}}/{{\cal D}}\lambda$ on the decomposition  of $f_{\mu\nu}$ leads to $$\begin{gathered}
 \frac{{{\cal D}}f_{\mu \nu}}{{{\cal D}}\lambda}
 = \frac{1}{2} \left( \frac{{{\cal D}}I}{{{\cal D}}\lambda} S_{\mu \nu}
 + I \frac{{{\cal D}}S_{\mu \nu}}{{{\cal D}}\lambda} \right)
 + \frac{{{\cal D}}P_{\mu \nu}}{{{\cal D}}\lambda}
 + \frac{{{\rm i}}}{2} \epsilon_{\alpha \mu \nu \beta}
 \left( \frac{{{\cal D}}V}{{{\cal D}}\lambda} e_{(0)}{}^\alpha n^\beta
 + V \frac{{{\cal D}}(e_{(0)}{}^\alpha n^\beta)}{{{\cal D}}\lambda} \right) \,.\end{gathered}$$ We then need to screen-project this equation in order to obtain the tensor-valued Boltzmann equation. Then the last two terms vanish. Indeed, first $$S_\mu{}^\alpha S_\nu{}^\beta \frac{{{\cal D}}S_{\alpha \beta}}{{{\cal D}}\lambda}
 = S_\mu{}^\alpha S_\nu{}^\beta \left[ \frac{{{\rm d}}p}{{{\rm d}}\lambda} \left(
 - \frac{p_\alpha e^{(0)}{}_\beta + e^{(0)}{}_\alpha p_\beta}{p^2}
 + 2 \frac{ p_\alpha p_\beta}{p^3} \right)
 + \frac{1}{p} \left( p_\alpha \frac{{{\cal D}}e^{(0)}{}_\beta}{{{\cal D}}\lambda}
 + \frac{{{\cal D}}e^{(0)}{}_\alpha}{{{\cal D}}\lambda} p_\beta \right) \right] = 0\,.$$ Second, from the normalization condition of $e_{(0)}{}^\mu$ and $n^\mu$, we can show that $$e^{(0)}{}_\mu \frac{{{\cal D}}e_{(0)}{}^\mu}{{{\cal D}}\lambda} = 0 \,, \qquad
 n_\mu \frac{{{\cal D}}n^\mu}{{{\cal D}}\lambda} = 0 \,.$$ This means that the derivative of $e_{(0)}{}^\mu$ is orthogonal to $e^{(0)}{}_\mu$ and the derivative of $n^\mu$ is orthogonal to $n_\mu$. We thus find that $$S_\mu{}^\alpha S_\nu{}^\beta \epsilon_{\gamma \alpha \beta \delta}
 \frac{{{\cal D}}(e_{(0)}{}^\gamma n^\delta)}{{{\cal D}}\lambda}
 = S_\mu{}^\alpha S_\nu{}^\beta \epsilon_{\gamma \alpha \beta \delta}
 \left( \frac{{{\cal D}}e_{(0)}{}^\gamma}{{{\cal D}}\lambda} n^\delta
 + e_{(0)}{}^{\gamma} \frac{{{\cal D}}n^\delta}{{{\cal D}}\lambda} \right) = 0 \,.$$ Finally, we obtain that the Boltzmann equation for the tensor-valued distribution functions can be split into the desired form as $$\begin{gathered}
 S_\mu{}^\alpha S_\nu{}^\beta \frac{{{\cal D}}f_{\alpha \beta}}{{{\cal D}}\lambda}
 = \frac{1}{2} \frac{{{\cal D}}I}{{{\cal D}}\lambda} S_{\mu \nu}
 + S_\mu{}^\alpha S_\nu{}^\beta \frac{{{\cal D}}P_{\alpha \beta}}{{{\cal D}}\lambda}
 + \frac{{{\rm i}}}{2} \frac{{{\cal D}}V}{{{\cal D}}\lambda}
 \epsilon_{\alpha \mu \nu \beta} e_{(0)}{}^\alpha n^\beta = C_{\mu \nu} \,.\end{gathered}$$

Expression of the Boltzmann equation for polarisation {#sapp:Boltzeq_pol}
-----------------------------------------------------

First of all, the Boltzmann equation for the circular polarisation is the same as that for the intensity, but with $\bar V=0$ and with a vanishing collision term since it is not generated by the Compton scattering. We shall not study further the equation dictating the evolution of $V$ since it should remain null at all time unless generated by other types of collisions.

As for the linear polarisation, let us write down the basic equations for the tetrad components as commonly done at first order. By moving to the tetrad components, one can rewrite the covariant derivative in the Liouville operator as $$\begin{aligned}
 e_{(a)}{}^\mu e_{(b)}{}^\nu L [{\bf P}\,]_{\mu \nu}
 &= S_{(a)}{}^{(c)} S_{(b)}{}^{(d)} \left( P_{(c) (d) | (e)}
 e^{(e)}{}_\mu \frac{{{\rm d}}x^\mu}{{{\rm d}}\lambda}
 + \frac{{{\partial}}P_{(c) (d)}}{{{\partial}}\ln q} \frac{{{\rm d}}\ln q}{{{\rm d}}\lambda}
 + D_{(i)} P_{(c) (d)} \frac{{{\rm d}}n^{(i)}}{{{\rm d}}\lambda} \right) \,,\end{aligned}$$ where $$\begin{aligned}
 f_{(a) (b) | (c)}
 &\equiv e_{(c)}{}^\mu {{\partial}}_\mu f_{(a) (b)}
 - w^{(d)}{}_{(a) (c)} f_{(d) (b)} - w^{(d)}{}_{(b) (c)} f_{(a) (d)} \,,\end{aligned}$$ and $w^{(a)}{}_{(b) (c)}$ is the Ricci rotation coefficient defined by $w^{(a)}{}_{(b) (c)} \equiv e^{(a)}{}_\mu \nabla_{(c)}
 e_{(b)}{}^\mu$. Here we also notice that only the spatial component has non-vanishing term because the projection of $S_\mu{}^\nu$ onto $e^\mu{}_{(0)}$ vanishes by construction. After the decomposition of the Boltzmann equation for polarisation, one obtains the following equation for the temperature part as in Eq. (\[eq:Boltz\_pol\_dis\]) up to the second order $$\begin{aligned}
 {{\cal L}}^P_{(i) (j)}
 &= {{\cal P}}_{(i) (j)}{}' + {{\cal P}}_{(i) (j), k} n^{(k)}
 + \frac{a^2}{q}\left[\left. \frac{{{\rm d}}\eta}{{{\rm d}}\lambda} \right|^{(1)} {{\cal P}}_{(i) (j)}{}'
 + \left. \frac{{{\rm d}}x^k}{{{\rm d}}\lambda} \right|^{(1)} {{\cal P}}_{(i) (j), k}
 + \left. \frac{{{\rm d}}n^{(k)}}{{{\rm d}}\lambda} \right|^{(1)} D_{(k)} {{\cal P}}_{(i) (j)}\right]
 \notag\\
 & \qquad
 - a\Bigl( w^{(k)}{}_{(l) (0)} + w^{(k)}{}_{(l) (m)} n^{(m)} \Bigr)^{(1)}
 \Bigl( S_{(i)}{}^{(l)} {{\cal P}}_{(k) (j)}
 + S_{(j)}{}^{(l)} {{\cal P}}_{(k) (i)} \Bigr)\,.\end{aligned}$$ As for the spectral distortion, the equation is given by $$\begin{aligned}
 {{\cal L}}^Y_{(i) (j)}
 &= Y_{(i) (j)}{}' + Y_{(i) (j), k} n^{(k)} \,.\end{aligned}$$ The complete expression of the collision term for the linear polarisation is given by $$\begin{aligned}
 {{\cal C}}_{(i) (j)}^P
 &= a \, \bar{n}_e \sigma_T \left(
 - {{\cal P}}_{(i) (j)} - \frac{3}{4} {{\cal T}}_{(i) (j)}{}^{(k) (l)}
 \Bigl[ \langle \Theta m_{(k) (l)} \rangle
 - 2 \langle {{\cal P}}_{(k) (l)} \rangle \Bigr]
 + v^{(k)} n_{(k)} {{\cal P}}_{(i) (j)} \right) \notag\\
 &\qquad
  + a \, \bar{n}_e \sigma_T \biggl(
 {{\cal T}}_{(i) (j)}{}^{(k) (l)} \Bigl[ {{\cal Q}}_{(k) (l)}^T
 - \Theta {{\cal Q}}_{(k) (l)}^{T ~~~ (1)} \Bigr] 
 - 3 {{\cal C}}^\Theta {{\cal P}}_{(i) (j)} + \delta_e {{\cal C}}^P_{(i) (j)} \biggr) \,, \\
 C^Y_{(i) (j)}
 &= a \, \bar{n}_e \sigma_T \left( - Y_{(i) (j)}
 - \frac{3}{4} {{\cal T}}_{(i) (j)}{}^{(k) (l)}
 \Bigl[ \langle y m_{(k) (l)} \rangle
 - 2 \langle Y_{(k) (l)} \rangle \Bigr]
 + {{\cal T}}_{(i) (j)}{}^{(k) (l)} {{\cal Q}}_{(k) (l)}^Y 
 - {{\cal C}}^\Theta {{\cal P}}_{(i) (j)} \right) \,,\end{aligned}$$ where ${{\cal T}}_{(i) (j)}{}^{(k) (l)}$ is a traceless projection operator with respect to $S_{(i) (j)}$ $$\begin{aligned}
 {{\cal T}}_{(i) (j)}{}^{(k) (l)}
 \equiv S_{(i)}{}^{(k)} S_{(j)}{}^{(l)} - \frac{1}{2} S^{(k) (l)} S_{(i) (j)} \,.\end{aligned}$$

Finally before closing this section let us explicitly write down the equation for the temperature part for the sake of completeness. It is of the form $$\begin{aligned}
 {{\cal L}}^P_{(i) (j)} = {{\cal C}}_{(i) (j)}^P \,,
 \end{aligned}$$ and using that at first order $$w^{(k)}{}_{(l) (0)} = \frac{1}{a}\beta^{[k}{}_{, l]} \,,\qquad
 w^{(k)}{}_{(l) (m)} =\frac{2}{a}  h_m{}^{[k}{}_{, l]} \,,$$ the explicit form of ${{\cal L}}^P_{(i) (j)}$ is $$\begin{aligned}
 {{\cal L}}_{(i) (j)}^P \equiv
 &{{\cal P}}_{(i) (j)}{}' + {{\cal P}}_{(i) (j), k} n^{(k)} \,,
 - \alpha {{\cal P}}_{(i) (j)}{}' - (\beta^k + h^k{}_l n^{(l)}) {{\cal P}}_{(i) (j), k} 
 \notag\\
 & \qquad
 - \Bigl( \alpha_{, k} - \beta_{l, k} n^{(l)} + h_{k l}{}' n^{(l)}
 + 2 h_{m [k, l]} n^{(l)} n^{(m)} \Bigr) D^{(k)} {{\cal P}}_{(i) (j)}
 \notag\\
 & \qquad \qquad
 - \Bigl( \beta^{[k}{}_{, l]} + 2 h_m{}^{[k}{}_{, l]} n^{(m)} \Bigr)
 \Bigl( S_{(i)}{}^{(l)} {{\cal P}}_{(k) (j)} + S_{(j)}{}^{(l)} {{\cal P}}_{(k) (i)} \Bigr) \,.\end{aligned}$$

Technical details about linear polarisation
===========================================

Gauge transformation of a tensor-valued distribution function {#AppGTtensor}
-------------------------------------------------------------

In this section, we investigate the transformation property of the distribution matrix up to the second order. The transformation properties of $I$ and $V$, which are scalar distribution functions, have already been investigated in the main text. In this section we focus on the case of the linear polarisation and set $I=V=0$ such that $f_{\mu\nu}=P_{\mu\nu}$. Since the polarisation matrix is at least first order, this means that we need only to keep terms which are first order in the coordinates transformation $(T,L^i)$.

First we must understand the transformation properties of the screen projector. As for the tetrads, it is not a purely geometric quantity and it is not invariant under a change of coordinates. Indeed, it is defined with respect to the time-like tetrad which depends on the choice of coordinates. We thus have the two related screen projectors $${\bm S}({\bm p})={\bm g}+{\bm e}^{(0)} \otimes {\bm e}^{(0)}-{\bm
  n}\otimes{\bm n}\,,\quad \tilde{\bm S}({\bm p})={\bm g}+\tilde {\bm
  e}^{(0)} \otimes \tilde {\bm e}^{(0)}-\tilde{\bm
  n}\otimes\tilde{\bm n}\,,\quad\text{with}\quad {\bm n}\equiv
\frac{{\bm p}}{{p}^\mu{e}^{(0)}_\mu}-{\bm e}^{(0)}\,,\quad \tilde{\bm n}\equiv
\frac{{\bm p}}{{p}^\mu \tilde {e}^{(0)}_\mu}-\tilde{\bm e}^{(0)}\,.$$ If we use the transformation rule at the same point  we find that at first order in the transformation the relation between these projectors is given by  [@Tsagas:2007yx] $$\label{RelSStilde}
\tilde S_{\mu \nu}\circ \tilde c{(\tilde{\bm x},q^{\widetilde{(a)}})}=
S_{\mu\nu}({\bm x},q^{(a)})+ 2 \left(e^{(0)}_{(\mu}+n_{(\mu}\right) S_{\nu)
  \alpha}({\bm x},q^{(a)}) V^\alpha\,,\quad \text{with}\quad {\bm V}
\equiv -\Lambda^{(i)}_{\,\,(0)} {\bm e}_{(i)}=-{{\partial}}^i T {\bm e}_i\,.$$ Note that if we use tetrad coordinates for the argument of the screen projectors, then $${\bm S}(q,n^{(i)})={\bm g}+{\bm e}^{(0)} \otimes {\bm e}^{(0)}-n_{(i)}
n_{(j)}{\bm e}^{(i)}\otimes {\bm e}^{(j)}\,,\qquad \tilde{\bm
  S}\circ \tilde c{(q,n^{(i)})}={\bm g}+\tilde{\bm e}^{(0)} \otimes \tilde{\bm
  e}^{(0)}-n_{(i)} n_{(j)}\tilde {\bm e}^{(i)} \otimes  \tilde {\bm e}^{(j)}\,.$$ This means that when the coordinates of the projectors are expressed in the tetrad basis associated with the corresponding coordinates system, we obtain $$\label{SijisSij}
S_{(i)(j)}(q,n^{(k)})=\delta_{ij}-n_{(i)}n_{(j)} = \tilde S_{\widetilde{(i)}\widetilde{(j)}}\circ \tilde c{(q,n^{(k)})}\,.$$ Thus in any gauge, the expression of the related screen projector in the tetrad basis is the same by construction, and $S_{(i)(j)}$ depends actually only on $n^{(k)}$. However it must remain clear that these two tensors are geometrically different since they are associated with different coordinates systems and indeed their relation is given at first order by .

The distribution tensor is also not geometrically invariant since for every observer used to define the screen projector, we must consider a different distribution tensor. However all the possible distribution matrices are related through projections and we find that the distribution tensors defined by the tetrad $\tilde{\bm e}_{(0)}$ and ${\bm e}_{(0)}$ are related at the same point of the tangent bundle by [@Tsagas:2007yx] $$f^{\tilde{\bm e}_{(0)}}_{\mu\nu}\circ
 \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 = \tilde S_\mu^\alpha \circ \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 \tilde S_\nu^\beta \circ \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 f^{{\bm e}_{(0)}}_{\alpha \beta}({\bm x},q^{(\imath)})\,.$$ This means that the only requirement is to project the distribution function so that it is projected with respect to the new observer and the new direction. Combining this transformation rule with  we obtain at first order the transformation rule as $$f^{\tilde{\bm e}_{(0)}}_{\mu\nu}\circ
 \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 = f^{{\bm e}_{(0)}}_{\mu \nu}({\bm x},q^{(\imath)})
 + 2 \left(e^{(0)}_{(\mu}+n_{(\mu}\right)
 f^{{\bm e}_{(0)}}_{\nu)\alpha}({\bm x},q^{(\imath)})V^\alpha\,.$$ If we project this expression onto the tetrad components, and noting $\tilde f_{(i)(j)}\equiv f^{\tilde{\bm e}_{(0)}}_{(i)(j)}$ and $f_{(i)(j)}\equiv f^{{\bm e}_{(0)}}_{(i)(j)}$, we obtain $$\label{Rulefijstrange}
 \tilde f_{\widetilde{(i)}\widetilde{(j)}}\circ
 \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 = f_{(i)(j)}({\bm x},q^{(\imath)})
 - f_{(i)(k)}({\bm x},q^{(\imath)})L_{\,\,\,,j]}^{[k}
 -f_{(k)(j)}({\bm x},q^{(\imath)})L_{\,\,\,,i]}^{[k}
 - 2 n_{((i)} f_{(j))(k)}({\bm x},q^{(\imath)}){{\partial}}^k T\,,$$ find the transformation rule under a gauge transformation, we need to expand the left hand side around $({\bm x},q^{(\imath)})$. At first order we get $$\begin{aligned}
\tilde f_{\widetilde{(i)}\widetilde{(j)}}\circ
 \tilde c{(\tilde {\bm x},q^{\widetilde{(\imath)}})}
 &\simeq\left(1+\xi^\mu \frac{{{\partial}}}{{{\partial}}x^\mu}
 +\delta q^{(i)}\frac{{{\partial}}}{{{\partial}}q^{(i)}} \right)
 \tilde f_{\widetilde{(i)}\widetilde{(j)}}\circ
\tilde c{({\bm x},q^{{(\imath)}})}+\cdots \nonumber\\
&\simeq \tilde f_{\widetilde{(i)}\widetilde{(j)}}\circ
 \tilde c{({\bm x},q^{{(\imath)}})}
 + \left(\xi^\mu \frac{{{\partial}}}{{{\partial}}x^\mu}+\delta q^{(i)}\frac{{{\partial}}}{{{\partial}}q^{(i)}}
 \right) f_{{(i)}{(j)}}+\cdots \,.\end{aligned}$$ Using the expansion (\[Expandderpi\]) we then obtain the gauge transformation rule for the tensor valued distribution function in tetrad coordinates. Noting $\widetilde{f_{(i)(j)}} \equiv \tilde
f_{\widetilde{(i)}\widetilde{(j)}}\circ \tilde c$ for simplicity, this reads $$\label{Gaugeruletensor}
\widetilde{f_{{(i)}{(j)}}} = f_{(i)(j)}-\left(\xi^\mu \frac{{{\partial}}}{{{\partial}}x^\mu}+\delta \ln q {{\frac{\partial}{\partial \ln q}}}+\delta n^{(i)}D_{(i)}  \right)f_{(i)(j)}-f_{(i)(k)}
L_{\,\,\,,l]}^{[k}S^{(l)}_{(j)} -f_{(k)(j)} L_{\,\,\,,l]}^{[k} S^{(l)}_{(i)}\,,$$ where it is implied that all quantities are evaluated either at ${\bm x}$ or at $({\bm x},q,n^{(i)})$.

For completeness we report the explicit form of the gauge transformation for ${{\cal P}}_{(i) (j)}$ which is obtained from the above transformation rule and the spectral decomposition (\[defyP\]) $$\begin{aligned}
\label{GT:calP}
 \tilde{{{\cal P}}}_{(i) (j)}
 &= {{\cal P}}_{(i) (j)}
 - L^{[k}{}_{, l]} \Bigl( S^{(l)}{}_{(i)} {{\cal P}}_{(k) (j)}
 + S^{(l)}{}_{(j)} {{\cal P}}_{(k) (i)} \Bigr) - \xi^\mu {{\partial}}_\mu {{\cal P}}_{(i) (j)}
 - \delta n^{(k)} D_{(k)} {{\cal P}}_{(i) (j)} \,.
 \end{aligned}$$ It is also found, as expected, that the spectral distortion $Y_{(i)(j)}$ part is gauge invariant since it vanishes on the background and at first order.

Gauge transformation for Liouville and Collision terms {#AppGTLC}
------------------------------------------------------

We deduce from the transformation rule (\[GT:calP\]) and the spectral decomposition (\[spdec:calP\_lhs\]) and (\[spdec:calP\_rhs\]) that the spectral distortion part, ${{\cal L}}^Y_{(i)(j)} $ and ${{\cal C}}^Y_{(i)(j)}$, must be gauge invariant. Concerning the temperature part, they should transform as (noting $\widetilde{{{\cal L}}^P_{(i)(j)}}
 \equiv {{\cal L}}^P_{\tilde{(i)}\tilde{(j)}} \circ \tilde c$ and $\widetilde{{{\cal C}}^P_{(i)(j)}}
 \equiv {{\cal C}}^P_{\tilde{(i)}\tilde{(j)}} \circ \tilde c$ ) $$\begin{aligned}
 \widetilde{{{\cal L}}^P_{(i)(j)}}
 &= {{\cal L}}^P_{(i)(j)} - L^{[k}{}_{, l]} \Bigl( S^{(l)}{}_{(i)} {{\cal L}}^P_{(k) (j)}
 + S^{(l)}{}_{(j)} {{\cal L}}^P_{(k) (i)} \Bigr) - \xi^\mu {{\partial}}_\mu {{\cal L}}^P_{(i)(j)}
 - \delta n^{(k)} D_{(k)} {{\cal L}}^P_{(i) (j)} \,, \\
 \widetilde{{{\cal C}}^P_{(i)(j)}}
 &= {{\cal C}}^P_{(i)(j)} - L^{[k}{}_{, l]} \Bigl( S^{(l)}{}_{(i)} {{\cal C}}^P_{(k) (j)}
 + S^{(l)}{}_{(j)} {{\cal C}}^P_{(k) (i)} \Bigr) - \xi^\mu {{\partial}}_\mu {{\cal C}}^P_{(i)(j)}
 - \delta n^{(k)} D_{(k)} {{\cal C}}^P_{(i) (j)} \,.
 \end{aligned}$$ Using the transformation rules derived in this paper, we checked that this is indeed the case when using the detailed form of the Liouville and collision operators.

Extraction of temperature and spectral distortion {#AppExtraction}
=================================================

The functions $y$ and $\Theta$ can be extracted thanks to the integrals of the type $${{\cal M}}_n[f]\equiv \frac{\int f q^{2+n} {{\rm d}}q}{(3+n)
 \int \bar I(q) q^{2+n} {{\rm d}}q} \,,$$ just by applying them order by order to $I(q)$, using that ${{\cal M}}_0[{{\cal D}}_q^2 \bar I]=0$. We then obtain

$$\begin{aligned}
\Theta^{(1)} &= {{\cal M}}_1 [I^{(1)}] = {{\cal M}}_0[I^{(1)}] \,,
 \label{ExtractTheta1} \\
 \frac{1}{2}\Theta^{(2)} &= \frac{1}{2} {{\cal M}}_0[I^{(2)}]-\Theta^{(1)2} \,, \\
 \frac{1}{2}y^{(2)} &= \frac{1}{2} \left( {{\cal M}}_1[I^{(2)}] - {{\cal M}}_0 [I^{(2)}]
 \right) - \frac{1}{2}\Theta^{(1)}{}^2 \,.\end{aligned}$$

Similarly to what can be done for the intensity part, the spectral components of polarisation can be extracted thanks to

$$\begin{aligned}
 {{\cal P}}^{(1)}_{\mu\nu}
 &= {{\cal M}}_1 [{\bf P}_{\mu\nu}^{(1)}] = {{\cal M}}_0 [P_{\mu\nu}^{(1)}] \,, \\
 \frac{1}{2} {{\cal P}}^{(2)}_{\mu\nu}
 &= \frac{1}{2} {{\cal M}}_0 [P_{\mu\nu}^{(2)}]
 - 3 \Theta^{(1)} {{\cal P}}^{(1)}_{\mu\nu} \,,\\
 \frac{1}{2}Y^{(2)}_{\mu\nu}
 &= \frac{1}{2} \Bigl( {{\cal M}}_1[P_{\mu\nu}^{(2)}]-{{\cal M}}_0[P_{\mu\nu}^{(2)}] \Bigr)
 - \Theta {{\cal P}}^{(1)}_{\mu\nu} \,.\end{aligned}$$

[10]{}

WMAP, E. Komatsu [*et al.*]{}, Astrophys. J. Suppl. [**192**]{}, 18 (2011), arXiv:1001.4538. Planck Collaboration, P. Ade [*et al.*]{}, (2013), arXiv:1303.5084. E. Komatsu and D. N. Spergel, Phys. Rev. [**D63**]{}, 063002 (2001), arXiv:astro-ph/0005036. J. M. Maldacena, JHEP [**05**]{}, 013 (2003), arXiv:astro-ph/0210603. D. Langlois and F. Vernizzi, Phys.Rev.Lett. [**95**]{}, 091303 (2005), arXiv:astro-ph/0503416.

D. Langlois and F. Vernizzi, Phys.Rev. [**D72**]{}, 103501 (2005), arXiv:astro-ph/0509078.

K. Enqvist, J. Hogdahl, S. Nurmi, and F. Vernizzi, Phys.Rev. [**D75**]{}, 023515 (2007), arXiv:gr-qc/0611020.

C. Pitrou and J.-P. Uzan, Phys.Rev. [**D75**]{}, 087302 (2007), arXiv:gr-qc/0701121.

G. L. Comer, N. Deruelle, D. Langlois, and J. Parry, Phys. Rev. [**D49**]{}, 2759 (1994). E. W. Kolb, S. Matarrese, A. Notari, and A. Riotto, Mod. Phys. Lett. [**A20**]{}, 2705 (2005), arXiv:astro-ph/0410541. D. H. Lyth, K. A. Malik, and M. Sasaki, JCAP [**0505**]{}, 004 (2005), arXiv:astro-ph/0411220.

N. Bartolo, S. Matarrese, and A. Riotto, JCAP [**0606**]{}, 024 (2006), arXiv:astro-ph/0604416.

N. Bartolo, S. Matarrese, and A. Riotto, JCAP [**0701**]{}, 019 (2007), arXiv:astro-ph/0610110.

C. Pitrou, Class. Quant. Grav. [**26**]{}, 065006 (2009), arXiv:0809.3036. C. Pitrou, Gen.Rel.Grav. [**41**]{}, 2587 (2009), arXiv:0809.3245.

C. Pitrou, J.-P. Uzan, and F. Bernardeau, JCAP [**1007**]{}, 003 (2010), arXiv:1003.0481. M. Beneke and C. Fidler, Phys. Rev. [**D82**]{}, 063509 (2010), arXiv:1003.1834. R. Durrer, Fund.Cosmic Phys. [**15**]{}, 209 (1994), arXiv:astro-ph/9311041.

C. Pitrou, Class. Quant. Grav. [**24**]{}, 6127 (2007), arXiv:0706.4383. Z. Huang and F. Vernizzi, (2012), arXiv:1212.3573. S.-C. Su, E. A. Lim, and E. Shellard, (2012), arXiv:1212.6968. G. W. Pettinari, C. Fidler, R. Crittenden, K. Koyama, and D. Wands, (2013), arXiv:1302.0832. P. Creminelli, C. Pitrou, and F. Vernizzi, JCAP [**1111**]{}, 025 (2011), arXiv:1109.1822. Planck Collaboration, P. Ade [*et al.*]{}, (2013), arXiv:1303.5081. A. Stebbins, (2007), arXiv:astro-ph/0703541. C. Pitrou, F. Bernardeau, and J.-P. Uzan, JCAP [**1007**]{}, 019 (2010), arXiv:0912.3655. R. L. Arnowitt, S. Deser, and C. W. Misner, (1962), arXiv:gr-qc/0405109. C. G. Tsagas, A. Challinor, and R. Maartens, Phys. Rept. [**465**]{}, 61 (2008), arXiv:0705.4397. S. Dodelson and J. M. Jubas, Astrophys. J. [**439**]{}, 503 (1995), astro-ph/9308019. M. Bruni, S. Matarrese, S. Mollerach, and S. Sonego, Class. Quant. Grav. [**14**]{}, 2585 (1997), arXiv:gr-qc/9609040. D. Hanson, K. M. Smith, A. Challinor, and M. Liguori, Phys.Rev. [**D80**]{}, 083004 (2009), arXiv:0905.4732. A. Lewis, A. Challinor, and D. Hanson, JCAP [**1103**]{}, 018 (2011), arXiv:1101.2234. E. Gourgoulhon, (2007), arXiv:gr-qc/0703035.

[^1]: We assume for the simplicity of the argument that one system of coordinates is enough to cover the entire manifold.

[^2]: Again here for simplicity, we assume that such coordinates system covers the whole manifold.

[^3]: Under our conventions, $\tilde q \equiv a p^{\widetilde{(0)}}$.
