Ray tracing and Hubble diagrams in post-Newtonian cosmology

TL;DR

This paper investigates how small-scale inhomogeneities and opaque objects in the universe affect light propagation and cosmological measurements, revealing potential biases in estimating dark energy parameters.

Contribution

It introduces a method to compute redshifts and luminosity distances in inhomogeneous cosmologies with large density contrasts, highlighting the impact of opaque objects on observational biases.

Findings

Opaque objects cause a 10% shift in distance measures.

Photon trajectories avoid high-curvature regions due to opaque objects.

Selection effects can bias estimates of cosmological parameters.

Abstract

On small scales the observable Universe is highly inhomogeneous, with galaxies and clusters forming a complex web of voids and filaments. The optical properties of such configurations can be quite different from the perfectly smooth Friedmann-Lema\^itre-Robertson-Walker (FLRW) solutions that are frequently used in cosmology, and must be well understood if we are to make precise inferences about fundamental physics from cosmological observations. We investigate this problem by calculating redshifts and luminosity distances within a class of cosmological models that are constructed explicitly in order to allow for large density contrasts on small scales. Our study of optics is then achieved by propagating one hundred thousand null geodesics through such space-times, with matter arranged in either compact opaque objects or diffuse transparent haloes. We find that in the absence of opaque objects, the mean of our ray tracing results faithfully reproduces the expectations from FLRW cosmology. When opaque objects with sizes similar to those of galactic bulges are introduced, however, we find that the mean of distance measures can be shifted up from FLRW predictions by as much as $10\%$. This bias is due to the viable photon trajectories being restricted by the presence of the opaque objects, which means that they cannot probe the regions of space-time with the highest curvature. It corresponds to a positive bias of order $10\%$ in the estimation of $\Omega_{\Lambda}$ and highlights the important consequences that astronomical selection effects can have on cosmological observables.

Tables

Table 1. Table 1 : Classification of various approaches to inhomogeneous cosmological modelling. Top-down versus bottom-up refers to whether an FLRW background is assumed or not. Screening refers to whether inhomogeneities can affect the entire space-time, or only their own locale.

	Top-Down	Bottom-Up
Screened	Swiss cheese models	Lindquist-Wheeler lattice
Not Screened	perturbation theory	Bruneton-Larena & black hole lattices
	exact solutions	post-Newtonian cosmology

Equations

$$\mathbf{R}[\bar{\mathbf{g}}]\ne \bar{\mathbf{R}}[\mathbf{g}],$$

$$\epsilon \equiv \frac{v}{c}\ll 1,$$

$${\partial }_t\phi \sim \epsilon {\partial }_{\mu }\phi \ll {\partial }_{\mu }\phi .$$

$$\mathrm{d}{s}^2=\left[-1+h_{tt}^{(2)}+h_{tt}^{(4)}\right]\mathrm{d}{t}^2+2h_{t\mu }^{(3)}\mathrm{d}t\mathrm{d}{x}^{\mu }+\left[{\delta }_{\mu \nu }+h_{\mu \nu }^{(2)}\right]\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }+...,$$

$$\rho ={\rho }^{(2)}+{\rho }^{(4)}+...\mathrm{and}p=p^{(4)}+...,$$

$${\left(\frac{\dot{a}}{a}\right)}^2=\frac{8\pi G}{3}({\rho }_{\text{m}}+{\rho }_{\text{r}})-\frac{k}{a^2}+\frac{\Lambda }{3}+\mathcal{B},$$

$$\mathcal{B}\approx -{\left[4\pi G{\rho }_{\text{m}}L(t)\right]}^2\left(1.50-0.80\frac{{\Omega }_{\text{r}}}{{\Omega }_{\text{m}}}+1.76\frac{{\Omega }_k}{{\Omega }_{\text{m}}}\right).$$

$$k^a{k}_a=0$$

$$\frac{\mathrm{D}{k}^a}{\mathrm{d}\lambda }\equiv k^b{\nabla }_b{k}^a=\frac{\mathrm{d}{k}^a}{\mathrm{d}\lambda }+\Gamma \text{indices}{}_b^a{}_c{k}^b{k}^c=0,$$

$$\mathbf{k}=\omega (\mathbf{u}+\mathbf{d}),$$

$$1+z=\frac{\omega }{{\omega }_0}.$$

$$\frac{\mathrm{d}z}{\mathrm{d}\lambda }=-k^a{k}^b{\nabla }_a{u}_b\equiv -(1+z{)}^2{H}_{||},$$

$$\frac{{\mathrm{D}}^2{\xi }^a}{\mathrm{d}{\lambda }^2}=R\text{indices}{}_b^a{}_c{}_d{k}^b{k}^c{\xi }^d,$$

$$\frac{\mathrm{D}{s}_A^a}{\mathrm{d}\lambda }=\frac{k^a}{\omega }\frac{\mathrm{D}{u}_b}{\mathrm{d}\lambda }{s}_A^b.$$

$$\frac{{\mathrm{d}}^2{\xi }_A}{\mathrm{d}{\lambda }^2}={\mathcal{R}}_{AB}{\xi }_B,$$

$$\mathcal{R}=\left(\begin{array}{cc}\mathcal{R}& 0\\ 0& \mathcal{R}\end{array}\right)+\left(\begin{array}{cc}-{\mathcal{W}}_1& {\mathcal{W}}_2\\ {\mathcal{W}}_2& {\mathcal{W}}_1\end{array}\right).$$

$$\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\\ {\dot{\xi }}_1\\ {\dot{\xi }}_2\end{array}\right)(\lambda )=\mathcal{W}(\lambda \leftarrow {\lambda }_0)\left(\begin{array}{c}{\xi }_1\\ {\xi }_2\\ {\dot{\xi }}_1\\ {\dot{\xi }}_2\end{array}\right)({\lambda }_0),$$

$$\frac{\mathrm{d}\mathcal{W}}{\mathrm{d}\lambda }=\left(\begin{array}{cc}{0}_2& 1_2\\ \mathcal{R}& 0_2\end{array}\right)\mathcal{W},$$

$$\mathcal{W}({\lambda }_2\leftarrow {\lambda }_0)=\mathcal{W}({\lambda }_2\leftarrow {\lambda }_1)\mathcal{W}({\lambda }_1\leftarrow {\lambda }_0),$$

$${\xi }_A(\lambda )=-{\mathcal{D}}_{AB}{\theta }_B,$$

$$D_{\text{A}}\equiv \sqrt{\frac{A_S}{{\Omega }_0}}.$$

$$D_{\text{A}}=\sqrt{\det \mathcal{D}}.$$

$$D_{\text{L}}\equiv \sqrt{\frac{L}{4\pi I}},$$

$$D_{\text{L}}=(1+z)\sqrt{\frac{A_0}{{\Omega }_S}}=(1+z{)}^2{D}_{\text{A}}.$$

$$\mathrm{d}{s}^2=-(1-2\Phi )\mathrm{d}{t}^2+(1+2\Psi ){\delta }_{\mu \nu }\mathrm{d}{x}^{\mu }\mathrm{d}{x}^{\nu }+\mathcal{O}({\epsilon }^3).$$

$$M={\rho }_0{L}_0^3=\frac{3H_0^2{\Omega }_{\text{m}}{L}_0^3}{8\pi }\approx 4.5\times 10^{10}{M}_{\odot },$$

$$\Phi (t,\mathbf{x})$$

