∎

11institutetext: Thomas Buchert1,4,‡ 22institutetext: Henk van Elst2,3
33institutetext: Asta Heinesen1 44institutetext: 1 Univ Lyon, Ens de Lyon,
Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F–69007
Lyon, France
2 parcIT GmbH, Erftstraße 15, D–50672 Köln, Germany
3 Cosmology and Gravity Group (CAGG), Department of Mathematics and Applied
Mathematics, University of Cape Town, ZA–7701 Rondebosch, South Africa
4 Fakultät für Physik, Universität Bielefeld, Postfach 100131, D–33501
Bielefeld, Germany
$\ddagger$ corresponding author.

44email<EMAIL_ADDRESS>$\cdot$<EMAIL_ADDRESS>$\cdot$
<EMAIL_ADDRESS>
# The averaging problem on the past null cone
in inhomogeneous dust cosmologies ††thanks: Work supported by ERC advanced
Grant 740021–ARTHUS.

Thomas Buchert and Henk van Elst and Asta Heinesen

(Received: 23 February 2022 / Accepted: 8 December 2022)

###### Abstract

Cosmological models typically neglect the complicated nature of the spacetime
manifold at small scales in order to hypothesize idealized general
relativistic solutions for describing the average dynamics of the Universe.
Although these solutions are remarkably successful in accounting for data,
they introduce a number of puzzles in cosmology, and their foundational
assumptions are therefore important to test. In this paper, we go beyond the
usual assumptions in cosmology and propose a formalism for averaging the local
general relativistic spacetime on an observer’s past null cone: we formulate
average properties of light fronts as they propagate from a cosmological
emitter to an observer. The energy-momentum tensor is composed of an
irrotational dust source and a cosmological constant – the same components as
in the $\Lambda$CDM model for late cosmic times – but the metric solution is
not _a priori_ constrained to be locally homogeneous or isotropic. This
generally makes the large-scale dynamics depart from that of a simple
Friedmann–Lemaître–Robertson–Walker solution through ‘backreaction’ effects.
Our formalism quantifies such departures through a fully covariant system of
area-averaged equations on the light fronts propagating towards an observer,
which can be directly applied to analytical and numerical investigations of
cosmic observables. For this purpose, we formulate light front averages of
observable quantities, including the effective angular diameter distance and
the cosmological redshift drift and we also discuss the backreaction effects
for these observables.

###### Keywords:

Relativistic cosmology Spacetime foliations Light propagation Cosmological
backreaction Dark Universe

††journal: GRG

###### Contents

  1. 1 Introduction
  2. 2 Local spacetime setting
    1. 2.1 Field equations and matter content
    2. 2.2 Past null cone
    3. 2.3 Foliation of the past null cone and light fronts
    4. 2.4 Propagation equations along the past null cone
    5. 2.5 Constraint equations on the light fronts
    6. 2.6 Photon conservation, distance measures, and redshift drift
  3. 3 Area-averaging over light fronts
    1. 3.1 Area-averaging operation and evolution equations for area-averaged variables
    2. 3.2 Average conservation laws and integral-geometric measures
    3. 3.3 Area-averaged observables
    4. 3.4 FLRW limit
  4. 4 Discussion
  5. A Light front-adapted spacetime metric

##  1 Introduction

Observational cosmology primarily relies on information that propagates along
our past null cone; cf. Ref. krisac1966 . Observational data is most commonly
interpreted within the Friedmann–Lemaître–Robertson–Walker (FLRW) class of
spacetime metrics. Perturbation theory at this background cosmology is then
employed to refine the description of the propagation of null signals along
the past null cone in an ‘almost FLRW’ cosmology that takes inhomogeneities
into account (see, e.g., Refs. bildhauer:backreaction2 ; sasaki1 ; sasaki3 ;
vanderveldetal ). The $\Lambda$CDM model (Cold Dark Matter including a
positive cosmological constant $\Lambda$) is a particular general relativistic
FLRW metric model, with ordinary matter, radiation, dark matter (CDM), and
dark energy (modelled with $\Lambda$) as energy-momentum sources. The
$\Lambda$CDM paradigm has proven successful in a largely consistent
interpretation of observational data; however, with a number of ‘tensions’
remaining FLRWtensions ; kids ; Riess2019 ; hubbletension1 ; hubbletension2 ;
Perivolaropoulos:2021jda . Importantly, the success of the $\Lambda$CDM
framework comes at the price of introducing dark components to the
cosmological energy budget planck .

The specification of a realistic inhomogeneous spacetime metric for a lumpy
model of the Universe and the description of realistic local properties for
the propagation of null signals is an involved task. Cosmological observations
involve the study of many sources and the propagation of null signals over
large spatial distances. We might thus aim at formulating average dynamical
equations – suitable for obtaining effective laws for the propagation of null
signals in a cosmological context – where exact knowledge of the local
spacetime metric might not be required. Averaging of Einstein’s field
equations in general relativity has been formulated in the context of a
$3$-dimensional volume-averaging of scalar-valued variables111Averaging over
tensor-valued variables introduces ambiguities in the volume-averaging
procedure (see Refs. ellis ; ellisbuchert ; buchert:review ; klingon ; singh2
and references therein). buchert:grgdust ; buchert:grgfluid . In this
approach, the spacetime metric is not specified and no approximation is
employed to restrict inhomogeneities. Instead, one concentrates on integral
properties of local spacetime variables such as the rest mass density, the
expansion rate, the shear rate and the spatial Ricci curvature scalar. The
resulting volume-averaged Einstein’s field equations can be viewed as balance
equations for effective macroscopic spacetime variables. The system of
equations governing the macroscopic variables is in general not closed and
additional constraints must be imposed in order to solve the system. (A
comprehensive discussion of the closure problem in $N$ dimensions, in
particular emphasizing topological constraints, may be found in Ref. GBC .)

In such average Universe descriptions, the structure of the underlying local
inhomogeneous distributions of the matter content and the spatial geometry
manifests itself globally through correction terms to Friedmann’s equations.
Such terms will be non-vanishing in general – also in cosmologies with a
notion of statistical spatial homogeneity and isotropy – and represent the
deviance of the ‘monopole state’ of a model universe with structure, i.e., an
appropriate lowest-order and spatially uniform solution, from the strictly
spatially homogeneous and isotropic FLRW solutions. The monopole state can,
nevertheless, be mapped to an effective FLRW cosmology by interpreting the
correction terms as effective source terms. For instance, volume-averaged
fluctuation terms have been discussed as geometrical candidates of kinematical
dark energy and dark matter sources; e.g. Refs. buchert:jgrg ; palle ; morphon
; David1 ; buchert:darkenergy ; kolbetal ; rasanen ; lischwarz ;
buchert:review ; Bolejko1 ; David2 ; wiegandbuchert ; RZA_2 ; roy:instability
; sussman1 ; Buchert2011 ; Rasanenreview ; buchertrasanen ; RZA_2 ; sussman2 ;
David3 ; zimdahl:LTB ; Bolejko2 ; quentin:darkmatter ; hubbletension3 .

Much work has been done in describing the past null cones of observers in
various inhomogeneous and anisotropic model universes (e.g. Refs. newpen1962 ;
dyeroe74 ; dyeroe75 ; weinberg ; linsenbuch ; maamat1994 ; mustapha1 ;
ellis:dahlem ; kantowski3 ; kantowski2 ; kantowski1 ; demetal03 ; pyne ;
Gourgoulhon2006 ; ruth:luminosity ; hellaby:mass ; luhellaby ; dominik:average
), including coordinate-independent considerations within gauge invariant
perturbation theory Marco:perturbations ; YooDurrer ; Yoo:2019qsl , and the
distortion of the average distance–redshift relation relative to that of a
reference FLRW cosmology in specific spatially inhomogeneous model universes
mustapha2 ; SikoraGlod ; LTBluminosity1 ; LTBluminosity2 ; LTBluminosity3 ;
Fleury:2013sna ; Fleury:2013uqa ; Koksbang:2017arw ; Adamek:2018rru . General
considerations on the propagation of null signals in statistically spatially
homogeneous and isotropic cosmologies with slowly evolving structure have been
carried out lightpropndust ; lightpropngeneral . Such considerations allow for
relating leading-order observable quantities to volume-averaged variables on
spacelike $3$-surfaces under the assumption that null signals sample the
Universe fairly in volume. Covariant null cone averaging formalisms suitable
for averaging on light fronts have been developed in Refs. lightconeav1 ;
Fanizza2019pfp , but have mainly been employed for the purpose of constraining
errors in the determination of FLRW background parameters in a $\Lambda$CDM
cosmology lightconeav2 ; BenDayan:2012 ; BenDayan:2013gc . Formalism that can
be used for model-independent analysis of standardizable objects and cosmic
drift effects, while making no assumptions about the form of the cosmological
spacetime metric, has been presented in Refs. Umeh:2013UCT ; Heinesen:2020bej
; Heinesen:2021nrc ; Heinesen:2021qnl  – see also Refs. Nez ;
Korzynski:2017nas ; Grasso:2018mei ; Korzynski:2021aqk ; carfora  for other
interesting formalisms which could be suitable for model-independent data
analysis.

In this paper, we employ an averaging operation similar to that introduced in
Refs. buchert:grgdust ; buchert:grgfluid ,222The same averaging operation
might be formulated in the notation of Refs. lightconeav1 ; Fanizza2019pfp .
but with the averaging domain adapted to an observer’s past null cone. More
precisely, we shall propagate $2$-dimensional inhomogeneous light fronts on
which we define area-averaged physical variables (of, e.g., matter density and
optical scalars), and invoke Einstein’s field equations to obtain propagation
equations for the same area-averaged variables. We do this without specifying
the spacetime metric. We shall characterize average cosmological distance
relations in terms of null cone averages over light fronts. We shall constrain
the local setting to that of an irrotational dust cosmology. This
approximation comes with limitations: on scales where gravitationally bound
structures are forming, the irrotationality assumption for matter breaks down.
At such scales, velocity dispersion arising from the internal motion effects
within the structures must generally also be included through an effective
anisotropic pressure term that violates the dust approximation. In the early
history of the Universe, where radiation constitutes a significant fraction of
the total energy density, pressure must also be included for an accurate
description. The framework developed here is thus applicable for scales above
the largest gravitationally virialized structures and for cosmic times well
after the epoch of radiation and matter equality. In order to include
modelling of smaller length scales and earlier cosmic times, the framework
must be generalized in terms of the assumed matter content, which may be done
by following similar approaches as in Ref. buchertetal:generalfluid , see also
buchertetal:foliations . We furthermore consider the idealized case where
caustics in the geodesic null congruence constituting an observer’s past null
cone can be ignored. This is consistent with the large-scale modelling of the
matter content of the Universe described above, where we neglect the physics
of strong gravitational lenses. For treatments including caustics, see, e.g.,
Refs. ellissolomons ; mustapha3 .

The general framework investigated in this paper offers the possibility of
employing non-perturbative and background-free approximations to describe the
effect of inhomogeneities on cosmological measurements.

Notation and conventions: We use units in which $c=G/c^{2}=1$. Greek letters
$\mu,\nu,\ldots$ label spacetime indices in a general basis. Summation over
repeated indices is understood. The signature of the spacetime metric,
$g_{\mu\nu}$, is taken to be $(-+++)$ dreimaennerbuch .

##  2 Local spacetime setting

Here we describe the assumptions made for the local cosmological spacetime. In
section 2.1 we describe the assumed matter content in the form of an
irrotational dust source. In section 2.2 we describe the observer’s past null
cone and in section 2.3 we introduce the foliation of the past null cone into
light fronts that will be used throughout this paper. In sections 2.4 and 2.5
we derive propagation equations and constraint equations, respectively, for
capturing the dynamics of relevant variables for the cosmological spacetime
and null signals travelling along an observer’s past null cone. Finally, in
section 2.6 we consider photon number conservation and effective cosmological
distance measures.

###  2.1 Field equations and matter content

We consider a cosmological spacetime which is dynamically described through
Einstein’s field equations ein1915 ; ein1917 ,

$R_{\mu\nu}-\frac{1}{2}\,g_{\mu\nu}R+\Lambda g_{\mu\nu}=8\pi T_{\mu\nu}\ ,$
(1)

where $R_{\mu\nu}$ is the $4$-Ricci curvature of the spacetime, $R$ is its
trace, and $\Lambda$ denotes the cosmological constant. We concentrate on a
matter-dominated epoch of the cosmological evolution and restrict the energy-
momentum tensor $T_{\mu\nu}$ to that of an irrotational continuum of dust,

$\displaystyle T_{\mu\nu}=\varrho u_{\mu}u_{\nu}\
;\quad\omega_{\mu\nu}:=h_{[\mu}{}^{\alpha}h_{\nu]}{}^{\beta}\nabla_{\alpha}u_{\beta}=0\
,$ (2)

where $u^{\mu}$ is the future-directed $4$-velocity of the matter congruence,
$\varrho$ is the matter’s rest mass density, $\nabla_{\mu}$ is the covariant
derivative, and square brackets denote anti-symmetrization in the indices
involved. We assume that the irrotational continuum of dust is transparent to
propagating null signals. From Eqs. (2) and the conservation of energy-
momentum it follows that $u^{\mu}$ is geodesic and can be expressed as the
gradient of a scalar-valued variable,

$\displaystyle u^{\nu}\nabla_{\nu}u^{\mu}=0\ ;\quad
u_{\mu}=-\,\partial_{\mu}\tau\ ,$ (3)

where $\tau$ can be interpreted as a proper time function along each
individual matter world line, since $u^{\rho}\partial_{\rho}\tau=1$. The
covariant $(1+3)$-decomposition of the covariant derivative of the
$4$-velocity of the matter congruence thus results in (cf. Refs. ehl1961 ;
ell1971 ; hve1996 ; wainwrightellis ; hveetal1997 ; hve1998 ; ellhve1999 ;
elletal2012 ; roy:1+3 )

$\nabla_{\mu}u_{\nu}=\Theta_{\mu\nu}=\frac{1}{3}\,\theta
h_{\mu\nu}+\sigma_{\mu\nu}\ ,$ (4)

where $\theta$ is the isotropic expansion rate describing the volume expansion
of the matter congruence, and $\sigma_{\mu\nu}$ is the volume shear rate
describing the volume-preserving deformation of the congruence. The covariant
derivative of the $4$-velocity is equal to its own projection
$\Theta_{\mu\nu}:=h_{\mu}{}^{\alpha}h_{\nu}{}^{\beta}\nabla_{\alpha}u_{\beta}$
because of the vanishing of the $4$-acceleration of $u^{\mu}$. The spatial
projection tensor

$h_{\mu}{}^{\nu}:=\delta_{\mu}{}^{\nu}+u_{\mu}u^{\nu}\ ,$ (5)

is orthogonal to the irrotational matter $4$-velocity field and plays the role
of the metric tensor on the spacelike $3$-surfaces orthogonal to $u^{\mu}$.

###  2.2 Past null cone

We now assume the presence of a non-gravitating, geodesic and irrotational
null congruence pervading the cosmological spacetime such that its future-
directed $4$-momentum $k^{\mu}$ obeys Seitz:1994xf

$\displaystyle k_{\mu}k^{\mu}=0\ ;\quad k^{\nu}\nabla_{\nu}k^{\mu}=0\ ;\quad
k_{[\mu}\nabla_{\nu}k_{\alpha]}=0\ .$ (6)

This implies that $k^{\mu}$ is proportional to the gradient of a scalar-valued
variable,

$\displaystyle k_{\mu}=\eta\,\partial_{\mu}V\ ;\qquad
k^{\rho}\partial_{\rho}V=0\ ;\quad k^{\rho}\partial_{\rho}\eta=0\ ,$ (7)

where the conditions $\eta>0$ and $n^{\rho}\partial_{\rho}V<0$ (with $n^{\mu}$
a future-directed timelike vector field) ensure the future-directed nature of
$k^{\mu}$. The function $V$ defines the direction of the $4$-momentum
$k^{\mu}$ and $\eta$ defines its normalization. The preservation of $V$ along
each individual ray of the geodesic null congruence follows from the vanishing
norm of $k^{\mu}$, and the preservation of $\eta$ follows from the affine
geodesic equation. Spacetime domains singled out by $\\{V=\text{constant}\\}$
are $3$-dimensional null surfaces.

As a special case of null surfaces we consider those associated with the past
null cones of an observer. We consider the case where the timelike world line,
$\gamma$, of this observer belongs to the matter congruence.333We might in
principle choose any world line to initialize a past null cone, but we shall
often be interested in observers comoving with the matter in the cosmological
spacetime. For each point $P$ of $\gamma$, we consider incoming geodesic null
rays of all spatial directions orthogonal to $u^{\mu}$ as initial conditions
for the direction of the $4$-momentum $k^{\mu}$. We extend this two-parameter
family of null directions away from $P$ via the geodesic equation to construct
the past null cone ${\cal C}^{-}(P)$ at $P$. From
$u^{\rho}\partial_{\rho}V<0$, we have that the scalar $V$ uniquely labels the
one-parameter family of past null cones ${\cal C}^{-}(P)$ along $\gamma$.
Since $\eta$ is a constant over each null cone, we have that $\eta$ is a
function of $V$, and can be absorbed into a redefinition of $V$. We thus set
$\eta=1$ without loss of generality. The normalization of the $4$-momentum
$k^{\mu}$ might for instance be chosen by fixing the energy function,

$E:=-u^{\mu}k_{\mu}=-u^{\rho}\partial_{\rho}V\ ,$ (8)

along $\gamma$. This choice of normalization specifies the frequency of
monochromatic null signals measured at the observer’s telescope. We might for
instance consider the normalization
$-u^{\mu}k_{\mu}\lvert_{\gamma}=\text{constant}$ that corresponds to the
situation where the central observer performs measurements at the same
radiation frequency at different instances of proper time.

