A new approach and code for spinning black holes in modified gravity

We discuss and implement a spectral method approach to computing stationary and axisymmetric black hole solutions and their properties in modified theories of gravity. The resulting code is written in the Julia language and is transparent and easily adapted to new settings. We test the code on both general relativity and on Einstein-Scalar-Gauss-Bonnet gravity. It is accurate and fast, converging on a spinning solution in these theories with tiny errors ($\sim \mathcal{O}\left(10^{-13}\right)$ in most cases) in a matter of seconds.


Introduction
In the last decade, with the observation of gravitational wave events by the LIGO Scientific Collaboration [1][2][3][4][5][6], and interferometry measurements of the centre of M87 and the Milky Way by the Event Horizon Telescope Collaboration [7][8][9], we have entered a new era of testing gravity, probing the nature of black holes and Einstein's theory of general relativity (GR) in the previously inaccessible strong field regime.
In GR, mathematical theorems guarantee that in (electro-)vacuum the gravitational field of stationary black holes is described uniquely by the Kerr(-Newman) metric [10]. As eloquently put by Subrahmanijan Chandrasekhar, the uniqueness theorems along with a set of other results dubbed no-hair theorems (see [11] for a review) assert that the Kerr metric provides "the absolute exact representation of untold numbers of massive black holes that populate the universe". While all strong regime observations are so far compatible with this "Kerr hypothesis", any eventual deviation would provide a much sought after smoking-gun for new physics.
Modification of the field equations describing gravity, however, naturally leads to an increase in their complexity such that analytic analysis becomes intractable. With closed-form solutions not available, one is forced to resort either to perturbation theory or numerical methods. In the strong-field regime, perturbative approximations may not be well-justified, leaving numerical studies as the most promising way forward. In this arena, the ever-increasing precision of our observations and measurements necessitates increasingly accurate solutions.
In this paper, we will describe a numerical method and code capable of solving with high accuracy a system of non-linear elliptic partial differential equations (PDEs), such as those that appear when analyzing stationary and axially symmetric spacetimes, and implement this in a publicly available code. A first version of our numerical implementation is available in the GitHub repository in Ref. [44]. The code is written in Julia and can be run with ease on laptop-class computers, with solutions being found in a matter of seconds. The Julia language is fast, memory efficient, and easy to manipulate, ensuring that implementing different modified gravity theories is not a difficult task.
Our code follows similar previous numerical solvers for these spacetimes, in particular the non-publicly-available FIDISOL/CADSOL solver [45][46][47] (which has been extensively used in the literature, see e.g. [13-15, 21, 23, 24, 26, 28, 48]) and a recent publicly available solver developed in Refs. [49,50]. We have several motivations for writing another code. First, in this work we show that pseudospectral methods * are ideally suited to solving the type of equations at hand. In our tailor-made implementation we therefore make use of such methods †, while both former codes utilise finite difference methods. In contrast to the first code mentioned above, our implementation is also open source, and moreover in our bench-marking we find our code to be far more accurate as detailed further below. The code of Refs. [49,50] is also significantly more accurate than that of the FIDISOL/CADSOL solver (though the documented accuracy is still less than our own when bench-marked on the Kerr solution) and is publicly available. This code is, however, written in C, and our use of Julia leads to simple code that can easily be adapted. Our code is also considerably faster. Our overall aim is a publicly available, accurate, well documented code that is transparent and easy to use code. Furthermore, our code provides a toolbox to explore several properties of the obtained black hole solutions, rather than being only a PDE solver.
This paper is organised as follows. In section 2 we introduce the reader to pseudospectral methods and the technical machinery that will be necessary to apply them in the context of black hole physics. Next, in section 3 we will describe how we can use the aforementioned methods to solve the stationary and axisymmetric field equations for gravity, discussing the boundary conditions, coordinate compactifications, and our numerical approach. We further discuss many of the properties that can be extracted from a spinning black hole solution. Finally, in section 4 we start by validating our methods and code against the Kerr black hole, which is known in closed form, and later use our machinery to obtain stationary and axisymmetric black holes in Einsteinscalar-Gauss-Bonnet gravity for linear and exponential couplings. We also discuss the accuracy of our code, and further compare with results from other codes in published literature. We work with units such that G = c = 1.

Spectral Methods
The idea behind spectral methods is to approximate a smooth solution to a system of differential or integral equations by a sum over a finite number of basis functions. In this section, we review how this works. Given that our aim is a clear and adaptable code, the presentation is relatively complete, and summarises that given in John P. Boyd's book on spectral methods [53], to which the reader can turn for full details (see also Ref. [51]).
For simplicity, we begin with the one dimensional case and illustrate how the method finds an approximation to the smooth solution to a differential equation, u(x), with the differential equation written in the form where R is called the residual of the system. The solution, u(x), can be approximated by a finite truncated series solution u N (x) such that where {φ n (x)} ∞ n=0 is a set of global and orthogonal basis functions, {α n } ∞ n=0 is the set of spectral coefficients, and N is the resolution. In this setup, u N (x) can be said to be a numerical solution of the system (1) if spectral coefficients are found such that the residual is below a certain prescribed tolerance. The method is therefore global rather than local, with an exponential convergence with N for problems with smooth solutions. Since black hole solutions are smooth, we expect exponential convergence when come to find such solutions using spectral methods. This is in contrast to the polynomial convergence rate of most other numerical methods, such as finite element or finite difference schemes. Furthermore, numerical solutions obtained via a spectral method provide an analytical approximation to the problem at hand (rather than a set of approximate numerical values at a discrete number of points).
As noted the basis functions must be orthogonal, which implies that where the brackets represent the inner product of two functions f (x) and g(x) with respect to the weight function, ω(x) > 0, on the interval [a, b] as The set of basis functions used should have a number of further properties: i) they should be easy to compute (e.g. trigonometric functions or polynomials); ii) the approximations built out of the basis functions should converge rapidly to the true solution as the resolution is increased; iii) they should be complete, which means that any solution can be represented to arbitrarily high accuracy by taking the resolution to be sufficiently high. Two commonly used sets of basis functions that obey these requirements are sines and cosines, as used in a Fourier series, and a special class of polynomials dubbed Chebyshev polynomials.

