SpinningBlackHoles.jl

SpinningBlackHoles is a pseudospectral solver tailor-made to obtain stationary and axisymmetric black hole solutions in modified theories of gravity.

If you use SpinningBlackHoles you are kindly asked to cite the following paper in any resulting works:

@article{SpinningBlackHoles,
    author = "Fernandes, Pedro G. S. and Mulryne, David J.",
    title = "{A new approach and code for spinning black holes in modified gravity}",
    eprint = "2212.07293",
    archivePrefix = "arXiv",
    primaryClass = "gr-qc",
    doi = "10.1088/1361-6382/ace232",
    journal = "Class. Quant. Grav.",
    volume = "40",
    number = "16",
    pages = "165001",
    year = "2023"
}

The user is also referred to this reference for details concerning the physics and mathematical setup involved.

This project (including all files in this repository) is licensed under the terms of the GNU General Public License Version 3 as published by the Free Software Foundation.

Installation

The first step is to install Julia on your machine, if not already installed. This can be done by following the instructions on the official website.

This package is not in the official registry. You can install it (and its dependencies) by first opening a Julia REPL

julia

and pressing the ] key to enter the Pkg REPL. Afterwards, you must specify the URL to the repository:

pkg> add https://github.com/pgsfernandes/SpinningBlackHoles.jl

and that is it, the package is installed. To exit the Pkg REPL simply press backspace, and to exit Julia you can type exit().

To check the source code and download the examples one can git clone the repository to a suitable folder

git clone https://github.com/pgsfernandes/SpinningBlackHoles.jl

To obtain the field equations for our metric ansatz, you can use whatever computer algebra system you prefer. In our examples we have used Mathematica together with the OGRe package. Instuctions for installing this package can be found at the link. The process of obtaining the field equations is time-consuming (for complicated field equations such as scalar-Gauss-Bonnet it could take about an hour or two), but needs only to be done once for every model one wants to study. To speed up the process, we recommend a quick adaptation of the OGRe package to use the Simplify function of Mathematica by default (rather than the FullSimplify function). This can be achieved by opening the OGRe.m file and replacing the three occurances of "FullSimplify" with "Simplify".

To export the field equations to a format compatible with our numerical infrastructure, we use again Mathematica. Using our exporting file (GenerateSystem.nb) is highly recommended, but not necessary. One can opt e.g. to adapt one of our examples manually or write their own exporter in their favourite language.

How to implement and run a model

We have two working examples of successfully implemented models in the Examples folder: general relativity and scalar-Gauss-Bonnet gravity. The first step in implementing a model is to obtain the field equations for our metric ansatz, which contains four free functions of $r$ and $\theta$: $f$, $g$, $h$, $W$. These field equations can be obtained from their covariant version, e.g. $G_{\mu \nu} = 0$, by computing the relevant geometrical quantities. The process is divided into two steps, both contained in the Deriving FEqs subfolder of each examples' directory. The first step is performed in the Mathematica file Derive_FEqs_Part_1.nb, where we provide our metric ansatz and compute the relevant geometrical quantities and effective stress-energy tensor. The equations of motion are then defined. By comparing the files for the two cases we have provided, it should be clear how different field equations could be implemented instead. For example there is an additional equation in the Gauss Bonnet case associated with the scalar field present. These quantities are then exported by the mathematica file to another file FEqs.m, such that later we can use them by simply importing the file rather than computing them again. Then, in the file Derive_FEqs_Part_2.nb, we import the file FEqs.m and manipulate the field equations by performing relevant linear combinations. Again comparing the files for the GR and Gauss-Bonnet cases is illustrative of how to adapt the code further. Finally, we perform a change of radial variable $x = 1-\frac{2r_H}{r}$ in order to use our spectral method, mapping $r \in [r_H, \infty[$ to $x\in [-1,1]$. Also, we re-label our angular coordinate $\theta \to y$. The final product is stored in the file EqsDiag.mx in the same folder. Note again that these notebooks will have to be adapted depending on the model you are considering.

