Perturbed redshifts from N-body simulations
Julian Adamek





julian.adamek@obspm.fr
Laboratoire Univers et Théories, Observatoire de Paris – PSL Research University – CNRS – Université Paris Diderot – Sorbonne Paris Cité,
5 place Jules Janssen, 92195 Meudon CEDEX, France



In order to keep pace with the increasing data quality of astronomical surveys the observed source redshift has to be modeled
beyond the well-known Doppler contribution. In this letter I want to examine the gauge issue that is often glossed over when
one assigns a perturbed redshift to simulated data generated with a Newtonian N-body code.
A careful analysis reveals the presence of a correction term that has so far been neglected. It is roughly proportional to the observed
length scale divided by the Hubble scale and therefore suppressed inside the horizon. However, on gigaparsec scales it can be comparable
to the gravitational redshift and hence amounts to an important relativistic effect.




Perturbed redshifts from N-body simulations
    Julian Adamek
    December 30, 2023
===========================================





§ INTRODUCTION



While the standard analysis of redshift space perturbations <cit.> only contains the leading Doppler term, recent and future
galaxy surveys have or will have sufficient statistical power to detect subleading terms such as gravitational redshift, weak lensing,
transverse Doppler shift, time delays etc. These include general relativistic effects
that are interesting targets for testing gravity at cosmological scales. For instance, a first detection of gravitational redshift in
galaxy clusters has been claimed in <cit.>. Perturbation theory can be used to understand these effects analytically
(e.g. <cit.>), and many results have been derived in longitudinal gauge. However, their application to N-body simulations
requires due care with respect to how the data are mapped to the relativistic spacetime.

One possibility is to run relativistic simulations directly in longitudinal gauge <cit.>, but the use of Newtonian N-body
codes is by far the more common practice. Coming from the Newtonian world it is not immediately clear how their results
should be interpreted in a relativistic context (e.g. <cit.>). In <cit.> it
was finally realized that one can specify a gauge, the so-called “N-body gauge,”[In fact, <cit.> had already
discovered “half of” the necessary coordinate transformation, leading to what they called “Newtonian matter gauge.” The spatial part
of the transformation is however different from the N-body gauge, leading to a non-vanishing volume perturbation.] in which the
relativistic equations are formally identical to the Newtonian ones, therefore providing a framework for a fully relativistic
interpretation of Newtonian simulations. The same authors went on to analyze initial conditions in this framework, showing that the usual
recipes are consistent with general relativity as well <cit.>.

The issue has so far mostly been discussed in relation to the matter power spectrum, and the implications for the perturbed redshift and
some other observables have therefore not yet been fully appreciated. In the recent literature the redshift formula of longitudinal gauge
is often directly applied to N-body simulation data, thereby disregarding the fact that
these are not provided in the appropriate coordinate system. For instance, <cit.> do not account for this aspect,
but I will show that the correction is fortunately very small on the scales they are interested in. While the gauge issue is noted in
<cit.> where a correction term is applied to the density perturbation (see their equation 17), the redshift of
individual sources is still computed without such a correction (see their equations 9 and 10). In the following I wish to clarify
this issue by studying the coordinate transformations involved, and finally by deriving the correct formulas for the perturbed redshift in
N-body gauge. I show that a small correction term due to the spatial coordinate transformation appears and should in principle be included
in the analysis. However, this term only becomes relevant on extremely large scales.



§ WEAK-FIELD METRIC



Astrophysical objects with high compactness exist on scales ≲ 0.01 parsec (the largest known supermassive black holes),
while on extremely large scales ≳ 100 megaparsec the Universe can be described entirely in terms of linear equations. In between
these two extremes there lies a vast range of scales where the distribution of matter can be very inhomogeneous but the gravitational fields
are weak, and the geometry is therefore only weakly perturbed. This empirical fact is confirmed every time we point a telescope at the
sky and observe that most rays seem to propagate along almost straight paths.

It is therefore possible to describe the dynamics on those scales in terms of the nonlinear evolution of matter on a
geometry with linearly perturbed metric. Of course, as for any geometry, there also exist many coordinate
systems in which the metric takes a wildly nonlinear form. A particularly well-known coordinate system in which the smallness of the
geometric perturbation is carried into effect is the one of the longitudinal gauge. For the purpose of a more general
discussion I adopt the notation introduced in <cit.> and write the generic line element for a metric with linear scalar
perturbations on top of a Friedmann-Lemaître model as

    ds^2 = a^2(τ) [-(1 + 2 A) dτ^2 + (1 + 2 H_L) δ_ij dx^i dx^j
    - 2 ∇_i B dx^i dτ - 2 (∇_i ∇_j - 1/3δ_ijΔ) H_T dx^i dx^j] ,

