# Towards Cosmography of the Local Universe

Julian Adamek  Institut für Astrophysik, Universität Zürich,
Winterthurerstrasse 190, 8057 Zürich, Switzerland  Chris Clarkson  School of
Physical and Chemical Sciences, Queen Mary University of London,
327 Mile End Road, London E1 4NS, UK
Department of Physics & Astronomy, University of the Western Cape,
Cape Town 7535, South Africa  Ruth Durrer  Département de Physique Théorique,
Université de Genève,
24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland  Asta Heinesen  Université
de Lyon 1 & ENS de Lyon & CNRS & Centre de Recherche Astrophysique de Lyon
UMR5574,
69007 Lyon, France
Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark  Martin Kunz  Département de Physique
Théorique, Université de Genève,
24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland  Hayley J. Macpherson
Kavli Institute for Cosmological Physics, The University of Chicago,
5640 South Ellis Avenue, Chicago, Illinois 60637, USA
NASA Einstein Fellow

###### Abstract

Anisotropies in the distance–redshift relation of cosmological sources are
expected due to large-scale inhomogeneities in the local Universe. When the
observed sources are tracing a large-scale matter flow in a general spacetime
geometry, the distance–redshift relation with its anisotropies can be
described with a geometrical prediction that generalises the well-known
Friedmann-Lemaître-Robertson-Walker result. Furthermore, it turns out that a
finite set of multipole coefficients contain the full information about a
finite-order truncation of the distance–redshift relation of a given observer.
The multipoles of the distance–redshift relation are interesting new
cosmological observables that have a direct physical interpretation in terms
of kinematical quantities of the underlying matter flow. Using light cones
extracted from $N$-body simulations we quantify the anisotropies expected in a
$\Lambda$ cold dark matter cosmology by running a Markov chain Monte Carlo
analysis on the observed data. In this observational approach the survey
selection implements an implicit smoothing scale over which the effective rest
frame of matter is fitted. The perceived anisotropy therefore depends
significantly on the redshift range and distribution of sources. We find that
the multipoles of the expansion rate, as well as the observer’s velocity with
respect to the large-scale matter flow, can be determined robustly with our
approach.

###### keywords:

Cosmology, Local Universe, $N$-body Simulations

##  1 Introduction

In modern cosmology, data are often analysed within a particular field theory
and with assumptions about the energy-momentum content, for instance general
relativity with ordinary matter, cold dark matter, radiation, and dark energy.
A complementary way to analyse data is through cosmography, where the field
theory of gravity is left unspecified. Cosmography thus presents an
independent approach to constraining the kinematics and the curvature of the
Universe, without relying on a particular dynamical theory (Weinberg, 1972;
Ellis et al., 1985).

Cosmographic approaches are particularly suitable for low-redshift
measurements, but can also be applied in the intermediate- to high-redshift
regime (see e.g. Cattoen & Visser, 2007; Vitagliano et al., 2010; Capozziello
et al., 2020). The low-redshift regime of cosmology is sensitive to the large-
scale structures within our cosmic neighbourhood. The anisotropies that are
expected in cosmological observables due to cosmic structures are commonly
modelled using peculiar velocity corrections as predicted by a $\Lambda$ cold
dark matter ($\Lambda$CDM) cosmology. This procedure has the disadvantage that
one must trust that the derived corrections indeed give a complete description
of the anisotropies in the observations. An alternative approach is to
incorporate the anisotropies into the cosmography itself, thereby fitting for
the anisotropies instead of making a model prediction for them. Such
cosmographic treatments can determine the recent expansion history of the
Universe in a model-independent way (Ellis et al., 1985; Clarkson et al.,
2012). Moreover, estimates of the impact of local structure on the
measurements on various scales might provide insights into cosmological
tensions and anomalies seen across a growing number complementary cosmological
data sets (Perivolaropoulos & Skara, 2022; Abdalla et al., 2022).

To formulate a cosmography that is valid for describing the anisotropy in our
cosmological measurements, we must go beyond the spatially maximally-symmetric
Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes and formulate
observables geometrically in the absence of metric symmetries. Such
formulations were pioneered by Kristian & Sachs (1966) and later developed by,
e.g., MacCallum & Ellis (1970); Clarkson & Umeh (2011); Clarkson et al.
(2012); Umeh (2013); Heinesen (2021a, b); Maartens et al. (2023). The
application of such frameworks to the analysis of data has so far been very
limited (though see Dhawan et al., 2023; Cowell et al., 2023) due to the high
quality of data and good sky coverage necessary to measure anisotropic
features in the observables robustly. However, as data sets continue to grow,
cosmographic analyses accounting for the imprint of large-scale structure in
observations are now within reach.

Cosmographic analyses can be done with a variety of astrophysical observables.
In this paper, we focus on angular diameter distance and analyse the
anisotropic distance–redshift data for observers within $\Lambda$CDM
cosmology. We use the UNITY simulations presented in Coates et al. (2021),
which were performed with the general-relativistic $N$-body code gevolution
(Adamek et al., 2016). Using relativistic ray tracing we calculate distances
and redshifts over the full skies of synthetic observers within the
simulations. We investigate the distance–redshift cosmography from Heinesen
(2021a) as inferred by the observers, presenting constraints on the dominant
multipoles in generalised cosmological parameters for different redshift
intervals (see Heinesen & Macpherson, 2022; Cowell et al., 2023, for estimates
of the expected dominant multipoles).

Our high-resolution simulations present a realistic scatter of data points
around any possible notion of a large-scale matter congruence. Here, a matter
congruence is defined as a set of timelike integral curves that provide a
threading of the observed patch of spacetime. To make such a notion more
precise, we run a simulation where small-scale perturbations are filtered out
in the initial conditions, leading to a situation where a smooth matter
congruence is uniquely defined. Additionally, in a separate set of simulations
we compare our results from gevolution to those found in the large-scale
numerical relativity fluid simulations presented in Macpherson & Heinesen
(2021).

In Section 2, we introduce the cosmography of Heinesen (2021a) and recast the
Taylor series expansion into an arguably more suitable Padé approximant, as
well as discuss fluctuations with respect to this smooth congruence model. In
Section 3 we discuss details of the UNITY simulations and ray tracing
procedure, and in Section 4 we outline our analysis of the synthetic
observations. We present our main results in Section 5 and conclude in Section
6.

##  2 Formalism

In this section, we reformulate the anisotropic cosmography of astrophysical
sources relative to an observer in an arbirtary spacetime without exact
symmetries. Following Heinesen (2021a), we first consider a situation where
overlapping streams of the sources (shell crossing) can be ignored at the
scales of interest, such that the astrophysical sources constitute a non-
singular congruence of world lines generated by a 4–velocity field $u^{\mu}$.
In this treatment we first assume that the observer is comoving with the
congruence of sources, and then generalise the result to arbitrary observer
frames. Furthermore, we assume that the photons emitted from the sources and
received at the observer are described in the geometric optics limit, and that
photon number is conserved, such that Etherington’s reciprocity theorem
applies.

###  2.1 Taylor series cosmographic expansion

When Etherington’s reciprocity applies the luminosity distance to the emitting
source can be written as $d_{L}=d_{A}(1+z)^{2}$ for a light beam in a general
geometry, where $d_{A}$ is the angular diameter distance (also known as the
area distance) to the emitting source and $z$ is the redshift.

Consider the case where $z$ is a parameter along each of the photon null rays,
such that the function $z\mapsto d_{L}(z,\boldsymbol{e})$ is well defined for
each direction $\boldsymbol{e}$ of incoming light as seen on the observer’s
sky. If this function is analytic, its Taylor series expansion

$d_{L}(z,\boldsymbol{e})=d_{L}^{(1)}(\boldsymbol{e})\,z+d_{L}^{(2)}(\boldsymbol{e})\,z^{2}+d_{L}^{(3)}(\boldsymbol{e})\,z^{3}+\mathcal{O}(z^{4})\,$
(1)

approximates the luminosity distance in some cosmic neighbourhood of the
observer. The coefficients of the series expansion can be represented as

$\displaystyle d_{L}^{(1)}(\boldsymbol{e})$
$\displaystyle=\frac{1}{\mathfrak{H}_{o}(\boldsymbol{e})}\,,\qquad
d_{L}^{(2)}(\boldsymbol{e})=\frac{1-\mathfrak{Q}_{o}(\boldsymbol{e})}{2\mathfrak{H}_{o}(\boldsymbol{e})}\
,$ (2) $\displaystyle d_{L}^{(3)}(\boldsymbol{e})$
$\displaystyle=\frac{-1+3\mathfrak{Q}_{o}^{2}(\boldsymbol{e})+\mathfrak{Q}_{o}(\boldsymbol{e})-\mathfrak{J}_{o}(\boldsymbol{e})+\mathfrak{R}_{o}(\boldsymbol{e})}{6\mathfrak{H}_{o}(\boldsymbol{e})}\,,$

where the effective cosmological parameters

