Numerical solutions for the f ( R ) -Klein-Gordon system

. We construct a numerical relativity code based on the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation for the gravitational quadratic f ( R ) Starobinsky model. By removing the assumption that the determinant of the conformal 3-metric is unity, we first generalize the BSSN formulation for general f ( R ) gravity theories in the metric formalism to accommodate arbitrary coordinates for the first time. We then describe the implementation of this formalism to the paradigmatic Starobinsky model. We apply the implementation to three scenarios: the Schwarzschild black hole solution, flat space with non-trivial gauge dynamics, and a massless Klein-Gordon scalar field. In each case, long-term stability and second-order convergence is demonstrated. The case of the massless Klein-Gordon scalar field is used to exercise the additional terms and variables resulting from the f ( R ) contributions. For this model, we show for the first time that additional damped oscillations arise in the subcritical regime as the system approaches a stable configuration.


Introduction
One of the most widely studied modifications of classical General Relativity (GR) are the scalar tensor f (R) theories, in which the Einstein-Hilbert gravitational action is generalized to be a nonlinear function of the Ricci scalar.While considerable progress has been made in understanding the implications of such theories [1,2,3,4,5], the technical nature of these theories limits analytical understanding to systems with highsymmetry or perturbations of such systems.In recent years, numerical relativity has contributed greatly to our knowledge of the dynamics of space-time in classical GR, and provides some hope of having similar applicability to questions in f (R) gravity.
The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [6,7] is a modification of the Arnowitt-Deser-Misner (ADM) formulation [8] which has become one of the most widely used formulations in the development of numerical relativity codes due to its high level of empirical stability and robustness [9].A corresponding BSSN formulation for f (R) gravity has also been shown to be strongly hyperbolic in certain regimes and may prove useful as a basis for numerical implementation [10,11].
The formulation presented in [10] is written with Cartesian coordinates in mind, which can lead to complications when the equations are applied, for instance, to problems with a radial coordinate [9,12].A main source of difficulty is the assumption that the conformal 3-metric has a unit determinant in Cartesian coordinates, and thus incompatible with the natural 3-metric developed using spherical polar coordinates.Another difficulty that arises is the singular nature of the initial data for the conformal connection function in spherical polar coordinates.Nevertheless, it has been shown in the context of GR that this difficulty may be treated by introducing a redefined field [9,12,13,14,15,16] whose initial data is zero.
Here, we present a modification of the BSSN formulation for f (R) gravity developed in [10], which allows for the determinant of the conformal 3-metric to be a general function rather than unity.Such a modification of the BSSN formulation is referred to as the "generalized-BSSN" (GBSSN) formulation and was constructed for the case of GR in [12,13] and examined further in [9].We have implemented the formulation in a numerical code, to demonstrate its potential for the study of scalar tensor f (R) theories in the strong-field regime.In particular, we focus our attention on the quadratic Starobinsky model [17,18,19] of f (R) gravity.This model has attracted increasing attention in recent years in the context of astrophysical applications, as well as having been widely used to provide a geometric origin to dark matter abundances [20].It has also been proposed to explain cosmic inflation [18,21].The modification to standard GR should be prominent in strong-gravity scenarios such as those studied in the present work.
For a specific numerical implementation, we study the Starobinsky-Klein-Gordon (SKG) system, i.e., the Starobinsky f (R) gravity together with a matter action describing a massless Klein-Gordon (KG) scalar field.As tests of the robustness and accuracy of the code, we first apply it to two scenarios which are well-studied in the context of standard GR: (i) the evolution of the Schwarzschild black hole solution (as has previously been done in [9,12,14,15]) using a fourth-order Runge-Kutta (RK4) scheme, and (ii) a dynamical gauge pulse in a Minkowski space-time, as studied in [9,14].While the standard set of GBSSN equations for GR together with the RK4 scheme do not allow for the stable evolution of a gauge pulse in flat space using spherical polar coordinates, it has been shown in [9] that by regularizing certain evolution variables the system can be evolved stably.Alternatively, it was also shown in [14] that stable evolutions are possible using a partially implicit Runge-Kutta (PIRK) scheme [22] without special regularizations.In the present work, we implement the latter and demonstrate the consistency and performance of our numerical approach through these test cases.
Finally, we apply the GBSSN formulation of f (R) gravity to study the evolution of dynamical models in SKG gravity.We develop techniques for determining constraint-satisfying initial data for the geometry given a KG scalar field configuration, and demonstrate evolutions in these theories using the PIRK scheme.In particular, we note that the techniques used in [9] to obtain initial data satisfying the Hamiltonian constraint are not suitable for the general f (R) case.This is because the Hamiltonian constraint in the quadratic f (R) case does not reduce to a linear differential equation and also contains higher-order spatial derivatives compared to the case of GR.We therefore require a different numerical approach to find suitable initial data.In addition, since the f (R) GBSSN equations include additional evolution variables compared to GR, the numerical implementations used in [9,14,16] are not immediately applicable since we also need to perform the time integration of the additional evolution equations.Here, we develop for the first time a numerical implementation which accomodates these additional f (R) GBSSN evolution variables.The results obtained when performing the evolution of the SKG system are then compared with those of the GR case, i.e., the Einstein-Klein-Gordon (EKG) system, which has been studied before using the GBSSN equations in [9,16].
This communication is organized as follows.Section 2 reviews the field equations for a general f (R) gravity theory within the metric formalism, i.e., when the metric tensor components are the only independent fields and the connection is the Levi-Civita one.In Section 3, we review the 3+1 decomposition of Riemannian space-time following the approach of [23,24,25].In Sections 4 and 5 we modify the BSSN formulation of the f (R) gravity given in [10] to accommodate curvilinear coordinates following [9,12,13] where this was done for GR.In Section 6 we consider the obtained system of evolution and constraint equations for the case of spherical symmetry.We emphasize that, up until this point, the discussion applies to general f (R) theories.When subsequently performing numerical simulations in Section 7 using these f (R) GBSSN equations, we limit our consideration to the quadratic Starobinsky model.We consider a number of test cases, including a Schwarzschild black hole, a dynamical gauge pulse in flat space, and an SKG model whose results are compared with those of the EKG system.Some technical details regarding the PIRK scheme that we have implemented are presented in the Appendices for the interested reader.To this end, we discuss explicitly and for the first time how we have treated the additional GBSSN evolution variables that are present in the quadratic f (R) case within the PIRK scheme in order to have stable and convergent numerical simulations.