The cosmological redshift, $z$, of a luminous astrophysical source comoving
with the matter congruence and located on the observer’s past null cone, as
measured at the point of observation $P$ (“here and now”), is given by (cf.
Ref. ellhve1999 )

$(1+z):=\frac{-u^{\mu}k_{\mu}}{-u^{\mu}k_{\mu}(P)}=\frac{E}{E(P)}\ ,$ (9)

where $E(P)$ denotes the energy of the incoming null ray as measured by the
observer at the point of observation $P$. The redshift function is independent
of the normalization procedure chosen for $k^{\mu}$.

We may decompose the $4$-momentum $k^{\mu}$ associated with the incoming
geodesic null congruence in terms of the $4$-velocity $u^{\mu}$ of the matter
congruence and a spatial unit vector $e^{\mu}$ orthogonal to $u^{\mu}$
according to

$k^{\mu}=E(u^{\mu}-e^{\mu})\ ,$ (10a) where $e^{\mu}$, as evaluated by the
observer at $P$, is multi-valued and takes values in the two-parameter family
of all possible spatial directions of pointings of the observer’s telescope.
For later reference, we introduce an auxiliary future-directed null
$4$-momentum, $l^{\mu}=E(u^{\mu}+e^{\mu})\ .$ (10b)

The null $4$-momentum $l^{\mu}$ points in the opposite spatial direction as
$k^{\mu}$, when viewed in the rest frame of the matter congruence. Note that
$l^{\mu}$ is not necessarily a geodesic null congruence (cf. footnote 5 for
further properties). We refer to geodesic null rays generated by $k^{\mu}$ as
incoming null rays and null rays generated by $l^{\mu}$ as outgoing null rays.
The field $l^{\mu}$ does not incorporate any new physical information about
the cosmological spacetime, but will be convenient for formulating constraint
equations in the covariant $(1+1+2)$-decomposition of the cosmological
spacetime considered later in this analysis (see Sec. 2.5).

###  2.3 Foliation of the past null cone and light fronts

In order to describe the evolution along a fixed past null cone ${\cal
C}^{-}(P):\\!\\{V=\text{constant}\\}$, we consider its foliation into a union
of spatial $2$-surfaces that we refer to as light fronts. One way to construct
light fronts is to consider the intersection of the past null cone with a
family of $3$-surfaces (timelike, spacelike or lightlike). An obvious choice
of an intersecting foliation is that defined from the irrotational
$4$-velocity $u^{\mu}$ of the matter congruence, so that the fronts of
intersection are simultaneously orthogonal to the null congruence and to the
matter congruence.

The resulting light fronts – also referred to as screen spaces, cf. Refs.
sac1961 ; wal1984  – are constant-level surfaces of $\tau$ and $V$. The local
screen basis propagated along the geodesics null rays is usually referred to
as Sachs basis sac1961 ; per2004 ; joretal1961 . In this setting, the full
observer’s past null cone can be considered to be the union of $2$-dimensional
light fronts with varying values of $\tau$.444Another option is to consider
level surfaces of constant affine parameter $\lambda$ of the geodesic null
congruence defined from the propagation requirement
$k^{\rho}\partial_{\rho}\lambda=1$. In order to uniquely define $\lambda$ as a
spacetime function, we must specify initial conditions. We might, for
instance, require setting $\lambda\lvert_{\gamma}=0$. From this it follows
immediately that $u^{\rho}\partial_{\rho}\lambda\lvert_{\gamma}=0$, and the
gradient of $\lambda$ is thus spacelike in the vicinity of $\gamma$, and
$\lambda=\text{constant}$-level surfaces define timelike cylinders in the same
vicinity. Far away from the vertex, the use of $\lambda$ as a meaningful
foliation scalar must be carefully re-assessed. We define the projection
tensor onto the light fronts:

$p_{\mu}{}^{\nu}:=\delta_{\mu}{}^{\nu}+u_{\mu}u^{\nu}-e_{\mu}e^{\nu}\ .$ (11)

We have
$p_{\mu}{}^{\nu}u_{\nu}=p_{\mu}{}^{\nu}e_{\nu}=p_{\mu}{}^{\nu}k_{\nu}=p_{\mu}{}^{\nu}l_{\nu}=0$,
$p_{\mu}{}^{\alpha}p_{\alpha}{}^{\nu}=p_{\mu}{}^{\nu}$, and
$p_{\mu}{}^{\mu}=2$. The tensor $p_{\mu}{}^{\nu}$ is thus orthogonal to the
space spanned by $k^{\mu}$ and $u^{\mu}$ while it acts as the metric tensor
for tensorial fields intrinsic to the screen space. The surface density on the
light fronts is given by

$\epsilon_{\mu\nu}:=u^{\alpha}p_{\mu}{}^{\beta}p_{\nu}{}^{\gamma}e^{\delta}\epsilon_{\alpha\beta\gamma\delta}=\epsilon_{[\mu\nu]}\
,$ (12)

where $\epsilon_{\alpha\beta\gamma\delta}$ is the spacetime permutation
tensor, equal to $\sqrt{-g}$ for even and $-\sqrt{-g}$ for odd permutations of
$0,1,2,3$, and $g$ is the determinant of the cosmological spacetime metric.
The surface density has the following properties:
$\epsilon_{\mu\nu}u^{\nu}=\epsilon_{\mu\nu}e^{\nu}=\epsilon_{\mu\nu}k^{\nu}=\epsilon_{\mu\nu}l^{\nu}=0$,
$\epsilon_{\mu\nu}\epsilon^{\alpha\beta}=2!\,p_{[\mu}{}^{\alpha}p_{\nu]}{}^{\beta}$
and $\epsilon_{\mu\nu}\epsilon^{\mu\nu}=2$. The area expansion tensor
associated with the light fronts can be decomposed as sac1961 ; joretal1961

$\displaystyle\hat{\Theta}_{\mu\nu}:=p_{\mu}{}^{\alpha}p_{\nu}{}^{\beta}\nabla_{\alpha}k_{\beta}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\,\hat{\theta}p_{\mu\nu}+\hat{\sigma}_{\mu\nu}\ ,$
(13a) where $\hat{\theta}:=p^{\mu\nu}\nabla_{\mu}k_{\nu}$ is the dimensionless
area expansion rate of a geodesic null ray bundle and
$\hat{\sigma}_{\mu\nu}:=p_{\mu}{}^{\alpha}p_{\nu}{}^{\beta}\nabla_{\alpha}k_{\beta}-\frac{1}{2}\hat{\theta}p_{\mu\nu}$
is the dimensionless area-preserving deformation or ‘shearing’ rate of the
same geodesic null ray bundle. The anti-symmetric part of the deformation
vanishes by construction via the requirement of irrotationality; see the last
condition in Eq. (6). We define an analogous area expansion tensor for the
auxiliary null $4$-momentum $l^{\mu}$ by
$\displaystyle\tilde{\Theta}_{\mu\nu}:=p_{\mu}{}^{\alpha}p_{\nu}{}^{\beta}\nabla_{\alpha}l_{\beta}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\,\tilde{\theta}p_{\mu\nu}+\tilde{\sigma}_{\mu\nu}\
,$ (13b) where we have labelled the kinematic variables associated with
$k^{\mu}$ by a hat and the kinematic variables associated with $l^{\mu}$ by a
tilde.555 As a result of the fact that $l^{\mu}$ does not in general generate
a geodesic null congruence, it does not in general satisfy the condition (6).
For the same reason, $l^{\mu}$ is not necessarily hypersurface-forming.
However, it is irrotational through its definition as a linear combination of
irrotational vector fields, Eq. (10b), after projection onto the screen space
normal to $k^{\mu}$ and $u^{\mu}$: $p_{\mu}^{\ \alpha}\;p_{\nu}^{\
\beta}\nabla_{[\alpha}l_{\beta]}=0$.

It shall furthermore be useful to also define the analogous expansion tensor
and kinematic variables associated with $e^{\mu}$ (labelled by a bar):

$\overline{\Theta}_{\mu\nu}=p_{\mu}{}^{\alpha}p_{\nu}{}^{\beta}\nabla_{\alpha}e_{\beta}=\frac{1}{2}\,\overline{\theta}\,p_{\mu\nu}+\overline{\sigma}_{\mu\nu}\
.$ (13c)

The trace component of $\overline{\Theta}_{\mu\nu}$ is given by
$\overline{\theta}:=p^{\mu\nu}\nabla_{\mu}e_{\nu}=\nabla_{\mu}e^{\mu}$, which
follows from the orthonormality of $u^{\mu}$ and $e^{\mu}$ and the geodesic
nature of $u^{\mu}$. The symmetric tracefree part of
$\overline{\Theta}_{\mu\nu}$ is
$\overline{\sigma}_{\mu\nu}:=p_{(\mu}{}^{\alpha}p_{\nu)}{}^{\beta}\nabla_{\alpha}e_{\beta}-\frac{1}{2}\,\overline{\theta}\,p_{\mu\nu}$.
The anti-symmetric component
$p_{[\mu}{}^{\alpha}p_{\nu]}{}^{\beta}\nabla_{\alpha}e_{\beta}$ vanishes due
to the requirement of irrotationality for both $u^{\mu}$ and $k^{\mu}$.

The evolution of the multi-valued spatial direction vector $e^{\mu}$ along the
geodesic null congruence generated by $k^{\mu}$ can be formulated as
Heinesen:2020bej

$\displaystyle E^{-1}k^{\rho}\nabla_{\rho}e^{\mu}=-E^{-1}\alpha
k^{\mu}-\left(\frac{1}{3}\,\theta
h_{\nu}{}^{\mu}+\sigma_{\nu}{}^{\mu}\right)e^{\nu}\ ,$ (14)

where we have made use of one of the two area-adapted longitudinal expansion
rate variables for the matter congruence,

$\alpha:=\frac{1}{3}\,\theta-(p^{\mu\nu}\sigma_{\mu\nu})\
,\qquad\beta:=\frac{1}{3}\,\theta+\frac{1}{2}\,(p^{\mu\nu}\sigma_{\mu\nu})\ ,$
(15)

which were introduced in Ref. steell1968 . In the present irrotational dust
cosmology setup the variable $\alpha$ determines the logarithmic rate of
change of the energy $E$ along the geodesic null rays, see Sec. 2.4, and
replaces the Hubble parameter in FLRW cosmology in the general cosmographic
representation of the luminosity distance; see Refs. Umeh:2013UCT ;
Clarkson:2011uk ; Heinesen:2020bej . Using Eq. (14), we can write the
‘acceleration’ vector of $e^{\mu}$ in the following way:

$\displaystyle
h_{\nu}{}^{\mu}e^{\alpha}\nabla_{\alpha}e^{\nu}=\kappa^{\mu}+p_{\nu}{}^{\mu}e^{\rho}\sigma_{\rho}{}^{\nu}\
,$ (16)

where $\kappa^{\mu}:=p_{\nu}{}^{\mu}u^{\rho}\nabla_{\rho}e^{\nu}$ is the
position drift describing the angular drift of the source on the observer’s
sky as seen in a non-rotating (Fermi-propagated) reference frame
Korzynski:2017nas . The acceleration vector (16) involves the volume shear
rate and the angular drift of the source on the observer’s sky and can thus be
thought of as a measure of the violation of isotropy around the observer’s
position. The acceleration vector (16) and $\kappa^{\mu}$ are both constructed
such that they are tangential to the $2$-dimensional screen space.

The dimensionless kinematic variables $\hat{\theta}$, $\hat{\sigma}_{\mu\nu}$,
$\tilde{\theta}$ and $\tilde{\sigma}_{\mu\nu}$ associated with the null
$4$-momenta $k^{\mu}$ and $l^{\mu}$ can be expressed algebraically in terms of
the kinematic variables of $u^{\mu}$ and $e^{\mu}$, and the energy $E$ of the
incoming null rays as follows:

$\displaystyle\frac{\hat{\theta}}{E}=2\beta-\overline{\theta}\
,\qquad\frac{\hat{\sigma}_{\mu\nu}}{E}=p_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}\sigma_{\alpha\beta}-\overline{\sigma}_{\mu\nu}\
,$ (17a) $\displaystyle\frac{\tilde{\theta}}{E}=2\beta+\overline{\theta}\
,\qquad\frac{\tilde{\sigma}_{\mu\nu}}{E}=p_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}\sigma_{\alpha\beta}+\overline{\sigma}_{\mu\nu}\
;$ (17b)

this result is obtained from the covariant decompositions (10a) and (10b), and
the definitions of the kinematic variables introduced in Eqs. (4), (13a),
(13b) and (13c). Angular brackets single out the symmetric and tracefree part
of the tensor in the involved indices.

The null kinematic variables $\hat{\theta}$, $\hat{\sigma}_{\mu\nu}$,
$\tilde{\theta}$ and $\tilde{\sigma}_{\mu\nu}$ inherit the kinematics of
$u^{\mu}$ and $e^{\mu}$ in the screen space. The longitudinal component of the
shear rate, $(p^{\mu\nu}\sigma_{\mu\nu})$, enters in the expansion rates
$\hat{\theta}$ and $\tilde{\theta}$ via $\beta$. In general, the kinematic
variables of $k^{\mu}$ and $l^{\mu}$ differ. This difference can be assigned
to the focusing of the null congruences towards the observer, and is for
instance present for radially propagating geodesic null congruences in FLRW
cosmology.

We may eliminate the kinematic variables of the matter congruence in Eqs.
(17a) and (17b) to obtain the useful relations:

$\tilde{\theta}=\hat{\theta}+2E\overline{\theta}\
,\qquad\quad\tilde{\sigma}_{\mu\nu}=\hat{\sigma}_{\mu\nu}+2E\overline{\sigma}_{\mu\nu}\
.$ (18)

The function $\overline{\theta}$ quantifies the departure of the
$2$-dimensional screen space from a minimal surface area666It is a well-known
result from calculus of variations in Riemannian geometry that surfaces
minimizing the area measure locally have zero trace of the extrinsic curvature
scalar of the embedding. within the spacelike $3$-surfaces orthogonal to
$u^{\mu}$. When $k^{\mu}$ is a generator of a past null cone, which is
narrowing towards the singularity at the vertex of the observer, we expect
$\overline{\theta}$ to be dominantly negative.777Since causal lines can only
leave a past null cone (and not enter), we indeed expect a negative
contribution to the overall expansion rate of the screen space from the drift
of the screen space boundaries relative to the matter congruence: the screen
space is sampling the cross section of fewer fluid elements as the vertex of
the past null cone is approached. However, local differential expansion of the
dust matter congruence and the spatial fluctuations in the rest mass density
$\varrho$ could potentially compensate this tendency locally.

###  2.4 Propagation equations along the past null cone

We now formulate the evolution equations along the observer’s past null cone
${\cal C}^{-}(P):\\!\\{V=\text{constant}\\}$ for the variables of our main
interest.

The propagation equations for the optical scalars associated with the incoming
geodesic null congruence $k^{\mu}$ can be derived from light front projections
of the Ricci identity
$2\nabla_{[\mu}\nabla_{\nu]}k^{\alpha}=-R_{\mu\nu\beta}{}^{\alpha}k^{\beta}$;
cf. Refs. (wal1984, , pp. 222–223) and joretal1961 . Evolving the energy-
rescaled variables of Eq. (17a), we obtain:

$\displaystyle E^{-1}k^{\rho}\partial_{\rho}E=-\,\alpha E\ ;$ (19a)
$\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\left(\frac{\hat{\theta}}{E}\right)=-\frac{1}{2}\left(\frac{\hat{\theta}}{E}-2\alpha\right)\left(\frac{\hat{\theta}}{E}\right)-2\left(\frac{\hat{\sigma}}{E}\right)^{2}-8\pi\varrho\
;$ (19b) $\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\left(\frac{\hat{\sigma}}{E}\right)^{2}=-2\left(\frac{\hat{\theta}}{E}-\alpha\right)\left(\frac{\hat{\sigma}}{E}\right)^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\ -\
2\left(\frac{\hat{\sigma}^{\mu\nu}}{E}\right)\left(p_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}E_{\alpha\beta}+\epsilon_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}H_{\alpha\beta}\right)\
,$ (19c)

where
$\hat{\sigma}^{2}:=\frac{1}{2}\hat{\sigma}_{\mu\nu}\hat{\sigma}^{\mu\nu}$ is
the (squared) optical shear scalar. In Eq. (19), the propagation of the
energy-rescaled optical shear scalar is sourced by the incoming radiative Weyl
curvature eigenfield defined by

$p_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}E_{\alpha\beta}+\epsilon_{\langle\mu}{}^{\alpha}p_{\nu\rangle}{}^{\beta}H_{\alpha\beta}\
,$ (20)

where
$E_{\alpha\beta}:=C_{\mu\nu\rho\sigma}u^{\mu}h_{\alpha}{}^{\nu}u^{\rho}h_{\beta}{}^{\sigma}$
and
$H_{\alpha\beta}:=-\frac{1}{2}\epsilon_{\rho\sigma\gamma\delta}C_{\mu\nu}{}^{\gamma\delta}u^{\rho}h_{\alpha}{}^{\sigma}u^{\mu}h_{\beta}{}^{\nu}$
are the electric and magnetic Weyl curvature tensors with respect to the
matter frame.

The propagation equations for the scalar-valued matter variables read:888With
Eq. (10a), when acting on scalar-valued variables, the operator identity
$E^{-1}k^{\rho}\partial_{\rho}=u^{\rho}\partial_{\rho}-e^{\rho}\partial_{\rho}$
applies.

