BiGONLight: light propagation with bilocal operators in numerical relativity​
​by Michele Grasso(Center for Theoretical Physics, PAS).
Email: grasso@cft.edu.pl
Source: https://github.com/MicGrasso/bigonlight

Summary:

BiGONLight is a Mathematica package that encodes the Bi-local Geodesic Operators formalism (BGO) to study light propagation in Numerical Relativity. With BiGONLight it is possible to compute 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. The user just needs to specify the metric
g
μν
and the source S and observer O four-velocities and four-accelerations: in the case that the metric is obtained from a numerical simulation, one can use as input the ADM quantities, i.e. the spatial metric
γ
ij
, the extrinsic curvature
K
ij
and the gauge functions lapse α and shift
i
β
obtained from the simulation. However, the code also can accept 4D metrics given in the analytical form, i.e. with the components of the metric as arbitrarily complicated functions of the spacetime coordinates, and use the powerful Mathematica symbolic algebra manipulation to perform the 3 + 1 splitting of the metric
g
μν
. The package contains a collection of functions to compute geodesics, parallel transported frames, and to solve the BGO equation, as well as other functions to calculate useful quantities in differential geometry like the Christoffel symbols and the Riemann tensor.​An example of how to compute observables with BiGONLight is described in the following.

Info on the notebook:

This notebook uses the BiGONLight package to compute observables in a numerically evolved dust universe (EdS model). The spacetime dynamics is obtained using the Einstein Toolkit (ET), see Frank Löffler et. al. (2012) Class. Quantum Grav. 29 115001, a community-driven collection of codes for relativistic simulations, together with the FLRWSolver thorn, see Hayley J. Macpherson et. al. (2017) Phys. Rev. D 95, 064028.The notebook is structured as follows: - model set-up: import the 3+1 quantities lapse α, shift
i
β
, spatial metric
γ
ij
, and extrinsic curvature
K
ij
from the numerical simulation,- geodesic: provide initial position and tangent vector
{
μ
x
,
μ
ℓ
}
as input to find null geodesics;- parallel transported frame: give the initial conditions for a frame
μ
ϕ
α
that will be parallel transported along the computed geodesic (we use the SNF, i.e.
μ
ϕ
α
={
μ
u
,
μ
ϕ
1
,
μ
ϕ
2
,
μ
ℓ
}
with
μ
u
the four-velocity and
μ
ϕ
A
two orthonormal vectors)- optical tidal matrix: compute the optical tidal matrix
μ
R
ρσν
ρ
ℓ
σ
ℓ
in the parallel propagated frame; - BGO: compute and solve the Geodesic Deviation Equation for the Bi-local Geodesic Operators along the null geodesic;- observables: compute the optical observables in the BGO framework and compare with known results. The observables we will compute are: redshift, redshift drift, angular diameter distance, luminosity distance, parallax distance and distance slip.​Notation: Greek indices (α, β, ...) run from 0 to 3, while Latin indices (i, j, ...) run from 1 to 3 and refer to spatial coordinates only. Latin indices (A, B, ...) run from 1 to 2. Tensors and bitensors expressed in a semi-null frame SNF are denoted using boldface indices: Greek boldface indices (α,β,...) run from 0 to 3, Latin boldface indices (a,b,...) run from 1 to 3 and capital Latin boldface indices (A, B, ...) run from 1 to 2.
Set up the current directory as working directory and load the package.
In[]:=
​​SetDirectory[NotebookDirectory[]];​​Needs["bigonlight`"]
In[]:=
prec=16​​$PreRead=(#/.s_String/;StringMatchQ[s,NumberString]&&Precision@ToExpression@sMachinePrecisions<>"`16."&);(*MachinePrecisions<>"`prec."&*)
Out[]=
16
a generic direction can be written as
μ

=
0

∂
0
+γ(sinχcosξ
∂
1
+sinχsinξ
∂
1
+cosχ
∂
3
)=
μ
0

(
∂
0
)
+γ
μ
N
with respect to the local tetrad
∂
μ
. In particular for a null direction
0

=±1
and
γ=1
, while for a timelike direction
0

=±cγ
and
ψ=cγβ
where
β=
vel
c
the normalized velocity and
γ=
1
1-
2
β
the Lorentz factor.
In[]:=
[Vt_,lorentz_,direction1_,direction2_]:=​​Module[{V0=Vt,γ=lorentz,χ=direction1,ξ=direction2},​​V0{1,0,0,0}+γ(Sin[χ]Cos[ξ]{0,1,0,0}+Sin[χ]Sin[ξ]{0,0,1,0}+Cos[χ]{0,0,0,1})​​]​​​​StandardPlotStyle[size_,sizelegend_,ylabel_,xlabel_,title_,legend_,position_]:=Module[{ft=size,ftL=sizelegend,yL=ylabel,xL=xlabel,pt=title,list=legend,pl=position},​​If[list{},{FrameLabel{{yL,""},{xL,""}},FrameTrue,FrameTicksTrue,FrameStyleDirective[FontSizeft],GridLinesAutomatic,GridLinesStyleDirective[Gray,Opacity[0.2],Dashed],PlotLabelStyle[pt,Bold,FontSizeft]},{FrameLabel{{yL,""},{xL,""}},FrameTrue,FrameTicksTrue,FrameStyleDirective[FontSizeft],GridLinesAutomatic,GridLinesStyleDirective[Gray,Opacity[0.2],Dashed],PlotLabelStyle[pt,Bold,FontSizeft],PlotLegendsPlaced[LineLegend[Table[Style[list[[m]],FontSizeftL],{m,1,Length[list]}]],pl]},Print["Error"]]];

FLRWSolver

Import data from spacetime simulation

In[]:=
​​iℋ=SetPrecision[Take[Cases[Import["Data_FLRWSolver/res_0.05/H.average.asc","Table"],{__?NumericQ}],All,-2],prec];
In[]:=
​​Do[If[iℋ[[i,2]]>
-8
10
,rel=i-1;Print[ToString[{i,iℋ[[i,1]],iℋ[[i,2]]}]];Break[],rel=i-1],{i,1,Length@iℋ}]
In[]:=
iα=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/alp.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];
In[]:=
​​iγxx=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gxx.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​iγxy=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gxy.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​iγxz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gxz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​iγyy=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gyy.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​iγyz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gyz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​iγzz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/gzz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];
In[]:=
ikxx=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kxx.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​ikxy=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kxy.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​ikxz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kxz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​ikyy=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kyy.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​ikyz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kyz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];​​ikzz=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/kzz.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];
In[]:=
​​iρ=SetPrecision[Take[Take[Cases[Import["Data_FLRWSolver/res_0.05/rho.average.asc","Table"],{__?NumericQ}],rel],All,-2],prec];
In this notebook we use geometric units G=c=1. This allows us to express mass, length and time as the same unit. For our convenience, we will express all quantities as dimension-less quantities in units of length L. Therefore, any physical quantity is expressed as
phys
Q
=
comp
Q
α
L
, where α is a certain exponent and
comp
Q
is dimensionless. L is arbitrary chosen: in the ET simulation we set up
-L≤{x,y,z}≤L
, with
Δx=Δy=Δz=0.05L
and
Δη=0.01L
. Within this set up, from Friedmann equation
2
ℋ
=
2

