GRFolres: A code for modified gravity simulations in strong gravity

The


Summary
Gravitational waves (GWs) are generated by the mergers of dense, compact objects like black holes (BHs) and neutron stars (NSs).This provides an opportunity to study the strong field, highly dynamical regime of Einstein's theory of general relativity (GR) at higher curvature scales than previous observations [1,2,3,4,5,6].It is possible that at such scales modifications to GR may start manifest.However, in order to detect such modifications, we need to understand what deviations could look like in theories beyond GR, in particular in the merger section of the signal in near equal mass binaries, which are key targets of the LIGO-Virgo-KAGRA network of detectors (and their future 3G successors).Such predictions necessitate the use of numerical relativity (NR), in which the (modified) equations of GR are evolved from an initial configuration several orbits before merger, through the merger period and the subsequent "ringdown", during which the gravitational wave signal can be extracted near the computational boundary.
Current waveforms are tested for consistency with GR by measuring parameterised deviations to the merger, inspiral and ringdown phases [7,8,9,10,11], and not by comparison to any particular theories.If we obtain predictions for specific models, we can check whether such parameterised deviations are well-motivated and consistent in alternative theories of gravity [1,12,13,14,15,16,17], and quantify our ability to constrain model parameters using GW observations.
There are many ways to modify GR, one of the simplest being to couple an additional scalar degree of freedom to gravity, which may (if certain conditions are satisfied) result in so-called "hairy" stationary black hole solutions; that is, black holes with a stable, non trivial configuration of the scalar field around them (see [18] for a review).An example of this is the class of Horndeski models [19].Cubic Horndeski theories have been studied in [20,21] and an implementation of this is included in GRFolres.Another more general example within the Horndeski models is the four-derivative scalar-tensor theory (4∂ST), which is the most general theory with up to fourth powers of the derivatives (but still second order equations of motion).Despite their relative simplicity, many models have lacked well-posed (and thus numerically stable) formulations until recently.
An important breakthrough was made in 2020 by Kovács and Reall, who showed that Horndeski theories are indeed well-posed in a modified version of the harmonic gauge [22,23] -a particular coordinate system already used in NR.Subsequently, several specific theories within these classes were probed in their highly dynamical and fully non-linear regimes [24,25,26,27].The extension of the results of [22,23] to the alternative "singularity avoiding" coordinates in [28,29,30] offers an alternative gauge in which to probe questions of hyperbolicity, and may even offer stability advantages for certain cases such as unequal mass ratios, as studied in [27].Numerical work on these theories is still in the early stages of development and many technical details on their numerical implementation need to be further investigated.Equally, many scientific questions, concerning our accurate understanding of binary black holes' phenomenology in alternative theories of gravity and their implications for tests of GR, also remain unanswered.
The goal of GRFolres is to meet this need for further research, and to provide a model code to help others develop and test their own implementations.The code is based on the publicly available NR code GRChombo [31,32], which itself uses the open source Chombo framework [33] for solving partial differential equations (PDEs).
In the following sections, we discuss the key features, motivations, and applications of the code.

Key features
GRFolres inherits many of the features of GRChombo and Chombo.Here we list the key features.
• Stable gauge evolution -The code implements the modified moving puncture gauge that ensures a well-posed evolution in the weak coupling regime, as proposed in [28].The precise form of the gauge and its parameters can be changed and the standard moving puncture gauge is safely recovered by setting certain parameters to zero.
• Modified gravity theories -The currently available theories in the code are 4∂ST and cubic Horndeski.The code is templated over the theory (in the same way that GRChombo is templated over a matter class) so that it can easily be changed without major code modifications.The code also provides an implementation of 4∂ST without backreaction onto the metric (but including the possibility of using the new gauge), to enable comparison with previous works using the decoupling limit approximation.
• Accuracy -The fields are evolved with a 4th order Runge-Kutta time integration and their derivatives calculated with the same finite difference stencils used in GRChombo (4th and 6th order are currently available).
• Initial Conditions -The current examples use solutions that approximately or trivially solve the modified energy and momentum constraints of the theory.An elliptic solver for more general configurations is under development, using a modified CTTK formalism [34,35].
• Diagnostics -GRFolres has routines for monitoring the constraint violation and calculating the energy densities associated with the different scalar terms in the action, as discussed in [28,29,30].Other diagnostics can be added as required.We also extract data for the tensor and scalar gravitational waveforms.
• C++ class structure -Following the structure of GRChombo, the GRFolres code is also written in C++ and uses object oriented programming (OOP) and templating.
• Parallelism -GRChombo uses hybrid OpenMP/MPI parallelism with explicit vectorisation of the evolution equations via intrinsics, and is AVX-512 compliant.
• Adaptive Mesh Refinement -The code inherits the flexible AMR grid structure of Chombo, which provides Berger-Oliger style [36] AMR with block-structured Berger-Rigoutsos grid generation [37].Depending on the problem, the user may specify the refinement to be triggered by the additional degrees of freedom, i.e. the scalar field, or those of the metric tensor.