$\displaystyle\mathfrak{H}(\boldsymbol{e})$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{E^{2}}\frac{{\rm d}E}{{\rm
d}\lambda}\,,\qquad\mathfrak{Q}(\boldsymbol{e})\equiv-1-\frac{1}{E}\frac{\frac{{\rm
d}\mathfrak{H}}{{\rm d}\lambda}}{\mathfrak{H}^{2}}\,,$
$\displaystyle\mathfrak{R}(\boldsymbol{e})$ $\displaystyle\equiv$
$\displaystyle
1+\mathfrak{Q}-\frac{1}{2E^{2}}\frac{k^{\mu}k^{\nu}R_{\mu\nu}}{\mathfrak{H}^{2}}\,,$
$\displaystyle\mathfrak{J}(\boldsymbol{e})$ $\displaystyle\equiv$
$\displaystyle\frac{1}{E^{2}}\frac{\frac{{\rm d^{2}}\mathfrak{H}}{{\rm
d}\lambda^{2}}}{\mathfrak{H}^{3}}-4\mathfrak{Q}-3,$ (3)

generalise the FLRW Hubble, deceleration, curvature, and jerk parameters,
respectively, in the distance–redshift cosmography. Here, $k^{\mu}$ is the
photon 4–momentum, $\lambda$ is an affine parameter of the geodesic defined by
$k^{\mu}\nabla_{\mu}\lambda=1$, and the operator $\frac{{\rm d}}{{\rm
d}\lambda}\equiv k^{\mu}\nabla_{\mu}$ is the directional derivative along the
incoming null ray, $E=-u^{\mu}k_{\mu}$ is the photon energy as measured by the
sources/observer, and $R_{\mu\nu}$ is the Ricci curvature of the spacetime.
The effective cosmological Hubble, deceleration, jerk, and curvature
parameters $\\{\mathfrak{H},\mathfrak{Q},\mathfrak{J},\mathfrak{R}\\}$ contain
information about local kinematics and curvature effects that impact the
luminosity distance–redshift relation and might for instance include
anisotropic expansion of space, bulk flow of sources and lensing of photons.

In this framework, the anisotropic Hubble parameter can be written as the
exact multipole decomposition

$\mathfrak{H}(\boldsymbol{e})=\frac{1}{3}\theta-e^{\mu}a_{\mu}+e^{\mu}e^{\nu}\sigma_{\mu\nu}\,,$
(4)

where $\theta$ is the volume expansion rate in the source frame,
$\sigma_{\mu\nu}$ is the anisotropic deformation rate of space as traced by
matter (shear), $a^{\mu}\equiv u^{\nu}\nabla_{\nu}u^{\mu}$ is the
4–acceleration of the source frame, and $e^{\mu}$ denotes the spatial unit
four-vector in direction $\boldsymbol{e}$ in that frame. The effective
deceleration parameter can similarly be written in terms of the exact
multipole decomposition

$\mathfrak{Q}(\boldsymbol{e})=-1-\frac{1}{\mathfrak{H}^{2}(\boldsymbol{e})}\biggl{(}\overset{0}{\mathfrak{q}}+e^{\mu}\overset{1}{\mathfrak{q}}_{\mu}+e^{\mu}e^{\nu}\overset{2}{\mathfrak{q}}_{\mu\nu}\biggl{.}\\\
\biggr{.}+e^{\mu}e^{\nu}e^{\kappa}\overset{3}{\mathfrak{q}}_{\mu\nu\kappa}+e^{\mu}e^{\nu}e^{\kappa}e^{\lambda}\overset{4}{\mathfrak{q}}_{\mu\nu\kappa\lambda}\biggr{)}\,,$
(5)

with coefficients

$\displaystyle\overset{0}{\mathfrak{q}}\equiv\frac{1}{3}\frac{{\rm
d}\theta}{{\rm
d}\tau}+\frac{1}{3}D_{\mu}a^{\mu}-\frac{2}{3}a^{\mu}a_{\mu}-\frac{2}{5}\sigma_{\mu\nu}\sigma^{\mu\nu}\,,$
$\displaystyle\overset{1}{\mathfrak{q}}_{\mu}\equiv-\frac{1}{3}D_{\mu}\theta-\frac{2}{5}D_{\nu}\sigma^{\nu}_{\;\mu}-\frac{{\rm
d}a_{\mu}}{{\rm
d}\tau}+a^{\nu}\omega_{\mu\nu}+\frac{9}{5}a^{\nu}\sigma_{\mu\nu}\,,\quad$
$\displaystyle\overset{2}{\mathfrak{q}}_{\mu\nu}\equiv\frac{{\rm
d}\sigma_{\mu\nu}}{{\rm
d}\tau}+D_{\langle\mu}a_{\nu\rangle}\\!+a_{\langle\mu}a_{\nu\rangle}\\!-2\sigma_{\kappa(\mu}\omega^{\kappa}_{\nu)}\\!-\frac{6}{7}\sigma_{\kappa\langle\mu}\sigma^{\kappa}_{\nu\rangle}\,,$
$\displaystyle\overset{3}{\mathfrak{q}}_{\mu\nu\kappa}\equiv-
D_{\langle\mu}\sigma_{\nu\kappa\rangle}-3a_{\langle\mu}\sigma_{\nu\kappa\rangle}\,,$
$\displaystyle\overset{4}{\mathfrak{q}}_{\mu\nu\kappa\lambda}\equiv
2\sigma_{\langle\mu\nu}\sigma_{\kappa\lambda\rangle}\,,$ (6)

where indices between angled brackets $\langle\rangle$ denote the traceless
and symmetric part of the tensor in the involved indices, and $D_{\mu}$ is the
spatial covariant derivative in the frame of the congruence. The operator
$\frac{{\rm d}}{{\rm d}\tau}\equiv u^{\mu}\nabla_{\mu}$ is the derivative
along the world lines of the source congruence, and $\omega_{\mu\nu}$ is the
vorticity tensor describing the rotation of the congruence. The higher-order
effective cosmological parameters, $\mathfrak{J}$ and $\mathfrak{R}$, can
similarly be written in terms of multipole expansions (see Heinesen, 2021a,
for their detailed expressions).

The radius of convergence of Eq. (1) and the magnitude of the remainder for
the truncation of the cosmography at third order depend on the particularities
of the spacetime and the astrophysical survey geometry, and must in principle
be assessed for each physical scenario of interest (see, e.g., Macpherson &
Heinesen, 2021, for estimates of the remainder terms within some numerical
examples).

###  2.2 Padé approximant cosmography

We may take Eq. (1) as a first ansatz for an approximation of the cosmic
neighbourhood of the observer. However, in the real Universe, the congruence
described by our time-like vector field $u^{\mu}$ can be understood as
capturing the coarse-grained dynamics at large scales. Individual
astrophysical sources would then have peculiar motion with respect to this
large-scale average, which would manifest as a scatter in the Hubble diagram.
The observed redshift therefore has particularly large random fluctuations
which renders it discontinuous and non-monotonic along the lines of sight in
the late-time, nonlinear Universe. On the other hand, a continuous and
monotonic observable — at least at low cosmological distances — is the angular
diameter distance. When considered as a function of the affine parameter along
the line of sight, $d_{A}$ is mainly perturbed by weak gravitational lensing
and by the relativistic aberration due to the peculiar motion of the observer.
These perturbations are very smooth as long as no caustics are encountered on
the light cone. In attempt to arrive at a framework which is useful for
observational constraints in the late-time Universe, we therefore reformulate
the distance–redshift relation using $d_{A}$ as the parameter in which the
series expansion is performed. The inverse series of Eq. (1) reads

$\displaystyle z(\boldsymbol{e},d_{A})$ $\displaystyle=$
$\displaystyle\frac{1}{d_{L}^{(1)}\\!(\boldsymbol{e})}d_{A}+\frac{2d_{L}^{(1)}\\!(\boldsymbol{e})-d_{L}^{(2)}\\!(\boldsymbol{e})}{\left[d_{L}^{(1)}\\!(\boldsymbol{e})\right]^{3}}d_{A}^{2}+\frac{5\left[d_{L}^{(1)}\\!(\boldsymbol{e})\right]^{2}\\!\\!+2\left[d_{L}^{(2)}\\!(\boldsymbol{e})\right]^{2}\\!\\!-6d_{L}^{(1)}\\!(\boldsymbol{e})\,d_{L}^{(2)}\\!(\boldsymbol{e})-d_{L}^{(1)}\\!(\boldsymbol{e})\,d_{L}^{(3)}\\!(\boldsymbol{e})}{\left[d_{L}^{(1)}\\!(\boldsymbol{e})\right]^{5}}d_{A}^{3}+\mathcal{O}(d_{A}^{4})\,,$
(7) $\displaystyle=$
$\displaystyle\mathfrak{H}_{o}(\boldsymbol{e})d_{A}+\frac{1}{2}\left[\mathfrak{Q}_{o}(\boldsymbol{e})+3\right]\left[\mathfrak{H}_{o}(\boldsymbol{e})d_{A}\right]^{2}+\frac{1}{6}\left[16+11\mathfrak{Q}_{o}(\boldsymbol{e})+\mathfrak{J}_{o}(\boldsymbol{e})-\mathfrak{R}_{o}(\boldsymbol{e})\right]\left[\mathfrak{H}_{o}(\boldsymbol{e})d_{A}\right]^{3}+\mathcal{O}(d_{A}^{4})\,.$

As suggested by Heinesen (2021a), it is convenient to define the following:

$\displaystyle\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})$ $\displaystyle\equiv$
$\displaystyle\mathfrak{H}^{2}_{o}(\boldsymbol{e})\left[\mathfrak{Q}_{o}(\boldsymbol{e})+1\right]\,,$
(8a) $\displaystyle\hat{\mathfrak{J}}_{o}(\boldsymbol{e})$
$\displaystyle\equiv$
$\displaystyle\mathfrak{H}^{3}_{o}(\boldsymbol{e})\left[\mathfrak{J}_{o}(\boldsymbol{e})-\mathfrak{R}_{o}(\boldsymbol{e})-1\right]\,,$
(8b)

which have a naturally-truncating, finite multipole decomposition within the
cosmographic framework. Each multipole of these parameters contains physical
information about the kinematical properties of the congruence $u^{\mu}$. For
the case of $\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})$, this can be seen in Eq.
(2.1) (see also the discussion in Umeh, 2013), while for
$\hat{\mathfrak{J}}_{o}(\boldsymbol{e})$ this is given in Appendix B of
Heinesen (2021a). With the new quantities in Eqs. (8), Eq. (7) simplifies to

$z(\boldsymbol{e},d_{A})=\mathfrak{H}_{o}(\boldsymbol{e})d_{A}+\left[\mathfrak{H}_{o}^{2}(\boldsymbol{e})+\frac{1}{2}\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})\right]d_{A}^{2}\\\
+\left[\mathfrak{H}_{o}^{3}(\boldsymbol{e})+\frac{\hat{\mathfrak{J}}_{o}(\boldsymbol{e})+11\mathfrak{H}_{o}(\boldsymbol{e})\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})}{6}\right]d_{A}^{3}+\mathcal{O}(d_{A}^{4})\,.$
(9)

The third-order Taylor series expansions discussed in this section are
expected to provide good approximations to the cosmographic model at very
small distances. However, at larger distances the truncation error can
increase very rapidly. An alternative series expansion is the Padé
approximant, which can often provide improved behaviour over a finite distance
interval. Padé approximants have already been proposed and implemented in the
place of Taylor series expansions for more accurate constraints from
distance–redshift data (e.g. Nesseris & García-Bellido, 2013; Sapone et al.,
2014; Demianski et al., 2017), and have been shown to have radii of
convergence beyond $z=1$ (Gruber & Luongo, 2014), which is the strict upper
limit for the radius of convergence in Taylor series expansions.

We propose the approximant

$z(\boldsymbol{e},d_{A})=\frac{\mathfrak{H}_{o}(\boldsymbol{e})d_{A}+\frac{1}{2}\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})d_{A}^{2}+\frac{1}{6}\hat{\mathfrak{J}}_{o}(\boldsymbol{e})d_{A}^{3}+\ldots}{1-\mathfrak{H}_{o}(\boldsymbol{e})d_{A}-\frac{4}{3}\hat{\mathfrak{Q}}_{o}(\boldsymbol{e})d_{A}^{2}-\ldots}\,,$
(10)

which is order $[3/2]$ and has the same Taylor series expansion as Eq. (9) up
to $\mathcal{O}(d_{A}^{4})$. Since we only fix the first three terms in the
Taylor series, other Padé approximants of order $[3/2]$ are possible, and our
choice is mainly guided by the simplicity of the expression in terms of our
preferred variables. In the FLRW limit, the truncated Taylor series (9) has an
error of $\approx 0.8\%$ with respect to the exact $\Lambda$CDM angular
diameter distance at redshift $z\approx 0.15$. Our expression (10) using a
Padé approximant with the same number of parameters reduces this error to
better than $0.2\%$ (see Figure 1).

Figure 1: We show the relative error on the distance–redshift relation as a
function of redshift for the third-order Taylor series expansion (9) and the
$[3/2]$-order Padé approximant (10) in red and blue, respectively. Both
approximations use the same expansion rate, deceleration parameter and jerk,
and were compared to an exact $\Lambda$CDM background solution with
$\Omega_{\mathrm{m}}=0.31$, $\Omega_{\mathrm{rad}}=10^{-4}$, and vanishing
curvature.

###  2.3 Fluctuations about the congruence

Eq. (10) can be understood as a model for the ‘cosmographic’ redshift of the
underlying coarse-grained matter congruence. The observed redshift of
astrophysical sources is more realistically given by the model plus some
stochastic component $\delta z$, representing the peculiar velocities of those
sources with respect to the inhomogeneous congruence, i.e.

$z_{\mathrm{obs}}=z(\boldsymbol{e},d_{A})+\delta z\,.$ (11)

This expression combined with Eq. (10) generalises the common interpretation
of observed redshifts in an FLRW cosmology to the inhomogeneous case.

Here we are still assuming that the observation is made in the rest frame of
the congruence. In reality, the observer can have some velocity with respect
to the large-scale congruence which may even need to be determined from the
same data. It is therefore more useful to describe the observation in an
arbitrary frame and allow for three additional parameters that characterise
the boost between the frames. These three parameters are included in our joint
parameter fit, and thus no knowledge about the velocity of the observer is
assumed.

We denote the peculiar three-velocity vector of the observer with respect to
the congruence $u^{\mu}$ as $\boldsymbol{\beta}$, with length $\beta=v/c$.
Special relativity implies that the observed redshift in the boosted frame,
$z^{\prime}_{\mathrm{obs}}$, is given by

$1+z^{\prime}_{\mathrm{obs}}=\frac{1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}}{\sqrt{1-\beta^{2}}}\left(1+z_{\mathrm{obs}}\right)\,,$
(12)

where $\boldsymbol{e}^{\prime}$ is the observed position in the boosted frame.
The latter is different from $\boldsymbol{e}$ in the rest frame of the
congruence due to relativistic aberration. They are related by

$\boldsymbol{e}=\frac{\sqrt{1-\beta^{2}}}{1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}}\boldsymbol{e}^{\prime}-\frac{1-\left(1-\sqrt{1-\beta^{2}}\right)\beta^{-2}\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}}{1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}}\boldsymbol{\beta}\,.$
(13)

The aberration also changes the angular distance which follows from
$d_{A}^{2}d\Omega=\left(d^{\prime}_{A}\right)^{2}d\Omega^{\prime}$, where
$d\Omega$ is the solid angle element. We therefore find the well-known result
that

$d_{A}=\frac{1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}}{\sqrt{1-\beta^{2}}}d^{\prime}_{A}\,.$
(14)

Putting everything together we have a model for the observed redshift in the
boosted frame. If the velocity is small, i.e. $v\ll c$, we find

$z^{\prime}_{\mathrm{obs}}\approx\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)z\left[\left(1+\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)\boldsymbol{e}^{\prime}-\boldsymbol{\beta},\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)d^{\prime}_{A}\right]\\\
-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}+\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)\delta
z\,.$ (15)

The aberration of the original anisotropy pattern is a second-order effect
since we may write

$z(\boldsymbol{e},d_{A})=z(\boldsymbol{e}^{\prime},d_{A})+(\boldsymbol{e}-\boldsymbol{e}^{\prime})\cdot\boldsymbol{\nabla}_{\boldsymbol{e}^{\prime}}z(\boldsymbol{e}^{\prime},d_{A})+\ldots\,,$
(16)

with
$\boldsymbol{e}-\boldsymbol{e}^{\prime}\approx(\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta})\boldsymbol{e}^{\prime}-\boldsymbol{\beta}$.
The second- and higher-order corrections in above equation are therefore
strongly suppressed for the low multipoles of the distance–redshift relation
we are interested in, allowing us to make the further approximation that
$z(\boldsymbol{e},d_{A})\approx z(\boldsymbol{e}^{\prime},d_{A})$. If we also
drop the term $\beta\times\delta z$ from Eq. (15), which is quadratic in
velocities, we finally arrive at

$z^{\prime}_{\mathrm{obs}}\approx\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)z\left[\boldsymbol{e}^{\prime},\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)d^{\prime}_{A}\right]-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}+\delta
z\,.$ (17)

This is the expression we use to constrain the anisotropies in our cosmography
model (10) using the observed redshifts, angular diameter distances, and sky
positions measured by the observers in the simulation. We describe the process
of constraining the parameters of the cosmography in Section 4.2.

##  3 Simulations

In order to assess how well the parameters encoding the anisotropy of the
distance–redshift relation can be measured when assuming the cosmographic
model presented in the previous section, we use synthetic source catalogues
that we generate from simulations. In Section 3.1, we describe the existing
UNITY simulation data we use for our main analysis. In Section 3.2 we
introduce new, low-resolution simulations, which represent ‘smooth’ analogues
of the UNITY simulations from our main analysis.

###  3.1 UNITY simulation suite