$\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\theta=-\frac{1}{3}\,\theta^{2}-2\sigma^{2}-4\pi\varrho+\Lambda-\theta^{\prime}\
;$ (21a) $\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\varrho=-\theta\varrho-\varrho^{\prime}\ ;$ (21b)
$\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\sigma^{2}=-\frac{4}{3}\,\theta\sigma^{2}-\sigma_{\mu}{}^{\nu}\sigma_{\nu}{}^{\alpha}\sigma_{\alpha}{}^{\mu}-(\sigma^{\mu\nu}E_{\mu\nu})-(\sigma^{2})^{\prime};$
(21c) $\displaystyle
E^{-1}k^{\rho}\partial_{\rho}(p^{\mu\nu}\sigma_{\mu\nu})=-\frac{2}{3}\,\theta(p^{\mu\nu}\sigma_{\mu\nu})-p^{\mu\nu}\sigma_{\mu\alpha}\sigma_{\nu}{}^{\alpha}$
$\displaystyle\qquad\qquad\qquad\qquad\quad\
+\frac{4}{3}\,\sigma^{2}-(p^{\mu\nu}E_{\mu\nu})-2\sigma_{\mu\nu}e^{\mu}\kappa^{\nu}-(p^{\mu\nu}\sigma_{\mu\nu})^{\prime}\
,$ (21d)

where the operator ${}^{\prime}:=e^{\mu}\nabla_{\mu}$ denotes the spatial
derivative along $e^{\mu}$, and $(p^{\mu\nu}E_{\mu\nu})$ is the longitudinal
projection of the electric Weyl curvature tensor, which accounts for tidal
effects in the evolution equations.

We can recast Eqs. (21a) and (21) into the following evolution equations for
the area-adapted longitudinal expansion rate variables $\alpha$ and $\beta$:

$\displaystyle E^{-1}k^{\rho}\partial_{\rho}\alpha$ $\displaystyle=$
$\displaystyle-\alpha^{2}-\alpha^{\prime}+(p^{\mu\nu}E_{\mu\nu})-\frac{4\pi}{3}\varrho+\frac{1}{3}\Lambda-(e^{\mu}\sigma_{\mu\alpha})(e^{\nu}\sigma_{\nu\beta})p^{\alpha\beta}$
(22a) $\displaystyle+\ 2\sigma_{\mu\nu}e^{\mu}\kappa^{\nu}\ ;$ $\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\beta$ $\displaystyle=$
$\displaystyle-\frac{1}{9}\,(\alpha+2\beta)(4\beta-\alpha)-\beta^{\prime}-\frac{1}{2}\,(p^{\mu\nu}E_{\mu\nu})-\frac{4\pi}{3}\varrho+\frac{1}{3}\Lambda$
(22b) $\displaystyle-\
\frac{1}{2}\,(e^{\mu}\sigma_{\mu\alpha})(e^{\nu}\sigma_{\nu\beta})p^{\alpha\beta}-\frac{1}{2}\,\sigma_{\mu\alpha}\sigma_{\nu\beta}p^{\mu\nu}p^{\alpha\beta}-\sigma_{\mu\nu}e^{\mu}\kappa^{\nu}\
.$

We may, in addition, formulate the evolution equation for the cosmological
redshift, defined in Eq. (9). Using Eq. (19a), we have (cf. Ref. (ell1971, ,
p. 147f))

$k^{\rho}\partial_{\rho}z=-\,E(P)\alpha(1+z)^{2}\ .$ (23)

The cosmological redshift is a monotonic function along the geodesic null
curves only when the area-adapted longitudinal expansion rate variable
$\alpha$ does not change sign. Thus, as pointed out by Räsänen lightpropndust
, cosmological redshift can in general not be used as a parameter along the
geodesic null congruence when scales of collapsing structures are considered;
see section 4 for a discussion on valid parametrizations of light fronts. For
the same reason, cosmological redshift cannot in general be thought of as a
cosmological time variable, since, typically, $u^{\rho}\partial_{\rho}z$ will
change sign along the world lines of the matter congruence. The proper time
function $\tau$ for the matter congruence, on the other hand, is a monotonic
function along the geodesic null congruence,

$k^{\rho}\partial_{\rho}\tau=-k^{\mu}u_{\mu}=E\ ,$ (24)

since $E$ is a positive-valued variable. We remark that level surfaces
$\\{\lambda=\text{constant}\\}$ – with $\lambda$ defined as an affine
parameter along the geodesic null congruence satisfying
$k^{\rho}\partial_{\rho}\lambda=1$ – do not in general coincide with spacelike
$3$-surfaces $\\{\tau=\text{constant}\\}$. $E$ is in general an inhomogeneous
function on the screen space, which renders $\lambda$ an invalid label for the
light fronts $\\{V=\text{constant},\tau=\text{constant}\\}$.

###  2.5 Constraint equations on the light fronts

The null kinematic variables given in Eqs. (17a) and (17b) are subject to
constraints from the embedding of the $2$-dimensional light fronts into the
cosmological spacetime; cf. Refs. steetal2003 , (wal1984, , p. 258) and
(ste1991, , Sec. 16.2). When expressed in terms of the covariantly
$(1+1+2)$-decomposed variables, the Gauss embedding constraint reads:

$0=K-(p^{\mu\nu}E_{\mu\nu})-\frac{8\pi}{3}\varrho-\frac{1}{3}\Lambda+\frac{1}{4}\,\frac{\hat{\theta}}{E}\,\frac{\tilde{\theta}}{E}-\frac{1}{2}\,\frac{\hat{\sigma}_{\mu\nu}}{E}\,\frac{\tilde{\sigma}^{\mu\nu}}{E}\
,$ (25)

while the two Codazzi embedding constraints are given by

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
p^{\mu\alpha}p^{\nu\beta}\mbox{D}_{\nu}\left(\frac{\hat{\sigma}_{\alpha\beta}}{E}\right)-\frac{1}{2}\,p^{\mu\nu}\mbox{D}_{\nu}\left(\frac{\hat{\theta}}{E}\right)+\left(\frac{\hat{\sigma}^{\mu\nu}}{E}-\frac{1}{2}\,\frac{\hat{\theta}}{E}\,p^{\mu\nu}\right)\frac{\mbox{D}_{\nu}E}{E}\
;$ (26a) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
p^{\mu\alpha}p^{\nu\beta}\mbox{D}_{\nu}\left(\frac{\tilde{\sigma}_{\alpha\beta}}{E}\right)-\frac{1}{2}\,p^{\mu\nu}\mbox{D}_{\nu}\left(\frac{\tilde{\theta}}{E}\right)+\left(\frac{\tilde{\sigma}^{\mu\nu}}{E}-\frac{1}{2}\,\frac{\tilde{\theta}}{E}\,p^{\mu\nu}\right)\frac{\mbox{D}_{\nu}E}{E}\
,$ (26b)

respectively. In Eq. (25) the extrinsic curvature of the embedding has been
formulated using the area expansion tensors associated with $k^{\mu}$ and
$l^{\mu}$. The scalar $K:=\frac{1}{2}{}^{(2)}\\!R$ is the Gaussian curvature
of the screen space, and ${}^{(2)}\\!R$ is the Ricci scalar of the same screen
space. The sectional curvature of the light fronts is given by the projection

$\frac{1}{2}\,R_{\mu\nu\alpha\beta}\,p^{\mu\alpha}p^{\nu\beta}=(p^{\mu\nu}E_{\mu\nu})+\frac{8\pi}{3}\varrho+\frac{1}{3}\Lambda\
,$ (27)

where $R_{\mu\nu\alpha\beta}$ is the Riemann curvature tensor for the
cosmological spacetime.

###  2.6 Photon conservation, distance measures, and redshift drift

In the present analysis we consider light rays propagating along the
observer’s past null cone ${\cal C}^{-}(P):\\!\\{V=\text{constant}\\}$ to be
associated with freely streaming test particles which are – apart from the
point of emission and the point of observation at the vertex $P$ of the past
null cone – non-interacting with the matter content of the cosmological
spacetime. In this setting we are entitled to assume the conservation of the
current density,

$J^{\mu}={\mathfrak{n}}\,k^{\mu}\ ,$ (28)

where $\mathfrak{n}$ is the number density of light front continuum elements
(photons) with respect to the area enclosed by the screen space. The
conservation of $J^{\mu}$ leads to

$\nabla_{\mu}J^{\mu}=0\;\Leftrightarrow\;k^{\rho}\partial_{\rho}{\mathfrak{n}}=-\hat{\theta}\,{\mathfrak{n}}\
,$ (29)

i.e. the number of radiation elements per area element is conserved. This
conservation law is fundamental for the definition of cosmological distance
measures, see, e.g., Refs. ell1971  and sasaki1 , and is an important
assumption in order to obtain the relativistically corrected luminosity
distances.

The observer area distance $d_{A}$, henceforth referred to as the effective
angular diameter distance,999In generic cosmological spacetimes the observer
area distance for a luminous astrophysical source does not necessarily
coincide with the linear size based angular diameter distance for the same
source due to potential Weyl curvature induced shearing along the observer’s
past null cone of an incoming geodesic null ray bundle. Only in an exact FLRW
cosmology do their conceptions become identical as a consequence of the
prevailing spatial isotropy, and thus vanishing Weyl curvature; cf. Refs.
(ste1991, , Eq. (25,27)) and (ellhve1999, , Sec. 4.5.2). for a luminous
astrophysical source is defined by (ell1971, , Eq. (6.27))

$d_{A}^{2}:=\frac{{\delta}A}{{\delta}\Omega}\ ,$ (30)

where ${\delta}A$ is the physical area covered by the source at its own
position in directions perpendicular to the spatial propagation direction of
the emanating geodesic null rays, and ${\delta}\Omega$ is the solid angle the
source subtends at the position $P$ of the observer (with
$0<{\delta}\Omega<4\pi$). The evolution of the effective angular diameter
distance along the geodesic null congruence $k^{\mu}$ is given by (cf. Ref.
(ell1971, , Eq. (6.19)))

$E^{-1}k^{\rho}\partial_{\rho}d_{A}=\frac{1}{2}\,\frac{\hat{\theta}}{E}\,d_{A}\
;\qquad
E^{-1}k^{\rho}\partial_{\rho}d^{-2}_{A}=-\,\frac{\hat{\theta}}{E}\,d^{-2}_{A}\
.$ (31)

This, when combined with Eq. (19b), yields the focusing equation (cf. Ref.
(per2004, , Eq. (44))):

$E^{-1}k^{\rho}\partial_{\rho}(E^{-1}k^{\sigma}\partial_{\sigma}d_{A})=\left(\frac{1}{2}\,\alpha\,\frac{\hat{\theta}}{E}-\left(\frac{\hat{\sigma}}{E}\right)^{2}-4\pi\varrho\right)d_{A}\
.$ (32)

Using the conservation law (29) and Etherington’s reciprocity theorem
Etherington ; ellhve1999 ; Nez , the luminosity distance between a source and
the observer at the vertex point $P$,

$d_{L}:=\sqrt{\frac{L^{\mathrm{bol}}}{4\pi I^{\mathrm{bol}}}}\ ,$ (33)

can be written in terms of the effective angular diameter distance as follows:

$d_{L}=(1+z)^{2}d_{A}\ ;$ (34)

$L^{\mathrm{bol}}$ is the bolometric luminosity of the source, and
$I^{\mathrm{bol}}$ is the bolometric flux of energy as measured by the
observer.

Let us finally consider the drift of cosmological redshift of a given source
in time as measured by an observer at the vertex point $P$. Let the source and
the observer be comoving with the irrotational $4$-velocity field $u^{\mu}$ of
the matter congruence and let the incoming geodesic null rays connecting them
be generated by $k^{\mu}$. The spacetime point of observation $P$ and the
spacetime point of the source $S$ are shifted along the matter congruence at
their respective locations, creating a change in cosmological redshift which
we measure in the observer’s proper time $\tau_{0}$ Heinesen:2020pms :

$\displaystyle\frac{dz}{d\tau_{0}}$ $\displaystyle=$ $\displaystyle
u^{\rho}\partial_{\rho}z\,\lvert_{P}+\frac{d\tau}{d\tau_{0}}\,u^{\rho}\partial_{\rho}z\,\lvert_{S}\,=\,(1+z)\,\alpha\,\rvert_{P}-(1+z)\,\frac{E^{\prime}}{E}\,\bigg{\rvert}_{P}-\alpha\,\rvert_{S}+\frac{E^{\prime}}{E}\,\bigg{\rvert}_{S},$
(35)

where $d\tau/d\tau_{0}=1/(1+z)$ is the Jacobian of the change of the proper
time measure between the source and observer world lines from the bijection
induced by the geodesic null congruence. The second equality can be obtained
by substituting $u^{\mu}=k^{\mu}/E+e^{\mu}$ and using Eq. (19a). We see that
the cosmological redshift drift (35) depends on the area-adapted longitudinal
expansion rate variable, $\alpha$, of the matter congruence (replacing the
Hubble parameter in the analogous FLRW expression for redshift drift) and on
the spatial gradient of the photon energy, $E$, along the spatial direction of
the incoming null ray.

##  3 Area-averaging over light fronts

We now consider the evolution of area-averaged physical quantities along the
observer’s past null cone. In section 3.1 we define the area-averaging
operation on the light fronts and construct evolution equations for the area
of the light fronts. In section 3.2 we consider conservation laws and the
Gauss–Bonnet theorem as constraints for the average dynamics of the light
fronts. In section 3.3 we derive evolution equations for area-averaged
cosmological redshift and distance measures, and in section 3.4 we check the
relations derived against the FLRW limit.

###  3.1 Area-averaging operation and evolution equations for area-averaged
variables

We shall now define the covariant area-averaging operation on the light
fronts. The area-averaging operation defined is invariant under coordinate
transformations (though we shall define useful area-adapted coordinates on the
$2$-surfaces of integration), and it can be formulated in terms of covariant
integrals in the cosmological spacetime, as detailed in Ref. lightconeav1 .

Consider an orientable and simply connected compact $2$-dimensional domain
${\cal F}$ in the light front which lies within a $2$-surface of intersection
between the observer’s past null cone ${\cal C}^{-}(P)$ and a spacelike
$3$-surface ${\cal S}(\tau)$, the latter labelled by the observer’s
corresponding value of proper time: ${{\cal F}}(\tau)\subseteq({\cal
C}^{-}(P)\cap{\cal S})(\tau)$. The total domain $({\cal C}^{-}(P)\cap{\cal
S})(\tau)$ is assumed to be topologically closed (as, e.g., the all-sky last
scattering surface of the Cosmic Microwave Background radiation). For the
general spatially inhomogeneous case, this assumption corresponds to a choice
of closed space form, e.g., a spherical topology in the simplest case. We
shall allow ${{\cal F}}(\tau)$ to be any compact $2$-dimensional subdomain of
$({\cal C}^{-}(P)\cap{\cal S})(\tau)$ with boundaries propagated along
$k^{\mu}$ such that we track a system of constant number of radiation
elements.

For the present purpose, it is convenient to use local spacetime coordinates
$x^{\mu}=(\tau,V,x^{A})$, where $x^{A}$, with indices $A=2,3$, are local
coordinates adapted to the screen spaces such that they are constant along the
geodesic null rays:101010In general, local coordinates preserved along
$k^{\mu}$ are not preserved along $u^{\mu}$, since
$k^{\rho}\partial_{\rho}x^{A}=0$ implies that
$k^{\alpha}\nabla_{\alpha}(u^{\rho}\partial_{\rho}x^{A})=\pounds_{\textbf{k}}(u^{\rho})\partial_{\rho}x^{A}$,
where the operator $\pounds_{\textbf{k}}$ is the Lie derivative along the
geodesic null rays associated with $k^{\mu}$. Thus, if the change of $u^{\mu}$
along the geodesic null rays has components tangential to the screen space
(which happens generically in spatially inhomogeneous cosmologies), then
$u^{\rho}\partial_{\rho}x^{A}=0$ cannot be satisfied globally on the past null
cone. $k^{\rho}\partial_{\rho}x^{A}=0$. Gasperini, Veneziano and their
collaborators introduced these coordinates as ‘Geodesic Light cone Coordinates
(GLC)’ with $\tau$ and $V$ exchanged; see Refs. lightconeav1 , Fleuryetal:GLC
and appendix A for details, and [chap.3 and references therein]elmardi:phd
for the relation to George Ellis’ observational ‘Perturbed Light cone Gauge’
(PLG) coordinates. In GLC coordinates the components of the cosmological
spacetime metric and its inverse read:

$g_{\mu\nu}=\left(\begin{array}[]{ccc}0&1/E&\vec{0}\\\
1/E&\,(1/E^{2})\\!+\\!U^{2}&U_{B}\\\ \vec{0}&U_{A}&p_{AB}\\\
\end{array}\right)\
;~{}~{}~{}~{}~{}~{}~{}g^{\mu\nu}=\left(\begin{array}[]{ccc}-1&\,E&-EU^{B}\\\
E&0&\vec{0}\\\ -EU^{A}&\vec{0}&p^{AB}\\\ \end{array}\right)\ ,$ (36)

and the screen space-forming vectors and associated one-forms read (see Fig.
1):

$\displaystyle k_{\mu}=(0,1,0,0)\ ;$ $\displaystyle k^{\mu}=(E,0,0,0)\ ;$
$\displaystyle l_{\mu}=(-2E,-1,0,0)\ ;$ $\displaystyle
l^{\mu}=(E,-2E^{2},2E^{2}U^{A})\ ;$ $\displaystyle u_{\mu}=(-1,0,0,0)\quad;$
$\displaystyle u^{\mu}=(1,-E,EU^{A})\ ;$ $\displaystyle e_{\mu}=(-1,-1/E,0,0)\
;$ $\displaystyle e^{\mu}=(0,-E,EU^{A})\ ,$ (37)

