3+1 non-linear evolution of Ricci-coupled scalar-Gauss-Bonnet gravity
Daniela D. Doneva Llibert Aresté Saló Stoytcho S. Yazadjiev

	
	
	
	daniela.doneva@uni-tuebingen.de
	Theoretical Astrophysics, Eberhard Karls University of Tübingen, Tübingen 72076, Germany
	INRNE - Bulgarian Academy of Sciences, 1784  Sofia, Bulgaria
	
	
	
	l.arestesalo@qmul.ac.uk
	School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London, E1 4NS, United Kingdom
	
	
	yazad@phys.uni-sofia.bg
	Theoretical Astrophysics, Eberhard Karls University of Tübingen, Tübingen 72076, Germany
	Department of Theoretical Physics, Faculty of Physics, Sofia University, Sofia 1164, Bulgaria
	Institute of Mathematics and Informatics, 	Bulgarian Academy of Sciences, 	Acad. G. Bonchev St. 8, Sofia 1113, Bulgaria
	

	Scalar-Gauss-Bonnet (sGB) gravity with an additional coupling between the scalar field and the Ricci scalar exhibits very interesting properties related to black hole stability, evasion of binary pulsar constraints, and general relativity as a late-time cosmology attractor. Furthermore, it was demonstrated that a spherically symmetric collapse is well-posed for a wide range of parameters. In the present paper we examine further the well-posedness through 3+1 evolution of static and rotating black holes. We show that the evolution is indeed hyperbolic if the weak coupling condition is not severely violated. The loss of hyperbolicity is caused by the gravitational sector of the physical modes, thus it is not an artifact of  the gauge choice. We further seek to compare the Ricci-coupled sGB theory against the standard sGB gravity with additional terms in the Gauss-Bonnet coupling. We find strong similarities in terms of well-posedness, but we also point out important differences in the stationary solutions. As a byproduct, we show strong indications that stationary near-extremal scalarized black holes exist within the Ricci-coupled sGB theory, where the scalar field is sourced by the spacetime curvature rather than the black hole spin.


3+1 non-linear evolution of Ricci-coupled scalar-Gauss-Bonnet gravity
    Stoytcho S. Yazadjiev
    ^1Division of Space Research and Planetary Sciences, University of Bern, Switzerland

    ^2Max Planck Institute for Astronomy, Heidelberg, Germany

    

==============================================================================================================================================================



 

§ INTRODUCTION

	The rapid advance of gravitational wave detectors gives us confidence that in the next decades we will have the necessary observations to allow for precise tests of gravity <cit.>. The availability of accurate theoretical gravitational waveforms both in general relativity (GR) and its modifications will be crucial for the correct interpretation of the detected gravitational wave events. While the accuracy of the former is on a steady path of improvement, numerical relativity simulations beyond GR are still in the development phase. The reason behind that is twofold – first, the field equations typically get increasingly more complicated compared to GR when we start modifying the original GR action. The second more fundamental obstacle is the question of well-posedness, which is one of the main building blocks if we want to be able to perform time evolution. Even though strong hyperbolicity is proven in certain formulations of the 3+1 field equations in GR <cit.>, it is by no means guaranteed that this will remain true in modified gravity <cit.>. Whether this is sourced by an intrinsic problem of the theory or just a gauge change is required, is a question that awaits answer in a number of GR modifications.
	
	Our focus in the present paper will be the scalar-Gauss-Bonnet gravity (sGB). It provides an important playground for studying the possible deviations from GR one can have in an effective field theory of gravity while keeping the field equations of second order. It gained particular attention because it was proven that black holes with scalar hair can exist within this theory <cit.>, including spontaneously scalarized ones <cit.>. Even though evolution in spherical symmetry can be hyperbolic for a weak enough coupling <cit.>, solving the full 3+1 field equations is much more subtle with the standard harmonic gauge proven to be non-well-posed <cit.>. Interestingly, loss of hyperbolicity is observed even at the level of linear perturbations <cit.>.  
	
	An important breakthrough was the proof that a modified harmonic gauge leads to well-posed 3+1 field equations in the weak coupling regime, i.e. when the contributions of the Gauss-Bonnet term to the field equations are smaller than the two-derivative Einstein-scalar field terms  
 <cit.>. This regime is exactly where sGB gravity can be considered as a viable effective field theory. This eventually allowed the development of 3+1 numerical relativity codes <cit.> and an extension to a well-posed modified puncture gauge <cit.>.
	
	The search for a well-posed formulation of sGB gravity also eventually ignited the interest in alternative ways to address the problem. An interesting approach is to “fix” the field equations, which can be regarded as providing a weak completion of the considered effective field theory <cit.> and is inspired by the dissipative relativistic hydrodynamics <cit.>. Even though it seems promising, further development is needed to have a self-consistent and robust 3+1 evolution. Another approach is to modify the original sGB action and a natural extension is to add a coupling between the scalar field and the Ricci scalar coupling <cit.>. It was shown that it can lead to hyperbolic evolution in the case of a spherically symmetric collapse of a scalar field <cit.>. As a matter of fact, this theory has other interesting features such as linear stability of black holes that are otherwise unstable in the standard sGB gravity and the possibility to evade binary pulsar constraints for certain ranges of parameters <cit.>. Furthermore, it also cures some of the problems encountered when treating scalar-tensor theories in a cosmological set-up. In particular, adding the Ricci scalar coupling is a minimal model which succeeds in having GR as a cosmological late-time attractor <cit.>, which is otherwise not true <cit.>. However, other problems occurring at early times persist within this model <cit.> and one would need to add extra operators <cit.> in order to cure them.  
	
 In the present paper, we aim to explore further the Ricci-coupled sGB gravity by investigating the well-posedness of the theory in the 3+1 formulation of the field equations using the modified gauge proposed in <cit.> in the puncture gauge approach <cit.>. We also compare the theory with certain subclasses of sGB gravity known to also lead to linearly stable black holes <cit.> and hyperbolic 3+1 evolution (for weak enough scalar fields) <cit.>. 
	
	The paper is organized as follows. In Section <ref> we define the Ricci-coupled scalar-Gauss-Bonnet theory considered in this work and give a brief overview of the modified CCZ4 formalism, together with the concepts of the effective metric and the weak coupling conditions, which are relevant in this manuscript. In Section <ref>, after motivating the specific form of the coupling functions that we are using and introducing the numerical set-up, we present our results regarding the non-linear evolution of rotating and non-rotating black holes and their comparison within two different set of coupling functions. Finally, in Appendix <ref> we include the equations of motion of the theory in the modified CCZ4 formalism and in Appendix <ref> we test the validity of the developed code.

 We follow the conventions in Wald's book <cit.>. 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.
	
	

