BiGONLight: light propagation with bi-local operators in Numerical Relativity

BiGONLight, Bi-local Geodesic Operators framework for Numerical Light propagation, is a new tool for light propagation in Numerical Relativity. The package implements the Bi-local Geodesic Operators formalism, a new framework for light propagation in General Relativity. With BiGONLight it is possible to extract observables such as angular diameter distance, luminosity distance, magnification as well as new real-time observables like parallax and redshift drift within the same computation. As test-bed for our code we consider two cosmological models, the $\Lambda$CDM and the inhomogeneous Szekeres model, and a simulated dust universe. All our tests show an excellent agreement with known results.


Introduction
Electromagnetic and gravitational radiation are the primary means by which cosmologists and astronomers try to learn about the structure and the evolution of the Universe. All the information we receive from distant objects is inferred from these signals which are affected by the presence of cosmic structures between the source and the observer. These effects are accumulated along the line of sight and they modify the perception of the observer which e.g. can receive the image of the source as if it were in a different position on the sky or with a distorted shape or measure a redshift in the source light spectrum. Actually, most of the information comes from these distortions and they will be measured with unprecedented precision in a wider range of scales and redshift by the next generation of galaxy survey and Cosmic Microwave Background experiments * . On top of that, thanks to the improvements in the modern experimental apparatus, nowadays we will be able of measuring small temporal changes of those effects, the so-called optical drift effects. They may provide important and new information about the Universe structure and evolution, marking the beginning of Realtime Cosmology, [1]. In order to make the most from this revolution in observational astronomy and cosmology, the same accuracy is required in the theoretical predictions and interpretations of what we measure. This tough task requires much effort from two sides. At one hand the most recent progress on cosmological dynamics are represented by general-relativistic simulations of cosmic structures with no assumed symmetries employing exact solutions of Einstein equations, [2,3,4,5,6,7,8], or approximated treatments, [9,10]. The common ambitious aim is to have a description valid from large to small scales which accounts for relativistic effects. On the other hand relativistic effects in the non-linear regime have recently become to be investigated with relativistic simulations in galaxy clustering and lensing observables, and the Hubble diagram, see e.g. [11,12,13,14,15,16,17,18,19] and refs. therein.
From the point of view of the basic theory of light propagation a new approach was presented in [20]. The key ingredients of this new formulation are the Bi-local Geodesic Operators (BGO) which constitute the map from the portion of spacetime occupied by the observer to that occupied by the source and give the full description of the distortion of light rays in between. The main advantage of the BGO formalism is that it provides a unified framework to compute all optical observables, namely those inferred from gravitational lensing effects, like e.g. magnification, shear, and angular diameter distance. The novelty is that source(s) and observer(s) are allowed to move and therefore the observables originated by the variation in time and the motion of sources and observers, like e.g. parallax, and redshift and position drifts, are automatically included, see [20]. In addition, by having the observables written in terms of the BGO, it is easy to disentangle the contributions of the spacetime curvature from those due to observer and source motion and construct specific new probes, e.g. for the curvature as the distance slip investigated in [21] in the ΛCDM cosmology.
In this paper we present BiGONLight * , Bilocal Geodesic Operators framework for Numerical Light propagation, a Mathematica package developed for extracting observables from numerically generated spacetimes using the BGO formalism. The principal aim of the package is to provide a unified procedure to calculate multiple observables in Numerical Relativity: this is guaranteed since, once that the BGO are determined from the output metric of a numerical simulation, a all set of observables are obtained within the same computation. In order to be compatible with the majority of the codes in Numerical Cosmology, BiGONLight encodes the BGO formalism in 3 + 1 form and it uses the Mathematica powerful symbolic algebra manipulation and precision control options.
The paper is organized as follows: in Sec. 2, we give the fundamentals of light propagation and its formulation in terms of Bi-local Geodesic Operators. In Sec. 3 we present BiGONLight and the equations to compute the BGO in 3 + 1 form encoded into the package. The recipe to compute observables in numerical relativity using BiGONLight and the expressions of the observables in terms of BGO are given in Sec. 4. The last section of this paper, Sec. 5, is dedicated to the code tests, which are performed in the following three cosmological models: ΛCDM (see Sec. 5.1), Szekeres (see Sec. 5.2), and a numerically evolved dust universe (see Sec. 5.3). We draw our conclusions in Sec. 6. ..), running from 0 to 7, denote indices for the components of the 8×8 BGO matrix W. A dot denotes derivative with respect to conformal time. Quantities with a subscript 0 are meant to be evaluated at present, whereas the subscript "in" indicates the initial time. Quantities with a subscript S (or O) are meant to be evaluated at the source (observer) position.