with the coordinate-independent projections
$k_{\mu}k^{\mu}=l_{\mu}l^{\mu}=u^{\mu}e_{\mu}=u_{\mu}e^{\mu}=0$,
$u_{\mu}u^{\mu}=-1$, $e_{\mu}e^{\mu}=1$,
$k_{\mu}u^{\mu}=k^{\mu}u_{\mu}=l_{\mu}u^{\mu}=l^{\mu}u_{\mu}=k_{\mu}e^{\mu}=k^{\mu}e_{\mu}=-E$,
$l_{\mu}e^{\mu}=l^{\mu}e_{\mu}=E$, $k_{\mu}l^{\mu}=k^{\mu}l_{\mu}=-2E^{2}$.
The $2$-dimensional tensor $p_{AB}$ is the induced metric on the screen space
and $p^{AB}$ is its inverse, while the $2$-dimensional vector
$U^{A}=u^{\rho}\partial_{\rho}(x^{A})/E$ defines the drift of the area-adapted
coordinates $x^{A}$ along the matter congruence, and
$U^{2}:=p_{AB}U^{A}U^{B}$.

Figure 1: The $4$-vector fields $\boldsymbol{k}$, $\boldsymbol{l}$,
$\boldsymbol{u}$ and $\boldsymbol{e}$ (right column of (3.1)) at a line of
intersection $\\{V=\text{constant},\tau=\text{constant}\\}$ of the past null
cone ${\cal C}^{-}(P)$, embedded into the cosmological spacetime with GLC
coordinates $\\{\tau,V,x^{A=2}\\}$, and a single spacelike $3$-surface ${\cal
S}(\tau)$. The third dimension of the family of spacelike $3$-surfaces
$\\{\tau=\text{constant}\\}$, as parametrized by $x^{A=3}$, is suppressed. The
wave vector $\boldsymbol{k}$ is tangent to the past null cone, while the wave
vector $\boldsymbol{l}$ represents an auxiliary null-like field pointing in
the opposite spatial direction of $\boldsymbol{k}$ if seen in the restframe of
the $4$-velocity $\boldsymbol{u}$, which is orthogonal to the spacelike
$3$-surfaces (containing the vector $\boldsymbol{e}$ as a tangent).
$\boldsymbol{u}$ is also, together with the wave vector $\boldsymbol{k}$,
orthogonal to the $2$-dimensional light fronts
$\\{V=\text{constant},\tau=\text{constant}\\}$ (blue dashed line). Recall the
vector relations (10a) and (10b), and see related figures in Ref.
Fleuryetal:GLC . (Past null cone in this figure inspired by Fig. 12 of Ref.
per2004 ; spacetime foliation in this figure: Pierre Mourier (priv. comm.).)

We define the area of a $2$-surface within a screen space by

$A_{{\cal F}}(\tau)=|{\cal F}(\tau)|=\int_{{\cal F}}\sqrt{p(\tau,x^{A})}\,{\rm
d}^{2}x\ ,$ (38)

where $p(\tau,x^{A}):=\det(p_{AB})$ is the Riemannian area element of the
$2$-surface. We define the area-average of a general scalar-valued function
$F(\tau,x^{A})$ as follows:

$\langle F\rangle_{{\cal F}}\,(\tau):=\frac{1}{A_{{\cal F}}(\tau)}\int_{{\cal
F}}F(\tau,x^{A})\sqrt{p(\tau,x^{A})}\,{\rm d}^{2}x\ .$ (39)

It shall also be convenient to define area-averaging of products of scalar-
valued functions $F$ and $G$ as evaluated at different points along the
geodesic null rays (parametrized in terms of the observer’s proper time
function $\tau$):

$\displaystyle\langle F(\tau,x^{A})G(\tau^{\prime},x^{A})\rangle_{{\cal F}}:=$
$\displaystyle\frac{1}{A_{{\cal F}}(\tau)}\int_{{\cal
F}}F(\tau,x^{A})G(\tau^{\prime},x^{A})\sqrt{p(\tau,x^{A})}\,{\rm d}^{2}x\ ,$
(40)

where the evaluation of the area element is at the unprimed coordinate time
$\tau$.

The commutation rule for the operations of area-averaging and evolution along
the observer’s past null cone is

$\displaystyle\partial_{\tau}\langle F\rangle_{{\cal
F}}-\left\langle\partial_{\tau}F\right\rangle_{{\cal F}}$ $\displaystyle=$
$\displaystyle\text{Cov}_{{\cal F}}\left(\frac{\hat{\theta}}{E}\,,\,F\right)\
,$ (41)

where $\partial_{\tau}:=E^{-1}k^{\rho}\partial_{\rho}$ is the directional
derivative along the null congruence as measured in units of $\tau$ (cf. Eq.
(24) for the Jacobian, $k^{\rho}\partial_{\rho}\tau=E$), and where, for
scalar-valued functions $X$ and $Y$:

$\displaystyle\text{Cov}_{{\cal F}}(X,Y):=\langle XY\rangle_{\cal F}-\langle
X\rangle_{\cal F}\langle Y\rangle_{\cal F}\ .$ (42)

The local evolution of a scalar-valued function per unit proper time is
naturally measured along the geodesic null rays of $k^{\mu}$, since the
boundaries of the spatial domain as well as the past null cone scalar $V$ are
comoving with $k^{\mu}$. In deriving Eq. (41), the Jacobi identity,

$E^{-1}k^{\rho}\partial_{\rho}\sqrt{p(\tau,x^{A})}=\partial_{\tau}\sqrt{p(\tau,x^{A})}=\frac{\hat{\theta}}{E}\,\sqrt{p(\tau,x^{A})}\
,$ (43)

has been used; cf. Refs. (ell1971, , Eq. (6.19)) and (sasaki1, , Eq. (3.14)).
We define the domain-dependent dimensionless area scale factor for the light
fronts:

$a_{{\cal F}}(\tau):=\left(\frac{A_{{\cal F}}(\tau)}{A_{{{\cal
F}}_{i}}}\right)^{1/2}\ ,$ (44)

where the reference area $A_{{{\cal F}}_{i}}:=A_{{\cal F}}(\tau_{i})$ is
evaluated at a reference light front domain with time label $\tau_{i}$. Since
the domain ${{\cal F}}(\tau)$ is defined as having no flow of geodesic null
rays across its boundary, it fails in general to preserve the number of matter
fluid elements. Thus, $a_{\cal F}(\tau)$ cannot be interpreted as an average
scaling of the coordinate distance between individual world lines of the
matter congruence, and so does not reduce to the FLRW scale factor in the case
of spatial homogeneity and isotropy on the spacelike $3$-surfaces orthogonal
to the matter congruence. We formulate the evolution along the observer’s past
null cone ${\cal C}^{-}(P):\\!\\{V=\text{constant}\\}$ of the area scale
factor as a first part of a theorem.

######  Theorem 4.a (Effective evolution equations for light fronts)

The evolution along the observer’s past null cone ${\cal
C}^{-}(P):\\!\\{V=\text{constant}\\}$ of the area scale factor is governed by
the screen space area expansion rate:

$\mathcal{H}_{{\cal F}}:=\frac{\partial_{\tau}a_{{\cal F}}(\tau)}{a_{{\cal
F}}(\tau)}=\frac{1}{2}\left\langle\frac{\hat{\theta}}{E}\right\rangle_{{\cal
F}}\ ,$ (45)

where we defined the $2$-surface area expansion functional $\mathcal{H}_{{\cal
F}}$. The second derivative of the area scale factor yields the ‘screen space
area acceleration law’:

$\frac{\partial^{2}_{\tau}a_{{\cal F}}(\tau)}{a_{{\cal
F}}(\tau)}=-4\pi\left\langle\varrho\right\rangle_{{\cal
F}}+\mathcal{Q}^{\mathbf{k}}_{{\cal F}}+{\mathcal{S}}^{\mathbf{k}}_{{\cal F}}\
,$ (46)

where we have used the optical evolution equations (19a) and (19b). The
function

$\displaystyle\mathcal{Q}^{\mathbf{k}}_{{\cal F}}$ $\displaystyle:=$
$\displaystyle\frac{1}{4}\,\text{\rm Cov}_{{\cal
F}}\left(\frac{\hat{\theta}}{E}\,,\,\frac{\hat{\theta}}{E}\right)-\left\langle\frac{\hat{\sigma}^{2}}{E^{2}}\right\rangle_{{\cal
F}}$ (47)

is a ‘screen space kinematic backreaction’ term arising from local spatial
inhomogeneity and anisotropy of the incoming geodesic null ray bundle
associated with $k^{\mu}$. The function

${\mathcal{S}}^{\mathbf{k}}_{\cal
F}:=\frac{1}{2}\left\langle\frac{\hat{\theta}}{E}\,\alpha\right\rangle_{{\cal
F}}=\frac{1}{2}\left\langle\frac{\hat{\theta}}{E}\,k^{\rho}\partial_{\rho}\left(\frac{1}{E}\right)\right\rangle_{{\cal
F}}=-\frac{1}{2}\left\langle\frac{\hat{\theta}}{E}\,\frac{\partial_{\tau}E}{E}\right\rangle_{{\cal
F}}$ (48)

arises from the re-parametrization of the incoming geodesic null congruence
and measures the failure of $\tau$ to be an affine parameter along this null
congruence.

Remarks to Theorem 4.a

Removing the term (48) is possible by averaging over
$\\{\lambda=\text{constant}\\}$ light fronts instead of
$\\{\tau=\text{constant}\\}$ light fronts. We stick to the
$\\{\tau=\text{constant}\\}$ foliation of the observer’s past null cone in the
present paper. The term ${\mathcal{S}}^{\mathbf{k}}_{\cal F}$ vanishes for
energy functions $E=k^{\rho}\partial_{\rho}\tau$, cf. Eq. (24), that are
constant along the integral curves of the geodesic null congruence, i.e., when
the geodesic null rays are not subject to cosmological redshift. The function
$1/E$ plays the role of an inhomogeneous lapse function along the geodesic
null congruence, cf. Eq. (36) and appendix A.111111For a discussion of the
degeneracy at the null cone of the standard ADM (Arnowitt, Deser and Misner)
$3+1$ slicing formalism ADM  and its generalization for lightlike foliations,
exemplified for double null foliations, see Refs. Israel:nullADM  and
Israel:doublenull .

We formulate the area-average of the Gauss embedding constraint (25) as a
second part of this theorem.

######  Theorem 4.b (Effective energy constraint on light fronts)

Averaging the Gauss embedding constraint (25) over the compact domain ${\cal
F}$ of the light fronts yields the ‘screen space area expansion law’:

$\displaystyle\mathcal{H}_{{\cal F}}^{2}+\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal F}}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{3}\left\langle\varrho\right\rangle_{{\cal
F}}+\frac{1}{3}\Lambda-\left\langle K\right\rangle_{{\cal
F}}-\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{{\cal F}}\ ,$ (49)

where

$\displaystyle\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{{\cal
F}}:=\frac{1}{4}\,\text{\rm Cov}_{{\cal
F}}\left(\frac{\hat{\theta}}{E}\,,\,\frac{\tilde{\theta}}{E}\right)-\frac{1}{2}\left\langle\frac{\hat{\sigma}_{\mu\nu}}{E}\,\frac{\tilde{\sigma}^{\mu\nu}}{E}\right\rangle_{{\cal
F}}-\left\langle p^{\mu\nu}E_{\mu\nu}\right\rangle_{{\cal F}}$ (50)

is a second backreaction term on the screen space, which may be substituted by
the linearly transformed backreaction variable, using the relations (18):

$\displaystyle\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{\cal
F}:=\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{{\cal
F}}-\mathcal{Q}^{\mathbf{k}}_{\cal F}=\frac{1}{2}\,\text{\rm Cov}_{{\cal
F}}\left(\frac{\hat{\theta}}{E}\,,\,\overline{\theta}\right)-2\left\langle\frac{\hat{\sigma}_{\mu\nu}}{E}\
\overline{\sigma}^{\mu\nu}\right\rangle_{{\cal F}}-\left\langle
p^{\mu\nu}E_{\mu\nu}\right\rangle_{{\cal F}}\ .$ (51)

Remarks to Theorem 4.b

Notice that the new backreaction term (50) contains products of the kinematic
variables associated with both the incoming and outgoing null congruences,
$k^{\mu}$ and $l^{\mu}$, contrary to the backreaction term (47).121212We
remark that this backreaction term appears to bear a relationship to the
Hawking–Hayward energies as defined in Refs. Hawking ; Hayward ; see Refs.
bengtsson:hawkingenergy ; Dennis1 ; Dennis2  for recent work. Moreover, in
both backreaction terms, (50) and (51), tidal effects enter through the screen
space projection of the electric Weyl curvature tensor. This is contrary to
the $3$-dimensional volume-averaging operation of Refs. buchert:grgdust ;
buchert:grgfluid , where Weyl curvature only enters implicitly in the large-
scale volume evolution.

Since Eq. (49) must be the integral of Eq. (46), we take the derivative of Eq.
(49) with respect to $\tau$ and re-insert Eqs. (49) and (46) to obtain the
following integrability condition on the light fronts that we formulate as a
third part of theorem 44.

######  Theorem 4.c (Integrability condition)

A necessary condition of integrability for Eq. (46) to yield Eq. (49) is given
by the integral constraint:

$\displaystyle\partial_{\tau}\mathcal{Q}^{\mathbf{k}}_{{\cal
F}}+4\mathcal{H}_{{\cal F}}\mathcal{Q}^{\mathbf{k}}_{{\cal
F}}+\partial_{\tau}\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{\cal
F}+2\mathcal{H}_{{\cal F}}\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{{\cal
F}}+\partial_{\tau}\left\langle K\right\rangle_{{\cal F}}+2\mathcal{H}_{{\cal
F}}\left\langle K\right\rangle_{{\cal F}}$ (52a) $\displaystyle+$
$\displaystyle\partial_{\tau}\left(\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal
F}}\right)+2\mathcal{H}_{{\cal F}}\left(\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal
F}}\right)+2\mathcal{H}_{{\cal F}}\left({\mathcal{S}}^{\mathbf{k}}_{\cal
F}-\frac{1}{3}\Lambda\right)$ $\displaystyle=$
$\displaystyle\frac{8\pi}{3}\left(\partial_{\tau}\left\langle\varrho\right\rangle_{{\cal
F}}+5\mathcal{H}_{{\cal F}}\left\langle\varrho\right\rangle_{{\cal F}}\right)\
,$ where it follows from averaging of the local law for rest mass flow Eq.
(21b) that
$\displaystyle\partial_{\tau}\left\langle\varrho\right\rangle_{{\cal
F}}+\left\langle\theta\right\rangle_{{\cal
F}}\left\langle\varrho\right\rangle_{{\cal
F}}=-\langle\varrho^{\prime}\rangle_{{\cal F}}-\text{\rm Cov}_{{\cal
F}}\left(\theta\,,\,\varrho\right)+\text{\rm Cov}_{{\cal
F}}\left(\frac{\hat{\theta}}{E}\,,\,\varrho\right)\ .$ (52b) The integrability
condition in Eq. (52a) can be written in compact form:
$\displaystyle\frac{1}{a_{\cal
F}^{4}}\partial_{\tau}\left[\mathcal{Q}^{\mathbf{k}}_{\cal F}a_{\cal
F}^{4}\right]+\frac{1}{a_{\cal
F}^{2}}\partial_{\tau}\left[\left(\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{{\cal
F}}+{\mathcal{S}}^{\mathbf{k}}_{\cal F}-\frac{1}{3}\Lambda+\left\langle
K\right\rangle_{{\cal F}}+\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal F}}\right)a_{\cal
F}^{2}\right]$ (52c) $\displaystyle=$
$\displaystyle\frac{8\pi}{3}\frac{1}{a_{\cal
F}^{5}}\partial_{\tau}\left[\left\langle\varrho\right\rangle_{{\cal F}}a_{\cal
F}^{5}\right]+\partial_{\tau}{\mathcal{S}}^{\mathbf{k}}_{\cal F}\ .$

Remarks to Theorem 4.c

Looking at Eq. (52c), we appreciate that separate conservation of screen space
kinematic backreaction implies the following scaling behaviour
$\mathcal{Q}^{\mathbf{k}}_{\cal F}\propto a_{\cal F}^{-4}$. Furthermore,
separate conservation of the area-averaged Gaussian curvature implies
$\left\langle K\right\rangle_{{\cal F}}\propto a_{\cal F}^{-2}$, being
inversely proportional to the square of the area scale factor; see section 3.2
below for a discussion on the Gauss–Bonnet theorem and area-averaged Gaussian
curvature. (Compare also Ref. GBC  for a discussion of the averaged equations
in a $(2+1)$-spacetime.) The scaling behaviours derived from requiring
separate conservation laws of the cosmological variables in Eq. (52c) can in
some cases provide some guidance, but will in some cases also be misleading.
In general we expect non-trivial couplings between the various macroscopic
variables.

Discussion of Theorem 44

The area-averaged equations (45), (46), (49) and (52) appear more involved
than the corresponding evolution equations for the volume of fluid-orthogonal
spacelike $3$-surfaces in an irrotational dust universe given in Ref.
buchert:grgdust ; cf. Ref. buchert:grgfluid  for the corresponding formulation
including non-zero pressure in the fluid description. This is not surprising,
since the consideration of generic cosmological dynamics and corresponding
observations along the observer’s past null cone introduces additional
variables. In particular, in a covariant $(1+1+2)$-decomposition we see the
effects induced by the Weyl curvature. Even though the incoming geodesic null
rays in the present description are non-gravitating test particles – and thus
do not collectively act as a source in Einstein’s field equations – they enter
in the description of the observer’s past null cone and in the evolution along
this past null cone of quantities to be measured.

