Interior dynamics of neutral and charged black holes in gravity
Abstract
In this paper, we explore the interior dynamics of neutral and charged black holes in gravity. We transform gravity from the Jordan frame into the Einstein frame and simulate scalar collapses in flat, Schwarzschild, and ReissnerNordström geometries. In simulating scalar collapses in Schwarzschild and ReissnerNordström geometries, Kruskal and Kruskallike coordinates are used, respectively, with the presence of and a physical scalar field being taken into account. The dynamics in the vicinities of the central singularity of a Schwarzschild black hole and of the inner horizon of a ReissnerNordström black hole is examined. Approximate analytic solutions for different types of collapses are partially obtained. The scalar degree of freedom , transformed from , plays a similar role as a physical scalar field in general relativity. Regarding the physical scalar field in case, when is negative (positive), the physical scalar field is suppressed (magnified) by , where is the coordinate time. For dark energy gravity, inside black holes, gravity can easily push to . Consequently, the Ricci scalar becomes singular, and the numerical simulation breaks down. This singularity problem can be avoided by adding an term to the original function, in which case an infinite Ricci scalar is pushed to regions where is also infinite. On the other hand, in collapse for this combined model, a black hole, including a central singularity, can be formed. Moreover, under certain initial conditions, and can be pushed to infinity as the central singularity is approached. Therefore, the classical singularity problem, which is present in general relativity, remains in collapse for this combined model.
I Introduction
The internal structure of black holes and spacetime singularities are key topics in gravitation and cosmology Burko_1997_book (); Brady_1999 (); Berger_2002 (); Joshi_2007 (); Henneaux_2008 (), and are great platforms to explore the connection between classical and quantum physics. It is widely believed that rotating black holes exist in reality. According to Price’s theorem, for a collapsing star, the gravitational radiation carries away all the initial features of the star’s gravitational field, except the mass, charge, and angular momentum parameters Price (). As a further step, it is natural to ask what the final state of the internal collapses might be.
Ever since the foundation of general relativity, people have been trying to go beyond it. This endeavor arises from unifying gravitation and quantum mechanics, and addressing some cosmological problems, including the singularity problem in the early Universe and the dark energy problem in the late Universe. Various modified gravity theories have been explored, including scalartensor theory, highdimensional theory, and gravity, etc. For a review of modified gravity theories, see Ref. Clifton_1106 (). For reviews of theory, see Refs. Sotiriou_0805 (); Felice_1002 (); Nojiri_1002 (); Capozziello_1108 ().
Static and spherically symmetric black hole solutions in gravity were explored in Ref. Cruz_0907 (). Charged BornInfeld black holes for theories were studied in Ref. Olmo_1110 (). Instabilities and (anti)evaporation of Schwarzschildde Sitter and ReissnerNordström black holes in modified gravity were discussed in Refs. Nojiri_1301 (); Sebastiani_1305 (); Nojiri_1405 ().
Gravitational collapses in some modified gravity theories have been studied numerically. Spherical collapse of a neutral scalar field in a given spherical, charged black hole in BransDicke theory was investigated in Ref. Avelino_2009 (). Spherical collapses of a charged scalar field in dilaton gravity and gravity were explored in Refs. Borkowska_2011 () and Hwang_1110 (), respectively. Spherical scalar collapse in gravity was simulated in Ref. Guo_1312 (). Asymptotic analysis was implemented in the vicinity of the singularity of a formed black hole.
i.1 Mass inflation
In the vicinity of the central singularity inside a Schwarzschild black hole, the tidal force diverges, and maximal globally hyperbolic region defined by initial data is inextendible. However, inside charged (ReissnerNordström) and rotating (Kerr) black holes, the central singularity is timelike. The globally hyperbolic region is up to the Cauchy horizon, and the spacetime is extendible beyond this horizon to a larger manifold. The ReissnerNordström inner (Cauchy) horizon is a surface of infinite blueshift, which in turn may cause the inner horizon unstable Simpson_1973 (). Furthermore, the strong cosmic censorship conjecture was proposed, which states that for generic asymptotically flat initial data, the maximal Cauchy development is future inextendible. For mathematical explorations of the internal structures of charge black holes, see Refs. Dafermos_2003 (); Dafermos_2014 (). For reviews on the Cauchy problem in general relativity and strong cosmic censorship, see Refs. Ringstrom_2015 (); Isenberg_2015 (), respectively.
The backreaction of the radiative tail from a gravitational collapse on the inner horizon of a ReissnerNordström black hole was investigated by Poisson and Israel Poisson_1989 (); Poisson_1990 (). It was shown that due to the divergence of the tail’s energy density occurring on the inner horizon, the effective internal gravitationalmass parameter becomes unbounded. This phenomenon is usually called mass inflation. These arguments were extended to the rotating black hole case in Ref. Barrabes_1990 ().
In Refs. Poisson_1989 (); Poisson_1990 (), approximate analytic expressions were obtained by considering a simplified model in which the perturbations were modeled by crossflowing radial streams of infalling and outgoing lightlike particles. To get more information, some numerical simulations in more realistic models have been performed. The dynamics of a spherical, charged black hole perturbed nonlinearly by a selfgravitating massless scalar field was numerically studied in Refs. Gnedin_1991 (); Gnedin_1993 (); Brady_1995 (); Burko_1997 (); Burko_1997b (); Hansen_2005 (). Under the influence of the scalar field, the inner horizon of a charged black hole contracts to zero volume, and the center becomes a spacelike singularity. The mass inflation phenomenon was observed. In Refs. Hod_1997 (); Oren_2003 (), with regular initial data, spherical collapse of a charged scalar field was simulated. An apparent horizon was formed. A null, weak massinflation singularity along the Cauchy horizon and a final, spacelike, central singularity were obtained. Spherical collapses in BransDicke theory, dilaton gravity, and gravity were investigated in Refs. Avelino_2009 (); Borkowska_2011 (); Hwang_1110 (). Mass inflation phenomena were also reported.
It is important to connect approximate analytic candidate expressions with numerical results. In Refs. Burko_1998 (); Burko_1999 (), the features of the Cauchy horizon singularity in charge scattering were studied. Analytic and numerical results were compared at some steps.
In this paper, we use the following notations:

Neutral collapse: neutral scalar collapse toward a black hole formation.

Neutral scattering: neutral scalar collapse in a (neutral) Schwarzschild geometry.

Charge scattering: neutral scalar collapse in a (charged) ReissnerNordström geometry. In this process, the scalar field is scattered by the inner horizon of a ReissnerNordström black hole.
i.2 New results
In this paper, we explore neutral scalar collapses in flat, Schwarzschild, and ReissnerNordström geometries in gravity, taking the wellknown HuSawicki model as an example Hu_0705 (). We seek approximate analytic solutions. A generalized MisnerSharp energy in gravity in the Jordan frame was defined in Ref. Cai_0910 (). In this paper, for ease of operation, we mainly work in the Einstein frame and compute the MisnerSharp mass function of the general relativity version, instead. Moreover, we will investigate a dark energy singularity problem.
We explore scalar collapses in both general relativity and gravity. For convenience, we transform the dynamical system in gravity from the Jordan frame into the Einstein frame. In the new system, we equivalently work in Einstein gravity to which a scalar degree of freedom and a physical scalar field are coupled. Basically, plays a similar role as what a physical scalar field does in Einstein gravity. While in gravity, the physical scalar field is suppressed (magnified) when is negative (positive), where is the coordinate time. For simplicity, the results in general relativity are presented in a separate paper Guo_1507 (), and we focus on collapse in gravity in this paper.
According to the strength of the scalar field, charge scattering can be classified into five types as follows:

Type I: spacelike scattering. When the scalar field is very strong, the inner horizon can contract to zero volume rapidly, and the central singularity becomes spacelike. The dynamics near the spacelike singularity is similar to that in neutral collapse.

Type II: null scattering. When the scalar field is intermediate, the inner horizon can contract to a place close to the center or reach the center. For each variable (the metric elements and physical scalar field), the spatial and temporal derivatives are almost equal. In the case of the center being reached, the central singularity is null. This type has two stages: early/slow and late/fast. In the early stage, the inner horizon contracts slowly, and the scalar field also varies slowly. In the late stage, the inner horizon contracts quickly, and the dynamics is similar to that in the spacelike scattering case.

Type III: critical scattering. This case is on the edge between the above two cases. The central singularity becomes null.

Type IV: weak scattering. When the scalar field is very weak, the inner horizon contracts but not much. Then the central singularity remains timelike.

Type V: tiny scattering. When the scalar field is very tiny, the influence of the scalar field on the internal geometry is negligible.
In this paper, we will explore the dynamics of Types I, II, and IV, and obtain approximate analytic solutions for the first two.
By comparing the dynamics in a ReissnerNordström geometry and charge scattering, we investigate the causes of mass inflation and seek further approximate analytic solutions with the following improvements. Usually, doublenull coordinates are used in studies of mass inflation in spherical symmetry. In the line element of doublenull coordinates, the two null coordinates and are present in the form of product . In the equations of motion, mixed derivatives of and are present quite often. In this paper, we use a slightly modified line element, in which one coordinate is timelike and the rest are spacelike. In this case, in the equations of motion, spatial and temporal derivatives are usually separated. This simplifies the numerical formalism and helps to obtain approximate analytic solutions. In addition, we compare numerical results and approximate analytic solutions closely at each step. We compare the dynamics for Schwarzschild black holes, ReissnerNordström black holes, neutral collapse, and charge scattering. We treat the system as a mathematical dynamical system rather than a physical one, examining the contributions from all the terms in the equations of motion.
In Ref. Hwang_1110 () where spherical charged scalar collapse in gravity was simulated, a singularity problem was reported. When a dark energy model is used, can be pushed to easily. Correspondingly, the Ricci scalar becomes singular. This singularity problem can be avoided when an model is used instead. However, the causes of this singularity problem were not explained. In this paper, we will consider a simpler case. Instead of simulating the collapse of a charged scalar field, we study neutral scalar collapse in a ReissnerNordström geometry. The same singularity problem is found. By analyzing the contributions from all the terms in the equations of motion for with , we interpret the causes for the singularity problem. Basically, near the inner horizon, in the equation of motion for , the scalar field and the geometry construct a positive feedback system. Depending on initial conditions, can be accelerated either in positive or in negative directions, until singularities are met. In the negative case, can be accelerated to negative infinity. Correspondingly, goes to zero as the central singularity is approached. However, in the positive case, can be pushed to zero in a short time. Correspondingly, and the Ricci scalar are pushed to and infinity, respectively. This is the cause of the singularity problem. Taking into account quantumgravitational effects at high curvature scale, one may obtain an additional term to the Lagrangian for gravity Starobinsky_1980 (); Vilenkin_1985 (). When this term is added to the function, a singular is pushed to regions where is also singular. Therefore, the singularity problem can be avoided Starobinsky_1980 (); Vilenkin_1985 (); Nojiri_0804 (); Bamba_0807 (); Capozziello_0903 (); Appleby_0909 (); Bamba_1012 (); Bamba_1101 ().
Although the dark energy singularity problem is avoided in the combined model (a combination of a dark energy model and the model), the classical singularity problem, which is present in general relativity, remains in collapse for this model. Under certain initial conditions, near the central singularity, can be positive. Then the positive feedback system in the equation of motion for can push and to positive infinity.
This paper is organized as follows. In Sec. II, we build the framework for charge scattering, including action for charge scattering, the coordinate system, and the model. In Sec. III, we set up the numerical formalism for charge scattering. In Secs. IV, V, and VI, scalar collapses in flat, Schwarzschild, and ReissnerNordström geometries will be explored, respectively. In Sec. VII, we consider weak charge scattering. In Sec. VIII, we discuss the causes and avoidance of the singularity problem. In Sec. IX, the results will be summarized.
In this paper, we set .
Ii Framework
In this section, we build the framework for charge scattering in gravity, in which a selfgravitating massless scalar field collapses in a ReissnerNordström geometry in gravity. Compared to general relativity, in this process, there is one extra scalar degree of freedom . For convenience, gravity is transformed from the Jordan frame into the Einstein frame. For comparison and verification considerations, we use Kruskallike coordinates, and set up the initial conditions by modifying those in a ReissnerNordström geometry with a physical scalar field, a scalar degree of freedom , and the potential for . The HuSawicki model is used as an example.
ii.1 Action
The action for charge scattering in gravity can be written as follows:
(1) 
with
(2)  
(3) 
, , and are the Lagrange densities for gravity, a physical scalar field , and the electric field for a ReissnerNordström black hole, respectively. is a certain function of the Ricci scalar , and is the Newtonian gravitational constant. is the electromagneticfield tensor for the electric field of a ReissnerNordström black hole.
The energymomentum tensor for the massless scalar field is
(4) 
The electric field of a ReissnerNordström black hole can be treated as a static electric field of a point charge of strength sitting at the origin . In the ReissnerNordström metric, the only nonvanishing components of are . The corresponding energymomentum tensor for the electric field is Poisson_2004 ()
(5)  
Although Eq. (5) is obtained in the ReissnerNordström metric, it is valid in any coordinate system, since as seen by static observers, the electromagnetic field should be purely electric and radial Poisson_1990 (); Poisson_2004 ().
ii.2 theory
The equivalent of the Einstein equation in gravity reads
(6) 
where and . The trace of Eq. (6) describes the dynamics of ,
(7) 
where is the trace of the stressenergy tensor . Defining a new variable by
(8) 
and a potential by
(9) 
one can rewrite Eq. (7) as
(10) 
The field equations for gravity (6) are somewhat different from the more familiar corresponding ones in general relativity. Therefore, for convenience, we transform gravity from the current frame, which is usually called the Jordan frame, into the Einstein frame, in which the formalism can be formally treated as Einstein gravity coupled to a scalar field.
Rescaling by
(11) 
one obtains the corresponding action of gravity in the Einstein frame Felice_1002 ()
(12)  
where , , , and a tilde denotes that the quantities are in the Einstein frame. The Einstein field equations are
(13) 
where
(14)  
(15) 
is the ordinary energymomentum tensor for the physical matter fields in terms of in the Jordan frame. With the expression for the energymomentum tensor for the scalar field in the Jordan frame, shown in Eq. (4), the corresponding expression in the Einstein frame can be written as
(16)  
which gives
(17) 
In the Jordan frame, in any coordinate system, the energymomentum tensor for the static electric field of a point charge of strength sitting at the origin can be expressed as [see Eq. (5)]
(18) 
Then we have in the Einstein frame,
(19)  
where and are the quantity in the Jordan and Einstein frames, respectively. Since we mainly work in the Einstein frame in this paper, we simply use for . We denote the total energymomentum tensor for the source fields as
(20) 
The equations of motion for and can be derived from the Lagrange equations as
(21) 
(22) 
Alternatively, Eqs. (21) and (22) can be obtained from the corresponding ones in the Jordan frame. Some details are given in the Appendix.
In the Einstein frame, the potential for can be written as
(23) 
Then we have
(24) 
ii.3 Coordinate system
In the studies of mass inflation, the doublenull coordinates described by Eq. (25) are usually used,
(25) 
where and are functions of the coordinates and . and are outgoing and ingoing characteristics (trajectories of photons), respectively. For convenience, in this paper, we use a slightly modified form described by Eq. (26), obtained by defining and Frolov_2004 (),
(26) 
This set of coordinates is illustrated in Fig. 1. Similar to the Schwarzschild metric, the ReissnerNordström metric can be expressed in Kruskallike coordinates Graves_1960 () (also see Refs. Poisson_1989 (); Poisson_1990 (); Poisson_2004 (); Reall_2015 ()). So for ease and intuitiveness, we set the initial conditions close to those of the ReissnerNordström metric in Kruskallike coordinates, taking into account the presence of a physical scalar field , a scalar degree of freedom , and the potential .
In the form of Ref. Reall_2015 (), the ReissnerNordström metric in Kruskallike coordinates in the region of can be written as
(27) 
where and are the locations and surface gravities for the outer and inner horizons of a ReissnerNordström black hole, respectively. is defined implicitly below Reall_2015 (),
(28) 
In this set of coordinates, as implied by Eq. (28), the exact inner horizon is at regions where and are infinite. However, it is found that, even when and take moderate values, still can be very close to the inner horizon, e.g., . (See Fig. 2.) Therefore, at such regions, the interaction between the scalar fields and the inner horizon still can be very strong, then we can investigate mass inflation numerically.
This formalism has several advantages as follows:

In the line element (26), one coordinate is timelike and the rest are spacelike. This is a conventional setup. It is more convenient and more intuitive to use this set of coordinates. For the set of coordinates described by Eq. (25), in the equations of motion, many terms are mixed derivatives of and ; while for the set of coordinates described by Eq. (26), in the equations of motion, spatial and temporal derivatives are usually separated.

We set initial conditions close to those in the ReissnerNordström metric. Consequently, with the terms related to the scalar fields being removed, we can test our code by comparing the numerical results to the analytic ones in the ReissnerNordström case conveniently. Moreover, by comparing dynamics for scalar collapse to that in the ReissnerNordström case, we can obtain intuitions on how the scalar fields affect the geometry.

The interactions between scalar fields and the geometry are local effects. In Refs. Gnedin_1991 (); Gnedin_1993 (), the space between the inner and outer horizons are compactified into finite space. This overcompactification, at least to us, makes it a bit hard to understand the dynamics. In the configuration that we choose, the space is partially compactified, and the picture of charge scattering turns out to be simpler.
ii.4 model
For a viable dark energy model, has to be positive to avoid ghosts Nunez (), and has to be positive to avoid the DolgovKawasaki instability Dolgov (). The model should also be able to generate a cosmological evolution compatible with the observations Amendola (); Guo_1305 () and to pass the Solar System tests Hu_0705 (); Justin1 (); Justin2 (); Chiba (); Tamaki_0808 (); Tsujikawa_0901 (); Guo_1306 (). Equivalently, general relativity should be restored at high curvature scale, and the model mainly deviates from general relativity at low curvature scale comparable to the cosmological constant. In this paper, we take a typical dark energy model, the HuSawicki model, as an example. This model reads Hu_0705 ()
(29) 
where is a positive parameter, and are dimensionless parameters, , and is the average matter density of the current Universe. We consider one of the simplest versions of this model, i.e., ,
(30) 
where is a dimensionless parameter. In this model,
(31) 
(32) 
As implied in Eq. (32), to make sure that the de Sitter curvature, for which , has a positive value, the parameter needs to be greater than . In this paper, we set to and set to or . Then, together with Eqs. (23) and (32), these values imply that the radius of the de Sitter horizon is about . Moreover, in the configuration of the initial conditions described in Secs. III.2 and VI, the radii of the outer apparent horizons of the formed black holes are about and , respectively. [See Figs. 6(f) and 11(f).] The potentials in the Jordan and Einstein frames are plotted in Figs. 3(a) and 3(b), respectively.
Iii Numerical setup for charge scattering
In this section, we set up the numerical formalisms for charge scattering in gravity, including field equations, initial conditions, boundary conditions, discretization scheme, and tests of numerical codes.
iii.1 Field equations
In doublenull coordinates (26), using
one obtains the equation of motion for ,
(33) 
where and other quantities are defined analogously. For simplicity, we define and integrate the equation of motion for , instead Frolov_2004 (). The equation of motion for can be obtained by rewriting Eq. (33) as
(34) 
provides the equation of motion for ,
(35) 
In doublenull coordinates, the equations of motion for (21) and (22) become, respectively,
(36) 
(37) 
where
(38) 
The and components of the Einstein equations yield the constraint equations
(39) 
(40) 
Via the definitions of and , the constraint equations can be expressed in coordinates. Equations and generate the constraint equations for the and components, respectively,
(41) 
(42) 
iii.2 Initial conditions
We set the initial data to be time symmetric:
(43) 
Therefore, in this configuration, the constraint equation (41) is satisfied identically. Note that, in this configuration, the values of and at are the same as those in the ReissnerNordström metric case.
We set the initial value for as
(44) 
In this paper, we give two sets of initial conditions:
(45)  
(46) 
where is de Sitter value, defined by . The initial value for is defined to be the same as the corresponding one in the ReissnerNordström case (27),
(47) 
where is defined by Eq. (28) with . We obtain the initial value for in charge scattering by combining Eqs. (33) and (42),
(48) 
We set at the origin as in the ReissnerNordström metric case. The definition domain for the spatial coordinate is . Then can be obtained by integrating Eq. (48) via the fourthorder RungeKutta method from to , respectively. The initial values of , , , and are shown in Fig. 11.
In this paper, we employ the finite difference method. The leapfrog integration scheme is implemented, which is a threelevel scheme and requires initial data on two different time levels. With the initial data at , we compute the data at with a secondorder Taylor series expansion as done in Ref. Pretorius (). Take the variable as an example,
(49) 
The values of and are set up as discussed above, and the value of can be obtained from the equation of motion for (37).
Up to this point, the initial conditions are fixed, with all the field equations being taken into account. The firstorder time derivatives of , , , and at described by Eq. (43) ensure that the constraint equation (41) is satisfied. The equation for at expressed by (48) implies that the constraint equation (42) is satisfied. Computations of , , , and at via a secondorder Taylor series expansion, as expressed by Eq. (49) for the case of , satisfy all the equations of motion.
Note that the region of is included in the initial conditions. This may not be physical. However, we are mainly interested in the interior dynamics of black holes, and then it is not important where the scalar field originally comes from. This setup makes it convenient for us to compare the results of charge scattering to the known solutions of the ReissnerNordström geometry. Therefore, we use this setup as a toy model.
iii.3 Boundary conditions
The values of , , , and at the boundaries of are obtained via extrapolations. In fact, since we are mainly concerned with the dynamics around , the boundary conditions will not affect the dynamics in this region, as long as is large enough.
iii.4 Discretization scheme
In this paper, we implement the leapfrog integration scheme, which is secondorder accurate and nondissipative. We let the temporal and spatial grid spacings be equal, .
The equations of motion for (36) and (37) are coupled. Newton’s iteration method is employed to solve this problem Pretorius (). With the illustration of Fig. 4, the initial conditions provide the data at the levels of “down” and “here”. We need to obtain the data on the level of “up”. We take the values at the level of “here” to be the initial guess for the level of “up”. Then, taking as an example, we update the values at the level of “up” using the following iteration:
where is the residual of the differential equation for the function , and is the Jacobian defined by
We do the iterations for the coupled equations one by one, and run the iteration loops until the desired accuracies are achieved.
iii.5 Tests of numerical code
To make sure that the numerical results are trustworthy, one needs to test the numerical code. We compare the numerical results obtained by the code with the analytic ones for the dynamics in a ReissnerNordström geometry, and examine the convergence of the constraint equations and dynamical equations in charge scattering.
The dynamics in a ReissnerNordström geometry is a special case for charge scattering, in which the contributions from the scalar fields are set to zero. This special case has analytic solutions expressed by Eqs. (27) and (28). Therefore, we can test our code by comparing the numerical and analytic results in the ReissnerNordström geometry. Set , , and . We plot the evolutions of and along the slice in Fig. 5(a). As shown in Fig. 5(a), numerical and analytic results match well at an early stage; while at a late stage where gravity and electric field become strong, the numerical evolutions have a time delay compared to the analytic solutions.
When the numerical results are obtained, we substitute the numerical results into the discretized equations of motion and the constraint equations, and find that they are well satisfied. See Figs. 13 and 15, for example. Moreover, the convergence of the constraint equations (41) and (42) is examined. We assume one constraint equation is thorder convergent: residual=, where is the grid size. Therefore, the convergence rate of the discretized constraint equations can be obtained from the ratio between residuals with two different step sizes,
(50) 
Our numerical results show that both of the constraint equations are about secondorder convergent. As a representative, we plot the results for the constraint equation (41) in Fig. 5(b) for the slice .
Convergence tests via simulations with different grid sizes are also implemented Sorkin (); Golod (). If the numerical solution converges, the relation between the numerical solution and the real one can be expressed by
where is the numerical solution. Then, for step sizes equal to and , we have
Defining and , one obtains the convergence rate
(51) 
The convergence tests for , , , and are investigated. They are all secondorder convergent. As a representative, the results for are plotted in Fig. 5(b) for the slice . The values of the parameters in charge scattering in this section are described at the beginning of Sec. VI. We use the spatial range of and the grid spacings of .
Iv Neutral scalar collapse
In this section, we consider neutral collapse in flat geometry in gravity and discuss the mass inflation which happens in the vicinity of the central singularity of the formed black hole.
iv.1 Numerical setup
The numerical setup in neutral scalar collapse in gravity is discussed in Ref. Guo_1312 (). The dynamical equations for , , , , and can be obtained by setting the terms related to the electric field in the corresponding equations in Sec. III.1 to zero:
(52) 
(53) 
(54) 
(55) 
(56) 
where .
In the equation of motion for (54), the term can create big errors near the center . To avoid such a problem, we use the constraint equation (39) alternatively Frolov_2004 (). Defining a new variable
(57) 
one can rewrite Eq. (39) as the equation of motion for ,
(58) 
In the numerical integration, once the value of at the advanced level is obtained, the value of at the current level will be computed using Eq. (57).
We set the initial data as
(59) 
The initial values for and are defined as
(60)  
(61) 
with , , , and . The parameters for the HuSawicki model (30) are set as and .
The local MisnerSharp mass is defined as Misner ()
(62) 
(See Ref. Hayward () for details on various properties of the MisnerSharp mass/energy in spherical symmetry.) Then on the initial slice , the equations for , , and are Frolov_2004 (); Guo_1312 ()
(63)  
(64)  
(65) 
Set at the origin . Then the values of , , and on the initial slice can be obtained by integrating Eqs. (63)(65) from the center to the outer boundary via the fourthorder RungeKutta method. The values of , , , and at can be obtained with a secondorder Taylor series expansion, as discussed in Sec. III.2. The value of at can be obtained using Eq. (57).
The range for the spatial coordinate is defined to be . At the inner boundary , is always set to zero. The terms in Eq. (55) and in Eq. (56) need to be regular at . Since is always set to zero at the center, so is . Then we enforce and to satisfy at . The value of at is obtained via extrapolation. We set up the outer boundary conditions at via extrapolation. The temporal and the spatial grid spacings are .
The numerical code is secondorder convergent and is the one developed in Ref. Guo_1312 ().