Numerical black hole solutions in modified gravity theories: Spherical symmetry case

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 and 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 nonlinear 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.