The f (R) field equations
We let U 4 denote four-dimensional Riemannian space-time which, following [26,27], is defined as the tuple (M, g, ∇) where M is the space-time manifold, g is the metric tensor and ∇ is the Levi-Civita connection.We consider a gravitational theory derived from U 4 together with the total action where ‡ is the gravitational action and S m is some matter action [1].We note that here we use geometrized units with G = c = 1 as well as the signature (−, +, +, +).The Lagrangian density in the gravitational action (3) consists of a general function f (R) of the Ricci scalar R. Varying the total action with respect to the metric tensor leads to the well-known equations of motion [1] f In the last expression, □ := g µν ∇ µ ∇ ν is the space-time covariant d'Alembertian operator and is the energy-momentum tensor.In addition, we use f R := df /dR and f RR := d 2 f /dR 2 §.
Taking the trace of the f (R) field equations (4) yields where T := T µ µ is the trace of the energy-momentum tensor.

The 3+1 decomposition
In order to re-cast the field equations as a system of evolution equations, we carry out a standard 3+1 decomposition of the four-dimensional space-time manifold M with metric g (see, for instance, [25,30,31]).Let t be a non-constant scalar field on M, ‡ We note that the gravitational action (3) may be written as [28,29] where one defines V (φ) := ξ(φ)φ − f (ξ(φ)) and φ := f R (ξ) where ξ is an auxiliary field.In such a formulation, φ describes the additional scalar degree of freedom that is present in general f (R) gravity theories.In its form above, the gravitational action should not be mistaken with the GR one, for which φ = 1 and V (φ) = 0.It is precisely due to the existence of φ ̸ = 1, i.e., f R ̸ = 1, that the dynamics of the subsequent equations of motion differs from its GR counterparts.In particular, when working in the so-called metric formalism such equations become of higher (greater than two) order although the undesired Ostrogradski instability is avoided.In the end, whether we consider the formulation (2) or (3), there will be one additional degree of freedom which needs to be accounted for in the subsequent 3 + 1 decomposition; resulting in additional evolution equations compared to the case of GR. § While it is common to denote derivatives of the f (R) function with respect to the Ricci scalar using a prime, we will use a prime later on in Section 6 to denote derivatives with respect to the radial coordinate.
and define the one-form dt normal to three-dimensional level surfaces of t, which we denote by Σ t .We would like t to serve as our time coordinate, and thus restrict our consideration to functions for which the dual of dt is timelike, so that Σ t is spacelike.We define and also define the normal vector field, n, to the hypersurface Σ t through The scalar α defined as is called the lapse function and gives the ratio between the proper time and coordinate time along integral curves of t.
Given g on M we can induce a 3-metric, dubbed γ, on each Σ t via We use D to denote the unique symmetric affine connection that is compatible with γ, i.e., it satisfies Dγ = 0 , and refer to this as the three-dimensional Levi-Civita connection.For any tensor field B of type (k, m), the three-dimensional covariant derivative DB is related to the spacetime covariant derivative ∇B through the projection into the hypersurface Σ t , The extrinsic curvature of a slice Σ t is the symmetric 2-tensor It can be shown that this is a purely spatial tensor satisfying n µ K µν = n ν K µν = 0.The intrinsic curvature of a slice Σ t is described by the three-dimensional Riemann tensor, 3 R, which is defined through the linear map together with the condition 3 R µνσρ n ρ = 0 , which ensures that 3 R lies entirely within Σ t .The space-time Riemann tensor R can be related to 3 R through projections onto the spacelike slices.Specifically, one can find that which are the Gauss, Codazzi and Ricci equations, respectively.In deriving the last of these, we make use of the fact that the components of the 4-acceleration satisfy Equations ( 16)-( 18) encompass all of the non-trivial contractions of R onto the spacelike slices Σ t and thus give a complete description of the four-dimensional Riemann tensor in terms of n and quantities that can be calculated within a slice Σ t using γ.
The Gauss, Codazzi and Ricci equations are sufficient to re-cast the fourdimensional Einstein field equations in terms of spatial quantities together with their derivatives along n, as originally carried out in the ADM formulation of [8].Similarly, they can be used to derive a 3+1 decomposition of the f (R) field equations (4) which was carried out in [32] and examined further in [10].
For the matter contribution, we follow the usual decomposition of the energymomentum tensor as carried out in, for instance, [25].Given the definition of the energy-momentum tensor (5), the various projections onto the spacelike slices of the foliation determined by n, yield the energy density, momentum density and stress-energy tensor, respectively.

An ADM formulation for f (R) theories
In this section we briefly review an ADM-like formulation for f (R) theories which was developed in [10,32].In order to construct such a formulation, we make use of the 3 + 1 decompositions of space-time and the energy-momentum tensor discussed in the previous section.A significant complication over standard GR results from the additional degrees of freedom embodied by the f (R) ̸ = R part of the gravitational action (3).The corresponding field equations (4) involve second derivatives of f , i.e., fourth derivatives when expanded in terms of the metric.In fact, the field equations (4) involve a d'Alembertian acting on f R which can be expanded as follows: As is done in [10,32], we introduce an auxiliary variable based on the Lie derivative of the 4-Ricci scalar∥ Substituting this definition into the right-hand side of (22) and substituting the trace (6) of the field equations into the left-hand side leads to the following result which can ∥ We note that, while in [10,32] the auxiliary field is defined as L n R, in the present work we have defined the auxiliary field (23) as be interpreted as an evolution equation for W : Meanwhile, the definition (23) itself can be re-cast as an evolution equation for the 4-Ricci scalar, which is promoted to the status of an independent evolution variable provided f RR ̸ = 0.