Chebyshev Polynomials
For non-periodic problems, Chebyshev polynomials are the most natural choice as the spectral series is guaranteed to converge exponentially fast (provided our domain is restricted to the interval x ∈ [−1, 1]). The nth Chebyshev polynomial (of the first kind) is defined as T n (x) = cos (nθ) , θ = arccos x, or equivalently by the three-term recurrence relation and hence form an orthogonal basis. Their derivatives are given by where U n (x) denotes the nth Chebyshev polynomial of the second kind, defined by the recurrence relation and with derivative Note that some derivatives require special care at the boundaries x = ±1, and must be computed as a well-defined limit, namely

Interpolation
Interpolation is the process by which a function is approximated by a finite sum of suitable basis functions. The idea is that the sum is constructed such that the aproximation agrees with the true function at the chosen set of interpolation points (also called collocation points). The objective is that the interpolant provides a good approximation to the true function also between those points. By virtue of the minimal amplitude theorem [53], Chebyshev polynomials are widely used in interpolations. The reason is twofold. First, when using the so called Chebyshev nodes (or Gauss-Chebyshev points), x n , as collocation points, the effect of the Runge phenomenon (numerical instabilities near the boundaries in the form of uncontrolled oscillations) is minimized. These points are the roots of the N th Chebyshev polynomial, and are given by Secondly, when Chebyshev polynomials are used as the basis for the interpolation, the interpolation error is distributed uniformly over the whole range. The algorithm to interpolate a smooth function u(x) using a truncated Chebyshev series written as * relies on finding the optimal spectral coefficients {α n }, and uses the discrete orthogonality relation of Chebyshev polynomials: where the x j are the points given in Eq. (12). These discrete relations imply that We present in Fig. 2 an illustrative example of a Chebyshev interpolation performed for several resolutions using the above expressions.

Trigonometric functions
For periodic problems, sines and cosines are the most suitable basis functions for a spectral series. These obey well known orthogonality relations, and form the basis for the Fourier series representation of a periodic function. As we will see, a finite sum of these functions can be used to generate a trigonometric interpolation to a periodic function. Moreover, we can often simplify further by taking into account symmetries. For example, considering the core problem considered in this paper, we note that stationary and axisymmetric black holes are solutions to a system of two-dimensional elliptic PDEs that depend on the radial coordinate and the zenith angle θ ∈ [0, π]. These solutions also often possess definite parity with respect to θ = π/2 (i.e. in most cases they are symmetric about θ = π/2), and therefore we need only to consider the range θ ∈ [0, π/2]. In this range, the following discrete orthogonality relations hold where θ n = (2n + 1)π 4N , n = 0, . . . , N − 1. Table 1 summarizes the parity properties of the functions appearing in the relations above, and together with these orthogonality relations we see that a function, u(θ), symmetric about θ = 0, π/2 can be interpolated using only even cosines such that with the spectral coefficients Fourier series Parity w.r.t. θ = 0 Parity w.r. Even Odd Even = 0 = 0 = 0 = 0 Table 1. Properties of the elements of a Fourier series of a function u(θ), depending on the parity symmetries, along with a scheme of its boundary values. Here, n ∈ N 0 . The entries on this table for θ = π would be equivalent to those of θ = 0.

Solving an ODE with a spectral method -a first example
So far we have seen how a known function can be approximated by a finite sum of suitable basis functions using interpolation. Now we turn to the problem of how to find such an approximation to an unknown function that is the solution to a given differential equation.
To understand how to solve differential equations using a spectral method, we will first consider a simple ordinary differential equation (ODE) example. Consider the one dimensional non-linear boundary value problem We will find an approximate solution to this boundary value problem in the form of a Chebyshev spectral series, and later compare our results with the known exact solution, given by To illustrate the calculations analytically, we will first consider a (very) low resolution approximate solution with N = 3, where Here there are N = 3 unknowns (α 0 , α 1 , and α 2 ), and N BC = 2 boundary conditions. Once we substitute our ansatz of Eq. (22) onto the residual given in Eq. (20) we obtain together with the boundary conditions To find the approximate solution we simply need to determine values for the three unknowns. Given that we have only N = 3 degrees of freedom, and the two boundary conditions provide two constraints, we need only one further equation to find the values. To get this constraint the idea is to evaluate the residual at the N −N BC = 1 collocation point given by Eq. (12) (with the N in that expression given by 1), which gives the point x = 0. With our resolution of N = 3, finding an approximate solution to the boundary value problem then reduces to solving three non-linear coupled algebraic equations for the spectral coefficients, given by the two boundary conditions of Eq. (24) together with the residual of Eq. (23) evaluated at x = 0. The solution to the system is By construction this is an interpolation to the exact solution (21). We note that had we chosen the resolution N = 4, the number of collocation points where we would have to evaluate the residual would be N − N BC = 2, and would be given by Eq. (12) as x = ± sin π/8 ≈ ±0.382683. This would give N = 4 coupled equations to find the four unknowns in this case. This then generalises to arbitrary N .
On Fig. 3 we plot the exact solution against the approximation obtained with N = 3, together with the absolute error, |1 − u N /u|, whose maximum can be seen to be O (10 −2 ) already for a very low resolution N = 3. We also plot the approximation to the solution of the boundary value problem, but for a resolution N = 24, where we observe that errors become of order machine precision (O (10 −16 )). In Fig. 4 we plot the behaviour of the maximum absolute error as a function of the resolution, where exponential convergence is observed. As a rule of thumb, the truncation error is typically the same order-of-magnitude as the last coefficient retained in the truncation series. An important point to make here is that even though for N = 3 the approximation system has a closed-form analytical solution for the spectral coefficients, once higher resolutions are considered, a numerical root-finding method (such as Newton-Raphson) has to be employed. To successfully employ a Newton-Raphson method, a good initial guess for the spectral coefficients is of the utmost importance. We will come back to this in the next section.
To conclude, spectral collocation methods, also known as pseudospectral methods, are powerful tools that can be used to find high accuracy numerical solutions to differential equations. They provide global analytical approximations for the solution, and handling any kind of boundary condition is straightforward.