We use data from existing simulations described in Coates et al. (2021): the
UNITY simulation suite consists of 34 independent realisations of a flat
$\Lambda$CDM Universe with cosmological parameters $h=0.67$,
$\Omega_{\mathrm{b}}=0.049$, $\Omega_{\mathrm{c}}=0.27$,
$T_{\mathrm{CMB}}=2.7255\,\mathrm{K}$ and a minimal-mass normal hierarchy of
neutrinos with $\sum m_{\nu}=0.06\,\mathrm{eV}$. Gaussian initial conditions
follow a power law for primordial perturbations with amplitude
$A_{s}=2.1\times 10^{-9}$ at the pivot scale of $0.05\,\mathrm{Mpc}^{-1}$ and
spectral index $n_{s}=0.96$. The simulated volume is a cube of
$4032\,\mathrm{Mpc}/h$ on each side, and the matter (baryons and cold dark
matter are treated as a single species) is sampled using $2304^{3}$ mass
elements that form a collisionless $N$-body ensemble. In the present analysis
we only use one of the realisations available, but the remaining ones could be
analysed in the future to generate an ensemble of observers.

The simulations were performed using the $N$-body code gevolution (Adamek et
al., 2016) which adopts a weak-field expansion of general relativity to
incorporate relativistic effects. The numerical output consists of data for
the matter distribution and metric perturbations on the past light cone of an
observer placed inside the simulation. For convenience, the data are
discretised in terms of concentric shells around the observer that are evenly
spaced in comoving coordinates at intervals of $1.75\,\mathrm{Mpc}/h$. Each
shell is pixelised using the HEALPix framework (Górski et al., 2005) at a
resolution adapted to the surface area of the shell in order to sample the
simulated matter and metric perturbations as close as possible to the
resolution of $1.75\,\mathrm{Mpc}/h$, which corresponds to the Cartesian mesh
used in the particle-mesh algorithm of the code.

The discretised data are post processed using a simple ray tracer that works
in the Born approximation (see Lepori et al., 2020, for details of the ray-
tracing procedure). Assuming that the lines of sight are unperturbed has the
benefit that observables can be computed very efficiently by essentially
performing weighted sums of pixel maps corresponding to consecutive shells
down the past light cone. The ray tracer computes two observables for each
pixel: the observed redshift $z_{\mathrm{obs}}^{\prime}$ and the angular
distance $d_{A}^{\prime}$, where here the primes indicate that the observation
is done in the rest frame corresponding to the time slicing of Poisson gauge,
which closely matches the rest frame of the cosmic microwave background.

Each observable includes all relativistic perturbations at first order. For
the observed redshift these are the Doppler shift due to peculiar motion of
the source, the gravitational redshift, the integrated Sachs-Wolfe effect and
the Shapiro time delay, while the angular distance is only perturbed by the
weak-lensing convergence. In order to assign the peculiar velocity we take a
mass-weighted average of the velocity of all particles in each pixel. It is
worth pointing out that, although all effects are added linearly by the ray
tracer, they arise from a completely non-perturbative underlying solution.

For this work we assume that astrophysical sources are good tracers of the
underlying matter, i.e. we sidestep the issue of bias and use a density of
sources that is directly proportional to the matter density. This also avoids
the issue of empty pixels in which no unique velocity can be assigned.

###  3.2 ‘Smooth’ simulations

The simulations described in the previous section contain small-scale
nonlinear structures at low redshift, thus, constructing a large-scale
congruence description of these simulations is non-trivial (e.g. Clarkson et
al., 2011). We attempt to do this by directly fitting the cosmographic
framework outlined in Section 2, which essentially introduces a ‘smoothing
length’ corresponding to the redshift range of the data we include. We also
perform the same fits using a simulation which should have a uniquely defined
large-scale congruence description. By comparing this fit to direct
calculations of the effective cosmological parameters (2.1) in the same
simulation (in a manner equivalent to those presented in Macpherson &
Heinesen, 2021), we can assess how well the fit is matching the ‘true’
congruence description of the simulation.

We take the UNITY simulation outlined in the previous section and remove all
perturbations with wavelengths shorter than $670\,\mathrm{Mpc}$ from the
initial data by applying a low-pass filter. We decrease the spatial resolution
of this simulation by a factor of four and the mass resolution by a factor of
eight — justified by the lack of small-scale power in the simulation. Aside
from this resolution change, we perform the simulation with parameters as
described above. The filter is applied only to the initial data, and naturally
some structures beneath this scale will form through the course of the
simulation. However, the power on small scales will be significantly reduced
at late times with respect to the original UNITY simulations. We thus consider
the scale of $670\,\mathrm{Mpc}$ as an approximate smoothing scale for the
late-time simulation, which is much larger than the resolution scale of
approximately $10\,\mathrm{Mpc}$. For this simulation, we use the analysis
code of Macpherson & Heinesen (2021) to calculate the effective cosmological
parameters (2.1) by taking numerical derivatives of dynamical quantities at
the observer position. We then use the data obtained by ray tracing to
constrain these same parameters and compare the results.

Macpherson & Heinesen (2021) performed their analysis on simulations from the
Einstein Toolkit (Löffler et al., 2012). To check for consistency, we
replicated one of those simulations with gevolution and repeated the analysis.
The Einstein Toolkit uses a fluid approximation for the matter field in
numerical relativity while gevolution uses an $N$-body particle description in
weak-field general relativity. We find good agreement between these two cases
(see Appendix A) and thus conclude that, under the right conditions, we can
find a large-scale congruence description of the $N$-body particle simulation
data.

It is worth pointing out that due to the non-commutation of averaging and time
evolution in nonlinear general relativity, we do not in general expect the
large-scale congruence of the UNITY simulation to coincide with that of the
simulation with smoothed initial data.

##  4 Parameter constraints

In this section we discuss details of our method for constraining the
anisotropic cosmography using synthetic observations from the UNITY
simulations. First, we describe the assumptions we make and how they reduce
the degrees of freedom of the model cosmography, and then we outline the
pipeline we use to constrain these degrees of freedom using simulation data.

###  4.1 Dominant multipoles

The effective cosmological parameters in the Padé approximant (10) of the
distance–redshift relation contain a hierarchy of multipoles. The multipole
decompositions of the Hubble and deceleration parameters are given in (4) and
(5), respectively, and the curvature and jerk parameter decompositions can be
found in Heinesen (2021a). This full hierarchy contains many degrees of
freedom (though see Heinesen & Macpherson, 2022, for a reduction in the
degrees of freedom by a factor of two), and constraining all of them with
current data is difficult. We can make some simplifying assumptions to reduce
the number of degrees of freedom and allow the best possible constraints. For
a congruence of collisionless matter the acceleration vector $a^{\mu}$
vanishes, which removes the dipole component of the anisotropic Hubble
parameter. In this work, we thus focus on constraining the monopole and
quadrupole of the Hubble parameter as well as the independent multipoles of
the deceleration parameter (monopole, dipole, quadrupole and octopole).
According to Eq. (2.1), the hexadecapole of the deceleration parameter is
fully determined by the quadrupole of the Hubble parameter, and we build this
constraint into our model. We fix $\hat{\mathfrak{J}}_{o}=0$ which is close to
the value for a flat FLRW universe with a realistic value of the present
energy density of radiation,
$\hat{\mathfrak{J}}_{o}/H_{0}^{3}=2\,\Omega_{\mathrm{rad}}\approx 10^{-4}$.
Note the “$-1$” that enters the definition of $\hat{\mathfrak{J}}_{o}$ in Eq.
(8) to see the connection to the canonical jerk parameter.

###  4.2 Markov chain Monte Carlo (MCMC) analysis

We constrain the components of each anisotropy tensor by modelling the
logarithm of the likelihood as

$\ln\mathcal{L}=-\frac{1}{2}\sum_{i}\left[\frac{\delta
z^{(i)}}{\sigma\left(d_{A}^{(i)}\right)}\right]^{2}-\sum_{i}\ln\sigma\left(d_{A}^{(i)}\right)\,,$
(18)

up to an irrelevant constant. In this formula, we compute the residual $\delta
z^{(i)}$ for each individual source from Eq. (17), i.e., we solve that
equation for $\delta z$ using the observed redshift
$z^{\prime}_{\mathrm{obs}}$, the observed direction $\boldsymbol{e}^{\prime}$
and the observed angular distance $d_{A}^{\prime}$ of each observed pixel $i$.
The model parameters are given by the boost velocity vector
$\boldsymbol{\beta}$ (three components) and all the parameters that enter Eq.
(10) describing the cosmography (22 independent parameters after fixing the
four-acceleration and the jerk). Furthermore, we add two nuisance parameters
to model the Gaussian scatter $\sigma(d_{A})$,

$\sigma\left(d_{A}\right)=\frac{\sigma_{0}}{1+s_{1}d_{A}}\approx\frac{\sigma_{0}}{1+s_{1}\left(1-\boldsymbol{e}^{\prime}\cdot\boldsymbol{\beta}\right)d^{\prime}_{A}}\,,$
(19)

in order to allow a certain evolution of the scatter with distance. This means
that the Gaussian scatter needs to be evaluated for each pixel separately, and
we try to find a global fit for $\sigma_{0}$ and $s_{1}$. The underlying
assumption of this likelihood model is that $\delta z$ is indeed a Gaussian
random variable with zero expectation value and a variance that mildly depends
on angular distance. Our results show that this is a reasonable assumption in
the presence of small-scale nonlinear structure, but not such a good
approximation for the case where we filtered the initial conditions. Details
are discussed in the next section.

