Light propagation through blackhole lattices
Abstract
The apparent properties of distant objects encode information about the way the light they emit propagates to an observer, and therefore about the curvature of the underlying spacetime. Measuring the relationship between the redshift and the luminosity distance of a standard candle, for example, yields information on the Universe’s matter content. In practice, however, in order to decode this information the observer needs to make an assumption about the functional form of the relation; in other words, a cosmological model needs to be assumed. In this work, we use numericalrelativity simulations, equipped with a new raytracing module, to numerically obtain this relation for a few blackhole–lattice cosmologies and compare it to the wellknown FriedmannLemaîtreRobertsonWalker case, as well as to other relevant cosmologies and to the EmptyBeam Approximation. We find that the latter provides the best estimate of the luminosity distance and formulate a simple argument to account for this agreement. We also find that a FriedmannLemaîtreRobertsonWalker model can reproduce this observable exactly, as long as a timedependent cosmological constant is included in the fit. Finally, the dependence of these results on the lattice masstospacing ratio is discussed: we discover that, unlike the expansion rate, the relation in a blackhole lattice does not tend to that measured in the corresponding continuum spacetime as .
a,b]Eloisa Bentivegna, c]Mikołaj Korzyński, d]Ian Hinder, d,e]and Daniel Gerlicher Prepared for submission to JCAP
Light propagation through blackhole lattices