Root-finding methods -Newton-Raphson
To numerically solve the system of algebraic equations for the spectral coefficients a root-finding method must, in general, be employed. In particular, we will utilise the well known Newton-Raphson method. In the one-dimensional case, the method attempts to solve the equation f (x) = 0 iteratively, starting with a initial guess, x 0 . Successive values of x are then generated until a value, x * , is reached at which the equation is approximately solved to a certain prescribed tolerance. The series of iterations takes the form For example, assume we want to find the root of the function f (x) = x 3 + x − 1, known to be x * ≈ 0.6823278 to eight decimal places. Starting with x 0 = 1 as our initial guess, using Eq. (25) we obtain thus converging to x * in four iterations to the prescribed tolerance of eight decimal places. Convergence, is however, not guaranteed, and particularly in more complicated settings an appropriate choice of starting point is extremely important, and must be chosen carefully.
The generalization of the method to N variables with N equations finds the root of a vector-valued function F : R N → R N , and amounts to solving the linear system at each iteration for the unknown x n+1 − x n , where J is the N × N Jacobian matrix of the system, defined as Constructing the Jacobian matrix of a given system is not always an easy task, but is relatively straightforward for the system of equations that arises when using the spectral method to solve an ODE. Another advantage of this method. As described, in this case, the system to be solved, F , will be composed of the residual R evaluated at the Gauss-Chebyshev points (12), and the boundary conditions, and will in general involve u, u x and u xx . Our unknowns are the spectral coefficients α j . Thus, to facilitate the computation of the Jacobian, we may use the chain rule and only then substitute the spectral expansions for the function u. This process is easily generalizable to a system of differential ODEs/PDEs (rather than a single one) * . As previously stated, when using a Newton-Raphson method, the choice of initial guess to the spectral coefficients is extremely important, because a non-appropriate choice will likely result in non-convergence of the algorithm. A good initial guess can sometimes be difficult to obtain, especially when dealing with systems of PDEs, where the number of coefficients is large (for our specific black holes problem, typically of O (10 3 ) coefficients). A good way of tackling this issue stems from a good understanding of the problem in question. For example, from an effective field theory point of view, a Kerr black hole is probably a good approximation to a black hole solution in modified theories. Therefore, since we have a closed-form expression for a Kerr black hole, an interpolation of this solution can used to generate an initial guess for the spectral coefficients in modified theories. * One must be careful when labelling the spectral coefficients of the different functions as it might be a source of errors.

Black Holes -Metric ansatz, The Kerr Solution, Boundary Conditions, and Connection with the Numerical Approach
We will now apply the methods described in the previous sections to black hole physics and approximately solve the coupled PDEs that arise when obtaining stationary solutions in a given theory of gravity. We will focus on a particular ansatz for the black hole spacetime written in quasiisotropic coordinates with line-element which is stationary, axisymmetric, and circular. Here f , g, h and W are dimensionless functions of the radial and angular coordinates r and θ, and where r H is the (coordinate) location of the event horizon. The spatial coordinates range over the intervals In order for the line-element to be a solution to the theory of gravity at hand, the functions, f , g, h and W must satisfy a set of PDEs that result from the field equations of the theory. The spacetime presented possesses two Killing vector fields, k = ∂ t and Φ = ∂ ϕ , and the linear combination where Ω H is the angular velocity of the horizon (to be defined below), is orthogonal to and null on the event horizon. This Lewis-Papapetrou form for the metric is motivated by the discussion of Ref. [54], which asserts that the above metric ansatz is consistent for a generic theory of gravity provided that its solutions can be obtained perturbatively from a solution in the general relativity limit. Note that our form of the metric functions on the line element of Eq. (29) differ somewhat from the standard form used in other works (see e.g. [13][14][15]21,23,24,26,28,48]). The reasons for this will become clearer once we make a connection to our numerical approach, and are related to numerical accuracy issues.

General Relativity -The Kerr Black Hole
To begin, let us consider the known Kerr black hole, which is the solution to the stationary and axisymmetric field equations of GR in vacuum. For completeness, we present its charged generalization, the Kerr-Newman solution of electrovacuum in Appendix A. The Kerr black hole solves the field equations where G µν is the Einstein tensor, which follow from the Einstein-Hilbert action where R is the Ricci scalar of the metric g µν . With the ansatz of Eq. (29) the Kerr black hole solution reads and M is the ADM mass of the black hole. The total angular momentum per unit mass, a, of the solution is related to M and r H via where we have defined the dimensionless spin The mass M and total angular momentum J can be read off from the metric components as r → ∞, where leading to Note that the Kerr black hole in the quasi-isotropic coordinate system presented in Eq. (29) can be obtained from the standard textbook Boyer-Lindquist coordinates solution with the radial coordinate transformation The inverse transformation is given by