Statement of Need
As far as we are aware there is currently no other publicly available code that implements the 4∂ST theory of modified gravity or the cubic Horndeski theory in (3+1)-dimensional numerical relativity.
There is at least one private code, based on the PAMR/AMRD and HAD [38,39] infrastructure, that was used in the first works to successfully implement the modified general harmonic gauge for 4∂ST [24,25,26,27].Since this code uses a Generalised Harmonic Coordinates (GHC) formulation, it necessitates excision of the interior of black holes, which can be difficult to implement in practice.As a consequence, many groups in the numerical relativity community have opted to use singularity avoiding coordinates such as the BSSN [40,41,42], Z4C [43,44] or CCZ4 [45,46] formulations in the puncture gauge [47,48], which do not require the excision of the interior of black holes from the computational domain.In GRFolres, we use the results of [28,29,30] to extend the well-posed formulations of modified gravity to singularity avoiding coordinates.This provides an alternative gauge to the modified GHC one used by other groups.Not only does this provide a valuable comparison to their work, but also eliminates the need for excision.
In spherical symmetry several codes have been developed that implement Einstein-scalar-Gauss-Bonnet (a subset of the 4∂ST theory that we include as an example in GRFolres).In particular, using the NRPy framework (http://astro.phys.wvu.edu/bhathome)[66] in [63], the private code of Ripley & Pretorius in [67,68,69,70,71], and a modification of the GR1D code [72,73] for calculating core collapse in Einstein-scalar-Gauss-Bonnet in [74].There is also the fully nonlinear spherical code developed in [75,76].Spherical codes provide a useful testing ground in which coordinate ambiguities can be avoided [67], and a well posed formulation is easier to obtain, but they lack the generality required to study objects with angular momentum, or binary mergers.

Research projects to date using GRFolres
So far the code has been used to study a range of fundamental physics problems, as listed here.
• The test field case was used in [49] to model the scalar waves produced during the ringdown stage of binary black hole coalescence in Einstein-scalar-Gauss-Bonnet, and quantify the extent to which current and future gravitational wave detectors could observe the spectrum of scalar radiation emitted.
• The regime of validity of effective field theory in collapse and binary evolutions in cubic Horndeski theories was studied in [20,21].It was found that the mismatch of the gravitational wave strain can be as large as 10%-13% in the Advanced LIGO mass range for such theories.
• In the work [28], the code was developed and tested, with waveforms for shift-symmetric theories of Einsteinscalar-Gauss-Bonnet gravity produced for equal mass binaries, as illustrated in Fig. 3.
• In the work [29], the studies were extended to binary mergers in theories with spin-induced scalarisation.The clouds formed are dumbbell-like in shape, as illustrated in Fig. 4.
• In the work [30], the dependence of the conditions for hyperbolicity and weak coupling were studied for spin-induced scalarisation, and the critical thresholds found for a number of cases, as illustrated in Fig. 5.  GB /M 2 = 0.0005, g 2 /M 2 = 0.001  Figure 5: The time evolution of the determinant of the effective metric in a case of spin-induced scalarisation.When the determinant is negative (in black) outside the apparent horizon (depicted with a dashed white line), the theory has become ill-posed.Taken from [30].

Figure 1 :
Figure 1: Contour plot of network signal-to-noise ratio (SNR) for the scalar ringdown of a binary black hole (BBH) in Einstein-scalar-Gauss-Bonnet gravity at 1 Gpc as observed by the Virgo, Livingston and Hanford network of detectors at design sensitivity.Taken from [49].

Figure 2 :
Figure 2: Energy density (in blue) of the scalar field surrounding the binary black holes for the Horndeski theory at a representative instant of time during the inspiral phase.The apparent horizon of the black holes is shown in orange.The region where the weak coupling conditions are larger than one is depicted in brown.Taken from [21].

Figure 4 :
Figure 4: The time evolution of the density of the scalar cloud that develops in Einstein-scalar-Gauss-Bonnet gravity with an exponential coupling, resulting in spin-induced scalarisation.Taken from [29].