§ THEORETICAL BACKGROUND

	
	We consider a scalar-Gauss-Bonnet theory with a Ricci coupling (which belongs to the Horndeski class) corresponding to the following action,
	
    S=1/16π∫ d^4x√(-g)(R+X-β(φ)R+λ(φ)/4ℛ^2_GB) ,

	where R is the Ricci scalar with respect to the spacetime metric g^μν, ℛ^2_GB is the Gauss-Bonnet invariant defined as ℛ^2_GB=R^2-4R_μνR^μν+R_μνρσR^μνρσ, φ is the scalar field with X=-1/2∇_μφ∇^μφ being its kinetic term. The Gauss-Bonnet and Ricci couplings are controlled by arbitrary functions of the scalar field with λ(φ) having dimensions of [length]^2 and β(ϕ) being dimensionless. Its equations of motion yield

     
    (1-β(φ))(R_μν- 1/2R g_μν) + Γ_μν= 1/2∇_μφ∇_νφ -  1/4g_μν∇_αφ∇^αφ +(g_μν- ∇_μ∇_ν)β(φ) ,
          ∇_α∇^αφ=  -  λ'(φ)/4 R^2_GB+β'(φ)R ,

	
	where   Γ_μν is defined as 	
	
    Γ_μν   =    -1/2RΩ_μν - Ω_α^ α(R_μν - 1/2R g_μν) 
    		+ 2 R_α(μΩ^ α_ν) - g_μν R^αβΩ_αβ 
    		+    R^β_ μανΩ^ α_β ,
  
	with Ω_μν= ∇_μ∇_ν λ(φ).
	
	

 §.§ Modified CCZ4 formalism

 The equations of motion that follow from varying (<ref>) in the modified harmonic gauge introduced by
<cit.> and supplemented by constraint damping terms are given by (<ref>) with the following replacement,
    
    R^μν-12R g^μν→ R^μν-12R g^μν+2(δ_α^(μĝ^ν)β-12δ_α^βĝ^μν)∇_βZ^α-κ_1[2n^(μZ^ν)+κ_2n^αZ_αg^μν] ,

    where ĝ^μν and g̃^μν are two auxiliary Lorentzian metrics that ensure that gauge modes and gauge condition violating modes propagate at distinct speeds from physical modes, as in <cit.> [Note that g̃^μν is hidden in the definition of the constraints Z^μ (see <cit.> for further details).]. They can be defined as
    
    g̃^μν=g^μν-a(x)n^μn^ν     ĝ^μν=g^μν-b(x)n^μn^ν ,

    where a(x) and b(x) are arbitrary functions such that 0<a(x)<b(x) and n^μ=1α(δ_t^μ-β^iδ_i^μ) is the unit timelike vector normal to the t≡ x^0 =const. hypersurfaces with α and β^i being the lapse function and shift vector of the 3+1 decomposition of the spacetime metric, namely
    
    ds^2=-α^2dt^2+γ_ij(dx^i+β^idt)(dx^j+β^jdt) .

   The damping terms in (<ref>), whose coefficients should satisfy κ_1>0 and κ_2>-22+b(x), guarantee that constraint violating modes are exponentially suppressed <cit.>.

    In Appendix <ref> we have written down the evolution equations for the 3+1 formalism. The versions of the 1+log slicing and Gamma-driver evolution equations that result in the modified puncture gauge are
    
    ∂_tα=β^i∂_iα-2α1+a(x)(K-2Θ) ,
    ∂_tβ^i=β^j∂_jβ^i+34Γ̂^i1+a(x)-a(x)α∂_iα1+a(x) ,

    where Θ=Z^0, K is the trace of the extrinsic curvature of the induced metric γ_ij, Γ̂^i=γ̃^klΓ̃^i_kl+2γ̃^ijZ_j, with Γ̃^i_kl being the Christoffel symbols associated to the conformal spatial metric γ̃_ij≡χγ_ij, where χ=(γ_ij)^-1/3.

	

 §.§ Effective metric


 The hyperbolicity of the equations of motion is held when its principal part is diagonalizable with real eigenvalues and a complete set of linearly independent and bounded eigenvectors that depend smoothly on the variables. The eigenvalues from the gauge sectors lie on the null cones of the auxiliary metrics, while the physical sector is described by a characteristic polynomial of degree 6 which factorises into a product of
