Interior dynamics of neutral and charged black holes in f(R) gravity

# Interior dynamics of neutral and charged black holes in f(R) gravity

Jun-Qi Guo    Pankaj S. Joshi Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India
10th July 2019
###### 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 Reissner-Nordström geometries. In simulating scalar collapses in Schwarzschild and Reissner-Nordström geometries, Kruskal and Kruskal-like 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 Reissner-Nordströ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 scalar-tensor theory, high-dimensional 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 Born-Infeld black holes for theories were studied in Ref. Olmo_1110 (). Instabilities and (anti)-evaporation of Schwarzschild-de Sitter and Reissner-Nordströ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 Brans-Dicke 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 (Reissner-Nordströ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 Reissner-Nordströ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 Reissner-Nordströ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 gravitational-mass 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 cross-flowing 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 self-gravitating 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 mass-inflation singularity along the Cauchy horizon and a final, spacelike, central singularity were obtained. Spherical collapses in Brans-Dicke 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:

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

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

3. Charge scattering: neutral scalar collapse in a (charged) Reissner-Nordström geometry. In this process, the scalar field is scattered by the inner horizon of a Reissner-Nordström black hole.

### i.2 New results

In this paper, we explore neutral scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in gravity, taking the well-known Hu-Sawicki model as an example Hu_0705 (). We seek approximate analytic solutions. A generalized Misner-Sharp 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 Misner-Sharp 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:

1. 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.

2. 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.

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

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

5. 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 Reissner-Nordström geometry and charge scattering, we investigate the causes of mass inflation and seek further approximate analytic solutions with the following improvements. Usually, double-null coordinates are used in studies of mass inflation in spherical symmetry. In the line element of double-null 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, Reissner-Nordströ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 Reissner-Nordströ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 quantum-gravitational 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 Reissner-Nordströ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 self-gravitating massless scalar field collapses in a Reissner-Nordströ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 Kruskal-like coordinates, and set up the initial conditions by modifying those in a Reissner-Nordström geometry with a physical scalar field, a scalar degree of freedom , and the potential for . The Hu-Sawicki model is used as an example.

### ii.1 Action

The action for charge scattering in gravity can be written as follows:

 S=∫d4x√−g[f(R)16πG+Lψ+LF], (1)

with

 Lψ =−12gαβψ,αψ,β, (2) LF =−FμνFμν4. (3)

, , and are the Lagrange densities for gravity, a physical scalar field , and the electric field for a Reissner-Nordström black hole, respectively. is a certain function of the Ricci scalar , and is the Newtonian gravitational constant. is the electromagnetic-field tensor for the electric field of a Reissner-Nordström black hole.

The energy-momentum tensor for the massless scalar field is

 T(ψ)μν≡−2√|g|δ(√|g|Lψ)δgμν=ψ,μψ,ν−12gμνgαβψ,αψ,β. (4)

The electric field of a Reissner-Nordström black hole can be treated as a static electric field of a point charge of strength sitting at the origin . In the Reissner-Nordström metric, the only nonvanishing components of are . The corresponding energy-momentum tensor for the electric field is Poisson_2004 ()

 T(F)μν ≡ −2√|g|δ(√|g|LF)δgνμ (5) = 14π(FμρFνρ−14δμνFαβFαβ) = q28πr4⋅diag(−1,−1,1,1).

Although Eq. (5) is obtained in the Reissner-Nordströ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 f(R) theory

The equivalent of the Einstein equation in gravity reads

 f′Rμν−12fgμν−(∇μ∇ν−gμν□)f′=8πTμν, (6)

where and . The trace of Eq. (6) describes the dynamics of ,

 □f′−2f−f′R3−8π3T=0, (7)

where is the trace of the stress-energy tensor . Defining a new variable by

 χ≡dfdR, (8)

and a potential by

 U′(χ)≡dUdχ=2f−f′R3, (9)