At this stage, the screen space area acceleration law (46) of theorem 44 for
the area scale factor $a_{\cal F}(\tau)$ in 4.a is subject to the area-
averaged constraint equation (49) in 4.b. These equations involve three
backreaction terms, $\mathcal{Q}^{\mathbf{k}}_{{\cal F}}$,
$\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{\cal F}$ and
$\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{\cal F}$, where the latter two are
trivially related to the first, $\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{\cal
F}-\mathcal{Q}^{\mathbf{k}\mathbf{e}}_{\cal F}=\mathcal{Q}^{\mathbf{k}}_{{\cal
F}}$. The integrability condition (52) in 4.c serves as a balance equation
among the various sources and backreaction terms, and is not an independent
equation. In addition, these equations involve the further variables
$\left\langle\overline{\theta}\right\rangle_{{\cal F}}$, the area-averaged
rest mass density $\left\langle\varrho\right\rangle_{\cal F}$, the area-
averaged Gaussian curvature scalar $\left\langle K\right\rangle_{\cal F}$, and
the non-affine term ${\mathcal{S}}^{\mathbf{k}}_{\cal F}$ (48). This forms a
set of two independent equations for seven unknown functions, in contrast to
the $3$-dimensional dust case in Ref. buchert:grgdust , where the counting
yields three independent equations (those that correspond to the two above,
but also the continuity equation for the averaged fluid density) for four
unknown functions. However, the present case is to be formally131313We still
refer to the dust case here; a generalization including fluid pressure with
general considerations on lapse and shift can be found in fanizza:pressure .
compared with the $3$-dimensional case in Ref. buchert:grgfluid , where the
counting also yields two independent equations for seven unknown functions
(that can be reduced by one via an equation of state). Unlike the dust case
buchert:grgdust , but similar to the more involved case buchert:grgfluid , the
evolution equations for
$\partial_{\tau}\left\langle\overline{\theta}\right\rangle_{{\cal F}}$ and
$\partial_{\tau}\left\langle\varrho\right\rangle_{{\cal F}}$, which can be
derived from Eqs. (21a) and (21b), i.e.
$\partial_{\tau}\left(\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal
F}}\right)+2\mathcal{H}_{{\cal F}}\left(\mathcal{H}_{{\cal
F}}\left\langle\overline{\theta}\right\rangle_{{\cal F}}\right)=...$,
$\partial_{\tau}\left\langle\varrho\right\rangle_{{\cal
F}}+5\mathcal{H}_{{\cal F}}\left\langle\varrho\right\rangle_{{\cal F}}=...$,
do not immediately serve to constrain the system further due to failure of the
light front domain to be comoving with the matter congruence, and due to the
appearance of further backreaction terms.

As detailed in the next section, invoking the Gauss–Bonnet theorem has the
potential to add a further constraint to the non-closed system of area-
averaged equations.

###  3.2 Average conservation laws and integral-geometric measures

Conservation laws can play an important role in the closure of averaged
evolution equations. In this section we investigate the possibility of
defining globally conserved quantities over the light fronts. The local photon
number conservation law (29) leads to the area-averaged equivalent:

$\partial_{\tau}\langle{\mathfrak{n}}\rangle_{{\cal F}}+2\mathcal{H}_{{\cal
F}}\langle{\mathfrak{n}}\rangle_{{\cal F}}=0\ \Leftrightarrow\
\langle{\mathfrak{n}}\rangle_{{\cal
F}}=\frac{\left\langle{\mathfrak{n}}\right\rangle_{{{\cal
F}}_{i}}}{a^{2}_{{\cal F}}(\tau)}\ .$ (53)

In the case when we consider a global area-average over a total light front
domain, a $2$-surface $\Sigma:=({\cal C}^{-}(P)\cap{\cal S})(\tau)$ which we
assume to be compact and topologically closed, we may invoke the Gauss–Bonnet
theorem for the area-averaged Gaussian curvature scalar $K$:

$\left\langle K\right\rangle_{\Sigma}=\frac{\left\langle
K\right\rangle_{{\Sigma}_{i}}}{a^{2}_{\Sigma}(\tau)}\ ,\qquad\left\langle
K\right\rangle_{{\Sigma}_{i}}A_{{\Sigma}_{i}}=2\pi\chi\ ,$ (54a) where
$\chi=2-2g$ is the Euler characteristic of $\Sigma$ with genus $g$. Thus, in
this case the area-averaged Gaussian curvature scalar obeys its own separate
conservation equation and decouples from the integrability condition (52a).
For light fronts $\Sigma$ of spherical topology we have $\chi=2$, and so
$\left\langle K\right\rangle_{{\Sigma}_{i}}A_{{\Sigma}_{i}}=4\pi$. (Level
$2$-surfaces of a past null cone must always have spherical topology, unless
caustics change the topology by destroying the bijection between the surfaces
ellissolomons ; mustapha3 ; Dennis1 .)

In the case of a light front with boundary, we may invoke the Gauss–Bonnet
theorem with boundary term, now considering a compact domain ${\cal
F}\subset\Sigma$:

$\left\langle K\right\rangle_{{{\cal F}}}A_{{\cal F}}\,+\int_{\partial{\cal
F}}\kappa\,{\mathrm{d}}{\ell}=2\pi\chi\ ,$ (54b)

where $\kappa$ is the extrinsic curvature of the $1$-dimensional boundary
$\partial{\cal F}$ of ${\cal F}$, and ${\mathrm{d}}{\ell}$ is the line element
for this boundary.

We notice from the above considerations the natural appearance of integral-
geometric measures of the light fronts, which can be viewed as generalizations
of the Minkowski functionals of convex sets in Euclidian space; cf. Ref.
beyond  and references therein: the $2$-surface area, $A$ (content), the
length of the circumference, $L:=\int_{\partial{\cal F}}{\mathrm{d}}{\ell}$
(shape), and the Euler characteristic, $\chi$ (connectivity). As the area of a
light front evolves, the shape of the light front also evolves, while its
Euler characteristic is preserved according to our assumption of caustic-free
evolution along the observer’s past null cone, irrespective of the choice of
compact submanifold of the total light front domain $({\cal C}^{-}(P)\cap{\cal
S})(\tau)$.

The evolution of the area-averaged Gaussian curvature scalar $\left\langle
K\right\rangle_{{\cal F}}$ is in general determined by the evolution of the
area of the $2$-surface $A_{{\cal F}}(\tau)$ and the boundary term
$\int_{\partial{\cal F}}\kappa{\mathrm{d}}{\ell}$. Whereas the evolution of
the area is given by the area-averaged expansion rate scalar of the geodesic
null congruence, implicitly given by Eq. (45), the evolution of the boundary
term involves shear degrees of freedom as well as Weyl curvature degrees of
freedom.

Let the boundary at each $2$-dimensional screen space be determined by the
level surfaces of value $B_{0}$ of a spacetime function $B$, such that

$n^{\mu}:=\frac{p^{\mu\nu}\partial_{\nu}B}{\mathcal{N}}\
;\qquad\mathcal{N}:=\sqrt{p^{\mu\nu}\partial_{\mu}B\,\partial_{\nu}B}\ ;\qquad
k^{\rho}\partial_{\rho}B=0\ ,$ (55)

defines the outward directed normal to the boundary.141414For an analogous
definition of domain boundaries in the $3+1$ slicing formalism, see Ref.
Gasperini:2009mu . The last condition ensures that the boundaries are constant
along the integral curves of the geodesic null congruence generating the
observer’s past null cone. The normalized tangent vector of the curve in the
light front plane, $q^{\mu}$, is determined by:

$q^{\mu}n_{\mu}=0\ ;\qquad q^{\mu}u_{\mu}=0\ ;\qquad q^{\mu}e_{\mu}=0\ ,$ (56)

and the induced metric on the boundary is given by

$q_{\mu}q^{\nu}=p_{\mu}{}^{\nu}-n_{\mu}n^{\nu}=\delta_{\mu}{}^{\nu}+u_{\mu}u^{\nu}-e_{\mu}e^{\nu}-n_{\mu}n^{\nu}\
,$ (57)

such that $\\{u^{\mu},e^{\mu},n^{\mu},q^{\mu}\\}$ constitute an orthonormal
vector basis for the cosmological spacetime. The extrinsic curvature
associated with the embedding of the boundary in the $2$-dimensional screen
space can now be expressed as follows:

$\kappa:=n^{\nu}q^{\mu}\nabla_{\mu}q_{\nu}\ .$ (58)

From the normalization requirement, $n_{\mu}k^{\nu}\nabla_{\nu}n^{\mu}=0$, and
the propagation rule, $k^{\rho}\partial_{\rho}B=0$, we have

$\displaystyle\frac{k^{\rho}\partial_{\rho}\mathcal{N}}{\mathcal{N}}=-\frac{1}{2}\,\hat{\theta}-n^{\mu}n^{\nu}\hat{\sigma}_{\mu\nu}\
.$ (59)

The evolution of the normal $n^{\mu}$ from one screen space to the next is
given by

$\displaystyle k^{\nu}\nabla_{\nu}n^{\mu}$ $\displaystyle=$
$\displaystyle-q^{\mu}\,q^{\alpha}n^{\nu}\hat{\sigma}_{\alpha\nu}-k^{\mu}\,n^{\alpha}e^{\nu}\sigma_{\alpha\nu}\
,$ (60)

which follows from the definition of $n^{\mu}$ in Eq. (55). Similarly, we may
compute the evolution of $q^{\mu}$ along the geodesic null congruence:

$\displaystyle k^{\nu}\nabla_{\nu}q^{\mu}$ $\displaystyle=$ $\displaystyle
n^{\mu}\,q^{\alpha}n^{\nu}\hat{\sigma}_{\alpha\nu}-k^{\mu}\,q^{\alpha}e^{\nu}\sigma_{\alpha\nu}\
,$ (61)

which follows from the orthogonality requirements (56). The evolution of the
boundary term is then given by

$\partial_{\tau}\int_{\partial{\cal
F}}\kappa\,{\mathrm{d}}{\ell}=\int_{\partial{\cal
F}}\frac{\nabla_{\mu}\left(\kappa\,\mathcal{N}\,k^{\mu}\right)}{\mathcal{N}E}\,{\mathrm{d}}{\ell}=\int_{\partial{\cal
F}}\left(E^{-1}k^{\rho}\partial_{\rho}\kappa+\frac{\kappa\,q^{\mu}q^{\nu}\nabla_{\mu}k_{\nu}}{E}\right){\mathrm{d}}{\ell}\
,$ (62)

where Eqs. (56) and (59) have been used, and where

$q^{\mu}q^{\nu}\nabla_{\mu}k_{\nu}=\frac{1}{2}\hat{\theta}+q^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu}$
(63)

is the projected expansion tensor along the tangent of the curve, describing
the expansion of the length of the boundary.

An important observation is that the intermediate equality in Eq. (62) shows
that the boundary term is conserved, if $(\kappa\,\mathcal{N}k^{\mu})$ is a
conserved current in the vicinity of the boundary, corresponding to $\kappa$
being a conserved density on the boundary domain. A further insight is
obtained when considering the evolution of the extrinsic curvature of the
boundary:

$\displaystyle k^{\rho}\partial_{\rho}\kappa$ $\displaystyle=$
$\displaystyle-\left(\frac{1}{2}\,\hat{\theta}+q^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu}\right)\left(\kappa-n^{\alpha}e^{\beta}\sigma_{\alpha\beta}\right)-n^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu}\left(q^{\alpha}e^{\beta}\sigma_{\alpha\beta}-\frac{q^{\rho}\partial_{\rho}E}{E}\right)$
(64) $\displaystyle+$ $\displaystyle
q^{\alpha}\nabla_{\alpha}(n^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu})-E\left(E_{\mu\nu}e^{\mu}n^{\nu}+\epsilon_{\mu\nu}{}^{\beta}H_{\alpha\beta}n^{\mu}q^{\nu}q^{\alpha}\right)\
.$

In deriving Eq. (64), the evolution equations (60) and (61) have been invoked
and it has been used that

$\displaystyle
n^{\mu}k^{\nu}q^{\alpha}q^{\beta}R_{\mu\alpha\nu\beta}=n^{\mu}k^{\nu}q^{\alpha}q^{\beta}C_{\mu\alpha\nu\beta}=-E\left(E_{\mu\nu}e^{\mu}n^{\nu}+\epsilon_{\mu\nu}{}^{\beta}H_{\alpha\beta}n^{\mu}q^{\nu}q^{\alpha}\right),$
(65)

where the first equality follows from the form of the energy-momentum tensor
(2), and the second equality follows from the decomposition of the Weyl
curvature tensor into its electric and magnetic parts with respect to the
matter frame Maartens:1997fg .

It can be seen from Eq. (64) that vanishing on the boundary of each of the
projected Weyl curvature terms of Eq. (65), the off-diagonal components of the
projected volume shear rate of the matter congruence,
$n^{\alpha}e^{\beta}\sigma_{\alpha\beta}$, and the incoming optical shear,
$n^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu}$,151515Any restriction imposed on the
volume shear rate of the matter congruence leads to constraints on the Weyl
curvature and on the spatial gradient of the expansion rate through the
geodesic deviation equation and constraint equations for the relevant
congruence Ellis:2011 . For example,
$n^{\mu}q^{\nu}\hat{\sigma}_{\mu\nu}=0\Rightarrow
k^{\alpha}n^{\mu}q^{\beta}k^{\nu}C_{\alpha(\mu\beta)\nu}=0$ in the present
set-up. By imposing such restrictions it must be checked whether there exist
non-trivial cosmological spacetimes fulfilling such conditions for the volume
shear rate. results in conservation of the boundary term. In this special
case, we have from Eq. (54b) that $\left\langle K\right\rangle_{{{\cal
F}}}\propto 1/A_{{\cal F}}$, with the constant of proportionality depending on
the initial boundary chosen. This is a scaling behaviour of the area-averaged
Gaussian curvature scalar equivalent to that of the boundary-free case. Given
the strict constraints which must be imposed in order to mimic the boundary-
free case, the scaling law $\left\langle K\right\rangle_{{{\cal F}}}\propto
1/A_{{\cal F}}$ is likely not satisfied in realistic cosmological spacetimes.

###  3.3 Area-averaged observables

We shall now perform the area-average over observable quantities for an
observer viewing multiple objects on the same $2$-dimensional screen space.
Using Eqs. (9) and (23), the evolution of the area-average of the logarithmic
cosmological redshift along the observer’s past null cone reads:

$\partial_{\tau}\langle\ln(1+z)\rangle_{{\cal
F}}=-\left\langle\alpha\right\rangle_{{\cal F}}+\text{Cov}_{{\cal
F}}\left(\frac{\hat{\theta}}{E}\,,\,\ln(1+z)\right)\ .$ (66)

Unlike the local cosmological redshift function, which, in general, cannot be
thought of as a time variable along the world lines of the matter congruence,
the area-averaged logarithmic cosmological redshift can be used as a time
label for the $2$-dimensional light fronts as long as the right-hand side of
Eq. (66) does not change sign.

We can, furthermore, compute the evolution along the observer’s past null cone
of the area-average of the logarithm of the effective angular diameter
distance, which, by Eq. (31), reads:161616We omit the normalization of the
dimensionfull quantities $d_{A}$ and $d_{L}$ inside the logarithm; the
corresponding equations hold for any choice of normalization.

$\displaystyle\partial_{\tau}\langle\ln(d_{A})\rangle_{{\cal
F}}=\mathcal{H}_{{\cal
F}}+\int^{\tau}\mathcal{M}(\tau,\tau^{\prime})\,{\mathrm{d}}\tau^{\prime}\ ,$
(67) $\displaystyle{\rm where}\ $
$\displaystyle\mathcal{M}(\tau,\tau^{\prime}):=\frac{1}{2}\text{Cov}_{{\cal
F}}\left(\frac{\hat{\theta}(\tau,x^{A})}{E(\tau,x^{A})}\,,\,\frac{\hat{\theta}(\tau^{\prime},x^{A})}{E(\tau^{\prime},x^{A})}\right)$
(68)

defines the memory function for the observer’s past null cone – cf. the
definition (3.1) – memorizing the spatial inhomogeneities and auto-correlation
properties in the expansion rate and in the energy function of the incoming
geodesic null rays encountered along the path from the source to the observer.
We have used that

$\displaystyle\ln(d_{A})=\frac{1}{2}\int^{\lambda}\hat{\theta}(\lambda^{\prime},x^{A})\,{\mathrm{d}}\lambda^{\prime}=\frac{1}{2}\int^{\tau}\frac{\hat{\theta}(\tau^{\prime},x^{A})}{E(\tau^{\prime},x^{A})}\,{\mathrm{d}}\tau^{\prime}\
,$ (69)

where integration is along each geodesic null ray labelled by comoving
coordinates $x^{A}$, together with Eqs. (31) and (24). The memory of shear and
Weyl curvature along the observer’s past null cone is encoded in the shape of
the boundary $\partial{\cal F}$ of the compact light front domain ${\cal F}$,
according to Eqs. (62) and (63).

The time-range of dependence is expected to be determined by the size of
typical matter structures in the cosmological spacetime. Differentiating Eq.
(67) and using Leibniz’ rule leads to the area-averaged logarithmic focusing
equation:

$\displaystyle\partial^{2}_{\tau}\langle\ln(d_{A})\rangle_{{\cal F}}$
$\displaystyle=$ $\displaystyle\partial_{\tau}\mathcal{H}_{\cal
F}+\mathcal{M}(\tau,\tau)+\int^{\tau}\partial_{\tau}\mathcal{M}(\tau,\tau^{\prime})\,{\mathrm{d}}\tau^{\prime}\
,$ (70a) with $\displaystyle\partial_{\tau}\mathcal{H}_{\cal
F}=-\mathcal{H}^{2}_{{\cal F}}-4\pi\left\langle\varrho\right\rangle_{{\cal
F}}+\mathcal{Q}^{\mathbf{k}}_{{\cal F}}\;+\;{\mathcal{S}}^{\mathbf{k}}_{\cal
F}\ .$ (70b)

The area-averaged logarithmic version of the reciprocity relation (34) reads:

$\langle\ln(d_{L})\rangle_{{\cal F}}=2\langle\ln(1+z)\rangle_{{\cal
F}}+\langle\ln(d_{A})\rangle_{{\cal F}}\ ,$ (71)