Light propagation and the bi-local geodesic operators
Let us start by considering an observer O and a source S separated by a large distance and moving freely along their time-like worldlines. Naming N O and N S the regions of the spacetime in which the observer and the source are moving, we assume that N O and N S are causally connected, so that any signal emitted by S is received by O at any later time. We also assume that the typical length scales of N O and N S are small compared to the distance between them. Within the geometric optics approximation, we can describe the signal as moving along null geodesics γ(λ) connecting the source and the observer, such that γ( observer's and source's positions, respectively. The curve γ is given by the geodesic equation where µ is the tangent vector, λ is the affine parameter spanning the geodesic γ, D Dλ ≡ σ ∇ σ is the covariant derivative along γ. A solution of Eq. (1) is specified by the initial position and the initial tangent vector. In astronomy, we are the only observer performing every measurement, thus it is natural to set the initial conditions (x µ O , µ O ) at the observer location and to trace the geodesic back to the source. Now, if the observer is displaced by δx µ O , a new geodesic connects O and S, and it is characterized by the new initial conditions ( Figure 1: When the observations are repeated over time, the observer and the source may change their positions. This implies that the null geodesic connecting O to S has to change its path accordingly. This situation is described with the deviations of the positions δx µ and of the tangent vectors ∆ µ at the two new extremes S and O. The deviations (δx µ O , ∆ µ O ) can be used to parametrize a family of null geodesics and they propagate according to the Geodesic Deviation Equation (GDE) with initial conditions provided that all the geodesics of the family stay close enough to γ * . In Eq. (3) R µ ν is the short-hand for the optical tidal matrix R µ ν ≡ R µ αβν α β . The standard procedure to solve the GDE (3) was introduced by Sachs in 1961, [22]. The key quantity of this approach is the deformation of the light bundle's cross section, i.e. the projection of the beam on a 2D screen spanned by two orthonormal vectors orthogonal to µ and parallelly transported along γ (the Sachs basis). The GDE is projected onto this basis and then re-written in terms of the deformation matrix, a 2 × 2 matrix containing the expansion and the shear: the GDE recasted this way gives the well-known Sachs equations, see [23] for a complete introduction to the Sachs formalism. Despite its great success, the Sachs formalism describes only momentary observations, so that it can not be used to take into account what happens when the observations are repeated in time and/or the observer and the source move. Obviously, this can be overcome by considering a new fiducial geodesic at some later time and solving again the Sachs equations but this procedure can be very complicated in some cases. In this paper we will use a different approach, which accounts also for these situations in a unified framework, based on the Bi-local Geodesic Operators (BGO), introduced in [20] and summarised below.
Let us start by noticing that the GDE (3) defines a linear, bijective map W(S, O) between its solutions (δx µ S , ∆ µ S ) and initial conditions ( assuming the optical tidal matrix R µ ν to be a smooth tensor field. The bi-tensors W XX , W XL , W LX , W LL acting from O to S are called bi-local geodesic operators. Equation (6) can then be written in the more compact form * Here the GDE is presented as a system of two first-order ODE. The reader may find more familiar the following form of the GDE as a second order ODE D 2 Dλ 2 ξ µ = R µ αβν α β ξ µ , where ξ ν = δx ν and D Dλ ξ ν = ∆ ν . In general, the GDE describes the deviation between any two infinity close null-like, time-like or space-like geodesics.
where W(S, O) is the 8 × 8 Wroński matrix of the GDE acting from O to S. By inserting (7) in the GDE (3), we obtain the propagation equation for the BGO with initial conditions The BGO W is a symplectic mapping, as firstly noticed in [24], since Eq. (8) for null geodesics can be formulated as a Hamiltonian system, and it satisfies the properties with p λ an arbitrary point on the fiducial geodesic γ. The usual set up is specified by giving the initial conditions at O and by integrating the GDE backward in time, up to the source S. The physical motivation is of course that every measurement is done from the observer position. The BGO W(p λ , O) obtained by integrating Eq. (8) backwards connect the observer with an arbitrary point p λ on the geodesic, ending at the source p λ S = S. Although it is natural to study light propagation this way, there are circumstances in which it is more convenient to think forward in time, namely to integrate Eq. (8) from the source to p λ , ending at the observer p λ O = O, and obtain W(p λ , S). To this second case belong all the simulations of cosmological dynamics, in which the Einstein equations are solved forward in time, with initial conditions given at the end of inflation. It follows that forward integration of light propagation, onthe-fly with the simulation of the spacetime dynamics, would be a cost-efficient and time-saving methodology instead of using the natural approach for light propagation in post-processing. These two opposing procedures find a meeting point within the BGO framework, which provides a relatively easy way to transform from forward integrated W(p λ , S) to backward integrated W(p λ , O). The transformation follows from the BGO properties and it is found by multiplying Eq. (11) by W −1 (S, p λ ) from the left and using Eq. (10) to obtain * We will give the explicit transformation rules for all the BGO components in Sec. 3.3.
Let us finally recall that the most important feature of the BGO formalism is that it provides a unified framework to compute optical observables, including also real-time ones, like e.g. parallax, and redshift and position drifts, [20]. In Sec. 4 we will give the explicit expressions of the observables considered in this work in terms of the BGO.

BiGONLight: light propagation with the BGO in Numerical Relativity
The BGO framework finds a natural application in studying light propagation in Numerical Relativity, since it allows to obtain multiple optical observables within the same calculation and this is easily adaptable to perform light propagation onthe-fly with a simulation of relativistic dynamics. For this purpose, we develop BiGONLight (Bi-local Geodesic Operator framework for Numerical Light propagation), a Mathematica package which simulates light propagation in Numerical Relativity within the BGO formalism, described in Sec. 2. It is publicly available at the repository https://github.com/MicGrasso/bigonlight.git and it works as an external library that can be called inside a Mathamatica notebook. Mathematica provides a large variety of numerical methods that can be used and customized to adapt them to the particular problem. This is exploited in the functions SolveGeodesic[], SolveEnergy[], PTransportedFrame[] and SolveBGO[] in which the user can choose the numerical methods used to solve the system of ODE * . Another useful feature is the Mathematica's precision control options, which allows the user to set the precision and accuracy of numerical result through the commands WorkingPrecision, SetPrecision and SetAccuracy. In the following we give a detailed description of the equations for light propagation and the procedure to obtain the BGO implemented in the code. We dedicate a separate section, Sec. 4, to the computation of the observables once the BGO are known.
The required input for BiGONLight are the specetime metric and the source/observer kinematics, which can be provided in two different ways: (i) from the analytical expression of the metric g µν and the emitter/observer four-velocities (u µ S , u µ O ) and four-accelerations (w µ S , w µ O ) for some exact model, and (ii) from the output of a relativistic numerical simulation of the spacetime dynamics. A large variety of numerical codes used in cosmology and astrophysics employs the 3+1 formalism to solve the Einstein equations and simulate full-GR dynamical systems. To be compatible with the numerical output generated by these codes, in BiGONLight we recasted the BGO formalism in the 3+1 form. In the following we summarise some basic definitions of the 3+1 formalism and report the 3+1 version of all the equations for light propagation used to compute the BGO. For comprehensive references of the 3+1 formalism see e.g. [25,26,27,28].
The procedure of splitting a four-dimensional spacetime (M, g µν ) into its 3+1 form, the so-called ADM formalism, was introduced by Arnowitt, Deser and Misner in [29]. It is constructed by foliating a four-dimensional manifold M with a family of 3D space-like hypersurfaces Σ t , labelled by a monotonic function t, such that t = const on each slice. This space-like foliation defines the time-like vector field n µ , which is orthonormal to * The user can choose between three different methods: "RK a 4 th order Runge-Kutta method, "A which lets Mathematica to decide what is the best method to use, and "SS denoting the Stiffness-Switching method, which allows to switch between implicit and explicit methods to resolve stiff problems. the slices and it can be regarded as the four-velocity of Eulerian observes. The geometry on each hypersurface is described by the following quantities: • the induced metric γ µν is defined as the covariant form of the orthogonal projector onto the slices γ µν = g µν + n µ n ν (13) and it is used to measure proper distances on Σ t ; • the covariant derivative on the slice which is written in terms of the 3D Christoffel symbol which represents the curvature of the 3D hypersurfaces with respect to the 4D embedding spacetime.
A natural choice for a coordinate system x µ is the one adapted to the foliation: the corresponding reference frame is such that the three vectors ∂ µ i = ∂ ∂x i µ are tangent to the hypersurface while ∂ µ 0 = ∂ ∂t µ is transverse to it. In particular, the time-like vector field ∂ ∂t µ is tangent to a congruence of world-lines of coordinate observers and it is given by ∂ ∂t where α is the lapse function, which measures the proper time of the Eulerian observers, and β µ is the shift vector, which quantifies the displacement on Σ t of the coordinate observer ∂ ∂t µ with respect to the Eulerian observer n µ . The components of the normal vector and the metric g µν in the adapted coordinates are written with the lapse, the shift and the induced metric as and where the Latin indices runs from 1 to 3. A generic 4D tensor is projected on the slice as BiGONLight is designed to accept as input directly the ADM quantities (α, β i , γ ij , K ij ) generated by a numerical simulation. However, the powerful symbolic algebra manipulation of the Wolfram language allows also to use the analytical form of the metric g µν as input. In this case, the ADM quantities are computed in BiGONLight by the function ADM[] accordingly to Eqs. (15), (17), (18). On top of that, BiGONLight contains other functions to calculate the Christoffel symbols, Christoffel[], and the Riemann tensor, Riemann[]. Note that in both cases, the input metric is provided in form of components within a specific coordinate system and not in full tensorial form. The package is designed to work in any gauge and any coordinate system, leaving the choice to the user. Once we have summarized the basics, we will now describe in more details and step by step the procedure for obtaining the BGO in the 3+1 decomposition: we report the 3+1 version of the geodesic equation following [30] and we derive all the 3+1 ingredients for the evolution equation of the BGO.

Geodesic equation in 3 + 1 decomposition
The first code solving the 3+1 geodesic equation was presented in [31] and it was used to map the event horizon in numerical simulations with black holes (like heads-on collision of two black holes), while a more recent formulation of the 3 + 1 geodesic equation was presented in [30]. Here, we briefly resume the procedure in [30], since we will use their approach throughout all our calculations.
The null geodesics representing light rays connecting source and observer is obtained by solving where µ is the tangent vector which obeys to the null condition It is 3 + 1 decomposed as: where E is defined by E = −n µ µ and V µ n µ = 0. In other words, EV µ is the component of µ tangent to Σ t and En ν is the orthogonal one. Substituting Eq. (22) in the geodesic equation (20) after a long but straightforward calculation, see [30], the geodesic equation decouples in two differential equations for the orthogonal and the tangent components, respectively. In the BiGONLight package the two functions EnergyEquations[] and GeodesicEquations[] give the differential equations (23) and (24), respectively. Next, we need to specify the initial conditions. The function InitialConditions[] is specifically constructed to get the initial conditions for V i and E consistent with the null condition (21) and the decomposition (22). The ODE (23) and (24)

Parallel transport equation in the 3 + 1 decomposition
The BGO map the changes of the deviations (δx µ , ∆ µ ) between the observer O and the source S along the photon geodesic γ, see Eq. (7). This is possible only if we introduce a frame parallel transported along γ, which allows us to compare quantities at the observation point and at source position. For our purposes, the definition and the parallel transport of the frame have to be split into the 3+1 form. Let us start with the parallel transport of a generic vector e µ along a given geodesic with tangent µ , which is governed by α ∇ α e µ = 0 .
The 3+1 decomposition of e µ reads where we have defined the orthogonal component Cn µ = −n α e α n µ and the tangent component E µ = γ µ α e α . The parallel transport equation (25) becomes Now, we make use of some 3 + 1 well-known relations: for the first term we have n µ ∇ µ C = 1 α L α n C, in the second bracket we substitute n α ∇ α n µ = D µ log α and the definition of the extrinsic curvature ∇ α n µ = −K µ α , and for the expansion of the last two terms we use the two identities: We then re-write Eq. (27) split into two parts, the one proportional to n µ and the other tangent to Σ t as where both must vanish individually in order to satisfy the parallel transport condition (25). The last steps consist in expanding the Lie derivative and converting partial derivatives with respect to the time t into total derivative via The final result for the 3 + 1 parallel transport equation is where we have only spatial indices i, j = 1, 2, 3 since all the quantities lie on Σ t . The system of equations (31) has to be solved for each vector of the frame we choose. Following [20], we choose the semi-null frame (SNF) composed by the tetrad of vectors where u µ is the matter four-velocity, µ is the tangent of the photon geodesic. The two vectors φ µ A are orthonormal to both u µ and µ , namely: with Q a real number. Each vector of the SNF is then decomposed in 3 + 1 form as and parallelly propagated by solving Eq. (31) for its orthogonal and tangent components, i.e. Φ α n µ = −n σ φ σ α n µ and F µ α = γ µ σ φ σ α . In BiGONLight this is demanded to the function PTransportedFrame[], which gives as output the components in (34) (24) and Eq. (31) we presented are valid for any type of geodesics and to parallelly transport any type of vectors along that geodesic. From now on, we will restrict the BGO formalism to the case of a SNF parallelly propagated along a null-like geodesic. The equation for the BGO has then to be projected onto the SNF. The result is simply given by where the only formal difference with respect to Eq. (8) is that the covariant derivative along the photon geodesic D/Dλ reduces to the total derivative d/dλ in the SNF. In Eq. (35), the optical tidal matrix is projected in the SNF and it is given by where the inverse of the induced metric of the frame is used to raise the indices. In general R µ ν is a 4×4 matrix with non trivial components. However, in the SNF it is easy to use the symmetries of the Riemann tensor to show that the components R 0 ν = R ν and R µ 0 = R µ vanish. Now, let us use Eq. (22) and Eq. (34) in Eq. (36) to write the optical tidal matrix in terms of 3 + 1 quantities After some tedious but straightforward calculations, we finally obtain where are the Gauss relation, the Codazzi relation and the Ricci relation, respectively (see e.g. [27] with initial conditions: The systems in Eq. (44) are solved separately in BiGONLight by SolveBGO[]. Getting the BGO with the procedure just described is one important part of BiGONLight. Let us recall that the only inputs required are the spacetime metric, the observer four-velocity components and the initial and ending points.
With initial conditions in Eq. (45), Eq. (44) gives the BGO W(p λ , O) integrated backward in time from the observer. The relation with the BGO W(p λ , S) integrated forward in time from the source can be explicitly written down for all the BGO components from Eq. (12). To this end we need the inverse matrix W −1 which is easily found from the symplectic property where Ω is the 8 × 8 non-singular, skew-symmetric matrix with Latin tilded indices running from 0 to 7 (ã, · · · = 0, 1, . . . , 7) and the Greek bold indices (α = 0, 1, . . . , 3) indicate the components in the SNF. By inverting Eq. (46) we find The transformation from forward to backward BGO in Eq. (12) finally reads These relations are coded in the section "Forward to backward transformation for W operators" of each sample notebook in the repository https://github.com/MicGrasso/ bigonlight.git.

Optical observables with BiGONLight
The recipe to obtain optical observables using BiGONLight can be summarized as follows. One needs to: (i) specify the spacetime metric g µν and the source S and observer O kinematics, namely four-velocity u µ and four-acceleration w µ . They can be given already in 3 + 1 components or as 4D quantities and the functions ADM[] and Vsplit[] will do the splitting of g µν , and u µ and w µ respectively; (ii) set up the initial photon geodesic using Vsplit[] for the null tangent µ , which gives its 3 + 1 components E and V i . Alternatively one can give E and a vector V i that has to be assigned by specifying the spatial direction V 2 and V 3 and use InitialConditions[] to get V 1 from the null condition; (iii) obtain the geodesic equations Eq. (23) and Eq. (24)  (iv) set up the initial conditions for the SNF according to Eq. (33), directly in 3 + 1 components or using SNF[], which is specifically designed to compute the SNF in 3 + 1. Then PTransportedFrame[] will give the SNF parallel transported along the light ray; * For the purposes of our paper, we need to solve only photon geodesics. However, the code can be used to trace any type of geodesics, namely time-like and space-like also, by specifying the appropriate initial tangent vector in the 3 + 1 splitting with Vsplit[].
(v) compute separately R µ ν projected into the SNF with OpticalTidalMatrix[]; (vi) obtain the ODE system for the BGO in Eq. (44)  Note that the components of the optical tidal matrix can have a very complicated expression, that may cause problems in solving the GDE (44). This can be overcame by using an interpolated form of R µ ν : the Mathematica Interpolation[] function allows to use different methods and reach an excellent precision, see Sec.5.
Step (vi) is the starting point to compute the optical observables, which are all given by different combinations and/or functions of the W components. It is important to remark that all the observables in the BGO formalism are written in terms of W(S, O), namely the map computed from the observer to the source. As we already recalled, this is obtained directly by integrating the GDE backward in time. However in some cases, e.g. if the spacetime model comes from a numerical simulation, it may be more convenient to get the inverse BGO map W −1 (S, O) and then use the transformations Eqs. (49)- (52) to obtain the W needed for the observables.
Here we list the four observables that we study in this paper. The redshift is simply given by its definition where σ is the tangent to the light ray, and u σ O and u σ S are the observer and source four-velocities. The same definition in 3 + 1 splitting reads: The angular diameter distance is formally given by where W XL A B is the map between the physical size of the source and the angle subtended in the sky, as measured at the observer position, namely Conversely, the parallax distance is related to the displacement of the observer position and the apparent angular shift of the source position, as measured from the observer The expression for the parallax distance is The last observable that we consider in this paper is the redshift drift ζ, given in terms of the BGO by, [32] ζ In the above expression τ O is the proper time of the observer, the first term represents the Doppler effect caused by the four-acceleration w σ of the observer and the source, and U is an 8 × 8 matrix given by the following combinations of the BGO Let us notice that, even if the BGO formalism is independent of the specific choice of the frame used, the observables are dependent on this choice, as evident by the explicit dependence on u µ O , u µ S , w µ O and w µ S in Eqs. (55)-(60). Indeed, it is possible to transform locally between two different frames using an appropriate Lorentz transformation, but this modifies the observables introducing special relativistic effects like Doppler effect or aberration.
The reader can find the derivation of the Eqs. (55)- (59) in [20,33,32]). All the expressions in Eqs. (55), (58) and (59) are new with respect to the standard approach, in the sense that these observables are expressed within a new, unified framework. However, while for D ang and D par there already exist analogous formulas, where instead of the BGO we have the magnification and the parallax matrix (see [21] for the comparison), it did not exist a general formula for the redshift drift: Eq. (59) looks the same for every spacetime model under consideration. Instead, in the standard approach the redshift drift is calculated by taking the derivative with respect to the time coordinate of the definition of the redshift, and this depends on the specific form of the metric tensor and null-geodesic, the latter depending in turns on the symmetries that one gives to the initial conditions. This means that the equations to get the redshift drift in the standard approach look different for each specific model and/or configuration of the light rays * .

Code tests
In this section we test the accuracy of BiGONLight within well-known cosmological models. The tests are performed by considering the following observables: redshift, * To be more precise, Eq. (68) is valid for the FLRW model only as well as Eq. (96) is valid for the Szekeres model and geodesics along the symmetry axis only. Instead, Eq. (59) looks the same in both cases. angular diameter distance, parallax distance, and redshift drift. We compare the results obtained with three different procedures, by defining the estimator ∆O(BGO, X) where O BGO refers to Eqs. (53), (55), (58), and (59). We consider the following three cases: (i) the ΛCDM model, where the specetime metric is the analytical input for BiGONLight to compute O BGO and for O X we use the analytical well-known solutions for all the four observables, see Sec. 5.1; (ii) the inhomogeneous Szekeres model [34], where the specetime metric is again the analytical input for BiGONLight but to obtain O X we solve numerically a specific differential equation for each observable, see Sec. 5.2; (iii) the Einstein-de Sitter model, where the input for BiGONLight are the 3 + 1 quantities coming from the Einstein Toolkit simulation and O X is obtained analytically, see Sec. 5.3. It worth noting that if O X is obtained analytically, then max |∆O(BGO, X)| represents the simulation error: in case (i) we have just the computational error from BiGONLight whereas in case (iii) the final error in the observables is the combined effect of both the Einstein Toolkit and BiGONLight finite precision. On the other hand, if O X is obtained numerically, as in case (ii), then ∆O(BGO, X) gives only an estimation of the accuracy of the two methods used.

The ΛCDM model
The first group of tests regards the study of light propagation in the flat ΛCDM model. This is an exact solution of the Einstein field equations representing an homogeneous and isotropic spacetime and the matter-energy content consists of irrotational dust of cold dark matter and a cosmological constant Λ. The line element is given by where η is the conformal time and a(η) is the scale factor, which is the solution of Einstein equations and describes the dynamics of the model. The explicit result is found to be [35] a(η) = 3 Ωm 0 where cn(y|r ) is the Jacobi elliptic cosine function, with y = 4 √ 4 . To test BiGONLight we consider two classical observables, namely the redshift z and the angular diameter distance D ang , and two interesting observables that are not yet measured in the cosmological context as they belong to the new research field named Real-time Cosmology, see Ref. [1]. One is the parallax distance which exploits the motion of the Solar System with respect to the Cosmic Microwave Background frame providing a baseline of 78 AU per year for the cosmic parallax * . Cosmic parallax was first proposed in 1986 in Ref. [37] and it is expected to be measured by the Gaia satellite, [38], in the next few years. For discussions and forecasts about the measurements of the cosmological parallax distance we refer to Refs. [39,40,41,1,36,42,43,44] and Refs. therein. The other is the redshift drift, i.e. the time variation of the redshift of a source. It was first derived for the FLRW models in Refs. [45,46]. Since then, and particularly in recent years, a lot of work has been done to investigate the measurability of the redshift drift in cosmology and the information gained, see e.g. Refs. [47,48,49,50,51].
The analytical expressions in the flat FLRW cosmologies for the four observables that we consider are where the Hubble parameter in the ΛCDM model is given by We normalize the today scale factor a 0 = 1 and we take H 0 = 67.36 km s −1 Mpc −1 , the matter parameter today Ω m 0 = 0.315, and the cosmological constant parameter Ω Λ = 0.685 from [52]. The integral in Eq. (66) and (67) can be solved analytically and the results is, [35] where F ξ(z)|r the elliptic integral of the first kind, with arguments r = 2+ (71) * For an exhaustive definition of the parallax see e.g. Ref. [36], in which the author distinguishes between the three different cases: one source observed by two observers separated by spacelike interval (classic parallax), two sources observed by two observers separated by spacelike interval and two sources observed by one observer at two different moments (cosmic parallax also known as position drift).
The fact that we have analytical expressions for the observables makes this model a perfect test-bed for the code. We compare the results from BiGONLight with the one in Eq. (65)-(68) by considering the variation As shown in Fig. 2, the numerical implementation of the BGO method is in excellent agreement with the analytical formulas in ΛCDM, the variation ∆O being 10 −22 ÷ 10 −32 .
The maximum value of ∆O, of the order of 10 −22 , represents the numerical error over the observables and we reached such small values by using the precision control options WorkingPrecision, PrecisionGoal and AccuracyGoal implemented in Mathematica.

The Szekeres model
The second group of code tests is performed considering an inhomogeneous cosmological model which is an exact solution of Einstein equations, firstly presented in [53]. The line element for the Szekeres spacetime is: with α and β functions of the spacetime coordinates (t, x 1 , x 2 , x 3 ) that are determined by solving the Einstein equations. We can distinguish two different families of Szekeres models depending whether ∂ x 3 β = 0 or ∂ x 3 β = 0: the first case defines the "class I" family, which is a generalization of LTB models (with non-concentric shells), while the case ∂ x 3 β = 0 corresponds to a simultaneous generalization of the Friedmann and Kantowski-Sachs models and it defines the "class II" family. For a cosmological formulation of the Szekeres spacetimes see e.g. [54]. For our tests, we will use a class II Szekeres model filled with dust and a positive cosmological constant as presented in [34]. For this model, the line element (73) is rewritten as: and, for the special case of axial symmetry around x 3 , we have the following form for the function Z: The function D is the growing mode solution of the first-order Newtonian density contrast defined as δ = ρ − ρ ΛCDM ρ ΛCDM and it is given by, see e.g. [55], with 2 F 1 (a, b, c, y) being the Gaussian (or ordinary) hypergeometric function. The term B in Eq. (75) is constant and given by (see App. C in [56]) where D in = a in for initial conditions set deeply in the matter-dominated era. The function β + specifies the spatial distribution of the first-order density contrast and it can be related to the peculiar gravitational potential φ 0 via the cosmological Poisson equation at present time For the tests, we will use a sinusoidal profile for the peculiar gravitational potential and amplitude A such that δ 0 = 0.1 for the density contrast today.
In the following, we will present the tests over the redshift, the angular diameter distance and the redshift drift. Contrary to the ΛCDM case, here the observables are obtained using two numerical methods and the difference is expressed by All the three tests are done considering that the observer O is located at the origin of the reference frame, with coordinates (η 0 , 0, 0, 0), and she/he sees the light coming from a comoving distant source S, with coordinates (η, 0, 0, x 3 ). The light beam is propagating along the x 3 axis, with tangent vector µ = ( 0 , 0, 0, 3 ).
The first observable we test is the redshift: for a photon travelling along the x 3 -axis it is where 0 is obtained by solving the geodesic equation In the above expressions Z is given in Eq. (75), a is the scale factor (64), and H the Hubble parameter (69) of the ΛCDM background. The variation ∆z(BGO, Sz) refers to the numerical solution of Eq. (81) as opposite to the 3+1 geodesic solved by BiGONLight.
The second observable under analysis is the angular diameter distance D ang . The standard procedure to compute D ang is solving the Sachs focusing equation simply vanishes, as shown in [57]. Using the Einstein equations to express R µν µ ν in terms of the density contrast and using conformal time instead of the affine parameter, the focusing equation assumes the form where the dots refers to derivatives respect to conformal time η and the initial conditions at the observation point are . In synchronouscomoving gauge the density contrast along the geodesic comes directly from the continuity equation and reads In our estimator Eq. (79) the angular diameter distance D Sz ang is obtained by integrating Eq. (83), whereas D BGO ang is obtained from Eq. (55) with BiGONLight. The last test concerns the calculation of the redshift drift, namely the secular variation of the redshift of the source. It was calculated for some inhomogeneous cosmological models, see e.g. [58,59,60,61,62], but to our knowledge there is no expression for the Szekeres model considered here, thus we give in the following a short derivation. Let us consider that during the proper time lapse δτ O the spacetime coordinates of the observer change from while the same quantities at the observer position Our final aim is to compute the redshift drift, namely Let us begin with the variation of the source redshift with respect to the observer proper time δz δτ O , where it is better to obtain first dδz dx 3 (Eq. (91)) and dδτ dx 3 (Eq. (94)) separately, and then combine them to get an ODE (Eq. (96)), whose solution gives the redshift drift. The starting point are the differentials dη dx 3 and dz dx 3 . The first is simply given by the null condition 0 = −Z 3 , and reads where we used Eq. (88), Eq. (81), and the fact that the redshift at the observer is fixed. Now we use Eq. (88) to differentiate Eq. (86) and similarly for the redshift we have where we keep first-order terms in δη and δz only. Rearranging terms and using Eq. (89) again, we get d dx 3 The final step is to express the variation of the conformal time δη in terms of the proper time at the observer δτ O . We use their relation, which is simply since δx i = 0. We need the derivative with respect to x 3 , which reads and, together with Eq. (89) and Eq. (93), the solution of the last equality is For our test, in the estimator Eq. (79), ζ Sz is the numerical solution of Eq. (96) and ζ BGO is the expression in Eq. (59) obtained with BiGONLight.
As for the ΛCDM model, also for the Szekeres model we have a very good agreement between the observables calculated using the standard procedure and the observables from BiGONLight. The smallness of all the variations ∆z(BGO, Sz), ∆ζ(BGO, Sz) and ∆D ang (BGO, Sz) shown in Fig. 3 means that our code could be a reliable tool for light propagation in inhomogeneous cosmologies, represented here with the Szekeres model, which is computationally more complicated than the homogeneous ΛCDM case.

A dust universe in Numerical Relativity
The main application of the BiGONLight package is the computation of observables from numerically simulated spacetimes. For our test we choose to use the FLRWSolver * , [5], which is a module (or thorn) of the Einstein Toolkit (ET), [2]: the ET is a collection of open-source codes, called thorns, based on the Cactus framework, [63], which allows to solve the Einstein equations in the BSSN formulation of the 3+1 splitting, [64,65]. The role of the FLRWSolver is to provide the initial conditions in the form of linear perturbations around the Einstein-de Sitter (EdS) background, which are then evolved with the ET. Here, we limit ourself to the EdS background model and set perturbations to zero. In other words, we consider a FLRW model in which the Universe is flat and contains only cold dark matter. The line element of the EdS model in conformal time is where a EdS (η) = η 2 is the scale factor. We carry out the simulation in a cubic domain −L ≤ {x 1 , x 2 , x 3 } ≤ L with periodic boundary conditions and spatial resolution ∆x = ∆y = ∆z = L 20 , where L is the simulation unit length in Mpc * . The initial data are given at η in = L and such that γ in ij = δ ij . The simulation runs with the ET up to η 0 = 33.2L, which corresponds to integrating from redshift z = 1100 to present time z = 0, and we choose a fixed temporal resolution ∆η = L 100 due to computational time convenience. To give an estimation of the simulation error we define which is the variation between the analytical scale factor in EdS, i.e. a EdS = η 2 and the scale factor from the numerical simulation a ET = det(γ ij ) where O BGO is computed numerically with BiGONLight using as input the EdS model simulated with the ET and O EdS is the analytical expression in the EdS model that  results are shown in Fig. 5. What we see is the error on ∆O(BGO, EdS) which has in principle two contributions: one from the simulation of the EdS model and the other from the simulation of light propagation. We have already isolated the second contribution coming from BiGONLight in the ΛCDM test. Indeed, since we use the analytical solution for ΛCDM both in the numerical and analytical computation for the observables, the error ∆O(BGO, ΛCDM) we find in Fig. 2 is entirely due to the simulation of light propagation and is of the order of 10 −22 ÷ 10 −31 . On the other hand, in Fig. 4 we see that the accuracy of the ET simulation we use is much larger, i.e. of the order of 10 −10 , which is of the same order of the one for ∆O(BGO, EdS) in Fig. 5. Therefore we can conclude that the final error on the observables we find in Fig. 5 is settled by the accuracy of the ET in simulating the EdS model. Let us finally remark that for this test we perform light propagation using forward integration, namely from the source S to the observer O. We also repeated the computation of ∆O(BGO, EdS) using backward integration, form O to S, in BiGONLight and we found the same results.

Conclusions
In this paper we present and test BiGONLight, a Mathematica package for relativistic light propagation in Numerical Relativity. Our new code implements the BGO formalism, a new approach to geometric optics in General Relativity, firstly introduced in Ref. [20], which we recasted here in the 3+1 form to be compatible with the structure of relativistic numerical simulations. The generality of the BGO framework allows BiGONLight to be suitable to study light propagation in different contexts: from small scales, e.g. inside our Galaxy, to the ultra-large scales of Cosmology, as long as light propagation can be treated within the assumptions of geometric optics. BiGONLight is extremely flexible also from the technical point of view, requiring as input only the spacetime metric and the observer and source kinematics: the user can choose any gauge and also the type of input to use, namely numerical from a relativistic simulation or analytical from an exact solution of the Einstein field equations. An analytical solution for the spacetime metric in any perturbative scheme can be also given as input but the code does not apply the perturbative scheme to light propagation. This is indeed the approach we adopted in Ref. [56]. We plan to implement perturbation theory in a future version of BiGONLight. A key feature of the BGO framework is that all optical effects are encoded in the bi-ejective W map, which can be written form the observer O to the source S or viceversa. In this paper we give the ready-to-use transformation relations in Eqs. (49)- (52). This property in BiGONLight means that the user is free to adopt different approaches to trace light propagation. On one hand, one can trace the photons from the observer to the source using the conventional backward integration, on the other hand one can choose to do the other way around and use forward integration. This second method is particularly convenient in Numerical Relativity because it allows us to perform light propagation on-the-fly with a simulation of spacetime dynamics, which employs integration forward in time by construction. We decide to test BiGONLight in the computation of four cosmological observables: the redshift, the angular diameter distance, the parallax distance, and the redshift drift, the last two being real-time observables that are new in the cosmological context. We perform three different kinds of tests. In the first we test the accuracy of the code against the analytical results in the ΛCDM model and we find an excellent agreement of the order of 10 −23 ÷ 10 −31 , see Fig. 2. We are able to reach this extremely good accuracy by using the precision control options and the many well-tested numerical methods to solve ODE implemented in Mathematica. For the second test we consider the Szekeres inhomogeneous cosmological model as presented in Refs. [57,56] and we compare the BGO formalism with the standard procedures used for the angular diameter distance and the redshift drift. Let us point out here that the BGO formalism provides a unified way to compute multiple observables within the same framework, contrary to the standard approach. In particular, for the angular diameter distance we compare the result obtained from the BGO, Eq. (55), with the one from the Sachs equation, Eq. (83). For the redshift drift we compare the BGO general formula, Eq (59), with the result of the standard calculation, Eq. (96), that we derive here for the specific set up for light propagation in the Szekeres model considered. Also in this case we find a very good agreement of the order of 10 −22 , see Fig. 3. For our final test we consider the EdS background simulated with the Einstein Toolkit together with the FLRWSolver. The numerical output of this simulation is used to compute the observables with the BGO in BiGONLight which are then compared with the usual analytical formulas. Our findings are shown in Fig. 5 where we see that the accuracy on the observables calculated with BiGONLight is ruled by the accuracy of the Einstein Toolkit simulation, which is limited to the order of 10 −10 for the specifications that we use in our paper. Contrary to the first two tests, here the BGO are calculated by integrating the GDE forward in time and then we use the BGO property in Eq. (12) to obtain the observables, see Sec. 4. The current version of BiGONLight is designed to perform light propagation in post-processing and not on-the-fly, with the advantage of being able to accept input from different numerical codes for cosmological dynamics. Within this choice, every physical quantity Q phys can be expressed as Q phys = Q comp L α , where Q comp is dimensionless, L is a arbitrary length to be chosen, and α is a certain exponent. This way of writing is particularly useful in numerical simulations, where all physical quantities are represented as dimensionless numbers and units are assigned when analysing the results. Usually L is fixed to be a length meaningful for the specific physical situation under consideration. For example, a common choice in numerical cosmological dynamics is to set L equal to a characteristic length of the simulation (e.g. N-body and GR hydrodynamics), such as the side of the simulated box. In cosmology it is usually chosen to set L equal to the conformal time in Mpc, i.e. L = η phys , thus η comp = 1, at some special instant like e.g. at initial time or today. This choice is particularly convenient, since conformal time is found by integrating the Friedmann equation and reads