Exact Black Hole Solutions in Modified Gravity Theories: Spherical Symmetry Case
Andrew Sullivan, Nicol\'as Yunes, Thomas P. Sotiriou
TL;DR
This paper introduces a new numerical code for computing stationary black hole solutions in various modified gravity theories, enabling detailed predictions of observable phenomena to compare with astrophysical data.
Contribution
A novel numerical method and code for obtaining stationary, spherically symmetric black hole solutions in a wide class of modified gravity theories.
Findings
Validated the code with known solutions in General Relativity.
Generated new solutions in scalar-Gauss-Bonnet gravity with different couplings.
Compared analytical and numerical results for physical observables like ISCO and photon sphere.
Abstract
Detailed observations of phenomena involving black holes, be it via gravitational waves or more traditional electromagnetic means, can probe the strong field regime of the gravitational interaction. The prediction of features in such observations requires detailed knowledge of the black hole spacetime, both within and outside of General Relativity. We present here a new numerical code that can be used to obtain stationary solutions that describe black hole spacetimes in a wide class of modified theories of gravity. The code makes use of a relaxed Newton-Raphson method to solve the discretized field equations with a Newton's polynomial finite difference scheme. We test and validate this code by considering static and spherically symmetric black holes both in General Relativity, as well as in scalar-Gauss-Bonnet gravity with a linear (linear scalar-Gauss-Bonnet) and an exponential…
| sGB | EdGB | ||
| 0 | 2 | ||
| 1 | 2 | ||
| 3 | 2 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 3 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 4 | ||
| 1 | 4 | ||
| 3 | 4 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 6 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 13 | 12 |
| sGB | EdGB | ||
| 0 | 2 | ||
| 1 | 2 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 3 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 4 | ||
| 1 | 4 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 6 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 11 | 10 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 11 | 12 |
| sGB | EdGB | ||
| 0 | 1 | ||
| 1 | 1 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 2 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 3 | ||
| 1 | 3 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 0 | 5 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 12 | 8 | ||
| ⋮ | ⋮ | ⋮ | ⋮ |
| 7 | 11 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsCosmology and Gravitation Theories · Black Holes and Theoretical Physics · Relativity and Gravitational Theory
Exact Black Hole Solutions in Modified Gravity Theories:
Spherical Symmetry Case
Andrew Sullivan
Department of Physics, Montana State University, Bozeman, MT 59717, USA
Nicolás Yunes
Department of Physics, Montana State University, Bozeman, MT 59717, USA
Thomas P. Sotiriou
School of Mathematical Sciences & School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD, United Kingdom
Abstract
Detailed observations of phenomena involving black holes, be it via gravitational waves or more traditional electromagnetic means, can probe the strong field regime of the gravitational interaction. The prediction of features in such observations requires detailed knowledge of the black hole spacetime, both within and outside of General Relativity. We present here a new numerical code that can be used to obtain stationary solutions that describe black hole spacetimes in a wide class of modified theories of gravity. The code makes use of a relaxed Newton-Raphson method to solve the discretized field equations with a Newton’s polynomial finite difference scheme. We test and validate this code by considering static and spherically symmetric black holes both in General Relativity, as well as in scalar-Gauss-Bonnet gravity with a linear (linear scalar-Gauss-Bonnet) and an exponential (Einstein-dilaton-Gauss-Bonnet) coupling. As a by-product of the latter, we find that analytic solutions obtained in the small coupling approximation are in excellent agreement with our fully non-linear solutions when using a linear coupling. As expected, differences arise when using an exponential coupling. We then use these numerical solutions to construct a fitted analytical model, which we then use to calculate physical observables such as the innermost stable circular orbit and photon sphere and compare them to the numerical results. This code lays the foundation for more detailed calculations of black hole observables that can be compared with data in the future.
I Introduction
The recent discovery of gravitational waves Abbott et al. (2016a, b, 2017a, 2017b, 2017c, 2017d) has inaugurated a new era of multi-messenger astrophysics that opens up an entirely new avenue to test Einstein’s theory of General Relativity (GR) in the extreme gravity regime Yunes and Siemens (2013). GR is thus far in excellent agreement with experiments and observations. However, this statement does rely on the assumption that dark matter or dark energy do not signal a deviation from GR. Moreover, GR has so far resisted quantisation attempts. These considerations have provided motivation for exploring modified theories of gravity.
Unsurprisingly, modified theories typically increase the complexity of the field equations to such a degree that the calculation of observable predictions becomes incredibly difficult and sometimes even seemingly impossible if working analytically. Historically, experimental tests that probe the weak-field regime prompted the calculation of solutions and observables through perturbation theory Will (2014), and more recently, this has been further justified through effective field theory arguments Cooney et al. (2009); Burgess (2004, 2007). A by-product of assuming that GR modifications are small relative to GR predictions is that analytical calculations become tractable again. But once one begins to probe the extreme gravity regime, perturbative techniques need not be well-justified, as they typically eliminate strong-field instabilities that could have observational consequences. The typical example of this is spontaneous scalarization, a process through which large modifications to GR arise in a class of scalar-tensor theories when considering non-linear solutions for neutron stars above a certain compactness Damour and Esposito-Farèse (1998); Sampson et al. (2014) or black holes within a certain mass range Silva et al. (2018); Doneva and Yazadjiev (2018).
The need for more precise solutions and observable predictions without the use of perturbation theory then becomes clear, and we here take steps in this direction. We present a numerical infrastructure to produce exact solutions that describe stationary black holes in a wide class of modified theories of gravity. This paper is focuses on describing the numerical setup and applying it in the simplest case of spherically symmetric solutions that can act as a benchmark. By “exact” we mean solutions obtained numerically and without the use of perturbative techniques, including small coupling expansions. The code is computationally efficient, converging to an answer of prescribed accuracy within only a few iterations, and thus allowing for calculations in laptop-class computers. We validate this infrastructure against known analytic solutions in General Relativity, as well as in a specific member of the class of modified theories to which this infrastructure is applicable.
The basic problem our numerical infrastructure solves is that of finding the solution to an elliptic system of nonlinear and coupled differential equations. The typical approach to solve this problem numerically is to discretize the differential equations through a finite difference scheme. This leads to a residual error on the solution of the system of equations that one wishes to minimize. Our algorithm employs a Newton polynomial method to discretize the equations, which appropriately allows for an easier analytical evaluation of the Jacobian of the system. By minimizing the residual, we can then iteratively converge to the true solution using a root-finding algorithm, such as the Newton-Raphson method. Our code uses an additional relaxation factor to improve convergence, as well as compactified coordinates to properly set boundary conditions, and adaptive mesh refinement near the boundaries. With all of this machinery, our code typically converges to the user prescribed tolerance in 1–3 iterations.
After constructing this numerical infrastructure, we implement it in a few different scenarios. We begin by validating the algorithm through a simple toy problem whose analytic solution is known, and then through the calculation of the Schwarzschild solution in General Relativity. Since in both cases an analytic solution is known, we can easily compare it to our numerical solutions point-wise in the entire domain. We find that our code converges to the correct (analytically known) solution in one and three iterations respectively to within the tolerance we specified.
With the validation complete, we then construct stationary and spherically symmetric black hole solutions in scalar-Gauss-Bonnet (sGB) gravity, a well motivated modified theory that is a member of the quadratic gravity class Yunes and Siemens (2013); Berti et al. (2015); Barack et al. (2018). Ours is of course not the first study of compact objects in sGB gravity. Mathematically, the theory evades the no-hair theorems of General Relativity, allowing black holes to have non-trivial scalar hair Campbell et al. (1992); Mignemi and Stewart (1993a); Kanti et al. (1996); Yunes and Stein (2011); Sotiriou and Zhou (2014a, b), while preventing neutron stars from having a monopole scalar charge Yagi et al. (2016), thus making it extremely difficult to place constraints with binary pulsar observations.
Compact objects in scalar-Gauss-Bonnet gravity are typically studied under two forms of the action which depend on the coupling function between the massless scalar field and the Gauss-Bonnet invariant. When the scalar field is coupled through an exponential to the Gauss-Bonnet invariant this is commonly referred to as Einstein-dilaton-Gauss-Bonnet (EdGB) gravity. In the regime where the scalar field is small, one can approximate the exponential as a linear coupling to the Gauss-Bonnet invariant which is commonly referred to as the linear scalar-Gauss-Bonnet gravity. This is the terminology that will be used throughout this paper namely: ‘linear sGB’ to refer to the linear coupling function and ‘EdGB’ to refer to the exponential coupling function.
Stationary black holes have been found in linear sGB analytically using the small coupling limit approximation, both in spherical symmetry Yunes and Stein (2011); Sotiriou and Zhou (2014a, b) and in axisymmetry using a slow-rotation approximation Campbell et al. (1992); Mignemi and Stewart (1993b); Ayzenberg and Yunes (2014); Pani et al. (2011); Maselli et al. (2015). Reference Sotiriou and Zhou (2014b) also studied numerical non-perturbative spherically symmetric black holes in linear sGB. In EdGB, numerical solutions are known in spherical symmetry Kanti et al. (1996); Torii et al. (1997); Alexeyev and Pomazanov (1997); Guo et al. (2008); Pani and Cardoso (2009) and in axisymmetry Kleihaus et al. (2016), but typically these numerical solutions are obtained from a proprietary code that has not been tested for black hole problems of this kind. It is also worth mentioning that, dynamical evolution of black holes and binaries in sGB gravity has also received attention recently Benkel et al. (2016, 2017); Witek et al. (2018); Ripley and Pretorius (2019).
Our numerical infrastructure is general enough to allow the tackling of any coupling function between the massless dilaton and the Gauss-Bonnet term in the action, and we here explore both the linear sGB and EdGB coupling functions. In the linear case, we find that the exact numerical solution agrees spectacularly well with the perturbative analytic solution that assumes weak-coupling, while differences arise in the exponential coupling case for large enough coupling. An example of this is shown in Fig. 1, where we present the relative fractional correction to the ADM mass (left) and the scalar monopole charge (right) as a function of the dimensionless sGB coupling parameter . Although the differences between the EdGB and linear sGB solutions appears large for large , the deviation is actually only appreciable near the horizon, as we will show later.
With these exact numerical solutions at hand, we construct analytical fitted models to allow for the rapid computation of physical observables, such as the location of the innermost circular stable orbit (ISCO) and the light ring. The fitting coefficients are available online by request though we provide a few examples in Appendix A. Figure 2 shows the fractional change in the location of the ISCO (left) and the light ring (right). We again find comparatively very good agreement between the analytic perturbative result and the exact numerical result of linear sGB for all s considered, and some disagreement with the exact numerical result of EdGB for large . Observe also that the observables computed with the exact numerical solution agree extremely well with those computed with our analytical fitted models.
Our paper differs from previous work in that its goal is not to study a particular theory of gravity, but rather to develop a computational infrastructure that is (i) free and open to the public, (ii) computationally accurate and efficient, (iii) taylor-made for finding black hole solutions in modified theories of gravity with vector or scalar non-minimal couplings, and (iv) able to provide enough accurate data to compute analytical fitted models. The algorithm presented here applies to spherically symmetric scenarios, but extensions to axisymmetry are straightforward and under way. The algorithm developed here can thus be used to calculate certain astrophysical observables, such as the electromagnetic emission of accretion disks around black holes Abramowicz and Fragile (2013), the shadows of black holes Ayzenberg and Yunes (2018), or the quasinormal modes of black hole mergers Blázquez-Salcedo et al. (2016) on the fly and with high precision, allowing for realistic data analysis investigations with Bayesian methods.
The remainder of this paper is organized as follows. Section II outlines the numerical algorithm as applied to a general system of equations. Section III validates the algorithm through a simple toy problem and a Schwarzschild black hole. Section IV applies the algorithm to sGB gravity and derives the results described above. Section V constructs a fitted analytical model from the numerical solutions and compares physical observables determined by the numerical solutions and the fits. Finally, Section VI summarizes our results and points to future directions. For the remainder of this paper we use the following conventions: Greek letters denote spacetime indices; the metric has the spacetime signature ; we use geometric units where .
II Numerical Methods
Our numerical algorithm extends the method presented in Schönauer and Weiß (1989) to build a partial differential equation solver that finds exact black hole solutions in an arbitrary modified theory of gravity. The algorithm is split into three main parts: the relaxed Newton-Raphson method, the discretization method, and the discretization error estimation. We begin by outlining the Newton-Raphson method in the continuum limit in Sec. II.1. The discretization method takes our generic system of nonlinear differential equations and recasts them for evaluation on a discrete domain through a Newton polynomial scheme outlined in Sec. II.2. This method naturally introduces errors that must be estimated and controlled, as outlined in Sec. II.4. Combining these results, we describe the application of the discrete Newton-Raphson method with a relaxation factor in Sec. II.5, which reduces the computational problem to solving a system of coupled linear equations through two iterative Krylov subspace methods, detailed in Sec. II.6.
Given the description above, clearly this section contains a considerable amount of descriptive detail about the numerical methods that constitute the computational infrastructure we have deployed to find numerical black hole solutions in modified gravity. Our hope is that these details will serve as a sort of tutorial for theoretical physics experts that are less familiar with the numerical methods we use here. Readers already familiar with these numerical methods may wish to skip this section altogether and jump directly to the validation of the algorithm and its application to black holes in sGB in Sec. III and beyond.
II.1 Iterative Relaxation Method in the Continuum Limit
We start with a system of , second-order, nonlinear elliptic partial differential equations. In the simplest 1-D (spherically symmetric and stationary) case, these equations are ordinary differential equations of some independent variable , such that in operator form
[TABLE]
where is an matrix of nonlinear differential operators and is the solution vector111For the remainder of this paper, the word vector stands for a standard Euclidean vector in flat space. with elements. The elements of are the fields in the problem, which we will sometimes denote with capital latin indices . We see then that the special solution vector is annihilated by the differential operator .
In general, however, one does not know the solution to the differential system, so a generic vector will not be annihilated by . Rather, for a generic vector that is not a solution to the differential system one has
[TABLE]
where the residual vector has elements and is a function of . The generic vector is a functional that spans the vector function space that contains the solution vector and any other vector that does not satisfy the differential system 222In particle physics, the set of functions that satisfy a differential system of equations is sometimes referred to as “on shell,” while those that do not are referred to as “off shell.” Thus spans the space that contains both on and off shell functions.. Clearly then, once one finds the correct vector , namely , then .
The vector is also subject to boundary conditions. We can express these in the form
[TABLE]
where is another matrix typically composed of a combination of real numbers and first-order differential operators, while is the boundary of the spatial domain. The boundary vector also has elements and is a set of constants defined on the boundary . Defining the boundary conditions in this way allows the imposition of Neumann boundary conditions for a subset of the components of and Dirichlet boundary conditions for other components.
Let us assume that we have a vector that we know is close the the actual solution vector . If so, the difference between it and the actual solution is some other small vector . Let us refer to the latter as the correction vector, which is mathematically defined via
[TABLE]
We begin by linearizing around through a first-order Taylor expansion, analogous to the familiar Taylor expansion of about . Doing so, we find
[TABLE]
We then evaluate this expression at the solution vector and obtain
[TABLE]
Using that and substituting the correction vector defined in Eq. (4), this can be simplified to
[TABLE]
where we have defined the Jacobian matrix
[TABLE]
We now wish to solve for the correction vector , which implies we must either invert the Jacobian or solve the linear system in Eq. (7). In practice, matrix inversion is typically more computationally expensive than linear system solving, so the latter is the method we use here, which we will explain in more detail in Sec. II.6. If one knew the correction vector exactly, one could then find the solution vector from Eq. (4). But in reality the linearizaton in Eq. (5) implies the correction vector we find by solving Eq. (7) is only an approximation to the true correction vector. This implies that to find the true correction vector we must apply this procedure iteratively.
Let us then describe the first couple of iterations of this procedure. We start with an initial guess that we know is close to the true solution, where the superscript in parenthesis is the iteration number. In this paper, the initial guess can be chosen to be either the GR solution, or an approximate solution for all fields in the modified theory. With this initial guess, we then find the initial residual vector
[TABLE]
Since the initial guess is not a solution to the differential system, the residual vector does not vanish, and we must thus correct the initial guess to find the first-iterated solution
[TABLE]
This requires the calculation of the zeroth-iterated correction vector , which we find by solving the linear system
[TABLE]
where the Jacobian is evaluated on the initial guess . This procedure then yields , and now it can be repeated until the nth-iteration to the solution is sufficiently close to , i.e. until the residual vector is below some specified tolerance .
II.2 Discrete Representation through Newton’s Polynomials
In order to numerically solve the differential system described above, one first needs to discretize it on a finite numerical grid. We here use Newton’s (centrally divided difference) interpolation polynomial method, which we describe next.
Newton’s interpolation polynomial provides a continuous local representation of a discrete function given by a set of data points. This procedure is a discrete analog of using a Taylor series to represent an approximation to a continuous function as a local polynomial about some point , namely
[TABLE]
where the primes denote derivatives with respect to the independent variable .
Now imagine that instead of a continuous function , we have a discrete function known only on a discrete collection of data points . Notice that the lower-case Latin subscript here does not denote the components of the of the previous subsection, but rather the element at which we evaluate the discrete function . For notational convenience, we identify with the discrete function evaluated at each point , namely
[TABLE]
How do we now approximate the discrete function in a neighborhood of some point in the spatial domain? We are tempted to use a Taylor expansion again, but because our function is discrete, we cannot take analytical derivatives as we did before, and instead we must use a finite difference approximation. Let us then define a neighborhood around some point as the region in the discrete spatial domain around which we wish to approximate our function. Clearly then, the value of the function at the point is simply . Let us further temporarily assume that the data points are equidistant, such that .
We are now almost ready to define our finite difference approximation, but first we must choose which discrete points in the neighborhood of we will use to approximate the function. For example, if we choose to use the points , then the Newton interpolation polynomial of our discrete function is
[TABLE]
The coefficients of the second and third terms in Eq. (14) are simply the first- and second-order forward finite difference approximation to the first and second derivative of the function in the limit, namely
[TABLE]
and
[TABLE]
respectively. Moreover, using that
[TABLE]
we notice that the term is the same polynomial that appears in the Taylor expansion, and that the additional term cancels with one of the in the denominator of Eq. (14). We can then combine this with the first-order term to obtain
[TABLE]
Let us now compare the above result to the first-order finite difference derivative in Eq. (15). First, the derivative now includes three points instead of two. Second, the accuracy has been increased to . This suggests that as we add increasingly higher-order polynomials to our discrete representation, they will automatically include corrections to each of the previous lower-order polynomials. Thus, the Newton interpolation polynomial is simply a discrete analog to the Taylor series, and they are formally equivalent in the continuous limit.
The full construction of the Newton interpolation polynomial requires certain Newton basis polynomials and certain divided differences coefficients, which are a generalization of the finite difference coefficients of the above simple example. The choice of points to include leads to the distinction between a forward and a backward divided difference, and this influences the directional finite difference that reduces to the Taylor derivative in the continuum limit. As a simple example, consider again the approximation of a discrete function in a neighbourhood around a point but this time using the points immediately surrounding : . Equation (14) then becomes
[TABLE]
which results in the first and second-order central finite difference equations in the continuum limit. Similarly if we use the points before only, namely , we then obtain the backwards finite difference equations in the continuum limit. The accuracy of the central finite difference equations is improved by a factor of over the forward or backwards finite differences.
Let us generalize the systematic application of these choices through the introduction of a stencil, namely a Euclidean vector whose elements are the labels of the points included in the evaluation of our Newton polynomial. For example, for a forward, central, and backwards Newton polynomial representation of a discrete function about a point , the stencil and respectively.
With this at hand we can write the general Newton interpolation polynomial of a discrete function about any point in the spatial domain as
[TABLE]
where recall that is given by Eq. (13) using stencil notation, and where
[TABLE]
and are the Newton basis polynomials, defined by
[TABLE]
If the grid is uniform and , one can check that this equation reduces to Eqs. (14) or (19) with the respective stencil.
The index in Eq. (20) is the order of the Newton polynomial and it indicates how many grid points are used and the maximum finite difference derivative order. In practice, as the order increases, we successively add points to either side of to keep the coefficients as central as possible, so our stencil has the form and will have elements. For even the stencil will be purely central, whereas for odd the stencil will be slightly forward. On the boundary of our domain, we must add points in a one-sided manner (either forward or backwards), and we must add an extra point beyond the number of points we keep away from the boundaries, to keep the accuracy of the represetation comparable to the central differences. Note that the specific sequence of the elements of does not change the resulting polynomial.
As an example of an application of this, let us calculate and compare the Newton interpolation polynomial representation and the comparable Taylor series expansion of a toy function
[TABLE]
on a uniform discretized grid where around the point . With a maximum order , the centralized stencil is . Figure 3 shows the Taylor series expansion and the Newton polynomial representation of our toy function order by order. The first three terms are shown in Eqs. (12) and (19), from which we see the order approximation in both cases is simply a constant . The order is a line whose slope is calculated by either the derivative of the function evaluated at , or the first-order finite difference evaluated at . Observe that the agreement between the Newton interpolation polynomial representation of the discrete function and the Taylor series expansion of the continuous function is very close but the agreement is better for even-orders of than for odd-orders. This is a result of the even-orders of being purely central for example for while the odd-orders are slightly forward for example which diminishes the accuracy. It is for this reason that we will restrict our choice in to be even.
II.3 Derivatives of Newton Polynomials and Discretization Error
One of the main advantages of using a Newton interpolation polynomial representation of a discrete function is that with it one can take analytic derivatives of the discrete function. For example, the derivatives can obtained by taking an analytic derivatives of the basis functions , which are the only dependent terms in Eq. (20). The two relevant derivatives for second order equations are
[TABLE]
and
[TABLE]
as one can verify by direct differentiation.
This property of Newton interpolation polynomial representations allows us to replace the action of any differential operator on our fields at any grid point with a finite difference coefficient coupled to the neighboring grid points, but this introduces some error. To quantify this error, let us introduce boldface vector notation to denote a discretized vector, i.e. the collection of data points of the Newton polynomial representation of a field evaluated on every point of the discretized grid. For example, for a given field . By construction, of course, the discretized field itself has no discretization error on any point of the discretized domain, but this is not the case for derivatives of the field. Let us then introduce the discretization error vector and the discretization operator , such that
[TABLE]
because by construction . These discretized errors will be incorporated into the Newton’s method of Sec. II.4 to maintain numerical accuracy.
Before moving on to the next subsection, let us here address some potential confusion due to the different choices of notation for a vector that we have employed. The boldface vector notation of the previous paragraph denotes a discretized vector, so its components are the values of the Newton polynomial representation of a function on each grid point of our spatial domain. The arrow vector notation of the previous section denotes a generic field vector, so its components are the different fields in the differential system. When we return to our general system of differential equations, we must generate a Newton polynomial for each component of , and then, we must discretize each of these representations on grid points. The quantity then becomes a 2-D matrix (see Eq. (27) below) with each element denoted as .
As we will see later, we will have to fold every discretized vector in our system of equations so that we can take proper element by element partial derivatives. To do this, we introduce a new index and fold the matrix into a single vector where , namely
[TABLE]
This folding replaces our system of continuous differential equations with a linear system of equations.
II.4 Minimization of the Discretization Error
Let us consider how the discretization error in the derivatives of the Newton polynomial representation of our discrete fields propagates into our system of differential equations. For a generic vector , our discretized differential system becomes
[TABLE]
where is the linearized differential operator resulting from the expansion of . As before, let us choose our vector to be close to the solution vector up to a small correction vector . Discretizing these relations we then have that
[TABLE]
where we have included a new correction vector to attempt to also minimize the discretization error. As in Sec. II.1, we then Taylor expand Eq. (28) to find
[TABLE]
where the discretized Jacobian is evaluated at , and we have defined the discretization error vector . The above equation can be satisfied by separately and simultaneously requiring that
[TABLE]
for the discretized correction vector and
[TABLE]
We recognize Eq. (31) as the discretized version of Eq. (7). By balancing both of these equations we can monitor the numerical accuracy of our discretization error.
The equation for the discretization error vector, Eq. (32), describes the change to our solution due to the discretized error vector . To control the discretized error, we wish to require that the relative correction is below a specified tolerance. To solve Eq. (32), we first must calculate which was defined as the discretized differential operator acting on :
[TABLE]
In practice, we will generally not have our system of nonlinear equations in operator form. For example, when calculating the Einstein field equations we will have the set of equations in the form of the residual vector as a functional of the fields themselves. If we wish to extract the linearized differential operator of Eq. (33) we must calculate it from the residual vector.
As a simple example, consider the following single continuous nonlinear differential equation in operator notation:
[TABLE]
If we wish to calculate we obtain,
[TABLE]
From this, we identify the linearized operator as
[TABLE]
and the general discretization error vector as
[TABLE]
where as before due to the Newton polynomial construction.
The derivative error can then be computed from Eq. (26), namely via
[TABLE]
but in general we do not know the analytical derivatives of any of the fields . Following Schönauer and Weiß (1989), we estimate the discrete derivative error from the difference between the discretized derivative at order and that at order . This is valid as long as the error decreases for increasing order, and if so, it provides a reasonable estimate as demonstrated in Schönauer et al. (1981). With this approximation, our discretized derivative error is
[TABLE]
Let us then describe the procedure we will follow to solve Eq. (30). First, before finding the field correction to the zeroth-iteration, we ensure the discretization error is under a specified tolerance by using Eq. (32) to determine the discretized grid size. This is achieved through a separate Newton-Raphson subroutine that we describe below. Once this is achieved, we then solve Eq. (31) to find the zeroth-iteration correction vector, which allows us to update the solution. This full procedure is then iterated, at each step ensuring that the discretization error is under control, until the residual in Eq. (31) is below a prescribed tolerance. We will detail this process in Sec. II.5.
Let us now describe the Newton-Raphson subroutine that we use to minimize the discretization error. First, we calculate from Eq. (37) using Eq. (39) for the discretization derivative error. We then find the correction to the solution by solving
[TABLE]
using one of the iterative methods in Sec. II.6. We stop the subroutine when the following condition is satisfied
[TABLE]
where is the maximum norm (also called infinity norm or supremum norm) among all fields and grid points. This condition forces the relative correction to the solution from the discretization error to be below some specified tolerance . In turn, this ensures that our final solution is within the discretized error tolerance.
If the above condition is not met for the chosen grid discretization, we adaptively adjust the grid step size until the relative discretization error correction is below the specified tolerance. To adjust the step size at each grid point defined as , we rescale the latter by the ratio of Eq. (41) with an additional “safety” factor of , and relax the ratio by the Newton polynomial order , such that
[TABLE]
where is the maximum norm among all fields evaluated at grid index .
Once the new step sizes are calculated, the entire grid must be adjusted and the condition in Eq. (41) must be rechecked. In the grid adjustment, we must ensure that the location of a new grid point remains within half an old step size of the old grid point. Then, each function is interpolated using the Newton polynomial representation calculated from the old grid point locations. From the new interpolated solution, and are recalculated and the condition of Eq. (41) is rechecked. If the condition is not satisfied, the process using Eq. (42) is repeated until the step sizes on the grid become sufficiently small that the condition is satisfied. Once the discretization error is sufficiently under control we are ready to finally apply the Newton-Raphson iterations to the solution vector.
II.5 Minimization of the Residual Vector
As long as we can ensure our discretization error is under control, we can safely ignore their contributions to Eq. (30) and we are left with the equation
[TABLE]
which we can now solve following the same steps outlined in Sec. II.1.
The full algorithm is then as follows. Using the system at iteration , , we first find the discretized residual vector
[TABLE]
and calculate the discretization error . We then check Eq. (41) to determine if we need to adjust the step size according to Eq. (42). Once the step size is sufficient, we solve Eq. (43) for the correction vector and the generic vector is updated with a relaxation factor (initially ) added
[TABLE]
Before the new solution is accepted, the residual of the system is re-calculated to ensure convergence, i.e. that the new residual is smaller than the old residual:
[TABLE]
If the convergence condition of Eq. (46) does not hold then the relaxation factor is reduced by half,
[TABLE]
and the new solution and residual is recalculated until Eq. (46) holds or until at which point the algorithm terminates. If the convergence condition does hold then is accepted and the relaxation factor grows by,
[TABLE]
Iterations are stopped when the maximum norm of the residual satisfies the condition,
[TABLE]
We analytically compute both the residual and the Jacobian from our set of field equations in the symbolic manipulation software Maple 2018, which are then automatically exported into the C programming language for evaluation by our algorithm. These are then used by an iterative linear solver described in the next section.
II.6 Linear System Solvers
The solution of linear systems of the form such as Eq. (43) is typically done by either direct or iterative methods. Direct methods find the solution through a finite number of steps. For example, LU decomposition uses gaussian elimination to decompose the matrix into a product of lower and upper triangular matrices which simplifies the computation of the solution. Iterative methods, on the other hand, gradually approach the solution by recursively minimizing the estimated error between the current solution and the exact solution. For very large matrices, direct methods become unwieldy and can become computationally cumbersome so we use two iterative methods here.
The most common type of iterative method used today are Krylov subspace methods. These methods generate a sequence of approximate solutions from the Krylov subspace of and the residual of the linear system defined as such that the corresponding residuals converge to the zero vector. The two popular forms of this type of methods used here are the generalized minimal residual (GMRES) Saad and Schultz (1986) and the biconjugate gradient stabilized (BiCGSTAB) van der Vorst (1992) methods. The error tolerance of both iterative methods was chosen as which effectively places a lower bound on our numerical accuracy.
The GMRES method uses the Arnoldi method to create a Krylov subspace from the Gram-Schmidt process. The method then computes the upper triangular matrix representation of in the Krylov basis. The solution is obtained from minimizing the norm of the residual of the system in this basis. BiCGSTAB is a modified variant of the conjugate gradient method which uses a method of gradient descent to find the minimum of the system of linear equations. BiCGSTAB is the generalized version of this method that applies to non-self-adjoint matrices and uses subroutine applications of GMRES to stabalize the conjugate gradient method.
III Validation
In this section we describe the steps we have taken to validate our algorithm. We first consider a simple ordinary differential equation. This toy problem is a useful example to demonstrate the structure of the Jacobian and to describe the detailed steps involved in applying the Newton-Raphson method. We then move on to General Relativity and consider static, spherically symmetric vacuum solutions to the Einstein equations. This demonstrates how the algorithm handles nonlinear coupled ordinary differential equations that represent black holes with coordinate singularities at the location of the horizon.
III.1 Toy Problem
Consider a simple second order ordinary differential equation in the operator notation of Eq. (1)
[TABLE]
The general solution to this differential equation can be found analytically to be
[TABLE]
If we choose the boundary conditions, and , the solution becomes
[TABLE]
Let us now apply our computational infrastructure to this simple problem. The system is a single ordinary differential equation (M=1) that we wish to solve on a uniform discretized grid of points. Let us choose the following initial guess for our generic vector :
[TABLE]
where the superscript stands for the iteration number, i.e. the initial guess is the iteration. The residual vector on the right-hand side of Eq.(2) can be analytically evaluated to
[TABLE]
which clearly does not vanish anywhere in the domain, i.e. in .
Since the residual does not vanish, we must correct the initial guess by some amount in the next (first) iteration. In order to find this correction, however, we must first find a discrete representation of in terms of a Newton polynomial on a uniform grid at points with . For this toy problem, we pick a Newton polynomial of order , and use a centralized stencil , where is the counter of the element that is closest to the value of .
With this at hand, Eq. (20) yields
[TABLE]
The basis is defined in Eq. (22), and since the order of the polynomial here is , the only basis that contribute are
[TABLE]
The coefficients , while the weighting factors that contribute at this polynomial order are
[TABLE]
Putting all of this together, we then have the discrete representation of shown below in Eq. (III.1).
The Newton polynomial representation of order is the second-order, discrete Taylor expansion of the function about the point . In practice, the only values that can take are in the set . Since this is a closed-form representation of the as a function of , we can evaluate the function and all derivatives at a point straightforwardly
[TABLE]
where , as must be the case since Eq. (III.1) is a discrete Taylor expansion.
[TABLE]
With this representation, we can now evaluate the discretized residual , which is simply given by
[TABLE]
Using the above expressions for the analytic derivatives of the discretized function , we then have that for the discretized residual at iteration
[TABLE]
and at the boundaries
[TABLE]
With the discretized residual, we can now evaluate the Jacobian. From Eq. (43), we have that
[TABLE]
in our toy problem, whose only non-vanishing components are
[TABLE]
which is a tridiagonal matrix. Given this, the discretization error calculated from Eq. (37) is
[TABLE]
with a slightly modified formula near the boundary.
With all of these quantities calculated, we now apply the Newton-Raphson method to invert Eq. (43). Let us then return to vector notation and re-introduce
[TABLE]
such that Eq. (43) becomes
[TABLE]
The inversion of this equation then yields the correction to our th solution, namely . Applying this algorithm, we find that our computational infrastructure converges to the tolerance required in a single iteration. In fact, after a single iteration, even after only requiring the tolerance of Eq. (49). Although the algorithm efficiently minimizes the residual of the system, the residual does not directly correlate with the error between the numerical solution and the exact solution. This is due to the additional errors that are introduced during the discretization procedure described in Sec. II.4. Therefore, the “true” error between the numerical solution and the exact solution is determined by a combination of the discretized residual and the discretization error.
In this toy problem, although the residual was minimized to , the relative discretization error is still from Eq. (41), so the “true” error after 1 iteration is still limited to (see the bottom left panel of Fig. 4). This implies that even though the residual converges far below our desired tolerance, the “true” error is only just below the desired tolerance because the “true” error is also indirectly influenced by the discretization error. The solution, error, and the residual are shown in Fig. 4.
III.2 Schwarzschild Black Hole
In the previous section we solved a simple ordinary differential equation using the method described in Sec. II. We would now like to apply it to solving the elliptic differential equations that arise from the vacuum Einstein equations in spherical symmetry and stationarity. The solution is the well-known Schwarzschild metric to which we can compare our numerical results. We therefore use this example as a benchmark of the performance of our algorithm.
The Einstein-Hilbert action in General Relativity in a vacuum is given by
[TABLE]
where is the Ricci scalar and is the determinant of the metric . Varying the action with respect to the metric gives the vacuum Einstein field equations
[TABLE]
where is the Einstein tensor.
Let us now consider a spherically symmetric and stationary metric ansatz in isotropic coordinates
[TABLE]
where is the isotropic radial coordinate, which is related to the Schwarzschild radial coordinate by
[TABLE]
where is the horizon radius in Schwarzschild coordinates and is the mass.
With this ansatz, the coupled system of ( in our computational infrastructure notation) differential equations that we wish to solve can be found from the and components of the Einstein tensor. These components are
[TABLE]
with primes standing for radial derivatives.
The Schwarzschild solution to these equations in these coordinates is
[TABLE]
The event horizon in isotropic coordinates is located at , as found from the condition .
The boundary conditions at the event horizon and at spatial infinity are determined by regularity and smoothness. At the event horizon we must have
[TABLE]
which follows from evaluation of the analytic solution at . At spatial infinity, asymptotic flatness requires that
[TABLE]
Our computational infrastructure allows us to not only find the numerical solution to Eq. (68), but also to find some observable global quantities that characterize the black hole spacetime. Asymptotically near spatial infinity, the leading order terms of the fields decay as
[TABLE]
where is the Arnowit-Deser-Misner (ADM) mass. Therefore, the ADM mass can be found from the coefficient of the expansion of the numerical solution near spatial infinity. Because of the high Newton polynomial order we use, we can interpolate our numerical solutions very close to spatial infinity to a high degree of accuracy, which allows us to extract the ADM mass from our numerical solutions easily.
The extraction of the ADM mass and the imposition of boundary conditions becomes more precise through the use of a compactified coordinate. We therefore introduce the coordinate , defined by
[TABLE]
and perform a coordinate transformation prior to solving our differential system. This changes our domain of integration from to the finite domain . In these compactified isotropic coordinates, the Schwarzschild solution has the form
[TABLE]
With the differential system in compactified coordinates, we then solve the problem numerically as specified in Sec. II. We begin by replacing each function and differential operator with a discretized Newton polynomial representation of order on a grid of points. We then initialize the numerical solver with a initial guess that is a small perturbation away from the Schwarzschild metric and that vanishes at the boundaries333For this toy problem, we know the exact analytic solution a priori, so we could initialize our solver with it. Doing so, however, would prevent us from validating our computational infrastructure.
[TABLE]
where . One can adjust to improve or worsen the initial guess, which in turn affects the number of iterations required to converge to a solution within the tolerance required.
Applying this algorithm, we find that our computational infrastructure converges to the desired tolerance in 3 iterations. The number of iterations is related to the initial guess, which in this case is controlled by the value of . In the limit as the initial guess becomes the exact solution, and the initial residual decreases below tolerance to within numerical precision.
Unlike in the toy problem from the previous section, the “true” error between the numerical solution and the exact solution is now of comparable order to the residual . This is due to the closed polynomial form of the Schwarzschild solution in compactified isotropic coordinates, which are very well approximated by our Newton polynomial. The comparison between the toy problem and the Schwarzschild application suggests that in problems where we do not have the exact solution to compute the “true” error we cannot use the residual as a direct measure of the error. In other words, even if we have a minimized residual far below the desired tolerance, the solution must be assumed to only be accurate to the desired tolerance. The numerical solution and residuals are shown in Fig. 5.
IV Spherically Symmetric Black Holes in Scalar-Gauss-Bonnet Gravity
In this section we solve the modified Einstein field equations in sGB gravity with both a linear coupling and an exponential coupling function, assuming a vacuum spacetime that is also spherical symmetric and stationary.
IV.1 Action and Field equations
The action in scalar-Gauss-Bonnet gravity in a vacuum is given by
[TABLE]
where is the Ricci scalar and is the determinant of the metric . The real dimensionless scalar field is coupled through a coupling constant which has dimensions of length squared and a function of the scalar field to the Gauss-Bonnet invariant
[TABLE]
We keep the coupling constant around in this section, but in all computation we set , as it can be eliminated through a redefinition of the scalar field and the coupling constant . Note this form of the action differs from that introduced in Yunes and Siemens (2013) by a factor of such that and .
By varying the action with respect to the metric and the scalar field we obtain two field equations. Variation with respect to the metric field yields
[TABLE]
where the scalar field stress-energy tensor is
[TABLE]
and
[TABLE]
Variation with respect to the scalar field yields
[TABLE]
The scalar field is subject to the following boundary conditions: it must be asymptotically flat, and its first derivative must vanish on the horizon in isotropic coordinates,444This follows from the regularity condition on the horizon Kanti et al. (1996); Sotiriou and Zhou (2014a, b). namely
[TABLE]
Utilizing the spherically symmetric metric ansatz in isotropic coordinates from before, the Einstein tensor is still given by Eqs. (71) and (72), and the two new sets of terms in the field equations are:
[TABLE]
where the primes denote derivatives with respect to the coordinate, e.g. and is given below in Eq. (89).
Let us end this subsection with a discussion of the coupling function . In full EdGB gravity, the coupling function is , but in the regime where is small, we can Taylor expand the coupling function as . The -independent term in this expansion is irrelevant as it leads to a theory that is identical to GR due to the Gauss-Bonnet invariant being a topological invariant. In this paper, we focus on numerical calculations for sGB gravity with both an exponential coupling and a linear coupling, namely for theories defined by
[TABLE]
We consider each of these cases separately in the next subsections.
[TABLE]
IV.2 Linear Scalar-Gauss-Bonnet Gravity
In the linear coupling theory, we can find an analytical perturbative solution to the field equations assuming the small-coupling limit, i.e. because is of the order of the curvature length of the system under consideration and we introduce the dimensionless coupling constant . Let us then use a deformed-Schwarzschild ansatz for the metric tensor
[TABLE]
where is a book-keeping parameter and is , and the following ansatz for the scalar field
[TABLE]
Both of these ansatz are assumed valid up to .
Inserting the ansatz in the field equations, we can analytically solve for the metric and the scalar field order by order in , imposing regularity on the horizon and asymptotic flatness at spatial infinity to fix any integration constants. At , and are just the Schwarzschild metric in isotropic coordinates of Eq. (73), while the scalar field vanishes due to asymptotic flatness. At , we find the metric perturbations vanish, while the scalar field perturbation is
[TABLE]
At , we find the first nontrivial correction to the metric tensor and we find that the scalar field perturbation at this order vanishes. Both are given below in isotropic coordinates in Eq. (93). We can then further express these analytic solutions in compactified coordinates for later use which are shown in Eq. (94). These results agree exactly with those found in Yunes and Stein (2011) and Sotiriou and Zhou (2014b) after performing the coordinate transformation from Schwarzschild to isotropic and compactified isotropic coordinates.
[TABLE]
[TABLE]
From these solutions we can obtain the ADM mass and scalar charge from an expansion as , namely
[TABLE]
where and is the bare mass of the black hole that appears in the background (Schwarzschild) metric. Just as the ADM mass is related to the coefficient of the term in an expansion of the metric about spatial infinity, the charge is related to the same coefficient but in the expansion of the scalar field about spatial infinity.
With this analysis in hand, let us now focus on numerically solving the field equations, Eqs. (82) and (85), simultaneously for both metric functions and and the scalar field without any approximations. In order to do so, we employ the computational infrastructure described in detail in Sec. II. In particular, we choose an initial grid of points and a Newton polynomial order . For the actual computation, we set which sets the bare mass of the black hole to . Different black hole masses can be obtained by scaling the radial coordinate appropriately. Note this renders our equations dimensionless and the correct units can be restored through a similar rescaling. The desired tolerance of the solution is which is both placed on the residual and on the relative tolerance of the discretization correction in Eq. (41). The tolerance in the iterative linear solvers described in Sec. II.6 is and places a lower bound on our numerical accuracy. We will also employ the compactified coordinate system defined in Eq. (77) to change our domain of integration to the finite domain . In these compactified isotropic coordinates, the perturbative corrections of Eq.(93) are shown below in Eq. (94).
We can now compare the analytic perturbative solutions to the full non-linear solutions. The top and middle panels of Fig. 6 show the sGB corrections to the metric components as a function of the compactified coordinate for different choices of the coupling constant . Observe that the numerical solution is almost indistinguishable from the analytic perturbative solution everywhere in the domain. The agreement is so remarkable that it is worthwhile exploring the difference between the analytic and the numerical solution, which we do in the bottom panel of Fig. 6. Observe that the difference is indeed very small, ranging from to almost depending on the metric component one studies, and it increases as the coupling strength is increased as expected. We have verified that the residual is always orders of magnitude smaller than this difference in our numerical solutions.
We observe similar behavior in the scalar field. The top panel of Fig. 7 shows the scalar field solved numerically and the analytic perturbative solution both as a function of the compactified coordinate. Observe again that the curves are right on top of each other. The difference between these curves is shown in the bottom panel of this figure, where we see clearly that the difference ranges from to for the largest couplings we considered. As before, we have verified that the residual of the scalar field equation of motion is smaller than this difference for all cases considered. Observe, however, that this time the difference between the numerical solution and the analytic perturbative solution is much larger than it was for either metric component. This makes sense because the modified field equations depend on the scalar field through its stress-energy tensor, which is quadratic in the scalar field. Thus, we expect a difference of for some to contribute a difference of to the metric components.
From these comparisons we can extract a few useful conclusions. Perhaps most importantly, we see that the non-linear corrections to the solution are truly small everywhere in the domain. This is sensible because the scalar field itself is small and the values of that we explore are small, so the corrections to GR can be treated perturbatively. Another interesting observation is that the largest deviations from GR manifest somewhere in the middle of the domain. This is in part due to the boundary conditions: at spatial infinity the metric must be asymptotically flat so the GR deformation must vanish at a suitable fall-off rate; near the horizon, the GR deformation must be regular, but due to its asymptotic form near the horizon, it must also vanish there. Our choice of quantities to compare also plays an important role. Note that we compare metric coefficients at different physical locations, as at the horizon by definition and the horizon can be at a different location in the two metrics we are comparing. It is also worth pointing out that observables may depend on derivatives of the metric or even integrals of combinations of metric functions in the spatial domain and, hence, agreement between metric coefficients in part of the domain does not necessarily imply that observables in linear sGB will be close to those in GR. We have seen this already in Fig. 1, and will return to this point in Sec. V.
IV.3 Einstein-dilaton-Gauss-Bonnet Gravity
Let us now consider the case of an exponential coupling function. The resulting field equations are Eqs. (82) and (85) with . In this case, a perturbative solution in small coupling does not exist, but we can still compare any numerical solutions to the small coupling approximation for the linear theory of the previous section. We find a numerical solution using the computational infrastructure of Sec. II, with the same choices for the grid spacing, Newton polynomial order, etc as in the previous subsection.
The top and middle panels of Fig. 8 show the EdGB corrections to the metric components as a function of the compactified coordinate for different choices of the coupling constant . In contrast to the linear sGB results of Fig. 6, in EdGB the corrections to the metric are immediately noticeable. When comparing the difference between the analytic and numerical solutions in the bottom panel of Fig. 6, we see they range from to , whereas the difference in the linear coupling case ranged from to . The results of this comparison are not very surprising because the analytic perturbative solution was found exclusively in the linear coupling case.
These same features can also be seen in the scalar field solution, as shown on the left panel of Fig. 9. The differences, shown on the right panel of this figure, range from to , while for comparison the differences in the linear coupling case ranged from to .
From these comparisons we can also extract a few useful conclusions. Unlike in the linear sGB case, we see here that the analytic perturbative solution found in that case does not agree with the fully nonlinear solution. This is expected. Perturbing in the coupling constant and perturbing in the scalar will formally yield the same field equations. This can be understood by the fact that one can always absorb the coupling constant in a scalar field redefinition. The range of validity of the two expansions, and hence their physical interpretation, can be different. Nonetheless, this argument clearly implies that the perturbative solution will provide a better approximation to the linear coupling than the exponential one. Finally, as in the linear sGB case, we find that the metric corrections vanish near the horizon and compactified infinity, and both deviations asymptote to each other near infinity, but derivatives of the metric potentials can be large.
V Properties of Solution
In this section we explore the physical properties of the numerical solutions found in the previous section. We begin by finding analytical models that we fit to the data to provide accurate, closed-form expressions that allow for the rapid computation of physical observables. We then use these fitted models and the numerical results to calculate the location of the innermost stable circular orbit and the light ring, and compare them with each other and with the analytical perturbative solutions.
V.1 Fitting Function
In the compactified coordinate system introduced in Eq. (77), the full nonlinear solutions can be expressed as
[TABLE]
where , , and are the analytical perturbed solutions of Eqs. (94) in the compactified coordinates and , , are nonlinear corrections that we wish to find.
Using the analytical perturbed solutions as an ansatz, we propose best fit models for the non-linear corrections of the form
[TABLE]
where we have set . We then fit these models to our numerical solutions to determine the constants on the grid domain and . For the analytical perturbative solution is indistinguishable from the full nonlinear solution to our specified tolerance. For we find pathologies in the numerical solution that we will describe in Sec. V.4. The fitting order of our models is determined by systematically increasing the polynomial order of each function until the residual between the numerical solution and the model saturates. Some of the best-fit coefficients are included in Appendix A, but they are available in a Mathematica file upon request.
A comparison between the numerical data and the analytical fitted models is presented in Fig. 10. Here we plot the field components (properly rescaled to fit all in the same figure) in each top plot and their corresponding residuals as a function of the compactified coordinate in each bottom plot for both the linear sGB and EdGB solutions and for two coupling values of . We find that the residual between the numerical data and the analytical fitted model is always below our desired tolerance on the numerical solution of . With this caveat, the fitted models can be treated as “exact” for practical applications (up to the accuracy mentioned above).
V.2 ISCO
The inner edge of accretion disks around black holes are typically characterized by the innermost stable circular orbit (ISCO) of massive test-particles Abramowicz and Fragile (2013). For our spherically symmetric ansatz, the marginal stable circular orbits are determined by solving for circular timelike geodesics for massive test-particles Barausse et al. (2011); Rezzolla and Zhidenko (2014) orbiting the black hole. This is equivalent to requiring conditions on the effective potential, , where this potential is given by
[TABLE]
where and are the energy and angular momentum per unit mass (for massive particles) respectively, defined from the conserved quantities corresponding to the temporal and azimuthal Killing vectors in a stationary and spherically symmetric spacetime.
As in Newtonian gravity, a stable or an unstable circular orbit occurs at local minima or maxima of the effective potential, such that and . In a Schwarzschild metric, there is both an unstable and a stable circular orbit, such that the unstable orbit is closer to the horizon, and the distance between these orbits is determined by the angular momentum . The innermost stable circular orbit is equivalent to finding the value of the angular momentum where these stable and unstable orbits coincide (because the unstable orbit will always be closer to the horizon than the stable orbit). By analogy, requiring these orbits coincide is equivalent to finding the saddle points of the effective potential, which are located at . By combining these three conditions, we find the generalized equation
[TABLE]
the solutions (there may be multiple) of which give the locations of the marginal stable circular orbits of the spacetime. The smallest of these solutions is identified as the ISCO.
In our compactified isotropic coordinates, the location of the ISCO in GR is which corresponds to the familiar when transformed to Schwarzschild coordinates. In sGB gravity, the ISCO location is shifted from this Schwarzschild value. We can find the ISCO shift using the perturbative solution of Eq. (94) to find
[TABLE]
which is identical to that of Yunes and Stein (2011) when converted to Schwarzschild coordinates. We can also find the ISCO shift for the numerical metric solving Eq. (99) with a Newton-Raphson algorithm. By taking the location of the ISCO in GR as our initial guess, we ensure that the converged root is the desired root, as we expect deviations to be comparably small.
We presented these results already on the left panel of Fig. 2 in Sec. I, where we saw that the ISCO shift is typically smaller than . We also saw there that the shift computed with the analytic perturbative solution in the linear sGB case [Eq. (100)] agrees well the shift computed with the numerical solution in linear sGB but disagrees in EdGB. Interestingly, the shift computed with the fitted models agree extremely well with the numerical solution in both cases.
V.3 Light Ring
The light ring or photon sphere is the surface generated by all unstable circular null geodesics of photons. The location of the light ring around black holes is important for observations with the Event Horizon Telescope Doeleman et al. (2009), which is imaging the black hole shadow of Sagittarius A*, i.e. the electromagnetically dark region caused by photons that cross the light ring and fall into the event horizon. Future observations of black hole shadows may be able to place constraints on the location of the light ring in other quadratic gravity theories Ayzenberg and Yunes (2018).
Similar to the ISCO calculation, the light ring can be found by requiring certain conditions on the effective potential. For massless particles, there is only a single unstable circular orbit and no stable circular orbits. Thus to find the unstable circular orbit, we need only require and , which leads to the equation
[TABLE]
As before, the smallest solution to this equation returns the location of the light ring around the black hole.
The location of the light ring in GR is simply , which reduces to in Schwarzschild coordinates. As in the ISCO case, the location of the light ring is shifted in sGB gravity. We can calculate this shift with the perturbative analytic solution to find
[TABLE]
a result that to the best of our knowledge had not appeared in the literature previously. In the full non-linear case, we must solve Eq. (101) numerically using a Newton-Raphson method with the GR shift as our initial guess.
The light-ring shift was already presented on the right panel of Fig. 2. The shift is comparable to the shift of the ISCO, typically smaller than . Interestingly, we do find a noticeable disagreement between the analytic perturbative solution and the linear coupling case for higher values of , which was not present in the other calculated observables. The comparison between these two solutions in Fig. 6 shows that the largest differences between them occur closer to the horizon. Therefore this difference is larger in the region around the location of the light ring () than in both the region around the location of the ISCO () and asymptotically far away (). Thus it is expected that an observable calculated in this region should have a comparatively magnified discrepancy between the analytic perturbative solution and the numerical linear sGB solution. We also find that the fitted models agree extremely well with the numerical solutions for both coupling cases.
V.4 Naked Singularity
Spherically symmetric black holes in sGB gravity have been shown to possess a minimum size for a given in both EdGB Kanti et al. (1996) and linear sGB Sotiriou and Zhou (2014b). This results from a consistency condition on the field equations, obtained by requiring the scalar field to be regular on the horizon. Physically, as one increases the coupling strength , the location of the curvature singularity inside the horizon grows while the location of the event horizon shrinks, until at some critical value of the two coincide. For values of larger than this critical value, the curvature singularity is outside the event horizon, leading to a naked singularity. Requiring that the latter do not exist yields a maximum value of the coupling strength (and a minimum size of the event horizon) for which sGB black hole solutions can exist.
Our numerical solutions confirm these results. In our numerical calculations, we impose boundary conditions on the horizon using compactified coordinates at . The transformation from Schwarzschild to compactified coordinates absorbs the horizon shift, so that physical the horizon is always located at in our numerical grid. This then implies that there is a maximum value of above which black hole solutions should not exist in our numerical code. Indeed, we find that for values of larger than roughly on our grid of points, our code ceases to converge to the required tolerance. This is because a curvature singularity sufficiently near (or inside) the computational domain induces large errors in the Newton polynomial representation of the solution near the horizon boundary, which then propagates through the entire domain in each iteration, preventing the algorithm from converging.
This is indeed what we see in our numerical calculations as we increase : the estimated discretization error on the horizon () begins to grow as the location of the curvature singularity approaches the event horizon boundary. Eventually, the curvature singularity is close enough to the horizon radius that the discretization error near the horizon becomes too large for the specified tolerance. In order to ensure that these results are not a numerical artifact, we implemented adaptive step size refinement on the computational grid that is triggered if the discretization error becomes too large. Even with this adaptive measure in place, the discretization error still grows near the horizon for sufficiently large , preventing the code from converging.
In order to further support these conclusions, we have computed the Gauss-Bonnet curvature invariant for different values of in both the linear sGB and EdGB theories. Figure 11 shows this invariant as a function of the compactified coordinate . Observe that as (near the horizon) the curvature invariant begins to grow to ever larger values as is increased. Observe that as increases, the correction to at orders larger than begin to quickly approach the correction at where this effect is not present.
VI Conclusions
We have here developed a new numerical framework to solve for stationary and spherically-symmetric spacetimes that represent black holes in a wide class of modified theories of gravity. This framework uses a Newton polynomial representation for the discretized functions, it then recasts the differential system as a linear algebra problem, and then solves the latter through a relaxed Newton-Raphson iterative method. Through the successive minimization of the residual, our framework is capable of controlling the maximum error in the final numerical solution. We have validated this framework through a toy problem consisting of a simple differential equation, through the Schwarzschild metric and by investigating black holes in sGB gravity.
With the sGB solutions at hand, we then investigated a series of physical properties of these spacetimes. First, we verified that the differences between exact numerical solutions and analytic perturbative solutions are very small when the coupling is linear. We then also verified that the exact numerical solutions in linear sGB differ quite significantly from those in EdGB. These similarities and differences manifest themselves not only through the metric tensor, but also through physical observables like the ADM mass, the scalar charge, the location of the ISCO, and that of the light-ring. We finally verified that sGB black holes do not exist beyond a critical value of the sGB coupling, as beyond this value a naked curvature singularity arises.
We then concluded our analysis by developing analytic fitting functions for the numerical solutions. These fitting functions are constructed through a combination of controlling factors (inspired by the analytic perturbative solution) and polynomials in the compactified coordinate. We verified that the fitting functions agree with the numerical solutions up to the numerical error in the latter. We then computed physical observables, like the ADM mass, the location of the ISCO and that of the light ring with the fitting functions and found excellent agreement between these results and those obtained from the exact numerical solutions.
The work we did here now opens the door to several further studies. The fitting functions described above, for example, could be used as the analytical background on which to study polar and axial perturbations. Such perturbations would then reveal the quasi-normal mode spectrum of sGB black holes for arbitrary values of the coupling. The quasi-normal mode frequencies could then be used to carry out spectroscopic tests of General Relativity with gravitational wave observations of merging black holes (provided the merger remnant has very small final spin).
Another interesting direction for future research is the extension of the methods developed here to axisymmetric black holes. The computational infrastructure we presented here is easily extendable to this case. The system of equations of course becomes more complicated not only because new metric functions must be solved for, but also because these functions will depend on both radius and polar angle. We have already extended the work presented in this paper to a two-dimensional grid that is capable of solving for the Kerr metric in General Relativity, and thus, we expect that extensions to modified gravity at this point should be straightforward.
Once such solutions are found, the fitting methodology developed here could be implemented in the axisymmetric case to find fully analytic approximations for all components of the metric tensor. Such a solution could then be used once more as a background on which to study the evolution of perturbations. The quasi-normal spectrum of these perturbations could then be used to place constraints on a variety of modified gravity theories through the future observations of gravitational wave ringdown modes with advanced detectors.
Acknowledgements.
We would like to acknowledge Hector Okada-Da Silva for useful comments and suggestions. A. S. and N. Y. would like to acknowledge support from the NSF CAREER grants PHY-1250636 and PHY-1759615, as well as NASA grants NNX16AB98G and 80NSSC17M0041. T. P. S. acknowledges partial support from the STFC Consolidated Grant No. ST/P000703/1 and networking support from the COST Action GWverse CA16104.
Appendix A
Here are a few fitting coefficients of Eq. (97). The full tables are available by request in Mathematica.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abbott et al. (2016 a) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116 , 061102 (2016 a), URL https://link.aps.org/doi/10.1103/Phys Rev Lett.116.061102 .
- 2Abbott et al. (2016 b) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116 , 241103 (2016 b), URL https://link.aps.org/doi/10.1103/Phys Rev Lett.116.241103 .
- 3Abbott et al. (2017 a) B. P. Abbott et al. (LIGO Scientific and Virgo Collaboration), Phys. Rev. Lett. 118 , 221101 (2017 a), URL https://link.aps.org/doi/10.1103/Phys Rev Lett.118.221101 .
- 4Abbott et al. (2017 b) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), The Astrophysical Journal Letters 851 , L 35 (2017 b), URL http://stacks.iop.org/2041-8205/851/i=2/a=L 35 .
- 5Abbott et al. (2017 c) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119 , 141101 (2017 c), URL https://link.aps.org/doi/10.1103/Phys Rev Lett.119.141101 .
- 6Abbott et al. (2017 d) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119 , 161101 (2017 d), URL https://link.aps.org/doi/10.1103/Phys Rev Lett.119.161101 .
- 7Yunes and Siemens (2013) N. Yunes and X. Siemens, Living Reviews in Relativity 16 , 9 (2013), ISSN 1433-8351, URL https://doi.org/10.12942/lrr-2013-9 . · doi ↗
- 8Will (2014) C. M. Will, Living Reviews in Relativity 17 , 4 (2014), ISSN 1433-8351, URL https://doi.org/10.12942/lrr-2014-4 . · doi ↗