whose evolution along the observer’s past null cone can be determined through
relations (66) and (67). There is a sense in which the logarithmic variables
$\ln(1+z)$, $\ln(d_{A})$ and $\ln(d_{L})$ are natural for discussing area-
averaged null signal propagation, as the logarithmic transformation simplifies
evolution equations and preserves the local form of the reciprocity theorem.
We note that the distance modulus is defined as a linear function of
$\ln(d_{L})$ by $\mu:=(5/\ln(10))\ln(d_{L}/10~{}\text{parsecs})$.

We now consider area-averaging of the cosmological redshift drift relation
(35) over sources sprinkled continuously and uniformly in volume over the
$2$-dimensional screen space:

$\left\langle\frac{{\mathrm{d}}z}{{\mathrm{d}}\tau_{0}}\right\rangle_{{\cal
F}}=\langle 1+z\rangle_{{\cal F}}\,\alpha\,\bigg{\rvert}_{P}-\langle
1+z\rangle_{{\cal
F}}\,\frac{E^{\prime}}{E}\,\bigg{\rvert}_{P}-\left\langle\alpha\right\rangle_{{\cal
F}}+\left\langle\frac{E^{\prime}}{E}\right\rangle_{{\cal F}}\ .$ (72)

Spatially inhomogeneous and anisotropic contributions enter the cosmological
redshift drift, also when area-averaged over many sources. In particular,
systematic effects are expected to be introduced through the terms evaluated
at the point of observation $P$, while non-cancelling effects from structure
along the individual null rays might also play a role; see Refs.
Heinesen:2020pms ; Heinesen:2021nrc  for detailed investigations of systematic
effects entering the cosmological redshift drift signal.

###  3.4 FLRW limit

It is worth considering the FLRW limit of the area-averaged equations (45),
(46), (49) and (52). We consider the FLRW spacetime metric written in hyper-
spherical local coordinates $\\{t,r,\theta,\phi\\}$ with line element,

${\mathrm{d}}s^{2}=-{\mathrm{d}}t^{2}+a(t)^{2}({\mathrm{d}}r^{2}+S^{2}_{k}(r){\mathrm{d}}\Omega^{2})\
,\qquad S_{k}(r):=\sqrt{-k}^{-1}\sinh(\sqrt{-k}r)\ ,$ (73)

where
${\mathrm{d}}\Omega^{2}={\mathrm{d}}\theta^{2}+\sin^{2}(\theta){\mathrm{d}}\phi^{2}$
is the solid angular element on the unit sphere, and $a(t)^{2}S^{2}_{k}(r)$ is
the adapted area measure on a sphere of proper radius $a(t)r$ contained in a
spacelike $3$-surface labelled by $t$. The dimensionless scale factor $a(t)$
is the conformal scaling of the static metric, and thus describes the temporal
dependence of the spatial volume in the comoving frame; we employ the usual
convention that $a(t_{0})=1$, where $t_{0}$ is the proper time elapsed since
the big bang singularity at the present epoch (“here and now”). The spatial
coordinates $\\{r,\theta,\phi\\}$ are comoving with a central observer, where
$\\{\theta,\phi\\}$ are angular coordinates describing directions on the
observer’s sky, and $r$ is a radial coordinate with dimension of length from
which we can get the proper geodesic distance $a(t)r$ away from the observer
at an instant of proper time $t$. The curvature of the 3-dimensional spatial
sections $k$ has dimensions of inverse length scale squared, and the spatially
flat FLRW model is obtained when $k\rightarrow 0$. The observer’s time
variable $t$ is synchronous with the proper time function $\tau$ associated
with the irrotational $4$-velocity $u^{\mu}$ of the matter congruence, of
uniform rest mass density $\varrho(t)$. The incoming geodesic null rays which
constitute the central observer’s past null cone are per construction of the
local coordinates propagating along paths with
$\\{\theta=\text{constant},\phi=\text{constant}\\}$. We thus have in this
reference frame:

$u^{\mu}=\delta_{t}{}^{\mu}\ ,\qquad
e^{\mu}=\frac{1}{a(t)}\,\delta_{r}{}^{\mu}\ ,\qquad
k^{\mu}=E\left(\delta_{t}{}^{\mu}-\frac{1}{a(t)}\,\delta_{r}{}^{\mu}\right)\
.$ (74)

The kinematic variables associated with $u^{\mu}$ and $e^{\mu}$ and the energy
function $E$ of the geodesic null congruence – obtained from the connection
associated with the line element (73) and the differential equation (19a) for
$E$ – are:

$\theta=3\,\frac{\dot{a}(t)}{a(t)}\ ,\quad\sigma_{\mu\nu}=0\
,\qquad\overline{\theta}=\frac{2}{a(t)}\frac{\partial_{r}S_{k}(r)}{S_{k}(r)}\
,\quad\overline{\sigma}_{\mu\nu}=0\ ,\qquad E=\frac{\mathcal{E}(L)}{a(t)}\ ,$
(75)

where $\dot{}:=\partial_{t}$, and
$L:=\int^{t}{\mathrm{d}}t^{\prime}/a(t^{\prime})+r$ is constant over each past
null cone of the central observer. The function $\mathcal{E}(L)$ is an initial
condition for the geodesic null congruence and might be chosen in accordance
with the observational frequency of interest; cf. section 2.2. The null scalar
variable of the geodesic null congruence is given by

$V(L)=-\int^{L}\mathcal{E}(L^{\prime})\,{\mathrm{d}}L^{\prime}\ ,$ (76)

which can be verified by computing its gradient and recovering $k^{\mu}$.

The kinematic variables associated with $k^{\mu}$ and $l^{\mu}$ follow from
plugging the expressions (75) into the relations (17a) and (17b), yielding:

$\displaystyle\frac{\hat{\theta}}{E}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\,\theta-\overline{\theta}=2\,\frac{\dot{a}(t)}{a(t)}-\frac{2}{a(t)}\,\frac{\partial_{r}S_{k}(r)}{S_{k}(r)}\
,\qquad\hat{\sigma}_{\mu\nu}=0\ ,$ (77a)
$\displaystyle\frac{\tilde{\theta}}{E}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\,\theta+\overline{\theta}=2\,\frac{\dot{a}(t)}{a(t)}+\frac{2}{a(t)}\,\frac{\partial_{r}S_{k}(r)}{S_{k}(r)}\
,\qquad\tilde{\sigma}_{\mu\nu}=0\ .$ (77b)

The contribution $\pm 2\partial_{r}(S_{k}(r))/S_{k}(r)/a(t)$ to the area
expansion rates accounts for the change of the area measure as the proper
radius of the $2$-dimensional screen space is decreased or increased,
respectively, when propagating the screen towards or away from the vertex
point $P$ of the observer’s past null cone along the incoming/outgoing null
rays. This contribution becomes singular at the vertex point $P$ at coordinate
value $r=0$, where the angular measure of the $2$-dimensional screen space
tends to zero.

We now consider the area-averaged equations (46, 49), in the limit of the FLRW
cosmologies. In this limit, the backreaction terms
$\mathcal{Q}^{\mathbf{k}}_{{\cal F}}$ and
$\mathcal{Q}^{\mathbf{k}\mathbf{l}}_{{\cal F}}$ and the right-hand side of Eq.
(52b) vanish. The Ricci scalar of the embedded $2$-surfaces is

$\displaystyle{}^{(2)}\\!R=\frac{{}^{(3)}\\!R}{3}+\frac{1}{2}\,\overline{\theta}^{2}\
;\qquad{}^{(3)}\\!R=\frac{6k}{a^{2}(t)}\ ,$ (78)

where ${}^{(3)}\\!R$ is the 3-Ricci curvature scalar associated with the
canonical FLRW spacelike $3$-surfaces ${\cal S}(t)$. Using these results,
along with Eqs. (45) and (77a), on both sides of the area-averaged Gauss
embedding constraint (49) gives:

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}+\frac{1}{4}\,\overline{\theta}^{2}-\frac{\dot{a}}{a}\,\overline{\theta}=\frac{8\pi}{3}\varrho+\frac{1}{3}\Lambda-\frac{k}{a^{2}}-\frac{\dot{a}}{a}\,\overline{\theta}+\frac{1}{4}\,\overline{\theta}^{2}\
,$

where it has been used that all area-averaged variables are constant over the
$2$-dimensional screen space. The terms involving $\overline{\theta}$ cancel,
and we arrive at the first of Friedmann’s equations:

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi}{3}\varrho+\frac{1}{3}\Lambda-\frac{k}{a^{2}}\
.$ (79)

The Raychaudhuri equation can be recovered by noticing that

$\displaystyle
E^{-1}k^{\rho}\partial_{\rho}\left(\frac{\hat{\theta}}{E}\right)$
$\displaystyle=$
$\displaystyle(u^{\rho}-e^{\rho})\partial_{\rho}\left(2\,\frac{\dot{a}}{a}-\overline{\theta}\right)=2\,\frac{\ddot{a}}{a}-2\left(\frac{\dot{a}}{a}\right)^{2}-\dot{\overline{\theta}}+\frac{1}{a}\,\overline{\theta}^{\prime}\
,$ (80a) $\displaystyle{\rm where}\qquad\
\dot{\overline{\theta}}=-\frac{\dot{a}}{a}\,\overline{\theta}\
,\qquad\overline{\theta}^{\prime}=-\frac{2k}{a}-\frac{a}{2}\,\overline{\theta}^{2}\
,$ (80b)

follows from the definition of $\overline{\theta}$, Eq. (75). This gives for
the acceleration equation (46):

$\displaystyle\frac{\partial^{2}_{\tau}a_{{\cal F}}(\tau)}{a_{{\cal
F}}(\tau)}$ $\displaystyle=$
$\displaystyle\frac{\ddot{a}}{a}-\frac{1}{2}\,\frac{\dot{a}}{a}\,\overline{\theta}-\frac{k}{a^{2}}\
.$ (81)

We also have from Eq. (46) that

$\displaystyle\frac{\partial^{2}_{\tau}a_{{\cal F}}(\tau)}{a_{{\cal
F}}(\tau)}=-4\pi\varrho+\left(\frac{\dot{a}}{a}\right)^{2}-\frac{1}{2}\,\overline{\theta}\,\frac{\dot{a}}{a}=-\frac{4\pi}{3}\varrho+\frac{1}{3}\Lambda-\frac{k}{a^{2}}-\frac{1}{2}\,\overline{\theta}\,\frac{\dot{a}}{a}\
,$ (82)

where it has been used that

$\displaystyle\frac{\hat{\theta}}{E}\,\frac{E^{-1}k^{\rho}\partial_{\rho}E}{E}=-2\left(\frac{\dot{a}}{a}\right)^{2}+\overline{\theta}\,\frac{\dot{a}}{a}\
.$ (83)

Finally, equating Eqs. (81) and (82), we arrive at Friedmann’s acceleration
law:

$\displaystyle\frac{\ddot{a}}{a}=-\frac{4\pi}{3}\varrho+\frac{1}{3}\Lambda\ .$
(84)

It is illustrative to derive the scaling of the Gaussian curvature scalar of
the FLRW screen spaces:

$\displaystyle
K=\frac{1}{2}{}^{(2)}\\!R=\frac{1}{a^{2}(t)S^{2}_{k}(r)}\left[\,kS^{2}_{k}(r)+\left(\partial_{r}S_{k}(r)\right)^{2}\,\right]=\frac{1}{a^{2}(t)S^{2}_{k}(r)}\
.$ (85)

This expression is manifestly positive when both $a(t)$ and $S_{k}(r)$ are
non-zero. The first equality in Eq. (85) follows from Eqs. (78) and (77a),
while the second follows from the definition of $S_{k}(r)$ and hyperbolic
identities. The FLRW expression for the Gaussian curvature scalar (85) is
nothing but the inverse area measure of the light fronts. The integral of Eq.
(85) over the full light front yields the area of the unit sphere, $4\pi$,
consistent with the Gauss–Bonnet theorem for closed surfaces (54a). For light
front subdomains, the integral of $K$ is simply the constant area covered by
the subdomain of the unit sphere.

##  4 Discussion

This paper offers an area-averaging formalism for describing the evolution of
dynamical variables and observable quantities along the observer’s past null
cone in the setting of generic irrotational dust cosmologies with a
cosmological constant. This formalism applies to all kinds of null signals,
e.g., those carried by photons, but also a future gravitational wave sky could
be examined using this formalism. The formalism is based on the assumption
that the cosmological dynamics progresses sufficiently smoothly so that no
caustics will form on the past null cone of the observer.

The area-averaged macroscopic system of equations on $2$-dimensional light
fronts derived in section 3.1 involves two independent backreaction terms. The
system of light front averaged equations comprises a larger set of global
variables than the analogous volume-averaged macroscopic system of equations
adapted to the spacelike $3$-surfaces given in Ref. buchert:grgdust . This is
to be expected, as additional dynamical variables are introduced when the
dynamics of the null cone generating congruence is included. Since the area-
averaging operation is adapted to the observer’s past null cone, there is a
flow of matter world lines across the averaging domain, which further
complicates the light front averaged equations. We note that caution is
required when invoking simplifying assumptions. A radical example would be to
get rid of all of the shear variables of the problem by simply setting the
area shear rate of the geodesic null congruence and the volume shear rate of
the matter congruence to zero. However, such approximations at the local level
of the cosmological spacetime turn out to be extremely restrictive; see Ref.
Ellis:2011  and references therein. If, for instance, the volume shear rate of
the matter congruence is required to be zero, $\sigma_{\mu\nu}=0$, then –
since its vorticity is already required to be zero – this implies that the
dust cosmology considered must be of the FLRW class ell1967 . Vanishing of the
area shear rate of the geodesic null congruence, on the other hand, implies
that the evolution of this null congruence is affected only by the matter
(hence, Ricci curvature) it encounters,171717The Weyl curvature has no
influence on the properties of the geodesic null congruence in this case. and
that this null congruence must be a principal null direction of the Weyl
curvature tensor; this fact might be viewed as a generalization of the
Goldberg–Sachs theorem for vacuum spacetimes Ellis:2011 ; GoldbergSachs .

The generality of the presented formalism leads, by construction, to a set of
balance equations that do not form a closed set. Quantification of the level
of backreaction is intimately related to the implementation of closure
conditions for the screen-space averaged variables. Concrete inhomogeneous
models at the level of area-averaged variables can be studied in order to
arrive at appropriate statistical descriptions for the propagation of geodesic
null ray bundles in spatially inhomogeneous dust cosmologies. One might for
instance employ exact solutions such as Lemaître–Tolman–Bondi solutions or the
more general classes of Szekeres solutions. They can be employed to design
‘Swiss cheese models’, within which light propagation has already been studied
in great detail; see Lavinto:2013exa ; Koksbang:2021zyi  for recent
investigations into observational backreaction effects in Swiss cheese models.
Generic structure formation models, based on perturbative assumptions and not
restricted by symmetries, are a natural next step of investigation. In the
present setup of the irrotational matter congruence, appropriate models have
been constructed that are based on relativistic perturbation theory combined
with exact averages, e.g. lischwarz , and the relativistic generalization of
Lagrangian perturbation theory, recently reviewed in Universe .181818See the
references therein, also to earlier work on other relativistic perturbation
theories, the history of Lagrangian perturbation theory, and e.g. RZA_2  for
their application as a closure condition for the $(3+1)-$averaged system.
Ongoing work investigates a non-perturbative generalization of these schemes,
controlled by the Szekeres class of solutions and that contain the most
general Szekeres solutions as the exact body GRZA .

Studies of exact scaling solutions of the area-averaged system could also
provide more insight, similar to those employed in the $3+1$ slicing formalism
roy:instability , and they can be explored to provide closure conditions for
the area-averaged equations. The scaling behaviours suggested by the
integrability condition (52c) at first glance are those where the relevant
variables are uncoupled, but we expect from Ref. roy:instability  that generic
scaling solutions dynamically couple the area-averaged variables. In this
context, topological closure conditions might be relevant: an insight from the
Gauss–Bonnet theorem applied to light fronts shows that there is no such
coupling of the scalar curvature of light fronts to other dynamical variables
for all-sky averages, unlike in the $3$-dimensional case: the area-averaged
scalar parts of Einstein’s field equations on the light fronts simplify
through the Gauss–Bonnet theorem when the total (all-sky) light front is
considered, or when the curvature of the embedding of the boundary is
associated with a conserved current, as detailed in section 3.2. Integral-
geometric properties analogous to the Minkowski functionals in Euclidian
spaces appear naturally when all-sky averaging over light fronts is performed.
Hence, it is for the case of all-sky averages, where we could expect
simplifications after evaluating global contributions on topologically closed
light fronts. We do not in general expect the light front backreaction terms
to vanish globally. However, in cases where those terms are covariant
divergences, a closed space would erase them, as in the case of flat light
fronts (corresponding insights have been developed in Newtonian cosmology
buchertehlers ).

Another arena of application of the formalism is that of general relativistic
numerical simulations. Cosmological simulations have already been investigated
in the context of backreaction and observables; cf. Giblin:2015vwq ;
Giblin:2016mjp ; Bentivegna:2015flc ; Adamek:2015eda ; East:2017qmk ;
East:2019chx ; Macpherson:2018btl  for recent works. In these works,
cosmological backreaction was generally found to be small at large scales
where observables were also found to agree well with the FLRW prediction,
whereas on smaller scales Newtonian simulations of structures combined with
relativistic ray tracing could account for distortions of the FLRW background
observations. The FLRW metric has thus proven itself robust towards the
initial perturbations employed within these simulated universe scenarios. This
robustness might seem surprising given the highly non-linear structures which
are allowed to develop in these codes. However, this may be related to the
fact that a foliation that inherits the properties of the longitudinal gauge
of standard perturbation theory is present throughout the cosmic evolution of
these simulations despite the presence of non-linear structures Giblin:2018ndw
; Clifton:2020oqx . Furthermore, the global architecture remains that of
Newtonian cosmology of a flat $3$-torus topology, which forces the global
model to evolve according to the assumed background cosmology. These are
conjectures that should be examined in future works, together with the
question of which type of initial perturbations and topological constraints
may cause the FLRW solution to be globally unstable as a model for the average
evolution, an instability that generically occurs in a dynamical system
analysis of the averaged $3+1$ system roy:instability , pointing to negative
average curvature on large scales buchertcarfora:curvature . The formalism
developed in this paper sheds light on the conditions that must be satisfied
for backreaction to be present in the context of observables, and may lead to
interesting developments on simulated spacetime scenarios with backreaction.