Now that we have our field equations, we are ready to export them to Julia code by using the Mathematica file GenerateSystem.nb. This file contains a function which receives the field equations we have just calculated as the input, and together with the appropriate boundary conditions which are defined in this file (and could be changed if needed for a new setting), outputs the necessary Julia code (including the Jacobian of the system). The function is:

PrintToJulia[name, dfeqs, bcH, bcI, bcHspin, bcIspin, vars, params, Funcs, rules]

where

• name is the name of the file to which we want to export the field equations. We typically use System.jl
• dfeqs are the field equations calculated in previous steps, contained in the file EqsDiag.mx
• bcH and bcI are respectively the boundary conditions for our fields/functions at the horizon and infinity, assuming that we want to provide the black hole horizon angular velocity as an input parameter. bcHspin and bcIspin are analogous, but for providing the dimensionless spin as the input
• vars are the fields present in our problem, such as our metric functions
• params are the all input parameters that our code uses, such as $r_H$...
• Funcs is related to coupling functions that might exist in our model. For example, in scalar-Gauss-Bonnet gravity we have a coupling function $F(\phi)$, which would be provided in this argument as {{F,p}}.
• rules are the string parsing rules to pass to the function in case the field equations contain some term which was not changed to Julia syntax. Example: if for some reason your field equations contain a term $\sqrt{h}$, it is likely that your generated System.jl file will contain a term Sqrt[h], which is incompatible with Julia syntax. To get rid of this term our rules parameter would have to be {"Sqrt[h]"->"h^(1/2)"}.

If everything was done correctly, you should now have a System.jl file in your folder. You should also have a Solver.jl file in the same folder, which is not generated by Mathematica, and must be adapted depending on the setting. Once again inspection of the Solver.jl files for the two implemented cases should make it clear how they can be adapted further.

Running an example

Once you have a System.jl file and a Solver.jl file adapted to the setting in question, you are now able to run an example. We will demonstrate now, as an example, how to obtain a black hole in the scalar-Gauss-Bonnet model.

The first step is to navigate on a terminal window to the folder containing your System.jl and Solver.jl files. Make sure that the resolution is set as you wish. The resolution can be changed in the beggining of the System.jl file

#DEFINE RESOLUTION IN x
global const Nx=40
#DEFINE RESOLUTION IN y
global const Ny=8

which by default is set to $N_x=40$ and $N_y=8$. Now enter the Julia REPL as before and to load the system run

julia> include("System.jl")

Spinning Black Holes package loaded!
Resolution: Nx=40, Ny=8

and finally, to obtain a black hole solution, run

julia> include("Solver.jl")

Solver initated with:
α/rh^2=1.0
χ=0.6

Iter     f(x) inf-norm    Step 2-norm 
------   --------------   --------------
     0     3.929702e+00              NaN
     1     7.231364e-03     4.183521e-04
     2     1.184027e-06     9.296671e-09
     3     3.623768e-13     3.904442e-16
     4     3.481659e-13     3.279524e-24

Success! Solution converged!
χ=0.6000000000000046
α/M^2=0.1597779841955325
Qs/M=0.0363666087161596
Smarr = 1.624256285026604e-13

The parameters with which you initialize your solver can be changed by altering the OneSolution function contained there

OneSolution(WBC=0.6,rh=1.0,A=1.0,spin=true,tol=1e-10,ToPrint=false,ergosphere=false,petrov=false,light_ring=false,isco=false,sphericity=false,linvel=false)

where WBC is the dimensionless spin of the solution, $\chi$, if spin is true or the horizon angular velocity if spin is false, rh is the event horizon radius that we typically set to 1.0, A is $\alpha/r_H^2$ (with $\alpha$ the Gauss-Bonnet coupling constant). You also have options to print the solutions, and information about the ergosphere, Petrov type, light rings, ISCO, sphericity and linear velocity of the horizon by setting the appropriate arguments to true. Note that if you want to change the resolution you will have to restart the Julia REPL by exiting and entering again. Also, the first time you run the code after entering the REPL will take longer because of Julia's just in time compilation system.