where a(τ) is the scale factor, τ and x^i are conformal time and comoving coordinates on the spacelike hypersurface,
respectively, and A, B, H_L, H_T are scalar functions describing perturbations.
Vector and tensor perturbations shall not be discussed here. A setting is said to be “weak-field” if there exists a coordinate system
for which all the above perturbation variables are small, i.e. |A|, |H_L|, √(∇_i B ∇^i B), √(∇_i ∇_j H_T ∇^i ∇^j H_T - (Δ H_T)^2/3)≪ 1. One can then find other such coordinate systems by making a small change of coordinates,
generated by two scalar fields T and L, so that τ→τ + T and x^i → x^i + ∇^i L. This shows that
the four scalar perturbations only contain two physical modes. The longitudinal gauge is obtained by choosing T
and L such that B = H_T = 0, leaving only the lapse perturbation A = Ψ and the volume perturbation H_L = Φ, where the two
potentials Ψ and Φ denote precisely the two physical modes when written as first-order gauge-invariant expressions. The typical
amplitude of these perturbations is ∼ 10^-5 in our Universe, except for the vicinity of black holes or neutron stars where the
weak-field description breaks down.

In longitudinal gauge the Hamiltonian constraint reads

    -ΔΦ + 3 ℋΦ' - 3 ℋ^2 Ψ = 4 π G a^2 δρ ,

where ℋ is the conformal Hubble rate, a prime denotes partial derivative with respect to τ, and the equation is
linearized in Φ and Ψ but not in δρ. In fact, the matter perturbation δρ can be very large, but it is
computed and evolved on a linearly perturbed geometry. This equation is however not the one used in a Newtonian N-body code. First, such
a code uses counting densities that are not corrected for perturbations of the volume element – Newtonian theory assumes Euclidean
geometry. Second, the Poisson equation lacks some of the terms featured on the left-hand side.

From now on I restrict the discussion to the case where gravitational fields are sourced exclusively by nonrelativistic matter.
In this case anisotropic stress can be neglected, implying Φ = -Ψ.



§ PERTURBED REDSHIFT IN LONGITUDINAL GAUGE



The effect of geometry (and perturbations thereof) on observables can be understood by studying the geodesics of photons that reach
the observation event. Let me denote the tangent vector of such a geodesic as k^μ.
The condition g_μν k^μ k^ν = 0 implies that k^i/k^0 = n^i (1 + 2 Ψ), where n^i is the unit
vector (δ_ij n^i n^j = 1) pointing in the direction the photon is traveling. The geodesic equation can then be expressed
as an evolution equation for the energy,

    dln k^0/dτ = - 2 ℋ - 2 n^i ∇_i Ψ ,

and an equation describing the deflection of the ray,

    dn^i/dτ = -2 (δ^ij - n^i n^j) ∇_j Ψ .