Turning to observables, the evolution equations for the area-averaged
effective angular diameter distance derived in section 3.3 are partly
expressed in terms of the area-averaged light front variables introduced in
section 3.1, but involve additional terms expressed through the memory
function $\mathcal{M}(\tau,\tau^{\prime})$, for which one needs to prescribe a
model to close the area-averaged macroscopic system of equations. Backreaction
in the area-averaged effective angular diameter distance might arise due to
‘direct backreaction’ through the memory function or
$\mathcal{Q}^{\mathbf{k}}_{{\cal F}}$, or through ‘indirect backreaction’
causing both $\mathcal{H}_{{\cal F}}$ and
$\left\langle\varrho\right\rangle_{{\cal F}}$ to evolve differently than in
the Friedmannian case. The evolutions along the observer’s past null cone of
the area-averaged cosmological redshift and of the area-averaged cosmological
redshift drift both induce additional terms which must be dealt with – for
detailed investigations of the expression for cosmological redshift drift in a
generic cosmological spacetime, and the effect of regional spatial
inhomogeneities and anisotropies on the cosmological redshift drift signal,
see Refs. Heinesen:2020pms ; Heinesen:2021nrc ; Heinesen:2021qnl ;
Korzynski:2017nas .

We emphasize that care must be taken when using the area-averaged relations to
interpret cosmological data sets. Some challenges for interpreting the
theoretically given area-averaged observable quantities include (i)
identifying spacelike $3$-surfaces of constant proper time on which the area-
averaging operation is formulated; (ii) relating the effective area-averaging
operation employed when observing many emitters over the sky to the volume-
averaging operation. Emitters might not sample the area of the emission screen
space fairly, and null signals will interact with the matter of the spacetime;
(iii) incorporating the properties of the observer’s local vicinity relative
to that of an observer with a volume-unbiased view of the Universe. Here
follows a discussion of these challenges.

  * (i)

Relative age measurements JimenezLoeb  might serve to distinguish spacelike
$3$-surfaces of constant proper time, although the errors associated with
current relative age measurements are of the order of 20% Ratsimbazafy . The
foliation of the observer’s $3$-dimensional past null cone is not unique, and
the choice to foliate it in terms of a family of spacelike $3$-surfaces ${\cal
S}:\\!\\{\tau=\text{constant}\\}$ of constant proper time with respect to the
irrotational matter congruence could be modified, which would change the
dynamical system of equations accordingly. It has been proposed to consider
level surfaces of constant cosmological redshift $z$ as an observationally
relevant foliation of the observer’s past null cone Fanizza2019pfp . Regions
of constant cosmological redshift are identifiable to a high precision, and
are therefore observationally preferred. However, when $z$ is not a monotonic
function along null paths, this destroys the otherwise observationally
intuitive notion of $z$ as a foliation scalar. In situations where the
effective fluid matter model is describing scales well above that of
collapsing structures, then the redshift function might be re-established as a
parametrization of light fronts. It must be assessed in the physical situation
at hand whether the cosmological redshift function is indeed a hypersurface-
forming scalar. In the general case, where the cosmological redshift function
is not hypersurface-forming, we might, nevertheless, average over regions of
constant cosmological redshift, without these regions being well-defined as
causally ordered and non-intersecting submanifolds.

  * (ii)

Luminous astrophysical sources emitting electromagnetic radiation or
gravitational wave signals are in general not uniformly distributed in volume
over the Universe. Such sources are concentrated in over-densities of the
matter distribution, and in this regard an observer should be observing a
mass-biased picture of the Universe. On the other hand, after being emitted
from a source, the null signals, which a given observer receives, tend to
propagate in empty space (since signals propagating in regions with matter are
likely to be scattered or absorbed by other particles, and, therefore, simply
do not arrive at the observer’s position). Such subtleties must ideally be
taken into account in an observationally matched averaging procedure. However,
correctly accounting for such biases would be involved and the correction
procedures would depend on the type of cosmological probe. As a first
approximation we might simply assume an ideal situation where sources and null
signals are probing the volume-averaged cosmological spacetime.

  * (iii)

All cosmological measurements are dependent on the properties that hold in the
immediate vicinity of the small segment of a single world line from which we
observe the Universe. Whereas cosmological observations are directly or
indirectly averaged over the emitter positions when performing statistical
analysis with many cosmological data points, the same does not apply to the
observer’s position. Some cosmological measurements might be sensitive to the
observer’s position, whereas others might be subject to less bias.

The area-averaging formalism provided in this paper may be combined with
appropriate statistical assumptions in order to formulate an average dynamical
theory for dynamics and observables on our past null cone.

Synthetic multi-variate cosmological data on the joint distribution of
luminous astrophysical sources and spatial geometry inside and on a modelled
observer’s past null cone may be generated by simulation and analyzed by means
of methods from inductive statistical inference; cf. the textbook
introductions, Refs. mce2020 ; geletal2014 , and hve2020 . This approach
allows for the possibility of factoring in all information that is relevant to
dealing systematically with the uncertainty pertaining to a quantitative
problem at hand. In this way (and when linked to exact cosmological dust
solutions to begin with), a set of interval estimates for area-averaged
observables such as, e.g., the cosmological distance measures may be
calculated for specific configurations and then compared to standard FLRW-
based values for these quantities. The amount of available observational data
to be integrated in this process is steadily increasing. On the technical
side, a long term project of this kind of cosmological data analysis may be
informed by and validated through supervised solution algorithms developed in
contemporary machine learning such as neural networks; cf., e.g., the textbook
introduction Ref. ng2018 , and examples on applications in Refs. gabetal2018 ;
machinelearning:lightcone .

Our analytical results can be profitably employed for numerical simulations of
null cones. As mentioned in the discussion above, it is important to examine
the regime of global stability of the FLRW solution for predicting observables
by means of numerical relativity. Relativistic effects in light propagation
have been computed in linearized weak-field approximations of general
relativity, while using (N-body or hydrodynamic) Newtonian simulations as
input to generate the matter distribution of the universe model
Borzyszkowski:2017ayl ; Breton:2018wzk . Light rays have also been tracked
fully relativistically through idealized but non-trivial inhomogeneous and
anisotropic model universes with backreaction Lavinto:2013exa ;
Koksbang:2019glb ; Koksbang:2020zej  and within post-Newtonian model universes
Sanghai:2017yyn ; Grasso:2021zra . Weak-field relativistic N-body simulations,
which allow to treat certain perturbative modes non-linearly, have been
employed in order to model both the matter distribution and the null cone of
an observer Adamek:2018rru . Null cones have been generated in general
relativistic simulations codes; however, only tracking the evolution in the
linear regime of density contrasts Giblin:2016mjp . See also East:2017qmk  for
a study comparing observables in Newtonian N-body codes and general
relativistic hydrodynamic simulations of structure formation, while employing
relativistic ray-tracing in both cases. Finally, ray tracing codes for the
general relativistic and non-linear universe simulations in Refs.
Macpherson:2016ict ; Macpherson:2018btl  have recently been developed
Macpherson:2022eve  with exciting applications to relativistic modelling of
cosmic observables Macpherson:2022 .

The framework developed in this paper is complementary to model-independent
cosmographic frameworks first considered in detail by krisac1966 ;
1985PhR…124..315E  and more recently by Clarkson:2011uk ; Heinesen:2020bej ;
Heinesen:2021nrc ; Heinesen:2021qnl , with applications to numerical
simulations in Macpherson:2021gbh ; Heinesen:2021azp . While the cosmography
developed in these works is powerful for extracting cosmological information
from data without the assumption of a cosmological metric or a dynamical
theory, and primarily apply to measurements made at low redshifts, the
framework considered in this paper considers the dynamical evolution and
theoretical average of light fronts as they propagate towards the observer
from astrophysical sources that may be distant. The cosmographic frameworks in
Heinesen:2020bej ; Heinesen:2021qnl  are furthermore concerned with the
angular dependence of observables over the sky of the observer, while the
framework developed here is concerned with the average over (a patch of) the
sky of the observer, which is applicable in general for finite-area screen
spaces with or without boundary. Future investigations may consider a link to
the analyses of galaxy catalog data on shells around the observer or Cosmic
Microwave Background radiation data that are often investigated globally on
topological all-sky support manifolds with the help of integral-geometric
measures like the Minkowski functionals.

###### Acknowledgements.

This work is part of a project that has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement ERC advanced grant 740021–ARTHUS, PI:
TB). The authors would like to thank Jürgen Ehlers, George Ellis, Pierre
Mourier, Dominik Schwarz and Nezihe Uzun for useful discussions, and Giuseppe
Fanizza, Syksy Räsänen and Dennis Stock for valuable comments on the
manuscript. We would in addition like to thank the anonymous referee for
constructive suggestions that helped improving the paper. This work has been
begun during a visit of TB in 2007 to the University of Bielefeld, Germany. TB
wishes to thank Dominik Schwarz for his invitation to hold a temporary
C4-chair at the department of physics. HvE acknowledges the generous
hospitality of the UCT Cosmology and Gravity Group during the period from July
to October 2010 when a significant share of the work underlying this paper was
accomplished.

##  Appendix A Light front-adapted spacetime metric

We shall consider the embedding of the light fronts into the cosmological
spacetime, and write the metric tensor in local coordinates that are adapted
to the screen space, $x^{\mu}=(\tau,V,x^{A})$. We require that the local
coordinates $x^{A}$, with $A=2,3$, satisfy the propagation law
$k^{\rho}\partial_{\rho}x^{A}=0$. The propagation laws for $\tau$ and $V$ are
also fixed through $k^{\rho}\partial_{\rho}\tau=E$ and
$k^{\rho}\partial_{\rho}V=0$ (see Eqs. (24) and (7)), and the full adapted
system of local coordinates $x^{\mu}=(\tau,V,x^{A})$ is thus specified
throughout the past null cone domain, once initial conditions for the local
coordinates are fixed at a screen space. In these local coordinates, we have
that

$\displaystyle k_{\mu}=(0,{1},0,0)\ ;$ $\displaystyle u_{\mu}=({-1},0,0,0)\ ;$
(A.1) $\displaystyle k^{\mu}=({E},0,0,0)\quad;$ $\displaystyle
u^{\mu}=({1,-E,EU^{A}})\ ,$ (A.2)

where $U^{A}:=u^{\rho}\partial_{\rho}(x^{A})/E$ defines the drift of the
screen space coordinates in the matter frame. The one-form components in (A.1)
follow directly from the definitions $k_{\mu}:=\partial_{\mu}V$ and
$u_{\mu}:=-\partial_{\mu}\tau$. The vector components in (A.2) follow from the
definition of the energy function $E:=-k^{\mu}u_{\mu}$ and the transport rules
(shift vectors) $k^{\rho}\partial_{\rho}(x^{A})=0$ and
$u^{\rho}\partial_{\rho}(x^{A})=EU^{A}$. We may now write the projection
tensor onto the light fronts (11) in the adapted local coordinate system
$x^{\mu}=(\tau,V,x^{A})$:

$p_{\mu\nu}=\left(\begin{array}[]{ccc}0&0&\vec{0}\\\ 0&\,\,U^{2}&U_{B}\\\
\vec{0}&U_{A}&p_{AB}\\\ \end{array}\right)\
;~{}~{}~{}~{}~{}~{}~{}p^{\mu\nu}=\left(\begin{array}[]{ccc}0&0&\vec{0}\\\
0&\,\,0&\vec{0}\\\ \vec{0}&\vec{0}&p^{AB}\\\ \end{array}\right)\ ,$ (A.3)

where $U_{A}:=p_{AB}U^{B}$ and $U^{2}:=p_{AB}U^{A}U^{B}$, and where the area-
adapted screen space metric $p_{AB}$ has inverse $p^{AB}$. The tensor
components in (A.3) follow from the orthogonality conditions
$p_{\mu\nu}k^{\nu}=p_{\mu\nu}u^{\nu}=0$ and
$p^{\mu\nu}k_{\nu}=p^{\mu\nu}u_{\nu}=0$, respectively. Note that in general
the values of the components $p_{11}$, $p_{1A}$ and $p_{A1}$ are non-zero,
which comes from generally non-zero values for $u^{1}$ and $u^{A}$ in (A.2).

Using the definitions (10a) and (11), we might formulate the metric tensor for
the cosmological spacetime as
$g_{\mu\nu}=k_{\mu}k_{\nu}/E^{2}-(k_{\mu}u_{\nu}+u_{\mu}k_{\nu})/E+p_{\mu\nu}$,
and insert Eqs. (A.1), (A.2), and (A.3) in this formulation to obtain

${g_{\mu\nu}=\left(\begin{array}[]{ccc}0&1/E&\vec{0}\\\
1/E&\,(1/E^{2})\\!+\\!U^{2}&U_{B}\\\ \vec{0}&{U_{A}}&p_{AB}\\\
\end{array}\right)\
;~{}~{}~{}~{}~{}~{}~{}g^{\mu\nu}=\left(\begin{array}[]{ccc}-1&\,E&{-}EU^{B}\\\
E&0&\vec{0}\\\ {-}EU^{A}&\vec{0}&p^{AB}\\\ \end{array}\right)\ .}$ (A.4)

Data Availability statement:

Data sharing not applicable to this article as no datasets were generated or
analysed during the current study.