one can rewrite Eq. (7) as

 □χ−U′(χ)−8π3T=0. (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

 κϕ≡√32lnχ, (11)

one obtains the corresponding action of gravity in the Einstein frame Felice_1002 ()

 SE = ∫d4x√−~g[12κ2~R−12~gμν∂μϕ∂νϕ−V(ϕ)] (12) +∫d4xLM(~gμνχ(ϕ),ψ,q),

where , , , and a tilde denotes that the quantities are in the Einstein frame. The Einstein field equations are

 ~Gμν=κ2[~T(ϕ)μν+~T(M)μν], (13)

where

 ~T(ϕ)μν =∂μϕ∂νϕ−~gμν[12~gαβ∂αϕ∂βϕ+V(ϕ)], (14) ~T(M)μν =T(M)μνχ. (15)

is the ordinary energy-momentum tensor for the physical matter fields in terms of in the Jordan frame. With the expression for the energy-momentum tensor for the scalar field in the Jordan frame, shown in Eq. (4), the corresponding expression in the Einstein frame can be written as

 ~T(ψ)μν = 1χ(∂μψ∂νψ−12gμνgαβ∂αψ∂βψ) (16) = 1χ(∂μψ∂νψ−12~gμν~gαβ∂αψ∂βψ),

which gives

 ~T(ψ)=~gμν~T(ψ)μν=−~gμν~∂μψ~∂νψχ=T(ψ)χ2. (17)

In the Jordan frame, in any coordinate system, the energy-momentum tensor for the static electric field of a point charge of strength sitting at the origin can be expressed as [see Eq. (5)]

 T(q)μdddν=q28πr\tiny JF4⋅diag(−1,−1,1,1). (18)

Then we have in the Einstein frame,

 ~T(q)μdddν ≡ ~gμα~T(q)αν (19) = gμαχ⋅T(q)ανχ = q28πχ2r\tiny JF4⋅% diag(−1,−1,1,1) = q28πr\tiny EF4⋅diag(−1,−1,1,1),

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 energy-momentum tensor for the source fields as

 ~T(total)μdddddddν=~T(ϕ)μdddν+~T(ψ)μdddν+~T(q)μdddν. (20)

The equations of motion for and can be derived from the Lagrange equations as

 ~□ϕ−V′(ϕ)−1√6κ~T(ψ)=0, (21)
 ~□ψ−√23\leavevmode\nobreak κ~gμν∂μϕ∂νψ=0. (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

 V(ϕ)=χR−f2κ2χ2. (23)

Then we have

 V′(ϕ)=dVdχ⋅dχdϕ=2f−χR√6κχ2. (24)

### ii.3 Coordinate system

In the studies of mass inflation, the double-null coordinates described by Eq. (25) are usually used,

 ds2=−4e−2σdudv+r2dΩ2, (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 (),

 ds2=e−2σ(−dt2+dx2)+r2dΩ2. (26)

This set of coordinates is illustrated in Fig. 1. Similar to the Schwarzschild metric, the Reissner-Nordström metric can be expressed in Kruskal-like 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 Reissner-Nordström metric in Kruskal-like 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 Reissner-Nordström metric in Kruskal-like coordinates in the region of can be written as

 ds2=r+r−k2+r2e−2k+r(rr−−1)1+k+|k−|(−dt2+dx2)+r2dΩ2, (27)

where and are the locations and surface gravities for the outer and inner horizons of a Reissner-Nordström black hole, respectively. is defined implicitly below Reall_2015 (),

 4uv=t2−x2=e2k+r(1−rr+)(rr−−1)−k+|k−|. (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:

1. 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.

2. We set initial conditions close to those in the Reissner-Nordströ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 Reissner-Nordström case conveniently. Moreover, by comparing dynamics for scalar collapse to that in the Reissner-Nordström case, we can obtain intuitions on how the scalar fields affect the geometry.

3. 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 f(R) model

For a viable dark energy model, has to be positive to avoid ghosts Nunez (), and has to be positive to avoid the Dolgov-Kawasaki 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 Hu-Sawicki model, as an example. This model reads Hu_0705 ()

 f(R)=R−R0D1RnD2Rn+Rn0, (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., ,

 f(R)=R−DR0RR+R0, (30)

where is a dimensionless parameter. In this model,

 f′=1−DR20(R+R0)2, (31)
 V′(ϕ)=R3√6κf′2(R+R0)2[1+(1−D)R0R(2+R0R)]. (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 double-null coordinates (26), using

 ~Gtt+~Gxx=8π[~T(total)tdddddddt+~T(total)xdddddddx],

one obtains the equation of motion for ,

 r(−r,tt+r,xx)−r2,t+r2,x=e−2σ(1−8πr2V−q2r2), (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

 −η,tt+η,xx=2e−2σ(1−8πr2V−q2r2). (34)

provides the equation of motion for ,

 −σ,tt+σ,xx+r,tt−r,xxr+4π(ϕ2,t−ϕ2,x+ψ2,t−ψ2,xχ−2e−2σV)+e−2σq2r4=0. (35)

In double-null coordinates, the equations of motion for (21) and (22) become, respectively,

 −ϕ,tt+ϕ,xx+2r(−r,tϕ,t+r,xϕ,x)=e−2σ[V′(ϕ)+1√6κ~T(ψ)], (36)
 −ψ,tt+ψ,xx+2r(−r,tψ,t+r,xψ,x)=√23κ(−ϕ,tψ,t+ϕ,xψ,x), (37)

where

 ~T(ψ)=~gμν~T(ψ)μν=−~gμν~∂μψ~∂νψχ=1χe2σ(ψ2,t−ψ2,x). (38)

The and components of the Einstein equations yield the constraint equations

 r,uu+2σ,ur,u+4πr(ϕ2,u+ψ2,uχ)=0, (39)
 r,vv+2σ,vr,v+4πr(ϕ2,v+ψ2,vχ)=0. (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,

 r,tx+r,tσ,x+r,xσ,t+4πr(ϕ,tϕ,x+ψ,tψ,xχ)=0, (41)
 r,tt+r,xx+2(r,tσ,t+r,xσ,x)+4πr(ϕ2,t+ϕ2,x+ψ2,t+ψ2,xχ)=0. (42)

### iii.2 Initial conditions

We set the initial data to be time symmetric:

 r,t=σ,t=ϕ,t=ψ,t=0dddatdt=0. (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 Reissner-Nordström metric case.

We set the initial value for as

 ψ(x,t)|t=0=a⋅exp[−(x−x0)2b]. (44)

In this paper, we give two sets of initial conditions:

 set 1: ϕ(x,t)|t=0=ϕ0, (45) set 2: ϕ(x,t)|t=0=ϕ0+ae−(x−x0)2, (46)

where is de Sitter value, defined by . The initial value for is defined to be the same as the corresponding one in the Reissner-Nordström case (27),

 e−2σ∣∣t=0=e−2σ∣∣RNt=0=r+r−k2+r2e−2k+r(rr−−1)1+k+|k−|, (47)

where is defined by Eq. (28) with . We obtain the initial value for in charge scattering by combining Eqs. (33) and (42),

 r,xx=−r,tσ,t−r,xσ,x+r2,t−r2,x2r−2πr(ϕ2,t+ϕ2,x+ψ2,t+ψ2,xχ)+12re−2σ(1−8πr2V−q2r2). (48)

We set at the origin as in the Reissner-Nordström metric case. The definition domain for the spatial coordinate is . Then can be obtained by integrating Eq. (48) via the fourth-order Runge-Kutta 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 three-level scheme and requires initial data on two different time levels. With the initial data at , we compute the data at with a second-order Taylor series expansion as done in Ref. Pretorius (). Take the variable as an example,

 ψ|t=Δt=ψ|t=0+ψ,t|t=0Δt+12ψ,tt|t=0(Δt)2. (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 first-order 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 second-order 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 Reissner-Nordströ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 second-order 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:

 ψnewup=ψup−G(ψup)J(ψup),

where is the residual of the differential equation for the function , and is the Jacobian defined by

 J(ψup)=∂G(ψup)∂ψup% .

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 Reissner-Nordström geometry, and examine the convergence of the constraint equations and dynamical equations in charge scattering.

The dynamics in a Reissner-Nordströ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 Reissner-Nordströ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 th-order 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,

 n=log2⎡⎢ ⎢⎣O(hn)O((h2)n)⎤⎥ ⎥⎦. (50)

Our numerical results show that both of the constraint equations are about second-order 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

 Freal=Fh+O(hn),

where is the numerical solution. Then, for step sizes equal to and , we have

 Freal=Fh2+O[(h2)n],
 Freal=Fh4+O[(h4)n].

Defining and , one obtains the convergence rate

 n=log2(c1c2). (51)

The convergence tests for , , , and are investigated. They are all second-order 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:

 r(−r,tt+r,xx)−r2,t+r2,x=e−2σ(1−8πr2V), (52)
 −η,tt+η,xx=2e−2σ(1−8πr2V), (53)
 −σ,tt+σ,xx+r,tt−r,xxr+4π(ϕ2,t−ϕ2,x+ψ2,t−ψ2,xχ−2e−2σV)=0, (54)
 −ϕ,tt+ϕ,xx+2r(−r,tϕ,t+r,xϕ,x)=e−2σ[V′(ϕ)+1√6κ~T(ψ)], (55)
 −ψ,tt+ψ,xx+2r(−r,tψ,t+r,xψ,x)=√23κ(−ϕ,tψ,t+ϕ,xψ,x), (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

 g≡−2σ−ln(−r,u), (57)

one can rewrite Eq. (39) as the equation of motion for ,

 g,u=4π⋅rr,u⋅(ϕ2,u+ψ2,uχ). (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

 r,tt=r,t=σ,t=ϕ,t=ψ,t=0dddatdt=0. (59)

The initial values for and are defined as

 χ(r)|t=0 =a⋅[1−tanh(r−r1)2]+χ0, (60) ψ(r)|t=0 =b⋅tanh(r−r2)2, (61)

with , , , and . The parameters for the Hu-Sawicki model (30) are set as and .

The local Misner-Sharp mass is defined as Misner ()

 gμνr,μr,ν=e2σ(−r2,t+r2,x)≡1−2mr. (62)

(See Ref. Hayward () for details on various properties of the Misner-Sharp mass/energy in spherical symmetry.) Then on the initial slice , the equations for , , and are Frolov_2004 (); Guo_1312 ()

 r,x =(1−2mr)eg, (63) m,r =4πr2[V+12(1−2mr)(ϕ2,r+ψ2,rχ)], (64) g,r =4πr(ϕ2,r+ψ2,rχ). (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 fourth-order Runge-Kutta method. The values of , , , and at can be obtained with a second-order 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 second-order convergent and is the one developed in Ref. Guo_1312 ().