Our data vectors consist of randomly selected pixels sampling the full sky at
redshifts below a chosen threshold. We set the probability for selecting a
pixel proportional to its comoving matter density. Since we want to model the
distance–redshift relation over a certain baseline distance, we want the
sources roughly evenly distributed across that distance. The volume factor
naturally favours objects at large distance, and we compensate this fact by
scaling the probability by an additional factor proportional to the solid
angle of the pixel. In this way, for any given solid angle of the
observations, the data are roughly evenly distributed in distance. This is of
course just one possible selection function, and we expect that different
selection functions may lead to different results depending on how much they
favour low-redshift sources over high-redshift ones, or vice versa. We use
different redshift thresholds to gain some first insights into this effect.

Our analysis does not require a lower redshift cutoff and we even include
sources with negative redshift. In fact, it may be important to keep such
sources for an accurate fit of the observer velocity. We vary the redshift
range of sources in our synthetic catalogs up to a maximum of $z=0.15$ to
ensure a small truncation error in our Padé approximant. Keeping the number of
sources at low redshift approximately constant, this results in a different
number of sources as we increase the maximum redshift. In the UNITY simulation
this number ranges from about $2000$ sources at $z=0.04$ up to about $7500$
sources at $z=0.15$. For the ‘smooth’ simulation we have fewer sources,
ranging between about $500$ at $z=0.025$ to about $2000$ at $z=0.15$.

##  5 Results

Figure 2: Multipole power
$C_{\ell}=\left(2\ell+1\right)^{-1}\sum_{m}\left|a_{\ell m}\right|^{2}$ for
the quadrupole ($\ell=2$) of $\mathfrak{H}_{o}/H_{0}$ which measures the shear
tensor of the coarse-grained matter congruence relative to the expansion, and
the independent multipoles of $\hat{\mathfrak{Q}}_{o}/H_{0}^{2}$ that
characterise the anisotropy in the deceleration parameter. The error bars
indicate $95\%$ confidence intervals. For all anisotropy parameters the power
drops nearly exponentially as the fitted range in redshift $z$ increases
(horizontal axis). The dashed lines show results from a simulation where
perturbations below the scale of 670 Mpc were filtered out in the initial
data, and hence a well-defined matter congruence exists. For this case, the
horizontal dotted lines show an alternative measurement using the properties
of null geodesics at the observer location ($z\rightarrow 0$).

Figure 2 shows the multipole power
$C_{\ell}=\left(2\ell+1\right)^{-1}\sum_{m}\left|a_{\ell m}\right|^{2}$ as a
function of the maximum redshift $z$ of the fitting interval. Solid curves
represent constraints in the fully nonlinear UNITY simulation and dashed
curves represent constraints in the simulation with filtered initial data. We
show the power in the quadrupole in $\mathfrak{H}_{o}/H_{0}$ (dark blue) as
well as the dipole (brown), quadrupole (red), and octopole (orange) in
$\hat{\mathfrak{Q}}_{o}/H_{0}^{2}$, where $H_{0}=\theta/3$ is the monopole of
the expansion rate $\mathfrak{H}_{o}$. The error bars indicate $95\%$
confidence intervals. The smaller error bars in the ‘smooth’ case are due to
the much reduced scatter in the data, owing to the lack of small-scale
perturbations. The horizontal dotted lines (without error bars) indicate the
multipole amplitudes found by taking numerical derivatives along null rays at
the observer location in the ‘smooth’ simulation, following the analysis of
Macpherson & Heinesen (2021) (see Appendix A for more details).

The ‘smooth’ simulations have structure below 670 Mpc removed from the initial
data, corresponding to the diameter of the light cone at $z\approx 0.08$.
Consequently there is relatively little evolution of the multipole amplitude
across the entire redshift range probed. The opposite is true for the
simulation that includes small-scale nonlinear structure: anisotropies are
very large if the fitting range includes only very low redshifts, and they
drop nearly exponentially as the fitting range is enlarged. However, the
anisotropy levels always remain larger than in the case of the smoothed
initial data, even beyond the smoothing scale. We reiterate that the initial
data of these simulations are otherwise identical; removing structure beneath
670 Mpc is the only difference.

The dotted lines in Figure 2, which are obtained by taking numerical
derivatives of the distance–redshift relation at the observer location in the
‘smooth’ simulation, are mostly in fair agreement with our MCMC analysis, with
the notable exception of the quadrupole of $\hat{\mathfrak{Q}}_{o}$ where the
measurements appear to disagree by more than two orders of magnitude. In this
particular case, the apparent disagreement is a result of projecting the
posterior onto $C_{\ell}$, which is a quadratic quantity. We show in Figure 7
(discussed in more detail below) that the posterior of the quadrupole
components of $\hat{\mathfrak{Q}}_{o}$ at low redshift encloses the zero
quadrupole value within its $1\,\sigma$-contours. However, $C_{\ell}$
effectively measures the square distance from zero in the five-dimensional
parameter space of quadrupole components. The steep increase in volume as one
moves away from zero initially overcomes the decrease in likelihood, giving a
large value for the marginalised $C_{\ell}$ but also an error bar of similar
order. This effect is problematic for multipoles whose full posterior overlaps
with the origin. On the other hand, if the entire posterior can be seen to be
far away from zero, the $C_{\ell}$ become a useful measure of the distance.

The physical reason for the expected smallness of the quadrupole term relative
to the dipole and octopole of $\hat{\mathfrak{Q}}_{o}$ can be seen in Eq.
(2.1): while the dipole and octopole are dominated by spatial gradients of the
matter flow field, the leading term of the quadrupole is the time derivative
of the shear (assuming a geodesic flow). For cosmological matter
perturbations, time derivatives are generally much smaller than spatial
gradients.

Cowell et al. (2023) placed an upper limit of $2.88\%$ on the amplitude of the
quadrupole in the Hubble parameter — relative to the monopole — using
Pantheon+ supernovae data. The data in this work lay in the redshift range
$z=0.023$ to $z=0.15$. The amplitude for a maximum redshift of $z=0.15$ that
we find from $C_{2}$ in Figure 2 is $\lesssim 1.6\%$ (at $95\%$ confidence),
which is consistent with the current upper limit from supernovae and suggests
that a detection should be possible in the near future. We note this amplitude
is similar for both the UNITY and ‘smooth’ simulations we use here.

Figure 3: 2D marginalised posteriors of $\sigma_{0}$, $H_{0}$ (in units of
$100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$) and $q$ for different
fitting ranges in redshift $z$ (different colours). The sampling is roughly
uniform in observed distance $d_{A}$ to the sources. Results from a simulation
where perturbations below the scale of 670 Mpc were filtered out in the
initial data are shown in orange and red. For this case, the red dashed lines
show an alternative measurement of $H_{0}$ and $q$ using numerical derivatives
along null geodesics at the observer location ($z\rightarrow 0$). The values
for the fiducial FLRW background model are indicated by grey dashed lines. The
inset labelled “Zoom-in” magnifies a region highlighted in the 2D plot for the
parameters $H_{0}$ and $q$.

Figure 3 shows two-dimensional marginalised posterior distributions of the
monopole of the Hubble parameter, $H_{0}$, and the monopole of the
deceleration parameter, $q\equiv-1-\overset{0}{\mathfrak{q}}/H_{0}^{2}$. The
exact monopole of $\mathfrak{Q}_{o}$ would receive small corrections due to
the anisotropy of $\mathfrak{H}_{o}$, but we neglect these higher-order terms
here. In other words, $q$ parameterises the exact monopole of
$\hat{\mathfrak{Q}}_{o}$ which has a direct interpretation in terms of
kinematical quantities of the cosmographic model as specified in Eq. (2.1).
Figure 3 also includes posteriors for the nuisance parameter $\sigma_{0}$
defined in (19). Grey to progressively darker blue contours show an increasing
maximum redshift of the sources as stated in the legend. Orange and red
contours show constraints in the ‘smooth’ simulations. Grey dashed lines in
the panels with $H_{0}$ and $q$ represent the true values of these parameters
in the FLRW limit of the simulations, while the red dashed lines indicate the
corresponding measurements based on numerical derivatives at the observer
location in the ‘smooth’ simulation.

As a general trend, the posterior contours tend to move through parameter
space and progressively shrink as the redshift range used in the fit
increases. While the tightening of contours is partly due to the increase of
the number of sources, the path of the best fit illustrates the dependence on
the sampled volume. At low redshift, $q$ is particularly poorly constrained,
as may be expected, but even at $z=0.15$ the contours have not yet converged
fully to the FLRW limit. For the ‘smooth’ simulation, the contours at low and
high redshift are roughly compatible, and both of them are in fair agreement
with the UNITY simulation when data up to $z=0.1$ or higher are included in
the fit, see the panel labelled “Zoom-in” in Figure 3. The fitted value for
$\sigma_{0}$ for UNITY is clustered around $8\times 10^{-4}$ which is
commensurate with the expected velocity dispersion. In the case with filtered
initial data the scatter is significantly smaller, of the order of few $\times
10^{-5}$.

The monopole values found via numerical derivatives (dotted lines in Figure 3)
should correspond closely to the constrained parameters in the ‘smooth’
simulation case. While we find good agreement in the case of $H_{0}$, the
deceleration parameter shows disagreement between these two calculations at
several $\sigma$.