quadratic and quartic polynomials <cit.>. The former is defined in terms of an “effective metric” and
is associated with a “purely gravitational” polarization, whereas the latter generically involves
a mixture of gravitational and scalar field polarizations. 

Even though the “fastest” degrees of freedom are associated with the quartic polynomial <cit.>, it is not necessarily the case
that hyperbolicity loss should occur first in their sector. Moreover, there is no simple way to study the hyperbolicity of that sector. This is why we have focused on the “purely gravitational” polarizations, which appear to coincide with the breakdown of the simulation as was seen in <cit.> and in this work.

In the Ricci-coupled sGB theory, the effective metric yields
	
    g_eff^μν=g^μν(1-β(φ))-Ω^μν ,

	and its determinant (normalized to its value in pure GR) can be expressed as
	
    (g^μν_eff)/(g^μν)   =   (11+Ω^⊥⊥-β(φ))^2{1χ[(γ^ij (1-β(φ)) - Ω^ij)(1 + Ω^⊥⊥- β(φ)) -2αΩ^⊥(iβ^j).
          . - (1-β(φ))Ω^⊥⊥β^i β^jα^2+ Ω^⊥ iΩ^⊥ j]} ,

	where Ω^ij=γ^i_μγ^j_νΩ^μν, Ω^⊥ i=-n_μγ^i_νΩ^μν and Ω^⊥⊥=n_μn_νΩ^μν. In the results shown later we will consider the normalized determinant
 
    G_ eff≡ (1+Ω^⊥⊥-β(φ))^2 (g^μν_eff)/(g^μν) ,

 which has no divergences when hyperbolicity is lost and is normalized to unity in the absence of any scalar field.
	
	

 §.§ Weak coupling condition


 One of the main reasons for the relevance of the sGB theory is that it accounts for the more general parity-invariant (up to field redefinitions)
scalar-tensor theory of gravity up to four derivatives, when considering a scalar field with no potential and neglecting the four-derivative scalar term, which
we see from our work in <cit.> is justified since it is always subdominant
to the effect of the Gauss-Bonnet term.

In this sense, we view these theories as effective field theories (EFTs) that arise as a
low energy limit of a more fundamental theory, whose terms are organised in a derivative expansion and appear multiplied by dimensionful coupling
constants that encode the effects of the underlying (unknown) microscopic theory.

