J
C
A
P03(2023)019
ournal of Cosmology and Astroparticle Physics
An IOP and SISSA journalJ
Cosmological distances with
general-relativistic ray tracing:
framework and comparison to
cosmographic predictions
Hayley J. Macpherson
Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Cambridge CB3 0WA, U.K.
E-mail: hayleyjmacpherson@gmail.com
Received October 3, 2022
Revised December 19, 2022
Accepted January 19, 2023
Published March 7, 2023
Abstract. In this work we present the first results from a new ray-tracing tool to calculate
cosmological distances in the context of fully nonlinear general relativity. We use this
tool to study the ability of the general cosmographic representation of luminosity distance,
as truncated at third order in redshift, to accurately capture anisotropies in the “true”
luminosity distance. We use numerical relativity simulations of cosmological large-scale
structure formation which are free from common simplifying assumptions in cosmology. We
find the general, third-order cosmography is accurate to within 1% for redshifts to z ≈ 0.034
when sampling scales strictly above 100 h−1 Mpc, which is in agreement with an earlier
prediction. We find the inclusion of small-scale structure generally spoils the ability of the
third-order cosmography to accurately reproduce the full luminosity distance for wide redshift
intervals, as might be expected. For a simulation sampling small-scale structures, we find a
∼ ±5% variance in the monopole of the ray-traced luminosity distance at z ≈ 0.02. Further,
all 25 observers we study here see a 9–20% variance in the luminosity distance across their sky
at z ≈ 0.03, which reduces to 2–5% by z ≈ 0.1. These calculations are based on simulations
and ray tracing which adopt fully nonlinear general relativity, and highlight the potential
importance of fair sky-sampling in low-redshift isotropic cosmological analysis.
Keywords: cosmological simulations, gravity, cosmic web
ArXiv ePrint: 2209.06775
c© 2023 The Author(s). Published by IOP Publishing
Ltd on behalf of Sissa Medialab. Original content from
this work may be used under the terms of the Creative Commons
Attribution 4.0 licence. Any further distribution of this work must
maintain attribution to the author(s) and the title of the work,
journal citation and DOI.
https://doi.org/10.1088/1475-7516/2023/03/019
J
C
A
P03(2023)019
Contents
1 Introduction 1
2 FLRW cosmography 2
3 General cosmography 3
3.1 Applications of cosmography 4
4 Simulations 5
5 Analysis software 6
5.1 Mescaline: Extracting interesting things from Cactus 6
5.2 Metric conventions and definitions 6
5.3 Assumptions and input 7
6 Ray tracing 8
6.1 Geodesics 8
6.1.1 Evolving the geodesic equation 9
6.2 Distances 9
6.2.1 Screen space 10
6.2.2 Jacobi matrix 11
6.2.3 Evolving the Jacobi matrix equation 12
6.2.4 Weyl tensor 13
6.3 Initial data 15
6.4 Time stepping, interpolation, and tests 16
6.5 Observers 16
7 Cosmography calculation 17
8 Results and discussion 18
8.1 Strictly large-scale simulations 18
8.2 A more realistic model universe 21
8.3 An even more realistic model universe 23
9 Caveats 28
10 Conclusions 29
A Propagation of the screen vectors 30
A.1 Observers at rest 32
A.2 Initial screen vectors 32
A.3 Test of initial screen vectors 32
B Boosting data to the hypersurface frame 35
C Spatial interpolation 36
– i –
J
C
A
P03(2023)019
D Transformation of vectors from local Fermi frame 37
E Tests of the mescaline ray tracer 38
E.1 Definition of error and convergence 38
E.2 EdS spacetime 39
E.2.1 Analytic solutions 39
E.2.2 Mescaline vs. analytic solutions 41
E.3 Linearly-perturbed EdS spacetime 42
E.3.1 Semi-analytic solutions 42
E.3.2 Mescaline vs. linearly-perturbed EdS 44
E.4 Errors in analysis of simulation data 45
E.5 Violation of the null constraint 46
E.5.1 Linearly-perturbed metric 47
E.5.2 Simulation test 48
1 Introduction
Standard cosmology commonly adopts the Copernican principle: the assumption that we
are not privileged observers in the Universe. Combined with primary cosmic microwave
background (CMB) radiation measurements, which strongly disfavour large deviations from
isotropy in the early Universe [1], this leads us to conclude that the Universe is statistically
homogeneous and isotropic on large scales. These conditions comprise the cosmological
principle which forms the backbone of the current standard Λ cold dark matter (ΛCDM)
cosmological model. This very simple standard model has enjoyed decades of success, providing
a consistent fit to most of our cosmological observations to date. However, as our cosmological
data becomes increasingly more precise we are beginning to uncover some disagreements
between theoretical predictions from ΛCDM and our observations [see 2, for a recent review].
The Universe is inhomogeneous and anisotropic on small scales. The transition to
statistical homogeneity in the galaxy distribution has been measured to occur at ∼ 70–80h−1
Mpc [e.g. 3, 4], however, some studies claim to have found cosmological structures above these
scales [e.g. 5, 6] — potentially violating the cosmological principle. Tests of isotropy using late
Universe large-scale structure probes do not offer a clear consensus — with some works yielding
consistency with the cosmological principle [e.g. 7–12] and some finding disagreement [e.g.
13–19]. See Aluri et al. [20] for an overview of observational tests and consistencies with the
cosmological principle.
Maintaining the cosmological principle, we might further assume that the transition to
statistical homogeneity and isotropy coincides with a transition to a Friedmann-Lemaître-
Robertson-Walker (FLRW) geometrical description. This coincidence is not necessarily obvious
either in the context of spatial averages [21–23] or light cone averages [24].
With the exception of several anomalies measured in CMB data (see Schwarz et al. [25],
section 6 of Planck Collaboration et al. [1], and section III of Aluri et al. [20]), we observe the
CMB radiation to be statistically isotropic. However, this condition alone is not sufficient
to guarantee a near-FLRW geometry in our vicinity. The latter also requires that the low-`
multipoles in the CMB — as well as their time and space derivatives — are small [26].
The FLRW metric assumption forms the basis of ΛCDM and is prevalent throughout
modern cosmological analysis. Particularly relevant to this work is the use of FLRW cosmog-
raphy: a theoretical approximation of the luminosity distance using a Taylor series expansion
– 1 –
J
C
A
P03(2023)019
in redshift, usually truncated at third order [see 27, and section 2]. This method allows
constraints of FLRW parameters without assuming a particular cosmological expansion history,
and is most commonly used in low-redshift supernova Type 1a (SN Ia) analysis [e.g. 28, 29].
Most state-of-the-art cosmological simulations implement ΛCDM by adopting an a
priori assumed FLRW background space-time which is unaffected by the nonlinear collapse
of structures. While these simulations have proven invaluable in developing our knowledge
of structure formation beyond perturbation theory, this assumption limits our ability to
study regimes beyond FLRW with these simulations. The weak field general-relativistic
N-body code gevolution offers an improvement by including the effects of perturbations
on the background expansion [30, 31]. A more recent development is the use of numerical
relativity (NR) in cosmological simulations [32–36], allowing the removal of the concept of a
background space-time altogether. Such simulations provide the ideal testbed for the FLRW
assumption in the late-time nonlinear Universe.
Recently, Heinesen [37] derived a generalised version of the luminosity-distance redshift
cosmography — incorporating all sources of inhomogeneity and anisotropy due to local
structures (see section 3). In Macpherson and Heinesen [38, hereafter MH21], the authors
studied the impact of local anisotropy on the generalised cosmographic parameters using
NR simulations. They found that such effects could bias an effective measurement of the
parameters if the observer did not fairly sample their sky. In appendix A of MH21, the authors
find that the general cosmography (as truncated to third order in redshift) may not converge for
large redshift intervals. In this work, we extend the analysis of MH21 using ray-tracing in fully
nonlinear GR. We will study the regime within which the general cosmography provides an
accurate description of the ray-traced luminosity distance, comparing to the estimate given in
MH21. We use the same set of observers in the very smooth simulations presented in MH21, as
well as studying the anisotropies in more realistic simulations with more small-scale structure.
Anisotropies in the luminosity distance have been studied in the context of perturbed
FLRW models [39–42, e.g.]. Scatter in the Hubble diagram of luminosity distances has been
studied using the gevolution code [43] as well as with NR simulations [44]. However, to the
best of our knowledge, an analysis of the anisotropic signatures in the luminosity distance in
simulations within the fully nonlinear context of GR has not yet been done [though see 45,
for an analysis of weak-lensing convergence in NR simulations].
In section 2 and 3 we provide a brief review of the FLRW and general cosmographies,
respectively, in section 4 we describe the NR simulations, and in section 5 we introduce our
post-processing analysis software. In section 6 we provide a detailed overview of the theory
behind the ray-tracing code, including the propagation of geodesics, distance calculations,
initial data, and some technical code details. In section 7 we briefly review the calculation
of the cosmographic parameters in MH21, and in section 8 we compare these calculations
to our ray-tracing analysis. We discuss some important caveats in section 9 and conclude in
section 10. Throughout this work, Greek indices represent space-time indices and take values
0 . . . 3, and Latin indices represent spatial indices and take values 1 . . . 3, and repeated indices
imply summation. We use geometric units with c = G = 1 throughout, where c is the speed
of light and G is the gravitational constant.
2 FLRW cosmography
For the low-redshift cosmological analysis of standarisable candles, we are able to use the
cosmographic relation between luminosity distance and redshift. Specifically, it is common to
– 2 –
J
C
A
P03(2023)019
adopt the FLRW cosmography of Visser [27], which first involves performing a Taylor series
expansion of the luminosity distance in redshift about the point of observation, namely,
dL(z) = d(1)L z + d
(2)
L z
2 + d(3)L z
3 +O(z4), (2.1)
as truncated at third order. For a strictly FLRW metric tensor, the coefficients d(i)L are
d
(1)
L,FLRW =
1
Ho
, d
(2)
L,FLRW =
1− qo
2Ho
, (2.2a)
d
(3)
L,FLRW =
−1 + 3q2o + qo − jo + Ωko
6Ho
, (2.2b)
which are expressed in terms of the Hubble (Ho), deceleration (qo), curvature (Ωko), and jerk
(jo) parameters. Here, the subscript o indicates the parameters are evaluated at the observer
location. These parameters are defined in terms of the FLRW metric components by
H ≡ a˙
a
, q ≡ − a¨
aH2
, j ≡
˙¨a
aH3
, Ωk ≡ − k
a2H2
, (2.3)
where a is the FLRW scale factor, k is the (constant) scalar curvature of spatial sections, and
an over-dot represents a derivative with respect to time.
The parameters (2.3) are derived without making any assumption on the form of the
equations specifying the evolution of the FLRW scale factor a(t), which makes it an attractive
formalism to use to infer cosmological parameters. This means, for example, we do not need
to assume a particular energy density of matter or dark energy to infer the parameters (2.3)
from low-redshift data. For this reason, such constraints are “model-independent” within the
context of the FLRW metric tensor. We should note that using a cosmographic expansion
with z as the expansion parameter is strictly only valid for redshifts z < 1, see section 3.1 for
further discussion on this.
3 General cosmography
The coefficients (2.2) are derived assuming an FLRW metric tensor. Recently, Heinesen [37]
derived these coefficients relaxing this assumption and thus arriving at a cosmographic
relation for the luminosity distance which is fully general in terms of the metric tensor and
field equations.1 Within this formalism, the cosmological parameters are replaced by their
“effective” namesakes which are dependent on both the observer’s position and the direction
of observation.
The coefficients of the expansion (2.1) in the generalised cosmography are
d
(1)
L =
1
Ho
, d
(2)
L =
1−Qo
2Ho
, (3.1a)
d
(3)
L =
−1 + 3Q2o +Qo − Jo +Ro
6Ho
, (3.1b)
1Some assumptions have still been made, including the existence of a time-like congruence of observers and
emitters as well as some regularity requirements for the Taylor series expansion, see Heinesen [37] for details.
– 3 –
J
C
A
P03(2023)019
where a subscript o again indicates that quantity is evaluated at the point of observation.
The effective Hubble, deceleration, curvature, and jerk parameters are defined as
H = − 1
E2
dE
dλ , (3.2a)
Q ≡ −1− 1
E
dH
dλ
H2
, (3.2b)
R ≡ 1 +Q− 12E2
kµkνRµν
H2
, (3.2c)
J ≡ 1
E2
d2H
dλ2
H3
− 4Q− 3 , (3.2d)
respectively. In the above, we are considering an observer moving with 4-velocity uµ observing
an incoming photon with 4-momentum kµ which has travelled along a geodesic parametrised
by λ. The energy of the photon, as seen by the observer, is E ≡ −uµkµ, the derivative along
the line of sight is ddλ = kµ∇µ, and Rµν is the Ricci tensor of the space-time. The set of
parameters (3.2) can be considered as inhomogeneous and anisotropic generalisations of the
FLRW parameters (2.3). Each of the effective cosmological parameters can be expressed as
multipole series expansions in the vector describing the direction of observation, eµ. This
might represent, for example, the direction of an object in a given cosmological survey. Such
an expansion allows us to pinpoint the anisotropic contributions to each effective cosmological
parameter. For the effective Hubble parameter, for example, we have
H(e) = 13θ − e
µaµ + eµeνσµν , (3.3)
where θ is the volume expansion of the congruence defined by uµ, aµ is its 4-acceleration,
and σµν is its shear tensor [see MH21 or 37, for explicit definitions of these quantities].
Similarly, the higher-order effective cosmological parameters can also be written as multipole
decompositions, however, we refer the reader to Heinesen [37] for their explicit forms.
3.1 Applications of cosmography
FLRW cosmography is perhaps most commonly adopted in the analysis of supernova Type Ia
(SN Ia) as standardisable candles to infer the local Hubble expansion [e.g. 28, 29]. Any Taylor
series expansion of luminosity distance using the redshift z as a parameter is divergent for
z > 1 [46]. This result holds for both the FLRW and general cosmographies outlined above
and we are therefore constrained to studying the low-redshift Universe with these methods.
Alternatives have been proposed to be able to study high-redshift data with cosmography
by, e.g., re-parameterising the redshift as y = z/(1 + z) [46] or using Padé approximations in
place of the Taylor expansion [47]. Such formalisms have been applied to study observational
samples with z > 1 such as quasars, gamma-ray bursts, baryon acoustic oscillations, and
SNe Ia [e.g. 48–50]. However, concerns have been raised regarding the convergence of these
techniques at high redshift [51].
The general cosmographic formalism outlined in section 3 contains a total of 61 indepen-
dent degrees of freedom [see 37] which is too many to constrain using current data. In Heinesen
and Macpherson [52], the authors used the ‘quiet universe’ dust models — which should
provide a good description to the late-time Universe — to bring this down to ∼ 30 degrees of
freedom. However, even this is still far too many to constrain with currently available data
sets. In MH21, the authors found the effective Hubble and deceleration parameters of the
– 4 –
J
C
A
P03(2023)019
general cosmography should be dominated by a quadrupole and dipole anisotropy, respectively.
Focusing on the anisotropic signals we might expect to be dominant, and neglecting all others,
is one way to drastically reduce the number of degrees of freedom in data analysis.
Recently, several works have constrained a dipole anisotropy in the cosmographic de-
celeration parameter from SN Ia data by adding a dipole term on top of the usual FLRW
cosmography [53–56]. Dhawan et al. [57] presented the first constraints of anisotropies in
SN Ia data which were theoretically motivated by the simplified general cosmography [MH21,
37], including the first constraints on a quadrupole in the Hubble parameter from SN Ia.
While we cannot claim a significant anisotropy in current data, new and improved data sets
will allow for much stronger constraints on individual multipoles — as well as more complete
hierarchies of multipoles — in near-future analyses [58].
In appendix A of MH21, the authors showed that the general cosmography as expanded
to third order in redshift may only accurately represent the true luminosity distance for narrow
redshift intervals. Specifically, for smoothing scales of 100h−1 Mpc, they predicted that the
cosmographic luminosity distance should be correct to within ∼ 1% for redshifts z ∼0.02–0.03.
For larger smoothing scales of 200h−1 Mpc, it was predicted to be accurate within ∼ 1% for
redshifts z ∼0.04–0.06. Determining the regime of applicability of the cosmographic luminosity
distance is vital in ensuring accurate results when applying it to data, however, the above
numbers are order-of-magnitude estimates. Methods to alleviate these potential issues could
involve allowing for a rapid decay of the anisotropies with redshift [as is considered in, e.g.
53, 55–57], or implementing some smoothing of observational data prior to the cosmological fit.
A key part of this paper is to validate the regime of validity of the cosmographic luminosity
distance, as estimated in MH21, using distances calculated with ray tracing.
4 Simulations
We are interested in studying the regime of validity of the general cosmographic framework.
Since this framework makes no assumption on the specific form of the metric tensor, we wish
to maintain this generality in our analysis as much as possible. Numerical relativity has
proven to be a viable tool for cosmological simulations without the need to make assumptions
on the existence of a fictitious global “background” metric tensor [32–34]. In particular, such
simulations have proven useful for studies of general-relativistic effects on observables in the
context of inhomogeneous cosmology [35, 38, 44, 59].
Here, we will use the same simulations as in MH21 in order to be able to directly
compare their predictions with our results. The simulations were performed using the Einstein
Toolkit2 [60, 61]; an open-source NR code based on the Cactus3 infrastructure. The ET has
been adapted and used for cosmological simulations in a number of works [e.g. 33, 34, 62, 63]
and has proven to be a valuable tool to study inhomogeneous cosmology in the nonlinear regime.
Our cosmological initial data was generated using FLRWSolver4 [34] under the assumption of
linear perturbations atop a flat FLRW background space-time. The initial density fluctuations
are assumed to be Gaussian random and are drawn from the matter power spectrum at
z = 1000 generated with the CLASS5 code [64, 65]. The initial power spectrum is generated
using h = 0.7 (and otherwise default parameters), which need only be used to translate quoted
2https://einsteintoolkit.org.
3https://www.cactuscode.org.
4https://github.com/hayleyjm/FLRWSolver_public.
5https://lesgourg.github.io/class_public/class.html.
– 5 –
J
C
A
P03(2023)019
length scales from the simulation and does not imply a particular Hubble expansion. The
simulation evolution is matter dominated (no cosmological constant is included in Einstein’s
equations) and thus we compare our results to the equivalent FLRW model, namely, the
Einstein-de Sitter (EdS) model. We adopt a continuous fluid description for the matter
content, with pressure such that P  ρ, and use periodic boundary conditions [see 66, for
more specifics on the simulation setup]. We note that once the simulation starts there is no
explicit constraint for the space-time to remain close to an FLRW background, however, we
do find that spatial averages of the model universe on the z ≈ 0 slice agree with the EdS
model to within a few percent [see MH21 and also 66].
We use two simulations with cubic domain lengths L = 12.8h−1 Gpc and 25.6h−1 Gpc,
both with numerical resolution N = 128 (where the total resolution is N3 grid cells). These
domain size and resolution combinations were chosen in order to explicitly exclude any
structure beneath 100 and 200h−1 Mpc, respectively, due to the requirements of the general
cosmography in sampling only large-scale (expanding) space-time regions [see 37, and MH21].
5 Analysis software
We use the new ray tracer built as a part of the mescaline post-processing analysis code [66]
to calculate the luminosity distance and redshift of geodesics in our simulations. In section 5.1
we present some key details of the code, in section 5.2 we introduce relevant conventions and
definitions, and discuss the assumptions made in the code in section 5.3.
5.1 Mescaline: Extracting interesting things from Cactus
Mescaline is a post-processing analysis code written specifically to study the effects of
nonlinear GR in cosmological simulations performed with the ET [66]. In Macpherson
et al. [66] and Macpherson et al. [67] the authors used mescaline to study the averaged
dynamics of inhomogeneous ET simulations on a variety of scales to address the backreaction
problem for the first time in a realistic cosmic web. As mentioned above, MH21 used
mescaline to study the anisotropies in the luminosity distance as predicted by the general
cosmography outlined in section 3. All works using mescaline thus far have considered
either spatial averaging or the study of local derivatives in the observer’s vicinity to predict
observable signatures.
In this paper, we make use of the new ray-tracing capabilities of mescaline for the first
time. This new feature of the code advances the geodesic equation and the Jacobi matrix
equation to calculate the angular diameter distance, luminosity distance, and redshift along
a geodesic in the numerical space-time (see section 6 for details of the equations evolved).
In appendix E we show that the ray tracer matches known analytic solutions and prove its
numerical convergence at the expected rate.
5.2 Metric conventions and definitions
The primary application of mescaline is the analysis of NR simulations. Therefore, we adopt
a 3+1 split of space-time, where the metric tensor gµν takes the form
ds2 = gµνdxµdxν , (5.1)
= −α2dt2 + γijdxidxj , (5.2)
where xµ = (t, xi) are the space-time coordinates of the simulation, α is the lapse function
describing the spacing between spatial surfaces in time, γµν ≡ gµν + nµnν is the induced
– 6 –
J
C
A
P03(2023)019
spatial metric of the surfaces, and we adopt a gauge condition for the shift vector of βi = 0
throughout. The normal vector describing the 3-dimensional spatial surfaces is defined as
nµ ≡ −α∇µx0, (5.3)
where ∇µ is the covariant derivative associated with gµν , and we have nµ = (−α,0) and
nµ = (1/α,0) for our chosen zero shift.
The extrinsic curvature, Kij , describes the embedding of the spatial surfaces in space-
time. It is defined as the covariant derivative of the normal vector projected onto the spatial
surfaces. We can also relate it to the time derivative of the spatial metric as
∂tγij = −2αKij , (5.4)
where ∂µ ≡ ∂/∂xµ. The Christoffel symbols associated with the metric tensor gµν are defined
as
Γµαβ =
1
2g
µν (∂αgνβ + ∂βgαν − ∂νgαβ) , (5.5)
and the spatial Christoffel symbols associated with γij are
(3)Γijk =
1
2γ
il (∂jγlk + ∂kγjl − ∂lγjk) , (5.6)
where in the special case of βi = 0 we have Γijk = (3)Γijk, and so we drop the superscript (3)
for the rest of this work. The Ricci tensor of the space-time is the contraction of the Riemann
tensor, namely Rµν ≡ Rαµαν , and R ≡ gµνRµν is its trace.
Mescaline calculates the spatial Christoffel symbols, the components of the 4-Ricci
tensor and its trace, some relevant components of the Weyl tensor Cµναβ (see section 6.2.4)
and its electric part in the fluid frame, the components of the 3-Ricci tensor (3)Rij and its trace
(3)R, and the fluid rest-frame 3-curvature scalar [see section 4 of 23]. The code also calculates
kinematic and dynamical quantities in the rest-frame of the fluid, namely, the expansion
scalar θ, the shear tensor σµν , the vorticity tensor ωµν , and the 4-acceleration aµ. Following
the averaging approach of Buchert et al. [23], mescaline also calculates fluid-intrinsic spatial
averages over a user-defined domain D and calculates the effective cosmological parameters
including the kinematical backreaction, average curvature, volume of the domain VD, and
effective scale factor, aD ≡ [VD(t)/VD(tini)]1/3, where tini is the initial simulation time.
5.3 Assumptions and input
Mescaline assumes a continuous fluid description of the matter content by taking input of
a smooth density field, ρ, and velocity field with respect to the Eulerian observer, vi, at
all points on the grid. It also takes as input the spatial metric tensor of the hypersurfaces,
γij , the extrinsic curvature, Kij , and the lapse function α. It enforces the gauge condition
βi = ∂tβi = 0 throughout, however, the lapse and its time derivative are kept general
(though the latter must be specified according to the gauge condition of the simulation).
The code assumes a uniform Cartesian grid, i.e. ∆x = ∆y = ∆z for the grid spacing in
each dimension xi = (x, y, z) when taking spatial derivatives. We adopt geometric units,
G = c = 1, throughout the code and enforce periodic boundary conditions in the calculation
of spatial derivatives.
The simulation frame is assumed to be generally separate from the frame of the fluid flow,
i.e., the hypersurface normal nµ is not equal to the fluid 4-velocity uµ. However, the observers
we define in both the cosmography and ray-tracing calculations are chosen to be co-moving
with the fluid flow by defining their 4-velocity to be that of the fluid at their location.
– 7 –
J
C
A
P03(2023)019
6 Ray tracing
In this section, we detail the ray-tracing calculation performed in mescaline. In section 6.1 we
introduce the geodesic equation and describe how it is evolved and in section 6.2 we introduce
the cosmological distances relevant in this work and describe the methods for calculating
these distances along the geodesic.
6.1 Geodesics
We consider the propagation of a light bundle, or ray, through the 3+1 space-time along null
geodesics. The photon 4-momentum is
kµ ≡ dx
µ
dλ
, (6.1)
where λ is the affine parameter of the geodesic, the total derivative is d/dλ ≡ kµ∂µ, and kµ
satisfies the null condition kµkµ = 0. The geodesic equation,
Dkµ
dλ
≡ dk
µ
dλ
+ Γµαβk
αkβ = 0, (6.2)
describes the evolution of kµ along the geodesic in the space-time described by gµν , where the
covariant total derivative D/dλ ≡ kµ∇µ. It is useful to decompose the photon 4-momentum
in the rest frame of an observer with 4-velocity uµ as
kµ = E(uµ − eµ), (6.3)
where E = −kµuµ is the energy of the photon as measured in the observer’s frame. In the
above, eµ is a unit vector describing the direction of the incoming photon on the observer’s
sky, which satisfies the constraints
gµνe
µeν = 1, eµuµ = 0, (6.4)
i.e., it is a space-like unit vector orthogonal to the observer’s 4-velocity.
A photon feels the effect of the expansion of the space-time and curvature via changes in
its energy, E, along the geodesic. This change in energy is measured as the photon redshift,
which in general is defined as the ratio of the photon energy when it was emitted at the
source, s, to when it was received at our detectors (at the observer), o, namely
1 + z ≡ Es
Eo
. (6.5)
The vector kµ fully specifies the geodesic itself, and only depends on the space-time
metric. The photon energy, and hence redshift, is defined once we choose an observer. In
mescaline, we choose our observers to be co-moving with the fluid flow. In practise, this
means that each observer’s uµ coincides with the fluid 4-velocity at their location in the
simulation. Usually, after correcting for local peculiar velocities (e.g., our motion around the
Sun and the galaxy or the peculiar motion of a nearby object) we assume our cosmological
measurements lie in the frame co-moving with the cosmic expansion. This frame might be
considered as “co-moving with the fluid flow”. However, connecting our observations to a
particular frame in the simulation is not straightforward. See section 6.5 for a discussion on
this topic.
– 8 –
J
C
A
P03(2023)019
6.1.1 Evolving the geodesic equation
We wish to numerically advance the geodesic equation (6.2) to trace the path of photons
through the simulated space-time. We can freely re-parametrise the geodesic in terms of the
simulation coordinate time, t, as follows:
dkµ
dt
= dλ
dt
dkµ
dλ
, (6.6)
= 1
k0
dkµ
dλ
, (6.7)
so long as the directional aspect of the derivative d/dλ is preserved. Namely, we must
ensure that we keep track of the change in position of the geodesic in time. In the above
formulation, we are re-casting the directional derivative along the geodesic to be spaced in
terms of coordinate time, t, instead of affine parameter, λ. Next, substituting the geodesic
equation (6.2) into the right hand side of (6.7) we arrive at the following system of equations
dk0
dt
= −Γ000k0 − Γ00jkj − kjΓ0j0 −
kk
k0
Γ0kjkj , (6.8a)
dki
dt
= −Γi00k0 − Γi0jkj − kjΓij0 −
kk
k0
Γikjkj , (6.8b)
dxi
dt
= k
i
k0
, (6.8c)
where the last equation describes the path of the geodesic through the simulation domain.
The time components of the Christoffel symbols appearing in the equations above, in terms
of 3+1 variables, are
Γ000 =
1
α
∂tα, Γ0j0 =
1
α
∂jα, Γ0jk = −
1
α
Kjk,
Γi00 = αγij∂jα, Γi0j = −αγikKjk, (6.9)
and Γijk is calculated using spatial derivatives of the metric tensor as defined in (5.6).
Substituting these into (6.8) we arrive at
dk0
dt
= −k
0
α
∂tα− 2
α
kj∂jα+
kk
αk0
Kkjk
j , (6.10a)
dki
dt
= −αk0γij∂jα+ 2αKijkj −
kk
k0
Γikjkj , (6.10b)
dxi
dt
= k
i
k0
, (6.10c)
which is the system solved in mescaline to advance the geodesic through the simulated
space-time. We describe the process of choosing initial data for kµ in section 6.3.
6.2 Distances
The geodesic equation tells us how the energy of one photon traverses through a given space-
time. This is not enough to study observables in cosmology, since many of our observations
consider distances. The most relevant distance for this work is the luminosity distance DL,
– 9 –
J
C
A
P03(2023)019
which can be defined as the ratio of the intrinsic luminosity of a source, L, to the flux received
on Earth, F , namely
DL =
√
L
4piF . (6.11)
The luminosity distance can also be defined in terms of the distance modulus µ ≡ m−M =
5log10
(
DL
10pc
)
, where m is the apparent magnitude of a source and M its absolute magnitude.
Another distance we are interested in is the angular diameter distance DA, which is defined as
D2A ≡
As
Ωo
, (6.12)
where As is the physical area of a source, and Ωo is its observed angular extent. Assum-
ing conservation of photons, these two distances are related by Etherington’s reciprocity
relation [which holds in any space-time and under any theory of gravity, see, e.g. 68]
DL = (1 + z)2DA. (6.13)
To calculate the angular-diameter distance we must track the evolution of a small bundle
of photons, or a light beam, as they traverse the space-time geometry. This beam might
be thought of as the collection of photons travelling from an extended source towards our
telescopes. We use a separation vector, ξµ, to describe the behaviour of two rays next to
one another in the bundle in order to eventually relate the physical area of the source to its
observed angular extent. Assuming the beam is narrow, we can track the evolution of the
bundle using the evolution of ξµ along the geodesic, which gives rise to the geodesic deviation
equation. We refer the reader to chapter 2 of Fleury [69] for an in-depth introduction to light
propagation and light beams, including the derivation of the geodesic deviation equation.
6.2.1 Screen space
The separation vector ξµ contains information about our infinitesimal ray bundle as it
traverses through the space-time along the geodesic. We are interested in further relating this
information to observable quantities. For this purpose, it is useful to introduce a 2-dimensional
plane (located in the 4-dimensional space-time) on which we project our light bundle. This
is known as the “screen space”, and is defined by the set of basis vectors sµA, for A ∈ [1, 2].
We can think of this 2-dimensional space as a flat plate that an observer holds perpendicular
to the light bundle, thus projecting the light onto the “screen” (much like we do with our
telescopes), at a point on the geodesic in order to measure the separation of rays within the
bundle. The equations we will solve (which we introduce in the next section) describe how
the separation between light rays projected onto this screen evolves as the bundle travels
along the geodesic. They describe how the curvature of space-time influences the shape and
size of the beam throughout its journey.
The vectors sµA form an orthonormal basis of the screen space within the space-time
defined by gµν , and must satisfy the constraints
sµAuµ = s
µ
Aeµ = 0, s
A
µ s
µ
B = δ
A
B, (6.14)
where δAB is the Kronecker delta. In words, the screen basis is orthogonal to the observers
4-velocity uµ, and the direction of observation eµ, and thus is also orthogonal to the direction
of propagation of the ray bundle (i.e. we will also have sµAkµ = 0 via (6.3)). They are also
orthogonal to each other and normalised to have unit length.
– 10 –
J
C
A
P03(2023)019
The screen vectors are essential for distance calculations along the geodesic, since the
beams morphology is calculated from the separation of geodesics in screen space. The basis
sµA must therefore also be propagated along the geodesic, however, since they are an arbitrary
basis they should not affect the physical evolution along the bundle. In appendix A, we
provide details of how we propagate the screen vectors and show that our physical results
remain unchanged with different initial sµA.
6.2.2 Jacobi matrix
Now that we are familiar with the concept of screen space, we will briefly describe the
connection of the separation vector to observables. We will only briefly touch on the derivation
of the equations involved, mainly for interest and physical intuition of the equations we solve
in mescaline.
Projecting the separation vector ξµ into the screen space gives the separation of rays in
the bundle in the 2-dimensional plane, namely ξA ≡ ξµsµA. In combination with the geodesic
deviation equation, ξA will allow us to track intrinsic properties of the beam and relate these
to quantities measured by observers in the space-time.
After this projection, the geodesic deviation equation gives rise to the Sachs vector
equation [see 69]
d2ξA
dλ2
= RABξB, (6.15)
where we have introduced the optical tidal matrix
RAB =
(
R 0
0 R
)
+
(
−Re(W ) Im(W )
Im(W ) Re(W )
)
, (6.16)
and the positioning of the indices A,B does not matter since they are raised and lowered with
δAB, however, summation is still implied over repeated indices. In (6.16), the Ricci lensing
scalar is
R ≡ −12Rµνk
µkν , (6.17)
which describes the focusing of the light beam due to matter in its path. The Weyl lensing
scalar is
W ≡ −12Cµναβσ
µkνkασβ, (6.18)
where σµ ≡ sµ1 − isµ2 (with i2 = −1), which describes the shearing of the light beam due to
nearby structures outside of the beam itself. This is natural given that the inherent nature of
the Weyl tensor describes the distortion of bodies (which in this case is the light beam) due
to external tidal forces. We discuss the method for calculating W in section 6.2.4 below.
The linearity of (6.15) implies that there exists a direct mapping between any solution
and its initial conditions. If we consider a beam of light that converges at an observer at
point o, and integrate (6.15) outwards to a source at point s, we can show that the physical
separation of two rays in the bundle at the source, ξAs , is
ξAs = DAB(s← o)
dξB
dλ
∣∣∣∣
o
, (6.19)
where s← o indicates the direction of integration. The derivative on the right hand side can
be written as [see 69]
dξB
dλ
∣∣∣∣
o
= −EoθBo , (6.20)
– 11 –
J
C
A
P03(2023)019
where θBo is the angular separation of the rays on the observer’s sky. Therefore, in (6.19),
DAB represents the Jacobi matrix mapping of the physical separation of rays in the bundle at
the source, ξAs , to the observed angular separation of the same rays at the observer, θBo . This
matrix therefore contains all information relating the physical attributes of a source to how it
appears to an observer. Especially relevant to this work is its determinant, which is the ratio
of the physical area of the source to its observed angular extent,
Eo det|D| = AsΩo , (6.21)
which is the definition of the angular diameter distance, DA, as given in (6.12). The Jacobi
matrix can also be decomposed in terms of quantities describing image shear and rotation due
to gravitational lensing [see section 2.2.2 of 69]. Evolving the Jacobi matrix equation directly
thus also allows for straightforward analysis of weak lensing observables using the same data.
We note that evolving a bundle of geodesics via the Jacobi matrix equation (or any
of its decompositions) is not the only viable method to calculate distances in simulation
data. One might instead choose to evolve the geodesic equation (6.2) individually for a set of
neighbouring light rays and calculate their separations ξA along the way [e.g. 70–74].
6.2.3 Evolving the Jacobi matrix equation
Interpreting the Jacobi matrix as relating observations of a source to its physical properties
relies on the evolution of the Sachs vector equation from an observer out to the source. We can
obtain the evolution of the Jacobi matrix along the geodesic by taking the second derivative
of (6.19) and inserting the Sachs vector equation (6.15), which gives its evolution in terms of
the optical tidal matrix [69]
d2
dλ2
DAB = RACDCB. (6.22)
As a next step, we might decompose the Jacobi matrix to define the Sachs optical scalars
θS and σS ,6 which describe the expansion and shear rates of the beam, respectively [75].
The expansion rate θS is singular at the point of observation, λ = λo, and so it is not ideal
for numerical integration. It is therefore common for ray-tracing in numerical data to cast
the evolution equations for the optical scalars into a different form to instead evolve, e.g.,
the area of the beam [e.g. 44] or the angular diameter distance directly [e.g. 76]. See also
Grasso et al. [77], Grasso and Villa [78] for light propagation methods using bi-local geodesic
operators [79]. However, here we do not take either of these routes. Instead, we directly
integrate the evolution equation (6.22) to obtain the full Jacobi matrix along the geodesic.
We now need to write (6.22) in a form suitable for numerical integration. First, as with
the geodesic equation, we will re-cast the derivative with respect to λ into a derivative with
respect to the time coordinate, t. Thus, we can write
d2
dλ2
= d
dλ
(
k0
d
dt
)
(6.23)
= dk
0
dλ
d
dt
+ k0 d
dλ
d
dt
(6.24)
= −Γ0αβkαkβ
d
dt
+
(
k0
)2 d2
dt2
(6.25)
6We add a subscript S here to avoid confusion with the kinematic fluid expansion and shear, however, these
scalars are usually referred to without this subscript.
– 12 –
J
C
A
P03(2023)019
where to get to the last line we have used the geodesic equation (6.2). Therefore, (6.22)
becomes
− Γ0αβkαkβ
d
dt
DAB + (k0)2
d2
dt2
DAB = RACDCB, (6.26)
and to simplify the system we define
PAB ≡
d
dt
DAB, (6.27)
such that we can re-write (6.26) to arrive at the following system of evolution equations
d
dt
PAB =
1
(k0)2
Γ0αβkαkβPAB +RACDCB, (6.28)
d
dt
DAB = PAB. (6.29)
Expanding the sum Γ0αβkαkβ in (6.28) for a 3+1 metric, we arrive at the full system of
evolution equations
d
dt
PAB =
[ 2
αk0
ki∂iα− 1
α(k0)2Kijk
ikj + 1
α
∂tα
]
PAB +
1
(k0)2R
A
CDCB, (6.30a)
d
dt
DAB = PAB, (6.30b)
and the initial conditions at the point of observation are [69]
DAB(λo) = 0, (6.31a)
PAB(λo) = δAB. (6.31b)
The system of equations (6.30) combined with the evolution equations for kµ in (6.10)
and the screen vector evolution in appendix A fully specifies the evolution of the Jacobi matrix
along the geodesic. From this system we can extract the angular diameter distance and the
redshift at all points along the geodesic for a given numerical metric tensor.
6.2.4 Weyl tensor
The Weyl tensor appears in the optical tidal matrix (6.16) in the Weyl lensing scalar (6.18).
It is therefore an important part of the evolution of the Jacobi matrix. The Weyl tensor is
defined as the trace-free part of the Riemann tensor, namely,
Cµναβ ≡ Rµναβ + 12
[
gµβRνα + gναRµβ − gµαRνβ − gνβRµα
]
+ 16
(
gµαgνβ − gµβgνα
)
R.
(6.32)
To calculate the Weyl curvature scalar W , we adopt a method similar to that described
in appendix B of Giblin et al. [44]. We begin by collapsing the space-time indices in (6.18)
into a new set of indices, namely
W = −12CΣΠX
ΣY Π, (6.33)
with XΣ ≡ σµkν and Y Π ≡ kασβ , and the index Σ therefore takes values over all combinations
of µν, and Π takes values over all combinations of αβ. Substituting the definition of the Weyl
– 13 –
J
C
A
P03(2023)019
tensor into (6.33), all terms in brackets in (6.32) are zero either due to the constraints on the
screen vectors (6.14) (i.e. after contracting with σµ ≡ sµ1 − isµ2 ), or because they cancel with
one another. The expression therefore simplifies to
W = −12RΣΠX
ΣY Π, (6.34)
where RΣΠ ≡ Rµναβ . This expression is symmetric in changes Σ↔ Π due to the symmetries
in the Riemann tensor. All twelve possible values for the indices are Σ,Π ∈ [01, 02, 03, 12,
13, 23, 10, 20, 30, 21, 31, 32], because combinations with repeated indices in either Σ or Π
(i.e., Σ = 00, 11, etc.) have zero Weyl tensor component. We define the ‘forward’ indices as
Σf ,Πf ∈ [01, 02, 03, 12, 13, 23], and the ‘backward’ indices as their reverse values Σb,Πb ∈
[10, 20, 30, 21, 31, 32]. The Riemann tensor RΣΠ is antisymmetric in exchanges Σf → Σb or
Πf → Πb. Expanding the sum in index Σ in (6.34) therefore gives
RΣΠX
ΣY Π = RΣfΠXΣfY Π +RΣbΠXΣbY Π, (6.35)
= RΣfΠ(XΣf −XΣb)Y Π, (6.36)
= 2RΣfΠX [Σf ]Y Π, (6.37)
since RΣbΠ = −RΣfΠ and we have denoted the antisymmetric part of XΣf as X [Σf ] ≡
1
2(XΣf −XΣb). Expanding the Π index in (6.37) then gives
2RΣfΠX [Σf ]Y Π = 2RΣfΠfX [Σf ]Y Πf + 2RΣfΠbX [Σf ]Y Πb ,
= 4RΣfΠfX [Σf ]Y [Πf ], (6.38)
following similar steps. The Weyl curvature scalar finally becomes
W = −2RΣΠX [Σ]Y [Π], (6.39)
where we have dropped the subscript f in the indices for brevity. We next need to split this
expression into its real and imaginary parts to calculate the tidal matrix components in (6.16).
The imaginary part comes from the vector σµ, and therefore is in XΣ and Y Π. We split these
into
XΣ = XΣ1 − iXΣ2 , (6.40a)
Y Π = Y Π1 − iY Π2 , (6.40b)
where
XΣ1 ≡ sµ1kν , XΣ2 ≡ sµ2kν , (6.41)
Y Π1 ≡ kαsβ1 , Y Π2 ≡ kαsβ2 . (6.42)
Substituting (6.40) into (6.39) gives
Re(W ) = −2RΣΠ
(
X
[Σ]
1 Y
[Π]
1 −X [Σ]2 Y [Π]2
)
, (6.43)
Im(W ) = 2RΣΠ
(
X
[Σ]
1 Y
[Π]
2 +X
[Σ]
2 Y
[Π]
1
)
, (6.44)
with Σ,Π ∈ [01, 02, 03, 12, 13, 23]. This calculation of the Weyl scalar is completely general
for any metric gµν , and only depends on the constraints on the screen vectors (6.14) and the
– 14 –
J
C
A
P03(2023)019
null condition for the photon 4-momentum (i.e., when either set of vectors are contracted
with gµν) as well as some symmetries of the Weyl tensor.
The components of the Riemann tensor we need for these values of Σ,Π can be written
in terms of 3+1 variables, for zero shift vector, as
R0i0j = α∂0Kij + α∂i∂jα− α∂k(α)Γkij + α2KkjγklKli, (6.45a)
R0ijk = α∂jKik − α∂kKij + αKljΓlik − αKlkΓlij , (6.45b)
Rijkl = γim
[
∂kΓmjl − ∂lΓmjk + ΓmnkΓnjl − ΓmnlΓnjk
]
+KikKjl −KilKjk. (6.45c)
6.3 Initial data
To evolve our system and trace the path of the geodesic back in time, we first need to specify
initial data for kµ, xµ,DAB, PAB, and the screen vectors sµA for each geodesic at the observer
position. The initial data for xµ is just the position of the observer and the initial data for
the Jacobi matrix was given in (6.31). The initial kµ will be different for each geodesic and is
dependent on what kind of sky sampling we are interested in. In general, we need to translate
a set of lines of sight as seen by the observer, e.g. the coordinates of a particular survey,
into an initial set of kµ. From the decomposition (6.3), the photon 4-momentum kµ is fully
specified by the observer 4-velocity, uµ, and the direction of the source on the sky, eµ, (up to
a scaling by the energy of the photon, E) both of which are arbitrary so long as the relevant
constraints are satisfied for eµ. We take our observers to be co-moving with the fluid flow
and thus uµ is specified using the 4-velocity of the fluid at the observers location. To set the
initial eµ, we first write the spatial components as
ei = Dmi, (6.46)
where mi is some spatial direction vector and D is a scaling factor to be determined from
the constraints. In the frame of the observer, eµ is purely spatial and thus e0 = 0. However,
the simulation is not necessarily performed in the rest frame of the observer and therefore
the coordinates {t, xi} are in general not adapted to the fluid flow. We will in general have
e0 6= 0, and we must also determine this from the constraints on eµ. We take the vector mi
to represent the {x, y, z} Cartesian coordinates of the direction of an incoming geodesic on
the observer’s sky. Given an initial mi, we determine D and e0 by substituting (6.46) into
the constraints (6.4), which gives
e0 = ±m
iui√
u20γijm
imj − α2mimjuiuj
, (6.47a)
D = ∓u0√
u20γijm
imj − α2mimjuiuj
, (6.47b)
where α, γij , and ui are evaluated at the observer’s position in the simulation [see also
appendix B3 of 38].
Mescaline takes a set of directions mi as defined in the observer’s frame representing,
e.g., a mock supernova survey or an even coverage in the direction of HEALPix7 indices [see
80]. To ensure our initial data is defined in the same frame used to advance the geodesics, we
first transform these mi into the simulation frame prior to calculating eµ (see appendix D).
7https://healpix.sourceforge.io.
– 15 –
J
C
A
P03(2023)019
After determining the components of eµ using (6.47), we then use the decomposition (6.3),
with an arbitrary choice of the photon energy at the observer E → Eo (since we are interested
in the ratio of the energies at observer and source), and the value of the 4-velocity at the
observer’s position, to determine the initial kµ.
We detail the process of generating initial data for the screen vectors in appendix A.2.
Their initial values are built from the constraints (6.14) after specifying the three remaining
degrees of freedom arbitrarily. This choice does not affect the final physical result, which we
show in appendix A.3.
6.4 Time stepping, interpolation, and tests
The equations are advanced in time using a Runge-Kutta 2nd order (RK2, or Heun) method,
and all spatial derivatives are calculated using fourth-order accurate stencils. The evolution
tracks the propagation of light beams through a discretised Cartesian grid at the specific
coordinate time intervals defined by the constant-t surfaces of the simulation output. The
position of the light beam at each time step will not necessarily be positioned on a spatial grid
cell. We therefore need to interpolate the quantities on the right hand side of the evolution
equations to the position of the beam at each time step. In mescaline, we adopt a 5th order
spatial, 3-dimensional B-spline interpolation (see appendix C).
In appendix E we present a variety of tests of the numerical accuracy of the mescaline
ray tracer. We test the redshift and distance calculations for analytic EdS and linearly-
perturbed EdS space-time metrics and confirm the error converges at the expected rate. We
also test the calculation on a set of controlled-mode NR simulations and perform a Richardson
extrapolation to estimate the error in our calculations for simulation data. Lastly, we test
the violation in the null condition kµkµ = 0 — which is not explicitly enforced anywhere in
mescaline and therefore serves as an additional numerical test — and confirm it reduces
with increasing resolution and is thus dominated by numerical error. Based on these tests, we
conclude the mescaline ray tracer is accurate and as precise as possible for a second-order
accurate integration scheme.
Each step of the ray-tracing calculation is explicitly propagating the light bundle onto
the next constant-t hypersurface of the simulation, and therefore each step for all observers
is associated with the same coordinate time. However, due to the inhomogeneous nature of
the simulated space-time, for each observer each step will correspond to a different redshift.
Further, not all lines of sight will have the same redshift at each time step. To isolate the
anisotropies in DL itself, we can interpolate each line of sight to the same z value after the
analysis is done. In this work, we use the mean value of z across the sky for each observer at
each step and linearly interpolate each line of sight to find DL at this redshift.
6.5 Observers
Quantities calculated via ray tracing are purely observable and are thus directly comparable
to distances and redshifts from cosmological surveys. However, they are dependent on our
choice of observer and so this choice requires some discussion. For the framework presented
here, the observers are moving with some general 4-velocity uµ and in mescaline we have
further assumed that this coincides with the fluid 4-velocity. This is a common choice in both
analytic and numerical calculations of observables which might manifest either by explicitly
taking both observer and emitter to be co-moving with the fluid flow [e.g. 40, 42, 81, 82], or
via the use of the synchronous co-moving gauge with observers at rest with respect to the
coordinates [e.g. 21, 44, 77, 83, 84].
– 16 –
J
C
A
P03(2023)019
Ideally, we want to mimic real cosmological observations as closely as possible with
our ray-tracing data. Relating the frame which is co-moving with the fluid flow within a
simulation to observational analysis is perhaps not straightforward. For example, in analysis
of observations we might assume that the rest frame of the CMB is the frame in which
the Universe is well-described by the FLRW models. The redshift of distant sources will
be dominated by cosmic expansion, whereas low-redshift sources can have a significant
contribution from peculiar velocities. For the latter, we might apply corrections to our
observed redshifts in an attempt to alleviate this and end up with a purely cosmological
redshift [see, e.g. 85]. Additionally, we as observers are also moving with respect to this
global cosmological frame and this must also be taken into account from observed redshifts.
Our motion is inferred from the kinematic dipole we observe in the CMB radiation [86],
though recently the purely kinematic origin of this dipole has been called into question [e.g.
14, 56, 87] [see also 88].
The use of the CMB frame could be motivated by the fact that the early Universe
was very close to homogeneous and isotropic, and therefore close to the reference FLRW
model which best fits the CMB radiation. The simplest way to define the rest frame of the
CMB in NR simulations is to use the dipole that each observer measures in their CMB map.
However, this would require ray-tracing the full sky to very high redshift which is beyond
the computational capabilities of this work. Since our simulations start close to an FLRW
model by construction, we might choose the simulation hypersurface to coincide with the
“CMB frame” for all observers. In fact, the perturbations to the metric tensor itself remain
small and our hypersurfaces thus remain close to the FLRW expansion chosen on the initial
slice. However, large perturbations in the matter develop during the simulation, and so the
frame co-moving with the fluid is not necessarily close to this FLRW model on all scales. We
thus might consider using the observer’s peculiar velocity with respect to the Eulerian grid
(i.e., the hypersurfaces) to boost their measured redshifts to this fictitious “CMB frame” for
analysis. We proceed using observers who are co-moving with the fluid flow for the remainder
of this work. However, in appendix B we outline the process by which we can correct our
simulation data using a single point-wise boost to the hypersurface frame. We also present
some of our main results after applying this boost for comparison.
For all calculations, observers are placed on the spatial hypersurface of the simulation
closest to z = 0. This slice is chosen using the effective scale factor of the simulation, aD,
and the initial redshift of the evolution. Observer positions are chosen (pseudo-)randomly in
space across the redshift zero slice.
7 Cosmography calculation
We calculate the luminosity distance as predicted by the general cosmography given in section 3
using the coefficients (3.1), using the effective cosmological parameters as calculated and
presented in MH21. We will briefly summarise the method used to calculate these parameters
in this section.
In MH21, the authors first calculated the effective Hubble parameter, using its multipole
decomposition given in (3.3) for each observer. The first step is to evaluate the volume
expansion θ, the 4-acceleration aµ, and the shear tensor σµν at the observer’s location (explicit
definitions of these quantities are given in appendix B of MH21). We then use these together
with the unique eµ for each line of sight (via the same process outlined in section 6.3) to
calculate H(e) across each observer’s sky. The effective deceleration parameter Q in (3.2b)
– 17 –
J
C
A
P03(2023)019
and jerk J in (3.2d) are then calculated using the first and second derivatives of H, respectively,
along the direction of the incoming null ray (see appendix B of MH21 for the exact expressions
of these derivatives). Finally, the effective curvature parameter R in (3.2c) is calculated using
the 4-Ricci tensor, at the observer’s location, contracted with kµ for each line of sight (kµ is
also calculated via the method given in section 6.3). All calculations were also done using
mescaline, though via a separate set of routines to the ray-tracing code.
In this work, we use the calculations of H,Q,R, and J from MH21 and substitute them
into the Taylor series expansion (2.1) to find the cosmographic dL(z) across the sky for all
100 observers (correct to third order in redshift). We calculate dL(z) for a set of redshift
values for each observer which coincide with the mean z across the sky as output from the ray
tracer (this is the same value which we interpolate DL(z) to create the ray-traced sky maps).
8 Results and discussion
In section 8.1, we compare the DL from ray tracing to the cosmographic dL for the simulations
presented in MH21. In section 8.2 we assess the accuracy of the cosmographic prediction in
a simulation containing more small-scale structure, but which still contains relatively low
density contrasts. In section 8.3, we further investigate the anisotropy in a simulation with as
much small-scale structure as we can reliably resolve with the simulation software we use.
In our results presented below, we often present the luminosity distances normalised by
the relevant FLRW luminosity distance. For our matter-dominated simulations, this is the
EdS model which has cosmological parameters ΩΛ = Ωk = 0 and Ωm = 1. This is the model
around which we apply small perturbations to set our initial data. We have
dL(EdS) =
2
H0
(
1 + z −√1 + z
)
, (8.1)
where H0 is the Hubble constant. We take the globally-averaged Hubble parameter in the
simulation Hall ≡ 〈θ〉all/3 as H0. Here, 〈θ〉all ≡
∫
all
√
γ θ d3X/Vall is the volume average
over the entire z = 0 slice, with Vall =
∫
all
√
γ d3X is the volume of the slice and γ is the
determinant of the spatial metric γij . As mentioned in section 4 above and in MH21, we
find this Hubble rate to coincide with the FLRW model prediction of H0 = 100h km/s/Mpc
(approximately 45 km/s/Mpc for the EdS model) to within one percent for all simulations we
use here.
Throughout our results, we generally distinguish between the ray-traced luminosity
distance and the cosmographic prediction for the luminosity distance (correct to third order
in redshift) using DL for the former and dL for the latter.
8.1 Strictly large-scale simulations
Here we compare the cosmographic predictions based on the simulations presented in MH21 to
our new calculations using ray tracing in full GR. These simulations enforce strict smoothing
scales by explicitly neglecting all structures beneath a chosen scale. There are two simulations
with resolution N = 128 and smoothing scales of 100 h−1 Mpc and 200 h−1 Mpc.
Figure 1 shows a comparison of the ray-traced luminosity distance DL (top row) to the
cosmographic dL (bottom row) for one observer in the simulation with a strict 200 h−1 Mpc
smoothing scale. Panels, left to right, show three redshift slices z ≈ 0.034, 0.07, and 0.107,
respectively. All distances are normalised by the EdS luminosity distance for that redshift. By
eye, the cosmographic dL appears to give a very good approximation to the exact DL at least
– 18 –
J
C
A
P03(2023)019
DL / dL(EdS) 1, z 0.034
-0.011 0.013
DL / dL(EdS) 1, z 0.070
-0.013 0.01
DL / dL(EdS) 1, z 0.107
-0.01 0.01
dL(3rdorder) / dL(EdS) 1, z 0.034
-0.011 0.013
dL(3rdorder) / dL(EdS) 1, z 0.070
-0.013 0.01
dL(3rdorder) / dL(EdS) 1, z 0.107
-0.01 0.01
Figure 1. Sky maps of the luminosity distance for one observer in the simulation with 200h−1
Mpc strict smoothing scale. Top row shows the ray-traced luminosity distance DL at three redshifts
z ≈ 0.034, 0.07, and 0.107 (left to right, respectively), normalised by the EdS luminosity distance at
that same redshift. Bottom row shows the cosmographic luminosity distance dL correct to third order
in z for the same three redshifts, also normalised by the EdS distance.
0.04 0.02 0.00 0.02 0.04 0.06
o
0.10
1.00
L 1
,d
iff
(D
L)
(%
)
z 0.034
RT vs cosmography
RT vs EdS
0.04 0.02 0.00 0.02 0.04 0.06
o
z 0.069
0.04 0.02 0.00 0.02 0.04 0.06
o
z 0.107
Figure 2. Percentage L1 difference between the ray-traced DL and the cosmographic dL (blue circles,
dashed line average), and the difference between the ray-traced DL and the EdS dL (red crosses, dotted
line average). We show the differences for 100 observers at three redshift slices (panels) as a function
of the density contrast at the observer, δo. All observers are randomly placed in the simulation with a
strict 200h−1 Mpc smoothing scale.
until z ≈ 0.07. To quantify the accuracy of the approximation, we calculate the difference
between the two distances along each line of sight for all 100 observers at three redshifts slices.
Specifically, we calculate the “L1 difference”, defined as
L1,diff(DL) ≡ 1
NLOS
NLOS∑
i=1
∣∣∣∣∣DiLdiL − 1
∣∣∣∣∣ , (8.2)
where DiL is the value of DL for line of sight i (with NLOS total lines of sight), and we
take either dL = dL(3rdorder) or dL = dL(EdS) to compare DL to both the cosmographic
prediction and that of the EdS model, respectively.
– 19 –
J
C
A
P03(2023)019
0.1 0.0 0.1 0.2
o
0.1
1.0
10.0
L 1
,d
iff
(D
L)
(%
)
z 0.034
RT vs cosmography
RT vs EdS
0.1 0.0 0.1 0.2
o
z 0.069
0.1 0.0 0.1 0.2
o
z 0.107
Figure 3. Same as figure 2 except for a different set of 100 observers in the simulation with the strict
100h−1 Mpc smoothing scale.
Figure 2 shows the percentage L1 difference between the cosmographic dL and the
ray-traced DL (blue circles) for 100 observers as a function of their local density contrast,
δo ≡ ρo/〈ρ〉all − 1. The red crosses show the difference between DL and the EdS distance for
the same observers. Horizontal lines show the mean over all observers for the relevant case.
Panels, left to right, show three redshift slices of z ≈ 0.034, 0.069, and 0.107, respectively. We
note the mean redshift on the slice differs slightly between observers. All observers shown
here are placed randomly in the simulation with a strict 200h−1 Mpc smoothing scale. The
cosmographic distance gives a sub-percent accurate representation of DL at all three redshifts
for all observers. For z . 0.07, the luminosity distance is captured to within ∼ 0.2% for all
observers. Due to the very low density contrasts in the simulation, the difference between DL
and the EdS distance is always ∼ 1%. However, the cosmographic dL still performs better
than the EdS relation for z . 0.1.
In appendix A of MH21, the authors predicted that the cosmography (as truncated at
third order in redshift) should represent the exact luminosity distance to within 1% at z ∼0.04–
0.06 for a strict smoothing scale of 200h−1 Mpc (as predicted based on this exact simulation).
We find the cosmographic prediction to actually perform better than this prediction, matching
within 1% even out to z ≈ 0.1 for all 100 observers.
The results so far have been based on the simulation in MH21 with strictly no structure
beneath 200 h−1 Mpc. Next, we will consider the second simulation presented in MH21: which
instead contains no structure beneath 100 h−1 Mpc, and is thus slightly less smooth. Figure 3
shows the percentage L1 difference between DL and the cosmographic dL (blue circles, dashed
line average) for a set of 100 observers in the simulation with a strict smoothing scale of
100h−1 Mpc. Again, we show differences as a function of the local density contrast at the
observer’s location. We note the higher density contrasts on the x-axis with respect to figure 2,
since this simulation has a smaller smoothing scale. Red crosses again show the difference
between DL and the EdS distance for the same three redshift slices as in figure 2. Here we
see the accuracy of the cosmography is better than 1% for all observers only to z ≈ 0.034,
reaching a few percent at z ≈ 0.069 and eventually ∼ 10 percent at z ≈ 0.107. We notice the
ray-traced distance DL is always within a few percent of the EdS distance for all observers out
to z ≈ 0.1. The cosmographic distance gives about an order-of-magnitude better prediction
than the EdS relation at low redshift, due to the incorporation of anisotropies induced by
the smaller-scale structures. However, as we approach higher redshifts the > O(z3) terms
– 20 –
J
C
A
P03(2023)019
0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
/
Figure 4. 2-dimensional slice of the density field at z ≈ 0 in a simulation with box size 3h−1Gpc and
resolution N = 256. All modes beneath ∼ 200h−1Mpc are removed from the initial data. We show
the matter density normalised by the average over the whole domain, ρ/ρ¯.
become important and the cosmographic relation (as truncated at third order) over-estimates
the anisotropy at z & 0.07 (which can also be seen in the right-most panel in figure 1 for the
case of the 200h−1 Mpc smoothing scale).
The prediction from MH21 for a 100 h−1 Mpc smoothing scale was that the cosmographic
dL should be accurate to within 1% for redshifts of z ∼0.02–0.03. We find our results are
consistent with, and even slightly better than, this prediction.
For both simulations we discussed in this section, the ray-traced distance is within a few
percent of the EdS distance for all redshifts we analysed. However, we stress again that all
structure beneath 100 h−1 Mpc or 200 h−1 Mpc has been strictly ignored (and thus cannot
affect the evolution of larger scales at all), which is not a realistic universe scenario. In the
next section, we perform a similar comparison for a simulation with more small-scale structure
to demonstrate the effects of local inhomogeneities on the ability of the cosmography (and
the EdS model) to capture distances.
8.2 A more realistic model universe
In this section, we will assess the ability of the cosmographic luminosity distance to capture
the ray-traced distance in a more realistic model universe. We use a simulation with box
length L = 3h−1Gpc and resolution N = 256, which was also performed using the ET with
initial data from FLRWSolver generated in a manner identical to that described in section 4.
However, for this simulation we have removed some small-scale structure from the initial
power spectrum to impose a “soft” smoothing scale. This means that initially our simulations
contain no structure beneath a certain scale. However, contrary to the previous section, we do
not strictly prevent structure beneath this scale from forming later in the simulation. In this
– 21 –
J
C
A
P03(2023)019
DL / dL(EdS) 1, z = 0.012
-0.027 0.073
DL / dL(EdS) 1, z = 0.024
-0.026 0.058
DL / dL(EdS) 1, z = 0.040
-0.025 0.042
dL(3rdorder) / dL(EdS) 1, z = 0.012
-0.034 0.068
dL(3rdorder) / dL(EdS) 1, z = 0.024
-0.067 0.067
dL(3rdorder) / dL(EdS) 1, z = 0.040
-0.25 0.073
Figure 5. Sky maps of the luminosity distance for one observer in a simulation with all structure
beneath 200h−1 Mpc removed from the initial data. Top row shows the ray-traced luminosity distance
DL at three redshifts z = 0.012, 0.024, and 0.04 (left to right, respectively), normalised by the EdS
distance. Bottom row shows the cosmographic luminosity distance dL correct to third order for the
same redshifts, also normalised by the EdS distance. Note that the colour bars between the top and
bottom rows differ.
section, we will study a simulation with all power below ∼ 200h−1Mpc initially removed. This
implies that the minimum mode in the initial data is sampled by about 17 grid cells. In practise
we remove this small-scale power by choosing a scale kcut = 2pi/λcut with λcut = 200h−1Mpc
and set P (k > kcut) = 0, with P (k ≤ kcut) left as the power spectrum output from CLASS
for k ≤ kcut. In section 8.3 below, we will study an identical simulation (same domain size,
resolution, and random realisation) with a slightly lower cut of λcut = 100h−1Mpc. Figure 4
shows a 2-dimensional slice of the density field (relative to the mean over the whole slice)
at redshift z ≈ 0 for the simulation with initial λcut = 200h−1Mpc. Here we see ∼ 10–20%
typical under-dense regions and ∼ 30–50% typical over-dense regions. Since this simulation is
matter-dominated, typical density contrasts on a given scale will in general be larger than in
the ΛCDM model due to the enhanced growth rate of structure. However, as also discussed in
MH21, we expect the qualitative anisotropic signatures to be the same as in a ΛCDM model,
though with generally larger amplitudes. Again, while all structures beneath ∼ 200h−1 Mpc
were removed from the initial data, the initial perturbations do collapse to scales smaller than
this by the end of the simulation. It is thus difficult to refer to a specific “smoothing scale” at
z ≈ 0, and so we refer to the simulations by their initial smoothing scales instead.
We place 25 observers and perform the ray tracing to obtain DL to z ≈ 0.1. We use fewer
observers for this simulation, relative to those in section 8.1, since it is higher resolution and
thus the calculations are more expensive. We also calculate the cosmographic dL for the same
set of observers using the same process as outlined in section 7 and in MH21. For both cases, we
calculate distances for 12×N2side lines of sight in directions of HEALPix indices with Nside = 32.
Figure 5 shows a comparison of the ray-traced DL (top panels) and the cosmographic
dL (bottom panels) at three redshift slices z = 0.012, 0.024, and 0.04 (left to right rows,
respectively) for one observer in the simulation with initial λcut = 200h−1 Mpc. We note the
– 22 –
J
C
A
P03(2023)019
0.1 0.0 0.1 0.2
o
0.1
1.0
10.0
L 1
,d
iff
(D
L)
(%
)
z 0.012
RT vs cosmography
RT vs EdS
0.1 0.0 0.1 0.2
o
z 0.024
0.1 0.0 0.1 0.2
o
z 0.040
Figure 6. Percentage difference between the ray-traced DL and the cosmographic dL (blue circles,
dashed line average), and the ray-traced DL and the EdS dL (red crosses, dotted line average), at three
redshift slices (panels). We show the difference for 25 observers as a function of the density contrast
at their location, δo. The observers are placed in the simulation with all power below 200h−1Mpc
removed from the initial data.
lower redshift slices here relative to the previous section, as the presence of more small-scale
structure spoils the ability of the cosmography to capture DL to higher redshift. We also
point out that the colour bars between the top and bottom rows of figure 5 are different
for DL and dL. Here we notice that the cosmography reproduces the anisotropic signatures
for this observer quite well to z ≈ 0.024, however, the amplitudes of these anisotropies are
mostly over estimated. We also notice that at the lowest redshift z ≈ 0.012 (left-most panels),
we have a quadrupole dominating the anisotropy in the luminosity distance, while as we
move to higher redshift of z ≈ 0.04, we see some higher-order multipoles mixed in as well
as a discernable dipole signature. Figure 6 shows the percentage L1 difference between the
ray-traced DL and the cosmographic dL (blue circles with average as a dashed line) and the
ray traced DL and the EdS dL (red crosses with average as a dotted line) for this case. For
z ≈ 0.012, the cosmography approximates the ray-traced DL to within 1% for all observers
and to within 0.5% for ∼ 70% of observers. The ray-traced distance is consistent with the
EdS distance to within ∼ 2–3% for most observers for z . 0.04. For redshift z ≈ 0.024 the
cosmography provides a ∼ 2% accurate prediction for DL on average across all observers. At
z ≈ 0.04, this has increased to ∼ 7%.
Even though we have included smaller scales than the simulations presented in section 8.1,
we are still neglecting a lot of small-scale structure in the simulations analysed in this section.
Due to the fluid approximation adopted in the code we use, we can only reliably simulate
physical scales above those of the largest galaxy clusters (∼8–10 h−1 Mpc). In the next
section, we will analyse a simulation sampling down to these smallest physically-reliable scales.
8.3 An even more realistic model universe
Here we assess the level of anisotropy in a simulation with even more small-scale structure than
that presented in the previous section. Namely, we remove all modes beneath ∼ 100h−1 Mpc
from the initial data (such that the minimum mode in the initial data is sampled by about 8
grid cells). The power spectrum and random seed used to generate the initial data (and all
parameters associated with setting initial data and evolution of the simulation) are identical
to the simulation in the previous section, and only the initial kcut differs between the two.
Figure 7 shows a 2-dimensional slice of the logarithm of the matter density field in the
– 23 –
J
C
A
P03(2023)019
0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0
log10( / )
Figure 7. 2-dimensional slice of the density field at z ≈ 0 in a 3 h−1 Gpc simulation with all modes
beneath ∼ 100h−1Mpc removed from the initial data. We show the logarithm of the matter density
normalised by the average over the whole domain, ρ/ρ¯.
simulation used in this section on the spatial slice at redshift z ≈ 0. Comparing to figure 4,
we can see the increase in small-scale structure and larger density contrasts by about an order
of magnitude.
In the previous section, we found that the cosmographic prediction was inaccurate at
the ∼ 10% level beyond z ≈ 0.04. We therefore expect it to provide an even less accurate
prediction for this simulation for these low redshifts. Further, the validity of the generalised
cosmography relies on the condition H 6= 0 at all points along the geodesic [see 37, for more
details and a discussion of the implications of this assumption]. Thus, the cosmography
applies only to large-scale cosmological models without collapsing structures. The simulation
shown in figure 7 contains non-linear structure, and thus collapsing regions. We therefore
expect H might change sign along some lines of sight, and the validity of the cosmographic
prediction breaks down. Accurately predicting the anisotropic features in this simulation using
the general cosmography therefore must involve some smoothing of small-scale structures in
order to build a large-scale model from this simulation, which is beyond the scope of this
work. In Macpherson and Heinesen [89, in prep.], we use these same simulations to assess
the ability of the general cosmography to predict the dipolar signature in DL for various
smoothing scales. Studying even smaller scales than those sampled in figure 7 — eventually
reaching a sampling which mimics our true observations — will require going beyond the fluid
approximation and using an N-body particle cosmological code such as gevolution [30].
Instead of directly comparing to the cosmographic prediction, we will present and assess
how the ray-traced luminosity distance varies across the sky in this simulations. We defer an
involved comparison with cosmographic predictions for specific multipole signatures to future
work [89, in prep.].
– 24 –
J
C
A
P03(2023)019
DL / dL(EdS) 1, z 0.032
-0.13 0.12
DL / dL(EdS) 1, z 0.070
-0.058 0.056
DL / dL(EdS) 1, z 0.101
-0.048 0.043
DL / dL(EdS) 1, z 0.301
-0.016 0.012
Figure 8. Sky maps of the ray-traced luminosity distance DL, relative to the EdS prediction, on four
uniform redshift slices z ≈ 0.031, 0.071, 0.101, and 0.301 (top left to bottom right, respectively) for one
observer in the simulation shown in figure 7.
We place 25 observers at the same locations as for the simulation in section 8.2, and ray
trace for the same lines of sight. Figure 8 shows all-sky maps of the ray-traced luminosity
distance DL, relative to EdS, for one observer in the simulation with structure beneath
∼ 100h−1 Mpc cut out of the initial data (as shown in figure 7). We show four uniform
redshift slices at z ≈ 0.031, 0.071, 0.101, and 0.301 in the top-left to bottom-right panels,
respectively. The maps here are just for one observer’s sky. The local environment of all
observers will yield different anisotropic signatures, since the effects themselves arise from
local kinematics such as differential expansion and shear. Figure 9 shows the ray-traced
luminosity distance as a function of redshift for a set of lines of sight for the same observer
as shown in figure 8. We show the deviation of DL from the EdS distance for 100 randomly
chosen lines of sight (coloured dashed curves), the mean over the whole sky (thick solid black
curve), and the mean over the 100 lines of sight shown here (thick dashed black curve). The
inset shows a Mollweide projection of the directions of the 100 lines of sight shown here on the
observer’s sky. The dark and light grey shaded regions show a ±2% and ±5% variation from
EdS, respectively. This observer is located in a region with density contrast of δo ≈ −0.33, i.e.,
a ∼ 30% under-density relative to the average over the spatial slice at z ≈ 0. This results in a
higher expansion in the observer’s vicinity relative to the “global” expansion, which presents
as smaller luminosity distances locally for z . 0.02 in figure 9.
Figure 10 shows the all-sky average luminosity distance, 〈DL〉sky = 1/NLOS∑iDiL,
relative to the EdS distance, as a function of redshift (mean z across the slice) for all 25
observers. This corresponds to the equivalent of the thick solid black curve in figure 9 for all
observers. The dark and light grey shaded regions again show ±2% and ±5% variation from
– 25 –
J
C
A
P03(2023)019
10 2 10 1
zmean
0.2
0.1
0.0
0.1
0.2
0.3
D
L
/d
L(
E
dS
)
1
Mean over whole sky
Mean over 100 LOS shown
o 0.33
Directions of LOS shown
Figure 9. Deviation in ray-traced luminosity distance from the EdS prediction as a function of
redshift. We show 100 lines of sight for the observer shown in figure 8 (coloured dashed curves), as well
as the mean over the whole sky (thick solid black curve) and the mean over the 100 lines of sight shown
(thick dashed black curve). The inset shows the directions of the lines of sight shown on the observer’s
sky, and the dark and light grey shaded regions show ±2% and ±5% deviations, respectively.
EdS, respectively. Each line segment is coloured according to the local density environment
nearby that observer for that redshift scale, δsmooth. Specifically, for each value of zmean we
use a Gaussian kernel to smooth the density field over the corresponding spatial scale, namely,
δsmooth ≡
∫
δ(xi)×W (xi, σ) dV, (8.3)
where δ(xi) is the density contrast defined on the z ≈ 0 simulation hypersurface, and
W (xi, σ) is a normalised Gaussian kernel centred on the observer’s position with width
σ = dC,EdS(zmean)/3. Here, dC,EdS(zmean) is the co-moving distance in the EdS model at
redshift zmean (the redshift at which 〈DL〉sky is calculated). We emphasize again that we use
the density field on the z ≈ 0 simulation hypersurface, so this will give an upper limit of the
density as smoothed on the observer’s past light cone. Further, we approximate the volume
element in the average (8.3) as dV ≈ d3X instead of the more general dV = √γ d3X where γ
is the determinant of the spatial metric of the hypersurface. This approximation is sufficient
to gain a rough insight to the local density contrast near the observer. We emphasize that
this smoothing method is used only for the purposes of visualisation in figure 10.
The inset in figure 10 shows a zoom in of the redshift range 0.05–0.3 such that the
absolute value of the smoothed density contrast (shown on a log scale) can be seen. We
expect this all-sky average to closely coincide with the monopole (isotropic) contribution
to the luminosity distance, since our distribution of “objects” fairly samples all points on
the sky (to within our HEALPix resolution). This kind of variance is sometimes referred to
as “cosmic variance”, and represents the variation we might expect in distances as we place
– 26 –
J
C
A
P03(2023)019
10 2 10 1
zmean
0.1
0.0
0.1
0.2
0.3
D
L
sk
y
/d
L(
E
dS
)
1
0.4
0.2
0.0
0.2
0.4
0.6
0.8
sm
oo
th
0.05 0.10 0.15 0.20 0.25 0.30
0.010
0.005
0.000
0.005
0.010
10 3
10 2
10 1
|
sm
oo
th
|
Figure 10. Mean luminosity distance relative to the EdS prediction as a function of redshift. We show
the all-sky average DL for 25 observers in the simulation shown in figure 7. Each curve is coloured
according to that observer’s local density contrast as smoothed to the scale corresponding to zmean,
namely, δsmooth. The dark and light shaded regions show 2% and 5% deviations from EdS, respectively.
The inset shows a zoom-in of the redshift range zmean = 0.05–0.3, with |δsmooth| on a log scale.
observers in regions with different local density contrast. In figure 10, we do notice some
correlation between the smoothed local density contrast and the deviation in DL at z . 0.01.
This correlation is no longer obviously present beyond z ≈ 0.02 and there is no discernible
correlation with variance in DL and δsmooth by z ≈ 0.05.
Aside from cosmic variance, the anisotropic variance of DL across the observer’s sky may
bias their cosmological inference if they do not fairly sample all directions. To estimate the
level of anisotropy across the skies of our 25 observers, we calculate the maximal sky-deviation
of the luminosity distance (in a manner similar to MH21),
∆DL ≡ D
max
L −DminL
2 , (8.4)
where DmaxL and DminL are the maximum and minimum values of DL on an observer’s sky,
respectively. Figure 11 shows the maximal sky variance (8.4) for all 25 observers in this
simulation. We show ∆DL relative to the EdS distance as a function of the mean redshift of
the slice. ∆DL here is computed using the uniform redshift slice of DL, namely, all lines of
sight are linearly interpolated to zmean (though zmean varies among observers). All observers
maintain a consistent maximal sky deviation of > 8% until z ≈ 0.03, which then remains > 5%
until z ≈ 0.06, and finally the variance is still above 1% out to the maximum redshift studied
here, z ≈ 0.3. At z ≈ 0.3 the range of variance we see is 1.2–2% across all 25 observers. In
appendix B we compare this result to the same calculation including a boost of the observers
to the “CMB frame” of the simulation (see also section 6.5).
The level of variance we see in figure 11 could bias cosmological inferences using data
which does not fairly sample the whole sky. The precise size of any effect on the resulting
– 27 –
J
C
A
P03(2023)019
10 2 10 1
zmean
10 2
10 1
100
D
L
/d
L(
E
dS
)
1
10%
1%
Figure 11. Maximal sky variance as a function of redshift for 25 observers in the simulation shown
in figure 7. We show ∆DL relative to the EdS distance at that redshift, dL(EdS). Dashed and dotted
horizontal lines show a 10% and 1% deviation, respectively.
cosmological parameters will be dependent on the particular survey geometry or sky coverage
in question as well as the redshift ranges considered (see MH21 for an approximate study of
the potential impact of this effect). We delay an in-depth investigation of the impact of these
effects — for example on low-redshift supernova data — to future work.
9 Caveats
The main caveats to our work are as follows:
1. Our simulations are matter dominated and contain no cosmological constant. The
simulations maintain a matter density Ωm ≈ 1 throughout the simulation on large
scales [i.e. there is no significant backreaction on the scale of the simulation domain,
see 66], and the structures that form will therefore grow at a faster rate. Ultimately,
this results in overall higher density contrasts in the simulation with respect to an
equivalent simulation containing a cosmological constant. In a model universe with
typically lower density contrasts, i.e. in ΛCDM, we would expect the amplitude of
anisotropic signatures to be reduced, however, the qualitative signatures will be robust
to these kinds of changes (see also MH21 for a discussion on this). Any universe with
cosmological structure will naturally contain these effects, however, the amplitude is
dependent on the particular cosmology in question.
2. All simulations we use here sample only large-scale structures, due to the limitation of
the continuous fluid approximation. We expect the anisotropic signatures to be different
with the inclusion of additional small-scale structure. In Macpherson and Heinesen [89,
in prep.], we investigate the impact of small-scale structure on the large-scale angular
features in the DL(z) relation by comparing the two simulations presented in section 8.2
and 8.3.
– 28 –
J
C
A
P03(2023)019
3. We have studied a limited number of observers in this work due to the computational
expense of performing all-sky ray tracing. While our observers span a wide range of
local density contrasts δo ≈ −0.4–1 (see figure 10), we have not nearly sampled all
possible local environments, nor are these environments especially like ours as observers
on Earth. Especially of relevance would be the study of targeted observers with similar
local environments to what we observe in galaxy surveys, rather than a sample of
randomly-placed observers. Our choice for using randomly-placed observers here was
motivated on gaining a broad understanding of the variance in these effects throughout
the simulation. This method certainly has benefits over studying only one observer,
however, the large variance we find among observers should further motivate studies of
these effects in constrained simulations of the local Universe [90, 91].
10 Conclusions
In this work we have presented an analysis of anisotropy in the luminosity distance at low
redshift with ray tracing. Our framework is based on fully nonlinear general relativity, with
no assumption of any background cosmology or perturbation theory at any point in the
calculation (beyond setting initial data for the simulations themselves).
We assessed the ability of the general cosmography [as truncated at third order in
redshift 37] to capture anisotropies in the luminosity distance as calculated using our new
ray-tracing tool. We found the cosmographic distance provides a < 1% accurate prediction
out to redshift z ≈ 0.1 for strict smoothing scales of 200h−1 Mpc, and to redshift z ≈ 0.034 for
strict smoothing scales of 100h−1 Mpc. We note that in the latter case, the smoothing scale
and the maximum redshift closely coincide with one another in terms of distance scale in an
EdS model approximation. In the former case, the maximum redshift represents slightly larger
scales than the smoothing scale in the EdS model. We therefore expect the cosmographic
prediction to be best at redshift scales closely coinciding with the smoothing scale itself, which
(as also stated in appendix A of MH21) may pose an issue for using the general cosmography,
as truncated at third order, for wide redshift intervals.
For simulations with small-scale structure beneath ∼ 200h−1 Mpc removed from the
initial data — rather than being strictly excluded for the whole simulation — we find that
the cosmography is accurate to within ∼ 1% to a redshift of z ≈ 0.012. While we cannot
explicitly state the size of the “smoothing scale” of this simulation at redshift zero, it is clear
that the presence of additional small-scale structure spoils the ability of the cosmography
(as expanded to third order) to capture DL to higher redshift. This breakdown might be
expected [38], and therefore some smoothing of the small-scale structure will likely be required
to yield an accurate large-scale prediction of anisotropies from the general cosmography as
truncated at third order [89, in prep.]. All simulations studied here are still very far from the
real Universe. The minimum scale we sample is ∼8–10 h−1 Mpc, which is the scale of large
galaxy clusters and many structures exist below this scale.
Therefore, we conclude that achieving a ∼ 1% accurate prediction of the luminosity
distance (including all of its anisotropic contributions) will be difficult for wide redshift ranges.
As mentioned in section 3.1, such issues can be alleviated by allowing for a decay of an
anisotropic signature included in a cosmological fit [such as in, e.g., 53–55, 57]. This would
ensure that the anisotropies themselves are being constrained by data within the narrow
redshift interval in which we expect the series to converge, while higher-redshift objects are
contributing to the background cosmological fit. While we cannot smooth the actual structure
– 29 –
J
C
A
P03(2023)019
in the Universe, as we can in simulations, some smoothing of the cosmological data itself may
be possible to allow anisotropic constraints [92, in prep.].
We further assessed the level of anisotropy in the luminosity distance in a simulation
including structures down to ∼ 100h−1 Mpc in the initial data. We found all 25 observers
saw a > 8% variance in DL across their sky to z ≈ 0.03 and > 1% variance at z ≈ 0.3. This
kind of variance could be especially important for analysis assuming an isotropic distance
law that does not fairly sample the full sky. However, further investigation is necessary to
estimate the size of the effect for specific survey geometries and redshift ranges. We also found
the “cosmic variance” of the ray-traced luminosity distance for z . 0.01 is ±10% (relative
to the EdS distance) with some correlation with the local density at the observer. Above
z ≈ 0.03, however, all 25 observers have a < 2% difference in their all-sky average (isotropic
contribution) relative to the EdS model. This “cosmic variance” is roughly consistent with
previous studies on the isotropic variance in the local Hubble parameter using Newtonian
N-body simulations [93–95] and NR simulations [67].
Acknowledgments
HJM would like to thank the anonymous referee whose careful reading and comments improved
the quality of this work. HJM would also like to thank Julian Adamek, Asta Heinesen,
Pierre Fleury, Jim Mertens, Sofie Koksbang, and Eoin Ó Colgáin for helpful discussions and
comments regarding this manuscript. HJM appreciates support received by the Herchel Smith
postdoctoral fellowship fund. HJM would also like to especially thank Pierre Fleury, Julian
Adamek, and Jim Mertens for very useful discussions on light propagation over the years that
inspired and formed the beginning of this work. This work used the DiRAC@Durham facility
managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and
STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. We
used GNU parallel in processing data for this work [96].
A Propagation of the screen vectors
In section 6.2.1 we introduced the notion of screen space, which is defined by the basis vectors
sµA. Not all of the components of s
µ
A are constrained by (6.14), with the orientation of the
basis remaining arbitrary. A family of observers lying along the geodesic co-moving with
uµ therefore may define their own basis vectors which satisfy the constraints but are not
necessarily aligned with one another. In this work, we will be considering a single basis per
observer per geodesic, and we therefore must maintain the same orientation of the basis along
the geodesic (although the initial orientation remains arbitrary, and our final results are
independent of this choice, which we show in appendix A.3 below). The screen vectors are
(by definition) orthogonal to the observer’s 4-velocity uµ, which is not necessarily parallel
transported along the geodesic. The basis vectors therefore obey a partial parallel transport
equation [see section 1.3.3 of 69, for a proof of this], which for convenience can be written
as [81]
DsµA
dλ
= k
µ
E
sαA
Duα
dλ
, (A.1)
where the covariant total derivative D/dλ ≡ kµ∇µ. The beams morphology itself is inde-
pendent of any observer, and so it can be shown that the evolution of the screen has no
– 30 –
J
C
A
P03(2023)019
physical effect on the beam [75]. It can therefore be safely assumed that the 4-velocity uµ is
parallel transported along the geodesic and thus so are the screen vectors. After making this
assumption, the screen basis is a “Sachs basis” and this can be useful in situations where a
clear matter rest-frame is not well-defined [e.g. 76]. In this work, we continue working with a
general uµ and thus partially parallel transport the screen vectors along the geodesic.
Similarly to our re-casting of the geodesic equation in section 6.1, we can re-parametrise
the geodesic in terms of coordinate time, t. After expanding the covariant total derivative
in (A.1), we can thus write
dsµA
dt
= k
µ
k0E
sαAk
δ∇δuα − 1
k0
Γµβγk
βsγA. (A.2)
From (A.2), we find the evolution equations for the time and space components of the screen
vectors to be
ds0A
dt
= s
0
A
E
{
k0∂tu0 + ki∂iu0 +
1
α
∂t(α)kiui − Evi∂iα+ αKkiukki
}
+ s
i
A
E
{
k0∂tui + kj∂jui +
1
α
∂i(α)kjuj + EKkivk − kjΓkijuk
}
, (A.3)
and
dsiA
dt
= s
0
A
E
{
ki∂tu0 +
kikj
k0
∂ju0 + Γki∂tα+
kikj
k0
αukKkj − αEγij∂jα+ αE
k0
Kijk
j
}
+ s
m
A
E
{
ki∂tum +
kikj
k0
∂jum + Γki∂mα− k
ikj
k0
Γkjmuk + αEKim −
E
k0
Γijmkj
}
, (A.4)
respectively. The derivatives of the 4-velocity, uµ, in the above can be expanded in terms of
derivatives of the velocity vi = ui/(αu0) as follows
∂tu0 =
u0Γ2
2 ∂t(v
ivi)− Γ∂tα, (A.5)
∂lu0 =
u0Γ2
2 ∂l(v
ivi)− Γ∂lα, (A.6)
∂tuj = Γγij∂tvi +
Γ3
2 γijv
i∂t(vkvk) + 2u0viKij , (A.7)
∂luj = Γγij∂lvi +
Γ3
2 γijv
i∂l(vkvk) + Γvi∂lγij , (A.8)
where Γ = 1/
√
1− vivi is the Lorentz factor, we have used u0 = −αΓ, and we have
∂t(vivi) = −2αKijvivj + 2γijvi∂tvj , (A.9)
∂l(vivi) = vivj∂lγij + 2γijvi∂lvj . (A.10)
For computational convenience, the evolution equations actually solved in mescaline
are written in terms of some quantities already calculated in other routines, namely
ds0A
dt
= s
0
A
E
{
u0Γ2
2
(
k0∂tv
2 + ki∂iv2
)
− E
α
∂tα− ∂iα
(
Γki + Evi
)
+ αKijkiuj
}
+ s
m
A
E
{
k0Γγim∂tvi +
k0Γ3
2 γimv
i∂tv
2 +Kimvi
(
2k0u0 + E
)
(A.11)
+ Γγlmki∂ivl +
Γ3
2 γlmv
lki∂iv
2 + Γkivl∂iγlm +
1
α
uik
i∂mα− kiΓjimuj
}
,
– 31 –
J
C
A
P03(2023)019
and
dsiA
dt
= s
0
A
E
{
u0Γ2
2
(
ki∂tv
2 + k
ikj
k0
∂jv
2
)
+ k
ikj
k0
(
αKkju
k−Γ∂jα
)
−αEγij∂jα+ αE
k0
Kijk
j
}
+ s
m
A
E
{
Γkiγlm
(
∂tv
l+ k
j
k0
∂jv
l
)
+ Γ
3
2 k
iγlmv
l
(
∂tv
2 + k
j
k0
∂jv
2
)
(A.12)
+2u0kivlKlm+Γki∂mα+αEγikKkm+
kikj
k0
Γvl∂jγlm− k
ikj
k0
Γkjmuk−
E
k0
Γijmkj
}
,
where ∂µv2 = ∂µ(vivi).
A.1 Observers at rest
We can show that s0A = 0 for any observer at rest with respect to the coordinates in the case
of zero shift. First, expanding the constraint that sµA must be orthogonal to the observer
4-velocity, i.e. gµνsµAuν = 0, we have
g00s
0
Au
0 + g0is0Aui + gi0siAu0 + gijsiAuj = 0. (A.13)
An observer at rest has 4-velocity uµrest = (1/α,0) in general, namely ui = 0 (which further
implies Γ = 1). In the above, we then have the remaining non-zero terms
g00s
0
Au
0 + gi0siAu0 = 0, (A.14)
and for zero shift we have gi0 = 0, implying g00s0Au0 = 0, and since either g00 or u0 cannot be
zero, we must have s0A = 0 for the constraint s
µ
Auµ to be satisfied.
A.2 Initial screen vectors
To calculate the initial values of the screen vectors, we begin with the constraints
sµAuµ = s
µ
Aeµ = gµνs
µ
As
ν
B = 0, (A.15)
which represent five equations for a total of eight unknowns (sµA). We can therefore freely
choose three components of sµA (i.e. the orientation of the basis) and solve the constraints for
the remaining five. We then must ensure our basis is normalised, i.e. that the following is
satisfied:
gµνs
µ
As
ν
A = gµνs
µ
Bs
ν
B = 1. (A.16)
We take our solutions to (A.15) and normalise them to unit length. This results in a set of
basis vectors which satisfies all of the required constraints.
A.3 Test of initial screen vectors
The screen vectors are an arbitrary set of basis vectors to describe the screen space, and any
set satisfying the constraints (for the same observer in the same environment) should give the
same physical result. To test this, we rotate the screen vectors by making a different choice
for our free three degrees of freedom and check the physical results are unchanged.
First, we perform a test run for linearly-perturbed EdS spacetime with perturbation
amplitude φ0 = 10−5 and resolution N = 48 (see appendix E.3 for further details on this test
– 32 –
J
C
A
P03(2023)019
4
2
0
2
D
L
/D
L,
E
dS
1
(×
10
4 )
init sA 1
init sA 2
Heun, N = 48
1.0
0.5
0.0
0.5
1.0
s0 1
(×
10
6 )
0.0
0.5
1.0
s1 1
(×
10
2 )
1
0
1
s2 1
(×
10
2 )
10 2 10 1 100
zEdS
2
1
0
1
2
s3 1
(×
10
2 )
Figure 12. Top panel: luminosity distance as a function of redshift in a linearly-perturbed EdS
metric for two test runs with different initial screen vectors sµA (dashed and solid curves). Bottom
four panels: screen vector components for the same two test runs (matching colours). Each of the ten
curves show an individual line of sight. The observer and all lines of sight are identical in both tests.
The screen vectors change (bottom four panels) but the physical result (DL, top panel) is the same.
setup). We run the same test twice with all parameters identical except for a different initial
set of screen vectors. Figure 12 shows the luminosity distance perturbation as a function of
redshift for the two tests with different sets of initial screen vectors (solid purple and dashed
cyan curves). In the bottom four panels we show each component of the screen vectors (for
A = 1) for the same tests, in the same colours. We can see the screen vector components
differ between the tests but the physical result, DL in the top panel, is consistent.
– 33 –
J
C
A
P03(2023)019
8
6
4
2
0
2
4
D
L
/D
L,
E
dS
1
(×
10
3 ) init sA 1
init sA 2
Heun, N = 64
4
2
0
2
4
s0 1
(×
10
6 )
5
0
5
s1 1
(×
10
4 )
10
5
0
5
s2 1
(×
10
4 )
10 2 10 1
zEdS
10
5
0
5
10
s3 1
(×
10
4 )
Figure 13. Top panel shows the ray-traced luminosity distance in an ET simulation sampling only
large-scale structures for two runs with different initial screen vectors (dashed and solid lines). Each of
the ten curves show a randomly drawn line of sight for the same observer. Bottom four panels show
the screen vector components for the same pair of runs.
We also perform this test with simulation output, using an ET simulation with resolution
N = 64, box size L = 640h−1 Mpc and all power beneath 100h−1 Mpc removed from the
initial data. We perform the ray-tracing to z ≈ 0.1 for two different sets of initial screen
vectors. The top panel of figure 13 shows the ray-traced luminosity distance DL, normalised
by the EdS distance, as a function of redshift for two different sets of initial screen vectors
(dashed cyan and solid purple curves). The bottom four panels show each of the screen vector
components for A = 1. We see the same result that the screen vectors differ between the two
runs but the physical result, DL in the top panel, is consistent.
– 34 –
J
C
A
P03(2023)019
We thus conclude from these tests that our choice of initial screen vectors, chosen via
the three arbitrary degrees of freedom, has no effect on our physical results.
B Boosting data to the hypersurface frame
In section 6.5, we introduced the possibility of approximating the CMB frame — as used in
observational cosmology — using the simulation hypersurfaces themselves. In this case, our
observers are still co-moving with the fluid flow, however, after measuring their distances
and redshifts in that co-moving frame, they boost their observations to the “CMB frame”.
Usually, such an observer would infer their velocity with respect to this frame using the
dipole they measure in the CMB [as in 86]. For this section, we instead assume the observer’s
velocity with respect to the fictitious “CMB frame” is their peculiar velocity with respect to
the Eulerian coordinates, namely, vµ. If we assume the “CMB frame” is well described by
an EdS expansion law, an observer moving with velocity vµ relative to this frame will see a
dipole modulation of their measured luminosity distance of the form [see, e.g. 41, 97]
dobsL (z) = dL,EdS(z) + vµeµ
(1 + z)
H(z) , (B.1)
where eµ is the direction of observation on the sky andH(z) is the conformal Hubble parameter.
Peculiar motion of the observer with respect to an FLRW frame induces a dipole in both the
measured redshift, z, and the luminosity distance, dL. However, both of these effects have
been accounted for in deriving the above expression [see 97, for more details], and thus dL(z)
and dL,EdS(z) are calculated at the same (observed) redshift z.
To apply a boost to each observer’s data in our simulations, we calculate dobsL above
using the peculiar velocity, vi, of the fluid at their location on the grid. We then subtract
the dipole component of dobsL dipole from the ray-traced luminosity distance DL to “correct”
for the motion of the observer. We confirm this process removes the majority of the dipole
contribution to the luminosity distance, however, some dipolar modulation still remains.
We would also like to emphasize that this process is not exactly the same as that done
with real observational data. In our observations, we measure the redshift but do not usually
have direct access to the luminosity distance. Peculiar velocity corrections are thus usually
applied to our redshift measurements. However, here we are working with simulated DL(z)
which has been interpolated along each line of sight, in order to arrive at a surface of constant
z. The intention of this type of data is to directly probe the DL(z) relation with respect to
the FLRW theory prediction. For this type of data, applying the boost to the luminosity
distance, using the method described above, allows us to estimate the impact of the observer’s
peculiar velocity on the anisotropy measured in DL. In future studies, we aim to work with
synthetic observational catalogues instead of shells of constant z, wherein we will consider
corrections which may be more in line with observational data analysis methods.
We might expect this correction to reduce the amount of anisotropy for our observers.
In figure 14 we compare the results for the maximal sky variance in DL calculated in the fluid
frame (as shown in figure 11, reproduced in the left panel) to the boosted observations in the
“CMB frame” (right panel). We see a reduction in the spread of variance across observers,
however, no significant reduction in the maximal sky variance is seen. We note again that
this “CMB frame” is not necessarily coincident with the way in which we define the CMB
frame in observational cosmology. The correct way to define such a frame would be to infer
the velocity with respect to the CMB that each observer would see. However, as mentioned in
– 35 –
J
C
A
P03(2023)019
10 2 10 1
zmean
10 2
10 1
100
D
L
/d
L(
E
dS
)
1
10%
1%
Fluid frame
10 2 10 1
zmean
10 2
10 1
100
D
L
/d
L(
E
dS
)
1
10%
1%
Boosted to 'CMB' frame
Figure 14. Maximal sky variance as a function of redshift for 25 observers in the simulation shown
in figure 7. We show ∆DL relative to the EdS distance at that redshift, dL(EdS). Left panel shows
the variance for observations in the frame co-moving with the fluid flow and the right panel shows
the same observations boosted to the “CMB frame” using each observers coordinate peculiar velocity.
Dashed and dotted horizontal lines in both panels show 10% and 1% deviation, respectively.
the main text, this would require ray tracing to very high redshift in order to obtain a map
to mimic the CMB radiation.
We note that the peculiar velocites of the “sources” (with respect to the simulation
hypersurface) considered here may be able to account for the anisotropies we see in the right
panel of figure 14. However, since we find that the simulation slice remains very close to the
background EdS model we choose for the initial slice [when averaging over large scales, see
66] this is perhaps not surprising. Usually, in cosmological data analysis we might also apply
some corrections for the peculiar velocity of low-redshift sources with respect to the CMB
frame [see 85]. Investigating the effect of such corrections on our simulation data falls beyond
the scope of this work. Connecting the observational CMB frame to ray-traced simulation
data is an important avenue to pursue if we wish to make predictions for future observations.
C Spatial interpolation
Mescaline reads in HDF5 data at a distinct set of time slices as output by the ET simulation.
The code itself is thus tracking the propagation of light beams through this simulated space-
time, at the specific coordinate times defined by the constant-t surfaces. The simulation is
also discretised in space and in general the light beam will not be positioned on a grid cell
at each time step; it will land somewhere in between. At every step, we therefore must use
spatial interpolation to calculate the simulation quantities (e.g. metric tensor, curvature)
appearing on the right-hand-side of the evolution equations at the position of the light beam.
For this purpose, we use the B-spline Fortran library,8 which provides one to six dimensional
B-spline interpolation on a regular grid. In mescaline, the order of the spatial interpolation
is left as a free variable. However, for all purposes in this paper we use fifth-order spline
interpolation. We have found that a higher-order interpolation offers very little improvement
in the accuracy of the results. Calculating all source terms at the fifth-order stencil is an
expensive part of the mescaline ray tracer, however, reducing to fourth order interpolation
only yields a ∼ 4% increase in speed.
8http://jacobwilliams.github.io/bspline-fortran.
– 36 –
J
C
A
P03(2023)019
D Transformation of vectors from local Fermi frame
Quantities output from the raytracer, such as redshift and luminosity distance, are calculated
as seen by observers moving with 4-velocity uµ. In this work, we take this to coincide with
the 4-velocity of the fluid flow. However, initial direction vectors we specify (eµ) are defined
in our coordinate system, which is not comoving with the fluid. We will generally be given
the directions of objects in the observer’s frame, e.g., the right ascension and declination
from a given survey. These coordinates are represented by the Cartesian direction vector
miˆ. We must transform these to the direction vector in our coordinate system, mi, which
is then used to calculate eµ. When plotting sky-distributions of redshift or distance we use
miˆ to ensure angular power spectra, direction of multipoles, etc., are correctly defined in the
observer’s frame.
In this section, we describe the transformation from a vector defined in the observer’s
local Fermi frame9 — also known as the local inertial frame, or proper reference frame of an
observer who is not necessarily geodesic [see pg. 327 of 98] — to the simulation coordinate
system. We consider the observer’s metric to be conformally Minkowski, with local “scale
factor” equal to γ1/3 — since the observer’s “scale factor” may not necessarily follow the
convention a(z = 0) = 1 within the simulation. We transform the local metric at the observer,
η˜αβ = γ−1/3ηαβ, to the metric tensor in our coordinate system using the following basis
(where αˆ represents the observer’s coordinate xαˆ)
gµν = eµαˆe
ν
βˆ
gαˆβˆ, (D.1)
= eµαˆe
ν
βˆ
γ−1/3ηαˆβˆ, (D.2)
such that
γ1/3gµν = eµ0ˆe
ν
0ˆη
0ˆ0ˆ + eµ
iˆ
eν
jˆ
ηiˆjˆ , (D.3)
= −uµuν + eµ
iˆ
eν
jˆ
δiˆjˆ , (D.4)
where in the last line we have used eµ0ˆ = u
µ, which is the observer’s time coordinate in our
coordinate system, i.e. the 4-velocity of the observer.
The transformation of miˆ is given by
mi = ei
jˆ
mjˆ , (D.5)
where we do not need to transform the time component since e0 is set via the normalisation
condition of eµ in mescaline. We therefore only need to solve six equations, namely
g˜ij ≡ γ1/3gij = −uiuj + eikejl δkl, (D.6)
(where we have now removed the hats for simplification) for the nine components of eij . We
therefore have three degrees of freedom remaining, and without loss of generality we can set
some components to zero:
eij =
exx 0 0exy eyy 0
exz e
y
z e
z
z
 . (D.7)
9The author would like to thank Jim Mertens for correspondence aiding with the calculations in this section.
– 37 –
J
C
A
P03(2023)019
From this choice, we find the remaining components are
ezz =
√
g˜zz + (uz)2, (D.8)
eyz =
g˜yz + uyuz
ezz
, (D.9)
eyy =
√
g˜yy + (uy)2 − (eyz)2, (D.10)
exz =
g˜xz + uxuz
ezz
, (D.11)
exy =
g˜xy + uxuy − exzeyz
eyy
, (D.12)
exx =
√
g˜xx + (ux)2 − (exy)2 − (exz )2. (D.13)
This transformation is implemented after determining the directions mjˆ but before calculating
the direction vector eµ, which is defined in the simulation coordinate system.
E Tests of the mescaline ray tracer
We must test the accuracy and precision of the mescaline ray tracer on known analytic
solutions. Tests such as these are necessary in order for us to trust results from simulation
output, where there is not necessarily a known analytic solution. Here we present the results
of ray tracing in two test metrics. In appendix E.2 we test within a pure EdS space-time: a
flat FLRW metric tensor with a matter-dominated source (no dark energy). In appendix E.3,
we test within a linearly-perturbed EdS space-time.
For all tests we use an RK2 (Heun) integration with Courant factor dt/dx= 0.1 and a
cubic domain with L = 4h−1 Gpc and three resolutions N = 12, 24, and 48. In all tests, we
feed an analytic form of the metric into mescaline before passing this into the regular (and
unchanged) ray-tracing routines. In the following sections, we will present the analytic metric
for both cases and the analytic solutions for the quantities we evolve, namely kµ, sµA, and DA
(and z in the case of pure EdS). We always compare three resolutions to the analytic solutions
to ensure the error is converging at the expected rate with a decrease in time step size.
E.1 Definition of error and convergence
We define the relative error with respect to the relevant analytic solutions, e.g., for the
luminosity distance in an EdS space-time, using
err(dL) ≡ dL,num
dL,EdS
− 1, (E.1)
where dL,num is the numerical luminosity distance and dL,EdS is the analytic luminosity
distance.
We calculate the rate of convergence by evaluating the error at all three resolutions with
grid spacing ∆x1, ∆x2, and ∆x3. For an integration method of order p, the error will be
O(∆xp), and the convergence rate of the error in (E.1) is given by
C = err∆x1 − err∆x2err∆x2 − err∆x3
, (E.2)
– 38 –
J
C
A
P03(2023)019
and the expected convergence rate is
Cexp = ∆x
p
1 −∆xp2
∆xp2 −∆xp3
. (E.3)
For the RK2 method we have p = 2 and therefore Cexp = 4 for the case of ∆x1 = 2∆x2 = 4∆x3.
E.2 EdS spacetime
For our first test we will use the EdS space-time, which has metric tensor
ds2 = −a2dη2 + a2δijdxidxj , (E.4)
where η is the conformal time, and a = ainit(η/ηinit)2 is the scale factor, and ainit and ηinit
are the initial values of the scale factor and conformal time, respectively.
E.2.1 Analytic solutions
We adopt the convention ainit = 1 for the tests, implying a(z = 0) = 1 + zinit. The scale factor
and redshift are then simply related by
z(η) = ainit
a(η) − 1. (E.5)
The luminosity distance in an EdS spacetime, with cosmological density parameters Ωm =
1,ΩΛ = Ωk = 0, is
dL,EdS(z) =
2
H0
(
1 + z −√1 + z
)
, (E.6)
where H(z) ≡ a′/a is the conformal Hubble parameter, where a dash is a derivative with
respect to conformal time, and H0 ≡ H(z = 0). Substituting the EdS metric into the evolution
equations (6.10) for k0 and ki we find
dk0
dt
= −2k0H, dk
i
dt
= −2kiH, (E.7)
which gives
k0 = k
0
ini
a2
, ki = k
i
ini
a2
(E.8)
for some arbitrary initial data, k0ini and kiini. Similarly, for the propagation of the screen
vectors (A.11) and (A.12) we find
dsiA
dt
= −siAH, (E.9)
and therefore
siA =
siA,ini
a
, (E.10)
for some arbitrary initial siA,ini. Since our observers are co-moving with the EdS matter
content, they have 4-velocity components ui = 0 and thus s0A = 0 and ds0A/dt = 0 (see A.1).
– 39 –
J
C
A
P03(2023)019
10 5
10 4
z/
z E
dS
1
N = 12
N = 24
N = 48
Heun
10 2 10 1 100
zEdS
0.975
1.000
1.025
C
/C
ex
p
Figure 15. Mescaline redshift, z, relative to the EdS redshift, zEdS, as a function of zEdS. We show
the relative error for three resolutions N = 12, 24, and 48 in the top panel, and the convergence rate,
C, relative to the expected (Cexp = 4 for RK2) in the bottom panel.
10 7
10 6
10 5
10 4
d L
/d
L,
E
dS
1
N = 12
N = 24
N = 48
10 2 10 1 100
zEdS
0.95
1.00
1.05
C
/C
ex
p
Figure 16. Ray-traced luminosity distance, dL, in an analytic EdS metric relative to the EdS analytic
solution, dL,EdS, as a function of zEdS. We show the relative error for three resolutions N = 12, 24,
and 48 in the top panel, and the convergence rate, C, relative to the expected value for RK2 in the
bottom panel.
– 40 –
J
C
A
P03(2023)019
10 7
10 6
10 5
10 4
k
/k
E
dS
1
k0
k1
k2
k3
N = 12
N = 24
N = 48
10 8
10 7
10 6
10 5
s A
/s
A,
E
dS
1
s11
s21
s31
s12
s22
s32
10 2 10 1 100
zEdS
0.99
1.00
1.01
C
/C
ex
p
10 2 10 1 100
zEdS
0.995
1.000
C
/C
ex
p
Figure 17. Photon 4-momentum, kµ, and screen vectors, sµA, as calculated with mescaline relative
to their analytic EdS evolution. The top left panel shows the relative error in each component of kµ
(different colours) for three resolutions (different line styles), and the bottom left panel shows the
convergence rate for all kµ components, relative to the expected rate for RK2. The top right panel
shows the relative error in siA for three resolutions, with the bottom right panel showing the respective
convergence rates.
E.2.2 Mescaline vs. analytic solutions
We set up the pure EdS metric (E.4) and send this into mescaline, using the same routines
which calculate z and DA for a general metric evolved using NR. We compare the redshift
z, luminosity distance dL, photon 4-momentum kµ, and screen vectors sµA to the analytic
solutions given in the previous section. The initial data for both kµ and sµA (k
µ
ini and s
µ
ini) is
dependent on the direction of observation and therefore is arbitrary. We also check the null
condition and the constraints on both eµ and sµA are satisfied to within round-off error after
generating initial data. We follow one line of sight only, since all directions are equivalent in
the EdS space-time, for each resolution until z ≈ 1.3.
Figure 15 shows the error in the redshift as calculated with mescaline relative to the
EdS solution (E.5) calculated using the conformal time. We show three resolutions (different
colours) in the top panel, and the respective convergence rate (E.2) in the bottom panel
relative to the expected Cexp = 4 for RK2. For our highest (yet still coarse) resolution of
N = 48 the error remains below 0.01%, and the convergence rate is as expected.
Figure 16 shows the error in the luminosity distance as calculated with mescaline for an
analytic EdS metric relative to the EdS solution. The top panel shows the relative error for
each resolution, and the bottom panel shows the convergence rate for the three curves, again
relative to the expected rate for RK2. Error in the luminosity distance lies below ∼ 0.005%
for N = 48 at redshift zEdS ∼ 1, and the error converges as expected.
Figure 17 shows the error in the photon 4-momentum, kµ, and the screen vectors, sµA,
with respect to their EdS analytic solutions (E.8) and (E.10), respectively, as a function of
EdS redshift. The top left panel of figure 17 shows the relative error in each component of
kµ (different colours) for three resolutions (different line styles), and the bottom left panel
– 41 –
J
C
A
P03(2023)019
shows the convergence rate for each set of three resolutions as different point styles (which all
closely coincide). The top and bottom right panels shows relative error and convergence for
the components siA. We note again that s0A = 0 for EdS, which we have checked is the case in
these tests. The error converges at the expected rate in all cases shown here.
E.3 Linearly-perturbed EdS spacetime
We now test the mescaline ray tracer in the presence of inhomogeneities. For this we will use
a linearly-perturbed EdS space-time, which (for scalar-only perturbations) has metric tensor
ds2 = −a2(1 + 2φ)dη2 + a2(1− 2φ)δijdxidxj , (E.11)
where φ = φ(xi) is a small perturbation with amplitude φ0  1. Since our observers are
co-moving with the fluid flow, they also must be perturbed with respect to the background.
For zero shift, the components of the fluid 4-velocity are ui = Γvi and u0 = Γ/α, where
Γ = 1/
√
1− vivi is the Lorentz factor and vi is the peculiar velocity of the fluid with respect
to the Eulerian observer. In the limit of linear perturbations, we therefore have ui = vi and
u0 = 1/α. The observer 4-velocity is thus uµ = (1/α, vi), where the velocity perturbation is
given in Macpherson et al. [66].
In the next section, we will derive the geodesic equation and the Jacobi evolution
equations for the linearly-perturbed EdS metric. We will then numerically solve this system
using an RK4 integrator to give a semi-analytic solution to compare with the mescaline
results. We use RK4 such that the analytic solution is closer to the true solution and thus
we can still check if our mescaline solutions (advanced with RK2) are converging at the
expected rate.
E.3.1 Semi-analytic solutions
The geodesic equation in the perturbed space-time with metric (E.11) is
dk0
dt
= −2Hk0 − 2k
j∂jφ
1 + 2φ , (E.12)
dki
dt
= −2Hki − k
0∂iφ
1− 2φ −
kk
k0
Γikjkj , (E.13)
which is correct to all orders in φ. The Christoffel symbols are
Γijk =
1
(1− 2φ)
[
δil∂l(φ)δjk − ∂j(φ)δik − ∂k(φ)δij
]
, (E.14)
which is also correct to all orders. The position of the photon is calculated via (6.8c), namely,
dxi/dt = ki/k0 with ki and k0 given as solutions to (E.12) and (E.13), respectively.
For the perturbed EdS metric (E.11), the evolution equations for the components of the
screen vectors sµA, to linear order in φ and all of its derivatives, become
ds0A
dt
= s
0
A
E
(
−Hk0α− aki∂iφ
)
+ s
m
A
E
[
a′δimvi
(
2ak0 − E)+ k0a2δim∂tvi + a2δlmki∂ivl],
(E.15)
– 42 –
J
C
A
P03(2023)019
and
dsiA
dt
= s
0
A
E
[
− ak
ikj
k0
(a′vkδkj + ∂jφ)− E∂iφ− EH k
i
k0
]
(E.16)
+ s
m
A
E
[
a2kiδlm
(
kj
k0
∂jv
l + ∂tvl
)
+ 2aa′kivlδlm + aki∂mφ− EHδim −
E
k0
Γijmkj
]
.
The Jacobi matrix equation is, also to linear order in φ and all of its derivatives,
d
dt
PAB =
(
2H+ 2
k0
ki∂iφ
)
PAB +
1
(k0)2R
A
CDCB, (E.17)
where we have used the null condition to substitute
δijk
ikj
(k0)2 =
α2
a2(1− 2φ) =
1 + 2φ
1− 2φ. (E.18)
We also must build the optical tidal matrix RAC for the linearly-perturbed space-time
to advance (E.17). The components of the Riemann tensor we need, correct to first order in
φ and all of its derivatives, are
R0i0j =
[
(a′)2 − aa′′
]
(1− 2φ)δij + a2∂i∂jφ,
R0xxy = −aa′∂yφ,
R0xxz = −aa′∂zφ,
R0yxy = aa′∂xφ,
R0yyz = −aa′∂zφ,
R0zxz = aa′∂xφ,
R0zyz = aa′∂yφ,
Rxyxz = a2∂y∂zφ,
Rxyxy = (a′)2(1− 6φ) + a2
(
∂2yφ+ ∂2xφ
)
,
Rxyyz = −a2∂x∂zφ,
Rxzyz = a2∂x∂yφ,
Rxzxz = (a′)2(1− 6φ) + a2
(
∂2zφ+ ∂2xφ
)
,
Ryzyz = (a′)2(1− 6φ) + a2
(
∂2zφ+ ∂2yφ
)
,
with R0xyz = R0yxz = R0zxy = 0. The components of the Ricci tensor are
R00 = −3H′ +∇2φ, (E.20)
R0i = 2H∂iφ, (E.21)
Rij =
[(
H2 + a
′′
a
)
(1− 4φ) +∇2φ
]
δij . (E.22)
In deriving all of the equations in this section, we have made use of the following relations:
∂iα =
a2∂iφ
α
, ∂tα =
a′α
a
, vj∂lγij = −2a2vj∂l(φ)δij .
– 43 –
J
C
A
P03(2023)019
10 5
10 4
D
A
/D
A,
(s
em
i
an
a)
1
N = 24
N = 48
N = 96
Heun
10 2 10 1 100
zEdS
0.995
1.000
1.005
C
/C
ex
p
Figure 18. Top panel: Relative difference between ray-traced DA and the semi-analytic solution for a
linearly-perturbed EdS spacetime. We show the average error over ten lines of sight for three resolutions
N = 24, 48, and 96 as a function of the background redshift, zEdS. Bottom panel: convergence rate C,
relative to the expected rate Cexp = 4 for RK2, for the three curves in the top panel.
E.3.2 Mescaline vs. linearly-perturbed EdS
We set up the linearly-perturbed EdS metric (E.11) in mescaline and send it into the ray
tracer. We use a functional form of φ(xi) of
φ(xi) = φ0
∑
i
sin
(
2pixi
L
)
, (E.23)
where L is the box length of the test domain and corresponds to the wavelength of the
perturbation. For the analytic solutions, we use the same metric and advance kµ, xµ, sµA,
and DAB according to the evolution equations given in the previous section. We store the
semi-analytic solutions for the angular diameter distance, photon 4-momentum, and screen
vectors at the same time steps as the ray-tracer output. We use three resolutions N = 24, 48,
and 96, a perturbation amplitude of φ0 = 10−8 to ensure second-order terms can be safely
neglected, and a 4 h−1 Gpc box size so that we can shoot rays to z ≈ 1 without sampling the
structure more than once. We place an observer at the centre of the domain, and shoot ten
randomly drawn lines of sight. We calculate the relative difference between the fully nonlinear
mescaline results (e.g. DA) and the semi-analytic solutions (e.g. DA,(semi−ana)) for each line
of sight, and average the resulting error over all lines of sight.
The top panel of figure 18 shows the ray-traced angular diameter distance relative to
the semi-analytic distance for the same linearly-perturbed EdS metric. Different colours show
resolutions N = 24, 48, and 96 and the bottom panel shows the convergence rate C for the
three curves relative to the expected for RK2 of Cexp = 4. The error is always below ∼ 0.003%
for resolution N = 96, and converges at the expected rate.
Figure 19 shows the relative difference between the photon 4-momentum components kµ
(top left) and screen vector components sµ1 (top right) and the semi-analytic solutions for the
– 44 –
J
C
A
P03(2023)019
10 7
10 6
10 5
10 4
k
/k
(s
em
i
an
a)
1
k0
k1
k2
k3
N = 24
N = 48
N = 96
10 8
10 7
10 6
10 5
10 4
10 3
10 2
s A
/s
A,
(s
em
i
an
a)
1
Heun
s01
s11
s21
s31
10 2 10 1 100
zEdS
0.95
1.00
1.05
C
/C
ex
p
k0
k1
k2
k3
10 2 10 1 100
zEdS
1.0
1.1
1.2
C
/C
ex
p
s01
s11
s21
s31
Figure 19. Top panels: Relative difference between kµ (left) and sµ1 (right) calculated with mescaline
and the semi-analytic solutions for a linearly-perturbed EdS space-time. We show the average difference
over ten lines of sight for three resolutions N = 24, 48, and 96 as a function of the background redshift.
Bottom panels: convergence rate C, relative to the expected Cexp = 4, for the curves in the respective
top panels.
linearly-perturbed EdS metric. The bottom panels show the convergence rate C relative to
the expected Cexp = 4 for RK2 for the respective top panels. All kµ and si1 components have
very small errors and converge at the expected rate. We note that we only show sµ1 (and not
sµ2 ) because both cases A = 1, 2 are evolved using the same source terms.
The error in s01 is a few orders of magnitude larger than the spatial components. For
the background we have s0A = 0, and so this component is perturbative only and thus has a
smaller magnitude than the spatial components, all of which have background contributions.
This further implies that a larger fraction of the error is made up of round-off error, and so
we might expect the relative error to be larger than the other components. Additionally, the
noise we see in the error in s01 can be attributed to a sign change in the error in some of the
lines of sight that make up the average shown here.
E.4 Errors in analysis of simulation data
We are also interested in the level of error in our ray-tracing calculations in a situation for which
we do not have an analytic solution. In such a case, we can estimate the error in our calculation
via a Richardson extrapolation. This is based on the expectation that our calculations should
approach the “true” solution as we increase our numerical resolution towards inifinty at a
rate which is determined by the order of accuracy of the implemented scheme.
We use three NR simulations with resolutions N = 64, 128, and 256 which are initialised
with a power spectrum of large-scale only perturbations within a L = 1.28h−1 Gpc domain.
Specifically, all power beneath ∼ 200h−1 Mpc is removed from the initial data, and thus
the simulations remain close to the linear regime for the extent of the matter-dominated
simulations from the initial redshift of z = 1000 until z = 0. The initial data is identical
between resolutions, and thus we expect minimal difference between the simulations at redshift
– 45 –
J
C
A
P03(2023)019
0.02 0.04 0.06 0.08 0.10 0.12
zmean
0.0
0.2
0.4
0.6
0.8
1.0
er
r(
D
L)
(×
10
5 )
Obs. mean
Figure 20. Sky-average of the relative error in DL for 10 observers (coloured dashed curves) as a
function of the mean redshift of the slice. The thick black curve shows the average over all observers.
zero for the large scales we sample. These simulations are thus suitable for a Richardson
extrapolation (for quantities calculated at individual positions in the domain), for which we
require derivatives to be the same between resolutions.
For all simulations, we place 10 observers randomly and propagate an isotropic sky of
12×N2side lines of sight for Nside = 16 until z ≈ 0.1 using an RK2 integrator. All lines of sight
and observer positions are identical for all simulations. We thus expect DL and z calculated
along each line of sight, and for each observer, to converge at a rate ∝ N−2. For each observer,
each iteration, and each line of sight, we fit a curve of the form DL(N) = DL,inf + X/N2,
where DL,inf is the N →∞ luminosity distance (or the “true” solution) and X is a constant.
We determine both DL,inf and X using the curve_fit function as a part of the SciPy10
Python package, using a similar relation for the redshift z. The relative error in DL (similarly
for z) is then
err(DL) ≡ DL,N=256
DL,inf
− 1, (E.24)
namely, the relative difference between our highest-resolution calculation and the “true”
solution at N →∞. We average the error over each observer’s sky at each iteration to track
the level of error for all observers at all redshifts, namely 〈err(DL)〉Ω(zmean) for each observer.
Figures 20 and 21 show the sky-averaged relative error in DL and z, respectively, as
a function of the mean redshift on that observer’s slice, zmean. The relative error in DL is
always of order 10−5 and the relative error in z is always of order 10−4.
E.5 Violation of the null constraint
In mescaline, the propagation of the components of the photon 4-momentum kµ, as well as
the setting of initial data, is done independently of the null condition requirement for photons,
10https://scipy.org.
– 46 –
J
C
A
P03(2023)019
0.02 0.04 0.06 0.08 0.10 0.12
zmean
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
er
r(
z)
(×
10
4 )
Obs. mean
Figure 21. Sky-average of the relative error in the redshift z for 10 observers (coloured dashed
curves) as a function of the mean redshift of the slice. The thick black curve shows the average over
all observers.
namely kµkµ = 0. We can thus use this condition as a numerical test of the accuracy of the
numerical evolution of kµ.
In this appendix, we assess the violation of the null condition in the mescaline ray
tracer. First, we test this using the analytic linear-perturbed metric tensor from appendix E.3,
and second we test this using a set of NR simulations similar to the one used in appendix A.3.
E.5.1 Linearly-perturbed metric
We use the test functionality of mescaline to pass the linearly-perturbed FLRW metric
tensor (E.11) into the ray tracing routines. For each case, we place a single observer in
the model universe and propagate 10 randomly-drawn lines of sight for a set time interval.
We repeat this process for three resolutions N = 12, 24, and 48 with the same spacing
in time (resulting in different numbers of iterations). We also use four different cases of
perturbation size, controlled by the amplitude of the metric perturbation φ0. We choose
φ0 = 10−2, 10−4, 10−6, and 10−8 for these tests.
Figure 22 shows the absolute value of the null condition violation, |kµkµ|, as a function
of redshift for four cases of φ0 (panels). Different colours represent different resolutions, as
indicated in the legend, and solid curves show the mean over all 10 lines of sight for this
observer and the shaded region represents the range from minimum violation to maximum
violation across all 10 lines of sight. All panels show the same observer, and the redshift on
the x-axis is the redshift averaged over all lines of sight. We note that for φ0 = 10−2, the final
redshift is z ≈ 0.27 whereas all other perturbation sizes give a final redshift of z ≈ 0.06 for the
same simulation time, due to the larger gravitational field the photon must travel through.
The violation in the null condition reduces as we increase resolution, as would be expected
of a numerical error. The ratio |kµkµ|/φ0 is approximately constant at ≈ 10−6–10−5 for
N = 48 across all cases. Since the violation reduces with resolution, we conclude that it is
dominated by numerical error.
– 47 –
J
C
A
P03(2023)019
10 2 10 1
zmean
10 10
10 9
10 8
10 7
10 6
10 5
10 4
|k
k
|
0 = 10 2
N = 12
N = 24
N = 48
10 3 10 2
zmean
10 12
10 11
10 10
10 9
10 8
10 7
10 6
|k
k
|
0 = 10 4
10 3 10 2
zmean
10 15
10 14
10 13
10 12
10 11
10 10
10 9
10 8
|k
k
|
0 = 10 6
10 3 10 2
zmean
10 16
10 15
10 14
10 13
10 12
10 11
10 10
|k
k
|
0 = 10 8
Figure 22. Absolute value of the violation of the null condition as a function of redshift (mean across
10 lines of sight) for a linearly-perturbed FLRW test metric tensor. Each panel shows a different
amplitude of perturbation, all with three resolutions N = 12, 24, and 48 (magenta, blue, and orange
curves, respectively). Solid curves show the mean over 10 lines of sight for one observer, and the
shaded regions show the range of minimum violation to maximum violation.
E.5.2 Simulation test
We also test the violation of the null condition in a set of NR simulations. In this case, we
randomly place 10 observers within the simulation z ≈ 0 hypersurfaces, and propagate 50
randomly-drawn lines of sight for each observer until z ≈ 0.06–0.09 using an RK2 integrator,
tracking the violation in kµkµ along the way. We use two simulations, with resolution N = 64
and 128, which have identical initial data and box size 640h−1 Mpc. The initial perturbations
are drawn from a CLASS power spectrum and have all modes below 100h−1 Mpc removed
to ensure modes are sampled with sufficiently many grid points. This set of simulations are
similar to those used in appendix E.4 except they sample smaller-scale structures initially.
Therefore, the physical structure in these simulations at z ≈ 0 differs between resolutions even
though the initial data is identical. We maintain the same spacing in iteration for the runs at
different resolution, to ensure we are studying convergence of both spatial and time derivatives.
– 48 –
J
C
A
P03(2023)019
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
zmean
10 6
10 5
|k
k
|
N = 64
N = 128
Figure 23. Absolute value of the violation of the null condition as a function of redshift for 10
observers placed in two NR simulations. Solid and dashed curves show the violation for 10 observers
as averaged over 50 lines of sight for resolutions N = 128 and 64, respectively. We show the violation
as a function of each observer’s mean redshift for that time slice, zmean.
Figure 23 shows the absolute value of the violation in the null condition for 10 observers
averaged across 50 lines of sight each for resolutions N = 128 (solid curves) and N = 64
(dashed curves). The violation in the null condition remains < 10−5 for all cases studied here,
and reduces by about a factor of four with increasing resolution N = 64→ 128, as is expected
based on a second-order accurate scheme.
References
[1] Planck collaboration, Planck 2018 results. VII. Isotropy and Statistics of the CMB, Astron.
Astrophys. 641 (2020) A7 [arXiv:1906.02552] [INSPIRE].
[2] L. Perivolaropoulos and F. Skara, Challenges for ΛCDM: An update, New Astron. Rev. 95 (2022)
101659 [arXiv:2105.05208] [INSPIRE].
[3] D.W. Hogg et al., Cosmic homogeneity demonstrated with luminous red galaxies, Astrophys. J.
624 (2005) 54 [astro-ph/0411197] [INSPIRE].
[4] M. Scrimgeour et al., The WiggleZ Dark Energy Survey: the transition to large-scale cosmic
homogeneity, Mon. Not. Roy. Astron. Soc. 425 (2012) 116 [arXiv:1205.6812] [INSPIRE].
[5] R.G. Clowes et al., A structure in the early universe at z ∼ 1.3 that exceeds the homogeneity
scale of the R-W concordance cosmology, Mon. Not. Roy. Astron. Soc. 429 (2013) 2910
[arXiv:1211.6256] [INSPIRE].
[6] I. Horváth, Z. Bagoly, J. Hakkila and L.V. Tóth, New data support the existence of the
Hercules-Corona Borealis Great Wall, Astron. Astrophys. 584 (2015) A48 [arXiv:1510.01933]
[INSPIRE].
[7] S. Sarkar, B. Pandey and R. Khatri, Testing isotropy in the Universe using photometric and
spectroscopic data from the SDSS, Mon. Not. Roy. Astron. Soc. 483 (2019) 2453
[arXiv:1810.07410] [INSPIRE].
– 49 –
J
C
A
P03(2023)019
[8] C.A.P. Bengaly, A. Bernui, J.S. Alcaniz and I.S. Ferreira, Probing cosmological isotropy with
Planck Sunyaev-Zeldovich galaxy clusters, Mon. Not. Roy. Astron. Soc. 466 (2017) 2799
[arXiv:1511.09414] [INSPIRE].
[9] D. Alonso et al., Homogeneity and isotropy in the Two Micron All Sky Survey Photometric
Redshift catalogue, Mon. Not. Roy. Astron. Soc. 449 (2015) 670 [arXiv:1412.5151] [INSPIRE].
[10] C. Marinoni, J. Bel and A. Buzzi, The Scale of Cosmic Isotropy, JCAP 10 (2012) 036
[arXiv:1205.3309] [INSPIRE].
[11] C. Gibelyou and D. Huterer, Dipoles in the Sky, Mon. Not. Roy. Astron. Soc. 427 (2012) 1994
[arXiv:1205.6476] [INSPIRE].
[12] C.M. Hirata, Constraints on cosmic hemispherical power anomalies from quasars, JCAP 09
(2009) 011 [arXiv:0907.0703] [INSPIRE].
[13] O. Luongo et al., Larger H0 values in the CMB dipole direction, Phys. Rev. D 105 (2022) 103510
[arXiv:2108.13228] [INSPIRE].
[14] N.J. Secrest et al., A Test of the Cosmological Principle with Quasars, Astrophys. J. Lett. 908
(2021) L51 [arXiv:2009.14826] [INSPIRE].
[15] K. Migkas et al., Cosmological implications of the anisotropy of ten galaxy cluster scaling
relations, Astron. Astrophys. 649 (2021) A151 [arXiv:2103.13904] [INSPIRE].
[16] P. Tiwari et al., Dipole anisotropy in sky brightness and source count distribution in radio NVSS
data, Astropart. Phys. 61 (2014) 1 [arXiv:1307.1947] [INSPIRE].
[17] S. Appleby and A. Shafieloo, Testing Isotropy in the Local Universe, JCAP 10 (2014) 070
[arXiv:1405.4595] [INSPIRE].
[18] M. Rubart and D.J. Schwarz, Cosmic radio dipole from NVSS and WENSS, Astron. Astrophys.
555 (2013) A117 [arXiv:1301.5559] [INSPIRE].
[19] A.K. Singal, Large peculiar motion of the solar system from the dipole anisotropy in sky
brightness due to distant radio sources, Astrophys. J. Lett. 742 (2011) L23 [arXiv:1110.6260]
[INSPIRE].
[20] P.K. Aluri et al., Is the Observable Universe Consistent with the Cosmological Principle?,
arXiv:2207.05765 [INSPIRE].
[21] T. Buchert, On average properties of inhomogeneous fluids in general relativity. 1. Dust
cosmologies, Gen. Rel. Grav. 32 (2000) 105 [gr-qc/9906015] [INSPIRE].
[22] S. Räsänen, Backreaction: directions of progress, Class. Quant. Grav. 28 (2011) 164008
[arXiv:1102.0408] [INSPIRE].
[23] T. Buchert, P. Mourier and X. Roy, On average properties of inhomogeneous fluids in general
relativity III: general fluid cosmologies, Gen. Rel. Grav. 52 (2020) 27 [arXiv:1912.04213]
[INSPIRE].
[24] T. Buchert, H. van Elst and A. Heinesen, The averaging problem on the past null cone in
inhomogeneous dust cosmologies, Gen. Rel. Grav. 55 (2023) 7 [arXiv:2202.10798] [INSPIRE].
[25] D.J. Schwarz, C.J. Copi, D. Huterer and G.D. Starkman, CMB Anomalies after Planck, Class.
Quant. Grav. 33 (2016) 184001 [arXiv:1510.07929] [INSPIRE].
[26] C. Clarkson and R. Maartens, Inhomogeneity and the foundations of concordance cosmology,
Class. Quant. Grav. 27 (2010) 124008 [arXiv:1005.2165] [INSPIRE].
[27] M. Visser, Jerk and the cosmological equation of state, Class. Quant. Grav. 21 (2004) 2603
[gr-qc/0309109] [INSPIRE].
– 50 –
J
C
A
P03(2023)019
[28] A.G. Riess et al., A Comprehensive Measurement of the Local Value of the Hubble Constant with
1 km s−1 Mpc−1 Uncertainty from the Hubble Space Telescope and the SH0ES Team, Astrophys.
J. Lett. 934 (2022) L7 [arXiv:2112.04510] [INSPIRE].
[29] W.L. Freedman et al., The Carnegie-Chicago Hubble Program. VIII. An Independent
Determination of the Hubble Constant Based on the Tip of the Red Giant Branch, Astrophys. J.
882 (2019) 34 [arXiv:1907.05922] [INSPIRE].
[30] J. Adamek, D. Daverio, R. Durrer and M. Kunz, General Relativistic N -body simulations in the
weak field limit, Phys. Rev. D 88 (2013) 103527 [arXiv:1308.6524] [INSPIRE].
[31] J. Adamek, D. Daverio, R. Durrer and M. Kunz, gevolution: a cosmological N-body code based on
General Relativity, JCAP 07 (2016) 053 [arXiv:1604.06065] [INSPIRE].
[32] J.T. Giblin, J.B. Mertens and G.D. Starkman, Departures from the
Friedmann-Lemaitre-Robertston-Walker Cosmological Model in an Inhomogeneous Universe: A
Numerical Examination, Phys. Rev. Lett. 116 (2016) 251301 [arXiv:1511.01105] [INSPIRE].
[33] E. Bentivegna and M. Bruni, Effects of nonlinear inhomogeneity on the cosmic expansion with
numerical relativity, Phys. Rev. Lett. 116 (2016) 251302 [arXiv:1511.05124] [INSPIRE].
[34] H.J. Macpherson, P.D. Lasky and D.J. Price, Inhomogeneous Cosmology with Numerical
Relativity, Phys. Rev. D 95 (2017) 064028 [arXiv:1611.05447] [INSPIRE].
[35] W.E. East, R. Wojtak and T. Abel, Comparing Fully General Relativistic and Newtonian
Calculations of Structure Formation, Phys. Rev. D 97 (2018) 043509 [arXiv:1711.06681]
[INSPIRE].
[36] D. Daverio, Y. Dirian and E. Mitsou, General relativistic cosmological N-body simulations. Part
I. Time integration, JCAP 10 (2019) 065 [arXiv:1904.07841] [INSPIRE].
[37] A. Heinesen, Multipole decomposition of the general luminosity distance ‘Hubble law’ — a new
framework for observational cosmology, JCAP 05 (2021) 008 [arXiv:2010.06534] [INSPIRE].
[38] H.J. Macpherson and A. Heinesen, Luminosity distance and anisotropic sky-sampling at low
redshifts: A numerical relativity study, Phys. Rev. D 104 (2021) 023525 [Erratum ibid. 104
(2021) 109901] [arXiv:2103.11918] [INSPIRE].
[39] N. Sugiura, N. Sugiyama and M. Sasaki, Anisotropies in Luminosity Distance, Prog. Theor. Phys.
101 (1999) 903.
[40] C. Bonvin, R. Durrer and M.A. Gasparini, Fluctuations of the luminosity distance, Phys. Rev. D
73 (2006) 023523 [Erratum ibid. 85 (2012) 029901] [astro-ph/0511183] [INSPIRE].
[41] C. Bonvin, R. Durrer and M. Kunz, The dipole of the luminosity distance: a direct measure of
H(z), Phys. Rev. Lett. 96 (2006) 191302 [astro-ph/0603240] [INSPIRE].
[42] I. Ben-Dayan, G. Marozzi, F. Nugier and G. Veneziano, The second-order luminosity-redshift
relation in a generic inhomogeneous cosmology, JCAP 11 (2012) 045 [arXiv:1209.4326]
[INSPIRE].
[43] J. Adamek et al., Bias and scatter in the Hubble diagram from cosmological large-scale structure,
Phys. Rev. D 100 (2019) 021301 [arXiv:1812.04336] [INSPIRE].
[44] J.T. Giblin, J.B. Mertens and G.D. Starkman, Observable Deviations from Homogeneity in an
Inhomogeneous Universe, Astrophys. J. 833 (2016) 247 [arXiv:1608.04403] [INSPIRE].
[45] J.T. Giblin, J.B. Mertens, G.D. Starkman and A.R. Zentner, General Relativistic Corrections to
the Weak Lensing Convergence Power Spectrum, Phys. Rev. D 96 (2017) 103530
[arXiv:1707.06640] [INSPIRE].
[46] C. Cattoën and M. Visser, The Hubble series: Convergence properties and redshift variables,
Class. Quant. Grav. 24 (2007) 5985 [arXiv:0710.1887] [INSPIRE].
– 51 –
J
C
A
P03(2023)019
[47] C. Gruber and O. Luongo, Cosmographic analysis of the equation of state of the universe through
Padé approximations, Phys. Rev. D 89 (2014) 103506 [arXiv:1309.3215] [INSPIRE].
[48] T. Yang, A. Banerjee and E.Ó. Colgáin, Cosmography and flat ΛCDM tensions at high redshift,
Phys. Rev. D 102 (2020) 123532 [arXiv:1911.01681] [INSPIRE].
[49] E. Lusso et al., Tension with the flat ΛCDM model from a high-redshift Hubble diagram of
supernovae, quasars, and gamma-ray bursts, Astron. Astrophys. 628 (2019) L4
[arXiv:1907.07692] [INSPIRE].
[50] S. Capozziello, Ruchika and A.A. Sen, Model independent constraints on dark energy evolution
from low-redshift observations, Mon. Not. Roy. Astron. Soc. 484 (2019) 4484
[arXiv:1806.03943] [INSPIRE].
[51] A. Banerjee et al., On problems with cosmography in cosmic dark ages, Phys. Lett. B 818 (2021)
136366 [arXiv:2009.04109] [INSPIRE].
[52] A. Heinesen and H.J. Macpherson, A prediction for anisotropies in the nearby Hubble flow,
JCAP 03 (2022) 057 [arXiv:2111.14423] [INSPIRE].
[53] W. Rahman et al., New constraints on anisotropic expansion from supernovae Type Ia, Mon.
Not. Roy. Astron. Soc. 514 (2022) 139 [arXiv:2108.12497] [INSPIRE].
[54] J. Colin, R. Mohayaee, M. Rameez and S. Sarkar, A response to Rubin & Heitlauf: “Is the
expansion of the universe accelerating? All signs still point to yes”, arXiv:1912.04257
[INSPIRE].
[55] D. Rubin and J. Heitlauf, Is the expansion of the universe accelerating? All signs still point to
yes a local dipole anisotropy cannot explain dark energy, Astrophys. J. 894 (2020) 68
[arXiv:1912.02191] [INSPIRE].
[56] J. Colin, R. Mohayaee, M. Rameez and S. Sarkar, Evidence for anisotropy of cosmic acceleration,
Astron. Astrophys. 631 (2019) L13 [arXiv:1808.04597] [INSPIRE].
[57] S. Dhawan, A. Borderies, H.J. Macpherson and A. Heinesen, The quadrupole in the local Hubble
parameter: first constraints using Type Ia supernova data and forecasts for future surveys, Mon.
Not. Roy. Astron. Soc. 519 (2023) 4841 [arXiv:2205.12692] [INSPIRE].
[58] D. Brout et al., The Pantheon+ Analysis: Cosmological Constraints, Astrophys. J. 938 (2022)
110 [arXiv:2202.04077] [INSPIRE].
[59] C. Tian et al., Question of measuring spatial curvature in an inhomogeneous universe, Phys. Rev.
D 103 (2021) 083513 [arXiv:2010.07274] [INSPIRE].
[60] F. Löffler et al., The Einstein Toolkit: A Community Computational Infrastructure for
Relativistic Astrophysics, Class. Quant. Grav. 29 (2012) 115001 [arXiv:1111.3344] [INSPIRE].
[61] M. Zilhão and F. Löffler, An Introduction to the Einstein Toolkit, Int. J. Mod. Phys. A 28 (2013)
1340014 [arXiv:1305.5299] [INSPIRE].
[62] E. Bentivegna, An automatically generated code for relativistic inhomogeneous cosmologies, Phys.
Rev. D 95 (2017) 044046 [arXiv:1610.05198] [INSPIRE].
[63] K. Wang, Numerical Relativity Investigation of the Effects of Gravitational Waves on the
Inhomogeneity of the Universe, Eur. Phys. J. C 78 (2018) 629 [arXiv:1801.08362] [INSPIRE].
[64] D. Blas, J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS)
II: Approximation schemes, JCAP 07 (2011) 034 [arXiv:1104.2933] [INSPIRE].
[65] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview,
arXiv:1104.2932 [INSPIRE].
[66] H.J. Macpherson, D.J. Price and P.D. Lasky, Einstein’s Universe: Cosmological structure
formation in numerical relativity, Phys. Rev. D 99 (2019) 063522 [arXiv:1807.01711] [INSPIRE].
– 52 –
J
C
A
P03(2023)019
[67] H.J. Macpherson, P.D. Lasky and D.J. Price, The trouble with Hubble: Local versus global
expansion rates in inhomogeneous cosmological simulations with numerical relativity, Astrophys.
J. Lett. 865 (2018) L4 [arXiv:1807.01714] [INSPIRE].
[68] G.F.R. Ellis, On the definition of distance in general relativity: I. M. H. Etherington
(Philosophical Magazine ser. 7, vol. 15, 761 (1933)), Gen. Rel. Grav. 39 (2006) 1047.
[69] P. Fleury, Light propagation in inhomogeneous and anisotropic cosmologies, Ph.D. thesis, Paris
University, VI, IAP, France (2015) [arXiv:1511.03702] [INSPIRE].
[70] C.J. Fluke, R.L. Webster and D.J. Mortlock, The ray bundle method for calculating weak
magnification by gravitational lenses, Mon. Not. Roy. Astron. Soc. 306 (1999) 567
[astro-ph/9812300] [INSPIRE].
[71] C.J. Fluke and P.D. Lasky, Shape, shear and flexion II — Quantifying the flexion formalism for
extended sources with the ray-bundle method, Mon. Not. Roy. Astron. Soc. 416 (2011) 1616
[arXiv:1101.4407] [INSPIRE].
[72] E. Bentivegna, M. Korzyński, I. Hinder and D. Gerlicher, Light propagation through black-hole
lattices, JCAP 03 (2017) 014 [arXiv:1611.09275] [INSPIRE].
[73] S.A. Ema, M.R. Hossen, K. Bolejko and G.F. Lewis, The density distributions of cosmic
structures: impact of the local environment on weak-lensing convergence, Mon. Not. Roy. Astron.
Soc. 509 (2021) 3004 [arXiv:2110.14089] [INSPIRE].
[74] M.-A. Breton and V. Reverdy, Magrathea-Pathfinder: a 3D adaptive-mesh code for geodesic ray
tracing in N-body simulations, Astron. Astrophys. 662 (2022) A114 [arXiv:2111.08744]
[INSPIRE].
[75] R.K. Sachs, Gravitational waves in general relativity. 6. The outgoing radiation condition, Proc.
Roy. Soc. Lond. A 264 (1961) 309 [INSPIRE].
[76] F. Lepori et al., Weak-lensing observables in relativistic N-body simulations, Mon. Not. Roy.
Astron. Soc. 497 (2020) 2078 [arXiv:2002.04024] [INSPIRE].
[77] M. Grasso, E. Villa, M. Korzyński and S. Matarrese, Isolating nonlinearities of light propagation
in inhomogeneous cosmologies, Phys. Rev. D 104 (2021) 043508 [arXiv:2105.04552] [INSPIRE].
[78] M. Grasso and E. Villa, BiGONLight: light propagation with bilocal operators in numerical
relativity, Class. Quant. Grav. 39 (2022) 015011 [arXiv:2107.06306] [INSPIRE].
[79] M. Grasso, M. Korzyński and J. Serbenta, Geometric optics in general relativity using bilocal
operators, Phys. Rev. D 99 (2019) 064038 [arXiv:1811.10284] [INSPIRE].
[80] K.M. Górski et al., HEALPix — A framework for high resolution discretization, and fast analysis
of data distributed on the sphere, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [INSPIRE].
[81] V.A.A. Sanghai, P. Fleury and T. Clifton, Ray tracing and Hubble diagrams in post-Newtonian
cosmology, JCAP 07 (2017) 028 [arXiv:1705.02328] [INSPIRE].
[82] S. Räsänen, Light propagation in statistically homogeneous and isotropic dust universes, JCAP
02 (2009) 011 [arXiv:0812.2872] [INSPIRE].
[83] N. Li and D.J. Schwarz, On the onset of cosmological backreaction, Phys. Rev. D 76 (2007)
083011 [gr-qc/0702043] [INSPIRE].
[84] E.W. Kolb, S. Matarrese and A. Riotto, On cosmic acceleration without dark energy, New J.
Phys. 8 (2006) 322 [astro-ph/0506534] [INSPIRE].
[85] E.R. Peterson et al., The Pantheon+ Analysis: Evaluating Peculiar Velocity Corrections in
Cosmological Analyses with Nearby Type Ia Supernovae, Astrophys. J. 938 (2022) 112
[arXiv:2110.03487] [INSPIRE].
[86] Planck collaboration, Planck 2013 results. XXVII. Doppler boosting of the CMB: Eppur si
muove, Astron. Astrophys. 571 (2014) A27 [arXiv:1303.5087] [INSPIRE].
– 53 –
J
C
A
P03(2023)019
[87] T.M. Siewert, M. Schmidt-Rubart and D.J. Schwarz, Cosmic radio dipole: Estimators and
frequency dependence, Astron. Astrophys. 653 (2021) A9 [arXiv:2010.08366] [INSPIRE].
[88] C. Dalang and C. Bonvin, On the kinematic cosmic dipole tension, Mon. Not. Roy. Astron. Soc.
512 (2022) 3895 [arXiv:2111.03616] [INSPIRE].
[89] H.J. Macpherson and A. Heinesen, Ray tracing in full general relativity and all-sky map of
cosmic distances, in prepaparation (2023).
[90] S. Heß, F.-S. Kitaura and S. Gottloeber, Simulating Structure Formation of the Local Universe,
Mon. Not. Roy. Astron. Soc. 435 (2013) 2065 [arXiv:1304.6565] [INSPIRE].
[91] S. Heß and F.-S. Kitaura, Cosmic flows and the expansion of the local Universe from non-linear
phase-space reconstructions, Mon. Not. Roy. Astron. Soc. 456 (2016) 4247 [arXiv:1412.7310]
[INSPIRE].
[92] J. Adamek, C. Clarkson, R. Durrer, A. Heinesen, M. Kunz and H.J. Macpherson, Towards
Cosmography of the Local Universe, in preparation (2023).
[93] H.-Y. Wu and D. Huterer, Sample variance in the local measurements of the Hubble constant,
Mon. Not. Roy. Astron. Soc. 471 (2017) 4946 [arXiv:1706.09723] [INSPIRE].
[94] I. Odderskov, S.M. Koksbang and S. Hannestad, The local value of H0 in an inhomogeneous
universe, JCAP 02 (2016) 001 [arXiv:1601.07356] [INSPIRE].
[95] R. Wojtak et al., Cosmic variance of the local Hubble flow in large-scale cosmological simulations,
Mon. Not. Roy. Astron. Soc. 438 (2014) 1805 [arXiv:1312.0276] [INSPIRE].
[96] O. Tange, Gnu Parallel 2018, Zenodo (2018) [DOI:10.5281/zenodo.1146014].
[97] L. Hui and P.B. Greene, Correlated Fluctuations in Luminosity Distance and the (Surprising)
Importance of Peculiar Motion in Supernova Surveys, Phys. Rev. D 73 (2006) 123526
[astro-ph/0512159] [INSPIRE].
[98] C.W. Misner, K.S. Thorne and J.A. Wheeler, Gravitation, Princeton University Press (2017).
– 54 –