$${\Phi }_0(\mathbf{x})=\{\begin{array}{ll}-\frac{M}{2R^3}({\mathbf{x}}^2-3R^2)& \text{if }|\mathbf{x}|\le R\\ \frac{M}{|\mathbf{x}|}& \text{if }|\mathbf{x}|\ge R.\end{array}$$

$${\left(\frac{\dot{a}}{a}\right)}^2=\frac{8\pi G\rho }{3}+\frac{\Lambda }{3}+\mathcal{O}({\epsilon }^4)\mathrm{and}\frac{\ddot{a}}{a}=-\frac{4\pi G\rho }{3}+\frac{\Lambda }{3}+\mathcal{O}({\epsilon }^4),$$

$$\mathbf{u}=\left(1+\Phi +\frac{1}{2}{H}^2{\mathbf{x}}^2\right){\partial }_t+H{x}^{\mu }{\partial }_{\mu }+\mathcal{O}({\epsilon }^3),$$

Full text

1 Introduction

The Universe we see around us appears to be close to homogeneous and isotropic on large scales, but is very inhomogeneous and anisotropic on small scales. The standard approach to modelling this situation is to assume the existence of a Friedmann-Lemaître-Robertson-Walker (FLRW) geometry that obeys Einstein’s equations for some averaged matter density. Using this solution as a background, one can then perform a perturbative expansion of both the metric and the energy-momentum tensor, in order to incorporate inhomogeneous astrophysical structures and their gravitational fields. A possible drawback of using this approach is that it assumes that the evolution and averaging operations commute. This is not true in Einstein’s equations, where we have terms such as

[TABLE]

due to the non-linearity of the Ricci tensor. This non-commutativity means that the original solution, used to model the background, may not stay a faithful description of the Universe on large scales, even if it started out as such at early times. In other words, there is the possibility in Einstein’s theory that small-scale structures could influence the large-scale cosmic expansion; a possibility referred to as “back-reaction” in the cosmology literature [1, 2, 3, 4].

One way to side-step the difficult mathematical problems associated with back-reaction is to construct bottom-up cosmological models in which the large-scale properties of the Universe emerge, rather than being specified from the beginning. Such models do not require averaging, in the way that is necessary in the standard top-down approach, and instead allow the expansion of the Universe to emerge naturally from the boundary conditions between neighbouring astrophysical bodies. A class of such models that has recently been developed by two of us has been dubbed “post-Newtonian cosmology”, due to the use of the well-known post-Newtonian expansions in its construction [5, 6, 7]. These models are described in more detail in section 2.1 below, and will form the basis of the study of cosmological redshifts and distance measurements that we perform in this paper. They are of particular interest for the study of ray optics that we perform here as they naturally allow for the presence on non-linear structures in a fully self-consistent and unambiguous way.

While the kinematical back-reaction effects in these models have been precisely quantified in refs. [6, 7], the optical properties have not until now been explicitly calculated. These latter properties are of great practical importance, as they are the direct observables upon which almost all astronomical probes are based. Above and beyond any questions involving the large-scale expansion, the influence of structure on the optical properties of a space-time are of very significant interest for the determination of cosmological parameters. If non-linear structures have any systematic effect on the propagation of rays of light, then this could potentially have significant influence on any inferences of redshifts and distance measures over cosmological scales, and could consequently bias the estimation of (for example) the amount of dark energy in the Universe. The appropriate formalism for investigating the optical properties of the Universe is geometric optics, which is outlined in section 2.2. We will apply this formalism to our post-Newtonian cosmological models, using direct ray tracing techniques, in order to determine the influence of non-linear structure on observables such as the distance-redshift relation.

Of course, the optical properties of inhomogeneous cosmological models have been extensively studied in the past, and various frameworks have been developed to try to model the general behaviours that are expected from the effects of a lumpy matter content. Of particular relevance for the present study are the Einstein-Straus Swiss cheese models in which regions of FLRW geometry are excised and replaced by Schwarzschild [8, 9, 10], the Lindquist-Wheeler models that construct an approximate space-time out of Schwarzschild directly [11], and the Bruneton-Larena model that creates an approximate model that is valid for a short period of cosmic time [12]. The optical properties of these models have all been studied in the past, and we will refer to them further below. In particular, the post-Newtonian cosmological models can be considered an improved version of the Lindquist-Wheeler model as the approximation scheme is under much better control. It could also be considered an improvement on the Swiss cheese and Bruneton-Larena models, as it removes the need for a background FLRW geometry and is valid for much longer periods of cosmic time.

In section 2 we discuss the physics of our cosmological model, and the optical equations necessary to calculate redshifts and distance measures. In section 3 we then outline the method we use to calculate the optics along very many lines of sight. Section 4 then gives a detailed account of the results of our numerical integrations. This is followed in section 5 by a discussion of these results in the context of the theorems of Weinberg, Kibble & Lieu, the stochastic approach to lensing, and the recent results that have been obtained using methods from numerical relativity. Finally, in section 6 we conclude. We work with geometrised units throughout, where $G=c=1$. Bold symbols can refer to four-vectors, spatial vectors, or matrices. Latin indices $a,b,c,\ldots$ are space-time indices and run from 0 to 3. Greek indices $\mu,\nu,\rho,\ldots$ are spatial indices and run from 1 to 3. We use the signature $(-+++)$ for the metric, and the gravitational potentials $\Phi,\Psi$ are taken to be positive. Appendices A and B give details of the components of geometric quantities within our perturbative expansion, and the numerical tests we performed on our code in de Sitter space-time.

2 Fundamentals

This section introduces the key ingredients for the present work; namely the post-Newtonian cosmological models with non-linear matter fields developed in refs. [6, 7], and the laws of geometric optics in curved space-times as discussed in refs. [13, 14, 15, 16].

2.1 Post-Newtonian cosmology

The post-Newtonian cosmological models are built block-by-block, by patching together regions of post-Newtonian space-time [5]. The large-scale dynamics of the resulting universe then emerges naturally from the junction conditions at the boundaries between these regions. The purpose of such a model is twofold. Firstly, it allows us to build a cosmological solution constructively, from the bottom up. This means that the expansion of such a universe does not need to be assumed from the outset, as is the case in most approaches to cosmological modelling, and that the effects of inhomogeneity on the expansion dynamics can be studied in an unambiguous way, without assuming anything about averaging or coarse-graining. Secondly, it makes explicit the link between the gravitational fields of the objects that exist within the Universe and the expansion of the Universe itself. This allows us, for example, to construct self-consistent parametrized frameworks that can be used to test gravity in both cosmology and the weak-field limit [17]. It also allows us to build Hubble diagrams from the gravitational fields associated with non-linear astrophysical bodies themselves, rather than in any putative averaged “background”. It is this last possibility that we study here.

2.1.1 The post-Newtonian formalism

The post-Newtonian formalism [18] is a perturbative expansion scheme for the equations of general relativity, which is intended to deal with weak gravitational fields generated by slowly moving sources. The expansion itself is performed in terms of the smallness parameter

[TABLE]

where $v$ is the typical velocity of sources of the gravitational field and $c$ is the speed of light. This implies that, for consistency, the typical time derivative of any field $\varphi$ must be much smaller than its typical spatial derivative

[TABLE]

The space-time metric and the matter fields are then expanded using this rule, and the field equations are constructed and solved, order-by-order in $\varepsilon$. For the metric this means writing $g_{ab}=\eta_{ab}+h_{ab}$, where $\bm{\eta}$ denotes the Minkowski metric and $\bm{h}\ll\bm{\eta}$ is a perturbation, so that the line-element becomes

[TABLE]

where a superscript in brackets indicate the order of smallness of the quantity in terms of $\varepsilon$, and the ellipsis denotes terms of higher order than those given explicitly. For matter fields being modelled as a perfect fluid we end up with

[TABLE]

for the energy density and isotropic pressure, respectively. The orders of accuracy given above are sufficient to model the gravitational physics associated with massive particles to first post-Newtonian order. For the post-Newtonian treatment of light, however, we can neglect $h_{tt}^{(4)}$, $h_{t\mu}^{(3)}$, $\rho^{(4)}$ and $p^{(4)}$. The reader will note that not every order in $\epsilon$ is required in every variable, in order to get a complete description.

The post-Newtonian formalism, as outlined above, is designed for use in the presence of astrophysical bodies like stars and planets. It is, however, unsuitable for direct application to cosmology. This is because on distances comparable to the Hubble radius the coordinate velocity of the cosmological fluid becomes close to the speed of light, and the metric coefficients differ very significantly from the assumed Minkowski background. The perturbative expansion therefore breaks down, and the entire formalism fails. This problem motivates the piecewise construction summarised below, which allows the post-Newtonian approach to be extended to cosmology in a fully self-consistent fashion.

2.1.2 A bottom-up cosmological model

Consider a region of space-time much, much smaller than the ball bounded by the Hubble sphere. At all points within this region, except in the immediate vicinity of compact objects like neutron stars and black holes, the post-Newtonian formalism should be expected to hold. Now consider a second comparable region, adjacent to the first one. A similar expansion can be performed there, except that the background will not be identical to that of the first region; they will both be regions of Minkowski space, but they will join at an angle to allow for curvature of space-time on cosmological scales (see the left-hand panel of figure 1). Thus the whole cosmological space-time can be obtained step-by-step, performing a post-Newtonian expansion in each small region independently, and then by gluing together resulting cells.

Any healthy junction between two space-time regions must fulfil the Darmois-Israel conditions [19, 20]. These state that, in the absence of surface layers, the intrinsic metric and the extrinsic curvature of the hypersurface that separates any two regions of space-time must be the same whether it is calculated on one side or the other. In order to explicitly ensure these junction conditions, refs. [6, 7] assumed that the Universe could be modelled as being periodic, and that neighbouring cells are symmetric with respect to the boundary that separates them (see the right-hand panel of figure 1). In this case, the whole lattice can be generated from a single cell by successive reflections.

Cosmic expansion naturally emerges from the boundary conditions between neighbouring regions, in this model. The large-scale cosmological dynamics can be studied by considering the length of a cell edge, as measured in the rest space of a set of observers that follows the trajectory of the boundary. The analogue of the cosmic scale factor is then given by the ratio $a(t)\equiv L(t)/L_{0}$ between the size of the cell edge at time $t$ and its value today. After extending the post-Newtonian formalism to include radiation and $\Lambda$, the equation that governs the evolution of this scale factor can be shown to be given by [7]

[TABLE]

where $\rho_{\text{m}}$ is the energy density in pressureless matter, $\rho_{\text{r}}$ is the energy density in radiation, $k$ is spatial curvature, and a dot denotes a derivative with respect to proper time measured by the comoving observers. Equation (2.5) is essentially the standard Friedmann equation, but with an additional source term $\mathcal{B}$ on the right-hand side that is of order $\varepsilon^{4}$. This new term encodes the back-reaction effect of structure within the cell on the large-scale expansion, and for a point-like mass at the centre of the cell takes the simple form

[TABLE]

This contribution to the large-scale expansion is purely relativistic, and in a universe dominated by dust and $\Lambda$ can be seen to be principally determined by a term that behaves precisely like an extra radiation component. Such a correction is reminiscent of the outcome from Green and Wald’s approach to backreaction [22, 23], and was also found in earlier work on reflection-symmetric subspaces [24]. The $\Omega_{\text{m}},\Omega_{\text{r}}$, and $\Omega_{k}$ parameters in eq. (2.6) denote the cosmological parameters associated with matter, radiation, and curvature respectively.

It is within the context of these lattice models that we will investigate the propagation of light, and the construction of Hubble diagrams. This will be done using opaque compact objects at the centre of each cell, as well as with dispersed transparent matter. We intend these configurations to roughly correspond to the type of matter that exists at the centre of galaxies, and the dark matter haloes within which galaxies are known to exist. The post-Newtonian approach to modelling gravitational fields in these constructions is ideal for addressing old questions about the effect of discretization of matter on light propagation in cosmology [25], as the matter fields are allowed to be highly non-linear and the intervening regions of vacuum are allowed to be entirely Ricci flat.

2.1.3 Comparison with other approaches

The post-Newtonian cosmological models described above can be contrasted with other approaches to inhomogeneous cosmology. Two of the most well-known of these are standard cosmological perturbation theory on an FLRW background [26], and the Swiss cheese models originally constructed by Einstein and Straus [8, 9]. The former of these approaches is of course extremely versatile, but is strictly only valid in the regime where density contrasts are small (of the same order as the expansion parameter). The latter allows for arbitrarily large density contrasts by modelling inhomogeneity as patches of either Schwarzschild [8, 9, 10, 27, 28, 29, 30, 31, 32], Lemaître-Tolman-Bondi [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45] or Szekeres solutions [46, 47, 48, 49, 50], but is somewhat less versatile due to its reliance on matter being well-modelled by these exact solutions. A common feature of both of these models, however, is that they can be considered to be top-down constructions, as they assume a global FLRW background geometry from the outset. This means that the large-scale cosmological behaviour is pre-specified, and cannot be affected by the formation of structure.

In comparison, post-Newtonian cosmological models do not contain any global assumptions at the level of the metric, but do assume that the matter distribution can be modelled as periodic (and hence statistically homogeneous). From this point of view they can be compared with a number of other recent articles on lattice cosmology. The first of these was based on the Lindquist-Wheeler construction [51, 52, 53], where Schwarzschild regions are glued together using approximate matching conditions in order to form a lattice of point masses [11, 54]. A second approach is the Bruneton-Larena lattice Universe [12, 55], which is based on perturbatively solving the Einstein field equations after performing a Fourier decomposition, and which results in a cosmological model that can be used for short periods of cosmic time. Further attempts at modelling the Universe as a lattice of point-particles were performed using Regge calculus by Liu and Williams [56, 57, 58], and by using geometrostatics and numerical relativity methods by a number of other authors [59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70]. The post-Newtonian method has the advantage of extra versatility, when compared to these other approaches.

As well as the dichotomy between top-down and bottom-up approaches, we could also classify inhomogeneous cosmological models according to whether or not inhomogeneities are ‘screened’ (in the sense that what happens in a given region of a model depends only on its immediate neighbourhood, or if it is affected by all matter that exists in the Universe). By construction, Swiss cheese models and Lindquist-Wheeler lattices belong to the first category, as within a Swiss cheese hole or a Lindquist-Wheeler cell space-time is entirely determined by the closest mass. By contrast, most other approaches result in individual inhomogeneities affecting all of the rest of space-time. The comparison of these approaches to inhomogeneous cosmology are summarised in Table 1. The reader may note that post-Newtonian cosmology combines the advantages of being background-free (bottom-up), as well as being realistic (not screened) and versatile.

2.2 Light propagation in curved space-time

The purpose of this article is to determine the Hubble diagrams that would be constructed by observers in the post-Newtonian cosmological models, described above. This requires us to understand how light propagates in curved space-times, so that we can determine both distance measures and redshifts. In this section we therefore provide a brief overview of the main results of geometric optics in the presence of gravitation. This will serve to introduce the reader to the subject, as well as to specify the conventions and notation we will use in the rest of the article. For more detailed descriptions see refs. [13, 14, 15, 16].

2.2.1 Rays of light

Light is an electromagnetic wave, and therefore propagates according to Maxwell’s equations. Assuming a minimal coupling between the electromagnetic and gravitational fields, and under the eikonal approximation of geometric optics, it can be shown that rays of light must follow null geodesics. The four-vector $\bm{k}$ tangent to each ray must therefore satisfy

[TABLE]

and

[TABLE]

with $k^{a}\equiv\mathrm{d}x^{a}/\mathrm{d}\lambda$ and ${\mathrm{d}k^{a}}/{\mathrm{d}\lambda}\equiv k^{b}k\indices{{}^{a}_{,b}}$, and where $\lambda$ denotes an affine parameter along the ray of light. If the rays of light are non-rotating then the four-vector components $k_{a}$ can also be written as the gradient of the wave’s phase, such that $k_{a}=\partial_{a}\phi$.

For an observer following an integral curve of the time-like four-velocity $\bm{u}$, the wave four-vector $\bm{k}$ can be split into a temporal part $\omega=-u_{a}k^{a}$ and a spatial part $e_{a}=(g_{ab}+u_{a}u_{b})k^{b}$. These quantities encapsulate the cyclic frequency and the direction of propagation, respectively. If we now write $\bm{e}=\omega\bm{d}$, then

[TABLE]

from whence it can be seen that $u_{a}d^{a}=0$ and $d^{a}d_{a}=1$. In writing eq. (2.9) we have chosen $\bm{k}$ to be future oriented, so that $\lambda$ increases to the future.

The redshift involved in most cosmological observations is defined as the fractional difference between the frequency $\omega$ emitted by a light source and the frequency $\omega_{0}$ at which it is observed, such that

[TABLE]

Now suppose that a collection of sources lie along the light ray, forming a field of four-velocities $\bm{u}(\lambda)$. In this case the evolution of $z$ with $\lambda$ can be shown to read

[TABLE]

where $H_{||}$ is the local rate of expansion of the family of sources in the direction of propagation of the ray of light, and where we have chosen units so that $\omega_{0}=1$. The function $z=z(\lambda)$ therefore contains information about the expansion of the Universe.

2.2.2 Narrow beams of light

To calculate other astronomical observables, such as measures of distance, requires a mathematical description of beams of light. For these purposes we will consider narrow beams of light, which can be modelled as bundles of null geodesics that are infinitesimally separated from each other. The distance between two neighbouring photons within the beam is then described by the separation vector $\bm{\xi}$, which is defined as being geodesic and orthogonal to $\bm{k}$. We therefore have $k_{a}\xi^{a}=0$ and

[TABLE]

where $R\indices{{}^{a}_{b}{}_{c}{}_{d}}$ are the components of the Riemann curvature tensor. This last expression is the geodesic deviation equation, which (after projection and integration) allows us to determine the angular diameter distance and luminosity distance along the beam.

In order to study the morphology of the beam it is convenient to project it onto a screen that is orthogonal to its direction of propagation. A set of orthonormal basis vectors that span such a screen can be written as $(\bm{s}_{1},\bm{s}_{2})=(\bm{s}_{A})_{A=1,2}$, where $u_{a}s_{A}^{a}=d_{a}s^{a}_{A}=0$. The basis vectors on screens at different positions along the beam can then be related by imposing the partial parallel-transportation condition $(\delta^{a}_{b}+u^{a}u_{b}-d^{a}d_{b})\mathrm{D}s_{A}^{b}/\mathrm{d}\lambda=0$, which for computational purposes can be more conveniently written as

[TABLE]

Transporting the basis vectors of the screen-space in this way prevents the beam’s morphology from spuriously rotating as one proceeds along its direction of propagation, and results in what is usually referred to as the Sachs basis.

The components of the separation vector in the Sachs basis can now be written as $\xi_{A}\equiv\xi_{a}s^{a}_{A}$. These new objects exist entirely within the plane of the screen, and give us direct information about the separation between photons within the beam (i.e. the separation that an observer would measure if he or she replaced the screen with a photographic plate). Using eq. (2.12), the projected separation vectors can be seen to obey the following evolution equation:

[TABLE]

where $\mathcal{R}_{AB}\equiv R_{abcd}s^{a}_{A}k^{b}k^{c}s_{B}^{d}$ is referred to as the optical tidal matrix, and the whole equation is known as the vector Sachs equation. The reader may note that the position of the indices $A,B,\ldots$ does not matter, as they are raised and lowered by $\delta_{AB}$. The optical tidal matrix can be decomposed into a scalar part and a trace-free part as

[TABLE]

The scalar part depends only on the Ricci curvature $R_{ab}$ of the space-time, with $\mathscr{R}=-(1/2)R_{ab}k^{a}k^{b}$, and causes the beam of light to be focussed without distortion. Meanwhile, the trace-free part depends only on the Weyl curvature $C_{abcd}$ of the space-time, with $\mathscr{W}_{1}\equiv C_{abcd}s_{1}^{a}k^{b}k^{c}s_{1}^{d}=C_{abcd}s_{2}^{a}k^{b}k^{c}s_{2}^{d}$, $\mathscr{W}_{2}\equiv C_{abcd}s_{1}^{a}k^{b}k^{c}s_{2}^{d}$, and causes distortion of the beam (by tidal gravitational forces).

Because the vector Sachs equation (2.14) is linear in $\xi_{A}$, any solution must be linearly related to its initial conditions at $\lambda_{0}$. More precisely, there must exist a $4\times 4$ Wronski matrix $\bm{\mathcal{W}}$ such that

[TABLE]

where an over-dot denotes $\mathrm{d}/\mathrm{d}\lambda$, and where $\bm{\mathcal{W}}(\lambda_{0}\leftarrow\lambda_{0})=\bm{1}_{4}$. From eq. (2.14), we can then write

[TABLE]

where $\bm{0}_{n}$ and $\bm{1}_{n}$ denote the $n\times n$ zero and unity matrices, respectively. This reduces the problem of finding $\xi_{A}$ along the beam to solving a set of first-order ordinary differential equations. It also allows us to use the following matrix multiplication rule:

[TABLE]

which is extremely advantageous when considering a cosmological model that is constructed in a piecewise fashion, as we are doing here. We will therefore use eq. (2.17) as the final form of the geodesic deviation equation when we numerically solve these equations in section 4.

Although we will integrate eq. (2.17) to find $\bm{\mathcal{W}}$ at all values of $\lambda$, the most interesting parts of this matrix are the $2\times 2$ cells in the top-right corner. These are collectively referred to as the Jacobi matrix, $\bm{\mathcal{D}}$, and are sufficient to determine the angular diameter distances of astrophysical objects that lie within the beam, as we will now discuss.

2.2.3 Distance measures

The first step in calculating the angular diameter distance of a source is recognising that the initial condition at the point of observation should be taken to be $\xi_{A}(\lambda_{0})=0$, as this is the point at which the beam converges. The projection of the separation vector $\xi_{A}(\lambda)$ is then linearly related to $\dot{\xi}_{A}(\lambda_{0})$ only, and must therefore be given by the Jacobi matrix $\bm{\mathcal{D}}=(\mathcal{D}_{AB})_{A,B=1,2}$. This gives

[TABLE]

where $\theta_{A}=-\dot{\xi}_{A}(\lambda_{0})$ is the angular separation between the two light rays separated by $\bm{\xi}$, as measured at the point of observation. In other words, the matrix $-\bm{\mathcal{D}}$ is the map from observed angular separations to spatial separations in the screen space at $\lambda$. This is exactly what is required to define the angular diameter distance.

The definition of the angular diameter distance to a source with cross-sectional area $A_{S}$, and that subtends the angle $\Omega_{0}$ on the observer’s sky, is

[TABLE]

This quantity is defined in analogy to the way that one would infer distance in a flat space-time, if the same source subtended the same angle. Using eq. (2.19), and recognising that the area of a parallelogram is given by the determinant of the matrix formed from the vectors that define it, one can directly deduce that

[TABLE]

Once the Wronskian matrix $\bm{\mathcal{W}}$ has been obtained, the angular diameter distance at all points along the beam can therefore be readily deduced by simple algebraic operations. This gives $D_{\text{A}}$ as a function of $\lambda$, and can be used to obtain $D_{\text{A}}$ as a function of redshift by using the solution of eq. (2.11).

From knowledge of the angular diameter distance, it is relatively straightforward to obtain an expression for the luminosity distance. This latter measure is defined as

[TABLE]

which would be the distance that one would infer in a flat space for a sources that has luminosity $L$, and is measured to have intensity $I$. In this case the beam of light must be taken to be focussed at the emitting source, such that $\xi_{A}(\lambda_{S})=0$. If the cross-sectional area of this beam is $A_{0}$ at the observer, and the angle it subtends at the source is $\Omega_{S}$, then by photon conservation the luminosity distance is simply given by

[TABLE]

The last of these equalities comes from Etherington’s reciprocity theorem, which states that $A_{0}\Omega_{0}=(1+z)^{2}A_{S}\Omega_{S}$ in any space-time [71]. The remarkable nature of this theorem means that $D_{\text{L}}(z)$ can now be obtained without integrating the light beam forward from the source, if one already has knowledge of the beam that is focussed at the observer.

3 Method

Let us now turn to the detailed implementation of light propagation in post-Newtonian cosmology using ray tracing techniques. In order to capture the leading-order post-Newtonian effects on the null geodesics that constitute the paths followed by individual rays of light we require the metric to be specified to order $\varepsilon^{2}$:

[TABLE]

The scalar gravitational potentials $\Phi$ and $\Psi$ are both objects of order $\varepsilon^{2}$ in the post-Newtonian expansion of the metric, and it can be noted that at this order there can exist no vector or tensor perturbations. The required geometric quantities associated with this metric are given in appendix A, up to order $\varepsilon^{2}$. While the form of these scalar potentials is similar to those used in standard cosmological perturbation theory, their precise meaning is different. In particular, the functions $\Phi$ and $\Psi$ in eq. (3.1) contain information about the global expansion, as well as the gravitational fields of nearby objects.

3.1 A universe in a cell

As the space-time we are considering is periodic, we only need to describe the properties of a single cell in order to get the geometry of the entire universe. The tessellation of space that we wish to consider here is based on a cubic cell, which we use to tile a flat three-dimensional reference space. In this case, the volume of the lattice cell can be taken to be $V_{0}\simeq L_{0}^{3}=1\>\mathrm{Mpc}^{3}$ today, where $L_{0}$ is the proper length of one of its edges at the present time. The condition that the lattice is constructed in a flat space means that the total rest mass within each cell must be exactly specified by [6, 7]

[TABLE]

where $H_{0}$ is the current value of the Hubble rate, $\rho_{0}$ the current mean density of matter (both dark and baryonic) in the Universe, and $\Omega_{\text{m}}=0.3$ its ratio with the critical density. In what follows we will assume this mass to be contained within a spherically symmetric static body at the centre of the cell, with radius $R\ll L_{0}$. This setup is illustrated in figure 2. One may note that other tessellations exist, in three-dimensional spaces of positive and negative curvature.

We intend the matter within each cell to represent either a galaxy, or a galactic dark matter halo. In the rest of this article we will therefore consider two different simulations:

(i) Galaxy simulations:

In this case the massive body at the centre of the cell models the compact bulge of a spiral galaxy. We take this object to have radius $R=3\>\mathrm{kpc}$, and assume it to be opaque. This last condition is considered to be largely true in the real Universe, where the density of luminous matter and gas in a bulge makes it almost impossible to view any objects that lie behind it. In practice, opacity is ensured in our models by artificially removing from the simulation any rays that enters into the region $r<3\>\mathrm{kpc}$. This setup is comparable to the Einstein-Straus Swiss cheese model [8, 9], and the black-hole lattice considered in ref. [72].

(ii) Halo simulations:

Here the central body represents a diffuse dark matter halo. We take this object to have radius $R=30\>\mathrm{kpc}$, and to be completely transparent. We assume for simplicity that the halo has uniform density $\rho_{\text{H}}=3M/4\pi R^{3}$. This is far from being fully realistic, but is expected to be sufficient to address the questions we pose in this article. More realistic density profiles will be implemented in the future. This setup is comparable to those that are often considered in the context of standard cosmological perturbation theory, as well as the LTB and Szekeres Swiss cheese models.

In both cases (i) and (ii) the gravitational potentials can be shown to be given by [6, 7]

[TABLE]

and $\Psi(t,\bm{x})=\Phi(t,\bm{x})-{\Lambda}/{4}$, where the origin of the spatial coordinate system lies at the centre of the cell, and where $\Phi_{0}(\bm{x})$ is given by

[TABLE]

Note that the case $r\leq R$ is relevant for the halo simulations only, as in the galaxy simulations light is not allowed to enter this region.

To the same order of accuracy, it can be shown that the global expansion in both cases (i) and (ii) is given by the following emergent Friedmann equations:

[TABLE]

where $\rho=M/L^{3}$ is the averaged energy density of matter in each cell (to order $\varepsilon^{2}$), and where $\Lambda$ is the cosmological constant. At this level of accuracy, these two equations are identical to the standard Friedmann equations of homogeneous and isotropic cosmological models. They are, however, derived from the Israel junction conditions applied at the boundaries between neighbouring cells, rather than by assuming that averaged energy densities and geometries can be substituted into Einstein’s equations directly. This means that they are valid for arbitrary distributions of matter within each cell, including the types of highly non-linear density contrasts required to describe galaxies and clusters of galaxies. Once these equations have been solved, the length of the cell edge is given by $L(t)=a(t)L_{0}$, where $a(t_{0})=1$.

The final piece of information required to compute observables such as redshift $z$ and angular diameter distance $D_{\text{A}}$ is a four-velocity field $\bm{u}$ that covers the whole of the space-time (or, equivalently, the whole of a cell). Once specified, $\bm{u}$ can be taken to define the rest frame of a hypothetical light source at any point in the space-time. We require observers who follow the integral curves of this field to be comoving with the boundary when they are coincident with it, and to be at rest at the centre of the cell if they are located at that position. In other words, we choose our hypothetical sources on the boundary to be comoving with the boundary. A vector field that obeys these conditions, and that has been continued into the interior of the cell, is given by [6]

[TABLE]

where $H\equiv\dot{a}/a$. The observers who follow the integral curves of this field are analogous to the fundamental (comoving) observers used in standard FLRW cosmology. We use them to define the frequency of light, the Sachs basis, and every other frame-dependent variable.

3.2 Initial conditions

The ray tracing procedure outlined in section 2.2 requires appropriate initial conditions to be set. Specifically, we need to specify the position of the observer in both space and time, as well as the direction on the observer’s sky in which the beam of light will be propagated. Once this information has been given, the tangent vector to the light ray, the Sachs basis vectors, and the Wronski matrix will be found by integrating eqs. (2.8), (2.13), and (2.17) backwards in time, from the observation event $O$ to its source at affine distance $\lambda$.

We choose the coordinates of $O$ by first performing a coordinate transformation, as follows:

[TABLE]

where hatted coordinates $\hat{x}^{a}$ are the analogue of comoving synchronous coordinates in the post-Newtonian cosmological framework [6], and where $\hat{r}_{0}\equiv\sqrt{\hat{x}_{0}^{2}+\hat{y}_{0}^{2}+\hat{z}_{0}^{2}}$. We want our observer to remain at fixed position with respect to these coordinates, so they are as similar as possible to the comoving observers used in FLRW cosmology. We then make the further choice that the time of observation is at

[TABLE]

which is the time at which an observer in this model exists in order to measure $H=H_{0}$, for any given values of $\Omega_{\Lambda}$ and $\Omega_{\text{m}}$. The remaining comoving spatial coordinates are chosen to be $\{\hat{x}_{0},\hat{y}_{0},\hat{z}_{0}\}=\{-0.4L_{0},0.1L_{0},0\}$. This places the observer in the bulk of the cell, far from the central mass, the cell edge, and all axes of discrete rotational symmetry [24].

If the observer at $O$ is comoving, then his or her four-velocity $\bm{u}_{0}$ is given by eq. (3.6) at $O$. The rest space of such an observer is then spanned by a triad $(\bm{e}_{\alpha})_{\alpha=1,2,3}$ with components

[TABLE]

so that $(\bm{u}_{0},\bm{e}_{\alpha})$ forms an orthonormal basis at $O$. Our choice of units of time is now such that the observed frequency of light at $O$ is given by $\omega_{0}\equiv-(u^{a}k_{a})_{0}=1$, and our specification of a spatial triad means the direction of incoming photons on the observer’s celestial sphere be written as

[TABLE]

where $\theta$ and $\phi$ are standard spherical coordinates. Once the pair of coordinates $(\theta,\varphi)$ have been chosen, the initial conditions for the four-vector tangent to the rays of light can then be seen to become

[TABLE]

For each light ray we randomly pick an observation direction, given by $(\theta,\varphi)$, in such a way that the observer’s celestial sphere is homogeneously covered; the associated probability density function thus reads $p(\theta,\varphi)={\sin\theta}/{4\pi}$ if $\theta\in[0,\pi]$ and $\varphi\in[0,2\pi)$, and zero otherwise. This fully specifies all of the initial conditions for our ray tracing experiment.

3.3 Reflection of light at the cell’s boundary

The periodic nature of our lattice means that instead of propagating light rays between neighbouring cells, which is the physical situation we wish to investigate, we can simply reflect our light rays off the cell boundaries and continue propagating it within the same cell. This will produce exactly the same result as propagating the light between cells, as we have reflection symmetry about each of our boundaries (see figure 3). It is also a more economical way of modelling an infinitely extended universe.

The reflection relations required at our cell boundaries can be derived using the fact that all optical quantities (wave four-vector, Sachs basis, Jacobi matrix, etc.) must be continuous. This fact is ensured to be true at our cell boundaries due to the satisfaction of the Israel junction conditions between neighbouring cells, and the fact that $\bm{u}$ has been taken to be comoving with the boundary. Thus, for any four-vector $\bm{v}$ attached to the light beam, its reflected counterpart $\bm{v}^{\prime}$ must read

[TABLE]

where $\bm{n}$ is the outward-pointing unit normal four-vector to the boundary, which for boundaries at coordinate positions $x^{\alpha}=\pm L(t)/2$ is given by

[TABLE]

The reflection equation (3.14) applies in particular to the wave four-vector $\bm{k}$, Sachs basis vectors $\bm{s}_{A}$, and separation four-vector $\bm{\xi}$. The reader may note, however, that by construction we have $\bm{u}^{\prime}=\bm{u}$, and hence $\omega^{\prime}=\omega$. A further consequence of eq. (3.14), applied to $\bm{s}_{A}$ and  $\bm{\xi}$ is that the screen components of $\bm{\xi}$ are unchanged by the reflection:

[TABLE]

This means the Jacobi matrix, the Wronski matrix, and all distance measures must be continuous at reflection, as required. It also means we can apply the Wronski multiplication rule (2.18) as if there were no reflection at all.

3.4 The ray-tracing code

To implement the integration of the optical equations, as outlined in section 2.2, we developed a ray-tracing code in C. Initial conditions for the numerical integrations we perform are set as described in section 3.2, and the reflection operations at cell boundaries are performed as outlined in section 3.3. The pipeline that implements this set of operations is sketched in figure 4.

The code is iterated to look-back time $t_{0}-t$ (in the coordinates of eq. (3.1)), with an evolving step $\Delta t_{i}\equiv t_{i}-t_{i+1}=a(t_{i})1.67\times 10^{4}\>\mathrm{yr}$. This choice is made so that the ratio between the time step and the cell size $L(t)$ remains essentially constant, at $\Delta t/L(t)=5\times 10^{-3}$. This ensures that the accuracy of the code is stable over time. At each time step, we evaluate the gravitational potential $\Phi$ and its derivatives using eq. (3.3). The sum over $\bm{p}$ in this equation is truncated to $\bm{p}\in\{-5,5\}^{3}\setminus\bm{0}$, so that only the $11^{3}=1,331$ closest masses are taken into account. The error that is made on the value of the potential due to this truncation is at most

[TABLE]

which essentially occurs in the region of space closest to the boundaries of the cell (where the potential is smallest). We expect this to have a negligible impact on our final results. The maximum value of the PN expansion parameter $\varepsilon$ is evaluated at the edge of the cell to be $\varepsilon\sim H_{0}L_{0}\sim 2\times 10^{-4}$; in terms of the gravitational potential, it is maximum at the surface of the mass clumps, with $\varepsilon\sim\sqrt{GM/R}\sim 10^{-3}$ in the galaxy simulation.

Meanwhile, the differential equations of light propagation (2.8), (2.13), and (2.17) are solved using a simple Euler integration with the above time step. The global numerical error on the redshift along an individual line of sight is estimated to be $N\times(\Delta t/L_{\text{H}})^{2}$, where $N$ is the number of steps, and $L_{\text{H}}$ is the Hubble radius (which is the typical distance over which $\bm{u}$ varies appreciably). As the integration is performed over cosmological distances we have $N\sim L_{\text{H}}/\Delta t$, so that

[TABLE]

We have tested our code against known exact expressions in de Sitter space-time, with the results detailed in appendix B. We find that the code exhibits an accuracy for $z(\lambda)$ in agreement with the above prediction. Estimating the global error on the Euler integration of the Jacobi matrix is more subtle. On the one hand, the local error is larger because the optical tidal matrix varies over distances comparable to $L\ll L_{\text{H}}$. On the other hand, its sign flips during the propagation through a cell—curvature increases while the photon is approaching the central clump, and then decreases while it moves away—so that local errors do not combine cumulatively. A conservative estimate can be obtained by considering a random local error with standard deviation $(\Delta t/L)^{2}$, which along an individual line of sight yields

[TABLE]

Numerical tests of the convergence of our numerical integrations indicate a slightly smaller error, below one part in a thousand. This error cannot be tested with the de Sitter case presented in appendix B, however, as the optical tidal matrix is exactly zero in that case.

We will propagate a given light beam within our lattice cell until it reaches a boundary, at which point it will be reflected according to the rules described in section 3.3. Due to the discrete nature of the time steps in the numerical integration, the actual intersection between the photon’s path and the world-sheet of the boundary generically occurs between two steps $i$ and $i+1$, so that $x^{\mu}_{i+1}$ lies outside of the cell. We correct for this by determining the actual point of intersection, and propagating the ray of light back to this point. Observables such as the redshift $z$, the Jacobi matrix $\bm{\mathcal{D}}$, and the associated angular distance $D_{\text{A}}$ are computed at each time step, and saved in an output file at every $z=0.1n$, with $n\in\mathbb{N}$.

Finally, the code stops along each ray of light if any one of three cases occurs:

The maximum redshift $z_{\text{max}}=1.5$ is reached; 2. 2.

A caustic occurs, at which point $\det\bm{\mathcal{D}}\leq 0$, and our approach breaks down; 3. 3.

The ray of light enters an opaque clump of matter (in the galaxy simulations only).

4 Results

Let us now present the results produced by our ray-tracing code in post-Newtonian cosmologies, before going on to compare these results with existing analytical and numerical approaches in section 5. Our code is constructed so that it can build up statistics associated with physical observables, such as redshift and angular diameter distance, by integrating along large numbers of individual lines-of-sight. Two examples of these individual paths are shown, for illustrative purposes, in figure 5. The left panel of this figure shows the deflection of light that occurs close to a compact object, while the right panel shows the first 20 reflections of a beam of light that stays far from the central mass. In both of these diagrams time flows from light to dark colours, while the arrows indicate the direction of numerical integration. Production of the first of these images required us to dramatically reduce the time step of the numerical integrations, in order to resolve the gravitational potential in the vicinity of the compact object111Note that any light ray passing this close to a mass point in our galaxy simulations would be excluded, due to our selection rules..

Let us now address the central question of the statistics of observables that are calculated from considering many such beams of light, and how these relate to the predictions of the homogeneous and isotropic FLRW models. We will do this by first calculating the redshift as a function of affine distance in section 4.1, before moving on to consider the angular distance-affine parameter relation in section 4.2. With both of these data sets in hand, we will then combine them to construct Hubble diagrams in section 4.3. This last quantity is the direct observable, and is what astronomers within our space-time should be expected to measure. In practice, for each simulation (galaxy and halo) we shoot $N_{\rm L}=10^{5}$ light beam in random directions on the observer’s celestial sphere. The statistical average of an observable $Q$ is then defined by

[TABLE]

which can be seen to correspond to a directional average, as the random set of beams evenly covers the sky.

4.1 Redshift

Let us start by considering the redshift along our rays of light. Figure 6 compares the fractional difference between the $z(\lambda)$ relation obtained in our simulation, using the direct methods outlined in section 2.2, with the FLRW relation $\bar{z}(\lambda)$ obtained by integrating the well-known expression

[TABLE]

We display results for both the galaxy simulation (with a compact opaque core), and the halo simulation (with a transparent and diffuse distribution of mass). It can be see that the mean deviation from the predictions of the corresponding FLRW models is less than one part in $10^{5}$, and decreases as the integration continues. While the numerical values obtained for $\left\langle z\right\rangle_{\Omega}/\bar{z}-1$ are larger than the error bars displayed, which account for the statistical deviation of numerical results obtained along $10^{5}$ different geodesics, they are not significantly larger than either the numerical error estimates obtained for our routine (see appendix B) nor the errors that should be expected from neglecting post-post-Newtonian gravitational potentials [6, 7]. Both of these latter sources of error are conservatively estimated to be at the level of about one part in $10^{5}$.

The results displayed in figure 6 lead us to conclude that the deviations of $z(\lambda)$ from their values in a corresponding FLRW universe are no greater than one part in $10^{5}$. This is an intrinsic limit to the accuracy that one could ever hope to obtain from a perturbative approximation of the gravitational field, as used in this study, and in this sense we have saturated the bound on deviations in redshift that can be induced from post-Newtonian gravity. To quantify the effect of inhomogeneity on cosmological redshift one must therefore go to post-post-Newtonian order (i.e. $\mathcal{O}(\epsilon^{4})$), and increase the precision of the numerical integrations to a comparable level of accuracy (around one part in $10^{10}$). This is beyond the scope of the present study.

These results should not, of course, be confused with the effects of peculiar motions on the $z=z(\lambda)$ relation. These motions can have a very considerable impact on redshifts [73, 74], as Doppler effects can contribute significantly to physical observables such as redshift-space distortions. By choosing a set of observers who are as close to comoving as possible, as specified by the four-vector field in equation (3.6), we have explicitly neglected such contributions. In a fully realistic model these effects can, and should, be added to the background cosmological redshift. This can be done in a relatively straightforward way by simply boosting either the observer or the source, but the effects of this are well studied, and would work in exactly the same way here as they do in the standard approach to cosmological modelling. We will therefore not consider them any further.

The negligible effect of inhomogeneity on background redshifts that we have found here can be compared to similar results in the Einstein-Straus Swiss cheese model [32], where deviations from $\bar{z}(\lambda)$ are due to the Rees-Sciama effect [75]. These models similarly neglect the contribution of peculiar velocities. On the other hand, standard cosmological perturbation theory, LTB, and Szekeres Swiss cheese models display much larger fluctuations in the observed redshifts. In each of these cases the cause of the difference is the peculiar motion about the large-scale average (although the exact models listed here do not require a background in order to be defined). Such effects could be incorporated into the general framework of post-Newtonian cosmological modelling, but would require a more realistic distribution of matter in each cell.

4.2 Angular distance

We now turn to angular diameter distance measurements. In this case, the FLRW relation between angular distance $D_{\text{A}}$ and affine parameter $\lambda$ is the solution of the following differential equation:

[TABLE]

with $\bar{D}_{\text{A}}(0)=0$ and $\mathrm{d}\bar{D}_{\text{A}}/\mathrm{d}\lambda(0)=-1$. In figure 7 we compare the output $D_{\text{A}}(\lambda)$ we obtain from ray tracing within our simulations to the $\bar{D}_{\text{A}}(\lambda)$ expected from in a comparable FLRW model. Contrary to $z(\lambda)$, the difference can be significant, and in the galaxy simulations is much greater than the estimated numerical error.

In the galaxy simulation the mean distance $\left\langle D_{\text{A}}\right\rangle_{\Omega}$, measured across the observer’s sky at a given $\lambda$, is systematically greater than its FLRW counterpart. In other words, light sources are systematically demagnified compared to observations in a homogeneous Universe. This demagnification reaches $\sim 10\%$ for values of $\lambda$ that correspond to $z\sim 1$. The reason for this difference is that light never crosses the opaque matter clumps in these simulations, but instead always propagates through regions of perfect vacuum. This means that the sole source of focussing in the beam is due to the Weyl curvature of the space-time, with no contribution from the Ricci curvature (which vanishes in vacuum). This is the exact opposite of what happens to a beam of light in an FLRW geometry, where Ricci focussing is always non-zero while the Weyl curvature vanishes. The point of relevance here is that the shear that is caused by the Weyl curvature is much less efficient than Ricci curvature at focussing beams of light (except when light passes very close to clumps of matter). It is this lack of focussing that is responsible for the behaviour seen in figure 7. The dispersion of our data, represented by the error bars in figure 7, is due to the fact that some light beams pass closer to galaxies than others, and are thus more sheared. Such behaviour is reminiscent of what occurs in Einstein-Straus models [32]. These issues will be discussed in more detail below.

In the halo simulation, on the other hand, light sources are not systematically demagnified. In fact, they are even slightly magnified on average (see section 5.1). Since the haloes are transparent, the deficit of Ricci focussing that occurs when light propagates through empty space is compensated when the haloes are crossed. As well as this, the fact that the haloes are less compact means that they tend to produce weaker Weyl curvature, meaning that the beams are less sheared. The dispersion of the data, which is slightly smaller than in the galaxy case, is due to the fact that some lines of sight pass through more haloes than others. Unlike redshifts, distance measures are not directly affected by peculiar matter flows, but rather due to gravitational lensing phenomena which is itself related to the local space-time curvature experienced by the beams of light.

4.3 Hubble diagram

Let us now combine the results from sections 4.2 and 4.1, to construct the Hubble diagram that an observer would produce within our model space-time. By convention, it is usual to present this information as a plot of the distance modulus $\mu$, which is related to the luminosity distance $D_{\text{L}}$ by

[TABLE]

In figure 8 we compare the Hubble diagrams generated from our simulations, within a post-Newtonian model, to those expected in a comparable FLRW model. As expected from section 4.2, the differences are most pronounced in the case of galaxy simulations. This follows immediately from the result that redshifts are largely unaffected by the inhomogeneity, leading to $\delta D_{\text{A}}(z)\equiv D_{\text{A}}[\lambda(z)]-\bar{D}_{\text{A}}[\bar{\lambda}(z)]\approx\delta D_{\text{A}}[\bar{\lambda}(z)]\,.$ Note, however, that this approximate equality is no longer true when second-order terms are taken into account. Such terms can bias the Hubble diagram in a non-trivial way, with the quantitive magnitude of the effect depending on the particular distance measure being considered (luminosity distance, luminous intensity, distance modulus, etc.) [76, 77, 78].

It is interesting to note that in the galaxy simulation the effect of inhomogeneity on the Hubble diagram acts to bias the inferred amount of dark energy. That is, if we fit the data from our galaxy simulation to a spatially flat $\Lambda$CDM model by minimizing the $\chi^{2}$ statistic,

[TABLE]

then we find that best-fitting model has $\Omega_{\text{m}}=0.22$, rather than the value of $0.3$ that would be obtained if one had direct access to the background rate of expansion. In performing this fit we have taken $z_{i}=0.1i$ and $\sigma^{2}(z_{i})=\sigma_{\text{int}}^{2}+\sigma_{\text{lens}}^{2}(z_{i})$, where $\sigma_{\text{int}}=0.1$ is the intrinsic dispersion that is expected to encapsulate the astrophysics of supernovae in real data [80] and $\sigma_{\text{lens}}$ is the dispersion due to gravitational lensing (corresponding here to the error bars in figure 8). This is comparable to similar results found in the Einstein-Straus Swiss cheese model [30]. It implies that if one fits the observed Hubble diagram in such a universe, by wrongly assuming that it is homogeneous and isotropic, then one overestimates its actual dark energy content by about $\Delta\Omega_{\Lambda}=0.08$. When applied to supernova data, this effect tends to improve agreement with Planck [31] on $\Omega_{\text{m}}$.

Finally, let us close this section with a remark on the scatter of supernova data due to gravitational lensing in the joint light-curve analysis (JLA) [79]. The results of the JLA, which assume $\sigma_{\text{lens}}(z)=0.055z$, are shown in figure 9 together with our numerical results from ray tracing. The JLA ansatz can be seen to lie somewhere between our halo model and galaxy model. The shape of the JLA ansatz is closer to the halo case, which starts to become linear around $z=0.4$, while in the galaxy case the behaviour is more complicated. This is not surprising as the JLA ansatz is motivated by standard cosmological perturbation theory, where lensing is mostly due to fluctuations in the Ricci curvature (as is the case in our halo simulation). This can be seen as a confirmation of the JLA approach, but also warns that more complicated effects can occur in a universe that contains high-density compact astrophysical structures. However, the reader may note that the precise values of $\sigma_{\text{lens}}$ in our simulations must be taken with some caution as they strongly depend on the parameters of the model, as discussed in section 5.2.

5 Discussion

In section 4 we presented the basic results from our ray-tracing code, in terms of redshifts and distance measures. The difference between the averages of some of these data sets, and the expectations from FLRW models, were sometimes seen to be significant. In addition, the variance within each data set has the potential to allow interesting information about the structure in the Universe to be extracted from (for example) supernova observations. It is therefore of considerable interest to be able to relate these results to the underlying features of the model, in order to be able to understand how different features of the each configuration affects cosmological observations, and to be able to extrapolate these results from our idealized models to the real Universe.

To this end, in this section, we will interpret the results from section 4 in terms of some of the most prominent frameworks that have been constructed to understand the effects of inhomogeneity on the propagation of light. In particular, we will investigate (i) the theorem developed by Weinberg and Kibble & Lieu that relates the mean of observables to FLRW expectations, (ii) the stochastic lensing formalism that aims to model the statistical distribution of observations made along many lines of sight, and (iii) the relationship with the empty beam approximation along special lines of sight, as recently studied using full numerical relativity. This will allow us to probe the regimes in which these constructions and approximations are valid, as well as provide ways to understand the results of the previous section within well-established frameworks.

5.1 Weinberg-Kibble-Lieu theorem

In 1976 it was argued by Weinberg that gravitational lensing due to inhomogeneity in the Universe does not change the average relation between luminosity distance and redshift [81]. However, it was later realised that this very general statement could not be true; Weinberg overlooked the importance of choosing the observable that one averages, as well as the importance of the choice of averaging procedure itself. That is, averaging luminosity distances is not equivalent to averaging magnitudes or luminous intensities, and averaging over angles is not equivalent to averaging over sources or area [82, 83, 84, 76, 77, 78]. A more detailed and accurate statement was proposed in 2005 by Kibble & Lieu [85], which stated that it is the directional average of the angular diameter distance squared on surfaces of constant $\lambda$ that is unaffected by the presence of matter inhomogeneities. Mathematically, this can be expressed as follows:

[TABLE]

Equivalently, one can say that the total area of surfaces of constant $\lambda$, $A(\lambda)=4\pi\left\langle D_{\text{A}}^{2}(\lambda)\right\rangle_{\Omega}$, is unchanged by gravitational lensing. We shall call this property the Weinberg-Kibble-Lieu (WKL) theorem. If one further assumes that the $z(\lambda)$ relation is not significantly affected by inhomogeneities, the WKL theorem can be re-written in terms of surfaces of constant $z$.

We can now test the validity of eq. (5.1) in our PN simulations, by comparing $\left\langle D_{\text{A}}\right\rangle_{\Omega}$ and $\left\langle D_{\text{A}}^{2}\right\rangle_{\Omega}$ to their FLRW counterparts. The results are shown in figure 10, for both the galaxy and halo simulations. The WKL theorem is clearly violated in the case of the galaxy simulation, while it appears an extremely good approximation in the case of the halo simulation. This difference is essentially due to the opacity of the galaxies, which means that the light rays cannot sample every point in space. The region that is excised in the galaxy simulation is the region of highest Ricci curvature, meaning that the average density of matter along the allowed lines of sight is lower than the cosmological average (or strictly zero, in the idealized simulation we have performed). This biases the average of the inferred measures of distance, and causes the discrepancy that can be seen in figure 10. The halo simulation, on the other hand, has no regions of space that the rays of light are forbidden from entering. The average density they experience is therefore much closer to the cosmological average. This is the essential reason why the WKL theorem works so well in the halo case222The reader may note that the accuracy of these results is much greater than the random error on any individual trajectory; we attribute this to the error on the mean decreasing as the inverse square root of the number of trajectories in the sample (which is $10^{5}$ here)..

The logic above is not only applicable to our post-Newtonian models alone, but should apply to every model in which opaque structures are present (including the real Universe). In particular, similar behaviour has been observed in the context of Einstein-Straus Swiss cheese models [32]. A further property that causes violation of the WKL theorem, that is not present in our models but may exist in more general setups, is statistical anisotropy [86]. This can cause additional Weyl lensing, which means that $\left\langle D_{\text{A}}^{2}(\lambda)\right\rangle_{\Omega}\neq\bar{D}_{\text{A}}^{2}(\lambda)$. Finally, the averaging over angles $\left\langle\ldots\right\rangle_{\Omega}$ that appears in eq. (5.1) is not always the relevant averaging procedure to perform for every cosmological observable; it is the natural notion of averaging that one would use for full sky observations such as the CMB [87, 45], but is not the type of averaging that one usually uses when constructing a Hubble diagram from observations [78]. For this latter case it would be more natural to perform an average over sources, which could lead to violations of the relations one might expect from naive extrapolation.

5.2 Comparison with stochastic lensing

A useful formalism for modelling the effects of small-scale inhomogeneities on the propagation of narrow light beams has recently been constructed by modelling lensing events as stochastic processes [88]. In a nutshell, this approach consists of modelling the fluctuations in the Ricci and Weyl focussing scalars as white noise. The moments of the distribution of angular diameter distance measures can then be calculated, and compared to observations and ray-tracing models. It is instructive to compare this approach to the results we obtained in section 4, in order to evaluate its efficacy, and in order to determine the features of the real Universe that it should be expected to faithfully capture.

The first step in the stochastic formalism of ref. [88] is to decompose Ricci focussing term from equation (2.15) as $\mathscr{R}=\left\langle\mathscr{R}\right\rangle+\delta\mathscr{R}$, where $\left\langle\mathscr{R}\right\rangle=-4\pi G(1+z)^{2}\left\langle\rho\right\rangle$ is the mean Ricci focusing associated with the mean matter density $\left\langle\rho\right\rangle$ encountered by the light beam, while $\delta\mathscr{R}$ encodes fluctuations about it. There is no need for such a decomposition of $(\mathscr{W}_{1},\mathscr{W}_{2})$, as in a statistically isotropic universe $\left\langle\mathscr{W}_{A}\right\rangle=0$. The stochastic processes, $\delta\mathscr{R}$ and $\mathscr{W}_{A}$, are then characterised by covariance functions $C_{\mathscr{R}}$ and $C_{\mathscr{W}}$, defined by

[TABLE]

where $\delta(\lambda-\lambda^{\prime})$ is the Dirac delta function and $\delta_{AB}$ is the Kronecker delta. From these assumptions, it is now possible to calculate the mean and variance of the angular diameter distance along different lines of sight.

In particular, the lowest-order contribution to the average of the angular distance $\left\langle D_{\text{A}}\right\rangle$ due to shear can be written as [88]

[TABLE]

where $D_{0}$ denotes the angular diameter distance in a homogeneous space with $\left\langle\mathscr{R}\right\rangle$ and $\mathscr{W}=0$. Next, the variance of angular diameter distances, $\sigma_{D_{\text{A}}}^{2}(z)\equiv\left\langle D_{\text{A}}^{2}(z)\right\rangle_{\Omega}-\left\langle D_{\text{A}}(z)\right\rangle^{2}_{\Omega}$, can be shown to satisfy the following differential equation [88]:

[TABLE]

where it has been assumed that $z(\lambda)$ is the same as in the corresponding FLRW model, and where a prime denotes a derivative with respect to $z$. It can be seen that both the Ricci term and an integrated Weyl term act as sources for $\sigma_{D_{\text{A}}}^{2}$.

It now remains to determine the covariance functions $C_{\mathscr{W}}$ and $C_{\mathscr{R}}$. In a universe randomly filled with spherical opaque clumps of matter with mass $M$ and size $R$, the associated Weyl covariance is found to be [88]

[TABLE]

where $g=GM/R^{2}$ is the gravitational acceleration at the surface of the clump. This result was originally derived in the context of Swiss cheese models, but can be easily generalised to the present case. Finally, it can be shown that for the galaxy simulations we have $C_{\mathscr{R}}=0$, while for halo simulations we have

[TABLE]

We now have all the information required to compute $\left\langle D_{\text{A}}(\lambda)\right\rangle$ and $\sigma_{D_{\text{A}}}^{2}$ for our two simulations, and to compare these to our results from ray tracing.

In figure 11 we compare the mean angular diameter distance from our ray tracing experiment to the results from the stochastic lensing formalism. We see that the prediction from equation (5.4) reproduce the results of our simulations very accurately. In the galaxy case we have also indicated the results that one would obtain from the empty-beam approximation, in which $\mathscr{R}=0$. We see that this empty beam approximation, although not very good, is still better than FLRW at modelling the results from ray tracing. This confirms that, in the case of our galaxy simulations, most of the departure from FLRW is caused by the deficit of Ricci lensing due to the opacity of the matter clumps. When the stochastic shear corrections are added to $D_{0}$, the numerical results are recovered to high accuracy. In the halo case, where the FLRW model is already a good model for $\left\langle D_{\text{A}}\right\rangle_{\Omega}$, the stochastic shear correction from eq. (5.4) can be seen to improve the fit even further.

Finally, in figure 12 we compare the variance of our ray tracing results with the results obtained from stochastic lensing. We see that equation (5.5) provides an excellent description of the halo simulations, but fails in the galaxy case. The inability of stochastic lensing to predict the correct variance of the angular distance in a universe where light passes through vacuum regions is caused by the strong non-Gaussianity of $\mathscr{W}$, which cannot be reasonably modelled as a white noise. This supports similar findings in the case of the Einstein-Straus Swiss cheese model [88].

5.3 Comparison with numerical relativity

A recent field of study, within the field of inhomogeneous cosmology, is the study of models created using the tools from numerical relativity [64, 65, 66, 67, 68, 69, 70, 89, 90, 91, 92, 93, 94, 95]. The prospect of creating high-accuracy inhomogeneous cosmological space-times without any continuous symmetries is a highly exciting one, and presents the prospect of testing a number of the ideas around the possible consequences of inhomogeneity in cosmology in a rigorous and well-defined way. Nevertheless, the construction of numerical cosmological models is not without its difficulties. Numerical artefacts in these solutions do not dissipate in the same way that they do in asymptotically flat space-times, and the amount of computational time required to calculate the space-time geometry (and observables within it) is not trivial.

These facts make comparison between numerical relativity and approximation schemes potentially useful for those who are trying to eliminate the effects of numerical errors and artefacts in their simulations. In the future, such a comparison may also provide a way for us to try and evaluate the extent to which approximation schemes such as ours can be used as valid approaches for trying to create more versatile frameworks for capturing the lead-order effects of inhomogeneity. To these ends, we will compare our ray-tracing results to the only results of calculating observables within a numerical cosmological model that are currently available; those of ref. [72]. Although the two approaches to cosmological modelling are formally quite different, the properties of the physical configuration they are trying to simulate are very similar: In both cases indeed the universe is periodic and made of compact matter clumps.

As in our models, the authors of ref. [72] considered a universe tiled with cubic lattice cells. They then consider two fiducial light beams, $A$ and $B$, as depicted in figure 13. These beams are very special, in that they are stable during their propagation; the periodicity of the model implies that neither of them can be deflected by the masses at the centre of the cells. This is guaranteed by the fact that reflection symmetric surfaces are always totally geodesic, and that each beam lies at the intersection of at least two reflection symmetric surfaces. Along each of their trajectories they then computed the $D_{\text{L}}(z)$ relation, which is a direct observable for anyone inside the lattice.

In ref. [72], the authors then compared the results of numerically integrating the optical equations in their numerical space-time with the following well-known cosmological models:

(i) Einstein-de Sitter (EdS):

An FLRW model with $\Omega_{\text{m}}=1$ and $\Omega_{K}=\Omega_{\Lambda}=0$;

(ii) $\Lambda$CDM:

The concordance FLRW model with $\Omega_{\text{m}}=0.3$, $\Omega_{K}=0$, and $\Omega_{\Lambda}=0.7$;

(iii) The Milne universe:

An FLRW universe with $\Omega_{\text{m}}=0$, $\Omega_{K}=1$, and $\Omega_{\Lambda}=0$;

(iv) The empty-beam approximation (EBA):

A cosmology with $D_{\text{L}}=-(1+z)^{2}\lambda$, and the EdS relation between redshift and affine parameter $z(\lambda)$.

They found that the last of these approaches is the best (and a good) fit to their numerical results. Initially, this may seem at odds with the results we presented in the left panel of figure 11, but we will show it is entirely consistent.

In figure 14, we present the results of considering the special light beams, $A$ and $B$, in our post-Newtonian simulations (without the cosmological constant, in order to facilitate a direct comparison with ref. [72]). Indeed, it can be readily observed that the EBA does indeed provide the best-fitting $D_{\text{A}}(z)$ along each of these beams in our simulations, just as it does in the numerical relativity simulations of ref. [72]. In fact, in the case of beam $A$ this result is exact, as along the edge of a cell tidal forces must be exactly zero due to discrete rotational symmetry. Furthermore, the beam does not encounter any form of matter along the cell edge, so both Ricci and Weyl lensing terms must vanish (as assumed in the EBA). For beam $A$, the apparent discrepancy with the results presented in figure 11 can now be seen to be entirely due to the very special choice of geodesics $A$, which is not representative of the experience of a typical ray of light within the model. In particular, beam $A$ experiences a drastically lower rate of shear than almost all other trajectories.

Let us now consider beam $B$, which is slightly focused with respect to the EBA. This is because along $B$ the tidal forces due to the masses above and below are not cancelled by the masses on the sides, which are further away. This beam is thus consistently sheared along the $z$ axis, which causes a correction with respect to the EBA. It is interesting to note that, even though shear is consistently produced along the same axis, the difference between $D_{\text{A}}(z)$ along the trajectory B and the EBA remains very small ($5\times 10^{-5}$ at $z\sim 1$). This can be compared with the left panel of figure 11, where the difference between $\left\langle D_{\text{A}}\right\rangle_{\Omega}$ and the EBA is at the level of one percent at the same redshift. We attribute this to the fact that the beam is always far from the central mass of each cell it passes. Consequently, this means that the majority of the shear corrections in $\left\langle D_{\text{A}}\right\rangle_{\Omega}$ must come from very close encounters with the central masses, and not from special trajectories where the contributions to shear are all cumulative. We conclude from this that the results presented in section 4 are not artefacts of the lattice structure of the model, but should be quite robust under changes of configuration.

6 Conclusions

In this paper we have ray traced through a set of models that were created to simultaneously model both large-scale expansion and non-linear structure. We have found that the redshifts that an observer in such a space-time would model are very close to those of an FLRW model with the same amount of total mass and dark energy. This result was found to be independent of the way in which matter is distributed within each of the primitive cells of our lattice model, and whether or not the matter is opaque or transparent. On the other hand, we have found that angular diameter and luminosity distances can differ considerably from the predictions of FLRW cosmology, as long as the mass is clumped into compact and opaque non-linear structures. In the extreme case, when all matter is clustered into high-density galaxy-bulge-sized objects we find that the difference can be as much as 10%, with objects appearing dimmer at the same redshift in the inhomogeneous cosmology. This difference drops dramatically when the matter is taken to be dispersed and transparent, however, and becomes less that 1% when dark matter halo-sized objects are being modelled.

When the results from considering dense opaque objects are fitted to FLRW cosmological models, we find that the best fitting spatially flat model has $\Omega_{\Lambda}=0.78$ when the actual value in the model is $\Omega_{\Lambda}=0.7$. This constitutes a bias of more than $10\%$ in the estimation of this crucial cosmological parameter, compared to the perfectly homogeneous and isotropic case, which is a considerable effect. Beyond the difference in the mean, the inhomogeneities also add a variance to the distribution around the mean. In the case of clumped opaque matter this is at the level of about $\pm 5\%$ for cosmologically interesting distances, while for diffuse transparent structures it is more like $\pm 3\%$. These are considerable deviations, that could in principle be used to extract information about inhomogeneity from the distribution of, for example, supernova observations.

By comparing our results to the expectations from the Weinberg-Kibble-Lieu theorem, we find that when the matter is transparent and diffuse the relevant averages in the inhomogeneous universe correspond to the FLRW values in a very precise way. When matter is opaque and compact, on the other hand, there is no such correspondence, and the averages of the relevant optical probes in the inhomogeneous universe do not correspond to FLRW values in any obvious way. We find that the recently developed stochastic approach to lensing is excellent at modelling the average of the luminosity distances that we calculate in both the opaque and transparent distributions of matter. This formalism is, however, rather less good at predicting the variance around the mean in the presence of opaque clusters of matter. Finally, we calculate observables along the special curves that have recently been considered in the context of numerical relativity simulations of inhomogeneous cosmological models. Our results are free from the oscillatory numerical artefacts that exist in these models, and could serve as a guide for the expected results that could be obtained once they have been brought under control.

The models that we have been considering should not, of course, be considered to be fully realistic representations of the actual Universe. This is primarily due to the periodicity that is assumed on large scales, in order to make their construction a tractable problem. While the real Universe does appear to display a homogeneity scale, and while similar periodicity assumptions are made in standard $N$-body simulations, this does not imply the stronger assumption of reflection symmetry that we have assumed. This means that the post-Newtonian cosmologies we have considered should be thought of as toy models. Nevertheless, these model do retain the distinct advantage of allowing relativistic effects to be consistently included in cosmology in the presence of non-linear structures, and in this regard they are much more flexible than other comparable approaches (such as Swiss cheese models or black hole lattices).

Due to the idealised nature of our model, we do not expect the $10\%$ bias in $\Omega_{\Lambda}$ discovered in our galaxy simulations to be fully representative of the value of the corresponding effect in the real Universe. This is because the real Universe contains both opaque galaxy bulges and dark matter haloes (as well as various other astrophysical matter). In this sense, the galaxy and halo simulations we have performed should be considered as extreme examples that provide upper and lower bounds on the magnitude of the effect that will occur in the real Universe. Given that we have no reason to expect the periodicity of our model to introduce any further biases beyond those discussed above, and given there is more dark matter than visible matter in the Universe, we estimate that in reality the astronomical selection effects should bias the estimation of $\Omega_{\Lambda}$ at about the level of $\sim 1\%$. This is a significant effect, in the age of precision cosmology.

In subsequent work we will consider the detailed statistical behaviour of the observables of gravitational lensing, such as optical shear and distance measures, as well as including more realistic distributions of matter. A first step in doing this would be to include opaque galaxy cores and, e.g., a Navarro-Frenk-White profile. A second step would be to arrange these haloes within a cell according to the actual matter power spectrum, and allowing for their peculiar motion under gravity. This should provide further insight into the biases that inhomogeneous structures can cause, as well as allowing interesting questions about the formation of caustics and the statistics of optical shear and distance measures to be addressed. If caustics form in any considerable number, then they could have profound consequences for both the construction of Hubble diagrams and CMB observations [96].

Acknowledgements

PF thanks Queen Mary University of London for its great hospitality at several stages of this project. VAAS would like to thank John Ronayne, Pedro Carrilho, Mike Cole, Cozmin Timis and Alex Owen for very helpful discussions and advice. We also thank Eloisa Bentivegna, Mikołaj Korzyński, Syksy Räsänen, Ian Hinder and Jean-Philippe Uzan for comments on the draft version of this article. VAAS and TC are both supported by the STFC. PF is supported by the Swiss National Science Foundation.

Appendix A Perturbed geometric quantities

At order $\varepsilon^{2}$ in the post-Newtonian expansion, the metric of space-time is given in equation (3.1). The perturbed Christoffel symbols associated with this geometry are then

[TABLE]

and the Riemann tensor reads

[TABLE]

These are all the components required for our numerical integrations.

Appendix B Tests in de Sitter space-time

In order to test the accuracy of our ray-tracing code, we considered the case of a de Sitter space-time, i.e. a universe with no matter but with a non-zero cosmological constant. In this case the potentials in the post-Newtonian metric (3.1) are given by

[TABLE]

We can now simulate observations within this geometry and compare them with the known exact expressions in de Sitter:

[TABLE]

which combine to give

[TABLE]

where $H_{0}=\sqrt{\Lambda/3}$. Recall that we chose the wave four-vector $\bm{k}$ to be future oriented, hence $\lambda\leq 0$ in the past. The relative difference between the output of our ray-tracing code and the exact results above is displayed in figure 15. It can be seen from this figure that the accuracy on the $z(\lambda)$ relation in our numerical implementation is at the level of about one part in $10^{6}$, which is even better than our estimates in section 3.4. The accuracy on the $D_{\text{A}}(\lambda)$ relation is two orders of magnitude better, and also converges as $\sim\Delta t/L$. Finally, the dispersion of the data is much smaller than the mean error. In other words, numerical errors appear like a systematic bias of the data more than a random process. We interpret this fact as being due to the Euler integration, for which the local errors are quadratic, and hence are always cumulative (as they have the same sign).