## References

  * (1) Adamek, J., Clarkson, C., Coates, L., Durrer, R., Kunz, M.: Bias and scatter in the Hubble diagram from cosmological large-scale structure. Phys. Rev. D 100, 021301 (2019). [arXiv:1812.04336]
  * (2) Adamek, J., Daverio, D., Durrer, R., Kunz, M.: General relativity and cosmic structure formation. Nature Phys. 12, 346 (2016). [arXiv:1509.01699]
  * (3) Aghanim, N. et al. (Planck Collaboration). Planck 2018 results. VI. Cosmological parameters. Astron. Astrophys. 641, A6 (2020). [arXiv:1807.06209]
  * (4) Arnowitt, R., Deser, S., Misner, C.W.: The dynamics of general relativity. In Gravitation, edited by L. Witten (Wiley, New York, 1962), 227. Reprinted: Gen. Relativ. Gravit. 40, 1997 (2008). [arXiv:gr-qc/0405109]
  * (5) Bagheri, S., Schwarz, D.J.: Light propagation in the averaged universe. J. Cosmol. Astropart. Phys. 10, 073 (2014). [arXiv:1404.2185]
  * (6) Beenakker, W. Venhoek, D.: A structured analysis of Hubble tension. [arXiv:2101.01372]
  * (7) Ben-Dayan, I., Gasperini, M., Marozzi, G., Nugier, F., Veneziano, G.: Backreaction on the luminosity-redshift relation from gauge invariant light-cone averaging. J. Cosmol. Astropart. Phys. 04, 036 (2012). [arXiv:1202.1247]
  * (8) Ben-Dayan, I., Gasperini, M., Marozzi, G., Nugier, F., Veneziano, G.: Average and dispersion of the luminosity-redshift relation in the concordance model. J. Cosmol. Astropart. Phys. 06, 002 (2013). [arXiv:1302.0740]
  * (9) Ben-Dayan, I., Gasperini, M., Marozzi, G., Nugier, F., Veneziano, G.: Do stochastic inhomogeneities affect dark-energy precision measurements? Phys. Rev. Lett. 110, 021301 (2013). [arXiv:1207.1286]
  * (10) Bengtsson, I.: The Hawking energy on photon surfaces. Gen. Relativ. Gravit. 52, 52 (2020). [arXiv:2001.08491]
  * (11) Bentivegna, E., Bruni, M.: Effects of nonlinear inhomogeneity on the cosmic expansion with numerical relativity. Phys. Rev. Lett. 116, 251302 (2016). [arXiv:1511.05124]
  * (12) Bildhauer, S., Futamase, T.: The cosmic microwave background in a globally inhomogeneous universe. Mon. Not. R. Astron. Soc. 249, 126 (1991).
  * (13) Biswas, T., Mansouri, R., Notari, A.: Nonlinear structure formation and apparent acceleration: an investigation. J. Cosmol. Astropart. Phys. 12, 017 (2007). [arXiv:astro-ph/0606703]
  * (14) Bolejko, K.: Volume averaging in the quasispherical Szekeres model. Gen. Relativ. Gravit. 41, 1585 (2009). [arXiv:0808.0376]
  * (15) Bolejko, K.: Emerging spatial curvature can resolve the tension between high-redshift CMB and low-redshift distance ladder measurements of the Hubble constant. Phys. Rev. D 97, 103529 (2018). [arXiv:1712.02967]
  * (16) Bonvin, C., Durrer, R., Gasparini, M.A.: Fluctuations of the luminosity distance. Phys. Rev. D 73, 023523 (2006). [arXiv:astro-ph/0511183]
  * (17) Borzyszkowski, M., Bertacca, D., Porciani, C.: LIGER: mock relativistic light-cones from Newtonian simulations. Mon. Not. R. Astron. Soc. 471, 3899 (2017). [arXiv:1703.03407]
  * (18) Brady, P.R., Droz, S., Israel, W., Morsink, S.M.: Covariant double-null dynamics: 2 + 2-splitting of the Einstein equations. Class. Quantum Gravity 13, 2211 (1996). [arXiv:gr-qc/9510040]
  * (19) Breton, M.A., Rasera, Y., Taruya, A., Lacombe, O., Saga, S.: Imprints of relativistic effects on the asymmetry of the halo cross-correlation function: from linear to non-linear scales. Mon. Not. R. Astron. Soc. 483, 2671 (2019). [arXiv:1803.04294]
  * (20) Brouzakis, N., Tetradis, N., Tzavara, E.: The effect of large-scale inhomogeneities on the luminosity distance. J. Cosmol. Astropart. Phys. 02, 013 (2007). [arXiv:astro-ph/0612179]
  * (21) Bruni, M., Sonego, S.: Observables and gauge invariance in the theory of non-linear spacetime perturbations. Class. Quantum Gravity 16, L29 (1999). [arXiv:gr-qc/9906017]
  * (22) Brunswic, L., Buchert, T.: Gauss–Bonnet–Chern approach to the averaged Universe. Class. Quantum Gravity 37, 215022 (2020). [arXiv:2002.08336]
  * (23) Buchert, T.: On average properties of inhomogeneous cosmologies. In: 9th JGRG Meeting, Hiroshima 1999, Y. Eriguchi et al. (eds.), _J.G.R.G._ 9, 306 (2000). [arXiv:gr-qc/0001056]
  * (24) Buchert, T.: On average properties of inhomogeneous fluids in general relativity: dust cosmologies. Gen. Relativ. Gravit. 32, 105 (2000). [arXiv:gr-qc/9906015]
  * (25) Buchert, T.: On average properties of inhomogeneous fluids in general relativity: perfect fluid cosmologies. Gen. Relativ. Gravit. 33, 1381 (2001). [arXiv:gr-qc/0102049]
  * (26) Buchert, T.: Dark Energy from structure: a status report. Gen. Relativ. Gravit. 40, 467 (2008). [arXiv:0707.2153]
  * (27) Buchert, T.: A cosmic equation of state for the inhomogeneous universe: can a global far-from-equilibrium state explain dark energy? Class. Quantum Gravity 22, L113 (2005). [arXiv:gr-qc/0507028]
  * (28) Buchert, T.: Toward physical cosmology: focus on inhomogeneous geometry and its non-perturbative effects. Class. Quantum Gravity 28, 164007 (2011). [arXiv:1103.2016]
  * (29) Buchert, T., Carfora, M.: Regional averaging and scaling in relativistic cosmology. Class. Quantum Gravity 19, 6109 (2002). [arXiv:gr-qc/0210037]
  * (30) Buchert, T., Carfora, M.: On the curvature of the present-day Universe. Class. Quantum Gravity 25, 195001 (2008). [arXiv:0803.1401]
  * (31) Buchert, T., Coley, A., Kleinert, H., Roukema, B.F., Wiltshire, D.L.: Observational challenges for the standard FLRW model. Int. J. Mod. Phys. D 25, 1630007 (2016). [arXiv:1512.03313]
  * (32) Buchert, T., Delgado Gaspar, I., Ostrowski, J.J.: On General-Relativistic Lagrangian Perturbation Theory and Its Non-Perturbative Generalization. Universe 8, 583 (2022). [arXiv:2209.13417]
  * (33) Buchert, T., Ehlers, J.: Averaging inhomogeneous Newtonian cosmologies. Astron. Astrophys.  320, 1 (1997). [arXiv:astro-ph/9510056]
  * (34) Buchert, T., Larena, J., Alimi, J.-M.: Correspondence between kinematical backreaction and scalar field cosmologies: The ‘morphon field’. Class. Quantum Gravity 23, 6379 (2006). [arXiv:gr-qc/0606020]
  * (35) Buchert, T., Mourier, P., Roy, X.: Cosmological backreaction and its dependence on spacetime foliation. Class. Quantum Gravity 35, 24LT02 (2018). [arXiv:1805.10455]
  * (36) Buchert, T., Mourier, P., Roy, X.: On average properties of inhomogeneous fluids in general relativity III: general fluid cosmologies. Gen. Relativ. Gravit. 52, 27 (2020). [arXiv:1912.04213]
  * (37) Buchert, T., Nayet, C., Wiegand, A.: Lagrangian theory of structure formation in relativistic cosmology. II. Average properties of a generic evolution model. Phys. Rev. D 87, 123503 (2013). [arXiv:1303.6193]
  * (38) Buchert, T., Räsänen, S.: Backreaction in Late-Time Cosmology. Annu. Rev. Nucl. Part. Sci. 62, 57 (2012). [arXiv:1112.5335]
  * (39) Carfora, M., Familiari, F.: A comparison theorem for cosmological lightcones. Lett. Math. Phys. 111, 53 (2021). [arXiv:2101.12698]
  * (40) Chirinos Isidro, E.G., Barbosa, R.M., Piattella, O.F., Zimdahl, W: Averaged Lemaître-Tolman-Bondi dynamics. Class. Quantum Gravity 34, 035001 (2017). [arXiv:1608.00452]
  * (41) Clarkson, C., Umeh, O.: Is backreaction really small within concordance cosmology? Class. Quantum Gravity 28, 164010 (2011). [arXiv:1105.1886]
  * (42) Clifton, T., Gallagher, C.S., Goldberg, S., Malik, K.A.: Viable gauge choices in cosmologies with nonlinear structures. Phys. Rev. D 101, 063530 (2020). [arXiv:2001.00394]
  * (43) Delgado Gaspar, I., Buchert, T., Ostrowski, J.J.: Beyond relativistic Lagrangian perturbation theory. I. An exact-solution controlled model for structure formation. Phys. Rev. D , accepted. [arXiv:2210.04004]
  * (44) Demianski, M., de Ritis, R., Marino, A.A., Piedipalumbo, E.: Approximate angular diameter distance in a locally inhomogeneous universe with nonzero cosmological constant. Astron. Astrophys. 411, 33 (2003). [arXiv:astro-ph/0310830]
  * (45) Di Valentino, E., Mena, O., Pan, S., Visinelli, L., Yang, W., Melchiorri, A., Mota, D.F., Riess, A.G., Silk, J.: In the realm of the Hubble tension—a review of solutions. Class. Quantum Gravity 38, 153001 (2021). [arXiv:2103.01183]
  * (46) Dyer, C.C., Roeder, R.C.: Observations in locally inhomogeneous cosmological models. Astrophys. J. 189, 167 (1974).
  * (47) East, W.E., Wojtak, R., Abel, T.: Comparing fully general relativistic and Newtonian calculations of structure formation. Phys. Rev. D 97, 043509 (2018). [arXiv:1711.06681]
  * (48) East, W.E., Wojtak, R., Pretorius, F.: Einstein–Vlasov calculations of structure formation. Phys. Rev. D 100, 103533 (2019). [arXiv:1908.05683]
  * (49) Ehlers, J.: Beiträge zur relativistischen Mechanik kontinuierlicher Medien. Akad. Wiss. Lit. Mainz, Abhandl. Math.–Nat. Kl. 11, 793–837 (1961). English translation: Contributions to the relativistic mechanics of continuous media. Gen. Relativ. Gravit. 25, 1225–1266 (1993).
  * (50) Einstein, A.: Die Feldgleichungen der Gravitation. Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.), 844–847 (1915). English translation: einsteinpapers.press.princeton.edu/vol6-trans/129.
  * (51) Einstein, A.: Kosmologische Betrachtungen zur allgemeinen Relativitätstheorie. Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.), 142–152 (1917). English translation: Cosmological considerations on the general theory of relativity. In The Principle of Relativity, (London, Methuen, 1923). Reprinted: (Mineola, Dover Publications, 1952), 175–188. einsteinpapers.press.princeton.edu/vol6-trans/433.
  * (52) Ellis, G.F.R.: Dynamics of pressure–free matter in general relativity. J. Math. Phys. 8, 1171–1194 (1967).
  * (53) Ellis, G.F.R.: Relativistic cosmology. In General Relativity and Cosmology, Proc. Int. School of Physics “Enrico Fermi” (Varenna), Course XLVII, edited by R. K. Sachs (Academic Press, New York, 1971), 104–182. Reprinted: Gen. Relativ. Gravit. 41, 581–660 (2009).
  * (54) Ellis, G.F.R.: Relativistic cosmology — its nature, aims and problems. In General Relativity and Gravitation, edited by B. Bertotti, F. de Felice and A. Pascolini, (D. Reidel Publishing Co., Dordrecht, 1984), 215–288.
  * (55) Ellis, G.F.R.: Shear free solutions in General Relativity Theory. Gen. Relativ. Gravit. 43, 3253 (2011). [arXiv:1107.3669]
  * (56) Ellis, G.F.R., Bassett, B.A.C.C., Dunsby, P.K.S.: Lensing and caustic effects on cosmological distances. Class. Quantum Gravity 15, 2345 (1998). [arXiv:gr-qc/9801092]
  * (57) Ellis, G.F.R., Börner, G., Buchert, T., Ehlers, J., Hogan, C.J., Kirshner, R.P., Press, W.H., Raffelt, G., Thielemann, F.-K., van den Bergh, S.: What do we know about global properties of the Universe? In Dahlem Workshop Report ES19 The Evolution of the Universe, Berlin 1995, edited by G. Börner and S. Gottlöber, (Wiley, Chichester, 1997), 51.
  * (58) Ellis, G.F.R., Buchert, T.: The universe seen at different scales. Phys. Lett. A 347 (Einstein Special Issue), 38 (2005). [arXiv:gr-qc/0506106]
  * (59) Ellis, G.F.R., van Elst, H.: Cosmological Models (Cargèse Lectures 1998). In Proc. of the NATO Advanced Study Institute on Theoretical and Observational Cosmology, Cargèse, France, August 17–29, 1998, edited by M. Lachièze–Rey, (Kluwer Academic, Boston, 1999), 1–116 and NATO Science Series C 541, 1–116 (1999). [arXiv:gr-qc/9812046]
  * (60) Ellis, G.F.R., Maartens, R., MacCallum, M.A.H.: Relativistic Cosmology, (Cambridge University Press, Cambridge, 2012).
  * (61) Ellis, G.F.R., Nel, S.D., Maartens, R., Stoeger, W.R., Whitman, A.P.: Ideal observational cosmology. Phys. Rep. 124, 315 (1985).
  * (62) Ellis, G.F.R., Solomons, D.M.: Caustics of compensated spherical lens models. Class. Quantum Gravity 15, 2381 (1998). [arXiv:gr-qc/9802005]
  * (63) Elmardi, M.Y.A.: Topics in Relativistic Cosmology: Cosmology on the Past Lightcone and in Modified Gravitation. Ph.D. thesis, University of Cape Town, South Africa (January 2018). [Inspire]
  * (64) van Elst, H.: Extensions and applications of $1+3$ decomposition methods in general relativistic cosmological modelling. Ph.D. thesis, University of London (1996). [NASA/ADS]
  * (65) van Elst, H.: $(1+3)$-covariant methods in general relativistic cosmology, (University of Cape Town, 1998). [Online resource]
  * (66) van Elst, H.: An introduction to inductive statistical inference: from parameter estimation to decision-making. [arXiv:1808.10173]
  * (67) van Elst, H., Uggla, C., Lesame, W.M., Ellis, G.F.R., Maartens, R.: Integrability of irrotational silent cosmological models. Class. Quantum Gravity 14, 1151–1162 (1997). [arXiv:gr-qc/9611002]
  * (68) Enqvist, K., Mattsson, T.: The effect of inhomogeneous expansion on the supernova observations. J. Cosmol. Astropart. Phys. 02, 019 (2007). [arXiv:astro-ph/0609120]
  * (69) Etherington, I.M.H.: On the definition of distance in general relativity. Phil. Mag. J. Sci. 15, 761 (1933). Reprinted: Gen. Relativ. Gravit. 39, 1055 (2007).
  * (70) Fanizza, G., Gasperini, M., Marozzi, G, Veneziano, G.: Generalized covariant prescriptions for averaging cosmological observables. J. Cosmol. Astropart. Phys. 02, 017 (2020). [arXiv:1911.09469]
  * (71) Fleury, P., Dupuy, H., Uzan, J.P.: Interpretation of the Hubble diagram in a nonhomogeneous universe. Phys. Rev. D 87, 123526 (2013). [arXiv:1302.5308]
  * (72) Fleury, P., Dupuy, H., Uzan, J.P.: Can all cosmological observations be accurately interpreted with a unique geometry? Phys. Rev. Lett. 111, 091302 (2013). [arXiv:1304.7791]
  * (73) Fleury P., Nugier F., Fanizza G.: Geodesic-light-cone coordinates and the Bianchi I spacetime. J. Cosmol. Astropart. Phys. 06, 008 (2016). [arXiv:1602.04461]
  * (74) Gabbard, H., Williams, M., Hayes, F., Messenger, C.: Matching matched filtering with deep networks for gravitational-wave astronomy. Phys. Rev. Lett. 120, 141103 (2018). [arXiv:1712.06041]
  * (75) Gasperini, M., Marozzi, G., Veneziano, G.: A covariant and gauge invariant formulation of the cosmological ’backreaction’. J. Cosmol. Astropart. Phys. 02, 009 (2010). [arXiv:0912.3244]
  * (76) Gasperini, M., Marozzi, G., Nugier, F., Veneziano, G.: Light-cone averaging in cosmology: formalism and applications. J. Cosmol. Astropart. Phys. 07, 008 (2011). [arXiv:1104.1167]
  * (77) Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B.: Bayesian Data Analysis, 3rd Edition, (Chapman & Hall, Boca Raton, 2014). [Online resource]
  * (78) Giblin, J.T., Mertens, J.B., Starkman, G.D.: Departures from the Friedmann–Lemaître–Robertson–Walker cosmological model in an inhomogeneous Universe: a numerical examination. Phys. Rev. Lett. 116, 251301 (2016). [arXiv:1511.01105]
  * (79) Giblin, J.T., Mertens, J.B., Starkman, G.D.: Observable deviations from homogeneity in an inhomogeneous Universe. Astrophys. J. 833, 247 (2016). [arXiv:1608.04403]
  * (80) Giblin, J.T., Mertens, J.B., Starkman, G.D., Tian, C.: Limited accuracy of linearized gravity. Phys. Rev. D 99, 023527 (2019). [arXiv:1810.05203]
  * (81) Goldberg, J.N., Sachs, R.K.: A theorem on Petrov types. Acta Phys. Polon. 22, 1323. Suppl. Reprinted as a Golden Oldie: Gen. Relativ. Gravit. 41, 433 (2009).
  * (82) Gourgoulhon, E., Jaramillo, J.L.: A $3+1$ perspective on null hypersurfaces and isolated horizons. Phys. Rep. 423, 159-294 (2006). [gr-qc/0503113]
  * (83) Grasso, M., Korzyński, M., Serbenta, J.: Geometric optics in general relativity using bilocal operators. Phys. Rev. D 99, 064038 (2019). [arXiv:1811.10284]
  * (84) Grasso, M., Villa, E., Korzyński, M., Matarrese, S.: Isolating nonlinearities of light propagation in inhomogeneous cosmologies. Phys. Rev. D 104, 043508 (2021). [arXiv:2105.04552]
  * (85) Hawking, S.: Gravitational radiation in an expanding universe. J. Math. Phys. 9, 598 (1968).
  * (86) Hayward, S.A.: Quasilocal gravitational energy. Phys. Rev. D 49, 831 (1994). [gr-qc/9303030]
  * (87) Heinesen, A.: Multipole decomposition of the general luminosity distance ’Hubble law’ – a new framework for observational cosmology. J. Cosmol. Astropart. Phys. 05, 008 (2021). [arXiv:2010.06534]
  * (88) Heinesen, A.: Multipole decomposition of redshift drift: Model-independent mapping of the expansion history of the Universe. Phys. Rev. D 103, 023537 (2021). [arXiv:2011.10048]
  * (89) Heinesen, A.: Redshift drift as a model independent probe of dark energy. Phys. Rev. D 103, L081302 (2021). [arXiv:2102.03774]
  * (90) Heinesen, A.: Redshift drift cosmography for model-independent cosmological inference. Phys. Rev. D 104, 123527 (2021). [arXiv:2107.08674]
  * (91) Heinesen, A., Buchert, T.: Solving the curvature and Hubble parameter inconsistencies through structure formation-induced curvature. Class. Quantum Gravity 37, 164001 (2000). (Focus issue on the Hubble constant tension); Erratum (2020). Class. Quantum Gravity 37, 229601 (2000). [arXiv:2002.10831]
  * (92) Heinesen, A., Macpherson, H.J.: A prediction for anisotropies in the nearby Hubble flow. J. Cosmol. Astropart. Phys. 03, 057 (2022). [arXiv:2111.14423]
  * (93) Hellaby, C.: The mass of the cosmos. Mon. Not. R. Astron. Soc. 370, 239 (2006). [arXiv:astro-ph/0603637]
  * (94) Hildebrandt, H., et al.: KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing. Mon. Not. Roy. Astron. Soc. 465, 1454 (2017). [arXiv:1606.05338]
  * (95) Israel, W.: Covariant double-null dynamics. Helv. Phys. Acta 69, 3 (1996).
  * (96) Jimenez, P., Loeb, A.: Constraining cosmological parameters based on relative galaxy ages. Astrophys. J. 573, 37 (2002). [arXiv:astro-ph/0106145]
  * (97) Jordan, P., Ehlers, J., Sachs, R.K.: Beiträge zur Theorie der reinen Gravitationsstrahlung. Akad. Wiss. Lit. Mainz, Abhandl. Math.–Nat. Kl. 1, 1–62 (1961). English translation: Contributions to the theory of pure gravitational radiation. Gen. Relativ. Gravit. 45, 2691–2753 (2013).
  * (98) Kantowski, R.: The Lamé equation for distance-redshift in partially filled beam Friedmann–Lemaître–Robertson–Walker cosmology. Phys. Rev. D 68, 123516 (2003). [arXiv:astro-ph/0308419]