Figure 4: Scatter of individual sources around the modelled distance–redshift
relation. The best-fit cosmography is shown in blue, and the scatter is
normalised according to the best-fit Gaussian inferred in the likelihood
analysis (dark blue line). The fiducial FLRW model is shown in red, and we
normalise the scatter by the standard deviation in this case, i.e. without
assuming a model. The histogram shows that the assumed Gaussian likelihood
model provides a fair description of the data, here for the case of
approximately 5000 sources with redshifts $z<0.1$.

Figure 4 shows a histogram of the observed scatter in redshift of the
individual sources in the UNITY simulation with respect to the best-fit
cosmographic model, Eq. (10), and the fiducial FLRW model (blue and red,
respectively). In the former case, the scatter of each source is normalised to
the inferred standard deviation $\sigma\left(d_{A}^{(i)}\right)$ of the best-
fit model, Eq. (19), whereas in the latter case we normalise the scatter to
the standard deviation of the entire distribution. We find an excellent
agreement with the Gaussian model assumed in the likelihood analysis (dark
blue curve), and the normality test using the Kolmogorov-Smirnov statistic is
passed easily with a p-value of $0.2$ (rising to $0.9$ if the sample mean is
subtracted). The scatter with respect to the fiducial FLRW model is also
Gaussian, yet with a standard deviation that is about $19\%$ larger than for
the best-fit cosmography. In addition, the central value of the distribution
is biased significantly — the Kolmogorov-Smirnov test rules out a zero-
centered Gaussian at more than $10\,\sigma$ in this case.

A corresponding histogram for the ‘smooth’ simulation is shown in Figure 5.
Since there is no intrinsic velocity dispersion in this model the small
residual redshift errors are dominated by non-Gaussian numerical noise due to
discretisation. We confirm this by conducting a resolution study, showing that
the scatter almost doubles when we decrease the resolution of the ‘smooth’
simulation by a factor of two. As a consequence, the Gaussian likelihood model
is not a good approximation in this case, and we caution that the
corresponding posteriors of our MCMC analysis should therefore be understood
as being only indicative.

Figure 5: Same as Fig. 4, but for the simulation with filtered initial
conditions and using approximately 2000 sources with redshifts $z<0.1$. The
Gaussian likelihood model is clearly not a very good description of the data
in this case.

It is interesting to note that the cosmographic model provides a better and
less biased description of the data than the fiducial FLRW model in all cases
shown in Figures 4 and 5.

Figure 6: 2D marginalised posteriors of the quadrupolar contributions to
$\mathfrak{H}_{o}$ for different fitting ranges in redshift $z$ (different
colours). The sampling is roughly uniform in observed distance $d_{A}$ to the
sources. The five independent components of $\sigma_{\mu\nu}$, which is the
shear tensor of the fitted coarse-grained matter congruence, govern the
quadrupole anisotropy of the Hubble parameter. Results from a simulation where
perturbations below the scale of 670 Mpc were filtered out in the initial data
are shown in orange and red. For this case, the red dashed lines show an
alternative measurement using the properties of null geodesics at the observer
location ($z\rightarrow 0$). All shown quantities would be zero in an exact
FLRW model (grey dashed lines).

Figure 6 shows two–dimensional posterior contours for the five independent
components of the quadrupole of $\mathfrak{H}_{o}$, the shear tensor. The
anisotropy in the expansion rate is typically of the order of a few percent in
the redshift ranges studied, and systematically decreases towards larger
redshift ranges. It is worth pointing out that the contours generally exclude
the isotropic FLRW model (grey dashed lines) with high significance, and that
the convergence towards isotropy does not follow a straight path in parameter
space. This means that the best-fit quadrupole shape depends quite strongly on
the redshift window over which it is inferred. The red dashed lines in Figure
6 indicate the quadrupole shape determined from numerical derivatives at the
observer location within the ‘smooth’ simulation. The corresponding contours
from the MCMC analysis (very tight orange and red contours) are consistent
with themselves over varying redshifts and agree extremely well with this
alternative measurement. We remind the reader that the confidence contours
should not be taken at face value in this case, since we know that the
Gaussian likelihood model is not a good description.

Figure 7: 2D marginalised posteriors of the quadrupolar contributions to the
distance–redshift relation for different fitting ranges in redshift $z$
(different colours). The sampling is roughly uniform in observed distance
$d_{A}$ to the sources. As in Figure 6, the five independent components of
$\sigma_{\mu\nu}$ govern the quadrupole anisotropy of the Hubble parameter,
whereas the five components of $\overset{2}{\mathfrak{q}}_{\mu\nu}$
characterise the quadrupole anisotropy of the deceleration parameter. Results
from a simulation where perturbations below the scale of 670 Mpc were filtered
out in the initial data are shown in orange and red. For this case, the red
dashed lines show an alternative measurement using the properties of null
geodesics at the observer location ($z\rightarrow 0$). All shown quantities
would be zero in an exact FLRW model (grey dashed lines). The panels
highlighted in light green show the degeneracies between pairs of parameters
that affect the same quadrupole anisotropy.

Figure 7 shows the same two-dimensional posteriors of the quadrupole
components of $\mathfrak{H}_{o}$ together with the quadrupole components of
the deceleration parameter $\hat{\mathfrak{Q}}_{o}$. The panels showing the
shear, in the lower-right part of the plot, are zoomed-in versions of Figure 6
to make the contours for the ‘smooth’ simulation more visible. The upper left
part of the plot shows the five independent components of the quadrupole of
$\hat{\mathfrak{Q}}_{o}$. As noted earlier, the posterior for the individual
quadrupole components in the ‘smooth’ simulation, when measured at low
redshift (orange contours), is compatible with zero at $1\,\sigma$ in all
cases. The large value of $C_{\ell}$ seen in Figure 2 therefore needs to be
interpreted as a projection effect from five-dimensional parameter space. The
five panels highlighted in light green in the lower left part of Figure 7 show
the degeneracy between $\mathfrak{H}_{o}$ and $\hat{\mathfrak{Q}}_{o}$ in the
anisotropic distance–redshift relation. The degeneracy is considerably weaker
in the ‘smooth’ simulation where the quadrupole of $\mathfrak{H}_{o}$ can be
measured very robustly. In the UNITY simulation the degeneracy remains
significant even at high redshift. Using prior knowledge about the expected
smallness of the deceleration quadrupole might therefore yield a more robust
constraint.

Figure 8: 2D marginalised posteriors of dipolar contributions to the
distance–redshift relation for different fitting ranges in redshift $z$
(different colours). The sampling is roughly uniform in observed distance
$d_{A}$ to the sources. The three components of
$\overset{1}{\mathfrak{q}}_{\mu}$ characterise the dipole anisotropy of the
deceleration parameter whereas the three components of the velocity vector
$\boldsymbol{\beta}$ parameterise the peculiar velocity of the observer with
respect to the fitted coarse-grained matter congruence. Results from a
simulation where perturbations below the scale of 670 Mpc were filtered out in
the initial data are shown in orange and red. For this case, the red dashed
lines show an alternative measurement using the properties of null geodesics
at the observer location ($z\rightarrow 0$). All shown quantities would be
zero in an exact FLRW model (grey dashed lines). The panels highlighted in
light green show the degeneracies between pairs of parameters that affect the
same dipole axis.

In Figure 8 we show two-dimensional constraints on the components of the
dipolar contributions to the distance–redshift relation. We assume from the
outset that the expansion rate has vanishing dipole as it would be
proportional to the acceleration at the observer, see Eq. (4). Therefore, the
possible sources of a dipole anisotropy are the dipole in the deceleration
parameter,
$(\overset{1}{\mathfrak{q}}_{x},\overset{1}{\mathfrak{q}}_{y},\overset{1}{\mathfrak{q}}_{z})$,
and the observer velocity with respect to the best-fit congruence,
$(\beta_{x},\beta_{y},\beta_{z})$. As previously, with increasing maximum
redshift we generally see the contours for the deceleration anisotropy
progress towards the isotropic limit. This progression is much slower for the
observer boost, which is still very far away from zero even for the largest
redshift interval that goes to $z=0.15$. For that case, the posteriors of the
components of $\boldsymbol{\beta}$ overlap with the ones measured in the
‘smooth’ simulation. For each of the three dipole axes, there is a noticeable
degeneracy between the amplitude of the deceleration dipole and the observer
velocity, as can be seen by the fact that the corresponding contours are
tilted. The degeneracy is largest at low redshift and gradually disappears as
the redshift interval grows — this is expected since the two effects act
differently across redshift. As a guide to the eye, we have highlighted the
three relevant panels in light green.

The dashed red lines in Figure 8 indicate (in the relevant panels) the
deceleration dipole measured in the ‘smooth’ simulation using numerical
derivatives at the observer location or the relative velocity of the observer
frame with respect to the matter at that same location. The MCMC analysis
recovers the observer boost very accurately from both redshift windows
(extremely tight orange and red contours), whereas the dipole of the
deceleration is not measured as well, especially when using only data at low
redshifts $z<0.025$.

##  6 Conclusion