Boundary Conditions
To solve the set of PDEs that result from the field equations in a particular theory of gravity, suitable boundary conditions should be imposed. These are obvious if an exact solution, such as the Kerr solution, is known by a trivial examination of the metric functions. However, in more intricate cases in modified gravity lacking an exact solution the boundary conditions must be found with a careful examination of the field equations and employing suitable expansions of the involved functions near the domain boundaries. For example if theories possess a GR limit when some parameter tends to zero, an expansion about the Kerr solution is possible. With this process, we find that in all cases to be discussed in this work within modified gravity theories, the metric functions must obey the same boundary conditions as the Kerr solution does. These conditions are summarized next.
(i) Axis boundary conditions: Axial symmetry and regularity of the solutions on the symmetry axis θ = 0, π, imply the following boundary conditions Moreover, the absence of conical singularities further imposes that on the symmetry axis h = 1, for θ = 0, π.
All solutions to be discussed in this work are are also symmetric with respect to a reflection on the equatorial plane θ = π/2. Therefore, as was discussed above, it is enough to consider the range θ ∈ [0, π/2] and one of the boundary conditions becomes (ii) Event horizon boundary conditions: The black hole solutions discussed here possess an event horizon located at a surface with constant radial variable r = r H . The boundary conditions that the metric functions f , g and h obey at r = r H are The reason for the Robin-type boundary conditions that the functions f and g obey comes from the inclusion of the N 2 factor in front of f in the coefficient that multiplies dt 2 in the metric ansatz, Eq. (29). This factor is chosen such that these functions do not contain a double-zero in a near-horizon expansion, allowing for more accurate solutions in this region, and therefore, a more accurate extraction of horizon physical quantities such as the area and temperature of the event horizon. We find that there are (at least) two possibilities for the condition that the function W should obey at the horizon, one of which must be chosen appropriately such that the number of input parameters is kept where Ω H is a constant interpreted as the angular velocity of the event horizon, which in the case of a Kerr black hole is given by (iii) Asymptotic boundary conditions: Requiring asymptotic flatness (i.e., that as r → ∞, our solution approaches the Minkowski spacetime), the functions f , g, and h obey lim Similarly to the boundary conditions at the event horizon, we find (at least) two suitable conditions for the function W lim or, from the asymptotic expansion of Eq. (39)