a
a
=
2
ℋ
0
a
a
0
η=
2
ℋ
0
1
2
a
a
0
, setting
η
in
=L
and using the values of
a
from the simulation, i.e.
a
in
=1and
1
2
a
0
=33.2
, we have
L=
2
ℋ
0
-
1
2
a
0
=268.11Mpc
,where we used
ℋ
0
=
67.36
299792.458
1
Mpc
. This value for L gives that the age of the EdS universe is
η
0
=8901.2Mpc
as aspected.From the expression for the redshift in FLRW,
z+1=
a

a
we have that the simulation runs from
z
in
=
a
0
a
in
-1=1101.2
to
z
0
=0
.

Units and Parameters

L=
2
1
2
iα[[Length@iα,2]]
1
H0
/.H0
67.36
299792.458
(*Mpc*);​​ρin=
3
2π(*G*)
(*
1
2
L
*);​​​​(*
ℋ
in
=

in
L
=
2
L
->*)in=2;​​(*
ℋ
0
=

0
L


0
=
ℋ
0
L=
2
√
a
0
*)0=
in
1
2
iα[[Length@iα,2]]
;​​

Model set-up

Geodesic

Initial conditions

Solve the equations numerically, given the functions in the metric explicitly

Test that the constraint conditions are satisfied for all the evolution

Plot the null geodesic

Parallel transport of a frame

set-up the initial conditions for the SNF

Parallel transport the SNF

Define the induced metric of the SNF “h”

Optical Tidal Matrix

Use the interpolation to make a better OPT

 operators

Forward to backwards transformation for  operators

​

Observables

Compute the redshift

Angular diameter distance

​
The comparison gives:

​
Luminosity distance and Parallax distance

The analytical expression of the parallax distance in EdS is

​
The distance slip μ

It is a new observable derived for the first time here, Phys. Rev. D 99, 064038 (2019), and calculated for the FLRW cosmologies here, Phys. Rev. D 101, 063506 (2020). Its peculiarity lies in the fact that its value is independent of the momentary motions of both the source and the observer and that for short distances it shows a tomographic property, being proportional to the amount of matter along the line of sight.

Redshift drift effects

In EdS, the analytical expression for the redshift drift is

Note:

BiGONLight can be used also to do light propagation on astrophysical scales, provided that the geometric optics approximation is valid. More examples can be found on the package’s repository https://github.com/MicGrasso/bigonlight

Additional information:

More information about the BGO framework can be found in Phys.Rev.D 99,064038 (2019) (arXiv:1811.10284) or Phys.Rev.D 101,063506 (2020) (arXiv:1912.04988). More information about the package can be found in Grasso, Villa (2021) (arXiv:2107.06306) and in Grasso, Villa, Korzynski, Matarrese (2021) (arXiv:2105.04552).