Evolution equations
The geometrical variables associated with the 3+1 decomposition are the 3-metric, γ, extrinsic curvature, K, the 4-Ricci scalar, R, and the new variable W defined in (23).Evolution equations for W and R can be derived from ( 24) and ( 25), respectively, while the evolution of the 3-metric is determined by the definition of the extrinsic curvature, (13), as usual.In order to derive an evolution equation for the extrinsic curvature, we first note that the field equations, (4) and their trace (6), provide us with the following expression for the Ricci tensor One can project this expression (26) with reference to the 3+1 slicing to yield Substituting the last projected expression (29) into the Ricci equation ( 18) yields an evolution equation for the extrinsic curvature tensor In the last expression, we have used the convention of denoting spatial components using Latin indices.Finally, we transform the right-hand side of the evolution equations into derivatives along the natural time coordinate by defining the vector field where β is spatial and is referred to as the shift vector.In terms of the lapse, shift and 3-metric, the line-element associated with the full metric tensor is [24] ds In summary, the following set of equations provide a prescription for the time evolution of the geometrical variables:

Constraints
Having stated the evolution equations in the previous subsection, we now turn our attention to deriving the constraint equations.Contracting the Gauss equation ( 16) twice leads to into which we can substitute (27) to give This is the Hamiltonian constraint for the f (R) theory.Meanwhile, taking the trace of the Codazzi equation ( 17) and substituting in (28) yields dubbed as the corresponding momentum constraint.As a consequence of the Bianchi identities, if these constraints (38) and ( 39) are satisfied on an initial spacelike slice then they are preserved by the evolution equations.

GBSSN evolution equations
In the case of standard GR, a direct implementation of the evolution system contained in equations ( 33) and ( 34) is known to lead to unstable numerical behaviour in 3 + 1 simulations.A prescription which has been successful in curing this problem is to introduce well-chosen auxiliary dynamical variables.Such a prescription introduces additional degrees of freedom but has the effect of re-casting the partial differential equation system in a strongly hyperbolic form.In particular, the so-called BSSN system [6,7] has proved to be a popular choice in 3+1 numerical relativity.In what follows, we briefly outline the BSSN formulation for f (R) theories following [10,24] and extend this to accommodate arbitrary coordinates.
To begin, we introduced the conformal 3-metric where ϕ is the conformal factor and is related to the determinants of the 3-metric and conformal 3-metric by The connection coefficients of the corresponding conformal three-dimensional Levi-Civita connection D, uniquely satisfying D γ = 0, are related to those of D by where γ ij denotes the inverse of the 3-metric tensor.
The conformal 3-metric γij and factor ϕ are the metric evolution variables of the new system.In standard BSSN, the additional degree of freedom resulting from the introduction of ϕ is fixed by placing a restriction on the determinant of the conformal metric, namely γ = 1.More generally [12,13], we can impose the evolution equation where v ∈ {0, 1}.Setting the parameter v = 0 corresponds to the so-called Eulerian condition while v = 1 corresponds to the Lagrangian description.It is at this point that the formulation presented here differs from that in [10].More specifically, in [10], Cartesian coordinates are assumed and the determinant of the conformal 3-metric is set to unity.Here, we do not impose such a condition on γ and instead assign the evolution equation ( 43); considering arbitrary coordinates as a result.While this was done for the case of GR in [12,13], and dubbed the GBSSN formulation, this has not be done before for a general f (R) theory.In what follows we derive, for the first time, the GBSSN formulation for f (R) gravity.
By making use of equations ( 41) and ( 43), one finds the following evolution equation for the conformal factor [12] where we have made use of the fact that which follows from equation ( 42).An alternative to ϕ which we adopt here, and is similar to what is used in [12,33], is to use the variable χ = e 2ϕ which evolves according to Following the standard BSSN prescription, we decompose the extrinsic curvature tensor K into its trace, K, and trace-free part and carry out a conformal rescaling consistent with that of the 3-metric to define An evolution equation for the conformal metric γ follows by substituting its definition (40) into the evolution equation for the 3-metric (33) and making note of (44).After some algebra, we arrive at Similarly, by substituting the definitions of the new variables into the evolution equation (34) of the extrinsic curvature we find the following evolution equations for K and the components of Ā: where the superscript TF is used to indicate the trace-free part of the corresponding tensor field.Finally, the GBSSN system introduces a further set of auxiliary variables given by components of the conformal Levi-Civita connection, namely Although they can be derived from the conformal 3-metric, these are treated as independent variables.An evolution equation follows from taking the time derivative of equation ( 52) and substituting in equations ( 39), ( 43) and (49): As noted in [9], when one imposes the Lagrangian condition in equation ( 43), i.e., setting v = 1, the GBSSN formulation of GR reduces to the BSSN formulation when specifying Cartesian coordinates.Similarly, here we note that when using Cartesian coordinates, the GBSSN formulation presented above for f (R) gravity reduces to the f (R) BSSN formulation given in [10] when we set v = 1.

Reduction to spherical symmetry
While the equations derived in the previous section apply to general geometries, it is also useful to consider the reduction to spherical symmetry both as a testing ground with fewer variables and simpler equations, as well as a platform for studying certain physical models such as irrotational collapse.For this purpose, we re-write the GBSSN system in terms of standard spherical polar coordinates (r, θ, ψ).In spherical symmetry, the conformal 3-metric (40) has only two independent components to be denoted as a = a(t, r) and b = b(t, r), and takes the diagonal form which was used in [9].The rescaled traceless extrinsic curvature Ā has a single independent component, which we denote A a = A a (t, r), and can be expressed as As suggested in [12], this choice of variable ensures that Āij is trace-free by construction.
For the conformal connection function, a straightforward coordinatization derived from the metric is However, this choice is awkward since its flat space initial data, 3 Γr = −2/r, is singular at the origin.Rather, following [9,13], we introduce the alternative variable ∆i : where 3 • Γ i jk are the connection coefficients associated with the Minkowski metric.These have the advantages that they are components of a true vector and have initial data ∆r = 0 in flat space.In terms of the conformal 3-metric components a and b, the regularized variable ∆r is [9] where the prime denotes differentiation with respect to the radial coordinate r.