Based on the cosmographic framework by Heinesen (2021a) for analysing the
distance–redshift relation in a model-independent way, we have developed a new
approach to constrain the anisotropy in the expansion rate, deceleration, and
higher-order differentials of the distance–redshift relation in the real
Universe that is inhomogeneous on small scales. These new observables
encapsulate information about the inhomogeneities within the local volume over
which the data are taken. The observational selection function introduces an
effective averaging procedure, and as the probed volume increases we see a
trend towards the FLRW limit, as expected in the simulated $\Lambda$CDM model
of our analysis. The departure from FLRW as a function of selection window can
be seen as a new observable that can be predicted for any cosmological model,
opening up a new way to test and constrain different possible scenarios. Our
analysis can be applied directly to real observations.

Here we used synthetic data, generated with $\Lambda$CDM $N$-body simulations
with relativistic ray tracing, to develop our analysis. Since a third-order
Taylor series expansion of the distance–redshift relation (including expansion
rate, deceleration, and jerk) has an error larger than $1\%$ already at
redshift $z=0.2$, we employ a Padé approximant that is better behaved over an
extended redshift range and limits our truncation error below $0.2\%$. We find
that the lowest multipoles in the expansion rate and deceleration can be
detected at high significance with a few thousand sources that sample the full
sky and have a scatter that is dominated by peculiar motion. Under such
optimal measurement conditions, the observed deviation from FLRW is still
significant even at redshift $z=0.15$, i.e. in a volume measuring more than 1
Gpc in diameter. In particular, the variation of the expansion rate can be as
large as $1\%$ across the sky, while the deceleration parameter can have
variations of order $10\%$.

We compare our new analysis to previous work where the multipoles of the
distance–redshift relation are computed directly from the kinematical
properties of the matter congruence at the observer location. This requires a
unique notion of matter congruence that only exists in a Universe where the
overlap of multiple streams of matter (often called ‘shell-crossing’) is
strictly avoided. We therefore run a simulation using initial conditions where
small-scale perturbations were artificially removed by applying a low-pass
filter. We find that the local definition of the expansion rate agrees well
with our new analysis when applied to the same simulation. The agreement is
indeed extremely precise for the monopole and quadrupole of the expansion
rate, where both methods give the same departure from the FLRW limit. The
inferred deceleration parameter and its multipoles, when measured using our
new analysis, do not resemble the local measurement that well. This may be
related to the fact that the latter measurement is based on numerical
derivatives of data, which are prone to numerical issues — such as particle
aliasing when performing a particle-to-mesh projection — especially for higher
derivative orders. We also note, however, that our MCMC analysis is based on a
Gaussian likelihood model, which does not well describe the true scatter of
the data when small-scale structure (otherwise the primary cause of scatter)
has been artificially suppressed.

While finalising this article, an interesting study by Kalbouneh et al. (2024)
appeared that is in many ways complementary to our work. Using a toy model for
a large inhomogeneity in the local Universe that still admits a unique matter
congruence, the anisotropies in the distance–redshift relation are discussed
in the cosmographic framework. Their discussion of the role of the observer
frame is consistent with our approach, showing that it can be included in the
analysis as we both suggest. While they advocate using a Taylor series
expansion in redshift for the distance–redshift relation (similar to earlier
work, including some of our own, e.g. Heinesen, 2021a; Heinesen & Macpherson,
2022), we believe that using the inverse series, i.e. an expansion in angular
diameter distance, has advantages when dealing with real data.111In practice,
there may be subtleties in using the distance as a parameter, since the
distances to astrophysical sources are often not a priori known and must be
fitted together with the cosmology, as is for instance the case in the
standardisation of supernovae of type 1a by Tripp (1998). In particular,
angular diameter distance is a smooth parameter along any null ray and is in
practice strictly monotonic over a significant portion of the local light
cone. The irreducible, intrinsic scatter of the data is instead linked mainly
to Doppler contributions that affect the observed redshift, rendering this
quantity discontinuous and badly non-monotonic. Our approach does not require
a binning in redshift and can effortlessly incorporate negative redshifts
which could be relevant for accurately measuring large observer boosts (it is
of course physically impossible to find negative values of the angular
diameter distance of a source).

In the realistic case where highly nonlinear structures form in the Universe,
the notion of a matter congruence becomes ambiguous. Our approach provides an
observational procedure to construct an effective congruence for a given
survey footprint, allowing us to quantify the expected deviations from the
FLRW limit. Since these deviations could be measured out to very large scales,
they could be used as new probes in the analysis of cosmic large-scale
structure, shedding a new light on the questions of how data should be
interpreted in an inhomogeneous Universe and how exactly the FLRW limit is
approached when taking summary statistics.

###### Acknowledgements.

The simulations used in this work were originally run on the DiRAC DIaL
system, operated by the University of Leicester IT Services, which forms part
of the STFC DiRAC HPC Facility (www.dirac.ac.uk). Reusing existing simulation
data greatly reduces the carbon footprint of our research activity. JA, RD and
MK acknowledge financial support from the Swiss National Science Foundation.
Support for HJM was provided by NASA through the NASA Hubble Fellowship grant
HST-HF2-51514.001-A awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. CC is supported by the UK Science &
Technology Facilities Council Consolidated Grant ST/T000341/1.

##  Appendix A Comparison with the Einstein Toolkit

(a)

(b)

Figure 9: 2D slices at redshift $z=0$ of the density, scalar expansion, shear
scalar, and acceleration scalar (top-left to bottom-right, respectively) for
gevolution data (a) and ET data (b) for runs with identical initial
conditions. The physical box length is $1.28\,h^{-1}\,\mathrm{Gpc}$ and all
modes with wavelength below $200\,h^{-1}\,\mathrm{Mpc}$ were removed from the
initial conditions.

In this Appendix, we consider two simulations with identical initial data
performed with different codes: gevolution and the Einstein
Toolkit222einsteintoolkit.org (ET; Löffler et al., 2012; Zilhão & Löffler,
2013). The purpose of this test is to ensure that our calculations of the
cosmographic parameters using $N$-body simulation data are consistent with a
simulation of a strict fluid congruence in a carefully controlled setup, which
we describe below.

The simulation performed with the ET mimics the setup used in Macpherson et
al. (2019); Macpherson & Heinesen (2021), namely that of a matter-dominated
universe initialised as a linearly-perturbed Einstein-de Sitter (EdS) model.
The initial conditions are set up with
FLRWSolver333https://github.com/hayleyjm/FLRWSolver_public (Macpherson et al.,
2017) as Gaussian random fluctuations following the matter power spectrum at
$z_{\mathrm{ini}}=1000$ from CLASS444http://class-code.net (Blas et al.,
2011). To replicate the same spacetime solution with gevolution we generate
the initial data for the ET simulation — which adopts a fluid approximation
for the matter content — and then translate these to particle displacement
fields and velocities for the $N$-body ensemble of gevolution using a method
described in Adamek et al. (2022).

The matter power spectrum we use to generate the initial data has initial
perturbations smaller than $\sim 200\,{\rm Mpc}/h$ removed, namely, we set
$P(k>k_{\rm max})=0$ for $k_{\rm max}=2\pi/\lambda_{\rm min}$ and
$\lambda_{\rm min}=200\,{\rm Mpc}/h$. Then the power at $k<k_{\rm max}$ is set
according to the output power spectrum from CLASS. We choose this particular
scale cut to ensure that the structures appear at large scales only, and thus
should remain well within the perturbative regime where the fluid description
is valid. In the presence of nonlinear structures, we might expect different
results when comparing simulations performed under a fluid approximation and
with an $N$-body description, especially when sampling small enough scales
such that we resolve shell-crossing, in which case the fluid approximation
breaks down completely.

Both simulations use a mesh of $256^{3}$ points, and the $N$-body simulation
is run with $1024^{3}$ particles to reduce noise from the particle
discretisation. We choose a physical box length of $L=1.28\,\mathrm{Gpc}/h$
for both simulations. The simulations are both initialised in the Newtonian
gauge, and the lapse and shift chosen in the ET simulation will maintain this
gauge so long as the perturbations remain in the linear regime (see Macpherson
et al., 2019, for more details on the gauge used). The relativistic $N$-body
simulation with gevolution uses the Newtonian gauge throughout.

Figure 9 shows two-dimensional slices of the fluid scalars in the final
snapshots of the two simulations run with gevolution (left panels) and the ET
(right panels). For each simulation, we show the density field, expansion
scalar, shear scalar, and acceleration magnitude (top-left to bottom-right on
each side, respectively). While the outputs of the two simulations show a good
overall agreement on large scales, the particle-to-mesh projection used to
construct the fields from the $N$-body ensemble leads to small-scale numerical
artefacts that are particularly noticeable in the bottom-right subplot showing
the acceleration.

The smooth density and velocity fields constructed from gevolution were then
passed through the same post-processing code as for the ET data (mescaline;
see Macpherson et al., 2019). For both simulations, we calculate the effective
cosmological parameters on the slice at redshift $z=0$ (see Macpherson &
Heinesen, 2021, for further details on the calculation itself). We calculate
the effective Hubble, deceleration, jerk, and curvature parameters first for a
set of 100 observers randomly placed in each simulation domain (with
equivalent positions between simulations), with lines of sight drawn in the
direction of $12\times N_{\rm side}^{2}$ HEALPix pixels with $N_{\rm
side}=32$.

(a)

(b)