Connection with the numerical approach
To recap, the field equations of a gravitational theory once applied to the line element of Eq. (29) will result in a set of non-linear coupled elliptic PDEs in r and θ subject to the boundary conditions described above. Our objective is to solve this system of PDEs numerically using a spectral method. For this we introduce the compactified radial coordinate mapping the range r ∈ [r H , ∞[ to With the compactified coordinate, the radial boundary conditions change, and we proceed to give the new conditions next.
Event horizon boundary conditions: The boundary conditions that the metric functions f , g and h now obey are for x = −1. For the function W , the first possibility (Eq. (46)) remains unchanged (W | x=1 = r H Ω H ), whereas the second becomes Asymptotic boundary conditions: The asymptotic boundary conditions the functions f , g, and h are now Asymptotically, function W now obeys either or at x = 1. With our compactified radial coordinate, and given the symmetries of our problem * , a suitable spectral expansion for the black hole metric functions (collectively denoted by F = {f, g, h, W }) is given by where N x and N θ are the resolutions in the radial and angular coordinates. Note that, as discussed above, the angular boundary conditions are automatically satisfied by this expansion (c.f. Table 1).
As mentioned previously, we will usually use the Kerr metric itself to set our initial guess when working with modified theories of gravity, and to do so we will need the expression for the spectral coefficients that follow from an interpolation of a two-dimensional function u(x, θ), which is given by where x k and θ l are given in Eqs. (12) and (17) respectively. Each Kerr black hole is uniquely described by two input parameters. For example, in the presentation given in Eq. (34), these are the location of the event horizon r H and the ADM mass M . We have seen, however, in expressions (36) and (48) that they are related to the dimensionless spin χ and the horizon angular velocity Ω H . Therefore, using the correct parametrization, the Kerr solution can be described by any input pair chosen from r H , Ω H , χ, and M . In the numerical approach, in a theory agnostic setting, one input parameter that must be used is r H because it enters directly the metric ansatz and the definition of our compactified coordinate x. We have, however, freedom in the choice of the other input parameter in the numerics. To the best of our knowledge, so far in the literature for similar problems [13-15, 21, 23, 24, 26, 28, 48], the other input parameter has always been chosen as the event horizon angular velocity Ω H . Using this input pair (r H , Ω H ), we find compatibility with the boundary conditions for the function W if we choose Eqs. (46) and (57) at the horizon and infinity, respectively. Then, in the case of a Kerr black hole, one finds that for a fixed value of r H , two branches of solutions exist, as shown in Fig. 5. This follows from inverting the relation (48). The first branch of solutions starts at a vanishing value of Ω H (for fixed r H ) and exists until at which point Then, a second branch appears, and Ω H tends backwards towards zero.
As Ω H → 0 on this second branch, extremal solutions are approached. The existence of two branches of solutions is not unique to Kerr, and is observed as well in the modified theories of gravity to be discussed in this work. We note that the numerical procedure gets rather difficult as near-extremal solutions are approached, as our metric ansatz with the described boundary conditions is not compatible with extremal solutions. A novel approach that we can also adopt is to choose the pair (r H , χ) as the input pair. This input pair is compatible with the W boundary conditions of Eqs. (55) and (58) while maintaining the number of input parameters at two. We often find it very convenient to use the dimensionless spin as an input parameter, for example when exploring domains of existence, or simply when working on a single solution where a certain χ is wanted. Our numerical spectral method is not only powerful because high accuracy solutions are produced, but also because highly non-linear boundary conditions can be handled with ease (which is the case of the boundary condition of Eq. (58)).
To solve the system of field equations subject to the discussed boundary conditions we must construct a suitable grid. This is done as follows. We assume a resolution N x × N θ . The discrete grid points in the x direction are chosen according to Eq. (12), where we take N = N x − 2, together with the boundary points x = −1 and x = 1, such that the total number of points in the x direction is N x * . In θ, our points are chosen as in Eq. (17), where we take N = N θ . The x and θ points together form the schematically shown in Fig. (6), in blue. Assuming there are a total number, N f uncs , of functions to solve for, there are N f uncs × N x × N θ degrees of freedom (spectral coefficients) in the problem, as seen in the spectral expansion of Eq. (59). For each value of θ in the grid at the x boundaries we impose for each function the horizon and asymptotic boundary conditions as discussed before. This gives us a total of N f uncs × 2 × N θ equations (Fig.  6, in red). The remaining N f uncs × (N x − 2) × N θ equations come from imposing the N f uncs residuals resulting from the field equations at each non-boundary x value, for each θ. The number of degrees of freedom is then equal to the number of equations to solve, as it should. A small caveat -the absence of conical singularities imposes that Eq. (43) must be obeyed (i.e. for our coordinate range, h = 1 at θ = 0). While we could leave this condition outside the numerical scheme and use it as another test to the code, we find that imposing it allows obtaining solutions with (much) higher accuracy. In our particular implementation, therefore, we have swapped the evaluation of one of the residuals at θ = 0 (for all interior values of x) * with the condition of Eq. (43), see Fig. 6 in yellow.

Numerical Approach: A summary
Here we summarize our numerical approach for clarity. To solve the field equations, some preliminary work must be done. First, we employ the metric ansatz of Eq. (29) which contains four unknown functions, f , g, h, and W . Plugging this metric ansatz onto the field equations of the theory, leads to a set of non-linear coupled PDEs that depend on the functions and their first and second derivatives (F, ∂ r F, ∂ 2 r F, ∂ θ F, ∂ 2 θ F, ∂ rθ F). The set of field equations is then expressed in terms of the compactified coordinate x defined in Eq. (52) and put in residual form (i.e., R (x, θ, ∂F) = 0). The same is done for the appropriate boundary conditions as discussed. This part of the process is usually done resorting to a computer algebra system such as Mathematica, Maple or SageMath †. Our code, which can be found at [44], includes detailed examples demonstrating how to derive the elliptic field equations for the two theories we will discuss: General Relativity and Einstein-Scalar-Gauss-Bonnet gravity. These examples are implemented using Mathematica. They serve as a valuable reference and can be easily adapted to different contexts. Due to their complexity, these elliptic equations can consist of hundreds or even thousands of independent terms, hence we won't present them here. The residuals (and appropriate Jacobian) are then exported to a Julia coding file in order to solve the problem using the developed numerical infrastructure. Each function is expanded in a spectral series given by Eq. (59) and the input parameters are then specified (depending on the chosen boundary conditions for the function W ). To successfully solve the field equations, a good initial guess must be provided to our Newton solver. For this, we interpolate the functions of the known Kerr solution using Eq. (60), obtaining appropriate spectral coefficients to be provided as a good initial guess. If new fields are present, as is the case with modified theories, we typically take advantage of perturbative solutions and interpolate them as a guess. Convergence is assumed once the norm difference between the spectral coefficients of two successive iterations is less than a certain prescribed tolerance.
To speed up the solver, the values of our basis functions and their first and second derivatives are calculated at all the grid points and stored, such that no repeated evaluations are performed. Another optimization that we found particularly impactful was to store the values on the grid of the trigonometric functions that typically appear in the residuals, sin θ and cos θ.
Once a solution is obtained, physical quantities can be extracted from it as we discuss in the next section, and the solution can be used for numerous investigations.

Physical Properties of Stationary and Axisymmetric Black Holes
Once a numerical stationary and axisymmetric black hole solution has been found using our code, we can extract important quantities of physical relevance. In this section, we review many of the quantities that one can extract from a solution, some of which can be used to test the accuracy of our code. We have implemented additional code to extract all these quantities from a numerical solution.

Quantities of interest
Starting with the asymptotic quantities, we have seen that the mass M and angular momentum J can be extracted from the asymptotic expansion of Eq. (38) or Eq. (39). In terms of the coordinate x these are given by We remark that such a simple expression for the extraction of J is the reason why we have defined the function W in this way -such that its decay is of the form ∼ 1/r, allowing for more accurate results. In a circular spacetime, the zeroth law of black hole mechanics holds, which means that the surface gravity is constant on the horizon of the stationary black hole. The surface gravity is defined as κ 2 = −1/2(∇ µ ξ ν )(∇ µ ξ ν ), where ξ was defined in Eq. (31). The Hawking temperature [56] can then be obtained from the surface gravity as The induced metric on the horizon is and from it we can compute several quantities of interest, the most important being the event horizon area Also of importance is the entropy, which is given in the Iyer-Wald formalism by [57] where µν is the binormal vector to the event horizon surface. In the case of a Kerr black hole the above expression reduces to the simple form S = A H /4. The horizon and asymptotic quantities are connected via the Smarr type relation [57][58][59][60] The Smarr relation is extremely important when studying numerical solutions as it provides a test to the code that relates physical quantities obtained on the horizon and asymptotic regions, allowing us to estimate the accuracy of the numerical method. Also of interest is the perimetral radius R which is a geometrically significant radial coordinate such that a circumference along the equatorial plane has perimeter 2πR. It is related to the coordinate r by To explore the horizon geometry, it is useful to define the horizon circumference along the equator and along the poles With these two quantities, we can define the sphericity For a Kerr black hole s ≥ 1, with s increasing with spin. That means that spin deforms the horizon towards oblateness. The linear velocity of the horizon quantifies how fast the null geodesic generators of the horizon spin relative to a static observer at infinity, and is given by where EllipticE denotes the complete elliptic integral of the second kind, and we also note that 2r H /M = 1 − χ 2 . The Kerr solution is Ricci flat, and thus the Lagrangian of GR vanishes on-shell and therefore so does the last term in Eq. (68).

Ergoregion
The ergoregion is defined as the domain outside the event horizon wherein the norm of the asymptotically timelike Killing vector k = ∂ t becomes positive, g µν k µ k ν > 0. It is bounded by the event horizon and by the surface where Within the ergoregion, an object cannot appear stationary with respect to a distant observer due to the intense frame-dragging. * Furthermore, ergoregions raise the possibility of extracting energy from a black hole via the Penrose process, or superradiant scattering [61]. Starting from the well-known result for the ergosphere of a Kerr black hole in Boyer-Lindquist coordinates and inverting the relation of Eq. (40) we obtain that in quasi-isotropic coordinates the ergosphere of a Kerr black hole is located at where the subscript "E" refers to "ergoregion". Due to the symmetries of our problem, we need only consider the range θ ∈ [0, π/2]. To visualize the ergoregion, we introduce the coordinates In Fig. 7 we observe the ergoregion of a Kerr black holes in the X − Z plane for several values of dimensionless spin.

Petrov type
The Petrov classification allows for a kinematic characterization of the gravitational field in a coordinate independent manner using algebraic properties of the Weyl tensor C µναβ , namely its number of distinct principal null directions. This classification is useful, for example, when searching for exact solutions, or for a Carter-like constant [62]. Using the Newman-Penrose formalism, the information is contained in five complex scalars known as the Weyl scalars. With the null tetrad {l µ , n µ , m µ , m µ }, where l µ and n µ are real, and m µ , m µ are complex conjugate satisfying the orthonormality conditions l µ n µ = 1, m µ m µ = −1 and all other products zero, the Weyl scalars are defined as With the above scalars, the following Lorentz invariant quantities can be constructed Given the above quantities, it is possible to determine the Petrov type of a given spacetime. The classification is summarized in Table 2 [63]. In particular, a spacetime is said to be algebraically special if D = 0. The Kerr(-Newman) spacetime is Petrov type D. In a numerical setup, we also find useful to introduce the speciality index defined as [64] With an appropriate choice of tetrad, following Ref. [64], it is possible to gauge away ψ 1 and ψ 3 to zero. Such a tetrad would be for example

Marginal Stable Circular Orbits: Light Rings and ISCO
The study of marginal stable circular orbits is highly relevant for the observational properties of black holes. The innermost stable circular orbit (ISCO) of massive particles is the smallest possible radius for a stable circular orbit and is often taken to mark the inner edge of an accretion disk around a black hole. Accelerated charged particles orbiting the black hole emit synchroton radiation whose physical properties are connected with the frequency of geodesics at the ISCO. Therefore, physical properties of an astrophysical black hole can be inferred via measurements of the ISCO through accretion disks.
Light rings are circular null geodesics, typically unstable, allowing light to encircle a black hole before being scattered to infinity or falling into the event horizon. From an observational point of view, they are important for observations made with the Event Horizon Telescope as they are intimately connected with the shadow of the black hole [65].
The two independent killing vectors of the spacetime, k µ = ∂ t and Φ µ = ∂ ϕ , have the associated conserved reduced energy E and angular momentum L where≡d/dλ. The above expressions can be rearranged in terms ofṫ andφ Considering orbits restricted to the equatorial plane, θ = π/2, the condition associated with the normalization of the four-velocity of the particles becomes with = {0, 1, −1} for a massless, massive and tachyon particle, respectively. We disregard = −1 from now on. Substituting the expressions of Eq. (84) in the above condition, and solving forṙ 2 , we can define the effective potential The conditions for a circular orbit areṙ = 0 andr = 0, from which follows that at the location of orbit. The dash denotes a derivative with respect to r. These conditions can further be rearranged into algebraic equations that must be satisfied simultaneously Light Rings For a light particle, = 0. In this case, calculations are simpler than in the massive case. Solving the first equation for L in (89) and substituting in the second we obtain which is to be evaluated on a radius r. The smallest root of the above equation is the location of the light ring. In Boyer-Lindquist coordinates the location of the circular photon orbits of a Kerr black hole are given by [66] r LR± BL = 2M 1 + cos where the plus sign refers to co-rotating photons, and the minus sign to counter-rotating photons. In quasi-isotropic coordinates the location of the circular photon orbits can be obtained using the inverse transformation in Eq. (41). ISCO For a massive particle, = 1. The ISCO is located at a saddle point of the effective potential, such that the condition U ef f = 0 should be imposed. This is equivalent to imposing E 2 g ϕϕ + 2ELg tϕ + L 2 g tt − g 2 tϕ − g tt g ϕϕ = 0, in addition to Eq. (89). To find the location of the ISCO, we first solve Eq. (89) for E and L as functions of the metric functions and their first derivatives, and later substitute these onto Eq. (92). Similarly to the light-ring case, we obtain a second order equation to be solved for r, the smallest root of which corresponds to the location of the ISCO. In Boyer-Lindquist coordinates the location of the circular massive particle orbits of a Kerr black hole are given by [66] r ISCO± where , and the plus sign refers to co-rotating particles, and the minus sign to counter-rotating particles. In quasi-isotropic coordinates the location of the circular orbits can be obtained using the inverse transformation in Eq. (41).

Orbital frequencies at the ISCO and Light Ring
The orbital angular frequency of particles both at the ISCO and light ring is given by where the above expression is to be evaluated at the location of the ISCO/light ring, ω + is the angular frequency of co-rotating particles and ω − is the angular frequency of counter-rotating particles. In the case of a Kerr black hole we have at the light ring, and at the ISCO. The orbital frequency at the ISCO is associated with the cut-off frequency of the emitted synchrotron radiation generated from accelerated charges in accretion disks, and the angular frequency at the light ring is related to the time-scale of the response of the black hole when it is perturbed (real part of the frequency of the black hole quasi-normal modes) [67].

Numerical spinning black hole solutions
In this section we first validate our numerical infrastructure against well-known results, namely the Kerr black hole, and then proceed to use it to obtain spinning black holes in a modified gravity theory, the Einstein-scalar-Gauss-Bonnet theory.

Validating the code against the Kerr black hole
To validate our numerical infrastructure we will solve the axisymmetric vacuum Einstein equations to numerically obtain the Kerr solution, and compare with analytical results. We choose to solve the the following combination of the field equations which diagonalize the Einstein tensor with respect to the operator ∂ 2 r + r −2 ∂ 2 θ : In Fig. 8 we present the results for the comparison of the metric functions obtained numerically with the analytically known ones for a Kerr black hole with χ = 0.6. Given that in this case the initial guess cannot be the Kerr metric itself, to obtain the results in Fig. 8 we used a Schwarzschild black hole with comparable r H * . The maximum observed error is of O (10 −13 ) for the metric function h, with all other metric functions being successfully obtained to machine precision. We also explored the whole domain of existence of Kerr black holes, comparing numerically obtained physically relevant quantities with analytical ones, see Fig. 9 below. These include the mass M , angular momentum J, horizon area A H and Hawking temperature T H of the black holes. Furthermore, we computed the (normalized) Smarr relation in Eq. (68). Overall, in all quantities we have found remarkable agreement between numerical and analytical results, with the Smarr relation providing accurate maximum error estimates. We also observe that errors are higher when the black holes approach the extremal case (χ → 1). This is because in the extremal limit, our setup is not valid and another metricänsatz is needed (see e.g. Ref. [14]).

Einstein-scalar-Gauss-Bonnet Gravity
Einstein-scalar-Gauss-Bonnet (EsGB) theories of gravity are a popular set of scalar tensor theories of gravity that have been extensively studied [16][17][18][19][20][21][22][23][24][25][26][27][28], and which admit black hole solutions different to those of GR. Here we use this set of theories to test our methods and code in a non-trivial, but previously studied setting. EsGB theories are described by the action where φ is a real scalar field that couples non-minimally to the Gauss-Bonnet term via the coupling function ξ(φ), and where α is a coupling constant with dimensions of length squared. No closed-form black hole solutions are known in these models, even in the static case. One is therefore forced to resort to numerical methods to study black holes in these theories. The field equations of the action (98) are where is the double-dual Riemann tensor (the square brackets denote anti-symmetrization). The scalar field equation is where the dot denotes differentiation with respect to the scalar field φ. In the stationary and axisymmetric setting, we find that the scalar field is subject to the boundary conditions [49,50] ∂ while the boundary conditions for the metric functions remain those given above. We therefore choose the same spectral expansion for the scalar field as we did for the metric functions. Black holes in the EsGB theory should obey the Smarr formula (68), which becomes where * and the entropy is given by Eq. (67) that in the EsGB case becomes where R (2) is the Ricci scalar of the induced metric on the horizon. We will focus on two coupling examples, the linear coupling and the exponential coupling We find that for the exponential coupling the Smarr relation takes a rather simple form where Q s is the scalar charge of the solution, appearing in the asymptotic expansion of the scalar field It can also be proved that for the linear coupling the following relation holds [68] Q s = 2παT H .
In what follows we use the relations in Eqs. (107) and (108) to address the accuracy of our numerical solutions for the exponential and linear couplings respectively. This is necessary as closed-form solutions are unknown. We use the same combination of field equations as in the Kerr case (Eq. (97)), along with the scalar field equation (100) to solve the system. To solve the system we use a comparable Kerr black hole as an initial guess for the metric functions, and for the scalar field we use the perturbative solution [49,50] φ ≈ α r 2 We present the accuracy estimate results (in a part of the domain of existence) using the Smarr relation for the exponential coupling and the relation in Eq. (108) for the linear coupling in Fig. 10. We observe that errors, as measured by the relations (107) and * This relation can also be written as provided the coupling does not obey ξ(φ) ∝ ξ (φ) and the scalar field asymptotically vanishes. This is advantageous from a numerical point of view because no second derivatives of the scalar field are required, increasing the accuracy in computing M s .  (108), are small and similar to those presented for the Kerr black hole in Fig. 9, despite a dramatic increase in the complexity and number of terms in the field equations. Our results also agree remarkably well with perturbative solutions, such as the ones obtained in Ref. [50].
As another test to the code, in Figure 11 we plot the accuracy as estimated by the Smarr relation (107) as a function of both resolutions N x and N θ . We observe exponential convergence, similarly to the toy model presented in Fig. 4. Note that the Smarr relation provides only an estimate of maximum error -recall the Kerr case, where most metric functions were actually obtained to a precision of ∼ O (10 −16 ) but the Smarr relation attained errors on the order of ∼ O (10 −13 ).
To further demonstrate the capabilities of our code, in the following we present some results for the physical properties of EsGB black holes. A plot of the ergoregion for a dilaton black hole with γ = 1, χ = 0.1 and α/M 2 = 1.15 can be found in Fig.  12. In Fig. 13 we plot |1 − S| as a function of x and θ, where S is the speciality index defined in Eq. (80), for the same EsGB black hole as before, where we can observe that the spacetime is not algebraically special, being Petrov type I. Spinning EsGB black holes were always observed to be Petrov type I. * The perimetral location and angular frequencies at the ISCO and light rings of EsGB dilaton black holes (γ = 1) are compared with those of a Kerr black hole (with the same χ and M ) in Fig. 14. Note that we have neglected any couplings between the dilaton and matter (see e.g. [48,69]). We have compared our results in the static and slowly rotating cases with those in Ref. [69], observing remarkable agreement (in the appropriate setup). From Fig. 14 we observe differences of a few percent in most cases, with the most drastic differences occurring for the location of the co-rotating light ring due to its proximity to the horizon. The qualitative behaviour is as follows: the perimetral radius of both the ISCO and the light ring decreases with α/M 2 , and the opposite happens for the angular frequencies * . Co-rotating orbits are most affected, and black hole spin enhances the differences of co-rotating orbits with respect to the Kerr case.

Comparison with other codes
Similar codes to the one we have developed in this chapter are scarce. Indeed, most of the numerical studies of spinning black holes in modified theories of gravity make use of the non-publicly-available FIDISOL/CADSOL solver [45][46][47], which implements a finite difference method together with the root finding Newton-Raphson method. The solver is written in Fortran and was first developed in the eighties. Works that use the FIDISOL/CADSOL solver can be found e.g. in Refs. [13-15, 21, 23, 24, 26, 28, 48]. Some of these works have applied the FIDISOL/CADSOL solver in studies of EsGB gravity, much like we did here. However, they report an error of order O (10 −3 ), as estimated by the Smarr relation. In the appendix of Ref. [70], the author gives a comprehensive overview of the FIDISOL/CADSOL solver, benchmarking it against the Kerr solution, with results again showing errors several orders of magnitude higher than those presented in Fig. 9. More recently, in Ref. [50] the authors developed the eXtreme Partial Differential Equations Solver (XPDES) code which is publicly available, to address similar problems. The code is written in C language, and implements a finite difference method to solve the field equations, similarly to the FIDISOL/CADSOL package. It makes use of the software Maple to export the field equations to many large C programming files. Ref. [50] χ=0.1, co-rotating does not discuss errors as estimated by Smarr relations, instead, they (also) benchmark their code against the Kerr solution, and compare their EsGB results to perturbative solutions, finding good agreement. They report typical maximum errors on obtaining the Kerr solution of O (10 −6 ), which represents a good improvement when compared with the FIDISOL/CADSOL package, especially given that the XPDES code is opensource and publicly available. Our code is written in the Julia programming language, which when compared with complied languages such as C code makes it logistically easier to use and adapt, and to implement new models. In our implementation the field equations and boundary conditions are written in a very simple way. For example, the boundary condition f − 2∂ x f = 0, is written as a residual in code language as f − 2 * df dx.
The code is memory efficient and fast, making use of pseudospectral methods as explained above, with solutions to the field equations being obtained in the order of a few seconds in laptop-class computers. In our (limited) comparisons with the XPDES code, we found that where our code took only a few seconds the XPDES code would take minutes to achieve a lower accuracy.
The results of this section, for example in Figs. 9 and 10, show that the accuracy of our code is many orders of magnitude better than the accuracy presented by either the FIDISOL/CADSOL package or even the XPDES code.
Once a solution to the field equations has been obtained, our code has built-in functions to compute all the physical properties of the black holes discussed in section 3.4, therefore allowing for a simple and comprehensive study of different models.

Conclusions
In this paper we have reviewed the spectral method for solving differential equations and subsequently argued that such methods are ideal for finding stationary and axisymmetric black hole solutions in modified theories of gravity. In particular, they allow complicated field equations and boundary conditions to be implemented in a straightforward manner. We showed how this can be done, and have implemented the method in a new code. To show it in action, and to benchmark its performance against other codes, we applied the code in the GR setting, and verified that the solution found is extremely close to the known Kerr black hole. We then applied it to a popular set of modified theories of gravity, Einstein-scalar-Gauss-Bonnet gravity, where it is known that black hole solutions different from Kerr exist. In this latter setting we verified the accuracy using analytical expressions that should hold identically. We found that even in the Gauss-Bonnet setting our code took just seconds to find accurate spinning black hole solutions.
Within the code we have also implemented many built in functions to calculate black hole properties of physical interest. In the future, obtained solutions together with these functions could be used to study a huge range of phenomena observational interest. Other possible studies include the quasi-normal modes of black hole mergers (hence permitting realistic data analysis with Bayesian methods), the electromagnetic emission from accretion disks, black hole shadows, and our code's solutions could also be used as seed solutions for numerical evolutions. Given that the code has been completed only recently, we have, however, not yet applied it widely. Although a first application in research work to EsGB theories is contained in Ref. [71]. In the future we hope to apply the code to other theories, such as the so called regularized 4D-Einstein-Gauss-Bonnet gravity theory [29][30][31][32][33][34][35][36][37] where thus far spinning black holes have not been found, and use it to further understand and constrain such theories using the physical properties described.
supported by a Royal Society University Research Fellowship.