Evolution equations
By imposing spherical symmetry and the Lagrangian condition (v = 1), equations ( 46) and (49) provide us with the following evolution equations for the metric variables: The two variables determining the extrinsic curvature, K and A a , evolve according to where the (r, r) component of the trace-free part of the Ricci tensor can be rewritten as In computing the right-hand side of equation ( 63), it is useful to note that for any smooth function E, the following expansion holds The use of the derivative of ∆r in the right-hand side of equation ( 64) has been shown to be important for the numerical stability of the evolution system [9,12,14,16].Additionally, as suggested in [9], we also introduce the field ∆r itself in the second line of equation ( 64) in the form of an additional term which evaluates to zero in the case that the constraint ( 58) is satisfied exactly.The evolution equation for ∆r follows from equations ( 53), ( 58) and (61): As noted above, in the case of f (R) ̸ = R theories we employ the two additional variables R and W .In spherical symmetry their evolution equations (35) and (36) reduce to assuming that f RR ̸ = 0 holds.In fact, in the upcoming sections we shall be considering classes of f (R) models such that f RR > 0 in order to avoid the so-called Dolgov-Kawasaki instability [34].
Finally, we note that the Hamiltonian constraint (38) in terms of the spherically symmetric variables is while the momentum constraint ( 39) is

Gauge conditions
In the present work, we make us of the Bona-Massó family of slicing conditions, i.e., we evolve the lapse function according to where h(α) is some function of the lapse [35,36].More specifically, we will use either the 1+log or harmonic slicing conditions, corresponding to h(α) = 2/α and h(α) = 1 respectively.For the shift vector, we introduce an auxiliary field B r and use the pair of equations In equation (73) µ is a constant parameter which we set either to µ = 0 or µ = 3/4 with the latter corresponding to the so-called Gamma-driver shift condition [37].

Klein-Gordon field
For the specifc f (R) models considered here, the non-GR related contributions vanish in vacuum space-times, and as such we need to impose a matter field to study any deviations from GR.As a simple model, we implement a massless KG field Φ = Φ(x) described by the action whose well-known associated energy-momentum tensor and equations of motion are respectively.In terms of the spatial metric and unit normal vector field, the equation of motion (76) is While we we will study the evolution of the GBSSN variables for the SKG system, the EKG system of equations has been studied numerically in spherical symmetry in [38,39] using the ADM formalism, and in [9,16,38] using the GBSSN formalism.The EKG system has also been studied for the case of axial symmetry in [40].
Following the development in [9] for GR, we define the variables and An evolution equation for Φ follows directly from the definition of Π, namely A spatial derivative of this expression leads to an evolution equation for Ψ Then a Lie derivative of (79) determines the evolution equation of Π, where a, as defined in equation ( 19), is the acceleration vector.Substituting the equation of motion (77) into (82) and relating the acceleration vector to the lapse function through a µ = D µ ln α gives In terms of the metric components, this expression reads ¶ The matter variables enter into the f (R) GBSSN evolution equations through the energy and momentum densities, ¶ We note that one could replace the first and second-order derivatives of Φ with Ψ and its first order derivative respectively in the square bracket of equation (84).However, for our code using the PIRK scheme, we have found that using the first and second derivatives of Φ leads to improved numerical behaviour.
respectively as well as the stress-energy tensor with components We also note that the (r, r) component of the trace-free part of the stress-energy tensor, used in (63), is given by

Numerical results
In this article, we perform numerical evolutions of three different scenarios: (i) the evolution of the Schwarzschild black hole solution; (ii) the evolution of flat (Minkowski) space-time endowed with non-trivial gauge dynamics; and finally (iii) the evolution of a massless KG scalar field.In the following, we shall consider these three scenarios in the context of the so-called Starobinsky f (R) gravity models for which the f (R) function in the gravitational action ( 3) is of the form [19] where ℓ is a constant of dimension length squared + .For the cases (i) and (ii) above, there is no formal difference between the Starobinsky and pure GR scenarios since the evolution variables R and W are identically zero.These two tests, then, demonstrate the robustness of our implementation in the case of standard GR and can be compared with similar results from the literature [9,12,14,15].For the case of a massless KG scalar field, we anticipate that the non-GR contributions from the Starobinsky model will be relevant since the space-time Ricci scalar has non-zero initial data for such a case.
Regarding this third case, we note that as a result of the additional degrees of freedom from the fourth-order derivatives of the metric, the numerical implementations discussed in [9,14,16] are not immediately applicable.Firstly, and as will be discussed in Section 7.3, the numerical methods employed in [9] to solve the Hamiltonian constraint in the initial hypersurface are not suitable for the general f (R) case.Secondly, it is necessary to introduce the evolution variables R and W into the numerical scheme according to the GBSSN equations presented in Section 6.
+ We note here that, for such a quadratic model, there will be non-trivial corrections to GR in the presence of matter.Referring the reader to the scalar-picture representation (2) for f (R) gravity, the field φ will indeed be dynamical.However, we note that for the case of a vacuum, the Schwarzschild solution is the only static and spherically symmetric solution [41].Also, it can be shown how in four dimensions quadratic curvature invariants such as R µν R µν , R µναβ R µναβ and ε µνσγ R µναβ R αβ σγ in a more general gravitational Lagrangian can be reduced to a term proportional to R 2 .See for instance [42] for further details.