Therefore, in order for our theory to be justified and valid as an EFT, one has to make sure that we are not beyond the threshold where the EFT breaks down and the higher derivative terms would become relevant. This is ensured as long as our theory is in the weak coupling regime throughout all its evolution. Namely, we require that the contributions of the Gauss-Bonnet term to the field equations are smaller than the two-derivative Einstein-scalar field terms, which can be expressed in the form of the following weak coupling condition <cit.>:
	
    √(|λ'(φ)|)/L ≪ 1 ,

	where L accounts for any characteristic length
scale of the system associated to the spacetime curvature
and the gradients of the scalar field, which can be computed as

    L^-1=max{|R_ij|^1/2,|∇_μφ|,|∇_μ∇_νφ|^1/2,|ℛ^2_GB|^1/4} .


	
	

§ RESULTS
 
	

 §.§ Coupling functions

	We will concentrate on the following forms of the coupling functions λ(φ) and β(φ):
	
    λ(φ)    =   λ_GB φ^2 + γ_GB φ^4, 
    β(φ)    =   β_Ricc φ^2.

	As it is well known, a pure φ^2 term in the λ(φ) function is enough to admit scalarized black holes <cit.>, but they are linearly unstable. The minimum modification to stabilize the solutions is to add a φ^4 term <cit.> or an alternative is to slightly change the theory, such as the introduction of a Ricci scalar coupling β(φ) <cit.> considered in the present paper. Thus, if γ_GB and β_Ricc are large enough by absolute value, the resulting black hole solutions are linearly stable.
	
	In the first part of the results presented below, we will consider the case of γ_GB=0 because our main goal is to examine the effect of the Ricci scalar coupling on the hyperbolicity. In the second part of the results, we will compare the effects of the γ_GB φ^4 term in λ(φ) on the one hand and the Ricci coupling on the other. The motivation is that these are the simplest modifications of pure sGB gravity with a λ(φ)=λ_GB φ^2 coupling that lead to a restoration of stability and share similar properties of the solutions.
	
	

 §.§ Numerical set-up and hyperbolicity loss treatment

	It was shown in <cit.> that the 1+1 non-linear evolution within the Ricci-coupled sGB theory is hyperbolic for an extensive region of the parameter space. It is natural to generalize these results to a 3+1 evolution where fixing a gauge is much more subtle. We have all reasons to believe that the modified gauge for sGB gravity <cit.> will also work for the considered theory with a Ricci scalar coupling. We also conjecture that, similarly to sGB gravity <cit.>, the loss of hyperbolicity, at least for the considered simulations, is related to the physical modes of the purely gravitational sector rather than the mixed scalar-gravitational one. This is based on the observation that the determinant of the effective metric (<ref>) turns negative right before the breakdown of the simulation. This is a signal that either the speed of these modes diverge or they become degenerate <cit.>.
	
	An important property of the sGB modified gauge proposed in <cit.> is that the mathematical proof for well-posedness is valid only in the weak coupling regime where the scalar field should be weak enough. As a matter of fact, in practice, hyperbolicity is also preserved when the weak coupling condition is slightly violated <cit.>. It is natural to assume that this will also remain true when we consider an additional coupling between the Ricci scalar and the scalar field. This is why we have performed a series of numerical relativity simulations to probe the hyperbolicity of the employed Ricci-coupled sGB theory. A newly developed modification of  <cit.> (based on  <cit.>) was implemented taking into account the Ricci scalar coupling in eqs. (<ref>), which are explicitly written down in our modified CCZ4 formalism in Appendix <ref>. Details about the code convergence and constraint violation are presented in  Appendix <ref>.
	
	We consider as initial data a Kerr black hole with a scalar field Gaussian pulse superimposed on it. In the theory we study, the Kerr black hole is a solution of the field equations but the parameters are chosen in such a way so that it is linearly unstable. As the evolution proceeds and the scalar pulse “hits” the black hole, the scalar hair starts developing  quickly until it reaches equilibrium or a loss of hyperbolicity occurs. The resolution that we mostly worked with is 128 points in each spatial direction with 6 refinement levels and a domain size of 256M.
	
	The employed puncture gauge enables us to evolve the spacetime also through the black hole horizon, thus no explicit excision of the horizon is being made. Still, in order to achieve a stable evolution one has to “turn-off” the Gauss-Bonnet coupling inside the black hole <cit.> and practically evolve GR in the interior. As the apparent horizon is being approached, the Gauss-Bonnet term is gradually turned on so that outside the apparent horizon we are solving the full field equations. As long as the turning on and off of the sGB coupling is performed entirely inside the apparent horizon, this does not affect the spacetime outside the black hole. As a matter of fact, there is a second important reason for switching off the Gauss-Bonnet terms in the black hole interior. Typically, a hyperbolicity loss develops first inside the black hole horizon and then it can emerge above it <cit.>. As long as this non-hyperbolic region is entirely inside the horizon, it is casually disconnected from the rest of the spacetime and it can be still accepted as a viable black hole solution. From the point of view of a numerical relativity code, though, as long as an elliptic region forms anywhere inside the computational domain it leads to unavoidable numerical divergences. Therefore, if we want to determine the threshold between hyperbolic and non-hyperbolic solutions (outside the apparent horizon) it is desirable to cure the black hole interior through the described procedure of switching on and off the Gauss-Bonnet terms.
	
	

 §.§ Hyperbolicity of black hole non-linear evolution in Ricci-coupled sGB theory

	

 First, we start presenting the results for the evolution of sequences of non-rotating black holes with fixed γ_GB=0 and increasing mass. Two different values of the Ricci coupling constant β_Ricc are considered, being adjusted in such a way that the resulting static black hole solutions are linearly stable. The mass and the scalar field at the black hole horizon are plotted in Fig. <ref> for sequences of models at the end state of the numerical relativity simulations of black hole scalarization (after the metric and scalar field stabilize and become nearly static). Red and black squares in the figure correspond to the two values of β_Ricc. Naturally, only the models where the evolution is hyperbolic are depicted because typically hyperbolicity is lost at early times of the scalar field development <cit.>. As evident in Fig. <ref>, we could reach higher maximum scalar fields at the horizon (before hyperbolicity is lost) for the smaller value β_Ricc=2. On the other hand, the range of values of M/√(λ_GB) where the black holes have a well-posed evolution enlarges with the increase of β_Ricc, which is consistent with the findings in spherical symmetry <cit.>.

     As a comparison, with solid and dashed lines in Fig. <ref> we plot the sequence of solutions resulting from solving the set of static field equations similar to <cit.>. Thus, the lines contain all asymptotically flat, regular, and linearly stable black hole solutions regardless of their hyperbolicity. The branches originate from <φ>_ AH=0, at the bifurcation point of the Schwarzschild solution, and they are terminated at some smaller M/√(λ_GB). As one can see in the figure, the lines match very well the points, which is a strong argument for the correctness of the developed extension of . As expected, the lines span a larger range of M/√(λ_GB) as compared to the models resulting from non-linear evolution. The reason for that is that black holes with larger <φ>_ AH  cannot be formed dynamically through a hyperbolic time evolution, i.e. the dynamical variables diverge before the scalar field settles to a constant value. Therefore, similarly to pure sGB gravity <cit.>, only the small scalar field black hole solutions are hyperbolic.
	
	
	

 	
 
	Let us also examine how the evolution of a single black hole looks in the event of a hyperbolicity loss. Several snapshots of the time evolution are depicted in Fig. <ref>, where the top panel represents the scalar field development while in the lower panel one can see the normalized determinant of the effective metric G_ eff defined by eq. (<ref>). The snapshots are adjusted in such a way that the first one is when the scalar field starts developing and is already non-negligible while the last snapshot is the time step just before the code “crashes”. 
	
	
	In the plots of the determinant G_ eff, negative values are depicted with black color. Let us remind the reader that inside the black hole horizon (the dashed white line) we have turned off the Gauss-Bonnet coupling, practically setting λ(φ)=0 and β(φ)=0 in the vicinity of the singularity. For the particular models in Fig. <ref> the cutoff is set at a coordinate radius of roughly r/M≅ 0.9 and, after that, the Gauss-Bonnet term is slowly turned on before the horizon is reached (at r/M≅ 1.07). This ensures that the Gauss-Bonnet term is turned on completely inside the horizon and that this transition region is far enough from the apparent horizon because, otherwise, some undesired numerical error might propagate outside it <cit.>.  Therefore, only the spacetime outside the apparent horizon is a self-consistent solution of the full field equations and G_ eff≅ 1 deep inside the black hole (i.e. the determinant of the effective metric (<ref>) is the same as in GR). 
 
    The most important fact that we observe in the graph is the development of a G_ eff<0 (black) region just before the evolution stops. This is a very strong argument that the breakdown of the code is caused by a hyperbolicity loss in the gravitational sector of physical modes (governed by the effective metric (<ref>)). Therefore, similarly to pure sGB gravity, it is unlikely that this can be improved by a gauge transformation.

    We point out that the loss of hyperbolicity in Fig. <ref> clearly happens inside the apparent horizon. Actually, what typically happens is that a non-hyperbolic region forms inside the black hole horizon, it grows and expands outside it <cit.>, rendering the solution non-hyperbolic[Note that if we have a non-hyperbolic region inside the horizon this does non necessarily mean that the solution is non-physical since this problematic region is casually disconnected from the rest of the spacetime.]. We cannot follow such growth, though, because the code crashes right after the determinant of the effective metric g^μν_ eff turns negative anywhere in the computational domain. Therefore, what is actually observed in simulations, including Fig. <ref>, is that hyperbolicity loss appears right above the region where we turn on the Gauss-Bonnet term even if this is below the apparent horizon. We have checked that, when moving the cutoff radius further inside or outside, hyperbolicity loss still happens for slightly shifted threshold values of the parameters. Nevertheless, the main qualitative features reported here remain unchanged.

	
	

 §.§ Rotating black holes with Ricci scalar coupling

	
 
 A natural question to ask is whether hyperbolicity is preserved for the models depicted in Fig. <ref> in case one includes rotation. 
	For that purpose, we have chosen a model from Fig. <ref> not far away from the point of hyperbolicity loss and performed evolutions for gradually increasing black hole angular momentum. 
	The time evolution of the scalar field on the pole is depicted in Fig. <ref>, which shows that we can perform stable evolution even for very rapidly rotating black holes.
	The maximum depicted value of a_0/M is 0.9. 
	Above that, i.e. close to the extremal limit, we could also perform evolution of the black hole scalarization but in these cases ending up with a stable hairy black hole requires a subtle adjustment of the auxiliary simulation parameters <cit.>.

	
	It is interesting to note that the domain of existence of scalarized rotating black holes in sGB gravity presented so far in the literature <cit.> (excluding the case of spin-induced scalarization <cit.>) seem to be vanishingly small at a moderate a_0/M due to the violation of the regularity condition. Our simulations suggest that, at least for the considered values of the parameters, scalarized black holes with non-negligible scalar field strength might exist up until (or close to) the extremal limit. Of course, we are considering different coupling functions and a Ricci scalar coupling compared to <cit.>. In addition, the time evolution we perform can not be a rigorous proof of the existence of stationary black hole solutions. Still, our results suggest that rotating scalarized black holes, where the scalarization is driven by the spacetime curvature rather than the spin of the black hole, exist up until close to the extremal limit. It is also highly likely that this is not attributed to the Ricci coupling alone, but perhaps a good choice of the coupling function in sGB gravity can lead to the same behavior.
	

	
	

 §.§ Comparison between φ^4 term and Ricci coupling

 	
 
	One of the most important effects of the Ricci coupling term on the spectrum of black hole solutions is that it manages to “stabilize” them and with the increase of β_Ricc the scalar field gets more and more suppressed. But  this is also exactly the effect that a φ^4 term has when added to the λ(φ) coupling function in (<ref>). Of course, the two theories are intrinsically different but it will be interesting to compare the loss of hyperbolicity for both of them. We have already pointed out that the loss of hyperbolicity is mainly controlled by the effective metric g^μν_ eff, which differs slightly in the case with and without a Ricci coupling. It might be interesting to ask how far away from the weak coupling condition one can deviate before the modified gauge <cit.> can no longer secure hyperbolic evolution and how strong the scalar field would be. 
 
 Such a comparison is made in Fig. <ref>, where the time evolution of the scalar field and the weak coupling condition defined by (<ref>) is plotted for models with fixed M/√(λ_GB). The simulations are performed for non-rotating black holes. In the upper panel, β_Ricc=0 and γ_GB is varied (thus we are in sGB gravity with a quadratic and quartic coupling) and in the lower panel γ_GB=0 while β_Ricc varies (Ricci-coupled sGB theory). The ranges of γ_GB and β_Ricc are chosen on the threshold of hyperbolicity loss. A star at the end of some lines marks hyperbolicity loss while for the rest we observe a saturation of the scalar field to a constant. As one can see, in the upper panels the behavior of the weak coupling condition is oscillatory at early times, which is an artifact of the changes in the scalar field gradient before it settles to an equilibrium value.
	
	In both cases, one can go beyond the weak coupling condition while still maintaining hyperbolicity, and the weak coupling condition defined by (<ref>) reaches the order of unity before hyperbolicity is lost. In the case with the Ricci coupling, one is able to reach a larger scalar field before hyperbolicity is lost but the maximum values for the two theories are still of the same order. Of course, this might change from model to model (e.g. when changing M/√(λ_GB)). However, basing ourselves on these results, one can conclude that both theories perform similarly in terms of hyperbolicity loss. 
	
 
	Because of these similarities, both in the evolution and the behavior of the spectrum of solutions, one can ask whether the two theories lead to black holes that can be distinguished through observations. For that purpose, we examined the radial profiles of the metric and the scalar field for two models in sGB gravity with and without Ricci scalar coupling in Fig. <ref>. The parameters of the model are adjusted in such a way that the masses and the scalar charge in the two theories are identical. As one can see, away from the horizon the two solutions look very similar but the differences close to the horizon can be significant. Of course, in this figure we examine only static solutions and it is yet unknown whether the non-linear dynamical will differ as well. Such a study is underway.

 	
 
	

	
    

§ CONCLUSIONS

    In the present paper, we have examined the 3+1 non-linear evolution of static and rotating black holes in scalar-Gauss-Bonnet gravity with an additional coupling between the scalar field and the Ricci scalar. The study was motivated by the recently discovered nice properties of this theory, such as having general relativity as a late-time cosmology attractor and being able to stabilize hairy black hole solutions that are otherwise unstable in certain flavors of pure sGB gravity <cit.>. Extending previous results on hyperbolic spherically symmetric scalar field collapse in sGB gravity with a Ricci coupling <cit.>, we explored in detail the well-posedness of the equations of motion in 3+1 evolutions. For that purpose, a modification of the  code (based on ) was developed in order to handle a self-consistent coupled evolution of the field equations.

    The results show that, as expected from the mathematical analysis, the modified gauge developed in sGB gravity <cit.> also leads to a hyperbolic evolution when adding a Ricci coupling as long as the weak coupling condition is satisfied. As a matter of fact, well-posedness is numerically preserved even slightly above the threshold corresponding to violation of the weak coupling condition. This applies to both static and rotating black holes. As a byproduct of our studies, we have discovered that rotating black holes with a scalar field sourced by the curvature of the spacetime exist for very large angular momenta, close to the extremal limit. This is in contrast with previous studies <cit.> in sGB gravity, where the domain of existence of black holes was getting really narrow as the extremal limit was approached due to a violation of the regularity condition at the horizon. Our results suggest that with a proper choice of the coupling function between the scalar field and the Gauss-Bonnet invariant, similar near-extremal scalarized black holes with a non-negligible scalar field also exist in pure sGB gravity. A systematic study of the stationary solutions in this case is underway.

    Finally, we have compared the results for the threshold of hyperbolicity loss in the Ricci-coupled sGB theory and also in sGB gravity with a more sophisticated coupling function, possessing both quadratic and quartic scalar field terms. Such a coupling has also a stabilization effect on the scalarized solution even for zero Ricci coupling. Our findings confirm that while the two theories are quite different, the threshold for hyperbolicity loss, in terms of scalar field strength and violation of the weak coupling condition, are very similar. In addition, the profiles of the spacetime metric and the scalar field in the black hole solutions are alike in both cases, even though important differences can be present in the near vicinity of the horizon. It is therefore interesting to study to what extent future observations, e.g. of gravitational waves emitted by merging black holes, will be able to distinguish between both.
         
    

§ ACKNOWLEDGEMENTS

    This study is in part financed by the European Union-NextGenerationEU, through the National Recovery and Resilience Plan of the Republic of Bulgaria, project No. BG-RRP-2.004-0008-C01. DD acknowledges financial support via an Emmy Noether Research Group funded by the German Research Foundation (DFG) under grant no. DO 1771/1-1. LAS is supported by an LMS Early Career Fellowship. We thank Nicola Franchini and Farid Thaalba for useful comments on the draft. We also thank Miguel Bezares and Thomas Sotiriou for useful discussions. We acknowledge Discoverer PetaSC and EuroHPC JU for awarding this project access to Discoverer supercomputer resources. We thank the entire  [] collaboration for their support and code development work.
    
	

§ EQUATIONS OF MOTION IN 3+1 FORM

	The 3+1 form of the Einstein equations in our formalism with an arbitrary T^μν yield
	
		
    ∂_⊥γ̃_ij   =    -2αÃ_ij+2γ̃_k(i∂_j)β^k-23γ̃_ij∂_kβ^k, 
    
    			∂_⊥χ   =   23χ(α K - ∂_kβ^k), 
    
    			∂_⊥K   =   -D^iD_iα +α[R+2 D_iZ^i +K(K-2Θ)] -3 κ_1(1+κ_2) α Θ+κ α2[S-3 ρ]
          -3 α b(x)/4(1+b(x))[R-Ã_ijÃ^ij+2/3K^2 -2κ_1(2+κ_2) Θ -2 κ ρ] ,
    
    			∂_⊥Θ   =   α/2[R-Ã_ij Ã^ij+2/3 K^2 
    			+2 D^iZ_i-2 Θ K]-Z_i D^iα-κ_1(2+κ_2) α Θ-κ α ρ
          -b(x)/1+b(x){α2( R-Ã_ij Ã^ij+2/3 K^2 )  -κ_1[2+κ_2] α Θ - κ α ρ} ,
    
    			∂_⊥Ã_ij   =    α[Ã_ij(K-2Θ)-2 Ã_ikÃ^k_ j]+2Ã_k(i∂_j)β^k - 23(∂_kβ^k)Ã_ij
          +χ[α(R_ij + 2D_(iZ_j)-κ S_ij)-D_iD_jα]^TF ,
    
    			∂_⊥Γ̂^i
    			   =   2 α(Γ̃^i_ i klÃ^kl-2/3γ̃^ij∂_jK-3/2 χ Ã^ij∂_jχ) -2 Ã^ij∂_jα-Γ̂^j∂_jβ^i + 2/3 Γ̂^i∂_jβ^j + 13 γ̃^ik∂_k∂_jβ^j + γ̃^jk∂_j∂_kβ^i 
          +2 α γ̃^ij(∂_jΘ- 1/α Θ ∂_jα - 23 K Z_j) -2 κ_1 α γ̃^ijZ_j -2 κ α γ̃^ijJ_j 
          -2α b(x)/1+b(x)[
    			D̃_jÃ^ij-(23)γ̃^ij∂_jK-32 χÃ^ij∂_jχ
    			+γ̃^ij(∂_jΘ-13 K Z_j) 
          - Ã^ijZ_j -κ_1 γ̃^ij Z_j-κ γ̃^ijJ_j
    			] ,

	
	where ∂_⊥=∂_t-β^i∂_i. Taking into account a scalar field with no potential, with an arbitrary β(φ)  coupling to the Ricci scalar and a non-zero contribution of λ(φ)=λ_GBf(φ), with an arbitrary coupling constant λ_GB and function f(φ), the equations become those above with the following decomposition of T_μν, 
	
		
    κ ρ   =   11-β(φ)(12ρ^φ-B+λ_GB ρ^GB) ,
    κ J_i   =   11-β(φ)(12J_i^φ-B_i+λ_GB J_i^GB) ,
    κ S_ij   =   11-β(φ)(12S_ij^φ-B_ij+γ_ij(B-B_nn)+λ_GB S_ij^GB) ,

	
	where the contribution from the kinetic term is given by
	
		
    ρ^φ   =   12(K_φ^2+(∂φ)^2) ,
     J_i^φ   =   K_φ ∂_iφ ,
     S_ij^φ   =   (∂_iφ)(∂_jφ)+12 γ_ij(K_φ^2-(∂φ)^2) ,

	
	with (∂φ)^2=γ^ij(∂_iφ)(∂_jφ) and K_φ=-1α∂_⊥φ. The elements B_nn, B_ij and B_i coming from the decomposition of B_μν=∇_μ∇_νβ(φ)=β'∇_μ∇_νφ+β”∇_μφ∇_νφ yield
	
		
    B_nn=n^μn^νB_μν   =   β”K_φ^2-β'α(D^kα D_kφ+∂_⊥K_φ) ,
    B_i=-γ^μ_ in^νB_μν   =   β”K_φD_iφ+β'(D_iK_φ-K^j_ iD_jφ) ,
    B_ij=γ^μ_ iγ^ν_ jB_μν   =   β”D_iφ D_jφ+β'(D_iD_jφ-K_φK_ij) ,

	
	with B=γ^ijB_ij. With regards to the Gauss-Bonnet sector, we define
	
		
    ρ^GB   =   Ω M2 - M_klΩ^kl , 
    
    			J^GB_i   =   Ω_iM2-M_ijΩ^j - 2(Ω^j_  [iN_j]-Ω^jkD_[iK_j]k) ,

	
	with
	
		
    M_ij   =   R_ij+1χ(29γ̃_ijK^2+13KÃ_ij-Ã_ikÃ_j^ k) , 
    N_i   =   D̃_jÃ_i^ j-32χÃ_i^ j∂_jχ-23∂_iK , 
    Ω_i   =   f'(∂_iK_φ-Ã^j_ i∂_jφ-K3∂_iφ)+f”K_φ∂_iφ ,
    Ω_ij   =   f'(D_iD_jφ-K_φK_ij)+f”(∂_iφ) ∂_jφ ,

	
	where N_i is the GR momentum constraint, and Ω_i and Ω_ij come from the 3+1 decomposition of 𝒞_μν=∇_μ∇_νf(φ), where λ(ϕ)=λ_GB f(φ). In addition, we have 
	
		
    M^TF_ij   ≡ M_ij-13γ_ijM ,
    Ω^TF_ij   ≡Ω_ij-13γ_ijΩ ,

	
	where M=γ^klM_kl is the GR Hamiltonian constraint and Ω=γ^klΩ_kl. Finally, the equations of the two additional degrees of freedom are:
	
		
    ∂_⊥φ   =    -α  K_φ, 
    ∂_⊥K_φ   =   α(-D^iD_iφ + KK_φ) - (D^iφ) D_iα+β'Z^β-λ_GB/4f'(ϕ)ℛ^2_GB .

		
	    
	where
	
    Z^β=R+K_ijK^ij+K^2-2α(∂_⊥K+D_iD^iα) .

	
	All the definitions above enable us to write down the 3+1 equations of γ̃_ij, χ, Θ, Γ̂^i and ϕ with a r.h.s. not depending on the time derivatives of the variables. The rest of the variables (Ã_ij, K and K_φ) include time derivatives in the r.h.s. and that's why we have to specify them with the following linear system,
	
    [       X^kl_ij          Y_ij             0;        X^kl_K           Y_K -3β'2(1-β(φ));      X^kl_K_φ         Y_K_φ             1 ][ ∂_t Ã_kl;     ∂_tK;   ∂_tK_φ ]=
    		[ Z_ij^Ã;    Z^K;  Z^K_φ ],

	where the elements of the matrix are defined as follows,
	
		
    X_ij^kl   =   γ_i^kγ_j^l(1-λ_GB3Ω)+2λ_GB(γ_(i^kΩ_j)^TF,l -γ_ij3Ω^TF,kl-λ_GB f'^2M^TF_ijM^TF,kl) ,
    
    			X_K^kl   =   λ_GB2χ(Ω^TF,kl-λ_GBf'^2M M^TF,kl) , 
    
    			X_K_φ^kl   =   λ_GB2χf'M^TF,kl , 
    
    			Y_ij   =   λ_GB3χ(-Ω^TF_ij-6f'β'M_ij^TF1-β(φ)+λ_GBf'^2M M^TF_ij) , 
    
    			Y_K   =   1-λ_GB3(Ω+3 f'β'M2(1-β(φ))-λ_GB4f'^2M^2 ) ,
    
    			Y_K_φ   =   2β'-λ_GB12f'M ,

	
	while the terms of the r.h.s. are
	
		
    Z_ij^Ã   =   χ[-D_iD_jα+α(R_ij + 2D_(iZ_j) -κ S̅_ij) ]^TF+β^k∂_kÃ_ij+2 Ã_k(i∂_j)β^k-23Ã_ij(∂_kβ^k)
          +α[Ã_ij(K-2Θ)-2Ã_ilÃ^l_ j] ,
    
    			Z^K   =   β^i∂_iK-D^iD_iα +α[R+2 D_iZ^i +K(K-2Θ)] -3 κ_1(1+κ_2) α Θ+κ α2(S̅-3ρ)
          -3 α b(x)/4(1+b(x))[R-Ã_ijÃ^ij+2/3K^2 -2κ_1(2+κ_2) Θ -2 κ ρ] ,
    
    			Z^K_φ   =   β^i∂_i K_φ+α(-D^iD_iφ+K K_φ)-(D^iφ)D_iα+α β' Z̅^β-λ_GB4 α f' ℛ̅^2_GB ,

	
	where the bar denotes that the terms depending on the time derivatives of Ã_ij, K and K_φ of the expressions Z^β, S_ij, S and ℛ^2_GB are subtracted, yielding
	
		
    κ S̅_ij   =   11-β(φ)(12S_ij^φ-B_ij+γ_ij(B-B̅_nn)+λ_GBS̅_ij^GB) ,
    ℛ̅^2_GB   =   -43 M[-1αβ^i∂_iK + 1αD_iD^iα-Ã_ijÃ^ij-K^23]-4  H 
          +8 M^TF,kl[1αD_kD_lα+1χ(Ã_kjÃ^j_ l-Θ̂_kl) ] ,
    Z̅^β   =   R+K_ijK^ij+K^2+2α(β^i∂_iK-D_iD^iα) ,

	
	with S̅=γ^ijS̅_ij, B̅_nn=β”K_φ^2-β'α(D^kα D_kφ-β^k∂_kK_φ) and
	
		
    S̅_ij^GB,TF   =   -13(Ω^TF_ij-λ_GBf'^2MM^TF_ij)[-1αβ^i∂_iK
    			+1αD_iD^iα-Ã_klÃ^kl-K^23]
          -M_ij^TF[Ω+f”(K_φ^2-(∂φ)^2)-β'f'Z^β-λ_GBf'^2H ]
          -13 Ω[1αD_iD_jα+1χ(Ã_imÃ^m_ j-Θ̂_ij) ]^TF-23 Ω_ij^TF(1αD_kD^kα-Ã_klÃ^kl) 
          +2 Ω_(i^TF,k[1αD_j)D_kα+1χ(Ã_j)mÃ_k^ m-Θ̂_j)k) ] +[N_(iΩ_j)]^TF
          -2 (13 γ_ij Ω^TF,kl+λ_GBf'^2M_ij^TFM^TF,kl)[1αD_kD_lα
    			+1χ(Ã_kmÃ^m_ l-Θ̂_kl) ]
          -2(D_kA_ij-D_(iA_j)k)Ω^k-γ_ ij (D^kA_kl) Ω^l+Ω_(iD^kA_j)k ,
    S̅^GB   =   23(Ω-λ_GB4f'^2M^2)[-1αβ^i∂_iK+1αD_iD^iα-Ã_ijÃ^ij-K^23]-2Ω^iN_i-Ω^TF,ij M^TF_ij
          +2  M(14 f”(K_φ^2-(∂φ)^2)-β'4f'Z^β-λ_GB4f'^2H+13 Ω) -ρ^GB
          +(Ω^TF,kl-λ_GBf'^2MM^TF,kl)(1αD_kD_lα+1χÃ_kmÃ_ l^m -Θ̂_klχ) ,

	
	where we have used Θ̂_kl=1αℒ_βÃ_kl+23(K-1α∂_iβ^i )Ã_kl with ℒ_βÃ_ij=β^k∂_k Ã_ij+2Ã_k(i∂_j)β^k and 
	
    H=
    			-43D_iK(N^i+D^iK3) +2 D_iA_jk(D^iA^jk-D^jA^ik)-2 N_i N^i .


	
	

§ CODE TESTING

	In this appendix we present the basic evidence for the validity of the developed extension of . The first test we have made is to verify that the late-time evolution of , namely when the black hole scalarizes and reaches a quasi-equilibrium state, agrees with the results from the solution of the static field equations. This is already presented in Fig. <ref>, where one can see a very good agreement between the masses and the scalar charges obtained by the modified  evolution and the static black hole solutions.

    The average value of the Hamiltonian constraint at the apparent horizon, as well as a convergence plot, are presented in Fig. <ref>. We observe that the convergence matches well to a fourth order, which is consistent with the order of the finite difference stencils, as was also shown in the pure sGB case <cit.>.
    
    
 
	
	
	
