=15.5pt empty

2028Second-order gauge-invariant formalism for the cosmological observables:
Complete verification of their gauge-invariance

1330Matteo Magi$^a$ and Jaiyul Yoo $^{a,b}$

$^a$ Center for Theoretical Astrophysics and Cosmology,
Institute for Computational Science

University of Zürich, Winterthurerstrasse 190,
CH-8057, Zürich, Switzerland

$^b$Physics Institute, University of Zürich,
Winterthurerstrasse 190, CH-8057, Zürich, Switzerland

August 27, 2024

<EMAIL_ADDRESS><EMAIL_ADDRESS>
Accounting for all the relativistic effects, we have developed the fully nonlinear gauge-invariant formalism for describing the cosmological observables and presented the second-order perturbative expressions associated with light propagation and observations without choosing a gauge condition. For the first time, we have performed a complete verification of the validity of our second-order expressions by comparing their gauge-transformation properties from two independent methods: one directly obtained from their expressions in terms of metric perturbations and the other expected from their nonlinear relations. The expressions for the cosmological observables such as galaxy clustering and the luminosity distance are invariant under diffeomorphism and gauge-invariant at the observed position. We compare our results to the previous work and discuss the differences in the perturbative expressions. Our second-order gauge-invariant formalism constitutes a major step forward
in the era of precision cosmology and its applications in the future will play a crucial role for going beyond the power spectrum and probing the early universe.

§ INTRODUCTION

Thanks to tremendous technological and observational advances in recent years,
we have witnessed an era of precision cosmology, and an unprecedented amount of
data in large-scale surveys allows us to measure important statistics and extract accurate information
about our universe. Furthermore, this trend will accelerate in the coming years,
making the near future of cosmology particularly promising. The upcoming large-scale
surveys on the ground include the Dark Energy Spectroscopic Instrument (DESI) <cit.>, the Vera
C. Rubin Observatory (formerly LSST) <cit.>, the Square Kilometre Array (SKA) <cit.>, and the next generation CMB experiment
CMB-S4 <cit.>, and several space-based missions will be operational such as
Euclid <cit.>, the Nancy Grace Roman Space Telescope (formerly WFIRST) <cit.> and the Laser Interferometer Space Antenna (LISA) <cit.>. These surveys will dramatically improve
our understanding of cosmology and fundamental physics, thanks to an unprecedented amount of high precision data. Cosmological perturbations on large scales originate
from random fields sourced by quantum fluctuations in the early universe. These
fields assume random values at each spatial point according to their corresponding
probability measure and cosmological information is encoded in the $N$-point
correlation functions. In order to take full advantage of the new opportunities
in the forthcoming era of precision cosmology, the theoretical predictions for
the $N$-point correlation functions should be made at least as accurate as
the level of precision set by the upcoming surveys.

The standard inflation model predicts a nearly Gaussian distribution in the early
universe <cit.>. While the non-Gaussianity we observe in large-scale structure largely
reflects the nonlinear evolution of gravity, a small amount of the primordial
non-Gaussianity might be hidden in our measurements (see, e.g., <cit.>).
To describe the deviation from Gaussianity, it is necessary to go beyond the
linear order in perturbation theory to capture statistics higher than the power
spectrum. The simplest measure of non-Gaussianity is the three-point correlation
function, which corresponds to the bispectrum in Fourier space. The most simple
type of non-Gaussianity is the local one (see, e.g., <cit.> for a review),
which peaks in the squeezed limit in Fourier space. Such a squeezed bispectrum is
particularly important because it breaks the degeneracy in galaxy bias <cit.> and
satisfies a particular consistency relation to the power spectrum <cit.>.
Since the standard inflationary model predicts a vanishingly small non-Gaussianity,
a detection of such primordial non-Gaussianity could potentially rule out a vast
class of inflationary models <cit.>. However, this exciting possibility can
only be realized, if we control the systematic errors in our theoretical
predictions below the level of signals we are interested in.

Cosmological perturbation theory has been extensively studied
in literature (see, e.g., <cit.>), and it solves the Einstein equation up
to the desired order in perturbations, after a gauge choice is often made
(see <cit.> for the calculations in a general metric representation).
In cosmological observations, there exists another step that relates the physical
solutions to the observations. Cosmological information is mostly carried by
photons (or gravitational waves), and they propagate along the past light cone to reach us in our rest frame.
For example, galaxy clustering is dominated by the matter density fluctuation,
which is a solution by the Einstein equation, but is not directly what we measure.
The light propagation and observations in the observer rest frame give rise
to extra contributions to galaxy clustering such as the redshift-space distortion,
gravitational lensing effects, the gravitational redshift, and other relativistic
effects (<cit.>, see <cit.> for a complete
treatment). A full relativistic description of galaxy clustering at the linear
order was developed <cit.>. For the same reason, the
relativistic effects are present in other cosmological probes such as the
luminosity distance <cit.>.
Given the relativistic nature of weak gravitational lensing and cosmic microwave
background (CMB) anisotropies, their theoretical descriptions were largely
complete from the inception for weak lensing <cit.> and for CMB <cit.> (see, e.g., <cit.> for a review).
However, extra relativistic contributions were found even at the linear order
in perturbations for weak gravitational lensing <cit.> and for CMB

Given these recent developments, it is therefore natural to develop the
second-order relativistic description of the cosmological observables. With the high
precision measurements of CMB anisotropies, the theoretical descriptions of CMB beyond
linear order have been developed in early days <cit.>, playing
a critical role in constraining the primordial non-Gaussianity and
deriving the lensing potential (see, e.g., <cit.>),
though only the linear-order deviation in the light propagation
from a straight path is needed for describing the second-order CMB anisotropies. The second-order description of light
propagation was employed in weak lensing formalism <cit.> and also
the luminosity distance <cit.>. The second-order relativistic description
of galaxy clustering and its applications to the bispectrum were recently
developed in <cit.>,
in which the second-order light propagation and
observations in the rest frame were carefully accounted for.

Despite all the efforts and developments in recent years, there exists little consensus in the second-order
calculations. Although the extension of the linear-order formalism to the
second-order one is conceptually straightforward, it involves remarkably complex expressions and subtleties, making it extremely challenging to track down the
differences and to make a meaningful comparison. A further complication arises
from gauge issues. The theory of general relativity is manifestly covariant;
however, splitting variables into a background part and perturbations breaks this
covariance by introducing a dependence on the choice of coordinates. It is common
practice to fix a gauge to work with more manageable expressions. However, by doing
so, one loses the ability of checking the gauge-invariance of the final results
for the cosmological observables, i.e., once gauge freedom is fully fixed,
all the resulting expressions are automatically gauge-invariant, not necessarily because
the expressions of the cosmological observables are correct, but simply because
gauge freedom was removed by hand. In this paper, we work with a general metric
representation without fixing a gauge. Only then, it is possible to exploit the geometric
nature of the defined quantities and check that they gauge-transform correctly
under an arbitrary change of coordinates. In the past this consistency test has
been performed only for the linear-order formalism.

Drawing on the previous work, here we develop a nonlinear theoretical framework
to describe the cosmological observables, present their second-order
perturbative expressions, and verify the gauge-transformation properties
of the second-order expressions. Adopting a geometric approach to describing
the light propagation and observations and using the tetrad formalism <cit.>,
we describe the observations in our rest frame and physical quantities in the
source rest frame, and their connection through light propagation was computed
in a nonlinear way. Thereby, we can clarify the separation
between the background quantities and the nonlinear perturbative quantities
in the theoretical descriptions of the cosmological observables.
By carefully expanding all the nonlinear expressions along the background light path and around the coordinate-independent reference positions of the observer and the
source parameterized by the observables, we derive the second-order expressions
of the cosmological observables that contain linear terms in the second-order
perturbation variables and quadratic terms of the linear-order perturbation
variables. The perturbation contributions at the observer position
are typically ignored in previous work, but their presence is crucial to ensure the
correctness of the expressions <cit.>. In particular, they couple to other perturbation
variables at second order in perturbations, acquiring the positional dependence.
Finally, we compute the gauge transformations of all the second-order expressions,
which are derived directly from the gauge transformation of their metric components
in the expressions. Utilizing the nonlinear expressions of those, we compute the
nonlinear (and second-order) expectations for the gauge transformations and verify
the consistency in all of our expressions. Again, we stress that these sanity checks
are a necessary step to establish the complete second-order formalism in the era of
precision cosmology and they constitute one of the main results in the present work.

The paper is structured as follows. Section $\ref{sec:nonlinear}$ introduces the nonlinear geometric formalism needed to treat all relativistic effects associated with light propagation in observation from a distant source to the observer. In Section <ref>, we present the second-order perturbative expressions of the cosmological observables in Section $\ref{sec:nonlinear}$. The main expressions we derive are the observed redshift fluctuations and the physical volume and area occupied by the source. In Section $\ref{sec:gauge}$ we explicitly verify that the perturbative expressions of Section $\ref{sec:second}$ gauge-transform according to the expectation derived from the definitions in Section $\ref{sec:nonlinear}$. Section $\ref{sec:compa}$ compares our expressions with previous work in literature, and in Section $\ref{sec:concl}$, we summarize and discuss our results. Moreover, Appendix $\ref{AppendixA}$ is reserved to further technical calculations related to the observer’s position, while Appendix $\ref{AppendixB}$ displays the second-order expressions in their complete form.

§ NONLINEAR EXPRESSIONS FOR COSMOLOGICAL OBSERVABLES

We present nonlinear expressions for cosmological observables obtained by the light propagation from the source to the observer position. In particular, we provide expressions for the observer and source positions, the observed redshift, and the area and volume occupied by the source; these will then be written in terms of nonlinear perturbative and exact background parts. The nonlinear perturbative part can later be expanded order-by-order. This section aims to lay the groundwork necessary to perform rigorous second-order perturbative calculations, which are the main focus of this work and are presented in Section $\ref{sec:second}$. Here we follow the general approach in <cit.> with nonlinear generalization in coordinates (see also <cit.>).

§.§ Geometric setup

We begin by introducing fundamental geometric quantities and stating our notation conventions. We denote with $\bar\mM$ the Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime with metric $\bar g$, it describes a homogeneous isotropic universe and is an excellent approximation to our universe on cosmological scales. We denote with $\mM$ the spacetime describing the real inhomogeneous universe with metric $g$, which deviates perturbatively from the FLRW background metric. The coordinate system we use is that of FLRW, in which the temporal coordinate is the cosmic time $t$, or more commonly the conformal time $\eta$, while the spatial coordinates are the comoving coordinates $x^\al$. We recall that with this choice, the fixed-time sections of $\bar\mM$ are maximally symmetric, and the comoving coordinates stretch according to the Hubble flow.

We adopt the following convention for the components of the background metric:
g̅_ηη(η)=-a^2(η)  ,g̅_(η)= a^2(η) _ ,
where $\de_{\al\be}$[We use Greek letters $\al, \be,\cdots$ as indices that can take values from one to three, while $\mu, \nu,\cdots$ are indices that can take values from zero to three.] is the Kronecker delta and $a(\eta)$ is the scale factor, containing the dynamics of the background universe. Note that the metric signature is mostly-plus, and we are restricting to the case of a flat universe with Cartesian comoving coordinates.
Our convention for the components of the metric of the real inhomogeneous universe is

g_ηη(η,x)   =   -a^2(η) ( 1 + 2Å(η,x) ) ,g_η(η,x)=-a^2(η) _α(η,x) ,
g_(η,x)   =   a^2(η)( _αβ+ 2 _αβ(η,x)) .
where spacetime fields $\AA,\, \BB_\al,\, \CC_{\al\be}$ are nonlinear perturbations to the background metric, and they are based on the flat 3-metric $\de_{\al\be}$. These perturbations represent the inhomogeneities of the real universe and if they vanish we recover $\bar g_{\mu\nu}$ the components of the background metric.