Schwarzschild solution
A standard test case for spherically symmetric evolution codes is the Schwarzschild solution which exercises the bulk of the evolution equations with non-trivial gauge dynamics and presents some challenges in regularizing the origin.As initial data, we impose a conformally flat 3-metric by setting a(0, r) = b(0, r) = 1, zero extrinsic curvature K(0, r) = A a (0, r) = 0 and also set ∆r (0, r) = 0.For the conformal factor χ we set where M is the mass of the black hole and has dimensions of length.Such initial data corresponds to the Schwarzschild solution in isotropic coordinates which possesses a wormhole topology.The evolution of such a configuration has been used as a test-bed in a number of studies [9,12,14].The initial lapse is set to α(0, r) = χ(0, r) and evolved using the 1+log slicing condition, corresponding to equation (71) with h(α) = 2/α.The shift vector is initially set to zero, i.e., β r (0, r) = 0, and evolved using the Gamma-driver condition (72)-( 73) with µ = 3/4.The spatial grid straddles the origin with gridpoints located at r n = (n − 1/2) ∆r with a base resolution of ∆r = 0.0125M and time-step of ∆t = ∆r/2.We use even boundary conditions at the origin for scalar functions and components of type (0, 2) tensor fields and odd boundary conditions for vectors.We evolve variables using the RK4 scheme.Spatial derivatives are calculated using second-order finite differences; using forward differences for the shift advection terms and central differences for all others.At the outer boundary, we freeze the data associated with the last two grid points; similar to what is done in [12].
In the left panel of Figure 1 we plot the natural logarithm of the root-mean-square (RMS) of the Hamiltonian constraint (69) over the entire numerical grid as a function of time, noting that it stabilizes at a value of less than 10 −5 at the base resolution.A spatial profile at t = 10M is shown to the right, and is evaluated at three different resolutions to demonstrate second-order convergence of the error.The solid blue, dotted orange and dashed green curves correspond to spatial resolutions of ∆r = 0.0125M , 0.025M and 0.05M respectively.The profiles corresponding to spatial resolutions of ∆r = 0.0125M and ∆r = 0.025M have been rescaled by factors of 16 and 4 respectively; consistent with second-order convergence.We note that for each spatial resolution used, we have kept the time-step as ∆t = ∆r/2.
Aspects of the apparent horizon evolution are plotted in Figure 2, including the apparent horizon position and mass evolution in the first and second panels respectively.The horizon radius r AH is computed by evaluating the expansion of outgoing null rays and numerically solving for Θ(r AH ) = 0 [43].Once the apparent horizon position has been obtained, we compute the mass using where A AH is the area of the apparent horizon.The deviation of the mass from its initial value is small; remaining within 0.2% of M throughout the simulation and within 0.04% at t = 200M .In the third panel, we show the convergence of the apparent horizon mass.Therein, we have plotted the natural logarithm of the difference, ∆M AH := M AH,2∆r − M AH,∆r where M AH,2∆r is the apparent horizon mass computed using a spatial resolution of 2∆r and M AH,∆r is the mass computed using a spatial resolution of ∆r and interpolated onto a grid with spacing 2∆r.The solid blue curve shows this difference rescaled by a factor of 4 for the case where ∆r = 0.0125M while the dashed orange curve shows this difference when ∆r = 0.025M .In the right panel, we have plotted the spatial profiles of the Hamiltonian constraint at t = 10M .The solid blue, dotted orange and dashed green curves show the Hamiltonian constraint for spatial resolutions of ∆r = 0.0125M (rescaled by a factor of 16), 0.025M (rescaled by a factor of 4) and 0.05M respectively.The rescalings used are consistent with second-order convergence.