Figure 10: Sky maps of the effective parameters for the same observer in the
simulation from gevolution (a) and ET (b) with domain size
$1.28\,h^{-1}\,\mathrm{Gpc}$ and all power beneath $200\,h^{-1}\,\mathrm{Mpc}$
removed from the initial conditions. Each parameter is shown normalised by its
EdS counterpart, with the exception of curvature.

Figure 10 shows a set of Mollweide projections of the effective Hubble,
deceleration, curvature, and jerk (top-left to bottom-right, respectively, in
each set of panels) for one observer placed at the same location in the
simulations from gevolution (left panels) and ET (right panels). With the
exception of curvature, all parameters are normalised by their EdS
counterpart. We notice both the anisotropic signatures and amplitudes of the
parameters are very close between the simulations, with the exception of the
jerk parameter (bottom right in both sets of panels), which shows visible
differences in anisotropic signature and a larger difference in amplitude.

Since this analysis is based on numerical derivatives of the fields at the
observer location, we expect that high-order derivatives — like the ones that
enter the calculation of the jerk — may become contaminated by the particle
discretisation, possibly requiring a higher-order particle-to-mesh projection
to mitigate (we used ‘cloud-in-cell’ projection). However, we do not expect
that our MCMC analysis in the main text suffers from similar shortfalls, as no
derivatives of data are taken at all in this case.

Next, we consider a set of 1000 observers (again placed in identical positions
between the two simulations) each with 300 randomly-drawn lines of sight
(lines of sight are unique for each observer but identical between
simulations). For each line of sight we calculate the four effective
cosmological parameters, and average the result over the sky for each
observer. This procedure was done in Macpherson & Heinesen (2021) to provide
an estimate on the observers “measurement” of these parameters, and to assess
how they might vary with observer position.

Figure 11: Sky averages for effective cosmological parameters for 1000
observers in ET (filled circles) and gevolution (crosses). Each point shows
the effective parameters relative to the EdS value for one observer averaged
over 300 randomly drawn lines of sight. Averages over all observers are shown
as dashed and dotted horizontal lines for the ET and gevolution data,
respectively.

Figure 11 shows these sky averages of the four cosmological parameters
(panels; normalised by each EdS counterpart) for all observers in the
simulations from gevolution (crosses) and ET (filled circles). We show the
mean over all observers as a dashed horizontal line for the ET and as a dotted
horizontal line for gevolution. For reference, the $\Lambda$CDM values of the
parameters are shown as black dot-dashed horizontal lines. Broadly, the
variance over observers matches between the two simulations. We note there are
some differences, however, we emphasise that we are comparing two independent
codes which have different numerical errors as well as distinct descriptions
of the matter field itself.

We conclude that the simulation run with gevolution provides a very good
description of the large-scale congruence in this controlled situation in
which we expect it to exist.

## References

  * Abdalla et al. (2022) Abdalla, E., et al. 2022, JHEAp, 34, 49, doi: 10.1016/j.jheap.2022.04.002
  * Adamek et al. (2022) Adamek, J., Calles, J., Montandon, T., Noreña, J., & Stahl, C. 2022, JCAP, 04, 001, doi: 10.1088/1475-7516/2022/04/001
  * Adamek et al. (2016) Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016, Nature Phys., 12, 346, doi: 10.1038/nphys3673
  * Blas et al. (2011) Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034, doi: 10.1088/1475-7516/2011/07/034
  * Capozziello et al. (2020) Capozziello, S., D’Agostino, R., & Luongo, O. 2020, Mon. Not. Roy. Astron. Soc., 494, 2576, doi: 10.1093/mnras/staa871
  * Cattoen & Visser (2007) Cattoen, C., & Visser, M. 2007, Class. Quant. Grav., 24, 5985, doi: 10.1088/0264-9381/24/23/018
  * Clarkson et al. (2011) Clarkson, C., Ellis, G., Larena, J., & Umeh, O. 2011, Rept. Prog. Phys., 74, 112901, doi: 10.1088/0034-4885/74/11/112901
  * Clarkson et al. (2012) Clarkson, C., Ellis, G. F. R., Faltenbacher, A., et al. 2012, Mon. Not. Roy. Astron. Soc., 426, 1121, doi: 10.1111/j.1365-2966.2012.21750.x
  * Clarkson & Umeh (2011) Clarkson, C., & Umeh, O. 2011, Class. Quant. Grav., 28, 164010, doi: 10.1088/0264-9381/28/16/164010
  * Coates et al. (2021) Coates, L., Adamek, J., Bull, P., Guandalin, C., & Clarkson, C. 2021, Mon. Not. Roy. Astron. Soc., 504, 3534, doi: 10.1093/mnras/stab1076
  * Cowell et al. (2023) Cowell, J. A., Dhawan, S., & Macpherson, H. J. 2023, Mon. Not. Roy. Astron. Soc., 526, 1482, doi: 10.1093/mnras/stad2788
  * Demianski et al. (2017) Demianski, M., Piedipalumbo, E., Sawant, D., & Amati, L. 2017, Astron. Astrophys., 598, A113, doi: 10.1051/0004-6361/201628911
  * Dhawan et al. (2023) Dhawan, S., Borderies, A., Macpherson, H. J., & Heinesen, A. 2023, Mon. Not. Roy. Astron. Soc., 519, 4841, doi: 10.1093/mnras/stac3812
  * Ellis et al. (1985) Ellis, G. F. R., Nel, S. D., Maartens, R., Stoeger, W. R., & Whitman, A. P. 1985, Phys. Rep., 124, 315, doi: 10.1016/0370-1573(85)90030-4
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, Astrophys. J., 622, 759, doi: 10.1086/427976
  * Gruber & Luongo (2014) Gruber, C., & Luongo, O. 2014, Phys. Rev. D, 89, 103506, doi: 10.1103/PhysRevD.89.103506
  * Heinesen (2021a) Heinesen, A. 2021a, JCAP, 05, 008, doi: 10.1088/1475-7516/2021/05/008
  * Heinesen (2021b) —. 2021b, Phys. Rev. D, 104, 123527, doi: 10.1103/PhysRevD.104.123527
  * Heinesen & Macpherson (2022) Heinesen, A., & Macpherson, H. J. 2022, JCAP, 03, 057, doi: 10.1088/1475-7516/2022/03/057
  * Kalbouneh et al. (2024) Kalbouneh, B., Marinoni, C., & Maartens, R. 2024.  https://arxiv.org/abs/2401.12291
  * Kristian & Sachs (1966) Kristian, J., & Sachs, R. K. 1966, Astrophys. J., 143, 379, doi: 10.1086/148522
  * Lepori et al. (2020) Lepori, F., Adamek, J., Durrer, R., Clarkson, C., & Coates, L. 2020, MNRAS, 497, 2078, doi: 10.1093/mnras/staa2024
  * Löffler et al. (2012) Löffler, F., et al. 2012, Class. Quant. Grav., 29, 115001, doi: 10.1088/0264-9381/29/11/115001
  * Maartens et al. (2023) Maartens, R., Santiago, J., Clarkson, C., Kalbouneh, B., & Marinoni, C. 2023.  https://arxiv.org/abs/2312.09875
  * MacCallum & Ellis (1970) MacCallum, M. A. H., & Ellis, G. F. R. 1970, Comm. Math. Phys., 19, 31, doi: 10.1007/BF01645496
  * Macpherson & Heinesen (2021) Macpherson, H. J., & Heinesen, A. 2021, Phys. Rev. D, 104, 023525, doi: 10.1103/PhysRevD.104.023525
  * Macpherson et al. (2017) Macpherson, H. J., Lasky, P. D., & Price, D. J. 2017, Phys. Rev. D, 95, 064028, doi: 10.1103/PhysRevD.95.064028
  * Macpherson et al. (2019) Macpherson, H. J., Price, D. J., & Lasky, P. D. 2019, Phys. Rev. D, 99, 063522, doi: 10.1103/PhysRevD.99.063522
  * Nesseris & García-Bellido (2013) Nesseris, S., & García-Bellido, J. 2013, Phys. Rev. D, 88, 063521, doi: 10.1103/PhysRevD.88.063521
  * Perivolaropoulos & Skara (2022) Perivolaropoulos, L., & Skara, F. 2022, New Astron. Rev., 95, 101659, doi: 10.1016/j.newar.2022.101659
  * Sapone et al. (2014) Sapone, D., Majerotto, E., & Nesseris, S. 2014, Phys. Rev. D, 90, 023012, doi: 10.1103/PhysRevD.90.023012
  * Tripp (1998) Tripp, R. 1998, Astron. Astrophys., 331, 815
  * Umeh (2013) Umeh, O. 2013, PhD thesis, University of Cape Town, Faculty of Science, Department of Mathematics and Applied Mathematics
  * Vitagliano et al. (2010) Vitagliano, V., Xia, J.-Q., Liberati, S., & Viel, M. 2010, JCAP, 03, 005, doi: 10.1088/1475-7516/2010/03/005
  * Weinberg (1972) Weinberg, S. 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity (New York: John Wiley and Sons)
  * Zilhão & Löffler (2013) Zilhão, M., & Löffler, F. 2013, Int. J. Mod. Phys. A, 28, 1340014, doi: 10.1142/S0217751X13400149