Geometrically, we describe the trajectories of the source and the observer with the curves $\mathcal{S}$ and $\OO$, respectively. We introduce a curve $\gamma$ to describe the light path intersecting the other two curves at $S$ and $O$. These points denote the position of the source at light emission time and the observer's position today. We affinely parameterize the curves $\mathcal{S}$ and $\OO$ with $\tau'$ and $\tau$, proper times of the source and the observer, such that $\mathcal{S}(\tau')$ is any point on the source world line and $\mathcal{O}(\tau)$ is any point on the observer world line. To parameterize the curve $\gamma$ we introduce the affine parameter $\Lambda$, which for the moment is arbitrary.
Moreover, we assume that $\mathcal{S}$ and $\OO$ are time-like curves, while $\gamma$ is a null geodesic.
We now consider a neighborhood $U$ of the image of $\gamma$ and assume that there exists $(U,\psi)$ chart of $\mM$. So, we can express the curves $\OO$, $\mathcal S$ (restricted to $U$), and $\gamma$ in coordinates as

x^μ(τ):=(ψ^μ∘)(τ) ,    x^μ(τ'):=(ψ^μ∘𝒮)(τ') ,    x^μ(Λ):=(ψ^μ∘)(Λ) .
Figure $\ref{fig:illustration}$ summarizes what has been introduced so far.

Sketch of the geometric configuration considered in the real inhomogeneous universe. The source follows the time-like world line $\mathcal{S}(\tau')$, and the observer follows the time-like world line $\OO(\tau)$, with $\tau'$ and $\tau$ proper times. Today the observer is at point $O$ and sees the source at light emission point $S$ via the null light path $\ga$. $U$ is the neighborhood in which we give coordinates to the curves with the chart $(U,\psi)$.

It is convenient to repeat the whole construction at the background level. For this purpose, we consider what we just introduced in the real universe $\mM$ and set the perturbations of the metric $g$ to zero. Hence, the real position of the observer today reduces to $\bar O$, and the real position of the source at emission time reduces to $\bar S$. Since the geometry of $\bar\mM$ is different from that of $\mM$, the curves introduced will also differ in the two manifolds; however, it is still possible to adopt the same parameterizations since these procedures do not depend on the target spaces.
By assigning coordinates to these background curves in the same way as before, we define the deviations of physical trajectories from the background trajectories as follows

x^μ(τ):=x^μ(τ)-x̅^μ(τ) ,
x^μ(τ'):=x^μ(τ')-x̅^μ(τ') ,
x^μ(Λ):=x^μ(Λ)-x̅^μ(Λ) ,
where the overbar denotes the coordinates of the curves defined in the homogeneous and isotropic spacetime $\bar\mM$.

We now narrow our attention to the light path $\ga(\Lambda)$. At each point of this path, through the chart $(U,\psi)$, we have a vector basis $\{\pa_\mu\}$, with which we can represent the components of any tensor. However, the components written in this basis have no physical meaning, as they transform each time we change coordinate system. The coordinate-induced vector basis is not suitable for describing what the observer is measuring. In fact, the measurements made by the observer occur in the privileged system of the laboratory or observatory, where the metric is Minkowski, and the outcome of such measurements does not depend on the coordinate choice of the FLRW manifold. One way to geometrically define the observer rest frame is to define a tetrad basis $\{e_a\}$, with Latin indices[
We use Latin letters to denote internal indices, while Greek letters are for spacetime indices. Hence, a quantity with Latin indices is not involved in FLRW coordinate transformation, while a quantity with Greek indices is. Our notation convention is: $\al,\be,\cdots$ and $i,j,\cdots$ run from one to three, while $\mu,\nu,\cdots$ and $a,b,\cdots$ run from zero to three.]
$a,b,\cdots$ running from zero to three. This frame is fixed in the tangent space of $\mM$, such that tensor components in this basis are independent of the coordinate system. Moreover, the tetrad basis is orthonormal by definition, and the components of the metric tensor in this basis are Minkowski $g(e_a,e_b)=\eta_{ab}$. We conclude that the tetrad basis has all the characteristics needed to describe the observer's measurements; hence we equip the point $O$ with such a frame. We also equip the point $S$ with a tetrad basis to define quantities in the source rest frame. Actually, we will need a tetrad basis at any point of $\ga(\Lambda)$ to describe fictitious observations made by fictitious observers along the light path.[To be more general one can imagine fictitious observers with tetrad vectors everywhere in the universe and describe the physical events as they could be measured by a fictitious observer in the rest frame.]
The physical specification of the tetrad basis will be presented in the following subsection because we have not yet introduced all the necessary ingredients.

Let us introduce our notation conventions. We indicate the coordinates of a path, either in the real or background universe, with the affine parameter as a subscript, e.g., $x^\mu_\Lambda$ and $\bar x^\mu_\Lambda$ stand respectively for $x^\mu(\Lambda)$ and $\bar x^\mu(\Lambda)$.
To indicate the coordinates of a specific point along a path, we display as a subscript the lowercase letter corresponding to the point, e.g., $x^\mu_o$ are the coordinates of the observer position $O$, parameterized for example by $\Lambda_o$. These notations also apply to deviations in a path, such that $\de x^\mu_\Lambda$ stands for $\de x^\mu(\Lambda)$ and $\de x^\mu_o$ are the coordinate perturbations at the observer position. For the coordinates of specific points along the background path, we will use roman subscripts to emphasize that these coordinates refer to the homogeneous and isotropic universe, e.g., $\bar x^\mu_{\rm o}$, but note that the point $\bar O$ is still parameterized by $\Lambda_o$. To denote the evaluation of a field, $k^\mu(x)$, at a specific point of a path, we use the same subscript of the coordinates of the point, e.g., $k^\mu_o$ stands for $k^\mu(x_o)$. If the evaluation is on the background path, we additionally put an overbar on the roman subscript to expressly refer to $\bar x^\mu$, e.g., $k^\mu_{\bobs}$ stands for $k^\mu(\bar x_{\rm o})$.
Finally, we use a prime to indicate partial derivatives with respect to conformal time and a dot to express partial derivatives with respect to cosmic time.

§.§ Observer four-velocity and space-like tetrad vectors

We now introduce a vector field $u$ on $\mM$, which we can express in coordinates via the chart $(U,\psi)$ as
$u=u^\mu(x)\pa_\mu$. This expression is valid over the whole set $U$, but we are interested in its restriction along the time-like trajectory of the observer $\OO$.
We want this vector field along the curve to describe the observer's four-velocity, so we impose that the following equations hold at each point of the observer's world line:

u=d/dτ=dx^μ(τ)/dτ/x^μ ,u_μu^μ=-1 .
If we further assume that the observer path is a geodesic, then the observer's acceleration $a^\mu$ vanishes:

^ρ^=0 ,
where $\nabla_\mu$ is the covariant derivative compatible with $g$ and $\Ga^\mu_{\rho\si}$ are the Christoffel symbols specifying the Levi-Civita connection. At the background level, the only non-vanishing component of the observer four-velocity is the temporal one, which is a consequence of the symmetries of the FLRW solution. In the real inhomogeneous universe, we can describe the deviations from the background by introducing the nonlinear perturbations $\de u$ and $\VV^\al$ to the observer four-velocity

u^μ=:1/a(1+u, ^)^μ ,
with $\VV^\al$ based on $\de_{\al\be}$ (see <cit.> for notation convention).
Given this expression, we readily derive

u_η   =   g_ημu^μ=-a(1+u+2Å+2Åu+_^) ,
u_   =   g_μu^μ=a(_-_-u_+2_^)=:a _α ,
where we defined the nonlinear perturbation $\VL_\al$.
The normalization condition of the time-like vector $u_\mu u^\mu=-1$ constrains $\de u$ and $\VV^\al$ to the metric perturbations:

-1/2_^(1+2_) .
Using Eq. ($\ref{uu}$) the geodesic condition in Eq. ($\ref{geou}$) can be explicitly written as
a^η=0=^η_ηη( 1+)^2+(2^η_η^+')( 1+)+^η_^^+^_,  ,
for the temporal component, and
a^=0=^_ηη( 1+)^2+( 2^_η^+'^-^ )( 1+)+^_^^+^^_,  ,
for the spatial components.
Among all the Christoffel symbols in the geodesic condition above, only $\Ga^\eta_{\eta\eta}$ has a non-vanishing background contribution, such that $a^\mu=0$ at the background. Consequently, the terms like $\Ga^\al_{\be\ga}u^\be u^\ga$ are at least third order in perturbations. With $\de\Ga^\eta_{\eta\eta}$ we denote the nonlinear perturbation $\de\Ga^\eta_{\eta\eta}:=\Ga^\eta_{\eta\eta}-\HH$, where $\HH=a'/a$ is the conformal Hubble parameter.

We are now in a position to resume the discussion about tetrads (see <cit.> for the tetrad formalism in cosmology). In the previous subsection, we stated that we want to equip the null light path $\ga$ with tetrad basis vectors at each point; however, there are many possible ways to do this, actually infinite. Indeed, if $\{e_a\}$ is a tetrad basis at a point and we make a change of basis (not a change of coordinates) to arrive at $\{e'_a\}$, then it can be proved that $\{e'_a\}$ is still a tetrad basis only if the transformation made belongs to the Lorentz group. We recall that the finite transformations in the Lorentz group are three independent boosts and three independent spatial rotations. The physical way to fix the freedom in boosts is to align the tetrad frame with the observer trajectory, namely $e_0\equiv u$, so that $e_0$ is time-like and $\{e_i\}_{i=1,2,3}$ are space-like. We fix the residual rotation freedom by imposing that at the background level the space-like tetrad vectors are aligned parallel to the FLRW spatial directions. Expressing the tetrads in terms of coordinates $e_a=[e_a]^\mu\pa_\mu$, with the components $[e_a]^\mu$ depending only on the coordinates of the point in which the tetrad is specified, we can fix

[e_0]^μ≡u^μ ,[e_i]^μ=:1/a( e^η_i,^_i+e^_i )^μ ,
where we introduced the nonlinear perturbations $\de e^\eta_i$ and $\de e^\al_i$ <cit.>.

We can express the metric $g$ both in the coordinate-induced basis and in the tetrad basis; therefore, by relating the components in the two bases, we derive the following equation
g_μν=η_ab[e^a]_μ[e^b]_ν ,
where $[e^a]_\mu$ is the inverse of $[e_a]^\mu$: $[e^a]_\mu [e_b]^\mu=\de^a_b~$ , $~[e^a]_\mu [e_a]^\nu=\de^\mu_\nu$. Combining these relations we obtain

g_μν[e_a]^μ[e_b]^ν=η_ab .
We can use this orthonormality condition to constrain the perturbations in the components $[e_a]^\mu$, but note that this will not constrain $\de e^\al_i$ completely. For this reason, we adopt the following exact splitting

e^_i =:-^_j(S^j_i+A^j_i) ,
with $S_{ij}$ and $A_{ij}$ nonlinear perturbations based on $\de_{ij}$, which are respectively symmetric and antisymmetric under the exchange of spatial indices.
We stress that only the symmetric part $S_{ij}$ is constrained in terms of metric perturbations via the orthonormality condition <cit.>.

§.§ Light wave vector and physical observables

As in the case of the observer four-velocity, we now define a vector field $k$, which we express in coordinates as $k=k^\mu(x)\pa_\mu$. We consider the restriction of this vector field along the light path $\ga$
k=k^μ(x_Λ)_μ ,
and we demand that it describes a null vector at each point:

k=d/dΛ=dx^μ(Λ)/dΛ/x^μ ,k_μk^μ=0 .
The requirement that $k$ is not any null vector, but the observed light wave vector, uniquely specifies the affine parameter $\Lambda$ and ensures that the components of $k$ are related to physical observables.

The observer in the rest frame makes cosmological measurements, which we describe geometrically in terms of the tetrad basis at point $O$. The components of the light wave vector in this reference are

k=k^a e_a ,

k^a≡ω(1,-n^i)^a ,_ijn^in^j=1 ,

where we suppressed the subscript $o$. With $n^i$ we denote the angular position of the source in the sky of the observer and the minus sign indicates that it points from the observer to the source. With $\omega$, or better $\omega_o$, we denote the angular frequency of the light received by the observer; in coordinates it reads $\omega=[e^0]_\mu k^\mu=-u_\mu k^\mu$. To obtain the components of $k$ in coordinates, we use the tetrad components, i.e., $k^\mu=[e_a]^\mu k^a$. Using the results from the previous subsection, we obtain that at the background level, the components of the wave vector are given by

k̅^μ:=dx̅^μ(Λ)/dΛ=ω/a( 1,-n^)^μ ,n^:=n^i^ _i .
Note that this expression is valid at any point of the background light path. However, only the angular frequency $\omega$ and the scale factor $a(\eta)$ change from point to point, meaning from observer to observer. While we align the internal frames of these fictitious observers such that $n^\al$ is the same at every point of the background light path in the background manifold $\bar\mM$.

In practice, it is more convenient to introduce a conformal wave vector $\hat k$, as it simplifies equations and does not alter physical conclusions. We will discuss the details of conformal transformations in the following subsection.

§.§ Conformal light wave vector

Let us introduce on $\mM$ an additional metric $\hat g$, which is related to the metric $g$ of the real inhomogeneous universe by the conformal transformation
g_μν(η,x)=a^2(η)ĝ_μν(η,x) .
Notice that the conformal metric $\hat g$ at the background level is the Minkowski metric.
In the rest of the work we will use $\mM$ as a shorthand for $(\mM,g)$ and $\hat\mM$ for $(\mM,\hat g)$.

The convenience in introducing a conformal transformation is that we can deal with a more straightforward metric, while preserving the causal structure of $\mM$. In fact, if a vector is time-like, space-like, or null with respect to the metric $g$, it will maintain this same property for the metric $\hat g$. One can say even more, namely, if two Lorentzian manifolds have the same causal structure, then the metrics must be related by a conformal transformation <cit.>.
However, what interests us most is yet another related property of conformal transformations that we now discuss.
Let $\nabla$ be the affine connection compatible with the metric $g$, and $\hat\nabla$ the affine connection compatible with the metric $\hat g$, then we consider a curve $c(\Lambda)$ and a tangent vector $V$. We assume that $c(\Lambda)$ is an affinely-parameterized geodesic on $\mM$ hence the tangent vector satisfies

∇_V V=0 .
Relating the two covariant derivatives we can see that, in general, $c(\Lambda)$ is not a geodesic in $\hat\mM$. There is an exception in the case of a null geodesic, for which the above equation reduces to <cit.>

∇̂_V V=(-2  V^μ∇_μlna) V ,
where the scale factor $a$ arises from the conformal transformation in Eq. ($\ref{confo}$).
Notice that this equation is of the form
∇̂_V V=f V ,f:=-2  V^μ∇_μlna=d/dΛln(a^-2)
which describes a path whose tangent vector changes parallelly to itself. This feature defines a geodesic, however the above equation looks different from Eq. ($\ref{geo}$) due to the parameterization of the curve. If we reparametrize $\Lambda\to \cc$ such that
d^2/dΛ^2=fd/dΛ ,
is satisfied, then we would recover the affinely-parametrized geodesic equation Eq. ($\ref{geo}$) for $c(\cc)$ with re-scaled tangent vector $(d\Lambda/d\cc)V$. In this sense, null geodesics are conformally invariant.

We now apply all this machinery to the light path $\ga(\Lambda)$ which is a geodesic on $\mM$ with tangent vector $k^\mu$. We then introduce a new parameter $\cc$, such that $\ga(\cc)$ is an affinely parameterized geodesic on $\hat\mM$ with rescaled tangent vector $\hat k^\mu$ defined as

k^μ=k̂^μd/dΛ ,

k̂^μ:=dx^μ()/d ,

where $\hat k^\mu$ satisfies $\hat k^\nu\hat\nabla_\nu \hat k^\mu=0$ and $\hat k_\mu\hat k^\mu=0$.
By direct substitution in Eq. ($\ref{nonaffine}$) we find that $d\cc/d\Lambda$ can be obtained up to a proportionality constant $\NC$ that cannot be constrained by the conformal transformation: <cit.>

dΛ/d=a^2 .
Notice that having physically fixed the parameterization $\Lambda$ does not fix the parameterization $\cc$. It is essential to understand that $k^\mu$ is a vector under diffeomorphism in the real inhomogeneous universe $\mM$ while $\hat k^\mu$ is not, exactly because the Jacobian of the line reparameterization is a function of spacetime coordinates. We will elaborate more on this point in Section $\ref{subsec:gauge}$.

We can now express Eq. ($\ref{kbg}$) in terms of the conformal wave vector in the background:

k̂̅̂^μ=aω( 1,-n^)^μ .
Our choice to fix the constant $\NC$ is to set unity

1≡(aω)_o ,
at the observer position.
The motivation behind this choice is that this combination has the same value at every point of the background light path. We explicitly show this by evaluating the combination $\NC(a\omega)_{\bobs}$ at the background observer position $\bar O$. First, we normalize the scale factor at the observer position in the background to unity $a(\bar\eta_{\rm o})\equiv1$, and the frequency redshifts to yield $\omega_{\bobs}=(a\omega)_{\bar{\rm p}}$, where $\bar P$ is any point along the background light path. Remember that $\omega$ is the "observed frequency" in the rest frame of a fictitious observer along the light path. Hence, we derived $\NC(a\omega)_{\bobs}=\NC(a\omega)_{\bar{\rm p}}$. In the presence of perturbations, the product $(a\omega)$ varies from point to point, while $\NC$ remains constant in $\mM$.
It proves convenient to account for changes in the product $(a\omega)$ at any point $P$ of the light path by introducing another nonlinear perturbation $\dhnu$, defined as follows:

(aω)_p=-(au_μk^μ)_p=-(û_μk̂^μ)_p=:1+_p .
We introduced the conformal four-velocity components $\hat u^\mu:=au^\mu$, and we used $\omega=-u_\mu k^\mu$ derived in the previous subsection. In these terms, our choice to fix the arbitrary constant $\NC$ at the observer position can be rephrased as $\dhnu_o\equiv 0$, independent of the coordinate system. We stress that the condition ($\ref{NC}$) is evaluated at the observer position $O$ in the real universe, not at $\bar O$ in the background universe. This conformal transformation has been widely used in literature <cit.>. Another choice is to impose the normalization condition at the source position $\dhnu_s\equiv0$.

We can now discuss the components $\hat k^\mu$ of the conformal wave vector. Combining Eq. ($\ref{NC}$) and Eq. ($\ref{kbarracappuccio}$) we introduce the nonlinear parameterization

^μ=:(1+ ,-n^-^)^μ ,
defining $\dnu$ and $\dea^\al$ nonlinear perturbations, with $\dea^\al$ based on $\de_{\al\be}$.
Using the above definition, together with the normalization condition in Eq. ($\ref{NC}$), and the relation between the two wave vectors in Eq. ($\ref{khatnohat}$), we obtain the useful nonlinear parameterization for the light wave vector:

k^μ=ω_oa_o/a^2( 1+,-n^-^)^μ .
At the observer position we can derive the perturbations $\dnu_o$ and $\dea^\al_o$ in terms of the tetrad components of Eq. ($\ref{tetrade}$). From $\hat k^\mu=\NC a^2 e^\mu_a k^a$ we obtain the boundary conditions: <cit.>

_o=_o-(e^η_i)_on^i ,    -^_o=^_o-n^i(e^_i)_o
=^_o+^_j(S^j_i+A^j_i)_on^i ,
where not only the symmetric part but also the antisymmetric part contributes. Mind again the notation convention
_o=(x_o)  ,^_o=^(x_o) .
Using the conformal metric $\CG$, we can compute the exact expressions for $\hat k_\mu=\hat g_{\mu\nu}\hat k^\nu$. By direct substitutions, we derive
_η   =
-(1+2Å+-+2Å-_^) ,_   =
-(n_+_+2_+_+ _+
2_^) ,
where we introduced the notation that each index contracted with $n^{\al}$ is replaced by $\para$ (parallel to $n^\al$). In turn, the exact null equation is

k̂_μk̂^μ=0   =   -+Å----2_^+2Å-_^+1/2^2

-1/2_^+Å^2-_^-_^^ .
We conclude by assigning coordinates to the null path parameterized by $\cc$ in the same way we did earlier in Eqs. ($\ref{chart}$) and ($\ref{deviations}$):
x^μ():=(ψ^μ∘)() ,x^μ():=x^μ()-x̅^μ() .

§.§ Light path geodesic equations

We have motivated that the curve $\ga$ parameterized by $\cc$ is a null geodesic in the conformal manifold$~\hat\mM$. Now we want to express the geodesic equations in coordinates.
The temporal and spatial components are
0   =   ^μ∇̂_μ^η=dd+^η_μν^μ^ν ,0=^μ∇̂_μ^=
-dd^+^_μν^μ^ν ,
where $\hat\nabla_\mu$ is the covariant derivative based on $\hat g$, $\hat\Ga^\mu_{\nu\rho}$ are the Christoffel symbols

^μ_νρ=1/2 ^μ(_ν,ρ+_ρ,ν-_νρ,  ) ,
of the corresponding Levi-Civita connection, and the derivative with respect to the affine parameter $\cc$ is obtained from the parameterization of $\hat k^\mu$ in Eq. (<ref>)

d/d=k̂^μ/x^μ=(1+)/η-(n^+^)/x^ .
It is clear from the definition of $n^\al=n^i\de^\al_i$ in Eq. ($\ref{kbg}$) that any coordinate derivative of $n^\al$ is zero as it does not depend on spacetime coordinates.

Since $\hat g_{\mu\nu}$ in the background is chosen to be the Minkowski metric $\eta_{\mu\nu}$ in Cartesian coordinates, it automatically follows that the Christoffel symbols $\hat\Ga^\mu_{\nu\rho}$ are vanishing in the background. Hence, we can define the following perturbations that source the geodesic equations ($\ref{kgeo}$):

+^η_(n^+^)(n^+^) ,

+^_(n^+^)(n^+^) .
We emphasize that $\hat{\Ga}^\al_{\be\ga}=0$ in the background given our choice of coordinates, otherwise $\GG^\al$ has extra contributions in proportion to $\hat{\bar{\Ga}}^\al_{\be\ga}$.
Next, we integrate the geodesic equations to derive a formal solution:

(x_)-(x_o)   =   -∫_0^d' ^η(x_') ,
^(x_)-^(x_o)=∫_0^d' ^(x_') ,
where we set the affine parameter $\cc$ equal to zero at the observer position $\ga(\cc_o)=O$, $\cc_o\equiv0$.

§.§ Time lapse and spatial shift at the observer position

In the real universe, the observer position $x^\mu(\tau)$ drifts away from the background position $\bar x^\mu(\tau)$ due to the metric perturbations and peculiar velocity $\VV^\al$.
We define the time lapse and spatial shift of the observer today by evaluating the first relation in Eq. ($\ref{deviations}$) at the observer position $O$ parameterized by $\tau_o$, i.e., $\OO(\tau_o)=O$:

t(τ_o):=t(τ_o)-t̅(τ_o) ,x^(τ_o):=x^(τ_o)-x̅^(τ_o) .
We will always work in cosmic time $dt=a d\eta$ when dealing with the time lapse and spatial shift (see <cit.>).

The relation between the vector field $u$ and the associated curve $\OO$ is given by the differential equations for the integral curve, which in coordinates read

d x^μ(τ)/dτ=u^μ(x_τ)  ,    ∀τ∈[0,τ_o] ,    x^μ(0)=0 ,
which we can formally integrate
x^μ(τ_o)=∫_0^τ_odτ u^μ(x_τ) ,
to obtain the time lapse and spatial shift at the observer

t(τ_o)=∫_0^τ_odτ u(x_τ) ,x^(τ_o)=∫_0^τ_odτ 1/a^(x_τ) .
The components of $u^\mu$ in cosmic time are given by

(1+u, 1/a^)^μ ,
and we used the solution in the background universe

dτ=τ_o=∫_0^∞dzH(z)(1+z) .
Combining the last equation with the definition of $u$ in Eqs. ($\ref{udef}$) and ($\ref{ucosm}$), the background path can be parameterized by the cosmic time as
x̅^μ_τ=(τ=t,0)^μ=:x̅^μ_t .
We stress again that Eq. ($\ref{lapseandshift}$) are formal solutions, as they require the knowledge of the whole trajectory $x^\mu(\tau)$ to obtain $\de x^\mu(\tau_o)$.

§.§ Observed redshift and parameterization of the source position

Using the tetrad basis at the source position $S$, we can define the emitted light frequency $\omega_s$, and compare it with the frequency $\omega_o$ that the observer measures in the rest frame. It is well known that these two frequencies are not equal because the cosmic expansion will stretch the wavelength of the emitted signal as it travels towards the observer. The observed redshift $z$ quantifies the discrepancy between these two frequencies:
1+z:=ω_s/ω_o .
We can combine this definition with $\dhnu$ in Eq. ($\ref{dhnu}$) for $P=S$, and with the normalization condition for $\NC$ in Eq. ($\ref{NC}$) to obtain the nonlinear expression for the observed redshift
1+z=a_o/a_s(1+_s) .
We define the nonlinear perturbation $\de z$ in the observed redshift as

1+z/a_s:=1+z .
At the background level, the redshift quantifies how much the universe has expanded since the light left the source at the given scale factor $a_s$ evaluated at the coordinate position $\eta_s$ of the source. While the observed redshift $z$ is well-defined, our choice of splitting the background part and the perturbation $\de z$ is not unique. Note that the coordinate position of the source is different in an inhomogeneous universe from that in the background, and it is also gauge-dependent. Furthermore, we remark that $z$ is a scalar, but $\de z$ does not have a well-defined geometric nature.

Combining the last two equations, we obtain the nonlinear expression for the perturbation $\de z$ in the observed redshift

z=a_o(1+_s)-1 ,
and the nonlinear expression for $\dhnu_s$ from its definition in
Eq. ($\ref{dhnu}$) is

_s   =   2Å+_-_+(1+2Å)(++)+^( _-_ )+^_(1+)

+2^(_+_^) .
We stress again that the right-hand side is evaluated at $x^\mu_s$ but there is also a dependence of $\dhnu_s$ on the reference position $x^\mu_o$ via $\dnu$ and $\dea^\al$; thus, we will use the notation $\dhnu(x_s;x_o)$. This argument also applies to $\dnu(x_s;x_o)$, $\dea^\al(x_s;x_o)$, $\de z(x_s;x_o)$ and $z(x_s;x_o)$.

To illustrate the significance of the observed redshift, we now turn our attention to the relation between $\hat k$ and $\gamma$, given by the differential equations for the integral curve

d x^μ()/d=k̂^μ(x_)  ,    ∀∈[0,_s] ,    x^μ(0)= x^μ_o .
For the moment, we consider the background level, so by integrating the differential equation, we obtain
η̅_-η̅_ o= ,x̅^_=-n^ ,
which shows how the affine parameter is related to the time coordinate of the background light path. However, since coordinates are defined up to diffeomorphisms, there exists ambiguity. We can use the observed redshift $z$ to define an unambiguous time coordinate $\bar\eta_z$, along the light path in the background: <cit.>

1+z=:1/a(η̅_z) ,
and this time coordinate is associated to the affine parameter $\cc_z$ defined by

η̅_z-η̅_ o=:_z=-∫_0^zdz'/H(z')=-_z ,    x̅^_z:=-_zn^=_z n^ ,
where $\rbar_z$ is the comoving distance out to redshift $z$ in the
background. Noting that the observed redshift and the position $\bar x^\mu_z$ are independent of coordinates (or gauge choice), we will use this physical parameterization to denote the reference position of the source in the real universe:

η_s=η̅_s+δη_s=:η̅_z+Δη_s ,    x^_s=x̅^_s+δx^_s=:x̅^_z+^_s ,
_s=:_z+_s ,
where the subscript $z$ indicates evaluation at $\cc_z$, and the subscript $s$ at $\cc_s$. Note that with lowercase delta $\de$ we indicate perturbations around the background source position parameterized by $\cc_s$, while we use uppercase delta $\Delta$ to indicate perturbations around the background source position parameterized physically by the observed redshift. These two types are in general different e.g. $\de\eta_s\neq\DT_s$.
We want to stress that the parameter $\cc_z$, unlike $\cc_s$, has been introduced to parameterize the light path in the background, where the source position $\bar S$ is identified precisely by the scalar $z$. If we want to use the parameterization $\cc_z$ on $\mM$, we need to account for the perturbation $\Dcc_s$ to identify $S$.

With these definitions, we can obtain useful nonlinear relations. At the source position, we derive

^μ_s=x^μ_s+x̅^μ_s-x̅^μ_z=x_s^μ+_s(1,-n^)^μ ,
where we used $\bar x^\mu_{\rm s}=\left( \cc_s+\bar\eta_{\rm o},-\cc_sn^\al \right)^\mu$ and $\bar x^\mu_z=\left( \cc_z+\bar\eta_{\rm o},-\cc_zn^\al \right)^\mu$.
While at the observer position, we obtain the boundary relations

(x_o;x_o)=0 ,    _o=0 ,    _o=_o≠0 ,    _o^=δx^_o≠0 .
Using the relation between $\de z$ and $\bar\eta_z$ in Eqs. (<ref>) and (<ref>), we derive the exact relation
between $\DT_s$ and $\dz$:

1+=a(η̅_z+_s)a(η̅_z) .
In addition to the observed redshift, we can use the observed angular position $n^i$, parameterized by the angles $\theta,\phi$, to fully specify the real location of the source in terms of quantities defined in the rest frame of the observer. These angles are related to $n^i$ by the following expression in Cartesian coordinates
n^i=(sinθcosϕ,sinθsinϕ,cosθ)^i .
We also define two unit vectors perpendicular to $n^i$ in the rest frame of the observer:

^i   :=    ∂/∂θn^i=( coscos,cossin,-sin)^i ,
^i   :=   1/sinθ ∂/∂ϕ n^i=(-sin,cos,0)^i ,
which can be defined on $\mM$ in the same way we did for $n^\al$, i.e., $\ttt^\al:=\de^\al_i\ttt^i$ and $\pp^\al:=\de^\al_i\pp^i$.
The parameterization we adopt for the source position is then <cit.>

x^μ_s=x̅^μ_z+^μ_s≡x^μ(z,θ,ϕ)   =   (η̅_z+Δη_s ,  (_z+)sin( +)cos( +) ,

(_z+)sin( +)sin( +) ,

(_z+)cos( +))^μ ,
with $\DT_s$, $\drr$, $\dtt$, $\dpp$ nonlinear perturbations that are functions of $(z,\theta,\phi)$. For the same reason above, the perturbations $(\de r,\dtt,\dpp)$ should be denoted as $(\Delta r_s,\Delta\ttt_s,\Delta\pp_s)$ but we keep this notation convention as there is no possibility of confusion in this case.
We stress again that $z,n^i,\ttt^i,\pp^i$ are observables, so their values do not depend on our coordinate-based description, and they are invariant under diffeomorphisms of $\mM$.

§.§ Distortion in the source position compared to the observed position

Let us consider again the relation between the conformal wave vector $\hat k$ and the light path $\gamma$ given in Eq. ($\ref{sourcepos}$).
Isolating the perturbative parts, it reads

d/dx^μ=(,-^)^μ ,
valid at any point along the light path.
Integrating this equation yields the formal solution

δx^μ_s   =   (_s , δx^_s)^μ=(η_o+∫_0^_s d  , x_o^-∫_0^_s d

=   (_o+λ_s _o-
∫_0^_s d(_s-)^η ,
x^_o-_s ^_o-∫_0^_s d(_s-)^)^μ ,
where we used the solution for $\dnu$ and $\dea^\al$ in Eq. ($\ref{kint}$).
Combining the last expression with Eq. ($\ref{Dxcc}$), we derive the following nonlinear expressions for the distortion in the source position with respect to $\bar x^\mu_z$:

_s   =   _o+_s+_s _o-
∫_0^_s d(_s-)^η(x_) ,

^_s   =   x^_o-_s n^-
_s ^_o-∫_0^_s d(_s-)^(x_) .
Solving Eq. ($\ref{De}$) for $\Dcc_s$, and utilizing the solution for $\dea^\al$ in Eq. ($\ref{kint}$), we derive a more convenient nonlinear expression for the spatial distortion at the source position:

_s^=x_o^-n^(_s-η_o)-∫_0^_sd( ^-n^) ,
where we kept the expression in a compact form by leaving unexpressed the perturbations $\dea^\al$ and$\,\dnu$. Note that $\cc_s$ in the integration limit also has a perturbative part $\Dcc_s$.

§.§ Physical volume and area occupied by the source

So far, we have discussed how the observer can infer the position of a point-like light source via the observed redshift $z$ and the observed angles $\ttt , \pp$. Now we turn our attention to the case of some extended light source subtended by the observed solid angle $\up d\Omega$, over the observed redshift range $\up dz$. Our interest is to provide expressions for the proper volume and area in terms of the observed quantities in the rest frame of the observer. The former is related to galaxy clustering, and the latter is related to the luminosity distance. Perturbative calculations have been performed in the past <cit.>

Although we are considering an extended source, it can still be geometrically described by the point $S$ on the manifold $\mM$. Indeed, we can define the volume and area in terms of a volume form $\Si$, i.e., a top-form defined at the source position $S$. Having introduced a tetrad basis $\left\{ e_a\right\}$ at the source position, we can define the volume form as the following exterior product:

:=e^0∧e^1∧e^2∧e^3 ,
where $\left\{ e^a\right\}$ is the dual basis of one-forms satisfying $e^a(e_b)=\de^a_b$. We use the tetrad components $[e^a]_\mu$ to recast this volume form at point $S$ in FLRW coordinates, and then using Eq. ($\ref{source}$), we relate the coordinates of the source position to the observed redshift and angles. Note that the parameterization in Eq. ($\ref{source}$) is valid at point $S$, when the proper time of the source is $\tau'=\tau_s$. Hence, we derive the volume form in coordinates <cit.>
=√(-g) ε_ρμν u^ρ_s ∂x^μ_s∂∂x^ν_s∂
∂x^_s∂ dτ_s∧d∧d∧d ,
where the Levi-Civita symbol is normalized as $\varepsilon_{\eta xyz}=1$, the (square root) of the metric determinant

√(-g):=√(-det g_μν)=:a_s^4(1+g) ,
arises to express the volume form $\Si$ in the FLRW basis $\left\{dx^\mu\right\}$ at the source position,
$u^\mu_s$ are the components of the source velocity $u_s=d/d\tau_s$, and $d\tau_s,dz,d\ttt,d\pp$ are the one-forms built from the corresponding parameters, and we defined the dimensionless perturbation $\de g$.
Considering small intervals $\up dz ,\up d\ttt, \up d\pp$ and fixed proper time $\tau_s$, we define the three-dimensional volume element $dV$ as the components of the volume form:

(τ_s,dzz,d,d)=:dV=√(-g) ε_ρμν u^ρ_s ∂x^μ_s∂∂x^ν_s∂
∂x^_s∂ dz d d ,
where we use the roman $\up d$ to denote intervals, while the italics $d$ denoted one-forms. Note that the expression above is a scalar under diffeomorphisms: $\up dz,\up d\ttt, \up d\pp$ are observables, so by construction, they are not affected by coordinate transformations, then $u_s^\mu$ and $x^\mu_s$ transform as components of four-vectors and their transformation is compensated by that of the scalar density $\sqrt{-g}$. Hence, Eq. ($\ref{volume}$) indeed describes the physical volume occupied by the source in its rest frame, that appears subtended by $(\up dz,\up d\ttt,\up d\pp)$ in the observer rest frame.

Given this definition of physical volume, there are many ways to introduce an area by considering different decompositions of $\up dV$. Here we define the physical area as the split of $\up dV$ into an element $\up dA$ perpendicular to the light propagation in the rest frame of the source. To derive the expression in FLRW coordinates of the light propagation, let us consider the light wave vector $k$ at the source position. Since at point $S$ there are both a tetrad basis $\left\{e_a\right\}$ and a coordinate-induced basis $\left\{\pa_\mu\right\}$, the following equality holds
k^ae_a=k=k^μ_μ ,
where we suppressed the subscript $s$. Substituting the expression for the internal components of the light wave vector $k^a$ in Eq. ($\ref{kobs}$), we obtain the propagation direction at the source position

n^μ_s:=-k^μ_sω_s+u^μ_s=1/a_s(-k̂^μ_s1+_s+û^μ_s) ,    (n^μn_μ)_s=1 ,     (n^μu_μ)_s=0 ,
where we additionally made use of the prescription at the source position $[e_0]^\mu\equiv u^\mu$, and via Eq. ($\ref{dhnu}$) we derived the second equality involving conformal four-vectors. The directional vector is normalized and spatial in the source rest frame. With this expression at hand, we define the physical area from Eq. ($\ref{volume}$):

dV   =   :dA (n_μx^μz)_sdz ,
dA   =   √(-g) ε_ρμν^ρ_s n^μ_s
d d .
The same comment about diffeomorphism invariance applies to this quantity, which indeed corresponds to the physical area occupied by the source in its rest frame. We emphasize again that both covariant expressions ($\ref{volume}$) and ($\ref{area}$) contain derivatives with respect to observed quantities $z, \ttt, \pp$, which affect directly only the observed angles $\ttt^\al, \pp^\al$, while FLRW coordinates are involved via the physical parameterization of the source position $x^\mu_s$ in Eq. ($\ref{source}$). The physical volume $\up dV$ additionally contains the derivative with respect to the observed redshift $z$, which describes the derivative along the line-of-sight of the observer. Finally, we stress that these are partial derivatives, so they are performed while keeping the other observables fixed.

We are now in the position to consider the perturbations in the physical volume and area due to the light propagation in the real universe $\mM$. We explicitly write $\up dV$ in Eq. ($\ref{volume}$)
using Eq. ($\ref{source}$) to expand the source coordinates, Eq. ($\ref{uu}$) for the components of $u^\mu_s$, and Eq. ($\ref{dzdz}$) to express the scale factor at the source position $a_s$ in terms of the observed redshift $z$ and the fluctuation $\de z$:

dV   =   (1+)^3(1+)^3(1+){
ε_ηαβγ(1+u) ∂(x̅_z^α+^) ∂

+ε_μνρ ^ ∂(x̅_z^μ+^μ)∂
}dz d d ,
where we suppressed the subscript $s$.
The background contributions to the homogeneous and isotropic universe indeed reduce to the background volume element
dV̅=_z^2H_z(1+)^3 dz  dΩ ,
and we define the dimensionless volume fluctuation $\de V$ as

dV=:_z^2H_z(1+)^3 dz  dΩ (1+)  ,
with $H_z=H(\bar\eta_z)$ and $\up d\Omega=\sin\ttt\,\up d\ttt\up d\pp$. From this equation we manifestly see that $\de V$ is a scalar under diffeomorphisms. Equations ($\ref{dV}$) and ($\ref{deV}$) provide the exact definition of the gauge-invariant fluctuation $\de V$ in the physical volume.

The relation between the physical area occupied by the source $\up dA$ and the observed solid angle element $\up d\Omega$ defines the angular diameter distance $\dD_A$:

dA=^2_AdΩ ,
which we can split into the background and the perturbative parts as

_A=:_A(1+) ,_A=_z/1+z ,
valid at any order.
Solving Eq. ($\ref{area}$) for the angular diameter distance using Eqs. ($\ref{splitDA}$) and ($\ref{lumD}$) we obtain

^2   =   _A^2
sin  (1+)(1+)^2[ε_ηαβγ(1+u)(n^+^1++^)

+ε_η ^(--1++u )

+ε_νρ ^(n^+^1++^)
=   _A^2(1+)^2=:_A^2(1+A) ,
where again we suppressed the subscript $s$, and we substituted the expression for the light propagation direction $n_s^\mu$ in Eq. ($\ref{propag}$) together with the equations mentioned for the expansion of $\up dV$. Moreover, we defined the perturbation $\de A$ in analogy with the volume fluctuation $\de V$.
We recall that the angular diameter distance is proportionally related to the luminosity distance $\dD_L$, and it can be decomposed in the same manner:
_L=(1+z)^2_A .
This relation suggests that the fluctuations in the angular and diameter distances are identical, and the perturbation in the physical area occupied by the source can be expressed as

A=^2+2 .
Again, we observe that this fluctuation, just like $\de V$, is a scalar under change of coordinates of $\mM$. Equation ($\ref{deA}$) provides the exact definition of the gauge-invariant fluctuations $\de\dD$ and $\de A$.

For later convenience, we introduce the following decompositions along the line-of-sight and observed angular directions
^   =   _n^+_θθ^+_ϕϕ^ ,    ^=_n^+_θθ^+_ϕϕ^ ,
with $\ttt^\al$ and $\pp^\al$ defined via Eq. ($\ref{angoli}$).

§.§ Gauge transformations of cosmological quantities

The study of gauge transformations of second-order perturbations represents the main effort of this work.
It will provide a validity check of our second-order expressions presented in Section $\ref{sec:second}$, and it will be carried out explicitly in Section $\ref{sec:gauge}$. Here we present a general introduction to gauge transformations in cosmology, and we establish the notation conventions relevant for the dedicated section (see <cit.> for a general discussion).

Recall that the theory of general relativity is based on the principle of general covariance, namely we are free to choose any coordinate system to describe physics. One consequence of this freedom is that quantities that have been geometrically defined as tensors must obey well-defined transformation laws under a change of coordinates.
Let us formalize the problem by considering a chart $(V,\varphi)$ of $\mM$, alternative to $(U,\psi)$, such that $V$ overlaps with $U$ and contains the light path $\ga$. Given a vector field $\zeta$, we can use the transformations that it induces to characterize the change of coordinates at a point $P$ in $U \cap V$:

x^μ_p(x_p)=e^ζ^ν_νx^μ|_p=x^μ_p+ζ^μ(x_p)+1/2 ζ^ν(x_p)_νζ^μ|_x_p+⋯ .
We use a tilde to denote the coordinate system given by the chart $(V,\varphi)$, defined in the same way as in Eq. ($\ref{chart}$). Note that we can continue to use the parameterization $\cc$, $\tau$, and $\tau'$ for null and time-like curves since the parameterization procedure is chart-independent.

Let us consider a tensor $K$ of arbitrary rank at point $P$. It is well known that the following exact transformation law holds

K^μ⋯ν_ρ⋯(x_p)=(∂x^μ/∂x^μ_1)_p⋯(∂x^ρ_1/∂x^ρ)_p⋯ K^μ_1⋯ν_1_ρ_1⋯_1(x_p) ,
which relates the components of $K$ in the new frame, indicated by a tilde, to the components of $K$ in the original frame. We stress that both sides of the equation are evaluated at point $P$, which in the tilde frame has coordinates $\widetilde x^\mu_p$, and in the original frame has coordinates $x^\mu_p$. If the tensor we are considering is a scalar, the above transformation law trivializes:

K(x_p)=K(x_p) .
In cosmology, we do not consider a generic coordinate transformation because we want to exploit the correspondence of the spacetime $(\mM,g)$ with the symmetric background solution given by FLRW $(\bar\mM,\bar g)$. Changing this correspondence between the real universe and the homogeneous isotropic universe is what we call gauge transformation. It is possible to relate a gauge transformation to a change of coordinates by evaluating both sides of Eq. ($\ref{tensor}$) at the same (numerical) value $x^\mu_p$ in the two frames, which corresponds to two different points of $\mM$. This procedure defines the gauge transformation of the tensor $K$, as generated by the vector field $\zeta$:

K^μ⋯ν_ρ⋯(x_p) = e^-_ζK^μ⋯ν_ρ⋯|_p=K^μ⋯ν_ρ⋯(x_p)- (_ζK^μ⋯ν_ρ⋯)(x_p) +1/2 (^2_ζK^μ⋯ν_ρ⋯)(x_p)+⋯ ,
with $\pounds_\zeta$ being the Lie derivative with respect to $\zeta$. Instead of working with the gauge vector field $\zeta$, we introduce an alternative vector field $\xi$ that simplifies the relation between new and old coordinates given in Eq. ($\ref{xzeta}$):

x^μ_p(x)=:x^μ_p+ξ^μ(x_p) .
The vector field $\xi$ will be our preferred choice to compute gauge transformations, and it is assumed to be a nonlinear perturbation.
We adopt the following parameterization for its components

ξ^μ=:(T,^)^μ ,
with $T$ and $\LL^\al$ nonlinear perturbations.
The conversion between the two gauge fields $\zeta$ and $\xi$ can be readily derived and is presented at the second perturbative order in Section $\ref{sec:gauge}$.

Note that the whole discussion above applies only to tensors. However, not all the quantities that we have introduced so far are tensors under diffeomorphism in $\mM$. The only tensor fields we defined are the metric $g$, the observer velocity $u$, the light wave vector $k$ and therefore all the tensor products that we can build. In particular, the redshift $z$, the perturbation in the volume occupied by the source $\de V$ and in the luminosity distance $\ddL$ are scalars under diffeomorphisms. Conformally rescaled quantities such as $\hat g$ or $\hat k$ instead, do not transform as tensors. To further illustrate this point let us consider the action of a generic diffeomorphism on the metric $g$:
g_μν(x)=x^ρ/x^μx^/x^νg_ρ(x) .
Given the conformal transformation $g=a^2\hat g$, we can derive from the equation above how the conformally rescaled metric $\hat g$ transforms:
ĝ_μν(x)=a^2(η)/a^2(η)x^ρ/x^μx^/x^νĝ_ρ ,
which is not how a rank-two tensor transforms. The same argument applies to $\hat k$ and in turn affects the perturbations $\dhnu$ and $\de z$.
Once again, the importance of computing gauge transformations in this work is that they provide constraints that our expressions must satisfy, and thus we use them as a theoretical consistency check. We will continue this discussion in more detail in the dedicated section.

§ SECOND-ORDER PERTURBATIVE CALCULATIONS

Given the definitions and exact relations of the previous section, we now perform perturbative calculations up to second order.
Let us consider some cosmological quantity $Q$ defined in the real universe $\mM$ and assume the split into background and perturbation as $Q=\bar Q+\delta Q$. We adopt the convention

Q =Q^(1)+Q^(2)+⋯ ,
where the superscripts represent the perturbative order.
For the sake of clarity, we will omit the superscripts for quadratic terms, i.e., terms made by the product of two first-order variables. The linear terms without superscripts, such as $\de Q$, represent first-order and second-order terms.

In the following, we will provide second-order perturbative expressions for the observer and source position, the observed redshift, and the area and volume occupied by the source; these will be obtained perturbing the expressions presented in the previous section. We will begin by computing the perturbations of geometric quantities to finally express all the second-order observables in terms of metric and observer velocity perturbations only. Note that we are not assuming Einstein equations or any other field equations that describe the dynamics of the universe. Hence, metric and observer velocity perturbations cannot be constrained at this stage and are treated as independent spacetime fields.

§.§ Metric-related quantities and observer four-velocity

A comprehensive second-order computation of the metric and the Einstein equations can be found in <cit.>.
We can compute the second-order expression of the inverse metric components perturbing the relation $g_{\mu\rho}g^{\rho\nu}=\de^\nu_\mu$. We obtain

g^ηη   =   1 a^2 ( - 1 + 2Å- 4 Å^2 + _α^α) ,      g^η=1 a^2 ( - ^α+ 2 Å^α+ 2 _β^αβ ) ,
g^   =   1 a^2 ( ^αβ - 2 ^αβ
-^α^β+ 4 ^α_γ^βγ ) .
The second-order perturbation in the determinant of the metric in Eq. ($\ref{detg}$) is[It can be computed perturbing the following exact equation
det g=(1/24)[(tr g)^4-6tr(g^2)(tr g)^2+3(tr(g^2))^2+8tr(g^3)tr g-6tr(g^4)] .

=Å+^α_α-12Å^2+1/2^α_α+Å ^α_α+12^α_α^_-
^αβ_αβ .
From the definition of $\hat\Ga^\rho_{\mu\nu}$ in Eq. ($\ref{gammahat}$) we derive the second-order Christoffel symbols defining the Levi-Civita connection based on $\hat g$
^η_ηη   =   Å'-2ÅÅ'+^('_-Å_,) ,
^η_η   =   Å_,-2ÅÅ_,-^('_-_[,]) ,
^η_   =   '_+_(,)-2Å( '_+_(,) )+^(_,-2_(,)) ,
^_ηη   =   -^'+Å^,+Å'^+2^('_-Å_,) ,
^_η   =   '^_-^[_,]+^Å_,+2^(_[,]-'_) ,
^_   =   -_^,+2^_(,)+^('_+_(,))
+2^(_,-2_(,)) .
Consequently, we can obtain $\GG^\eta$ and $\GG^\al$ at second order from Eq. ($\ref{GGt}$) and Eq. ($\ref{GGa}$)

^η   =
Å'-2+_,+'+2(Å' - )
- 2Å(Å'-2+_,

-_(α,) -_α')
-^(Å_,α-'_α+2_[,] -2'_
+ 2 _α,
-_,) ,
^   =   Å^,α-'^α
-  _^  ,α + ^α_  ,
- 2 '^α_+
2 ^α_, - _^    ,α
+(2Å^,α - 2'^α-_^  ,α
+ ^α_  ,- 2 '^α_ )
-^(_^  ,α
- ^α_  ,+ 2 '^α_-4 ^α_(,γ)+2_γ^    ,α
)+ ^α(Å'+ _, + '_- 2 Å_,)

- 2^α(  2 _,
- _,+Å_,
-'_ -2_[,γ]
- 2_γ')
where round brackets around indices represent the symmetrization and square brackets the antisymmetrization with respect to those indices, e.g., $2\BB_{(\al,\be)}=\BB_{\al,\be}+\BB_{\be,\al}\,$ , $\,2\BB_{[\al,\be]}=\BB_{\al,\be}-\BB_{\be,\al}$.

To obtain the perturbed expression of the observer four-velocity we can use the time-like normalization $u_\mu u^\mu=-1$ in Eq. ($\ref{delu}$). We derive up to second order <cit.>

=- Å+3/2 Å^2 + 12
^α_α-_^ ,
and the covariant components $u_\mu$ by using Eq. ($\ref{delubasso}$):

_η= - a ( 1 + Å-1/2Å^2
+ 12 ^α_α) ,    _α= a ( _α- _α+Å_α+ 2 ^β_αβ )=a_ .
As our approach is purely geometric, the observer velocity $\VV^\al$ and the metric perturbations $\AA,\, \BB_\al,\, \CC_{\al\be}$ cannot be constrained.
All our perturbative expressions presented in this section will ultimately depend on these independent spacetime fields.

§.§ Subtleties in second-order computations

Second-order perturbative calculations are more challenging than linear order calculations because of the many subtleties that arise.
Let us take the time lapse at the observer position as an example to clarify this point (see, e.g. <cit.>).
We begin by considering the differential equation ($\ref{curvint}$) for the integral curve. This equation conveniently relates $\de t$ to $\de u$ so that we can exploit the second-order expression in Eq. ($\ref{du2}$) to compute the lapse in terms of metric perturbations and observer velocity only.
To perturb this differential equation we make the following consideration: $\de t(\tau)$ and $\de u(x_\tau)$ are well-defined quantities, and their perturbative expansion can be generically expressed according to Eq. ($\ref{exp}$):
t(τ)   =    t^(1)(τ)+t^(2)(τ) +⋯ ,
u(x_τ)   =   u^(1)(x_τ)+u^(2)(x_τ) +⋯ .
However, the perturbative expansion of $\de u(x_\tau)$ is somewhat ambiguous, as there exist different ways to define it. This quantity is in fact a functional of the observer world line, so it can be expanded around the background path $\bar x^\mu_\tau$ as
u(x_τ=x̅_τ+x_τ)=u(x̅_τ)+x^μ_τ(/x^μu)|_x̅_τ+⋯ ,
such that the distortion in the path appears explicitly. We now make the following identification between perturbative orders, that we assume to be valid for any field on $\mM$

,u^(2)(x_τ)=u^μ(2)(x̅_τ)+x^μ(1)_τ(/x^μu^(1))|_x̅_τ .
Note that at linear order the deviations from the background path $\bar x^\mu_\tau$ do not come into play as $\de u$ is already a first-order perturbation. We can then write the formal perturbative solutions of Eq. ($\ref{curvint}$) as

t_o=∫_0^t̅_odt u(x̅_t) ,
at linear order, and

t_o=∫_0^t̅_odt [u(x̅_t)+x^μ_t(_μu)|_x̅_t ] ,
at second order. To emphasize the fact that the integrations are performed along the background path we use the parameterization $t$ instead of $\tau$. The reason why we decide to expand each field around the background is that this allows us to better keep track of the perturbative orders, and in this way, we can actually obtain the solution $\de x^\mu_\tau$ at each perturbative order, which is then fed into a higher order solution. This choice will also prove convenient when we will compute gauge transformations.

Again, we take the time lapse at the observer position to illustrate another subtle point.
We have repeatedly stated that $\de t(\tau)$ is not a field but rather a (perturbed) time coordinate.
However, it is also possible to define the time lapse as a field along the observer world line. By specifying a point on the trajectory $x^\mu_\tau$ we are also specifying a value of the field $u^\mu(x_\tau)$ and from Eq. ($\ref{curvint}$) this also means specifying the value $\de t(\tau)$. This procedure defines a map $\de t(x_\tau)$ from the coordinates of the path to the corresponding time lapse such that $\de t(x_\tau):=\de t(\tau)$ for all $\tau$ in $[0,\tau_o]$.[In fact, we can think of the time lapse for another trajectory of a fictitious observer. In this way, it can be defined as a field on $\mM$.] Using this field definition of the time lapse, we can replace $\de u$ in Eq. ($\ref{dt2}$) by $\dot{\de t}$, as the partial derivatives are now justified. We can now integrate by part, thus obtaining

t_o=x_o^μ(_μt)|_+∫_0^t̅_odt [u(x̅_t)-ẋ^μ_t(_μt)|_x̅_t ] ,
where the initial conditions in Eq. ($\ref{curvint}$) have been used. We further remark that the integral in the equation above is indeed the second-order expression for $\de t(\bar x_{\rm o})$, i.e., the time lapse for another observer at $\bar x^\mu_{\rm o}$ in an inhomogeneous universe. Since in the real universe the observer is located at $O$, and not at $\bar O$, we will need the expression for $\de t_o$ only, and not at all for $\de t_{\bobs}$. Therefore, $\de t_o$ will not be expanded around $\bar x^\mu_o$ in any of the following equations while the perturbative expression in the right-hand side of Eq. ($\ref{dt2nogeo}$) is expanded around the background path. The same arguments apply to $\de x^\al_o$, $\dnu_o$, $\dea^\al_o$, $\dnu_s$, $\dea^\al_s$ and all the other quantities we will compute. We conclude by stressing that we will only use the field-definition of path deviations $\de x^\mu$ to derive simpler equations for formal solutions, particularly to perform integration by part. On any other occurrence, they are treated simply as functions of the affine parameter denoting the path.

§.§ Time lapse and spatial shift at the observer position

Resuming the discussion, we now devote ourselves to computing perturbative expressions for the time lapse and spatial shift at the observer position. We will first derive the linear-order expression and then the second-order one. We adopt this structure throughout the whole section.

Linear-order lapse and shift
Combining the expression for $\de u$ in Eq. ($\ref{du2}$) and the linear-order solution for the time lapse at the observer in Eq. ($\ref{dt1}$) we obtain at linear order

δt_o=-∫_0^t̅_odt Å(x̅_t) ,δx^_o=∫_0^t̅_o dt 1/a^(x̅_t) .
Note that since the integrand is already first order in perturbations, the fields are evaluated at $\bar x^\mu_t$, such that the integration is along the time coordinate $t$ without any spatial shift $\bar x^\al_t=0$.
We will use the parameterization $dt$ to emphasize that the integration is performed along the background path.

If the observer experiences no acceleration, the geodesic equations ($\ref{geou}$) provide extra linear-order conditions

Å=/t(av) ,v_=0 ,
where we introduced the Helmholtz decomposition $\VL_\al:=v_\al-v_{,\al}$, with $\pa^\al v_\al=0$.
Using these geodesic conditions, the linear-order time lapse at the observer is further simplified as <cit.>

t_o=-v(x̅_o) .
Moreover, we derive the linear-order relations that are valid at any point along the background path:

/tδt=-Å(x̅_t) ,/tδx^=1a^(x̅_t) .

Second-order lapse and shift
Combining the second-order expression for $\de t_o$ in Eq. ($\ref{dt2nogeo}$) with $\de u$ in Eq. ($\ref{du2}$) and Eq. ($\ref{dt1nogeo}$) we derive

To compute the spatial shift at the observer we can follow the same steps to obtain at the second order

x^_o=x_o^μ(_μx^)|_+∫_0^t̅_odt 1/a(^+Å ^-^x^_,)|_x̅_t
In both expressions, any occurrence of $\de x^\mu$ is intended as the linear-order expression in Eq. ($\ref{dxlin}$), even if we do not perform the explicit substitution. In this sense the partial derivatives $\pa_\nu\de x^\mu$ are well-defined.

The second-order expression for the time lapse can be further simplified by assuming that the observer follows a geodesic. In particular, the linear-order condition in Eq. ($\ref{dt1geo}$) at any point of the background path reads $\de t=-av$. Using this equation we can rearrange Eq. ($\ref{dt2met}$) as
(Åv-δx^v_,)_ .
Moreover, we derive the second-order geodesic equations from the definition in Eq. ($\ref{geou}$):
a_=0   =   Å_, + 1/a[a(_- _)]^'+Å(-2Å_,+2'_+2_-'_-_)

+ Å^'_+2'^__+2^_( -Å_,+'_+_ )+ ^( _, - 2_[,])  ,
a^η=0   =   (^-^)a_ .
Noting that the linear-order geodesic condition in Eq. ($\ref{lindt}$) can be recast as $\AA_{,\al}=-(1/a)\left( a\VL_{\al} \right)^\prime$, we can rewrite the previous geodesic condition by using $\VL_\al$ in Eq. ($\ref{ubasso}$):
a_=0=Å_,+1/a(a_)^'( 1+Å)+^( _,+_, ) .
Furthermore, by considering $\VV^{\be}\BB_{\be,\al}=\VV^{\be}\VV_{\be,\al}-\VV^{\be}\VL_{\be,\al}$ we manipulate the equation
0   =Å_,+1/a(a_)^'+2^_[,]-ÅÅ_,+^_,
=(a_)^·+[ Å-1/2Å^2+1/2^_ ]_,,
which can be readily solved by the Helmholtz decomposition as
Å-1/2Å^2+1/2^_-( av )^·=0 ,v_=0 .
Finally, we manipulate the integrand of Eq. ($\ref{dt2met}$)
-1/2^_ ,
such that by combining this equation with the previous one, we derive the second-order time lapse at the observer with geodesic motion

t_o=-(v-Åv+δx^v_,)_ .
This result is consistent with the one derived in <cit.>.
Due to our use of conformal time, it proves more convenient the observer time lapse $\de\eta_o$ in terms of the conformal time. Given the definition of the cosmic time $dt=a d\eta$, we integrate it to obtain the matching condition for $\eta_o=\bar\eta_o+\de\eta_o$:

t_o=t̅_o+t_o=∫_0^η_odη a(η)   =   t̅_o+η_o+1/2(η̅_o)η^2_o+⋯ ,
and this yields the second-order expression for the time lapse in conformal coordinates

η_o=t_o-1/2_t_o^2 .

§.§ Distortion in the conformal wave vector

Having set up the observer position at second order, we proceed to compute the perturbations in the conformal wave vector $\hat k^\mu$ along the light path using the geodesic equations ($\ref{kint}$).

Linear-order solution
By substituting $\GG^\eta$ and $\GG^\al$ in Eq. ($\ref{GG2}$) into the formal geodesic solutions in Eq. ($\ref{kint}$), we obtain at linear order

|^s_o   =   ( -2Å)|^z_-∫_0^_z d(Å--)' ,
_|^s_o   =   -( _+2_ )|^z_-∫_0^_z d(Å--)_, .
Recalling the arguments presented in Subsection ($\ref{subsec:sub}$), the integration is performed along the background path, parameterized by the comoving distance $\rbar$ instead of the affine parameter $\cc$. In fact, from Eq. ($\ref{dr}$) the following relation holds at the background
d/d=-d/d=η-n^x^ .
We can obtain the linear-order expressions for $\hat k_o^\mu$ using the tetrad basis at the observer position.
From the definition in Eq. ($\ref{tetrade}$) and the orthonormality condition in Eq. ($\ref{orto}$) we derive the space-like tetrad components at linear order
e_i^η=^_i(_-_) ,S_ij=^_i^_j_ ,
where $S_{ij}$ is the symmetric part of the spatial components $\de e_i^\al$ in Eq. ($\ref{split}$). To derive the antisymmetric part $A^j{}_i$, we cannot rely on the orthonormality condition, but we can find constraints based on gauge transformation properties <cit.>. The derivation is shown in Appendix $\ref{AppendixA}$, here we state the result:

_jA^j_i=C_[,]^_i+ε_ij^j ,
where $C_\al$ is the vector part of $\CC_{\al\be}$ in its scalar-vector-tensor decomposition, $\varepsilon_{ijk}$ is the Levi-Civita symbol with $\varepsilon_{123}=1$, and $\Om^i$ is a gauge-invariant vector describing the misalignment of the space-like tetrad basis against a FLRW coordinate <cit.>.
Recall that in Eq. ($\ref{tetrade}$) we fully specified the tetrad basis only at background level, where we required the components of the space-like tetrads to be fully aligned with FLRW spatial directions. In the real universe $\mM$, we fixed the three degrees of freedom associated with boosts, but we still have three degrees of freedom associated with spatial rotations that are specified by $\Om^i$, but not determined by metric perturbations. Combining these results with the definition of $\hat k^\mu_o$ in Eq. ($\ref{ko}$) we derive the linear-order distortions at the observer

_o=-( Å+_-)_ ,(_)_o=-(_+_+C_[,]+ϵ_i jn^i^j)_ .

Second-order solution
To obtain the second-order solution we need to expand the formal expression in Eq. ($\ref{kint}$) around the background. This implies expanding the fields as functionals of the light path, but also expanding the upper limit of integration $\cc_s$ around affine parameter $\cc_z$ via $\cc_s=\cc_z+\Dcc_s$:

|^s_o=∫_0^_zd[^η(x̅_)+x^μ_(_μ^η)|_x̅_]-_s^η(x̅_z) .
We can simplify this equation using the geodesic equation ($\ref{kint}$) $\frac{d}{d\cc}\dnu=-\GG^\eta$, and then integrate by part.
Using the perturbations for $\de\hat k^\mu$ in Eq. ($\ref{dxder}$) and for $\Delta x_s^\mu$ in Eq. ($\ref{Dxcc}$) we arrive at

|^s_o=[^μ_(_μ)|_x̅_]^z_+∫_0^_zd(^η+'-^_,)|_x̅_ ,
where the boundary terms on the right-hand side appear as the difference between the terms inside the square brackets evaluated at $\rbar=\rbar_z$ and $\rbar=\rbar_{\rm o}$.
We can repeat the same steps for $\dea^\al$ and derive

^|^s_o=[^μ_(_μ^)|_x̅_]^z_-∫_0^_zd(^-'^+^^_,)|_x̅_ .
To obtain these expressions in terms of perturbations of the metric and velocity, we should substitute the linear-order solutions for $\dnu(\bar x)$ and $\dea^\al(\bar x)$ in Eq. ($\ref{dk1s}$), then the second-order expressions for $\GG^\eta(\bar x)$ and $\GG^\al(\bar x)$ in Eq. ($\ref{GG2}$), as well as the expression for $\DX^\mu_s$ (which we will derive later) and $\de x^\mu_o$ from the previous subsection. The complete expression is presented in Appendix $\ref{AppendixB}$.

Again, we can compute the components of $\hat k^\mu_o$ using the tetrad at the observer in Eqs. ($\ref{tetrade}$) and ($\ref{split}$). Employing the orthonormality condition ($\ref{orto}$) we derive at second order

e^η_i   =   _i-_i+Å(2_i-_i)+_i(^+^)+ε_ij(^-^)^j+(^-^)C_[i,] ,
S_ij   =   _ij+1/2_i_j-1/2_i_j-1/2_i_j+1/2_ij^k_k-3/2^_i_j-^_(iε_j)k^k-1/2C_,(i^_j)

-1/4C^_,(iC_j), ,
where spacetime indices have been contracted with Kronecker delta involving internal indices. From these relations, we derive the second-order expression for $\dnu_o$ defined in Eq. ($\ref{ko}$)

_o   =   (- Å-_++3/2Å^2-_( ^+^)+Å( _-2)-^_
+1/2_^+(^-^)C_[,]+ε_i j( ^-^)n^i^j)_+x^μ_o[_μ(-Å-_+_)]|_ ,         where the last term is the expansion needed to have the right-hand side evaluated at the background observer position $\bar x^\mu_o$. The second-order expression for $\dea^\al_o$ can also be computed from the definition in Eq. ($\ref{ko}$), taking into account the above expression for $S_{ij}$

-^_o   =   ^_o+(A^_i+^_i+1/2_i^-1/2_i^-1/2_i^+1/2^_i^k_k-3/2^_i^_-^_(iε^)_k^k

-1/4C^_,(iC^)_, )_n^i
+x^μ_o[_μ(^+^_in^i+C^[_,i]n^i+ε^_ijn^i^j)]|_ .
Note that the second-order expression for the antisymmetric part $A^j{}_i$ is not presented here, as it does not contribute to any of the second-order quantities we explicitly present in Appendix $\ref{AppendixB}$.

§.§ Fluctuations in the observed redshift

The observed redshift $z$ in Eq. ($\ref{dzdz}$) is affected by the cosmic expansion set by $a_s$ at the source position and the fluctuations $\de z$ that include the Doppler effect and the Sachs-Wolfe effect. Here we compute the second-order expression for $\de z$ by using Eq. ($\ref{dz}$).
Following the steps in Section $\ref{subsec:sub}$, we expand the fields on the right-hand side in Eq. ($\ref{dz}$) around the background reference positions $\bar x_z$ and $\bar x_{\rm o}$.
The scale factor at the observer position is
a(η_o=η̅_o+η_o)=1+_η_o+1/2(^2_+'_)η_o^2 ,
and we expand the perturbation $\dhnu$ around the observer position $\bar x_z$ as
We choose not to expand the reference position around $\bar x^\mu_o$, keeping it fixed at the real observer position $O$.
Combining these results, we obtain the second-order expression for the fluctuation in the observed redshift

z(x_s;x_o)   =   _η_o+(x̅_z;x_o)+1/2(^2_+'_)η_o^2
+^μ_s(_μ)|_x̅_z .        In order to express $\de z$ in terms of the metric perturbations and the observer velocity, we first need to derive a second-order expression for $\dhnu$. This is readily obtained from the exact relation in Eq. ($\ref{dhnuexact}$):

(x̅_z;x_o)   =   +Å+_-_+Å-1/2Å^2+Å_+_( ^-^)
+1/2 ^_+2_^ ,        where both sides are evaluated at $\bar x^\mu_z$, with reference $x^\mu_o$. The full expression for $\de z$ is given in Appendix $\ref{AppendixB}$ (see also <cit.>).

§.§ Distortion in the source position compared to the observed position

In Eq. ($\ref{xz}$) we introduced two ways to define deviations in the source position, $\de x_s^\mu$ and $\DX^\mu_s$. The first is a deviation compared to the unphysical background position $\bar x^\mu_s$, while the latter is compared to the background observed position $\bar x^\mu_z$. Our interest is indeed in computing the second-order perturbative expression for $\DX^\mu_s$, since it contains the entire chart dependence of the source position $x^\mu_s$.

Linear-order source position
To compute the temporal distortion $\DT_s$, we expand up to first order the exact relationship with the redshift $\de z$ in Eq. ($\ref{dzrel}$). We obtain

_s=z/_z=1/_z(_η_o++Å+_-_) ,
where we substituted the linear-order expression for $\de z$ computed in the previous subsection.
For the spatial components of the distortion $\DX^\al_s$, from the exact relation in Eq. ($\ref{DXal}$) we derive instead

_s^=x_o^-n^(_s-η_o)+∫_0^_zd( ^-n^) ,
as well as the following linear-order relation valid at any point along the background light path

d/d( ^+n^)=(^-n^) .
We stress that the left-hand side of this equation is evaluated at a given comoving distance $\rbar$, while the right-hand side at the background position $\bar x^\mu_{\rbar}$. In analogy with the observer lapse discussed in Subsection $\ref{subsec:sub}$, we can define the coordinate distortion $\DX^\mu(\bar x_{\rbar})$ as a field along the background path, instead of a function of affine parameter.

Given the physical parameterization of the source position in Eq. ($\ref{source}$), we can use the previous results to compute the radial and transverse distortions:

=   n_α^α_s=n_x_o^-_s+∫_0^_zd _   =
x_o +_o-_z+∫_0^_zd (
Å--) ,

_z   =   _α^α_s=_x_o^+∫_0^_zd  _^   =
-∫_0^_zd [θ_α(^α+2 ^α_)

+_x_o^+_z_α(^+^α+2 ^α_)_ ,

_zsin   =   _α^α_s=_x_o^+∫_0^_zd  _^   =
-∫_0^_zd [_α(^α+2 ^α_)

+_x_o^+_z_α(^+^α+2 ^α_)_ ,
where we used $\bar x^\al_z=\rbar_z n^\al(\ttt,\pp)$ to recast the partial derivatives.

Second-order source position
Expanding up to second order the exact relation between $\DT_s$ and $\de z$ in Eq. ($\ref{dzrel}$) we obtain
z=_z_s+1/2(^2_z+'_z  )^2_s ,
again, we used the fact that $\de z$ has functional dependence only on $x^\mu_s$, not on $x^\mu_o$.
Substituting the linear-order expression for $\DT_s$ we can solve the above equation for the second-order temporal distortion:

_s=z/_z-1/2^2_z+'_z/_z(z/_z)^2 .
For the spatial distortion $\DX^\al_s$ we perturb the exact definition in Eq. ($\ref{DXal}$) in the same way we did for $\dnu$ and $\dea^\al$, i.e., we expand around the background both the integrand and the integration limit. We obtain the following result
_s^   =   ∫_0^_zd[ ^-n^+x^μ__μ( ^-n^)|_x̅_]

+x_o^-n^(_s-η_o)-_s( ^-n^)|_x̅_z .
We can further simplify this expression using Eq. ($\ref{dx1}$) to replace $\dea^\al-n^\al\dnu$ evaluated along the background path:

_s^   =   ∫_0^_zd[ ^-n^+x^μ__μd/d(^+n^)

+x_o^-n^(_s-η_o)-_sd/d(^+n^)|_x̅_z ,
where the field definition of $\de x^\mu$ has been used.
We can now integrate by part and, making use of Eqs. ($\ref{Dxcc}$) and ($\ref{dxder}$), we can manipulate the above equation to finally obtain

_s^   =   ∫_0^_zd[ ^-n^+('^+n^')-^(^_,+n^_,)

+x_o^-n^(_s-η_o)+[_^μ _μ( ^+n^)|_x̅_]^z_ ,
where the partial derivatives of $\DX^\mu$ act on the linear-order expression, and the square brackets represent the difference between the terms evaluated at $\rbar_z$ and $\rbar_{\rm o}$.

We can proceed to compute the radial and transverse distortions as before. By a direct computation we find

=   n_^_s+_z/2 ( ^2+sin^2 ^2 ) ,
_z   =   _^_s-+_z/2cossin ^2 ,
_zsin   =   _^_s-sin -_zcos  ,
where again we do not make any substitution, but we derived explicit expressions in terms of metric and observer velocity perturbations for each term appearing.

§.§ Distortion in the physical volume and area

We conclude the section by presenting the second-order expressions for the fluctuations in the volume and area occupied by the source in its rest frame compared to the observed volume and area. The relevant nonlinear equations are Eq. ($\ref{deV}$) for $\de V$ and Eq. ($\ref{deA}$) for $\de A$. The distortion in the observed proper area of the source $\de A$ is related to the perturbations $\de\dD$ in the luminosity distance or angular diameter distance via Eq. ($\ref{dA}$).

For convenience we introduce the distortion matrix $\DD$, usually employed in the standard formalism of gravitational lensing to compare the source angular position $(\ttt+\dtt,\pp+\dpp)$ with the observed angular position $(\ttt,\pp)$ <cit.>:
:=∂(θ+θ,ϕ+ϕ)/∂(θ,ϕ) .
Given this definition, we can compute the determinant of the distortion matrix

det 𝔻=sin(+)/sin( 1+++-)=:1+ ,
where the nonlinear perturbation $\de\DD$ has been introduced. Here we introduced the distortion matrix $\DD$ merely to simplify the second-order expression for $\de V$ and $\de A$, which would be too lengthy otherwise. At linear order, the perturbation $\de\DD$ reads
𝔻=(+)+=: -2κ ,
with $\kappa$ being the standard gravitational
lensing convergence describing the linear-order lensing magnification. Already apparent at the linear order, $\de\mathbb{D}$ is gauge-dependent. At second order we obtain

𝔻   =   (+)++-+ (+)
-1/2^2 .To derive the second-order expressions for $\de V$ and $\de A$ we first need to expand up to second-order the physical parameterization of the source position in Eq. ($\ref{source}$), which together with
zη̅_z=-1/H_z ,z_z=1/H_z ,
obtained from the definition of the comoving distance $\rbar_z$ in Eq. ($\ref{dr}$), allows us to recast the terms in Eq. ($\ref{dV}$) for $\up dV$.
From direct computation we obtain
ε_ηαβγ ∂x^α∂
=   ^2_zsinH_z{1+2 _z+𝔻
+H_z∂∂ +^2^2_z+2 _z
(H_z∂∂ -2 κ)
-2κ H_z∂∂-H_z∂∂  ∂∂ -H_z∂∂  ∂∂ } ,

ε_μνρ ^ ∂x^μ∂
=   ^2_zsinH_z{_(1+2 _z
-2 κ-H_z∂∂ )+_ +_sin
(+)} ,
where we suppressed the subscript $s$. Substituting these equations in the definition of $\de V$ in Eq. ($\ref{deV}$) yields

=    3++2_z++H_zz++_+3 +3 ^2
-Å(2 _z-2κ+H_zz)

+(+3 )
+^μ_μ( 3 +-Å+_) ,
where the left-hand side is evaluated at the source position $\de V(x_s)$, while the right-hand side is evaluated at the background reference position $\bar x^\mu_z$. The last bracket expresses the expansion of the fields around these coordinates; once again, we point out that $\drr$, $\dtt$, $\dpp$, and $\DT$ are not fields, so they do not contribute to such expansion. To keep the expression as concise as possible, we did not express all the second-order terms and most of the linear-order ones.
A further complication in the expression above is that spacetime fields and observed angles mix; both derivatives concerning observables and spacetime coordinates also appear. These aspects will be treated in more detail when we consider gauge transformations in the next section.

We can proceed in the same way to compute the distortion in the observed proper area element of the source. From direct computation, we derive
ε_μνρ ^(a^μ)
∂x^ν∂∂x^ρ∂   =
(_∂∂+_sin∂∂)} ,

∂x^β_s∂∂x^γ_s∂   =   _z^2
_+-+2 _z++2 _z (_+_-)

-2 κ(2 _z+_+_-)
+^2^2_z-1/2^2-12(sin )^2

)_z+(-_)} ,          which can be further simplified using the linear-order null condition in Eq. ($\ref{null}$) and the linear-order expression for $\dhnu$ in Eq. ($\ref{dhnu2}$). Substituting these results in Eq. ($\ref{deA}$) for $\de A$ we obtain

A   =    2++2_z++_++_--Å(-2κ+2/_z-)

-1/_z( _+_/sin )-2κ(2/_z-)-2/_z+^2/^2_z-1/2^2-1/2(sin )^2

+1/_z+sin (_+_)+1/_z+z^2+2 z
+( 2 z+)(  -Å+2/_z-2κ-)
+^μ_μ( 2 z+-Å+_+_-) ,        The same remarks made for $\de V$ apply equivalently to $\de A$. Similar expressions are derived in <cit.>. In Section $\ref{sec:compa}$ we compare our results with previous work.

§ SECOND-ORDER GAUGE TRANSFORMATIONS

To test the validity of the second-order expressions derived in the previous section, we will now turn our attention to the gauge transformation properties of each of those expressions.
Resuming the discussion of Subsection $\ref{subsec:gauge}$, we can compute the gauge transformation of a tensor using Eq. ($\ref{gheigi}$), which at second order reads

K^μ⋯ν_ρ⋯(x_p=x_p) = K^μ⋯ν_ρ⋯(x_p)- (_ζK^μ⋯ν_ρ⋯)(x_p) +1/2 (^2_ζK^μ⋯ν_ρ⋯)(x_p)  .
Recall that we evaluate both sides of the equation at value $x_p^\mu$ so that the background contributions of the tensor $ \bar K^{\mu\cdots \nu}_{\rho\cdots\si}(t_p)$ are identical in the two frames.
As anticipated, the gauge field $\zeta(x)$ is treated as a perturbation, and it is expanded according to Eq. ($\ref{exp}$):
ζ^μ(x)=ζ^μ(1)(x)+ζ^μ(2)(x)+⋯ ,
the same kind of expansion holds for the gauge field $\xi^\mu(x)$, which will be our preferred choice since it makes computations simpler. From Eqs. ($\ref{xzeta}$) and ($\ref{xcsi}$), the second-order relation between the two gauge fields is given by

ζ^μ=ξ^μ-1/2 ξ^νx^νξ^μ ,
which can be expressed in components using Eq. ($\ref{xicomp}$):
ζ^η=T-1/2 T'T  -1/2 T_, ^
,ζ^=^-1/2 '^T -1/2 ^_, ^ .
With these relations we can readily recast the second-order gauge transformation in Eq. ($\ref{ancoragheigi}$) in terms of our preferred gauge field $\xi^\mu$. Since we will deal only with scalars, vectors, and symmetric tensors with two covariant indices, the only transformation laws we need to consider are:

K   =   K-ξ^μK_,μ+ξ^μξ^ν_,μK_,ν+1/2ξ^μξ^νK_,μν
K^μ   =   K^μ-_ξK^μ+1/2^2_ξK^μ+1/2ξ^ρξ^ν_,ρK^μ_,ν-1/2(ξ^ρ_,νξ^μ_,ρ+ξ^ρξ^μ_,νρ)K^νK_μν   =   K_μν-_ξK_μν+1/2^2_ξK_μν+1/2ξ^ξ^ρ_,K_μν,ρ+K_ρ(μ[ ξ^_,ν)ξ^ρ _,+ ξ^ξ^ρ_,ν)]
where the Lie derivative is expressed in terms of $\xi^\mu$ and both sides are evaluated at the same coordinate value $x^\mu_p$. Notice that while a scalar field does not transform under a change of coordinates, it indeed transforms non-trivially under a gauge transformation. In general, we are interested in the gauge transformation of the perturbation variables. Splitting a tensor $K$ into the background and the perturbation $K=\bar K+\de K$ we derive the gauge transformation law for the perturbations:
K   =   K-K̅'T-ξ^μK_,μ+K̅'T_,μξ^μ+1/2K̅” T^2
K^μ   =   K^μ-K̅'^μ(T-T_,νξ^ν)+K̅^ν(ξ^μ_,ν-ξ^ρξ^μ_,νρ)+ ξ^μ_,νK^ν-ξ^νK^μ_,ν-K̅'^νTξ^μ_,ν+ 1/2K̅”^μT^2
K_μν   =   K_μν-2K̅_ρ(μξ^ρ_,ν)-TK̅_μν-ξ^ρK_μν,ρ-2K_ρ(μξ^ρ_,ν)+2K̅_ρ(μ[ξ^ρ_,ν)ξ^+ξ^_,ν)ξ^ρ_,]
+K̅_ρξ^ρ_,μξ^_,ν+2K̅'_ρ(μξ^ρ_,ν)T+1/2 K̅”_μνT^2+K̅'_μνT_,ρξ^ρ ,
where background quantities are purely time-dependent due to the symmetries of the FLRW solution.

In Eq. ($\ref{xz}$) we parameterized the source position in terms of the observed position $\bar x^\mu(\cc_z)$ parameterized by the observed redshift and angles
x̅^μ(_z)=(η̅_z,_z n^)^μn^=(sincos, sinsin,cos)^ ,
and its deviation $\DX^\mu$ around the observed position.
Here we emphasize why this is more advantageous than the splitting around $\bar x^\mu(\cc_s)$. As we have extensively discussed, the observed redshift $z$ and the observed angles $\ttt,\pp$ are independent of our FLRW coordinate description; hence they are scalars. From Eqs. ($\ref{etaz}$) and ($\ref{dr}$) we note that these scalars fully determine the background source position $\bar x^\mu_z$. Therefore, the reference coordinates $\bar x^\mu_z$ will not be affected by diffeomorphisms of $\mM$. On the contrary, the coordinates $\bar x^\mu_s$ of the source will be affected by the coordinate transformation, so the parameterization of the source position in terms of deviations around these coordinates is not a convenient choice to deal with gauge transformations. Given the
boundary relations in Eq. ($\ref{z0}$), the background observer position $\bar x^\mu(\cc_o)=(\bar\eta_{\rm o},0)^\mu$ also remains unaffected by a change of coordinates. With these remarks, we can readily see that the transformations induced by the gauge field $\xi$ in Eq. ($\ref{xcsi}$) will involve only the deviations $\DX^\mu_s$ and $\de x^\mu_o$, such that we can readily derive the second-order gauge transformation

x^μ_o=x_o^μ+ξ^μ(x̅_o)+x^ν_o(/x^νξ^μ)|_x̅_o ,    where the fields are expanded as usual around the background reference points.

In the rest of the section, we will use gauge transformations as a consistency check of the second-order expressions derived in Section $\ref{sec:second}$. First, we compute their expected gauge transformations employing the nonlinear relations and exact definitions given in Section $\ref{sec:nonlinear}$. Second, we directly compute the gauge transformations of the second-order expressions in Section $\ref{sec:second}$ in terms of metric perturbations and observer velocity. The consistency of these two independent procedures provides a strong validity test for our second-order expressions in Section $\ref{sec:second}$. For instance, the gauge transformation of $\DX^\mu_s$ in Eq. ($\ref{trasfx2}$) is the expected transformation from its definition. We can then use its expression in Eqs. ($\ref{deta2}$) and ($\ref{dx2}$) to check the consistency of the gauge transformation.

§.§ Metric and observer velocity perturbations

The gauge transformations of the metric perturbations and the observer velocity will play a crucial role in the consistency check because all the second-order expressions in Section $\ref{sec:second}$ finally depend on these quantities. We first apply the second-order gauge transformation equation ($\ref{gauge2}$) to the metric $g_{\mu\nu}$ in Eq. ($\ref{metrica}$) and the observer velocity $u^\mu$ in Eq. ($\ref{uu}$). After straightforward computations, we obtain the second-order gauge transformations of the metric perturbations <cit.>

Å   =    Å-(T^' +T )- Å^'T
- 2 Å(T'+T)+32T'^2
+TT^'' +1 2
( 2^2+') T^2
+ 3T^'T-Å_,α^α-_α'^α
+^α(T^'_,α+T_,α )
- 12'^α ^'_α ,

_α   =   _α-T_,α+ ^'_α- 2ÅT_,α- (^'_α+ 2 _α) T
- _αT^'+2 T^' T_,α
+T (T^'_,α+ 2T_,α )
-_α,β^β-_β^β_  ,α+ 2_αβ '^β
+ T_,β^β_  ,α+^T_,α-T ( ^''_α+2^'_α)

_(,)'^-^'_, ,

_αβ    =    _αβ-T_αβ-_(,β)
+_(α T_,β)-(^'_αβ+
2 _αβ )T
+ 1 2 (
2^2+')  T^2 _αβ  -12T_,αT_,β

+ TT^'_αβ -_αβ,γ ^γ- 2_γ(α ^γ_  ,β) + ^'_(α
T_,β)+^γT_,γ_αβ + 2_(,β)T+
'_(,β) T

+^_,( _),+^_(,β)+1/2^_,_, ,
where we used Eq. ($\ref{xicomp}$) to express the components of the gauge field $\xi^\mu$. Remember that we chose rectangular coordinates to simplify the calculations. In the same way, we obtain the second-order gauge transformation of the observer velocity

u   =   u+T+T'+(𝒜'  -  𝒜 -T”)T+ 12 (ℋ^2  -  ℋ' )T^2 -  𝒜T'   + 𝒰^ T_,β

+ ^( 𝒜_, -  ℋ T_,β -  T'_,β)  ,
^   =   ^+ '^+(T-Å) '^+(^-
”^-𝒰'^ )T   + 𝒰^ ^_,β

-  ^β( 𝒰^_,β+'^_,β) .
We can readily check that Eq. ($\ref{deutilde}$) is consistent with the transformation obtained by directly substituting the gauge transformations of the metric perturbations derived above in the second-order expression for $\de u$ in Eq. ($\ref{du2}$).

Note that both sides of each equation presented in this subsection are evaluated at the same coordinate position $x^\mu$. For later convenience, we will assume these coordinates to be either on the background observer path $\bar x^\mu_t$ or on the background light path $\bar x^\mu_{\rbar}$.

§.§ Time lapse and spatial shift at the observer position

We begin the consistency checks by computing the second-order transformations of the time lapse $\de t_o$ and spatial shift $\de x^\al_o$ at the observer position. On the one hand, we compute the transformations from the definitions of the lapse and shift in Eq. ($\ref{lapseshift}$). On the other hand, we compute the transformations of the second-order expressions in Eqs. ($\ref{dt2met}$) and ($\ref{dx2met}$) in terms of metric and observer velocity perturbations. The former procedure will be denoted second-order expectation, while the latter second-order direct computation. We employ this terminology throughout the whole section. By checking that the two procedures give the same result, we ensure the sanity of our second-order expressions. We will only show the second-order gauge-transformation, from which the linear-order gauge transformation can be trivially read off.

Second-order expectation
Let us write in components the transformation of the time lapse and spatial shift in Eq. ($\ref{trasfx2}$) in the observer position:

η_o=η_o+T(x̅_o)+x^μ_o(/x^μ T)|_x̅_o ,
x^_o=x_o^+^(x̅_o)+x^μ_o(/x^μ^)|_x̅_o .    To derive the expectation for the time lapse in proper time, let us use the relation between $\de t_o$ and $\de\eta_o$ in Eq. ($\ref{lapseconf}$). Using the gauge transformation of $\de\eta_o$ above, we derive the expectation
t_o   =   η_o+1/2_η_o^2=t_o+T_+x^μ_o[x^μ (aT)]|_+1/2_T^2_ .

Second-order direct computation
We start the direct computations from the spatial shift at the observer position in Eq. ($\ref{dx2met}$). Directly substituting the gauge transformations of the metric perturbations and the observer velocity presented in Subsection $\ref{subsec:4.1}$, and using the linear-order transformations of the coordinate perturbations in Eq. ($\ref{trasfx2}$), we obtain
δx^_o   =    δx_o^μ(
_μδx^)|_+∫_0^t̅_odt 1/a
(^+Å ^-^δx^_,)
=   δx^_o+x^μ_o(_μ^)|_+ξ^μ_[_μ( x^+^)]|_

+∫_0^t̅_odt1/a[/η(^-T^-T'^-^^_,)-^^_,-'^x^_,]|_x̅_t ,
where the boundary terms in the first line come exclusively from the gauge transformation of $\de x_o^\mu(\pa_\mu \de x^\al)\big\rvert_{\bobs}$.
From the linear-order condition in Eq. ($\ref{dt1nogeo}$), valid at any point of the background observer path $\bar x^\mu_t$, we can rewrite the last two terms in the above integral as
-^^_,-'^x^_,=-/η(^x^_,  ) .
With this result and the relation $dt=a\,d\eta$, we notice that the integrand above is a partial derivative with respect to cosmic time $t$, and the integration yields boundary terms evaluated at the background observer position. Most of the terms canceled in the resulting expression to yield
δx^_o=δx^_o+^_+x_o^μ(x^μ^)|_ ,
which indeed coincides with the second-order expectation.

We can proceed in the same way for the time lapse at the observer position. We have derived two expressions for the second-order lapse, the general expression in Eq. ($\ref{dt2nogeo}$) and the simplified expression in Eq. ($\ref{dt2geo}$) under the assumption that the observer motion is geodesic. We will compute the direct transformation of the general expression:
δt_o   =    δx^μ_o(_μδt)|_+∫_0^t̅_odt
=   δt_o+x^μ_o[_μ(aT)]|_+ξ^μ_[_μ(t+aT)]|_

+∫_0^t̅_odt [1/a/η( aT+aÅT-a^T_,
-aTT'-1/2 aT^2 )+^Å_,
-1/a'^t_, ] .      Again, we can use the linear-order condition in Eq. ($\ref{dt1nogeo}$), valid on the observer background path, to rearrange the last two terms in the integral
-1/a'^t_,=-1/a/η( ^t_,  ) ,
such that the integration yields boundary terms at the observer position. Simplifying the resulting expression, we obtain
t_o=t_o+T_+x^μ_o [x^μ(aT)]|_+1/2_T^2_ ,
in agreement with the expected result (see <cit.>).

§.§ Distortion in the conformal wave vector

Here we want to compute the second-order gauge transformation of the components of the conformal light wave vector $\hat k^\mu$ defined in Eq. ($\ref{hatk}$). In the real universe $\mM$, we have argued that the conformal light wave vector is not an actual vector, and its components $\dnu$ and $\dea^\al$ depend on the position on the light path $x^\mu_\cc$ and the reference position $x^\mu_o$. As far as gauge transformations are concerned, we are considering an active change of coordinates only for the position on the light path, while the reference is fixed at the physical point $O$ in both frames.[Note that if we also expanded $\widetilde x^\mu_o$ around $x^\mu_o$, this position would not correspond to the point $O$ in the frame denoted by a tilde.] For this reason, in this section, we will consider the perturbations as a function of the light path only: $\dnu=\dnu(x_\cc)\,,\, \dea^\al=\dea^\al(x_\cc)$.

Second-order expectation
To compute the gauge transformation of $\hat k^\mu$, we first compute the gauge transformation of $k^\mu$ parameterized as in Eq. ($\ref{kmu}$), since the light wave vector $k^\mu$ is a four-vector in $\mM$ while $\hat k^\mu$ is not. Mind the following subtle point: the left-hand side of Eq. ($\ref{kmu}$) is evaluated at a point $x^\mu_\Lambda$ in the real universe, while the right-hand side depends on $x^\mu_\cc$ and $x^\mu_o$. Consequently, a change of coordinates in a point of the real universe implies a change of coordinates in two points of the conformal universe.

To compute the second-order expectation of $\widetilde\dnu$ and $\widetilde\dea{}^\al$ we employ the transformation law for the components of a tensor evaluated at a point $P$, given in Eq. ($\ref{tensor}$). This equation provides the expression for a tensor of arbitrary rank; in the case of a four-vector, it reads

k^μ(x_p)=(∂x^μ/∂x^ν)_pk^ν(x_p) .
In principle, we could also apply Eq. ($\ref{gauge2}$) directly to derive the sought gauge transformations, but with this procedure we will derive intermediate relations that will be crucial for the consistency check. Substituting in the above equation the parameterization of $k^\mu$ in Eq. ($\ref{kmu}$) we derive

1+(x)   =   a^2(η)/a^2(η)a(η_o)/a(η_o)∂η/∂x^ν(1+(x) , -n^-^(x))^ν   =   a^2(η)/a^2(η)a(η_o)/a(η_o)[(1+T')(1+)-T_,(n^+^)](x) ,
for the zero-th component $\mu=\eta$, and

-n^-^(x)   =   a^2(η)/a^2(η)a(η_o)/a(η_o)∂x^/∂x^ν(1+(x) , -n^-^(x))^ν   =   a^2(η)/a^2(η)a(η_o)/a(η_o)['^(1+)-(^_+^_,)(n^+^)](x) ,
for the spatial components $\mu=\al$. In both equations, we suppressed the subscript $p$, and we used $\widetilde\omega_o=\omega_o$, as the observed frequency does not depend on our coordinate description. Note that these equations are valid at any perturbative order at this stage and $n^\al$ is the observed angle in the rest frame of the observer, not a function of position $x^\mu$. By expanding $\widetilde x^\mu$ in Eq. ($\ref{xcsi}$) around $x^\mu$ in the right-hand side, we derive at the second order:

(x)    =   -d/dT +(2T-_oT_o)(1+-d/dT) + (2 ℋ^2 +  ℋ')T^2

+T'- T_,^-2 ℋ ℋ_o T T_o +12 ( ℋ_o^2 - ℋ'_o)T_o^2

^(x)   =   ^+d/d^+(2T-_oT_o)(n^+^α+d/d^) -'^+^_,^
+ (2 ℋ^2 +  ℋ' ) T^2n^α+12 ( ℋ_o^2 -  ℋ'_o)T_o^2n^- 2 ℋ ℋ_o  T T_o n^α ,
which will be used to check the validity of our expressions. Note that the left-hand side is evaluated at $\widetilde x^\mu$, while the right-hand side at $x^\mu$. To obtain the gauge-transformation relation at the same value, we first expand the left-hand side around $x^\mu$, move the quadratic terms to the right-hand side, and finally use the linear-order relation:

=   -d/dT + Td/dT' +(2T-_oT_o)(1+-d/dT)

-  ℒ^ (_,+2 ℋT_,-d/dT_,  )+ (2 ℋ^2 -  ℋ')T^2 -  2ℋTT'

+T'-T_,^-T'-2 ℋ ℋ_o T T_o +12 ( ℋ_o^2 - ℋ'_o)T_o^2

^   =   ^+d/d^-Td/d'^+(2T-_oT_o)(n^+^α+d/d^)

-  ^β( ^_,+2 ℋ  T_,βn^α+d/d^_, )-2TT'n^+ (2 ℋ^2 -  ℋ' ) T^2n^α
-'^+^_,^-T'^+12 ( ℋ_o^2 -  ℋ'_o)T_o^2n^- 2 ℋ ℋ_o  T T_o n^α .
As a concluding remark we stress again that the coordinates of the reference position have not been expanded as they refer to the physical point $O$. In practice, this means that in the left-hand side the equations read $\widetilde\dnu(\widetilde x=x;\widetilde x_o)\,,\,\ \widetilde\dea^\al(\widetilde x=x;\widetilde x_o)$, while the right-hand sides contains only fields evaluated at $x^\mu$ and $x^\mu_o$.

Second-order direct computation
Our second-order expressions for $\dnu$ and $\dea^\al$ are given in Eq. ($\ref{dnuint}$) and
Eq. ($\ref{deaint}$) respectively. We first check the consistency of the expression for $\dnu$ by computing the gauge transformation

(x_s)-(x_o)   =   [^μ_(_μ)|_x̅_]^z_+∫_0^_zd(^η+'-^_,)|_x̅_
^μ_μ(+ 2T-d/dT )+T'+^_,-Td/dT'

+∫_0^_zd(^η+'-^_,)|_x̅_ ,
where only the boundary terms have been expanded, for which we used the gauge transformation of $\DX^\mu$ in Eq. ($\ref{trasfx2}$) and the linear-order expression of $\dnu$ in Eq. ($\ref{dk1s}$). We present the gauge transformation of the integral separately to avoid lengthy expressions:

∫_0^_z d ( ^η+'-^_, )|_x̅_=

=∫_0^_z d ( ^η+'-^_, )|_x̅_
+{-d/dT+Td/dT'+2T( 1+-d/dT )
-2TT'-^( _,+2T_,-d/dT_, )
-_oT_o( -d/dT )+( 2^2-' )T^2
-T'+T'-T_,^-2_oTT_o}|^z_ ,
where we used $\GG^\eta$ in Eq. ($\ref{GG2}$), $\dnu$ and $\dea^\al$ in Eq. ($\ref{dk1s}$), and we stress that all the fields inside the curly brackets are evaluated along the background light path $\bar x^\mu_{\rbar}$. Combining the two previous equations ($\ref{dnutilde1}$) and ($\ref{dnutilde2}$), we obtain after straightforward simplifications:

(x_s)-(x_o)   =   (x_s)-(x_o)+{-d/dT+2T+^μ_μ( 2T-d/dT )

-_oT_o)( -d/dT )+(2^2+')T^2

-2_oTT_o}|^z_ .
The second-order expectation for this quantity is given by evaluating Eq. ($\ref{dnuPtilde}$) at the source position $P=S$ and the observer position $P=O$. The difference of these two contributions yields

(x_s)-(x_o)   =   (x_s)-(x_o)+{-d/dT+2T
+(2T-_oT_o)( -d/dT )+(2^2+')T^2
-2_oTT_o}|^s_o .
There seems to be a discrepancy between the two equations but note the different evaluation positions. The fields in the curly brackets in Eq. ($\ref{direct}$) are evaluated at coordinates $\bar x^\mu_z$ and $\bar x^\mu_o$ along the background paths $\bar x^\mu_{\rbar}$. On the other hand, the fields in the curly brackets of the expectation in Eq. ($\ref{expect}$) are evaluated at coordinates $x^\mu_s$ and $x^\mu_o$ along the real path $x^\mu_\cc$.
When we discussed subtleties at second order in Subsection $\ref{subsec:sub}$, we argued that evaluating a field either on the real path or on the background path makes a difference at second order, and the discrepancy is given in Eq. ($\ref{2bar}$). Hence, we need the following second-order relations to compare the two equations properly:
d/dT|_x=d/dT|_x̅+^μ[_μ(d/d T)]|_x̅ ,    (T)(x)=(η̅)T(x̅)+^μ_μ(T)|_x̅ ,  where we do not show any subscript as these equations are valid both for source and observer positions. Moreover, with $\DX^\mu$ we refer both to $\DX^\mu_s$ and $\de x^\mu_o$.
Substituting these relations in Eq. ($\ref{expect}$), we readily check that the expectation and the direct computation agree and indeed the apparent discrepancy was just a consequence of different evaluation positions.

We can repeat the same procedure for the spatial perturbation $\dea^\al$ given in Eq. ($\ref{deaint}$). The transformation of the boundary terms is given by
[^μ_(_μ^)|_x̅_]^z_   =
^μ_μ(^+ 2Tn^+d/d^)+T'^+^^_,+Td/d'^
+2'T^2n^+2TT'n^+2^T_,n^}|^z_ ,
while for the integral we derive
∫_0^_z d ( ^-'^+^^_, )|_x̅_=
=∫_0^_z d ( ^-'^+^^_, )|_x̅_
-{d/d^+2T( n^+^+d/d^)
-Td/d'^-2TT'n^+( 2^2-' )T^2n^-^( ^_,+2T_,n^+d/d^_, )
-_oT_o( ^+d/d^)-T'^-'^+^_,^-2_oTT_on^}|^z_ .
The substitutions employed to obtain these equations are the same mentioned in the previous case of $\dnu$, and we do not repeat them here.
The sum of these two contributions leads, after straightforward simplifications, to
^(x_s)-^(x_o)   =   ^(x_s)-^(x_o)+{d/d^+2Tn^+^μ_μ( 2Tn^+d/d^)

-_oT_o)( ^+d/d^)-2_oTT_on^
+(2^2+')T^2n^}|^z_ .
As before, the expectation for this equation can be read from Eq. ($\ref{dnuPtilde}$), and we find no discrepancy once the following second-order relations are taken into account
d/d^|_x=d/d^|_x̅+^μ[_μ(d/d ^)]|_x̅ ,    (T)(x)=(η̅)T(x̅)+^μ[_μ(T)]|_x̅ .

§.§ Fluctuations in the observed redshift

The observed redshift $z(x_s;x_o)$ is a cosmological observable, and as such, does not depend on the choice of FLRW coordinates we use to describe it. Geometrically, this means that $z$ is a scalar under diffeomorphisms of $\mM$.
We will make use of this fact to compute the gauge transformation of the perturbation $\de z(x_s;x_o)$, which again is not a tensor. As discussed in the previous subsection, we will suppress the dependence of $z$ and $\de z$ from the reference position $x^\mu_o$, since the evaluation at the observer position $O$ does not participate in gauge transformations.

Second-order expectation
We recall the transformation law of a scalar under change of coordinates in Eq. ($\ref{scalar}$), therefore the observed redshift satisfies
z ( x_s)=z(x_s) ,
where the reference coordinates in the left-hand side are $\widetilde x^\mu_o$, and in the right-hand side $x^\mu_o$. Using our definition for the splitting in Eq. ($\ref{dzdz}$)
we can recast the previous equation in terms of $\de z$ as

z( x)=a(η)/a(η)(1+z(x))-1 ,
where we suppressed the subscript $s$. This expression is valid at any perturbative order, and we can obtain a second-order expression expanding $\widetilde x^\mu$ in Eq. ($\ref{xcsi}$) around $x^\mu$. The expansion of the scale factor is
Substituting this result in Eq. ($\ref{dzimp}$) we obtain at second order

z(x)=z+T(1+z) +1/2(^2+'  )T^2 ,
which is the expression that we will need for the consistency check. Combining the three equations above, we derive the second-order gauge transformation of $\de z$:

z   =   z+T(1+z) -TT'-^T_,+1/2(^2-'  )T^2
-Tz'-^z_, ,
where all the quantities are evaluated at the coordinate value of the source position $x^\mu_s$, and the dependence on the observer reference position is only implicit.

Second-order direct computation
To compute the gauge transformation of our second-order expression of the distortion in the observed redshift given in Eq. ($\ref{dz2}$), we need to compute the gauge transformations of the perturbations $\de\eta_o$ and $\dhnu$. The perturbation $\de\eta_o$ is related to the time lapse at the observer position via Eq. ($\ref{lapseconf}$), and we have already checked the consistency of our second-order expression for the time lapse at the observer position. We only need to check the gauge transformation of $\dhnu$ directly from Eq. ($\ref{dhnu2}$). Note that this expression contains perturbations of the conformal wave vector and perturbations of the metric. Since we have already checked the consistency of perturbations in the conformal wave vector, we can use the results of the previous subsection, together with the results in Subsection $\ref{subsec:4.1}$, to derive the second-order gauge transformation
=      +(ℋ T -ℋ_oT_o)(1+)+ 12(ℋ^2  -  ℋ' ) T^2+ 12 (ℋ_o^2 -  ℋ_o') T_o^2- 'T

-  ℒ^ (_, +  ℋ  T_,) -  ℋ ℋ_o TT_o  -TT'  .
This expression is evaluated at $\bar x^\mu_z$ as in Eq. ($\ref{dhnu2}$) with the reference observer position $O$. For consistency, we expand the second-order fields around the background observer position $\bar x^\mu_{\rm o}$ via Eq. ($\ref{2bar}$):

Substituting the previous two equations in the expression for the redshift perturbation, together with the transformation of $\de\eta_o$ in Eq. ($\ref{detao}$), we derive the second-order
z(x)=z(x)+_zT_z(1+z) +1/2(^2+' )T^2+^μ_s[x^μ( T )]|_x̅_z ,
where the subscript $z$ indicates the evaluation of the second-order fields at $\bar x^\mu_z$.
Again, the fields in the expectation in Eq. ($\ref{dztildetilde}$) are evaluated at the source position $x^\mu_s$, while the expression above is evaluated at $\bar x^\mu_z$. Once we expand $(\HH T)(x_s)$ around $\bar x^\mu_z$ as in Eq. ($\ref{HT}$) we recover the consistency between the expectation and direct computation <cit.>.

§.§ Distortion in the source position compared to the observed position

Here we will compute the transformation of the distortion in the source position around the background $\bar x^\mu_z$ parameterized by the observed redshift and angles. The computations will be performed both in rectangular coordinates $\DX^\mu_s$ and in spherical coordinates $\drr\,,\,\dtt\,,\,\dpp$. The expression in spherical coordinates around observable quantities will be needed to compute the gauge transformations in the next subsection.

Second-order expectation
The expectation for the second-order transformation of $\DX_s^\mu$ has already been derived in Eq. ($\ref{trasfx2}$), here we write the components explicitly:

_s   =   _s+T(x̅_z)+_sT'|_x̅_z+^_sT_,|_x̅_z ,
^_s   =   ^_s+^(x̅_z)+_s'^|_x̅_z+^_s^_,|_x̅_z .
With the transformation of the spatial distortion $\DX_s^\al$, we can derive the transformation laws satisfied by the deviations in spherical coordinates. By direct substitution, we obtain
=   n_^_s+_z/2 ( ^2+sin^2 ^2 )
=   +_+^μ_s_μ_+1/2_z(^2_+^2_)+_ +_sin  ,
_z   =   _^_s- +_z/2cossin ^2
=   _z+_+^μ_s_μ_-/_z_-_(+_/_z)+1/2_z^2_
+_cos  ,
_zsin    =   _^_s-sin  -_zcos     =   _zsin +_+^μ_s_μ_-/_z_-_/_z( _zsin +_)
-__z( _zcos -_)
where we defined $\LL_\para:=\LL^\al n_\al\,$, $\LL_\ttt:=\LL^\al \ttt_\al\,$, $\LL_\pp:=\LL^\al \pp_\al\,$ and expressed each field in the right-hand side at the reference position $\bar x^\mu_z$.
Note that, once we have verified the consistency of our second-order expressions for $\DT_s$ and $\DX_s^\al$ in Eqs. ($\ref{deta2}$) and ($\ref{dx2}$) with the expectations in Eq. ($\ref{gtdx2}$), then the consistency of our second-order expressions for the distortion in spherical coordinates in Eq. ($\ref{dangoli2}$) is automatically recovered because $\drr,\dtt,\dpp$ are simply reparameterizations of $\DX^\al_s$.

Second-order direct computation
We start by checking the consistency of $\DT_s$ given in Eq. ($\ref{deta2}$). Our second-order expression depends only on $\de z$, of which we verified the consistency in the previous subsection. We can therefore substitute Eq. ($\ref{dztildetilde}$) for the gauge transformation of $\de z$. After straightforward manipulations, we recover the expected result
_s   =   _s+T(x̅_z)+_sT'|_x̅_z+^_sT_,|_x̅_z .
We now turn our attention to the spatial components $\DX_s^\al$ of the distortion in the source position. From Eq. ($\ref{Dx2}$), we need to compute the gauge transformations of boundary terms and line of sight integration. To avoid cumbersome equations, we present the terms separately. For the boundary terms, we derive
x^_o   =   x^_o+^_+x^μ_o(_μ^)|_ ,

-n^(_s-η_o)   =   -n^(_s-η_o+ T_z-T_+^μ_s(_μT)|_z
-x^μ_o(_μT)|_) ,

^μ_μ( ^+n^)   =   ^μ_μ( ^+n^+^+n^T)
+T( '^+n^'
+'^+n^T' )+^( ^_,+n^_,+^_,+n^T_, ) .
To compute the transformations of the distortion in the observer and source position, $\de x^\mu_o$ and $\DT_s$, we used Eq. ($\ref{trasfx2}$) as we already verified the consistency of these expressions. From the results of previous subsections, we can readily see that also the transformation of the linear-order perturbation $\DX_s^\al$ in Eq. ($\ref{Dx1}$) satisfies the linear-order expectation in Eq. ($\ref{trasfx2}$).
The gauge transformation of the integral yields
∫_0^_zd[ ^-n^+('^+n^')-^(^_,+n^_,)

∫_0^_zd{ ^-n^+('^+n^')-^(^_,+n^_,)

-T( '^+n^'+'^+n^T' )-^( ^_,+n^_,+^_,+n^T_, )]}|_x̅_ ,  where for the gauge transformations of the components of the conformal wave vector we used Eqs. ($\ref{dnuP}$) and ($\ref{deaP}$) since we have already checked the consistency of our second-order expressions.
Summing the contributions of the boundary terms and the integral, most of the terms cancel, and we recover
^_s   =   ^_s+^(x̅_z)+_s'^|_x̅_z+^_s^_,|_x̅_z ,
in agreement with the expectation.

§.§ Distortion in the physical volume and area

Second-order expectation
In Subsection $\ref{subsec:vol}$ we showed that the fluctuations $\de V$ and $\de A$ in Eqs. ($\ref{deV}$) and ($\ref{deA}$) are invariant under diffeomorphisms of $\mM$, which is consistent with the requirement that they are observed quantities. Consequently, the following relations hold under a change of coordinates

V(x_s)=V(x_s) ,A(x_s)=A(x_s) ,
i.e. invariant under diffeomorphism <cit.>.
These are valid at any perturbative order, while the second-order gauge transformations can be read from Eq. ($\ref{scalar}$):
V(x̅_z)=V(x̅_z)-TV'|_x̅ _z-^V_,|_x̅ _z
,  A(x̅_z)=A(x̅_z)-TA'|_x̅ _z-^A_,|_x̅ _z ,

Second-order direct computation
Given the intricacy of our second-order expressions for volume and area fluctuations in Eqs. ($\ref{DV2}$) and ($\ref{DA2}$), the direct computation of the gauge transformation of the perturbations $\de V$ and $\de A$ is a very laborious task. Beyond having to make many substitutions, the main complication is that we need to relate partial derivatives with respect to spacetime coordinates $x^\mu$ to derivatives with respect to the observed quantities $z,\ttt,\pp$. Instead of explicitly showing the consistency between the expectation and the direct computation, which would involve many equations, we will provide the critical relations used to prove the expected consistency so that the interested reader can reproduce the calculations.

* Relations between partial derivatives:
/z=-1/H_zη+1/H_zn^x^≡1/H_z d/d ,    /=_z^/x^ ,    /=_zsin ^/x^ .

* Identities with unit vectors
(+	/)^+1/sin/^   =    -2n^ ,
/^= ^cos ,    /^= -n^ ,
n^sin+^cos   =   -/^ .

* Relations involving partial derivatives with respect to Cartesian coordinates and spherical coordinates
^   =   _n^+_^+_^ ,

^_,=_,+_,^+_,^   =   _,+1/_z^_
+1/_zsin^_ ,

1/2^_,^_,-1/2^_,^_,   =   _,_,^+_,_,^+_,^_,^-_,_,^
-_,_,^-_,^_,^ ,
which also apply to $\VV^\al$, or any other field.

* Relations involving derivatives of spatial deviations
^   =   n^+^_z +^_zsin  ,

z^   =   n^z+1/H_z^+^_zz+1/H_z^sin +^_zsinz ,

^   =   n^+^+_z^+^_z+^_zcos +^_zsin ,

^   =   n^+^sin +_z^+^_z+_zsin ^
+^_zsin .

Employing these relations, and directly substituting in Eqs. ($\ref{DV2}$) and ($\ref{DA2}$) the gauge-transformation laws derived previously in the section we obtain
V(x_s)=V(x̅_z)+^μ_s(x^μV)|_x̅_z ,
A(x_s)=A(x̅_z)+^μ_s(x^μA)|_x̅_z ,  in agreement with the expectation in Eq. ($\ref{legge}$) once the right-hand side is expanded up to second order around $\bar x^\mu_z$.

§ COMPARISON WITH PREVIOUS WORKS

In this section, we compare our main results with previous work in the literature. The physical volume and area, as well as the observed redshift and source position, have already been computed at the second order in perturbations, and many applications have been proposed <cit.>.
However, in none of these previous works have the gauge transformation properties of the main results been explicitly verified. Given the complex nature of the second-order relativistic expressions of the cosmological observables, complete verification of the gauge-transformation properties in our work provides strong support for the validity of our results that are likely to have less (or no) flaws. Due to the different notations and intricacy of the
expressions, it is not straightforward to make a detailed comparison between our results and those of other groups. Hence, we can only conduct a qualitative comparison at this stage, and any quantitative comparison is reserved for future works.
Therefore, we will list previous efforts in literature, highlighting only the main differences with our results.

* Yoo, J., & Zaldarriaga, M. (2014). Beyond the linear-order relativistic effect in galaxy clustering: second-order gauge-invariant formalism <cit.>.

This article can be regarded as a prior work of what we have set out here. In fact, except for a few minor details and notation conventions, the approach and the methodology coincide but here we focus more on the exact nonlinear relations and definitions, separating perturbative expressions. In <cit.>, the general relativistic description of galaxy clustering was extended to the second order in metric perturbations without assuming any gauge conditions or adopting any restrictions on vector and tensor perturbations. Moreover, second-order gauge-invariant formalism was discussed, computing second-order gauge-invariant variables.
However, no explicit check of gauge transformation was made, over which we improve here. Furthermore, with the recent development of the tetrad formalism <cit.>, we have clarified the relation of the physical observables to coordinate-dependent quantities.

* Di Dio, E., Durrer, R., Marozzi, G., & Montanari, F. (2014). Galaxy number counts to second order and their bispectrum <cit.>.

The authors present the galaxy number counts to the second order in perturbations
and compute the angular bispectrum. They identify the dominant contributions
to the angular bispectrum of the galaxy number counts,
which originate from the density fluctuation, the redshift space distortion, and
gravitational lensing. Using the geodesic light-cone coordinate <cit.>, the authors
derive the nonlinear expressions for the cosmological observables without solving
for the Einstein equation (see, e.g., <cit.>). Then they transform their expressions to
the second-order Poisson gauge, neglecting vector and tensor contributions
to the metric perturbations but without assuming vanishing anisotropic stress.
Their expression in Eq. (3.23) shows the relation of perturbation variables
to the observed redshift:
1+z_s=a(η_o)/a(η_s)[1+z^(1)+z^(2) ] ,
which is different from ours given in Eq. (<ref>), but we stress that the
difference arises only from a different choice of splitting the background and the
perturbation. However, their explicit expression in Eq. (3.23) in terms of metric
perturbations contains only the contributions evaluated at the source position and
along the-line-of-sight direction, while no contributions at the observer position
are present. This equation can be compared to Eq. (<ref>) in our work, where
there exist not only the contributions at the observer position, but also
the quadratic contributions that involve one at the observer position and the
other at the source position. Similar differences are present in their expressions
for the angular position in Eq. (3.15), which can be compared to Eq. (<ref>) in our
work. These differences persist in their subsequent expressions based on the source
position such as the fluctuation in the physical volume and the area.

* D. Bertacca, R. Maartens, and C. Clarkson. Observed galaxy number counts on the light cone up
to second order: II. Derivation <cit.>.

This paper aims to develop the second-order framework for studying
the relativistic effects beyond the linear order, in the same spirit
of our current work. The final results are presented both in the Poisson gauge and
in a general metric representation without assuming a gauge condition,
but the vector and tensor perturbations at the linear order are ignored.
Due to the different notation convention, it is not easy to conduct
a detailed comparison, but the calculational steps are quite similar to ours.
In <cit.> (BMC hence forth) they use the comoving distance as an affine parameter $\chi$, which corresponds
to our affine parameter $-\lambda$ with negative sign and gives rise to
the sign difference in the conformal wave vector. The perturbative deviation of their affine parameter $\chi$ around the background comoving distance $\bar\chi$ is equivalent to ours in Eq. ($\ref{xz}$).
The main difference lies in their notation convention for the conformal wave vector, which is given in their Eq. (11):
k^μ(χ̅)=dx^μ/dχ̅(χ̅)=( -1+, n^i+^i)^μ(χ̅) .
The corresponding quantity $\hat k^\mu$ in our notation is defined in Eq. ($\ref{khatnohat}$) and it their notation is
k^μ(χ)=dx^μ/dχ(χ)=( -1-, n^i+^i)^μ(χ) .
Therefore, what we call $\dnu$ and $\dea^i$ in our papers are different quantities, where the difference amounts to another choice of the parameterization of the light path.
For instance, the difference in what we call $\dnu$ is given by
_MY=-_BMC-( 1-_BMC)d/dχχ .
Notice also that their geodesic equation for the conformal light wave vector above Eq. (86) looks different from what we use to derive Eqs ($\ref{dnuint}$) and ($\ref{deaint}$), namely
d/dk̂^μ+^μ_νρ|_x̅_k̂^νk̂^ρ+x^_^μ_νρ|_x̅_k̂^νk̂^ρ=0 .
However, further expanding this equation around $\bar x(\cc_z)$ and $\cc_z$, and using the linear-order relation
Δx^μ=x^μ+_z k̂̅̂^μ ,
we find consistency with their expression.
The use of $\bar\chi$ instead of $\chi$ yields their definition of $\de x^{0}$ and $\de x^i$ in Eqs. (16) and (17), which amounts to
adopting the Born approximation. However, the authors do include post-Born terms (see for instance their Eq. (113)) which make
a difference at the second order in which we have to account for the evolution
of perturbations on a perturbed light path. The observed redshift is defined in their Eq. (27), and the normalization condition for the conformal wave vector in their Eq. (28) corresponds to our
condition in Eq. ($\ref{NC}$). However, adopting $a(\eta_o)=1$, they ignored
the contributions at the observer position. Indeed, as apparent in our Eq. ($\ref{lapseandshift}$)
the observer position in an inhomogeneous universe is not uniquely set, but stochastically
set in a given realization of the inhomogeneous universe. With such normalization, all
the background equations with $a(\bar\eta_{\rm o})$, which are mostly not shown explicitly,
result in extra perturbative quantities, because of the normalization condition.
Similarly, the spatial shift of the observer position $\delta x^\alpha_o$ in our Eq. ($\ref{dxlin}$)
is ignored, which matters starting at the second order in perturbations.
Finally, the antisymmetric part of the tetrad vectors in our Eq. ($\ref{asim}$) is not
considered in BMC, which again matters in the second-order expressions for the
cosmological observations.
No verification of the gauge-transformation properties of their expressions are made.

* O. Umeh, C. Clarkson, and Roy Maartens. Nonlinear relativistic corrections to cosmological
distances, redshift and gravitational lensing magnification: II. Derivation <cit.>.

This paper studies light propagation in a generic spacetime to derive a nonlinear expression for the area distance. This is computed up to second order in perturbations in the Poisson gauge, and their results are expressed using the observed redshift as parameterization.
It is not straightforward to directly compare our results for the area distance with those obtained by the authors. These expressions are lengthy, the notation conventions are different, and the authors fix a gauge. For these reasons, we will compare the calculations in the intermediate steps rather than the final but lengthy result (see also <cit.> for similar
derivations based on the geodesic light cone gauge that include the contributions
at the observer position).
The authors consider two spacetimes: quantities defined on the physical spacetime are denoted by a hat, while quantities defined on a conformally related manifold have no hat. They specify the conformal transformation in their Eq. (21), where the only difference with our derivation is that they fix the constant $\NC=1$. A consequence of this choice is that quantities like the light wave vector or the perturbation in the observed redshift are different from our expressions. However, this should not impact the final expressions for the cosmological observables. Another difference between the two approaches is that the authors fix $\de k^\mu_o=0$ in the light wave vector  at the observer position, both at first and second order. In our case, such quantity satisfies the boundary conditions in Eq. ($\ref{ko}$) involving tetrad vectors at the observer. Such contributions at the observer position appear in our final result for the area distance but are neglected in the paper we are discussing. Moreover, in their definition of the redshift perturbation in Eq. (69), they adopt $a(\eta_o)=1$, ignoring contributions at the observer position.
Notice that the authors obtain their expression for the angular diameter distance by solving the differential equation
d^2^nD_A/d^2=^nS() ,
with boundary conditions
^nD_A(_o)=0 ,d^nD_A/d|_o=-^nE_o ,
where $\de^n S$ contains the area distance up to order $n-1$ and, in our notation, $E_0=-(\hat u^\mu\hat k_\mu)_o=1+\dhnu_o$. Given our choice of the constant $\NC$ we fixed $\dhnu_o=0$ while for the authors $E_o$ contain also a perturbative part. As a consequence, the boundary condition in Eq. (11) are different in the two approaches, hence also the solution for the angular diameter distance will be different.
There is another difference present for the relation between the affine parameter and the observed redshift. In their Eq. (126)
=ν++1/2^2 ,
where the definition of $\cc$ coincides with our $\cc_s$, while $\nu$ coincides with our $\bar\eta_z$. In Eqs. (132) and (133) they derive their expressions for $\de\cc$ and $\de^2\cc$. Such expressions coincide with what we derived for $\DT_s$ in Eq. ($\ref{deta2}$). So in our notation we could rewrite Eq. (12) as
_s=η̅_z+_s ,
which actually differs from our expression in Eq. ($\ref{xz}$) of the terms $-\bar\eta_{\rm o}-\de\eta_s$ and would then impact all the expansions in the mapping from affine parameter to observed redshift.

§ SUMMARY AND DISCUSSION

In this paper we have developed the nonlinear gauge-invariant formalism for describing the cosmological observables and presented complete second-order expressions associated with light propagation and observations without choosing a gauge condition. Our second-order expressions in general-metric representations allow us to explicitly check the validity by directly computing their gauge transformations and by comparing to the expectations. Our main results are summarized as follows:

* Adopting a geometric approach to describing the light propagation and observations, we have developed the nonlinear gauge-invariant formalism that is fundamental to computing all the general relativistic effects in cosmological observables. In particular, we have clarified the separation between the background and the perturbation quantities in a fully nonlinear way. Drawing on the previous work, we have introduced
tetrad vectors at the observer and the source positions to describe the observations recorded in the rest frame of the observer and to define the physical quantities in the source rest frame that are measured by the observer. The observable quantities are the redshift and the angular position of the source and their spectroscopic information; the physical quantities are the physical size and shape of the source and their spectrum in the rest
frame. These quantities are connected by light propagation and observations, and expressed in terms of intermediate quantities such as the coordinate position of the source or the distortion in the observed redshift. We have presented fully nonlinear expressions for the physical volume and area occupied by the source that can be measured in galaxy clustering and the luminosity distance. For future work we have left a nonlinear computation of the distortion of the source shape that can be measured in weak gravitational lensing.

* By expanding the nonlinear expressions, we have computed the second-order expressions for the cosmological observables and their intermediate quantities in terms of metric perturbations and the velocities of the source and the observer. Any nonlinear quantities can be perturbatively expanded and truncated in many different ways, though the resulting expressions should be all consistent, if the perturbation order is defined in the same way. Paying particular attention to these subtleties, we have expanded the second-order quantities along the
background photon path and around the coordinate-independent reference points of the observer and the source parameterized by the observables. The reason behind this expansion is that first we can keep perturbative orders under control by expanding each quantity around background paths. More importantly, by having employed coordinate-independent parameterization for physical positions, the entire coordinate-dependence is contained in the deviations $\de x^\mu_o$ and $\DX^\mu_s$. The main results are the second-order gauge-invariant expressions of the physical volume and area for galaxy clustering and the luminosity distance.