Dipartimento di Fisica e Astronomia
Università degli Studi di Catania
Via S. Sofia 64, 95123 Catania
Italy 
INFN
Sezione di Catania
Via S. Sofia 64, 95123 Catania
Italy 
Center for Theoretical Physics
Polish Academy of Sciences
Al. Lotników 32/46, 02668 Warsaw
Poland 
MaxPlanckInstitut für Gravitationsphysik
AlbertEinsteinInstitut
Am Mühlenberg 1, D14476 Golm
Germany 
Technische Universität München
Boltzmannstrasse 15, D85748 Garching
Germany
Keywords: cosmological simulations, ray tracing, gravity, GR black holes
Contents
1 Introduction
General relativistic spaces filled with black holes have recently been under scrutiny as exact cosmological models with a discrete mass distribution which is, in some sense, uniform on large scales. The construction of these spaces in numerical relativity has enabled the investigation of several questions without approximations, such as how such configurations evolve in time and what their global physical properties are [1, 2, 3, 4, 5]. At the same time, the numerical simulations have been complemented by insight coming from analytical studies, which have illustrated some general features of these spacetimes such as the behaviour of special submanifolds [6, 7, 8], the conditions under which they behave like the FriedmannLemaîtreRobertsonWalker (FLRW) models [9], and the link between their behaviour and the validity of Gauss’s law in a generic theory of gravity [10].
In this work, we use numerical spacetimes representing blackhole lattices (BHLs) to probe a different aspect of inhomogeneous cosmologies, namely their optical behaviour. As is well known, null geodesics are the bedrock of cosmological observations: light from distant sources is the primary tool for measuring the Universe’s density parameters, equation of state, and perturbations. Increasing the accuracy of models of light propagation and identifying the biases introduced by various approximation frameworks is thus critical.
Modelling light propagation in inhomogeneous cosmologies is a longstanding effort, which has followed two complementary courses: approximation schemes on one hand, and toy models on the other. The bestknown approach in the former class is the EmptyBeam Approximation (EBA) of Zeldovich [11], later generalized by Dyer and Roeder [12, 13]. This approach is based on the idea that different effects are at play when light propagates in a perturbed fluid or through discretelydistributed point masses, as different components of the curvature become dominant in either regime (this is sometimes referred to as the RicciWeyl problem [14]). This approach provides an excellent estimate of light propagation in SwissCheese models, and can be used to constrain the fraction of voids in a cosmological model [14, 15, 16]. The notion that discreteness may affect light propagation more than inhomogeneity itself has also appeared in other studies, such as those on light propagation through Schwarzschildcell universes [17, 18, 19].
The existing literature points in a number of common directions: first, examining individual geodesics, one concludes that the effective value of the cosmological constant (the one obtained fitting the spacetime to an FLRW model with the same matter density) is higher than its microscopic value (the one appearing in the gravitational action). Second, a statistical average of photon trajectories usually leads to a partial suppression of this difference. A suppression is also obtained by considering the perturbative solution corresponding to a regular arrangements of objects of equal mass, at least until the perturbative condition is respected [20].
Though consistent on many aspects, these studies are limited by the conditions imposed on the underlying model: most of the discretemass studies are either based on sphericallysymmetric building blocks or on the requirements that the objects be not too compact. It is presently not clear what the optical properties of a more generic space would be.
To investigate this issue and test the generality of the existing results, in this paper we compute the photon redshift and luminosity distance along null geodesics running through a BHL spacetime, constructed exactly through the integration of Einstein’s equation, nonperturbatively and in three dimensions. First we compare the result to some reference models from the FLRW class, to the Milne cosmology, and to a generic universe in the Empty Beam Approximation (EBA) [11, 12, 13]. We find that the latter provides the closest approximation to light propagation on the BHL, and derive a simple argument to explain this result, which in some sense extends the reasoning of [16] to completely vacuum spacetimes. We then turn to the question of whether it is possible to tune the cosmological parameters in the FLRW class to improve the fit. We find, in particular, that one can reproduce the luminositydistance–to–redshift relationship of a BHL with that of an FLRW model with the same average matter density and a fictitious, timedependent cosmological constant , and provide the first measurement of this running in our base configuration. Finally, we study how this behaviour depends on the BHL inhomogeneity parameter [4], which roughly corresponds to the ratio between the central mass and the lattice spacing, and in particular we analyse the continuum limit of .
An important factor in this discussion is the choice of light sources and observers, as the photon frequencies and number counts will depend on the reference frame in which they are measured. In FLRW models there is an obvious option: the comoving sources and observers. In inhomogeneous spaces, on the other hand, identifying a “cosmic flow” is more tricky (when at all possible) and relies on the somewhat arbitrary split between global cosmological evolution and “local effects” sourced by nearby gravitational structures. For the purpose of this work, we sidestep this question by noticing that, for a given geodesic, the angular and luminosity distances can be obtained by applying a certain linear operator to the fourvelocity of the observer, with no dependence whatsoever on the motion of the light source. It is therefore straightforward to quantify the effect of different observer prescriptions on these observables.
Section 2 introduces the formalism of light propagation and justifies the approach we take in our analysis, providing some examples in simple spacetimes. Section 3 provides an approximate description of light propagation in a BHL via a perturbative analysis. We present the numerical results in section 4 and in section 5 we comment on them. We provide tests of the geodesic integrator, used for the first time in this study, in the appendix. We use geometric units everywhere.
2 Fundamentals of light propagation
Let us start by considering a null ray emanating from a light source and reaching an observer : this curve can be described as an affinelyparametrized null geodesic , with and as end points corresponding to the affine parameter values and :
(2.1)  
(2.2) 
The curve is described by the geodesic equation:
(2.3) 
where:
(2.4) 
is the tangent vector to . In order to measure distances with null rays, however, we need more than a single geodesic: we need to consider a whole beam of rays [21], centred on , and study the evolution of its crosssectional area as it makes its way from to .
The time evolution of a beam’s cross section is described by the geodesic deviation equation (GDE). Let be the separation vector between the fiducial geodesic and a neighbouring one, called . It satisfies
(2.5) 
The GDE is a second order ODE for the 4–vector , or equivalently a first order ODE for and . It is valid for any neighbouring geodesic, but since in the geometrical optics we are only interested in null geodesics, we impose a restriction on the solution of the form:
(2.6) 
which ensures that is null. Note that if the equation above is satisfied at one point, then it is automatically satisfied along the whole of because of equation (2.5).
Let us now restrict the geodesics under consideration to those which lie on the same wavefront as , i.e. for which the separation vector satisfies
(2.7) 
The condition above means that, for a given observer at a given time, the photon corresponding to the geodesic and the one corresponding to lie on the same 2plane perpendicular to the direction of propagation (see Figure 1). This condition is Lorentzinvariant, meaning that if it is satisfied in one reference frame then it is valid in all frames. Moreover, for null geodesics it propagates along , i.e. if it is satisfied at one time it is satisfied along the whole of . This follows easily from (2.6) and (2.5).
The reason why we are interested only in geodesics which lie on the same wavefront is that we want to study geodesics which cross at one point, either the emission point or the observation point . If this is the case, then at either or , so that (2.7) is trivially satisfied there and thus also everywhere on .
By imposing (2.6) and (2.7) we have effectively reduced the number of degrees of freedom from four to three. It turns out that a further reduction is possible. Note that at every point we are free to add a vector proportional to to both and . The former corresponds to using a different point of the same geodesic in the definition of the separation vector , while the latter is just a rescaling of the affine parametrization of . Neither transformation affects the physical content of the equations, as long as we are in the regime of geometrical optics. As a matter of fact, it is easy to see that equations (2.5)–(2.7) are insensitive to these transformations as well:
(2.8)  
(2.9)  
(2.10) 
It follows that (2.5)–(2.7) can be reinterpreted as equations on the space , consisting of vectors orthogonal to and divided by the relation . We shall denote the equivalence class corresponding to a vector in as . The space is two–dimensional and inherits the positivedefinite metric from via the relation , where and are any vectors in the tangent space corresponding to the equivalence classes and , respectively. It can be thought of as the space of null geodesics lying in the neighbourhood of on the same wavefront, without any specification of which point on we assign to which point of . It is straightforward to verify that the covariant derivative can also be defined as an operator on .
In the standard formalism due to Sachs [22, 23], we then introduce a frame with two spatial, orthonormal screen vectors and , both orthogonal to and to a timelike observer . Notice that this is not strictly necessary: all that matters in geometrical optics are the equivalence classes and , which turn out to be entirely independent. More precisely, for any other choice of the observer and the corresponding and perpendicular to , the classes and are related to and via a simple spatial rotation.
The image distortion of a distant object and its angular distance can now be calculated by finding the Jacobi matrix of the GDE in the space
(2.11) 
with the initial data of the form
(2.12)  
(see [23] for its geometric definition and the discussion of its properties). Note that the initial data depends on the choice of parametrization of the null geodesic , because if we rescale , the tangent vector rescales accordingly via . Thus is parametrizationdependent. Nevertheless, the tensor product is parametrizationindependent and is therefore an intrinsic property of the light cone centred at the observation point . In practice the equations (2.11)–(2.12) are solved by first introducing a Sachs frame and then using the corresponding screen vectors and as a basis in .
The image distortion seen by the observer with 4velocity at the observation point is finally:
(2.13) 
while the angular distance is
(2.14) 
(see also [23] and references therein). Note that the result does not depend on the 4velocity of the source, while the dependence on the 4velocity of the observer is quite simple. For instance, it is easy to prove that, on an FLRW spacetime, observers boosted with respect to the comoving frame measure smaller angular distances, because the quantity decreases as the boost parameter is increased. One can therefore use equation (2.14) to work out which observers (if any) would measure a specified angular distance for an object in a given spacetime.
The luminosity distance is defined using the total energy flux from the source through a fixed area at the observation point. In the formalism above it can be expressed as
(2.15) 
where satisfies (2.11), but with the initial conditions (2.12) imposed at the source rather than at the observer, and is the relative change in the photon frequency as it moves along the geodesic, also known as its redshift:
(2.16) 
The fundamental result by Etherington [24] relates these quantities: the reciprocity relation reads
(2.17) 
It follows easily that
(2.18) 
Relation (2.17) allows one to calculate both distances by solving the GDE with the initial conditions (2.12) imposed either at the source or at the observation point.
In this paper we have found it much simpler to impose the initial conditions at the location of the source, and to integrate the equations forward in time. Moreover, instead of solving the GDE directly, we simply use the geodesic tracker and follow directly two additional null geodesics and , slightly perturbed with respect to the principal one, which we denote with . We specify the initial conditions for them at the source:
(2.19)  
(2.20)  
(2.21) 
where are the coordinates of geodesic and is its 4momentum. We can then compute by using the fact that:
(2.22) 
where is the determinant of at the geodesic initial location. This is the approach we take in the computations described in Section 4.
2.1 Homogeneous cosmologies
This formalism takes on a particularly simple form in the exactly homogeneous and isotropic cosmological models (the FLRW class), defined by the line element:
(2.23) 
where is the line element of one of the three threedimensional constantcurvature spaces of Euclidean signature. In this case, geodesics can move along coordinate lines and be parametrized by the coordinate time. In the flat case, for instance, we can choose as the geodesic direction (so that and , where is the scale factor). The matrix is then given by:
(2.24) 
where is the coordinate distance travelled along the geodesic at time :
(2.25) 
Given the initial normalization , equation (2.15) becomes:
(2.26) 
Noticing that, in an FLRW model, the redshift only depends on the ratio between the scale factor at the time of detection and the scale factor at the time of emission:
(2.27) 
it is easy to show that equation (2.26) coincides with the usual textbook expression for , which we quickly recall. We first need to calculate the comoving distance covered by a photon between and :
(2.28) 
with
(2.29) 
and
(2.30) 
Notice that the reference values for all quantities are those at the source: , , and are the model’s scale factor, Hubble rate, and density parameters at the time the photon is emitted, respectively. As is customary, we also define the curvature parameter by:
(2.31) 
Notice that referring to the initial values of these parameters rather than the final ones changes our expressions from the standard textbook treatment. It is straightforward to show that the usual formulae are recovered if one expresses all quantities at the source in terms of the corresponding ones at the observer.
Having found an expression for , we can use it to derive the apparent luminosity of an object of intrinsic luminosity (for details, see e.g. [25]):
(2.32) 
Since the apparent luminosity is defined as:
(2.33) 
we finally obtain:
(2.34) 
This can be easily identified, on a flat background, with (2.26). In homogeneous and isotropic cosmologies, therefore, the luminosity distance only depends on the redshift, and is parametrized by global quantities such as the matter density and the curvature of spatial slices. In the Einsteinde Sitter (EdS) model, simply reduces to
(2.35) 
2.2 Inhomogeneous cosmologies
The propagation of light in lumpy spacetimes has been studied since the 1960’s with various approaches, starting with the EBA proposed in [11] and later generalized in [12, 13]. The key idea inspiring these studies is that, in cosmological models where the matter is distributed in lumps, a large fraction of the light beams would not contain matter, and would therefore not be affected by the Ricci focusing characteristic of their FLRW counterparts.
Other limitations of the FLRW approximation and the related physical effects were subsequently analysed, both in approximate scenarios and in exact cosmological models (typically belonging to the SwissCheese family) [16, 21, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. A few robust features of these studies, that do not depend on the details of the models used, include that:

Light sources appear reduced in size and dimmer in a lumpy spacetime than in a homogeneous one with the same mean density;

The angular distance does not have a maximum, but keeps growing all the way to the cosmic horizon;

The actual deceleration parameter is larger than in the case where the same data is analysed with an FLRW model with the same mean density.
Later, when we measure the relationship in BHL spacetimes, we will use these features as guidelines for what to expect. Many of them do indeed hold for such highly nonlinear spacetimes too.
In fact, as discussed at length in Section 4, the luminosity distance in a BHL follows rather closely the EBA [11], which we report for completeness:
(2.36) 
In Section 3 we will explain why the EBA is a good approximation of the redshift–luminosity distance in a BHL, and point out that it is equivalent to neglecting the Ricci term in the standard geodesic deviation equation.
2.3 Geodesics and observer classes
As with many other quantities of interest that can be calculated in inhomogeneous cosmologies, the calculation of requires the choice of a time coordinate. In general, representing the spacetime in the geodesic gauge will lead to coordinate observers which are diversely affected by neighbouring gravitational structures, and may experience, e.g., light redshifting which has nothing to do with a global, suitably defined expansion rate (an observational cosmologist would call these local effects).
A study of light propagation in inhomogeneous spaces, especially one that is targeted at the comparison with the FLRW class, is then left with two possibilities: a statistical approach in which observers and sources are distributed stochastically throughout the spacetime, and a single relationship is obtained by averaging over their locations and fourmomenta; or the construction of one or more classes of cosmological observers, based on geometryinspired considerations such as following the geodesics of the average gravitational field, or geodesics with minimal deviation. We find the latter approach more likely to yield insight on the different gauge choices and related effects, and therefore use it in the remaining of this paper. Statistical reasoning is, however, also an important ingredient, as the observational data is arguably to be modelled through a mix of different observer and source states of motion. As statistical analyses are a tricky endeavour in cosmology, we leave this task for future work.
Notice that the second strategy is particularly difficult to deploy on vacuum spacetimes, as the sources of the gravitational field are only perceived through their effect on the metric tensor, and not through the presence of matter, so singling out a “local” component of the gravitational field will in some cases not even be well defined (for a discussion of this point, see e.g. [40] or [23] and references therein). We will however exploit the existence of global (albeit discrete) symmetries in our BHLs and only turn our attention to geodesics which are by construction least affected by local effects: these include, for instance, the geodesics running along the edges of the fundamental periodic cell constituting the lattice.
3 Light propagation in BHLs
In this section, we build an approximate model for the propagation of light in a BHL, based on a perturbative expansion in the BHL compactness parameter. This will serve as a qualitative analysis of the physics of the propagation of light and as a support in the interpretation of the numerical results presented in section 4. Note that an expansion in a similar parameter has already appeared in the context of BHLs [20], although the details are different.
3.1 A perturbative expansion in the compactness parameter
Let denote the characteristic size of a lattice cell, such as its initial geodesic length, and let be a characteristic mass, i.e. the total mass contained in a cell. As in [4], we can introduce the dimensionless parameter
(3.1) 
measuring the lattice compactness. If we additionally introduce the characteristic mass density , we can see that
(3.2) 
i.e. it goes down to zero as we decrease the size of a cell keeping the mass density of the corresponding FLRW model fixed. Note that is related to the curvature scale of the Friedmann model via , so can be reinterpreted as the separation of scales between the size of an individual lattice cell and the radius of curvature of the FLRW model:
(3.3) 
Note that the definition of involves a certain vagueness: we may take for the mass scale the ADM mass of the black hole measured at the other end of the EinsteinRosen bridge, but also some other related parameter. Also the choice of the length scale involves a certain arbitrariness. At the leading order we expect this ambiguity to be irrelevant.
We will now show how can be used to find a perturbative approximation for the metric tensor of the lattice model. The approximation is different from the standard perturbative approximation on an FLRW background, in the sense that it does not require the density contrast of the matter perturbation to be small. Obviously the problem of BHLs lies beyond the validity regime of the cosmological perturbation theory, because in a BHL we are dealing with everywhere.
We begin by introducing a coordinate system on a single cell. Let denote the background FLRW metric and be the Riemannian normal coordinate system around any point . The metric takes the form of
(3.4) 
Since is the curvature scale of the metric, the coefficients in the expansion above are of order (the flat metric), (the Riemann tensor), (the next term involving ), and so on. The Taylor expansion in the Riemannian normal coordinates becomes thus the expansion in negative powers of . We now introduce the rescaled coordinates and the rescaled Riemann tensor at point in coordinates .
(3.5) 
Both and are in the expansion in , at least within a single lattice cell around . The metric can be expressed in the new coordinates. The metric tensor components in those coordinates will be denoted by , i.e.
(3.6) 
Its expansion in takes the form of
(3.7) 
The first of the remaining higherorder terms contains the covariant derivative , which is as we noted before. Therefore the whole term in question can be reexpressed as , where we have defined by analogy the rescaled derivative of the curvature , which again is in . We see that the whole term turns out to be . Similar reasoning can be applied to all higher terms, yielding higher powers of the dimensionless parameter . We thus see that
(3.8) 
i.e. in the rescaled coordinates the expansion in the negative powers of turns in a natural way into an expansion in powers of , valid in a region of size around .
We can explain the physical meaning of the expansion above in the following way: if the background metric has the curvature scale of , then in an appropriately picked, quasiCartesian coordinate system it has the Taylor expansion in which the terms are of increasing order in . If we then pick a domain of size , then the metric in this domain, again in appropriate coordinates, has the form of the flat metric plus perturbations from the curvature and its derivatives. A simple way to obtain a perturbation of this kind is to use the Taylor expansion we mentioned before and rescale the coordinates by , which yields an expansion in powers of .
Now we can add the perturbation due to the discrete matter content. We assume the full metric to be
(3.9) 
with the perturbation of order in . Note that the dependence on means that the characteristic physical size of the perturbation is the size of a cell, i.e. . The Einstein tensor of the metric above is
(3.10) 
where is the linearisation of the Einstein tensor around a flat metric . In particular, in the harmonic gauge it is simply . We now return to the original, unrescaled coordinate system, where this equation takes the form of
(3.11) 
i.e. the perturbation of the Einstein tensor is , just like the Einstein tensor of the FLRW metric. It means that this approximation works even if the density perturbation is of the order of the background energy density. We may therefore use to cancel the stressenergy tensor of the underlying FLRW metric everywhere except on a single worldline.
Recall that , where is the cosmic fluid 4velocity. We impose the linear PDE on the metric perturbation:
(3.12) 
with periodic boundary conditions and with the constant chosen so that the RHS integrates out to zero over one cell. The solution can be obtained using Appell’s function [41]. It diverges at the centre, where the approximation fails, but near the cell’s boundary it is likely to work well. The resulting approximate metric is vacuum everywhere and periodic.
3.2 The continuum limit
Let us now consider the metric (3.9) along with its Christoffel symbols and Riemann tensor. It is straightforward to see that
(3.13)  
(3.14)  
(3.15) 
We can now go back to the original, unrescaled coordinates and obtain
(3.16)  
(3.17)  
(3.18) 
plus higher order terms in . Consider now the limit , i.e. where the size of the perturbations decreases in comparison to the curvature scale of the background FLRW model, or the limit where the compactness vanishes. Obviously we see that the metric tensor and the Christoffel symbols converge to the FLRW values in this case, while the curvature does not. This is due to the fact that the metric is that of a vacuum spacetime for all positive , while the FLRW one is not. This is a key observation in the study of the optical properties of a BHL, which are determined by the GDE and are therefore sensitive to the form of the Riemann tensor.
To illustrate this point, consider first a null geodesic. It follows from equations (3.16)–(3.18) above that its equation has the form of a perturbed FLRW geodesic
(3.19) 
where the tilde denotes the FLRW solution without the inhomogeneities. The parallel transport of a frame along the geodesic has a similar expansion in :
(3.20) 
We can now rewrite the GDE in the parallelpropagated frame along the geodesic
(3.21) 
We see that, already at the leading order in , we must take into account the full physical Riemann tensor instead of the simple FLRW one. In particular, since the BHLs are vacuum spacetimes, we need to solve the Riccifree GDE and possibly take into account the nonvanishing Weyl tensor along the way in order to calculate the angular and luminosity distance. Neglecting the Ricci tensor in the GDE is equivalent to the EBA (for a discussion of this point, see e.g. [16]). We may thus expect the redshift–luminosity relations for BHLs in the continuum limit to be close to the EBA.^{1}^{1}1We neglect here the finitebeamsize effects which would become large when becomes very small: the beam may at some point become wide enough to encompass a large number of black holes. In this situation the interaction between the beam and the black holes becomes quite complicated as we cannot use the GDE approximation any more.
At the order we may expect additional corrections to and due to higherorder contributions to the geodesic equation as well as to the GDE equation. Additionally, at this order we need to take into account the impact of the inhomogeneities on the observers in their free fall. In this work, we will not concern ourselves with a quantitative analysis of these effects, but we will signal their appearance to the reader when appropriate.
4 Results
In order to compute the relationship between redshift and luminosity distance on the spacetime of an expanding BHL, we carry out the numerical integration of the geodesic equation (with null tangent), along with the integration of Einstein’s equation required to obtain the metric tensor. The latter operation is performed by a code generated with the Einstein Toolkit, based on the Cactus [42] software framework along with modules such as Carpet [43, 44], McLachlan [45, 46], and CT_MultiLevel [47], as already presented in [2, 4, 8]. The geodesic integrator, on the other hand, is a new Cactus module that we have written. It implements a 3+1 decomposition of the geodesic equation in the form given in [48] and we have verified it against several exact solutions, as reported in Appendix B.
4.1 Initial data and evolution
As in [1, 4] we first construct an initialdata configuration by solving the Hamiltonian and momentum constraints on the cube with periodic boundary conditions. In particular, we choose free data corresponding to conformal flatness:
(4.1) 
and set the trace of the extrinsic curvature to zero around the origin and to a negative constant near the boundaries, with a transition region starting at a distance from the origin:
(4.2)  
(4.3) 
where we choose and . We represent the traceless part of the extrinsic curvature as:
(4.4) 
and the conformal factor as:
(4.5) 
where is the bare mass of the central black hole, and solve the constraints for and . For our basic configuration, we use and as in [4].
We then proceed to the time evolution of and using a variant of the BSSN formulation, implemented in the McLachlan module, and to the concurrent integration of the geodesic equation (2.3).
4.2 Computation of geodesics
In order to compute geodesics in a 3+1 numerical spacetime, we first perform a 3+1 decomposition of the geodesic equation (2.3),
(4.6) 
We decompose the geodesic tangent vector into its components along and orthogonal to the unit hypersurface normal , which we call and , respectively: . The vector is spatial, i.e. , and . We use an affine parametrisation, and is normalized as
(4.7) 
with for null geodesics. The spatial coordinates, covariant components of the tangent vector and affine parameter of the geodesic, (, , ) satisfy
(4.8)  
(4.9)  
(4.10) 
where
(4.11) 
is the time component of in the foliationadapted coordinate basis. Note that the derivative is with respect to coordinate time , not the affine parameter . These equations are the same as those given in [48], and a derivation is outlined in Appendix A.
Given at a time , equations (4.8)–(4.10) determine their evolution along a single geodesic. The right hand sides of eqs. (4.8)–(4.10) are computed by interpolating the metric quantities , , from the evolution grid to the point using fourthorder Lagrange interpolation, and is integrated using a fourthorder RungeKutta method using the Cactus MoL component. Additionally, the metric and various other quantities of interest are interpolated to , and all quantities are output as curves parametrised by for use in any subsequent analysis once the simulation is complete.
We implement the above prescription in two new Cactus components Geodesic and ParticleUtils. The former contains the equations themselves, and the latter provides librarytype functionality for integrating systems of equations along curves. A few validation tests are provided in Appendix B.
We now face the crucial task of selecting which geodesics to track. Let us notice that, on a space filled with periodic cells, symmetry reasons imply that an obvious class of cosmological observers is that formed by observers sitting at the cell vertices. Due to the symmetry, these observers do not exhibit any proper motions on top of the cosmic expansion, and the ratio of the proper distances between arbitrary pairs of observers is constant at all times. For this study, we construct and analyse two geodesics from this class (which we will denote and ), starting at the vertex , with initial tangents equal to and respectively. and are chosen by the geodesic integrator to ensure that the geodesics are null. The two geodesics are plotted in Figure 2.
In order to measure the luminosity distance along geodesics and , we evolve two further pairs of geodesics, with spatial directions given by:
(4.12)  
(4.13) 
and
(4.14)  
(4.15) 
with , representative of two narrow beams close to each original geodesic. We can then construct the redshift and luminosity distance along the two beams. Again, we emphasize that, since we keep the source parameters fixed and observe the time evolution of each geodesic, this setup is different (but essentially equivalent) to the one usually adopted in cosmology, where the observer is fixed and sources with different parameters are considered.
As in [4], we run this configuration on a uniform grid with three different resolutions (corresponding to , , and points per side) in order to estimate the numerical error. All results presented below are convergent to first order, consistently with the convergence order reported for the geometric variables in [4]. All curves represent the Richardson extrapolation, at this order, of the numerical data. The corresponding truncation error (when visible) is indicated by a shaded region around each curve.
4.3 Smallredshift behaviour
For small distances from the source, we expect the photon redshift and luminosity distance to behave respectively like
(4.16)  
(4.17) 
where is related to the first time derivative of the local volume element at the source location:
(4.18) 
(see [26]). Figure 2 shows that this expectation is confirmed by our computation. For large , however, both quantities grow larger than the linear order. Furthermore, the redshift clearly exhibits a nonmonotonic behaviour engendered by the inhomogeneous gravitational field. This is easy to explain as a small, periodic redshift due to the photons climbing a potential hill near the vertices (away from the nearest black holes) and falling into wells near the edge or diagonal midpoints (closer to the black holes). Naturally, the two geodesics are affected in different ways as they trace different paths through the gravitational field.
4.4 Luminosity distance
Due to numerical error, the geodesics deviate from the cell edge and face diagonal during the evolution, but remain quite close to them (the coordinate separation is less than after three cell crossings, in both cases). We can compare the relationship for geodesics and to the same quantity calculated according to four reference models:
All models are fitted according to two prescriptions: the initial scale factor is always set according to
(4.20) 
while the initial expansion rate is set to either (i) the initial time derivative of the proper length of the domain edge (say, the one between and )), which we call a global fit, and is the same procedure as [4]; or (ii) equation (4.18) (which we call a local fit).
Figure 3 shows all the resulting curves. We first recall that the expansion of the BHL, measured by the proper distance of one of its cell edges, could be fitted quite well by an EdS model with the same initial expasion, as shown in [4]. The two models, however, exhibits markedly different optical properties. For geodesic , the relative difference reaches by redshift . This is not surprising: the conditions under which these light rays propagate in a BHL and in an EdS model are substantially different. In the former case, for instance, null geodesics infinitesimally close to or accelerate away from, rather than towards, them.
We notice that the EBA provides the best estimate for in a BHL. We conjecture that this result is due to the fact that this approximation can capture both the largescale geometrical properties of a nonempty universe and the smallscale behaviour of light rays in vacuum. None of the other models satisfies both these conditions. Note also that, for longer times, the EBA works better for the geodesic (along the edge) than for geodesic (along the face diagonal). This is easy to explain if we notice that, because of the 4fold discrete rotational symmetry around the edge, there are no Weyl focusing effects on and therefore the GDE with the Ricci tensor neglected and no Weyl contribution is likely to be a good approximation for the propagation of the neighbouring light rays. On the other hand along the face diagonal we may expect a nonvanishing Weyl lensing around the midpoint area due to the tidal distortion of the rays. Such an effect is not taken into account in the EBA.
4.5 Fitting the FLRW class
It is tempting to consider an FLRW cosmology with the same matter content and initial expansion as the reference EdS, plus an additional stressenergy contribution coming from a cosmological constant, and attempt to tune its value to reproduce the luminosity distance in the BHL.
The left panel of Figure 4 shows a plot of the required at each , for values of in . The right panel shows a cross section of this surface with the planes and , where is the initial proper length of a cell edge. Notice, however, that none of these models would reproduce the expansion history of the BHL spacetime, which follows closely that of a region of an EdS model ( and ) with the same and , as discussed in [4]. This is the core of the fitting problem: the mapping between different properties of an inhomogeneous spacetime to the FLRW class will be different, and in general it will not be possible to identify a single FLRW counterpart capable of reproducing all of the dynamical and optical aspects of an inhomogeneous cosmology.
In Figure 5, we show the constant models which best fit the curves for geodesics and . They are obtained for and , respectively. The relative difference between these models and the exact solution is largest around , where it reaches .
Notice that essentially all quantities discussed so far are affected by oscillations with a substantial initial amplitude, which is subsequently damped. Similarly to the oscillations in the redshift, we conjecture that these features are due to the inhomogeneous gravitational field, and in particular to radiative modes which likely originate in the oversimplified initialdata setup we employed. In a space without an asymptoticallyflat region, it is of course difficult to test (or even formulate) this conjecture rigorously. The compactness of the spatial hypersurfaces, furthermore, means that one cannot simply ignore this initial transient as is customary in, e.g., binaryblackhole simulations, as the waves cannot escape from the domain (although their amplitude is significantly attenuated by the expansion). The presence of this unphysical component of the gravitational field, which we could barely notice in the length scaling we measured [4], affects very prominently, on the other hand, the BHL optical properties, and in particular the photon redshift. Better initialdata constructions which are free from these modes are an interesting field of investigation which goes beyond the purpose of this work.
Finally, it is worth observing that, as mentioned in Section 2, different observers would measure a different luminosity distance on the same spacetime, thereby potentially bringing the BHL result closer to the EdS curve. A boost with respect to the lattice would, for instance, lower the value of the distance, according to equation (2.14). So would a stronger gravitational field, as would be the case if an observer was located closer to the centre of a lattice cell.
4.6 Continuum limit
Finally, it is instructive to study how this behaviour depends on how tightly packed the BHL is, as represented by the quantity introduced in Section 3. For simplicity, here we use the bare mass of the central black hole as an estimate of , and the coordinate size of a cell edge as . In order to keep constant at the value of our base configuration (which had and ), we need to have . As representative masses we choose ; various properties of this BH family are illustrated in Table 1.
0.010  2.15  2.73  0.0046 
0.125  5.00  6.28  0.0250 
0.500  7.94  9.84  0.0630 
1.000  10.00  12.26  0.1000 
5.000  17.10  21.77  0.2924 
We plot the luminosity distance as a function of in Figures 6 and 7. We observe, in particular, that the difference between the luminosity distance in a BHL and in an appropriately fitted EdS does not tend to zero as . The EdS model, therefore, can reproduce the largescale expansion history of a BHL (as illustrated numerically in [3, 4], and deduced analytically in [9]), but is unable to fit its optical properties, even in the limit .
The numerical result is in agreement with the result of the perturbative analysis of Section 3, where we identified differences in the GDE of a BHL with respect to that of an FLRW model. This indicates that cosmologicaldistance estimates of a lumpy spacetime based on a fit with the FLRW class will exhibit a systematic error, regardless of how lumpy the spacetime is. These effects are substantially, but not exhaustively, captured by the EBA, as already observed in the case of other inhomogeneous spacetimes [34, 35, 16].
An important remark is that we observe that the tensor modes discussed in Section 4 intensify as , affecting the smaller BHLs to the point that it becomes impossible to identify a monotonic trend in the luminosity distance for large . For this reason, we are forced to limit our study to very small .
5 Discussion and conclusions
We have investigated the propagation of light along two special curves in the spacetime of a BHL, constructed by numerically integrating Einstein’s equation in 3+1 dimensions. In particular, we have measured the redshift and luminosity distance along these curves, and compared them to the estimates of these observables obtained in suitably fitted homogeneous models and in the EBA. The comparison shows that the latter approximation is the one most capable of reproducing the exact behaviour; we have built a heuristic argument to explain this finding, based on the analysis of the different curvature terms in the GDE. Our finding is congruous with the conclusions of similar studies in other inhomogeneous spacetimes [16]; in our case, however, the models are not backreactionfree by construction, so that we can measure all the relevant contributions to the GDE.
We have also fitted the relationship from the FLRW models with both a constant and a dependent to the data, finding that a value of approximately equal to reproduces the optical properties of the BHL better than the corresponding models with . In other words, in the BHL spacetime the luminosity distance for a redshift is larger than in the corresponding EdS model (the correspondence being based on the same initial proper size and expansion rate). This is also in line with the conclusions of previous studies [18], and arguably equivalent to the finding that fitting alternatively leads to a smaller value for this parameter [15].
Finally, we have examined a family of BHLs with varying BH masses and separations, in order to estimate how our result changes as . In this limit, it was proven in [9] that the expansion history of a BHL tends to that of a flat FLRW model with the same average density. Here, however, we find that the optical properties of a BHL exhibit a finite deviation from the corresponding FLRW model, which reaches by . Given a considerable pollution by tensor modes, which we conjecture originate in our initialdata construction, the luminosity distance is oscillatory, and we are unable to evaluate the continuum limit for larger .
Building a picture of the mechanisms involved in these results, as well as generalizing it to inhomogeneous spacetimes with different matter content and density profiles, is a particularly intriguing but hardtoapproach task. We can start to tackle it by comparing our results to a recent study [50], which also measured the effects of light propagation in an inhomogeneous model which, unlike the ones considered in this work, was filled with dust. In that investigation, percentlevel deviations were detected from the homogeneous Hubble law, which are about an order of magnitude smaller than the deviations reported here. From the arguments presented in this paper, we infer that the discrepancy is largely due to the different representation of the matter filling the two models. The quantitative formulation of this statement is a problem which we reserve for further study.
Acknowledgements
MK and EB would like to thank the Max Planck Institute for Gravitational Physics (Albert Einstein Institute) in Potsdam for hospitality. The work was supported by the project “The role of smallscale inhomogeneities in general relativity and cosmology” (HOMING PLUS/20125/4), realized within the Homing Plus programme of Foundation for Polish Science, cofinanced by the European Union from the Regional Development Fund, and by the project “Digitizing the universe: precision modelling for precision cosmology”, funded by the Italian Ministry of Education, University and Research (MIUR). Some of the computations were performed on the Marconi cluster at CINECA.
Appendix A Geodesic equation 3+1 decomposition
The tangent to a geodesic satisfies
(A.1) 
or, to simplify the following derivation,
(A.2) 
The covariant derivative is expanded in terms of the partial derivative and the Christoffel symbol of the spacetime metric,
(A.3) 
and the LHS is recognised as the derivative along the curve of the component with respect to the curve parameter,
(A.4) 
The Christoffel symbol is expressed in terms of derivatives of the metric,
(A.5) 
and we note that the first and third terms in parentheses are antisymmetric in and , whereas is symmetric, so these terms sum to zero, giving
(A.6) 
as the geodesic equation for the covariant component of the tangent vector.
We now wish to express the RHS in terms of the 3+1 quantities available in a numerical relativity simulation. We summarise briefly the standard 3+1 decomposition of a spacetime (see, e.g. [51]). A foliation of constanttime hypersurfaces is represented by a oneform , which locally can be written as the differential of the coordinate time, . The lapse function is defined as , and the timelike unit hypersurface normal as , so that . The spatial metric on the hypersurface is
(A.7) 
with . A vector is described as spatial if . Since the direction of time evolution, , is not necessarily aligned with the normal to the hypersurface, we express it in terms of a normal component, and the spatial shift vector ,
(A.8) 
To simplify the expressions, we will work in the standard coordinate basis in which the timelike basis vector is , and hence has components . In such a basis, we have , , , , where is any spatial vector, and lower case Latin indices from the middle of the alphabet () indicate spatial components (i.e. ). These relations will simplify the derivation.
We now need to express equation (A.6) in terms of partial derivatives of the spatial quantities available in an NR simulation. We first decompose into a timelike and spatial part,
(A.9) 
where , and aim to find an equation for the evolution of . Note that since