A reference clock is specified through a unit timelike vector
u^μ = a^-1 (1-Ψ, v^i) where v^i = dx^i / dτ is the peculiar (coordinate) velocity of the clock's rest frame.
In such a frame, the measured photon energy is -g_μν u^μ k^ν. A first integral of
eq. (<ref>) yields the following expression for the observed redshift between source (src) and observer (obs):



    1+z = g_μν u^μ k^ν|_src/g_μν u^μ k^ν|_obs = 
    a_obs/a_src(1 + n_i v^i_obs - n_i v^i_src + Ψ_obs - Ψ_src - 2∫_src^obsΨ' dχ)

Here dχ is a conformal distance element along the photon path. The redshift perturbations are easily identified as the Doppler
shift due to peculiar motion, the gravitational redshift due to time dilation, and the Rees-Sciama effect (also known as integrated
Sachs-Wolfe effect in the context of linear theory).

The boundary terms in this expression have to be evaluated at the coordinate time at which the photon geodesic actually intersects the
world lines of source and observer. The time of flight of the photon is affected by the Shapiro delay, which for a coordinate distance
χ between source and observer is given by the following relation:

    χ = ∫_src^obs dχ = ∫_src^obs(1 + 2 Ψ) dτ = τ_obs - τ_src + 2 ∫_src^obsΨ dχ





§ FROM LONGITUDINAL TO N-BODY GAUGE



Given a linearly perturbed geometry specified in longitudinal gauge I now make a small change of coordinates to set
H_L = 0. Considering the Lie derivative of the metric tensor the required transformation is generated by T, L
that satisfy

    ℋ T + 1/3Δ L = Ψ .

I furthermore choose L' = 0 such that velocities are not transformed.
In this new coordinate system one finds B = T and H_T = -L. It is important to verify that these new perturbations do not
violate the weak-field conditions[For instance, had I made a change of coordinates that sets A = B = 0 instead, like in
a synchronous gauge, I would have found L' = T = -a^-1∫^a (Ψ/ℋ)dã and that H_L receives a contribution
∼ΔΨ/ℋ^2 ∼δρ/ρ. The volume perturbation therefore does not remain small everywhere in such a
coordinate system, rendering the weak-field treatment inconsistent.]. As explained below, L shall be chosen such that Δ L
is of order Ψ, and hence the same is true for the perturbations generated by H_T. One can then easily convince oneself that
∇_i T is of the order of a peculiar velocity — in fact, it coincides with the Zel'dovich approximation thereof — and therefore
√(∇_i B ∇^i B)∼ v ≲ 10^-3 at low redshift.
The shift perturbation in the new coordinate system is therefore substantially larger than Ψ, but still comfortably
within the weak-field regime.
The lapse perturbation becomes

    A = Ψ + ℋ T + T' ,

and as explained shortly I arrange that it vanishes at leading order.

The Hamiltonian constraint becomes

    Δ(ℋB - 1/3Δ H_T) + 3 ℋ^2 A = 4 π G a^2 δρ ,

and as a consequence of H_L = 0 the density perturbation δρ can be obtained simply by counting the mass elements per
coordinate volume, in accordance with the procedure relevant for Newtonian codes.
For nonrelativistic particles the geodesic equation
reduces to

    dv^i/dτ + ℋ v^i = ∇^i(ℋ B + B' - A) .

Since I assume that matter is nonrelativistic and hence
the pressure perturbation can be neglected, the spatial trace of Einstein's equations yields

    ℋ A' - (ℋ^2 - 2a”/a) A = 0   ⇒   A ∝1/aℋ^2 .

This shows that an appropriate choice of boundary conditions will set A = 0.
Furthermore, with such a choice eqs. (<ref>), (<ref>) are formally identical to the ones of Newtonian
gravity if eqs. (<ref>), (<ref>) are used,

    ΔΨ = 4 π G a^2 δρ ,    dv^i/dτ + ℋ v^i = -∇^i Ψ .


Keeping in mind that the above equations remain valid even if δρ/ρ becomes large I now want to set the boundary conditions
at early times when matter perturbations are still linear. According to eq. (<ref>) the condition A = 0 is satisfied when

    T = -1/a∫^a Ψ/ℋ dã .

The linear solution of Ψ is constant in matter domination, and with
eq. (<ref>) the corresponding choice of L is given by Δ L = 5 Ψ_in. Here I introduce Ψ_in
to denote the linear initial condition for Ψ in matter domination.

With this choice one can see that ∇^i B = v^i in the linear regime, which (together with H_L = 0) is the original gauge condition
used in <cit.> for the N-body gauge. So even though my gauge condition H_T' = -L' = 0 is different, the resulting
coordinate system is the same whenever only nonrelativistic matter is present. The advantage of my condition is that it does not
explicitly refer to a matter perturbation and can hence be easily extended into the nonlinear regime of matter as long as gravitational
fields remain weak. In a slight abuse of terminology I shall therefore always call this system of coordinates the one of N-body gauge, as
the state of a Newtonian N-body simulation is given in precisely these coordinates.

I shall now discuss the repercussions of this change of coordinates from longitudinal to N-body gauge. The null condition is solved in
N-body gauge by k^i/k^0 = n^i + ∇^i B + n^j (∇_j ∇^i - δ^i_j Δ/3)H_T, and the photon geodesic
equation can be written as

    dln k^0/dτ = -2ℋ - n^i n^j ∇_i ∇_j B ,

and

    dn^i/dτ = -(δ^ij - n^i n^j)∇_j (n^k ∇_k B - 1/3Δ H_T) .

Considering how the coordinate transformation acts on u^μ one sees that u^μ = a^-1 (1, v^i) in the new coordinates. Thus,
a first integral of eq. (<ref>) gives the following new expression for the observed redshift:

    z+1 = a_obs/a_src(1 + n_i v^i_obs - n_i v^i_src + Ψ_obs - Ψ_src
    
     +ℋB|_obs - ℋB|_src - 2∫_src^obsΨ' dχ)

In order to recover all the terms of eq. (<ref>) I used (ℋ B)' = Ψ' and ℋB + B' = -Ψ, but
evidently a new boundary term ℋB appears. Noting that B = T this boundary term can be understood as the
result of the change of coordinates acting on the background term a_obs/a_src. In other words, the term has to appear
because the equal-time hypersurfaces in N-body gauge do not coincide with the ones of longitudinal gauge. The coordinate time
of a Newtonian N-body simulation is the one of N-body gauge, and hence this boundary term needs to be taken into account.

Let me now inspect the time of flight for the photon in N-body gauge,

    χ = ∫_src^obs(1 + n^i ∇_i B + n^i n^j ∇_i ∇_j H_T - 1/3Δ H_T)dτ
    
     = τ_obs - τ_src + B_obs - B_src + n^i ∇_i H_T |_obs
    - n^i ∇_i H_T |_src + 2 ∫_src^obsΨ dχ ,

where I again use the gauge conditions to recover the terms known from longitudinal gauge. Compared to eq. (<ref>) there are
two new boundary terms. These are expected from the gauge transformation, since the coordinate time transforms as τ→τ + T = τ + B, and the coordinate distance transforms as χ→χ + n^i ∇_i L|_obs - n^i ∇_i L|_src = χ - n^i ∇_i H_T|_obs + n^i ∇_i H_T|_src.
For a photon trajectory with fixed endpoints, the coordinate time of emission and observation therefore transforms such that the
boundary terms due to the shift perturbation in eqs. (<ref>) and (<ref>) cancel exactly.
The other boundary term in eq. (<ref>) gives precisely the change in the coordinate distance due to the spatial transformation
between longitudinal and N-body gauge. Therefore, the perturbed redshift for the trajectory remains invariant.



§ DISCUSSION



The precedent analysis clarifies that in order to use the longitudinal gauge for computing the perturbed redshifts with N-body simulation
data one should transform the coordinates appropriately. In particular, the coordinate distance between observer and sources changes
according to a spatial transformation that is independent of time. This has already been pointed out in <cit.>, and I
explicitly show how to recover this result in the relativistic framework provided by the N-body gauge. Alternatively the computation of
the perturbed redshift can also be carried out directly in N-body gauge. In this case the coordinates of sources are directly taken from
the simulation, and the effect appears as a modification of the Shapiro delay.

In order to estimate the amplitude of the correction, let me compute the typical change δχ = n^i ∇_i L|_obs - n^i ∇_i L|_src of the coordinate distance. Using the relation Δ L = 5 Ψ_in the variance of δχ
is given by

    ⟨δχ^2⟩ = 50 ∫_0^∞dk/k^3[1/3 - 1/k χj_1(kχ) + j_2(kχ)]Δ_k^Ψ ,

where Δ_k^Ψ is the dimensionless power spectrum of Ψ_in. Unfortunately the integral has an infrared divergence
for nearly scale invariant spectra, but in practice this divergence is regulated by the finite size of a simulation. Imposing a cutoff
close to the Hubble scale one finds that δχ / χ∼ 10^-4 almost independent of scale χ and precise value of the
cutoff. Considering how a change in the coordinate distance affects the time of flight one sees that the typical correction to
the redshift due to this coordinate effect is δ z / (1 + z) ∼ℋδχ. Therefore the effect becomes of the order
of the gravitational redshift for trajectories at or above the gigaparsec. At these extreme scales the gauge correction is of the same
order as all other relevant terms and should be taken into account.

As suggested in <cit.> the gravitational redshift may be measured statistically by looking for “excess redshift” of
the brightest galaxies at the center of clusters when compared to the fainter galaxies in the outskirts. For simulating such a measurement
the relevant correction is given by the difference in time of flight between the different sources, and therefore the effect is suppressed
by the small ratio between the scale of the cluster and the Hubble scale. One expects that in this case the correction is typically less
than 1% of the signal and can therefore safely be neglected.






I thank R Durrer, C Rampf and Y Rasera for comments on the manuscript as well as C Fidler and T Tram for many insightful discussions about 
the N-body gauge. I further enjoyed a valuable correspondence with D Bertacca and C Porciani about <cit.>. Prior to
submission, C Fidler et al. also kindly shared their manuscript <cit.> on a closely related topic.


utcaps


