Solving the initial conditions problem for modified gravity theories

,


I. INTRODUCTION
Recent breakthroughs in well posed formulations [1,2] have resulted in an expansion of the class of modified theories of gravity to which numerical relativity (NR)numerical simulations that solve the Einstein Equations as a time evolution problem -can be applied.This has allowed for the simulation of strong gravity spacetimes in these theories, including the fully non-linear backreaction of the additional curvature terms onto the metric [3][4][5][6][7][8], building on previous works that neglected such effects [9][10][11][12][13][14][15][16][17][18].Such simulations are of particular interest for the construction of binary merger waveforms, but also have relevance to more general questions about the hyperbolicity of such theories in the strongly coupled regime [3][4][5][19][20][21], and potentially in early Universe cosmology [22].However, in works to date in 3+1 spacetime dimensions, the challenge of satisfying the more complicated constraint equations has limited the physical scenarios that can be investigated.Most have begun with a trivial solution that satisfies the regular constraint equations of general relativity (GR) and set the additional scalar degree of freedom to zero (it is then evolved to a non trivial configuration over time) [3,[5][6][7].Alternatively, studies of spin-induced scalarisation that require a scalar "seed" have tolerated a certain level of constraint violation in the initial data arising from a non trivial profile of the scalar, and relied on constraint damping to remove it during the initial stages of the evolution [4,8,20].
The Arnowitt-Deser-Misner (ADM) [23] decomposition of the Einstein Field Equations provides the basis for NR simulations, posing the system as a Cauchy problem, with a set of evolution and constraint equations for the 3-dimensional spatial metric γ ij and the extrinsic curva-ture, K ij , of the spatial slices.The ADM approach gives rise to four independent constraint equations, which must be satisfied on each time slice, and in particular in the initial data that starts an NR simulation.These constraints are given by where K = γ ij K ij is the trace of the extrinsic curvature (also known as the mean curvature), and R and D i are the Ricci scalar and covariant derivative associated with the 3-dimensional spatial metric.The energy density ρ and momentum densities S i are functions of the matter fields and their derivatives.These can be derived directly from the matter action as projections of the stress-energy tensor T µν .In the case of modified gravity theories like EsGB, the additional curvature terms can be treated as further effective matter sources that depend on higher order derivatives of the other metric variables.
Including the lapse α and shift β i , there are 16 components that must be specified initially, but only four constraint equations.Whilst four of the extra components relate to physical degrees of freedom, the remainder are gauge.Values that represent free data and those that are fully determined by the physical scenario are not easy to separate, and so some must be chosen somewhat arbitrarily.The result is that there will often be a large number of possible ways to set the free data that appear to still meet the physical requirements.However, poor choices can result in uniqueness and existence problems in the initial constraints that prevent unique solutions being found.
Many different approaches have been developed for solving the constraint equations (for reviews see the stan-dard NR texts [24][25][26]).These approaches vary in the degrees of freedom that have their values imposed, and those that are solved for.The metric and extrinsic curvature can be decomposed in various ways (e.g., through a conformal decomposition), allowing for the freely chosen variables to be more closely identified with a particular physical interpretation.This can also leave the equations in a form more amenable to numerical solutions.
One such approach is known as the Conformal-Transverse-Traceless (CTT) method.In this approach the extrinsic curvature is decomposed into its trace K and a traceless tensor A ij .The 3-metric γ ij is decomposed into a conformal metric with determinant one γij , and a conformal factor ψ, as Āij , and the conformal Āij is further decomposed into a transverse-traceless part ĀT T ij and a vector potential W i (see equation 3).The conformal metric γij , the mean curvature K, and the transverse-traceless extrinsic curvature ĀT T ij are then imposed to have some values (usually trivially those of a flat spacetime), and the constraints are used to solve for the conformal factor ψ and the vector potential W i .Āij can then be reconstructed as, In this method, the term in the energy density ρ is generically problematic, due to the nature of the resulting elliptic equation for ψ.It is normally only possible to specify a conformally-related density ρ, as this allows for the term in ρ to appear with the "right sign" to guarantee unique solutions of the elliptic equation for ψ [27].
A modification to the CTT technique, known as CTTK, was recently proposed to address this limitation, which is more acute for problems involving fundamental fields as opposed to fluids [28].In this approach the variables are decomposed in the same way, but the elliptic equation for ψ is split into an algebraic equation for K in terms of ρ, and ψ is then only required to satisfy Laplace's equation with a source term in Āij .The conformal metric γij and ĀT T ij still need to be chosen, but the mean curvature K is now solved for as part of the algorithm and this simplifies the solution for ψ.Overall the method was found to be highly robust, with unique solutions being found by the solver over a wide range of scenarios.
In this work, the CTTK method has been modified to solve the full constraint equations of the most general four derivative, parity invariant scalar-tensor theory of gravity (4∂ST ).It has been shown that the CTT formulation of the elliptic equations for these theories has unique solutions in the weak coupling limit [29], and one can expect this to carry over to the CTTK case.To achieve this, as suggested in [29], the contributions from the higher-derivative curvature terms in the action are treated as additional sources that are split between the available degrees of freedom.The CTTK technique is modified in two different ways, depending on the boundary conditions of the problem at hand, to ensure that solutions exist.We describe the methods in more detail below, and demonstrate their effectiveness in practise by showing that the solutions converge as expected.
We follow the conventions in Wald's book [30].Greek letters µ, ν . . .denote spacetime indices and they run from 0 to 3; Latin letters i, j, . . .denote indices on the spatial hypersurfaces and they run from 1 to 3. We set G = c = 1.

A. Additional terms
In 4∂ST gravity, the Einstein-Hilbert action is modified to include a scalar field coupled non-trivially to gravity, through the Gauss-Bonnet term L GB , and g 2 (ϕ) and λ(ϕ) are smooth functions of ϕ.This is the most general parity-invariant scalar-tensor theory of gravity that includes up to four derivatives of the metric, and is an example of the wider class of Horndeski theories [31].The local magnitude of the coupling λ(ϕ) controls the deviation of the solutions from the GR case, with the g 2 (ϕ) term being generally subdominant for similar coupling constants and the same field amplitude.Deviations are also amplified in regions of high curvature (i.e., with a larger Gauss-Bonnet invariant).The effects of these additional terms can be included in ρ and S i , so that the constraints remain as they are in Equations 1 and 2 but with an effective ρ and S i given by Here ρ GB and S GB i include higher derivative contributions from the scalar and tensor metric variables, and ρ X and J X i contain the contributions of the g 2 (ϕ) term.Explicitly: with where Π is the conjugate momentum of ϕ, N i is the GR momentum constraint, Ω = γ ij Ω ij , and Ω i and Ω ij come from the 3 + 1 decomposition of the Weyl tensor.
The contributions from the kinetic and potential terms are given by the usual minimally coupled, real scalar terms, As mentioned in the introduction, previous work in such theories has required that ρ GB + ρ X ≪ ρ SF and S GB i + S X i ≪ S SF i everywhere.This reduces the problem to (approximately) satisfying the regular GR constraint equations, and still allows the solution to evolve away from that of GR during the evolution.However, here we are able to treat the deviations from GR fully nonperturbatively.

B. Choice of components to solve for
As discussed in the introduction (and in more detail in [28]), in contrast to the CTT method that imposes a spatially constant value of the mean curvature K, and solves for the conformal factor ψ = (det γ)1/12 , the CTTK method allows both K and ψ to have a spatially varying profile.The momentum constraint (now with an additional source term from the spatial variation of K) is solved as in the CTT method, for the vector potential W i , which is further decomposed into a scalar U and a vector V i .This additional flexibility in K can be used to absorb some of the problematic matter source terms, resulting in a more robust method for general field configurations.The same approach will be used here, but we will see that in general the 4∂ST terms are better absorbed into the elliptic equation for ψ, rather than in the mean curvature K.
The transverse-traceless part of the extrinsic curvature, ĀT T ij , is usually set to zero, which roughly corresponds to a spacetime containing no gravitational waves.The conformal metric is also assumed to be δ ij for simplicity (i.e., conformal flatness is assumed).It should be possible to apply the same techniques with more general choices of ĀT T ij and γij , such as those that are required for highly spinning black hole initial data, which would simply result in additional source terms in the equations.For now we maintain these choices for simplicity.

C. Black hole spacetimes
With the 4∂ST terms included and an appropriate form chosen for U (see the appendix of [28] for a discussion), we write the Hamiltonian and momentum constraints in CTTK as follows The GR case with ρ = ρ SF and S i = S SF i is derived in [28] -in what follows we explain the motivation for the positioning of the additional terms.
In regions of physical relevance, we can demand from an effective field theory (EFT) perspective that ρ GB + ρ X < ρ SF .However, in some regions of high curvature this may not hold.Moreover, the value of ρ GB + ρ X is not positive definite, and in practise for most cases of interest it varies between positive and negative regions in different parts of the spatial slice.For this reason, ρ SF and ρ GB + ρ X have been separated between equations ( 12) and (13).By including only the positive definite part in Eq. ( 12), we avoid the possibility of K 2 having a negative source, even in regions where the additional curvature contributions dominate.However, since it is also not negative definite, the ρ GB + ρ X term in equation ( 13) is likely to appear with the "wrong sign" in a linearised expansion of ψ in certain regions, which violates the conditions of the maximum principle and removes the guarantee of unique solutions [27].In all our test cases (see Section III), this has not caused issues, and no further modification has been necessary for convergence to a solution.We speculate that this may be because we also include the contribution from Āij Āij in the elliptic equation for ψ, which tends to stabilise the solutions 1 .
If the spacetime contains one or more black holes, divergences will appear in ψ and Āij .This can be dealt with by decomposing Āij (and therefore U and V i ) into a black hole part, Ābh ij , which contains the divergences, and a regular part, Āreg ij , which is solved for.The conformal factor ψ can also be decomposed into a sum of isolated black-hole solutions, ψ bh , and the regularised remainder, ψ reg .The Poisson equations for ψ and V i can then be solved directly, or by linearising around the previous solution and solving for the error.For example, the results given in Section III A are calculated by fully solving for V i at each step, and iterating for ψ as ψ = ψ 0 + δψ.
The forms of A bh ij and ψ bh are given in [28], and are those proposed by Bowen and York in [32,33].

D. Cosmological spacetimes
For simulations of cosmological spacetimes, periodic boundary conditions are often used (see [34][35][36][37][38][39][40][41][42][43][44], although other approaches are possible [45]).In the CTTK method without a 4∂ST term, the method for solving the constraints with periodic boundaries is very similar to that for black hole spacetimes.Again, the constraints reduce to Equations ( 12) to (14).However, the elliptic equations now impose a set of integrability conditions, as discussed in [39,46].In cosmological spacetimes, the differential operator in equation ( 14) has a non-trivial kernel.Therefore, by the Fredholm alternative, solutions to this equation (which are necessarily non-unique) will only exist if the source is in the orthogonal complement of the kerneli.e., in the adjoint.Since the kernel includes constants, it is necessary (but not sufficient) for the right hand side of the Poisson equation to equal zero when integrated over the entire grid.A similar requirement applies to the right hand side of equation ( 13) when it is treated as a constant source.With no 4∂ST term, this simply corresponds to requiring a periodic distribution of the scalar field, and ensuring that there is no net momentum flux through the box in any direction.However, with λ(ϕ) ̸ = 0, this is no longer sufficient, as there is no guarantee of the source terms including ρ GB , ρ X , S GB i and S X i averaging to zero across the grid, and no obvious way of choosing them such that this is always the case.
In simple cases where only the elliptic equation for ψ is problematic (e.g., with a simple sinusoidal profile for ϕ and no conjugate momentum Π), this can be solved by further dividing K into two parts, one constant and one spatially varying, which we call respectively K C and K GR .The 'GR' part of K satisfies equation (12), as in the GR case, sourced by ρ SF , that is The constant part of K, K C , can then be used to compensate for any violation of the integrability condition for ψ.This is achieved by setting K C to a value satisfying the condition As a result of these choices, the equation for ψ becomes Equation ( 16) will always have real solutions, unless ρ GB dominates over the combined contributions of ρ GR and Āij Āij when averaged across the grid2 .The integrability condition for the momentum constraint can also be spoiled by the presence of various nonlinear terms in S GB i .With a shift-symmetric or quadratic coupling and Π = 0, many of these terms are simplified.A sinusoidal profile for ϕ then gives a contribution to S GB i that satisfies this constraint.For a more general S GB i , removing the assumption of conformal flatness would provide another source in the constraints, potentially allowing for this contribution to be cancelled out by a judicious choice.In our tests below we restrict ourselves to showing that the method works in the simpler case, and leave more general conditions to future work.
Even once the source in equation ( 14) is in the adjoint and the right hand side of ( 17) integrates to zero, the solution to the elliptic equations at each step suffers from non uniqueness -the equations have multiple solutions, where solutions differ by the addition of a constant or linear term in the equation for the conformal factor3 , and one or more Killing vectors of the conformally flat metric in the case of V i (see [46]).This is addressed by solving for their values perturbatively, starting with an initial guess and solving for a small correction (δψ, δV i ) at each step.This perturbative treatment naturally generates a linear term in the equation for δψ that prevents the constant and linear modes from growing.The freedom in δV i can also be eliminated by adding a small linear coefficient to the Poisson equation for δV i .The addition of the conformal Killing vectors in V i is unimportant as they do not change the resulting value of the extrinsic curvature K ij [46].However, any non uniqueness in ψ has a physical consequence -its value at the start of the iteration picks out the final uniquely chosen solution, and this in turn determines physical properties of the field -e.g., the density of the gradient terms measured by normal observers.

III. TESTS
The methods described above for both black hole and cosmological spacetimes have been tested with a modified version of the CTTK solver used in [28].This solver is constructed using the open-source Chombo [47] code for finite difference solution of PDEs with Adaptive Mesh Refinement -in particular here we adapt their multigrid solver for elliptic equations.The results are imported FIG.1: Plots of the L 2 norm across the grid of the Hamiltonian constraint violation H and the momentum constraint violation M = √ M i M i , normalised to their initial values H 0 and M 0 , against non-linear iterations of the elliptic solver for the black hole spacetime with a dumbbell scalar field configuration shown in Fig. 2(b).We see that the solver converges to a good solution within ten iterations.The behavior is similar for the cosmological initial data.
into GRFolres [48], the modified version of the NR evolution code GRChombo [49,50], to check the constraint violation using the methods verified in [7,8].
Here we show typical results for both overall convergence to a solution and convergence to the true zeroconstraints solution with increasing resolution.The Chombo multigrid solver is designed to be second-order in all derivatives, so with the assumption that the errors are dominated by errors in the derivatives, the solver errors should also converge at that rate as the resolution is increased.This means that doubling the resolution should reduce the constraints by a factor of four4 .These tests have been conducted with a variety of coupling functions and potentials, although (as described above) the possible scalar field configurations are restricted in the case of periodic boundaries.Figure 1 shows the convergence of the constraint violations to zero (with a fixed resolution) as the number of non-linear iterations increases for the black hole spacetime described in Section III A. This shows good convergence to a solution of the full non-linear problem within ∼ 10 iterations.In practise, this global measure is rather crude and ignores the fact that in some regions with small errors (often nearer the boundaries in BH spacetimes) the solver takes longer to show good local convergence.More information can be gained by checking the spatial profiles.The convergence tests with increasing resolution for two particular cases are shown below in Figures 2 and 3.These plots show that the solver is consistently displaying the desired second-order convergence, but in the black hole case the solver was run for approximately 1000 non-linear iterations.Given that the solver is much less costly to run than the evolution code, such a high number of steps is not prohibitive -with three levels of refinement, this takes a few hours with ∼ 100 CPU cores on a typical computing cluster.In the cosmological case, where only one level is used, it takes under an hour.

A. Tests of black hole method
The black hole method described in Section II C has been tested with a dumbbell-shaped scalar field and momentum configuration around a central black hole with dimensionless spin coefficient a/M = 0.5 and a spin-axis along the z direction, where M is the total bare mass.The potential and coupling functions are both chosen to be quadratic, i.e., λ(ϕ) = 1 2 λ 2 GB ϕ 2 and V (ϕ) = 1 2 m 2 ϕ 2 , with m = 1 and λ GB /M = 1, a scalar field amplitude of 0.1, and a momentum amplitude of 0.01.In this test g 2 (ϕ) was set to zero -its impact on the solutions was found to be negligible.This configuration was inspired by the scalar field configuration that the field has been found to settle into in spin induced cases -see for example the plots in [8,20]-although in such cases the field is massless.We do not attempt to match the stationary configuration exactly since this is simply a proof of principle, and we have tested other spatial configurations that show similar results to the ones here.Figure 2a shows the amplitude of the scalar field across a slice of constant y-coordinate, passing through the singularity.Figure 2b then shows the convergence as the grid spacing is halved.The results with a finer grid are approximately secondorder across the grid, other than at grid boundaries.In the Hamiltonian constraint the errors are not fully dom-  FIG.3: Scalar field configuration and convergence plots for a sinusoidal scalar field configuration and periodic boundary conditions.inated by the derivatives, which is likely to cause the small deviation from exact second-order convergence in the central region.The convergence towards a solution for a fixed resolution is also shown, in Figure 1.
This method has also been tested with initial data for a black hole binary.The same coupling functions and amplitudes as the dumbbell test were used for two black holes with Bowen-York masses m 1,2 = 0.5M , momenta |P | = 2M perpendicular to their separation of 12M , and dimensionless spin parameters a/M = 0.6.In this case we set g 2 (ϕ)/M 2 = 1, although as above this only has a small impact on the solutions.A Gaussian scalar field with momentum was imposed over each black hole, as shown by the contours in Figure 4, along with the change in the conformal factor over that of the bare punctures.

B. Tests of cosmological spacetimes method
A similar test was conducted with periodic boundary conditions, and a sinusoidal profile for ϕ in each direction.The same coupling and potential functions were used, with m = 0.5 and the other parameters unchanged.As explained above, we needed to set the momentum of the field Π to be zero in this case.Figure 3 shows the scalar field configuration and tests of convergence with increasing resolution in the periodic case, and also displays second-order convergence across the grid.The speed of convergence with non-linear iterations is similar to that in the black hole case.

IV. DISCUSSION AND FUTURE WORK
In this work the CTTK method was adapted to solve the constraint equations of 4∂ST gravity, treating the additional curvature terms as another source to the Poisson equations of the GR case, as suggested in [29].Whilst such a treatment is a sensible first guess given that the new terms should be small in an effective theory (and given that in the weakly coupled regime a unique solution has been shown to exist [29]), it is far from clear that such an approach will work in practise given the high non-linearity of the problem, especially as the coupling is made large.We have demonstrated that this is the case, and that the method is robust, provided cer-tain choices are made about the split of the new terms between the available degrees of freedom.In fact, the method appears to be robust up to the strongly coupled regime of the theory, with coupling parameters of order 1.Beyond this regime the well-posedness of the evolution is no longer guaranteed, and in practise it will break down [3][4][5][19][20][21]; therefore the solutions beyond this point are not of particular physical interest.
This adapted method has been tested for both black hole and cosmological spacetimes, and shows appropriate convergence in both cases.Techniques have also been described for guaranteeing the existence and uniqueness of solutions with either asymptotically flat or periodic boundary conditions.This is the first method that has been demonstrated to work for fully satisfying the constraint equations for generic initial data in 3+1 dimensions, and it therefore expands the possibilities for numerical investigation of these theories.Future applications of this technique may include simulations of black holes with non-trivial initial scalar field configurations, and tests of inflation in 4∂ST gravity with inhomogeneous initial conditions.It would also be useful to extend the method to permit non conformally flat spacetimes (which could address the issues with integrability of the momentum constraints that we encountered in the cosmological case) and to check that a similar approach works in the alternative extended conformal thin sandwich (XCTS) method [51,52], in which one solves also for the lapse and the shift functions to achieve a specific initial time evolution of the variables.XCTS is more widely used in practise than CTT and offers a number of advantages, especially for identifying equilibrium initial data.We anticipate that the approach developed here should work equally well in such cases, but we hope to see this demonstrated in future work by groups working with such solvers.
FIG.2: Scalar field configuration and convergence plots for a growing dumbbell scalar field configuration around a spinning black hole.

( a )
Sinusoidal scalar field configuration with periodic boundaries.Here the conjugate momentum of the field is set to zero.(b) Local values of the Hamiltonian and momentum constraint violations along a line for two different resolutions,showing second order convergence.

FIG. 4 :
FIG. 4: Correction to the conformal factor of the bare punctures for a black hole binary, with a superimposed contour plot showing the Gaussian scalar field over each black hole.