Gauge pulse in flat space
Non-trivial gauge dynamics in the evolution of flat (Minkowski) space-time can be induced in the metric variables by using the following initial data for the lapse function [9] α(0, r and M AH,∆r denote the horizon mass computed using a spatial resolution of ∆r and 2∆r respectively.The solid blue curve, which is rescaled by a factor of 4, and the dashed orange curve show this difference when the spatial resolution is ∆r = 0.0125M and ∆r = 0.025M respectively.This bottom panel shows that the errors in the mass estimate converge away at second-order. where d and σ have dimensions of length and are the position and width of the initial profile respectively.For our consideration, we take d/σ = 5.By inserting such initial data into the numerical code described in the previous subsection, the GBSSN variables K, A a and ∆r start to diverge early on and we end up with NaN (not a number) values; indicating the presence of instability.We note that we have experimented with both 1+log slicing and harmonic slicing for the lapse function while using either a zero shift vector or the Gamma-driver condition and, regardless of the combination of these gauge choices, the aforesaid instability is noticed * .The unsuitability of the numerical scheme outlined in the previous subsection to evolve a gauge pulse in flat space using spherical polar coordinates has been noted in [9].To perform the stable evolution of a gauge * While the instability occurs for any of the combination of gauge choices mentioned above, we notice that K diverges slightly faster when using harmonic slicing and ∆r diverges faster when using the Gamma-driver condition.
pulse, a number of approaches are possible [9,14].One such approach would be to perform a regularization by introducing new evolution variables through a redefinition of fields [9,38].Another approach involves the application of some partially implicit numerical scheme as opposed to the RK4 method [14,16].In the present work, we use the latter approach by implementing the PIRK method used in [14].Thus, by following [14], we evolve a, b, α, β r and χ explicitly.Subsequently, we then evolve the extrinsic curvature variables K and A a partially implicitly using updated values of the conformal 3-metric components.We then evolve ∆r using updated values of the conformal 3-metric components and the extrinsic curvature, and finally evolve the auxiliary field B r .For a detailed outline of the PIRK scheme used here we direct the reader to Appendix A. Although this procedure allows for non-vanishing shift conditions, as in the Gamma-driver condition, we set the shift to zero for our gauge pulse evolutions.We also note that we make use of 1 + log slicing for our gauge pulse evolutions.
In Figures 3 and 4 we give the results of the evolutions for the case of flat spacetime endowed with non-trivial gauge dynamics.Figure 3 shows spatial profiles of the extrinsic curvature trace K at times t = 0, 5σ, 10σ and 15σ; produced using a spacestep of ∆r = 0.0125σ and a time-step of ∆t = ∆r/2.In Figure 4, we plot the natural logarithm of the RMS of the Hamiltonian constraint, in the left panel, up to t = 200σ; demonstrating the long-term stability of our code for this case.In the right panel, we show the spatial profiles of the Hamiltonian constraint for three resolutions.The solid blue, dotted orange and dashed green curves are produced using spatial resolutions of ∆r = 0.0125σ, ∆r = 0.025σ and ∆r = 0.05σ respectively while keeping ∆t = ∆r/2 for each case.For the cases where ∆r = 0.0125σ and 0.025σ, we rescale the profiles by 16 and 4 respectively which correspond to second-order convergence.It is evident from Figure 4 that these profiles are almost indistinguishable; indicating the second-order convergence of our numerical code for the case of flat space non-trivial gauge dynamics.

Scalar field evolution
Having verified that our numerical code correctly simulates the evolution of both the Schwarzschild solution as well as non-trivial gauge dynamics in flat space, we now turn our attention to the evolution of a massless KG scalar field in a dynamical space-time as generated by the quadratic f (R) Starobinsky model.While the space-time Ricci scalar had initial data of zero for the two previous cases, R will be non-zero for the case of a massless KG field.As a result, there will be non-trivial effects arising from the non-GR contribution of the f (R) Starobinsky model.Indeed, we set the initial data for the KG scalar field using where p is the dimensionless initial amplitude of the Gaussian profile while d and σ have dimensions of length and give, respectively, the position and width of the profile.The  Results corresponding to the evolution of a gauge pulse (92) in the lapse function in flat space using 1 + log slicing with a vanishing shift vector.Here, we plot the evolution of the trace, K, of the extrinsic curvature tensor.The first, second, third and fourth panels show the spatial profiles of K at t = 0, 5σ, 10σ and 15σ respectively.
To produce these plots, we have used a space-step of ∆r = 0.0125σ and a time-step of ∆t = ∆r/2.
initial data for the variable Ψ is computed analytically by taking the spatial derivative of Φ(0, r).For all of the subsequent numerical considerations, we set the initial position at d/σ = 5.We assume a conformally flat initial geometry, setting a = b = 1 and ∆r = β r = K = A a = W = 0.In the initial slice, the space-time Ricci scalar and the threedimensional Ricci scalar coincide, i.e., we have R(0, r) = 3 R(0, r).In terms of the conformal factor, the initial Ricci scalar yields Initial data for χ, and hence for R, requires that we solve the Hamiltonian constraint with the given initial matter profile (93).We note that the momentum constraint (70) is identically zero in the initial slice when substituting in the aforementioned initial data for the extrinsic curvature, matter field and the new field W . From the f (R) Hamiltonian constraint equation (69) we have Figure 4: Results corresponding to the evolution of a gauge pulse (92) in the lapse function in flat space using 1 + log slicing with a vanishing shift vector.Here, we have plotted the spatial profile of the Hamiltonian constraint at t = 10σ for three different spatial resolutions: ∆r = 0.0125σ, 0.025σ and 0.05σ.The profiles corresponding to spacesteps of ∆r = 0.0125σ and 0.025σ have been rescaled by factors of 16 and 4 respectively; consistent with second order convergence.To obtain these results, we have used a timestep of ∆t = ∆r/2 for each choice of spatial resolution.The fact that the lines in this figure lie on top of each other demonstrates the second-order convergence of our numerical relativity code for the case of the evolution of a gauge pulse in flat space.
in the initial slice.For the Starobinsky gravity model (88) equation (95) gives us where the Ricci scalar is related to the conformal factor through equation (94).In the case of GR (ℓ = 0), and as mentioned in [9], one can perform a redefinition of fields based on 1/ √ χ so that equation ( 96) is a linear differential equation in the aforesaid redefined variable.For the Starobinsky model such that ℓ ̸ = 0, this resulting equation is, however, non-linear.Here we have used a Newton's method approach to iteratively solve (96) for the conformal factor; using second-order finite differencing for the derivatives.Upon obtaining the initial conformal factor, the initial Ricci scalar is determined by (94).We now turn our attention to the evolution of this initial data.We once again take the shift vector to be zero throughout the simulation while evolving the lapse function using the harmonic slicing condition, i.e., we use equation (71) with h(α) = 1.In the case of ℓ = 0, i.e., for GR, we evolve the standard GBSSN variables consisting of the metric components, the trace and trace-free part of the extrinsic curvature tensor, the regularized conformal connection function and the three matter fields.To evolve this system using the PIRK method, we need to introduce the matter variables Φ, Ψ and Π into the scheme.Here, still in the pure GR scenario, we follow the procedure of treating the fields Φ and Ψ explicitly, and therefore evolve these in the same way as the metric components, and evolve the field Π partially implicitly after the extrinsic curvature terms and before the regularized conformal connection function.For a detailed discussion on how the matter fields are evolved using the PIRK scheme, the reader is directed to Appendix B.
In order to evolve the SKG system numerically, we need to evolve the fields R and W within the PIRK scheme.It is now a question of whether we treat these fields explicitly or partially implicitly.We have found that, when treating both R and W explicitly, i.e., updating these along with the metric components, and evolving the SKG system, the variable R quickly diverges followed by the divergence of A a , K, ∆r and then the remaining evolution variables; resulting in unstable evolution.However, when updating R explicitly and updating W partially implicitly along with the extrinsic curvature variables, we obtain the stable evolution of the SKG system.In Appendix C we discuss in detail how we introduce the variables R and W into the PIRK scheme for the non-trivial f (R) scenario.
We now turn our attention to demonstrating the stability and convergence of our numerical code for the case of a massless KG scalar field and a non-zero value for ℓ in (88).As a test case, we specify the Starobinsky gravity model by setting ℓ/σ 2 = 10 −4 and specify the initial data for the matter by setting p = 0.01 in (93).We consider two different spatial resolutions: ∆r = 0.025σ and ∆r = 0.05σ while keeping ∆t = ∆r/2 and track the RMS of the Hamiltonian constraint.After obtaining the initial data for the conformal factor when ∆r = 0.025σ, we subsequently interpolate the obtained profile to find the initial data for the geometry when ∆r = 0.05σ.In doing this, we ensure that the initial data is the same for both spatial resolutions.
In Figure 5, we plot the natural logarithm of the RMS of the Hamiltonian constraint for the two specified spatial resolutions; rescaling the RMS associated with the higher resolution by a factor of 4 which corresponds to second-order convergence.In computing the RMS of the Hamiltonian constraint, we ignore the first two points next to the origin for the ∆r = 0.025σ case while ignoring the first grid point next to the origin for the ∆r = 0.05σ case.Therefore, the results of Figure 5 indicate that the numerical code is second-order convergent for r/σ > 0.05.This further indicates that there is a significant amount of numerical error contained very close to the origin as a result of the 1/r coordinate singularity; similar to what occurs in the vacuum cases of the previous subsections.We also note that we have evolved the system to t = 100σ to demonstrate the stability of our code.
Next, we compare the evolutions for a fixed p and different values for ℓ, i.e., we use the same initial data for the matter and consider different realizations within the class of Starobinsky f (R) models.We expect that the most notable differences in the evolutions will occur for high curvatures in subcritical evolutions, i.e., no apparent horizon is detected.We set p = 0.04 and carry out evolutions for the GR case (ℓ = 0) as well as for a number of ℓ ̸ = 0 cases.It is only for the non-GR (ℓ ̸ = 0) cases that R and W are used as evolution variables since the numerical scheme requires the condition f RR ̸ = 0 to be satisfied and where in fact we have assumed f RR > 0 in order to avoid the so-called Dolgov-Kawasaki instability being present as mentioned above.For the GR (ℓ = 0) case we evolve the standard set of GBSSN evolution variables that excludes R For the initial data, we have set p = 0.01 in (93).We have plotted the RMS of the Hamiltonian constraint using spacesteps of ∆r = 0.025σ and ∆r = 0.05σ; rescaling the former by a factor of 4 which is consistent with second-order convergence.
and W .For these evolutions we use a spatial resolution of ∆r = 0.01σ and a time-step of ∆t = ∆r/2.
In Figure 6 we plot the central value of the lapse function α 0 and the trace of the extrinsic curvature tensor K 0 over time for the cases of ℓ = 0 and ℓ/σ 2 = 10 −3 .The dashed blue curves correspond to the case where ℓ = 0 (GR) and the solid orange curves correspond to the case where ℓ/σ 2 = 10 −3 .The profiles are almost identical for the early portion of the evolution.However, the extrema, which occur at around t ∼ 8σ and t ∼ 9σ, are greater in magnitude for the case where ℓ/σ 2 = 10 −3 .Also for the case of ℓ/σ 2 = 10 −3 , we note the appearance of damped oscillations beginning at around t ∼ 10σ.
Let us now study the deviation of the results for ℓ ̸ = 0 from the observed behaviour in the GR case.That is, we wish to analyse how the central value for the lapse function is affected by the non-GR contributions of a Starobinsky f (R) model.In order to do this, we consider the cases where ℓ/σ 2 = 10 −5 , 10 −4 , 5 × 10 −4 and 10 −3 and plot the central value of the lapse function for each value of ℓ, dubbed α 0 , after subtracting the central value for the lapse in the case of GR, dubbed α GR 0 .These are shown in Figure 7.The top-left panel shows α 0 − α GR 0 from t = 0 until t = 20σ.The remaining three panels show zoomed-in pieces of the plots contained in the aforementioned topleft panel.The dash-dotted blue, dotted orange, dashed green and solid red curves correspond to ℓ/σ 2 = 10 −5 , 10 −4 , 5 × 10 −4 and 10 −3 respectively.It is clear from Figure 7 that, as the value for ℓ is decreased, the profiles approach zero as expected.For the cases where ℓ/σ 2 = 5 × 10 −4 and 10 −3 , damped oscillations starting around t ∼ 10σ To produce these plots, we have used a space-step of ∆r = 0.01σ and a time-step of ∆t = ∆r/2.The dashed blue curves are associated with the case where ℓ = 0 while the solid orange curves are associated with that of ℓ/σ 2 = 10 −3 .In both figures, we see that the extrema are greater in magnitude when ℓ/σ 2 = 10 −3 .In addition, around t ∼ 10σ, we notice the occurrence of damped oscillations for the case where ℓ/σ 2 = 10 −3 .The top-left panel shows the plots from t = 0 until t = 20σ while the remaining panels show zoomed-in pieces of these plots.produce these plots, we have set ∆r = 0.01σ, ∆t = ∆r/2 while taking the initial amplitude of the KG field to be p = 0.04.The profiles converge to zero as ℓ decreases and we note that damped oscillations appear in the cases where ℓ/σ 2 = 5 × 10 −4 and ℓ/σ 2 = 10 −3 ; with the amplitudes associated with the latter being greater in magnitude.
are seen while no such oscillations are clearly seen for the smaller values of ℓ/σ 2 = 10 −5 and ℓ/σ 2 = 10 −4 .It is possible that similar oscillations may appear in the cases of ℓ/σ 2 = 10 −5 and ℓ/σ 2 = 10 −4 for larger subcritical values of p, however, this would require further investigation which is not covered in this work.

Conclusions
We have for the first time generalized the BSSN formulation of the f (R) gravity presented in [10] to the case where the determinant of the conformal 3-metric is not necessarily unity, in particular for application to models in spherical symmetry.This modification, which was suggested in [12,13] for the case of GR, is referred to as the GBSSN formulation.The extension of this formalism to f (R) theories of gravity in the metric formalism has allowed us to construct a numerical relativity code that we have applied to a system of one minimally coupled KG scalar field together with the Starobinsky f (R) gravity.We note that, while the GBSSN formalism discussed here is for a general f (R) theory, the aforesaid numerical relativity code is specific to the quadratic Starobinsky model.The latter leads to the usual Einstein term plus a quadratic contribution in the Ricci scalar which in recent years has attracted considerable attention in both cosmological and astrophysical studies.This system is of particular interest since it can be easily generalized, on the one hand, to other scalartensor higher-order gravity theories, and on the other hand, to massive scalar fields with non-minimal couplings.The formalism and its numerical implementation was verified through a number of standard numerical tests in the case of GR, demonstrating consistent and convergent results comparable to those found in the literature.In the case of a Schwarzschild black hole, previously used as a test case in [9,12,14,15], we determined an apparent horizon mass which remained within 0.2% of its expected value for the duration of the evolution.As a second test, we evolved a gauge pulse in order to induce dynamics in the metric components of a flat space-time [9,14].For this test, we used a partially implicit (PIRK) scheme following [14] and demonstrated second-order convergence of the implementation.
Finally, we performed tests of the full SKG system, which in addition to the standard GBSSN variables also incorporates the Ricci scalar and a new variable, W , defined in equation (23).The PIRK treatment of the GBSSN equations was extended to accommodate these additional variables.As a diagnostic, we examined the central value for the lapse function as an indication of the strength of the local field, finding that a larger gravitational R 2 -contribution lead to larger variations in this value.For a scalar field initially positioned at r/σ = 5, damped oscillations begin at around t ∼ 10σ in a subcritical scenario.The amplitude for these oscillations are greater when increasing the R 2 -contribution in the gravitational action.The appearance of oscillatory behaviours in the evolution of the Starobinsky model, as well as in other power-law f (R) models, has been reported in several contexts, ranging from the exterior solutions for realistic static neutron stars [44] to the evolution of the matter power spectrum [45].In many of them, the effect has been understood as a competing balance between the so-called curvature effective fluid, as generated by non-GR f (R) terms, and the standard matter fluid, the latter usually divided by f R ̸ = 0 (see for instance this decomposition in [46], Sec. 2).When the field equations are adequately rewritten, the addition of those two fluids is covariantly conserved, whereas separately none of them are.
Although the GBSSN formulation of f (R) gravity can be applied to a variety of f (R) models satisfying f RR ̸ = 0, we have only considered the Starobinsky gravity theory here given its importance and relevance in strong-gravity regimes.We have shown that the PIRK scheme with R being treated explicitly while treating W partially implicitly allows for the stable evolution of the SKG system.Nonetheless, this explicit-implicit approach may not be convenient to study the evolution of a KG scalar field in the presence of other f (R) models.Therefore, a possible extension of this work would be to consider f (R) models with a non-negligible contribution in the strong-gravity regime and see if the same PIRK scheme would allow for stable evolution.Complementarily, our code could be directly used to study the phenomenon of critical collapse for Starobinsky gravity using the GBSSN formulation.We also note that we have only considered the f (R) gravity here when discussing the GBSSN formulation.Therefore, one would need to construct modifications to this formulation in order to consider more general higherderivative gravitational theories.Nevertheless, the fact that we had to treat one of the additional fields compared to GR in the SKG system partially-implicitly may indicate that any numerical implementation of the GBSSN equations for other higher-derivative theories using the PIRK scheme may require some partially-implicit treatment of new variables.
With the above definitions in mind, let us turn our attention to implementing the evolution of the matter fields in the PIRK scheme given in Appendix A. Below, we refer to the steps stated therein.
To introduce a massless KG scalar field into the scheme presented in Appendix A, we first compute the following in Step 1 of the scheme Φn = Φ n + ∆tF Φ (α n , β and F K,1 while all terms depending on Π are contained in F K,2 and do not appear in F A,2 .On the other hand, for the treatment of the regularized conformal connection function, all terms depending on the relevant matter fields are contained in F ∆r ,1 .

Appendix C. PIRK with f (R) contributions
In this section, we specify how the f (R) (f RR ̸ = 0) GBSSN variables R and W are introduced into the PIRK scheme presented in Appendix A. As mentioned in Subsection 7.3, we treat R explicitly while treating W partially implicitly.Let us define F W and F R as being the right-hand sides of equations ( 67) and (68) respectively.In addition, we define F W,2 := αKW + β r W ′ , (C.1) Since R is treated explicitly, we evolve it together with the metric components and, when including a massless KG scalar field, with Φ and Ψ.That is, in Step  It now remains to state how we treat the terms containing R and W when updating the extrinsic curvature and the regularized conformal connection function.For K and A a , we place the terms containing R into F K,1 and F A,1 while there is no dependence on W .
For ∆r , all terms depending on R and W are placed in F ∆r ,1 .

Figure 1 :
Figure 1: Hamiltonian constraint (69) for the evolution of Schwarzschild wormhole initial data.The left panel shows the natural logarithm of the RMS of the Hamiltonian constraint when we have used a spatial resolution of ∆r = 0.0125M and a time-step of ∆t = ∆r/2.In the right panel, we have plotted the spatial profiles of the Hamiltonian constraint at t = 10M .The solid blue, dotted orange and dashed green curves show the Hamiltonian constraint for spatial resolutions of ∆r = 0.0125M (rescaled by a factor of 16), 0.025M (rescaled by a factor of 4) and 0.05M respectively.The rescalings used are consistent with second-order convergence.

Figure 2 :
Figure2: For the evolution of Schwarzschild initial data, the horizon is an effective diagnostic.In the top panel we see that the horizon radius grows from its initial value of r = M/2 due to gauge dynamics, but within 25M has settled to an effectively constant position.The horizon mass shown in the centre panel similarly experiences a transient oscillation, though remains within 0.2% of M throughout the evolution and within 0.04% at the final time of t = 200M .In the bottom panel, we plot the difference ∆M AH := M AH,2∆r − M AH,∆r where M AH,2∆r and M AH,∆r denote the horizon mass computed using a spatial resolution of ∆r and 2∆r respectively.The solid blue curve, which is rescaled by a factor of 4, and the dashed orange curve show this difference when the spatial resolution is ∆r = 0.0125M and ∆r = 0.025M respectively.This bottom panel shows that the errors in the mass estimate converge away at second-order.

Figure 3 :
Figure3: Results corresponding to the evolution of a gauge pulse (92) in the lapse function in flat space using 1 + log slicing with a vanishing shift vector.Here, we plot the evolution of the trace, K, of the extrinsic curvature tensor.The first, second, third and fourth panels show the spatial profiles of K at t = 0, 5σ, 10σ and 15σ respectively.To produce these plots, we have used a space-step of ∆r = 0.0125σ and a time-step of ∆t = ∆r/2.

Figure 5 :
Figure 5: Natural logarithm of the RMS of the Hamiltonian constraint for the evolution of a massless KG scalar field in a dynamical space-time generated by the quadratic Starobinsky f (R) model with ℓ/σ 2 = 10 −4 in (88).For the initial data, we have set p = 0.01 in (93).We have plotted the RMS of the Hamiltonian constraint using spacesteps of ∆r = 0.025σ and ∆r = 0.05σ; rescaling the former by a factor of 4 which is consistent with second-order convergence.

Figure 6 :
Figure 6: Central values of the lapse function α 0 and the trace of the extrinsic curvature tensor K 0 for the f (R) parameter values ℓ = 0 and ℓ/σ 2 = 10 −3 and we have set p = 0.04.To produce these plots, we have used a space-step of ∆r = 0.01σ and a time-step of ∆t = ∆r/2.The dashed blue curves are associated with the case where ℓ = 0 while the solid orange curves are associated with that of ℓ/σ 2 = 10 −3 .In both figures, we see that the extrema are greater in magnitude when ℓ/σ 2 = 10 −3 .In addition, around t ∼ 10σ, we notice the occurrence of damped oscillations for the case where ℓ/σ 2 = 10 −3 .

Figure 7 :
Figure7: Plots of the central value for the lapse function for various values of ℓ after subtracting off that of GR (ℓ = 0).The dash-dotted blue, dotted orange, dashed green and solid red curves correspond to ℓ/σ 2 = 10 −5 , 10 −4 , 5 × 10 −4 and 10 −3 respectively.The top-left panel shows the plots from t = 0 until t = 20σ while the remaining panels show zoomed-in pieces of these plots.produce these plots, we have set ∆r = 0.01σ, ∆t = ∆r/2 while taking the initial amplitude of the KG field to be p = 0.04.The profiles converge to zero as ℓ decreases and we note that damped oscillations appear in the cases where ℓ/σ 2 = 5 × 10 −4 and ℓ/σ 2 = 10 −3 ; with the amplitudes associated with the latter being greater in magnitude.