Interior Dynamics of Neutral and Charged Black Holes in f(R) Gravity

Abstract

In this paper, we explore the interior dynamics of neutral and charged black holes in f(R) gravity. We transform f(R) 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 f′ 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 f′, plays a similar role as a physical scalar field in general relativity. Regarding the physical scalar field in f(R) case, when dΦ/dt is negative (positive), the physical scalar field is suppressed (magnified) by Φ, where t is the coordinate time. For dark energy f(R) gravity, inside black holes, gravity can easily push f′ to 1. Consequently, the Ricci scalar R becomes singular, and the numerical simulation breaks down. This singularity problem can be avoided by adding an R² term to the original f(R) function, in which case an infinite Ricci scalar is pushed to regions where f′ 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, f′ and R 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.

1. Introduction

The internal structure of black holes and spacetime singularities are key topics in gravitation and cosmology [1,2,3,4,5], 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 [6]. 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 f(R) f(R) gravity, etc. For a review of modified gravity theories, see Reference [7]. For reviews of f(R) f(R) theory, see References [8,9,10,11].

Static and spherically symmetric black hole solutions in f(R) f(R) gravity were explored in Reference [12]. Charged Born-Infeld black holes for f(R) f(R) theories were studied in Reference [13]. Instabilities and (anti)-evaporation of Schwarzschild-de Sitter and Reissner-Nordström black holes in modified gravity were discussed in References [14,15,16].

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 Reference [17]. Spherical collapses of a charged scalar field in dilaton gravity and f(R) f(R) gravity were explored in References [18] and [19], respectively. Spherical scalar collapse in f(R) f(R) gravity was simulated in Reference [20]. Asymptotic analysis was implemented in the vicinity of the singularity of a formed black hole.

1.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 [21]. 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 References [22,23]. For reviews on the Cauchy problem in general relativity and strong cosmic censorship, see References [24,25], 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 [26,27]. 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 Reference [28].

In References [26,27], 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 References [29,30,31,32,33,34]. 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 References [35,36], 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 f(R) f(R) gravity were investigated in References [17,18,19]. Mass inflation phenomena were also reported.

It is important to connect approximate analytic candidate expressions with numerical results. In References [37,38], 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.

1.2. New Results

In this paper, we explore neutral scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in f(R) f(R) gravity, taking the well-known Hu-Sawicki model as an example [39]. We seek approximate analytic solutions. A generalized Misner-Sharp energy in f(R) f(R) gravity in the Jordan frame was defined in Reference [40]. 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 f(R) f(R)singularity problem.

We explore scalar collapses in both general relativity and f(R) f(R) gravity. For convenience, we transform the dynamical system in f(R) f(R) 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 ϕ(≡3/2ln^f′/8πG) ϕ(≡3/2ln^f′/8πG) and a physical scalar field ψ are coupled. Basically, ϕ plays a similar role as what a physical scalar field ψ does in Einstein gravity. While in f(R) f(R) gravity, the physical scalar field ψ is suppressed (magnified) when dϕ/dt dϕ/dt is negative (positive), where t is the coordinate time. For simplicity, the results in general relativity are presented in a separate paper [41], and we focus on collapse in f(R) f(R)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 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 u and v are present in the form of product dudv dudv. In the equations of motion, mixed derivatives of u and v 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 Reference [19] where spherical charged scalar collapse in f(R) f(R) gravity was simulated, a singularity problem was reported. When a dark energy f(R) f(R) model is used, ^{f′ f′} can be pushed to 1 easily. Correspondingly, the Ricci scalar R becomes singular. This singularity problem can be avoided when an ^{R2 R2} 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 ϕ(≡3/2ln^f′/8πG) ϕ(≡3/2ln^f′/8πG) with ^f′≡df/dR ^f′≡df/dR , 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, ^{f′ f′} goes to zero as the central singularity is approached. However, in the positive case, ϕ can be pushed to zero in a short time. Correspondingly, ^{f′ f′} and the Ricci scalar R are pushed to 1 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 ^{R2 R2} term to the Lagrangian for gravity [42,43]. When this ^{R2 R2} term is added to the f(R) f(R) function, a singular R is pushed to regions where ^{f′ f′} is also singular. Therefore, the singularity problem can be avoided [42,43,44,45,46,47,48,49].

Although the dark energy f(R) f(R) singularity problem is avoided in the combined model (a combination of a dark energy f(R) f(R) model and the ^{R2 R2} 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, dϕ/dt dϕ/dtcan be positive. Then the positive feedback system in the equation of motion for ϕ can push ϕ and R to positive infinity.

This paper is organized as follows. In Section 2, we build the framework for charge scattering, including action for charge scattering, the coordinate system, and the f(R) f(R) model. In Section 3, we set up the numerical formalism for charge scattering. In Section 4, Section 5 and Section 6, scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries will be explored, respectively. In Section 7, we consider weak charge scattering. In Section 8, we discuss the causes and avoidance of the singularity problem. In Section 9, the results will be summarized.

In this paper, we set G=c=4π_ϵ0=1 G=c=4π_ϵ0=1.

2. Framework

In this section, we build the framework for charge scattering in f(R) f(R) gravity, in which a self-gravitating massless scalar field collapses in a Reissner-Nordström geometry in f(R) f(R) gravity. Compared to general relativity, in this process, there is one extra scalar degree of freedom ^f′≡df/dR ^f′≡df/dR . For convenience, f(R) f(R) 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 ^{f′ f′} , and the potential for ^{f′ f′}. The Hu-Sawicki model is used as an example.

2.1. Action

The action for charge scattering in f(R) f(R)gravity can be written as follows:

S=∫^d4x−gf(R)16πG+_Lψ+_LF

with

_Lψ=−12^gαβ _{ψ,α ψ,β}

_LF=−_Fμν ^Fμν4

f(R)/(16πG) f(R)/(16πG) , _{Lψ Lψ} , and _{LF LF} are the Lagrange densities for f(R) f(R) gravity, a physical scalar field ψ, and the electric field for a Reissner-Nordström black hole, respectively. f(R) f(R) is a certain function of the Ricci scalar R, and G is the Newtonian gravitational constant. _{Fμν Fμν}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μν=_{ψ,μ ψ,ν}−12_gμν ^gαβ _{ψ,α ψ,β}

The electric field of a Reissner-Nordström black hole can be treated as a static electric field of a point charge of strength q sitting at the origin r=0 r=0 . In the Reissner-Nordström metric Equation (120), the only nonvanishing components of _{Fμν Fμν} are _Ftr=−_Ftr=−q/^r2 _Ftr=−_Ftr=−q/^r2 . The corresponding energy-momentum tensor for the electric field is [50]

_^T(F)νμ≡−2|g|δ(|g|_LF)δ_gμν=14π^Fμρ _Fνρ−14_δνμ ^Fαβ _Fαβ=^q28π^r4·diag(−1,−1,1,1)

Although Equation (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 [27,50].

2.2. f(R) f(R)Theory

The equivalent of the Einstein equation in f(R) f(R)gravity reads

^f′ _Rμν−12f_gμν−_{∇μ ∇ν}−_gμν□^f′=8π_Tμν

where ^f′≡df/dR ^f′≡df/dR and □≡_∇α ^∇α □≡_∇α ^∇α . The trace of Equation (6) describes the dynamics of ^{f′ f′},

□^f′−13(2f−^f′R)−8π3T=0

where T is the trace of the stress-energy tensor _{Tμν Tμν}. Defining a new variable χ by

χ≡dfdR

and a potential U(χ) U(χ)by

^U′(χ)≡dUdχ=13(2f−χR)

one can rewrite Equation (7) as

□χ−^U′(χ)−8π3T=0

The field equations for f(R) f(R) gravity Equation (6) are somewhat different from the more familiar corresponding ones in general relativity. Therefore, for convenience, we transform f(R) f(R)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χ

one obtains the corresponding action of f(R) f(R) gravity in the Einstein frame [9]

_SE=∫^d4x−g˜12^κ2R˜−12^g˜μν _∂μϕ_∂νϕ−V(ϕ)+∫^d4x_LMg˜μνχ(ϕ),ψ,q

where κ=8πG κ=8πG , _g˜μν=χ·_{gμν g˜μν}=χ·_gμν , V(ϕ)≡(χR−f)/(2^{κ2 χ2}) V(ϕ)≡(χR−f)/(2^{κ2 χ2}), and a tilde denotes that the quantities are in the Einstein frame. The Einstein field equations are

_G˜μν=^κ2_T˜μν(ϕ)+_T˜μν(M)

where

_T˜μν(ϕ)=_∂μϕ_∂νϕ−_g˜μν12^g˜αβ _∂αϕ_∂βϕ+V(ϕ)

_T˜μν(M)=_Tμν(M)χ

_{Tμν(M) Tμν(M)} is the ordinary energy-momentum tensor for the physical matter fields in terms of _{gμν gμν} in the Jordan frame. With the expression for the energy-momentum tensor for the scalar field ψ in the Jordan frame, shown in Equation (4), the corresponding expression in the Einstein frame can be written as

_T˜μν(ψ)=1χ_∂μψ_∂νψ−12_gμν ^gαβ _∂αψ_∂βψ=1χ_∂μψ_∂νψ−12_g˜μν ^g˜αβ _∂αψ_∂βψ

which gives

^T˜(ψ)=^g˜μν _T˜μν(ψ)=−^g˜μν _∂˜μψ_∂˜νψχ=^{T(ψ) χ2}

In the Jordan frame, in any coordinate system, the energy-momentum tensor for the static electric field of a point charge of strength q sitting at the origin r=0 r=0 can be expressed as (see Equation (5))

_Tdddν(q)μ=^q28π^_rJF4·diag(−1,−1,1,1)

Then we have in the Einstein frame,

_{T˜dddν(q)μ}≡^g˜μα _T˜αν(q)=^gμαχ·_Tαν(q)χ=^q28π^{χ2 _rJF4}·diag(−1,−1,1,1)=^q28π^_rEF4·diag(−1,−1,1,1)

where _{rJF rJF} and _{rEF rEF} are the quantity r in the Jordan and Einstein frames, respectively. Since we mainly work in the Einstein frame in this paper, we simply use r for _{rEF rEF}. We denote the total energy-momentum tensor for the source fields as

_{T˜dddddddν(total)μ}=_{T˜dddν(ϕ)μ}+_{T˜dddν(ψ)μ}+_{T˜dddν(q)μ}

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

□˜ϕ−^V′(ϕ)−16κ^T˜(ψ)=0

□˜ψ−23κ^g˜μν _∂μϕ_∂νψ=0

Alternatively, Equation (21) and Equation (22) can be obtained from the corresponding ones in the Jordan frame. Some details are given in Appendix A.

In the Einstein frame, the potential for ϕ can be written as

V(ϕ)=χR−f2^{κ2 χ2}

Then we have

^V′(ϕ)=dVdχ·dχdϕ=2f−χR6κ^χ2

2.3. Coordinate System

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

d^s2=4^e−2σdudv+^r2d^Ω2

where σ and r are functions of the coordinates u and v. u and v are outgoing and ingoing characteristics (trajectories of photons), respectively. For convenience, in this paper, we use a slightly modified form described by Equation (26), obtained by defining u=(t−x)/2 u=(t−x)/2 and v=(t+x)/2 v=(t+x)/2 [51],

d^s2=^e−2σ(−d^t2+d^x2)+^r2d^Ω2

This set of coordinates is illustrated in Figure 1. Similar to the Schwarzschild metric, the Reissner-Nordström metric can be expressed in Kruskal-like coordinates [52] (also see References [26,27,50,53]). 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 ^{f′ f′} , and the potential U(^f′) U(^f′).

Figure 1. Initial and boundary conditions for charge scattering. Initial slice is at t=0 t=0 . Definition domain for x is [−xbxb] [−xbxb] . u=(t−x)/2 u=(t−x)/2 and v=(t+x)/2 v=(t+x)/2 .

In the form of Reference [53], the Reissner-Nordström metric in Kruskal-like coordinates in the region of r>_r− r>_r−can be written as

d^s2=_{r+ r−k+2} ^{r2e−2_k+r r_r−−11+_k+|_k−|}(−d^t2+d^x2)+^r2d^Ω2

where _r±(=m±^m2−^q2) _r±(=m±^m2−^q2) and _k±[=(_r±−_r∓]/(2_r±2)] _k±[=(_r±−_r∓]/(2_r±2)] are the locations and surface gravities for the outer and inner horizons of a Reissner-Nordström black hole, respectively. r(t,x) r(t,x) is defined implicitly below [53],

4uv=^t2−^x2=^e2_k+r1−r_r+^{r_r−−1−_k+|_k−|}

Some details for this set of coordinates are given in Appendix B. In this set of coordinates, as implied by Equation (28), the exact inner horizon is at regions where uv uv and (^t2−^x2) (^t2−^x2) are infinite. However, it is found that, even when uv uv and (^t2−^x2) (^t2−^x2) take moderate values, r still can be very close to the inner horizon, e.g., r=(1+¹⁰⁻¹⁰)_r− r=(1+¹⁰⁻¹⁰)_r− . (See Figure 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 Equation (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 Equation (25), in the equations of motion, many terms are mixed derivatives of u and v; while for the set of coordinates described by Equation (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 References [29,30], 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.

Figure 2. Contour lines for r defined by Equation (28) in a Reissner-Nordström geometry with m=1 m=1 and q=0.7 q=0.7 . Although the exact inner horizon is at regions where uv uv and (t2−x2) (t2−x2) are infinite, r can be very close to the inner horizon r=r− r=r− even when uv uv and (t2−x2) (t2−x2) take moderate values.

2.4. f(R) f(R)Model

For a viable dark energy f(R) f(R) model, ^{f′ f′} has to be positive to avoid ghosts [54], and ^{f′′ f′′} has to be positive to avoid the Dolgov-Kawasaki instability [55]. The model should also be able to generate a cosmological evolution compatible with the observations [56,57] and to pass the Solar System tests [39,58,59,60,61,62,63]. Equivalently, general relativity should be restored at high curvature scale, and the f(R) f(R) model mainly deviates from general relativity at low curvature scale comparable to the cosmological constant. In this paper, we take a typical dark energy f(R) f(R) model, the Hu-Sawicki model, as an example. This model reads [39]

f(R)=R−_R0D1 ^Rn_D2 ^Rn+_R0n

where n is a positive parameter, _{D1 D1} and _{D2 D2} are dimensionless parameters, _R0=8π_ρ¯0/3 _R0=8π_ρ¯0/3 , and _{ρ¯0 ρ¯0} is the average matter density of the current Universe. We consider one of the simplest versions of this model, i.e., n=1 n=1,

f(R)=R−D_R0RR+_R0

where D is a dimensionless parameter. In this model,

^f′=1−D_R02^(R+_R0)2

^V′(ϕ)=^R36κ^{f′2 (R+_R0)2}1+(1−D)_R0R2+_R0R

Figure 3. Potentials for f(R) f(R) models. χ≡f′ χ≡f′ and χp≡χ/(1+χ) χp≡χ/(1+χ) . U(χ) U(χ) and V(χ) V(χ) are the potentials in the Jordan and Einstein frames and can be obtained from Equation (9) and Equation (23), respectively. (a) and (b) are for the Hu-Sawicki model Equations (30), f(R)=R−DR0R/(R+R0) f(R)=R−DR0R/(R+R0) , while (c) and (d) are for the combined model Equations (79), f(R)=R−DR0R/(R+R0)+αR2 f(R)=R−DR0R/(R+R0)+αR2 . D=1.2 D=1.2 , R0=10−5 R0=10−5 , and α=1 α=1 .

As implied in Equation (32), to make sure that the de Sitter curvature, for which ^V′(ϕ)=0 ^V′(ϕ)=0 , has a positive value, the parameter D needs to be greater than 1. In this paper, we set D to 1.2 1.2 and set _{R0 R0} to ^{10−5 10−5} or ^{10−6 10−6} . Then, together with Equation (23) and Equation (32), these values imply that the radius of the de Sitter horizon is about 1/_R0∼¹⁰³ 1/_R0∼¹⁰³ . Moreover, in the configuration of the initial conditions described in Section 3.2 and Section 6, the radii of the outer apparent horizons of the formed black holes are about 2.1 2.1 and 3.7 3.7 , respectively. (See Figure 6f and Figure 11f). The potentials in the Jordan and Einstein frames are plotted in Figure 3a,b, respectively.

3. Numerical Setup for Charge Scattering

In this section, we set up the numerical formalisms for charge scattering in f(R) f(R)gravity, including field equations, initial conditions, boundary conditions, discretization scheme, and tests of numerical codes.

3.1. Field Equations

Details on the Einstein tensor and the energy-momentum tensor for a scalar field are given in Appendix C. In double-null coordinates Equation (26), using

_G˜tt+_G˜xx=8π_{T˜dddddddt(total)t}+_{T˜dddddddx(total)x}

one obtains the equation of motion for r,

r(−_r,tt+_r,xx)−_r,t2+_r,x2=^e−2σ1−8π^r2V−^{q2 r2}

where _r,t≡dr/dt _r,t≡dr/dt and other quantities are defined analogously. For simplicity, we define η≡^r2 η≡^r2 and integrate the equation of motion for η, instead [51]. The equation of motion for η can be obtained by rewriting Equation (33) as

−_η,tt+_η,xx=2^e−2σ1−8π^r2V−^{q2 r2}

_Gθθ=8π_{T˜dddddddθ(total)θ Gθθ}=8π_{T˜dddddddθ(total)θ}provides the equation of motion for σ,

−_σ,tt+_σ,xx+_r,tt−_r,xxr+4π_ϕ,t2−_ϕ,x2+_ψ,t2−_ψ,x2χ−2^e−2σV+^{e−2σq2 r4}=0

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′(ϕ)+16κ^T˜(ψ)

−_ψ,tt+_ψ,xx+2r(−_{r,t ψ,t}+_{r,x ψ,x})=23κ(−_{ϕ,t ψ,t}+_{ϕ,x ψ,x})

where

^T˜(ψ)=^g˜μν _T˜μν(ψ)=−^g˜μν _∂˜μψ_∂˜νψχ=1χ^e2σ(_ψ,t2−_ψ,x2)

The {uu} {uu} and {vv} {vv}components of the Einstein equations yield the constraint equations

_r,uu+2_{σ,u r,u}+4πr_ϕ,u2+_ψ,u2χ=0

_r,vv+2_{σ,v r,v}+4πr_ϕ,v2+_ψ,v2χ=0

Via the definitions of u=(t−x)/2 u=(t−x)/2 and v=(t+x)/2 v=(t+x)/2 , the constraint equations can be expressed in (t,x) (t,x) coordinates. Equations (40) − (39) and (40) + (39) generate the constraint equations for the {tx} {tx} and {tt}+{xx} {tt}+{xx}components, respectively,

_r,tx+_{r,t σ,x}+_{r,x σ,t}+4πr_{ϕ,t ϕ,x}+_{ψ,t ψ,x}χ=0

_r,tt+_r,xx+2(_{r,t σ,t}+_{r,x σ,x})+4πr_ϕ,t2+_ϕ,x2+_ψ,t2+_ψ,x2χ=0

3.2. Initial Conditions

We set the initial data to be time symmetric:

_r,t=_σ,t=_ϕ,t=_ψ,t=0dddatdt=0

Therefore, in this configuration, the constraint Equation (41) is satisfied identically. Note that, in this configuration, the values of _{r,t r,t} and _{σ,t σ,t} at t=0 t=0are 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

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

set 1:d_ϕ(x,t)|t=0=_ϕ0

set 2:d_ϕ(x,t)|t=0=_ϕ0+a^{e−(x−_x0)2}

where _{ϕ0 ϕ0} is de Sitter value, defined by ^V′(_ϕ0)=0 ^V′(_ϕ0)=0 . The initial value for σ is defined to be the same as the corresponding one in the Reissner-Nordström case described by Equation (27),

^e−2σ _|t=0=^e−2σ _|t=0RN=_{r+ r−k+2} ^{r2e−2_k+r r_r−−11+_k+|_k−|}

where r is defined by Equation (28) with t=0 t=0. We obtain the initial value for r in charge scattering by combining Equations (33) and (42),

_r,xx=−_{r,t σ,t}−_{r,x σ,x}+_r,t2−_r,x22r−2πr_ϕ,t2+_ϕ,x2+_ψ,t2+_ψ,x2χ+12r^e−2σ1−8π^r2V−^{q2 r2}

We set _r,x=_σ,x=0 _r,x=_σ,x=0 at the origin (x=0,t=0) (x=0,t=0) as in the Reissner-Nordström metric case. The definition domain for the spatial coordinate x is [−_xbxb] [−_xbxb] . Then _{r(x,t)|t=0 r(x,t)|t=0} can be obtained by integrating Equation (48) via the fourth-order Runge-Kutta method from x=0 x=0 to x=±_xb x=±_xb , respectively. The initial values of r, σ, ^{f′ f′} , and ψ are shown in Figure 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 t=0 t=0 , we compute the data at t=Δt t=Δt with a second-order Taylor series expansion as done in Reference [64]. Take the variable ψ as an example,

_{ψ|t=Δt =ψ|t=0}+_{ψ,t |t=0Δt+12ψ,tt|t=0} ^(Δt)2

The values of _{ψ|t=0 ψ|t=0} and _{ψ,t |t=0 ψ,t |t=0} are set up as discussed above, and the value of _{ψ,tt |t=0 ψ,tt |t=0}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 r, σ, ϕ, and ψ at t=0 t=0 described by Equation (43) ensure that the constraint Equation (41) is satisfied. The equation for _{r,xx r,xx} at t=0 t=0 expressed by Equation (48) implies that the constraint Equation (42) is satisfied. Computations of r, σ, ϕ, and ψ at t=Δt t=Δt via a second-order Taylor series expansion, as expressed by Equation (49) for the case of ψ, satisfy all the equations of motion.

Note that the region of x<0 x<0is 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.

3.3. Boundary Conditions

The values of r, σ, ϕ, and ψ at the boundaries of x=±_xb x=±_xb are obtained via extrapolations. In fact, since we are mainly concerned with the dynamics around x=0 x=0 , the boundary conditions will not affect the dynamics in this region, as long as _{xb xb}is large enough.

3.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, Δt=Δx Δt=Δx.

The equations of motion for ϕ (36) and for ψ (37) are coupled. Newton’s iteration method is employed to solve this problem [64]. With the illustration of Figure 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:

_ψupnew=_ψup−G(_ψup)J(_ψup)

where G(_ψup) G(_ψup) is the residual of the differential equation for the function _{ψup ψup} , and J(_ψup) J(_ψup)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.

Figure 4. Numerical evolution scheme.

3.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 Equations (27) and (28). Therefore, we can test our code by comparing the numerical and analytic results in the Reissner-Nordström geometry. Set m=1 m=1 , q=0.7 q=0.7 , and Δx=Δt=¹⁰⁻⁴ Δx=Δt=¹⁰⁻⁴ . We plot the evolutions of r and σ along the slice (x=3×¹⁰⁻⁴,t=t) (x=3×¹⁰⁻⁴,t=t) in Figure 5a. As shown in Figure 5a, 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 Figure 13 and Figure 15, for example. Moreover, the convergence of the constraint Equations (41) and (42) is examined. We assume one constraint equation is nth-order convergent: residual= O(^hn) O(^hn), where h 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=_log2O(^hn)O^h2n

Our numerical results show that both of the constraint equations are about second-order convergent. As a representative, we plot the results for the {tx} {tx} constraint Equation (41) in Figure 5b for the slice (x=x,t=0.65) (x=x,t=0.65).

Figure 5. Tests of numerical code for charge scattering. (a) Numerical vs analytic results for a Reissner-Nordström black hole. m=1 m=1 , q=0.7 q=0.7 , and Δx=Δt=10−4 Δx=Δt=10−4 . The slice is for (x=3Δx,t=t) (x=3Δx,t=t) . This is a special case of charge scattering with contributions of scalar fields being set to zero. Numerical and analytic results match well at an early stage, while at a later stage gravity and electric field become stronger, the numerical evolutions have a time delay, compared to analytic solutions; (b) Numerical tests for the {tx} {tx} constraint Equation (41) and the evolution of ψ on the slice (x=x,t=0.65) (x=x,t=0.65) . They are both second-order convergent.

Convergence tests via simulations with different grid sizes are also implemented [65,66]. If the numerical solution converges, the relation between the numerical solution and the real one can be expressed by

_Freal=^Fh+O(^hn)

where ^{Fh Fh} is the numerical solution. Then, for step sizes equal to h/2 h/2 and h/4 h/4, we have

_Freal=^Fh2+O^h2n,ddd_Freal=^Fh4+O^h4n

Defining _c1≡^Fh−^Fh2 _c1≡^Fh−^Fh2 and _c2≡^Fh2−^Fh4 _c2≡^Fh2−^Fh4, one obtains the convergence rate

n=_{log2c1 c2}

The convergence tests for η≡^r2 η≡^r2 , σ, ϕ, and ψ are investigated. They are all second-order convergent. As a representative, the results for ψ are plotted in Figure 5b for the slice (x=x,t=0.65) (x=x,t=0.65) . The values of the parameters in charge scattering in this section are described at the beginning of Section 6. We use the spatial range of x∈[−1010] x∈[−1010] and the grid spacings of h=Δx=Δt=0.02 h=Δx=Δt=0.02.

4. Neutral Scalar Collapse

In this section, we consider neutral collapse in flat geometry in f(R) f(R)gravity and discuss the mass inflation which happens in the vicinity of the central singularity of the formed black hole.

4.1. Numerical Setup

The numerical setup in neutral scalar collapse in f(R) f(R) gravity is discussed in Reference [20]. The dynamical equations for r, η, σ, ϕ, and ψ can be obtained by setting the terms related to the electric field in the corresponding equations in Section 3.1 to zero:

r(−_r,tt+_r,xx)−_r,t2+_r,x2=^e−2σ(1−8π^r2V)

−_η,tt+_η,xx=2^e−2σ(1−8π^r2V)

−_σ,tt+_σ,xx+_r,tt−_r,xxr+4π_ϕ,t2−_ϕ,x2+_ψ,t2−_ψ,x2χ−2^e−2σV=0

−_ϕ,tt+_ϕ,xx+2r(−_{r,t ϕ,t}+_{r,x ϕ,x})=^e−2σV′(ϕ)+16κ^T˜(ψ)

−_ψ,tt+_ψ,xx+2r(−_{r,t ψ,t}+_{r,x ψ,x})=23κ(−_{ϕ,t ψ,t}+_{ϕ,x ψ,x})

where ^T˜(ψ)=^e2σ(_ψ,t2−_ψ,x2)/χ ^T˜(ψ)=^e2σ(_ψ,t2−_ψ,x2)/χ.

In the equation of motion for σ (54), the term (_r,tt−_r,xx)/r (_r,tt−_r,xx)/r can create big errors near the center x=r=0 x=r=0 . To avoid such a problem, we use the constraint Equation (39) alternatively [51]. Defining a new variable g

g≡−2σ−ln(−_r,u)

one can rewrite Equation (39) as the equation of motion for g,

_g,u=4π·r_r,u·_ϕ,u2+_ψ,u2χ

In the numerical integration, once the value of r at the advanced level is obtained, the value of σ at the current level will be computed using Equation (57).

We set the initial data as

_r,tt=_r,t=_σ,t=_ϕ,t=_ψ,t=0dddatdt=0

The initial values for χ[≡exp(2/3κϕ)] χ[≡exp(2/3κϕ)] and ψ(r) ψ(r)are defined as

_χ(r)|t=0=a·[1−tanh^(r−_r1)2]+_χ0

_ψ(r)|t=0=b·tanh^(r−_r2)2

with a=0.2 a=0.2 , b=0.1 b=0.1 , _r1=_r2=4 _r1=_r2=4 , and ^U′(_χ0)=0 ^U′(_χ0)=0 . The parameters for the Hu-Sawicki model described by Equation (30) are set as D=1.2 D=1.2 and _R0=¹⁰⁻⁶ _R0=¹⁰⁻⁶.

The local Misner-Sharp mass m is defined as [67]

^gμν _{r,μ r,ν}=^e2σ(−_r,t2+_r,x2)≡1−2mr

(See Reference [68] for details on various properties of the Misner-Sharp mass/energy in spherical symmetry.) Then on the initial slice (x=x,t=0) (x=x,t=0) , the equations for r, m, and g are [20,51]

_r,x=1−2mr^eg

_m,r=4π^r2V+121−2mr_ϕ,r2+_ψ,r2χ

_g,r=4πr_ϕ,r2+_ψ,r2χ

Set r=m=g=0 r=m=g=0 at the origin (x=0,t=0) (x=0,t=0) . Then the values of r, m, and g on the initial slice (x=x,t=0) (x=x,t=0) can be obtained by integrating Equations (63)–(65) from the center x=0 x=0 to the outer boundary x=_xb x=_xb via the fourth-order Runge-Kutta method. The values of r, σ, ϕ, and ψ at t=Δt t=Δt can be obtained with a second-order Taylor series expansion, as discussed in Section 3.2. The value of g at t=Δt t=Δt can be obtained using Equation (57).

The range for the spatial coordinate is defined to be x∈[020] x∈[020] . At the inner boundary x=0 x=0 , r is always set to zero. The terms 2(−_{r,t ϕ,t}+_{r,x ϕ,x})/r 2(−_{r,t ϕ,t}+_{r,x ϕ,x})/r in Equation (55) and 2(−_{r,t ψ,t}+_{r,x ψ,x})/r 2(−_{r,t ψ,t}+_{r,x ψ,x})/r in Equation (56) need to be regular at x=r=0 x=r=0 . Since r is always set to zero at the center, so is _{r,t r,t} . Then we enforce ϕ and ψ to satisfy _ϕ,x=_ψ,x=0 _ϕ,x=_ψ,x=0 at x=0 x=0 . The value of g at x=0 x=0 is obtained via extrapolation. We set up the outer boundary conditions at x=20 x=20 via extrapolation. The temporal and the spatial grid spacings are Δt=Δx=0.005 Δt=Δx=0.005.

The numerical code is second-order convergent and is the one developed in Reference [20].

4.2. Black Hole Formation

The evolutions of r, σ, ϕ, and ψ are plotted in Figure 6a–d, respectively. On the apparent horizon of a black hole, the expansion of the outgoing null geodesics orthogonal to the apparent horizon is zero [69]. Then in double-null coordinates, on the apparent horizon, there is [70]

^gμν _{r,μ r,ν}=^e2σ(−_r,t2+_r,x2)=1−2mr=0

Using this property, we locate the apparent horizon and plot it in Figure 6e–f. As shown in Figure 6a, the central singularity is also approached in the collapse. Therefore, a black hole is formed.

Figure 6. Evolutions in neutral collapse for the Hu-Sawicki model Equation (30), f(R)=R−DR0R/(R+R0) f(R)=R−DR0R/(R+R0) . In (a)–(d), the time interval between two consecutive slices is 10Δt=0.05 10Δt=0.05 . (e) and (f) are for the apparent horizon and the singularity curve of the formed black hole.

4.3. Asymptotic Dynamics in the Vicinity of the Central Singularity of the Formed Black Hole

We focus on the dynamics in the vicinity of the central singularity of the formed black hole. We will discuss that, in the vicinity of the singularity, due to the backreaction of the scalar fields on the geometry, the Misner-Sharp mass diverges. In other words, in addition to the inner horizons of Reissner-Nordström and Kerr black holes, mass inflation also happens in the vicinity of the central singularity of a Schwarzschild black hole.

In the vicinity of the central singularity, the field equation can be reduced to the following forms [20]:

r_r,tt≈−_r,t2

_σ,tt≈_r,ttr+4π_ϕ,t2

_ϕ,tt≈−2r_{r,t ϕ,t}

_ψ,tt≈−2r_{r,t ψ,t}+23κ_{ϕ,t ψ,t}

The asymptotic solutions to Equations (67)–(69) are [20]

r≈A^ξβ

σ≈Blnξ+_σ0≈[β(1−β)−4π^C2]lnξ+_σ0

ϕ≈Clnξ

The variable ξ is defined as ξ≡_t0−t ξ≡_t0−t , where _{t0 t0}is the coordinate time on the singularity curve.

As implied in Equation (70), ψ is suppressed (magnified) by ϕ when _{ϕ,t ϕ,t} is negative (positive). Due to the complex competition between gravity and dark energy, an approximate analytic expression for ψ is not obtained. Substituting Equation (71) into Equation (67) yields

(1−β)^ξ2(β−1)≈β^ξ2(β−1)

Then we have

β≈12

Using Equations (11) and (73), there is

χ≈^ξ23κC

So as shown in Figure 6c, when C is positive, χ also approaches zero as ξ and r approach zero. Then the transformation between the Jordan and Einstein frames, _gμν(EF)=χ·_{gμν(JF) gμν(EF)}=χ·_gμν(JF), breaks down. Actually this is not a serious problem, because numerical simulation stops anyway when the central singularity is approached. Moreover, in this paper, we discuss the dynamics in the vicinity of the central singularity rather than on the central singularity.

4.4. Mass Inflation

In the vicinity of the singularity, in the equation of motion for σ (68), because of the contribution from ϕ, σ(x,t) σ(x,t)is greater than the corresponding value in the Schwarzschild black hole case. This makes the mass function divergent near the singularity as will be discussed below.

Near the singularity, using Equation (62), Equation (71)–Equation (74), the mass function can be written as

m=r2[1+^e2σ(_r,t2−_r,x2)]≈18(1−^K2)^{A3 e2_σ0ξ−8πC2}≈18(1−^K2)^{A3+16πC2 e2_σ0r−16πC2}

where K≡|_r,x/_r,t| K≡|_r,x/_r,t| . In the Schwarzschild black hole case, C=0 C=0 . The mass function is always constant and is equal to the black hole mass. In neutral collapse, the parameter β does not change much and remains about 1/2 1/2 . However, the parameter C is not zero. Then the metric quantity σ is modified. (See Equation (72)). As a result, the delicate balance between r and ^e2σ(−_r,t2+_r,x2) ^e2σ(−_r,t2+_r,x2) is broken. Consequently, as implied in Equation (76), near the singularity, the mass function diverges: mass inflation occurs.

During the collapse, before the black hole is formed, the energy of the scalar fields accumulates in the central region. As a result, the scalar fields near x=0 x=0 are stronger than those at large-x regions. Next we discuss three consequences. As a support, we examine the dynamics in the vicinity of the singularity via mesh refinement that was implemented in References [20,71], and plot two sample sets of results on the slices (x=1,t=t) (x=1,t=t) and (x=2,t=t) (x=2,t=t) in Figure 7 and Figure 8, respectively.

Figure 7. (color online) Dynamics on the slice (x=1,t=t) (x=1,t=t) in neutral collapse for the Hu-Sawicki model Equation (30). (a)–(c): dynamical equations for r, η, and σ. (b) Near the central region, due to accumulation, the scalar field ϕ is strong. σ is positive. As a result, in the equation of motion for η (53), the terms 2e−2σ 2e−2σ and 16πe−2σr2V 16πe−2σr2V are negligible. The equation is reduced to η,tt≈η,xx η,tt≈η,xx . (d) ln(mEF)=alnr+b ln(mEF)=alnr+b , a=−1.6438±0.0008 a=−1.6438±0.0008 , b=−0.998±0.003 b=−0.998±0.003 . ln(mJF)=alnr+b ln(mJF)=alnr+b , a=−2.385±0.002 a=−2.385±0.002 , b=−3.443±0.009 b=−3.443±0.009 .

(1)

Values of σ. Due to the backreaction of the scalar fields on the geometry, σ in small-x regions is greater than in large-x regions. In fact, σ is positive in small-x regions, while negative in large-x regions. (See Figure 6b). Our numerical results of the parameter C in ϕ≈Clnξ ϕ≈Clnξ are C≈0.18 C≈0.18 at x=1 x=1 and C≈0.07 C≈0.07 at x=2 x=2 . Then we have 4π^C2≈0.41>1/4 4π^C2≈0.41>1/4 and σ>0 σ>0 at x=1 x=1 , while 4π^C2≈0.06<1/4 4π^C2≈0.06<1/4 and σ<0 σ<0 at x=2 x=2. (See Equation (72)).

(2)

Equation of motion for η (53). For positive σ, the terms 2^e−2σ 2^e−2σ and 16π^{e−2σ r2}V 16π^{e−2σ r2}V in Equation (53) are negligible, compared to the other two. Then Equation (53) is reduced to _η,tt≈_{η,xx η,tt}≈_η,xx . (See Figure 7b). However, for negative σ, the term 2^e−2σ 2^e−2σ is important, and Equation (53) is reduced to _η,tt≈−2^e−2σ _η,tt≈−2^e−2σ . (See Figure 8b).

(3)

Growth of the mass function. As implied in Equation (76), the mass function grows faster in the strong scalar field case than in the weak one. (See Figure 7d and Figure 8d).

Figure 8. (color online) Dynamics on the slice (x=2,t=t) (x=2,t=t) in neutral collapse for the Hu-Sawicki model. (a)–(c): dynamical equations for r, η, and σ. (b) At large-x regions, the scalar field ϕ is weak, and σ is negative. As a result, in the equation of motion for η (53), the term 2e−2σ 2e−2σ is important. The equation is reduced to η,tt≈−2e−2σ η,tt≈−2e−2σ . (d) lnm=alnr+b lnm=alnr+b , a=−0.2673±0.0008 a=−0.2673±0.0008 , b=0.217±0.003 b=0.217±0.003 .

Since f(R) f(R) gravity is defined in the Jordan frame, it is interesting to examine the mass function in the Jordan frame. A generalized Misner-Sharp energy in f(R) f(R) gravity in the Jordan frame was defined in Reference [40]. However, due to the complexity of some integrals, an explicit quasi-local form is usually not available with the exceptions of a Friedmann-Robertson-Walker universe and static spherically symmetric solutions with constant scalar curvature. In this paper, for simplicity, we remain to use the conventional format of the Misner-Sharp function. Considering the transformation that we used, _gμν(EF)=χ·_{gμν(JF) gμν(EF)}=χ·_gμν(JF) , there are ^e2σ _{|JF=χ·^e2σ|EF} ^e2σ _{|JF=χ·^e2σ|EF} and _rJF=^χ−1/2·_{rEF rJF}=^χ−1/2·_rEF . Considering that near the central singularity [20,41]

_KJF≡|_{r,x r,t|JF}≈|_{ϕ,x ϕ,t}|

one can obtain that _KJF≈_{KEF KJF}≈_KEF. Then the mass function in the Jordan frame can be written as

_mJF=_rJF21+^e2_σJF(_r,t2−_r,x2)_|JF≈18^1−23κC2(1−^K2)^{A3 e2_σ0ξ−8πC2−16κC}≈18^1−23κC2(1−^K2)^{A3+16πC2+23κC e2_σ0r−16πC2−23κC}

In the case of C>0 C>0 , due to the factor ^{r−2/3κC r−2/3κC} , _{mJF mJF} is greater than _{mEF mEF} . The _{mJF mJF} along the slice (x=1,t=t) (x=1,t=t) is plotted in Figure 7d.

For stationary black holes (e.g., Schwarzschild and Reissner-Nordström), the Misner-Sharp mass function is always equal to the black hole mass. In spherical symmetry, at spatial infinity, the mass function describes the total energy/mass of an asymptotically flat spacetime [68]. In gravitational collapse case, it means the total mass of the collapsing system.

In the vicinities of the central singularity of a Schwarzschild black hole and the inner horizon of a Reissner-Nordström or Kerr black hole, the dynamics and some quantities are local. The mass function is just a parameter which varies at each point, not giving global information on the black hole mass.

5. Neutral Scalar Scattering

In this section, we consider neutral scattering, in which a neutral scalar field collapses in a Schwarzschild geometry in f(R) f(R) gravity. The numerical formalism is a simpler version of the one in charge scattering that has been constructed in Section 3, and it can be obtained by removing the electric terms in the field equations presented in Section 3.1 and replacing the Reissner-Nordström geometry with a Schwarzschild one.

5.1. A Dark Energy f(R) f(R)Singularity Problem

In neutral scattering, for usual initial conditions of the scalar degree of freedom ^{f′ f′} , ^{f′ f′} asymptotes to zero as the central singularity is approached, which is similar to what happens in the neutral scalar collapse discussed in Section 4. Details are skipped here. On the other hand, when the initial velocity or acceleration of ^{f′ f′} is large enough, ^{f′ f′} can become 1 before the central singularity is approached. As implied in Equation (29), for dark energy f(R) f(R)models, this means that the Ricci scalar becomes infinite, and the simulation breaks down. In this section, we will focus on this singular circumstance.

Figure 9. (color online) Results for neutral scalar collapse in a Schwarzschild geometry for the Hu-Sawicki model. (a) and (b): evolutions of r and f′ f′ . The time interval between two consecutive slices is 30Δt=0.15 30Δt=0.15 . (c) evolutions of r and f′ f′ on the slice (x=−2,t=t) (x=−2,t=t) . (d) dynamical equation for ϕ on the slice (x=−2,t=t) (x=−2,t=t) . ϕ is mainly accelerated by the spatial derivative ϕ,xx ϕ,xx and the geometrical term −2r,tϕ,t/r −2r,tϕ,t/r . Eventually, ϕ and f′ f′ approach 0 and 1, respectively. Then the Ricci scalar becomes singularity, and the simulation stops.

The parameters are set as follows:

(1)

Schwarzschild geometry: m=1 m=1.

(2)

Physical scalar field: _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b , a=0.08 a=0.08 , b=1 b=1 , and _x0=4 _x0=4.

(3)

f(R) f(R) model: f(R)=R−D_R0R/(R+_R0) f(R)=R−D_R0R/(R+_R0) , D=1.2 D=1.2 , and _R0=¹⁰⁻⁵ _R0=¹⁰⁻⁵.

(4)

Scalar degree of freedom: _ϕ(x,t)|t=0=a+b·exp−^(x−_x0)2/c _ϕ(x,t)|t=0=a+b·exp−^(x−_x0)2/c , a=−0.02 a=−0.02 , b=−0.1 b=−0.1 , c=1 c=1 , and _x0=−2 _x0=−2.

(5)

Grid. Spatial range: x∈[−1010] x∈[−1010] . Grid spacings: Δx=Δt=0.005 Δx=Δt=0.005.

The results in this circumstance are plotted in Figure 9. We plot the dynamics of ϕ on the slice (x=−2,t=t) (x=−2,t=t) in Figure 9d, from which one can see that ϕ is mainly accelerated by the spatial derivative _{ϕ,xx ϕ,xx} and the geometrical term −2_{r,t ϕ,t}/r −2_{r,t ϕ,t}/r . Eventually, ϕ and ^{f′ f′} approach 0 and 1, respectively. Then the Ricci scalar becomes singularity, and the simulation stops, as shown in Figure 9a.

5.2. Avoidance of the Singularity Problem

In fact, it has been argued that such a singularity problem can also be caused in cosmology and compact stars [72]. This problem can be avoided by adding an ^{R2 R2} term to the dark energy f(R) f(R) model [42,43,44,45,46,47,48,49]. In this paper, we add the ^{R2 R2}term to the Hu-Sawicki model,

f(R)=R−D_R0RR+_R0+α^R2

In this combined model, at high curvature scale, ^f′≈2αR ^f′≈2αR . So a singular R is pushed to regions where ^{f′ f′} and ϕ are also singular. Then the singularity problem is avoided. The parameters take the same values as in the last subsection. In addition, the new parameter α in the f(R) f(R) model Equation (79) is set to 1.

The numerical results for neutral scattering for this modified model are plotted in Figure 10. As shown in Figure 10b, for this model, ^{f′ f′} can cross 1 without difficulty. What are the limits that ^{f′ f′} and R can reach? As shown in Figure 10d, as the central singularity is approached, there is

_ϕ,tt≈−2r_{r,t ϕ,t}

Figure 10. (color online) Results for neutral scalar collapse in a Schwarzschild geometry in a combined f(R) f(R) model Equation (79). (a) and (b): evolutions of r and f′ f′ . The time interval between two consecutive slices is 30Δt=0.15 30Δt=0.15 . (c) evolutions of r and f′ f′ on the slice (x=−2,t=t) (x=−2,t=t) . (d) dynamical equation for ϕ on the slice (x=−2,t=t) (x=−2,t=t) .

Since _{r,t r,t} is negative, the above equation describes a positive feedback system of _{ϕ,tt ϕ,tt} and _{ϕ,t ϕ,t} : |_ϕ,t| |_ϕ,t| produces more of |_ϕ,tt| |_ϕ,tt| , and in turn |_ϕ,tt| |_ϕ,tt| produces more of |_ϕ,t| |_ϕ,t| . Since _{ϕ,t ϕ,t} is positive, ϕ can be rapidly accelerated to positive infinity. Then ^{f′ f′}and R will go to infinity as the central singularity is approached. This feature surely deserves some attention as discussed below.

The Starobinsky model, f(R)=R+α^R2 f(R)=R+α^R2 , was obtained by taking into account quantum-gravitational effects. It could cause inflation in the early Universe and is conventionally believed to be singularity-free [42]. The combined model is reduced to the Starobinsky model at high curvature scale. It is found that, in neutral scattering for the combined model, a new black hole, including a new central singularity, can be formed. Near the central singularity, gravity dominates other terms, including the potential related to the ^{R2 R2} term, such that the Ricci scalar R can be pushed to infinity by gravity. We also simulate scalar collapse for the Starobinsky model in flat geometry. Similar results are obtained [73]. Therefore, the classical singularity problem, which is present in general relativity, remains in collapse for these models. Further details are skipped.

6. Results for Charge Scattering

In this section, we explore charge scattering: scalar collapse in a Reissner-Nordström geometry. We study the evolutions of the metric components and scalar fields and obtain approximate analytic solutions. We closely compare the dynamics in Schwarzschild black holes, Reissner-Nordström black holes, neutral scalar collapse, and charge scattering.

In this section, the parameters are set as follows:

(1)

Reissner-Nordström geometry: m=1 m=1 , and q=0.7 q=0.7.

(2)

Physical scalar field: _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b , a=0.08 a=0.08 , b=1 b=1 , and _x0=4 _x0=4.

(3)

f(R) f(R) model: f(R)=R−D_R0R/(R+_R0) f(R)=R−D_R0R/(R+_R0) , D=1.2 D=1.2 , and _R0=¹⁰⁻⁵ _R0=¹⁰⁻⁵.

(4)

Scalar degree of freedom: _ϕ(x,t)|t=0=_{ϕ0 ϕ(x,t)|t=0}=_ϕ0 , with ^V′(_ϕ0)=0 ^V′(_ϕ0)=0.

(5)

Grid. Spatial range: x∈[−1010] x∈[−1010] . Grid spacings: Δx=Δt=0.005 Δx=Δt=0.005 for Section 6.1 and Section 6.2, and Δx=Δt=0.0025 Δx=Δt=0.0025 for Section 6.3 and Section 6.4.

6.1. Evolutions

6.1.1. Outline

In this subsection, we describe the evolutions of r, σ, ^{f′ f′} , and ψ that are plotted in Figure 11. Examining the equations of motion (33)–(37) and the numerical results plotted in Figure 13 and Figure 15, one can see that in the charge scattering dynamical system, there are three types of quantities as follows:

(1)

Metric components: r and σ. They contribute as gravity.

(2)

Scalar fields: ϕ and ψ. They contribute as self-gravitating fields.

(3)

Electric field and V(ϕ) V(ϕ) . As implied in Equation (33), they are repulsive forces. However, the numerical results show that in charge scattering, compared to the contributions from other quantities, the contribution from V(ϕ) V(ϕ)is negligible.

Figure 11. Evolutions in charge scattering. (a)–(d): evolutions of r, σ, f′ f′ , and ψ. The time interval between two consecutive slices is 30Δt=0.15 30Δt=0.15 . (e) and (f) are for the apparent horizon and the singularity curve of the black hole. When the scalar field is strong enough (around x=2 x=2 ), the inner horizon can be pushed to the center, and the central singularity becomes spacelike. When the scalar field is weak enough (e.g., −4<x<−1 −4<x<−1 ), the inner horizon does not change much. At the intermediate state (e.g., 0<x<1.5 0<x<1.5 ), the inner horizon contracts to zero, and the central singularity becomes null. The results for the inner horizon, especially for x>2 x>2 , are not that accurate. We are aware that the inner horizon is actually at infinity, while r still can be very close to r− r− when x and t take moderate values.

Furthermore, these quantities can be separated into two sides: the gravitating side (r, σ, ϕ, and ψ) and the repulsive side [electric field and V(ϕ) V(ϕ)]. The dynamics in charge scattering consists mainly of how these variables interact and how the gravitating and repulsive sides compete.

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. Sample slice: (x=1.5,t=t) (x=1.5,t=t) in Figure 11. See Section 6.2.

(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 spatial and temporal derivatives are almost equal. In the case of the center being reached, the central singularity becomes 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 case. Sample slice: (x=0.5,t=t) (x=0.5,t=t) in Figure 11. See Section 6.3 and Section 6.4.

(3)

Type III: critical scattering. This case is on the edge between the above two cases. When the central singularity is reached, it becomes null. Sample slice: (x=1.4225,t=t) (x=1.4225,t=t) in Figure 11. Due to the similarity to that in general relativity discussed in Reference [41], details on this type of scattering in f(R) f(R)gravity are skipped in this paper.

(4)

Type IV: weak scattering. When the scalar field is very weak, the inner horizon does not contract much. Sample case: Figure 19. See Section 7.

(5)

Type V: tiny scattering. When the scalar field is very tiny, the influence of the scalar field on the geometry is negligible. Sample slice: (x=−3,t=t) (x=−3,t=t) in Figure 11.

In this paper, we will discuss Types I, II, and IV.

6.1.2. Causes of Mass Inflation and Evolutions

The local Misner-Sharp mass for a charge black hole hole is

^gμν _{r,μ r,ν}=^e2σ(−_r,t2+_r,x2)=1−2mr+^{q2 r2}

In a Reissner-Nordström geometry, in Kruskal-like coordinates expressed by Equation (27), near the inner horizon, although σ asymptotes to positive infinity, (_r,t2−_r,x2) (_r,t2−_r,x2) is much less than ^{e−2σ e−2σ} . Consequently, ^e2σ(_r,t2−_r,x2) ^e2σ(_r,t2−_r,x2) approaches zero. Then as implied in Equation (81), the mass function m takes a finite value and is equal to the black hole mass. For more details, see Reference [41].

In charge scattering, the equations of motion for η≡^r2 η≡^r2and σ are

−_η,tt+_η,xx=2^e−2σ1−8π^r2V−^{q2 r2}

−_σ,tt+_σ,xx+_r,tt−_r,xxr+4π_ϕ,t2−_ϕ,x2+_ψ,t2−_ψ,x2χ−2^e−2σV+^{e−2σq2 r4}=0

At the beginning, |_ϕ,t| |_ϕ,t| and |_ψ,t| |_ψ,t| may be less than |_ϕ,x| |_ϕ,x| and |_ψ,x| |_ψ,x| , respectively. However, as r decreases toward the central singularity, gravity becomes stronger. Then |_ϕ,t| |_ϕ,t| and |_ψ,t| |_ψ,t| become greater than |_ϕ,x| |_ϕ,x| and |_ψ,x| |_ψ,x| , respectively. (See Figure 15). As a result, in this case, the repulsive “force”, 4π[_ϕ,t2−_ϕ,x2+(_ψ,t2−_ψ,x2)/χ]+^{e−2σ q2}/^r4 4π[_ϕ,t2−_ϕ,x2+(_ψ,t2−_ψ,x2)/χ]+^{e−2σ q2}/^r4 , is greater than the corresponding one, ^{e−2σ q2}/^{r4 e−2σ q2}/^r4 , in a Reissner-Nordström geometry. This makes σ accelerate faster than in the Reissner-Nordström geometry. Consequently, the repulsive force from 2^e−2σ(^q2/^r2−1) 2^e−2σ(^q2/^r2−1) for η≡^r2 η≡^r2 is much weaker than the corresponding value in the Reissner-Nordström geometry. As a result, near the inner horizon, |_r,t| |_r,t| is much greater than the corresponding one in the Reissner-Nordström case. Then (_r,t2−_r,x2) (_r,t2−_r,x2) moves from extremely tiny values in the Reissner-Nordström metric case to moderate values, and r crosses the inner horizon r=_r− r=_r− for the given Reissner-Nordström geometry. With Equation (81), the mass parameter grows dramatically: mass inflation takes place. In other words, regarding the causes of mass inflation in charge scattering, the scalar fields’ backreaction on r is more important than that on σ. The evolutions of r, σ, ^{f′ f′} , and ψ are plotted in Figure 11a–d.

Because of gravity from the black hole and the contribution from the physical scalar field, for some configurations, ^{f′ f′} can decrease to zero as the central singularity is approached. (See Equation (36) and Figure 11, Figure 13, and Figure 15). Nothing is wrong with this circumstance. However, for certain configurations, ^{f′ f′} can be pushed to 1, as shown in Figure 21. Correspondingly, the Ricci scalar R becomes singular. We will discuss the former case in this section and the latter case in Section 8.

The evolution of ψ is plotted in Figure 11d. In the configurations that we consider, as the central singularity is approached, because of the strong suppression from the dark energy scalar ϕ, ψ asymptotes to constant values.

6.1.3. Locations of Horizons

We locate the outer and inner horizons using the following equation,

^gμν _{r,μ r,ν}=^e2σ(−_r,t2+_r,x2)=1−2mr+^{q2 r2}=0

The results are plotted in Figure 11e,f. Due to the absorptions of the energies of ϕ and ψ, the outer horizon increases from the original value of 1.7 1.7 to 3.7 3.7 . Note that the results for the inner horizon, especially at regions where x>2 x>2 , are not that accurate. We are aware that the inner horizon is actually at infinity, while r still can be very close to _{r− r−}even when x and t take moderate values.

6.2. Spacelike Scattering

In spacelike scattering, the scalar field is so strong, such that the inner horizon can contract to zero volume rapidly, and the central singularity converts from timelike into spacelike. Taking the slice (x=1.52,t=t) (x=1.52,t=t) as an example, we plot the evolutions of r, σ, ^{f′ f′} , ψ, and m on this slice in Figure 12. The mass function remains equal to the mass of the original Reissner-Nordström black hole, _m0=1 _m0=1 , until r is very close to _{r− r−} . By then mass inflation takes place. We examine the dynamics near the central singularity via mesh refinement and plot the terms in the dynamical equations in Figure 13.

The strongness of the scalar field causes several consequences as below.

(1)

Motion of r. The quantity r does not decelerate much when it crosses the inner horizon of the given Reissner-Nordström black hole, and it can approach the center. (See Figure 12).

(2)

Nature of the central singularity. The central singularity converts from timelike into spacelike.

(3)

The dynamics in spacelike scattering is similar to that in strong, neutral scalar collapse. The quantity σ takes large positive values, such that in the vicinity of the central singularity, compared to other terms, the term ^{e−2σ q2}/^{r2 e−2σ q2}/^r2 in the equation for r Equation (33) and the term ^{e−2σ q2}/^{r4 e−2σ q2}/^r4 in the equation for σ Equation (35) are negligible. As a result, in the vicinity of the central singularity, the dynamics is similar to that in strong, neutral scalar collapse as expressed by Equations (67)–(70). (See Figure 7 and Figure 13). Then the quantities r, σ, ϕ, and m take similar forms as those in neutral collapse. (See Figure 14).

In both strong, neutral scalar collapse and spacelike scattering, the equation of motion for η is reduced to _η,tt≈_{η,xx η,tt}≈_η,xx . (See Figure 7b and Figure 13b).

Since _{ϕ,t ϕ,t} is negative, the term 2/3κ_{ϕ,t ψ,t} 2/3κ_{ϕ,t ψ,t} in Equation (37) functions as a friction force for ψ. Consequently, compared to ϕ, ψ grows slowly. As the singularity is approached, ψ even approaches constant values. (See Figure 14c,d). As shown in Figure 13e, near the singularity, two major sets of terms can be expressed as

_{r,t ψ,t}≈_{r,x ψ,x},dd_{ϕ,t ψ,t}≈_{ϕ,x ψ,x}

Alternatively,

_{ψ,t ψ,x}≈_{r,x r,t}≈_{ϕ,x ϕ,t}

(4)

Growth of mass function.

In spacelike scattering, the equation of motion for ϕ can be simplified as

_σ,tt≈4π_ϕ,t2

Consequently, with the results obtained in Section 4.3, σ has the following asymptotic solution:

σ≈Blnξ+_σ0≈−4π^C2lnξ+_σ0

with ϕ≈Clnξ ϕ≈Clnξ . Then similar to the mass function Equation (76) in neutral collapse, with Equaution Equation (81), the mass function in spacelike scattering can be written as

m=r21+^{q2 r2}+^e2σ(_r,t2−_r,x2)≈18(1−^K2)^{A3 e2_σ0ξ2B−12}≈18(1−^K2)^{A4(−B+1) e2_σ0r4B−1}≈18(1−^K2)^{A4(4πC2+1) e2_σ0r−16πC2−1}

Numerical results show that, for the sample slice (x=1.52,t=t) (x=1.52,t=t) , near the central singularity, the slope of the singularity curve K is about 0.04 0.04 . As shown in Figure 14d, we linearly fit the numerical results of m via

lnm≈alnr+b

obtaining

a=−16.937±0.003,ddb=−29.04±0.01

Fitting numerical results for σ according to Equation (87) and combining Equations (88) and (89), we obtain

_aanalytic=4B−1=17.0±0.8

_banalytic=ln18(1−^K2)^{A4(−B+1) e2_σ0}=−28±4

Similarly, fitting numerical results for ϕ according to ϕ≈Clnξ ϕ≈Clnξ, we obtain

_aanalytic=−16π^C2−1=18.06±0.01

_banalytic=ln18(1−^K2)^{A4(4πC2+1) e2_σ0}=−28±4

One can see that the above three sets of results match well.

Figure 12. Evolutions along the slice (x=1.52,t=t) (x=1.52,t=t) in one spacelike scattering process. (a) evolutions of r and σ. (b) evolutions of f′ f′ , ψ, and m. In this case, the scalar fields are so strong, such that r can rapidly cross the line of r=r− r=r− and approach zero. Other variables (σ, f′ f′ , ψ, and m) also evolve rapidly after r has crossed the line of r=r− r=r− .

6.3. The Late/fast Stage of Null Scattering

When the scalar fields are less strong, the inner horizon may still contract to zero. However, in this case, the central singularity becomes null rather than spacelike. The equations of motion remain null: they have similar forms as free wave equations, e.g., _ϕ,tt≈_{ϕ,xx ϕ,tt}≈_ϕ,xx.

In the Reissner-Nordström black hole case, near the center, the repulsive (electric) force dominates gravity, and the central singularity is timelike. In spacelike scattering as discussed in the last subsection, gravity from the scalar field and the background geometry dominates the repulsive force. As a result, the central singularity is spacelike. At the late stage of null scattering that will be studied in this subsection, the scalar field is less strong, and the central singularity is null. Because of this, one may say that null scattering is a critical case of the competition between repulsive and gravitational forces, in which case the two types of forces have a balance.

Figure 13. (color online) Dynamics along the slice (x=1.52,t=t) (x=1.52,t=t) near the spacelike, central singularity. Near the singularity, we can approximately rewrite the original equations of motion for r, η, σ, ϕ, and ψ as follows. (a) rr,tt≈−r,t2 rr,tt≈−r,t2 . (b) η,tt≈η,xx2 η,tt≈η,xx2 . (c) σ,tt≈4πϕ,t2 σ,tt≈4πϕ,t2 . (d) ϕ,tt≈−2r,tϕ,t/r ϕ,tt≈−2r,tϕ,t/r . (e) r,tψ,t≈r,xψ,x r,tψ,t≈r,xψ,x and ϕ,tψ,t≈ϕ,xψ,x ϕ,tψ,t≈ϕ,xψ,x . Then we have ψ,t/ψ,x≈r,x/r,t≈ϕ,x/ϕ,t ψ,t/ψ,x≈r,x/r,t≈ϕ,x/ϕ,t .

Figure 14. Solutions along the slice (x=1.52,t=t) (x=1.52,t=t) near the spacelike, central singularity. (a) lnr≈alnξ+b lnr≈alnξ+b , a=0.5032±0.0001 a=0.5032±0.0001 , b=−0.3033±0.0008 b=−0.3033±0.0008 . σ≈alnξ+b σ≈alnξ+b , a=−4.0±0.2 a=−4.0±0.2 , b=−10±2 b=−10±2 . (b) ϕ≈alnξ+b ϕ≈alnξ+b , a=0.5827±0.0002 a=0.5827±0.0002 , b=−1.294±0.002 b=−1.294±0.002 . lnm≈alnr+b lnm≈alnr+b , a=−16.937±0.003 a=−16.937±0.003 , b=−29.04±0.01 b=−29.04±0.01 . (c) and (d): ψ(x,t) ψ(x,t) asymptotes to constant values as the central singularity is approached.

Take the slice (x=0.5,t=t) (x=0.5,t=t) as a sample slice, we plot the terms in the field equations for r, σ, ϕ, and ψ in Figure 15, and the evolutions of r, σ, ϕ, ψ, m, and |1−^K2| |1−^K2| in Figure 16. We investigate the dynamics in the vicinity of the central singularity via mesh refinement and plot the results in Figure 17.

The null scattering has two stages: early/slow and late/fast. As shown in Figure 15 and Figure 16, at the beginning of charge scattering, because of the repulsive force from the electric field, r, σ, ϕ, and ψ evolve slowly. As a result, the mass function m also grows slowly. We call this stage the early/slow stage. Later on, as the center is approached, gravity becomes very strong. Then these quantities evolve faster. We call this stage the late/fast stage.

As shown in Figure 15, when r is very small, the equations of motion for r (33), σ (35), and ϕ (36) can be rewritten as

−r_r,tt≈−r_r,xx≈_r,t2≈_r,x2

_σ,tt≈_σ,xx≈4π_ϕ,t2≈4π_ϕ,x2

_ϕ,tt≈_ϕ,xx≈−2r_{r,t ϕ,t}≈−2r_{r,x ϕ,x}

Figure 15. (color online) Dynamics along the slice (x=0.5,t=t) (x=0.5,t=t) in null scattering. (a)–(f): dynamical equations for r, η, σ, ϕ, and ψ.

Figure 16. Evolutions along the slice (x=0.5,t=t) (x=0.5,t=t) in null scattering. (a)–(f): evolutions for r, σ, f, ψ, m, and |1−K2| |1−K2| . At the early stage of mass inflation, 1.2<t<2 1.2<t<2 , r varies slowly; while at the late stage, t>2 t>2 , r varies rapidly toward zero.

Figure 17. Solutions near the null central singularity along the slice (x=0.5,t=t) (x=0.5,t=t) in null scattering. (a) lnr≈alnξ+b lnr≈alnξ+b , a=0.4421±0.0008 a=0.4421±0.0008 , b=−0.887±0.002 b=−0.887±0.002 . σ≈alnξ+b σ≈alnξ+b , a=−5.15±0.01 a=−5.15±0.01 , b=27.26±0.02 b=27.26±0.02 . ϕ≈alnξ+b ϕ≈alnξ+b , a=0.583±0.002 a=0.583±0.002 , b=0.465±0.005 b=0.465±0.005 . (b) ψ approaches a constant value 0.4 0.4 . lnm≈alnr+b lnm≈alnr+b , a=−25.2±0.1 a=−25.2±0.1 , b=17.9±0.3 b=17.9±0.3 .

Since the above three equations have some similarities to the corresponding ones in spacelike scattering, it is natural to guess that the quantities r, σ, ψ, and m may have expressions similar to those in spacelike scattering. In fact, this guess is verified by the numerical results plotted in Figure 17. Then we have

r≈A^ξ12

σ≈Blnξ+_σ0≈−4π^C2lnξ+_σ0

ϕ≈Clnξ

m≈18(1−^K2)^{A3 e2_σ0ξ2B−12}

≈18(1−^K2)^{A4(−B+1) e2_σ0r4B−1}

≈18(1−^K2)^{A4(4πC2+1) e2_σ0r−16πC2−1}

As shown in Figure 16f, at the late stage, (1−^K2) (1−^K2) is around ^{10−4 10−4}. Due to the similarity to spacelike scattering that we discussed in the last subsection, a comparison between numerical and analytical results for the fast stage of null scattering is skipped.

6.4. The Early/slow Stage of Null Scattering

As shown in Figure 15, at the early/slow stage of null scattering, the equations of motion for r (33), σ (35), ϕ (36), and ψ (37) are reduced as follows:

_r,tt≈_r,xx,dd_r,t2≈_r,x2

_σ,tt≈_σ,xx,dd_ψ,t2≈_ψ,x2,dd_r,tt≈_r,xx,dd_ϕ,t2≈_ϕ,x2

_ϕ,tt≈_ϕ,xx,dd_{r,t ϕ,t}≈_{r,x ϕ,x}

_ψ,tt≈_ψ,xx,dd_{ϕ,t ψ,t}≈_{ϕ,x ψ,x},dd_{r,t ψ,t}≈_{r,x ψ,x}

The above equations are like free scalar wave equations in flat spacetime. The derivatives of one variable (r, σ, ϕ, and ψ) are independent from the derivatives of another. Arbitrary functions of (t+x) (t+x) or (t−x) (t−x) can satisfy the above equations, and in principle, the initial conditions right after the collision between the scalar fields and the inner horizon will decide which function each variable can take. On the other hand, we find that, as shown in Figure 18, the constraint Equations (41) and (42) provide some useful information on the connections between some variables at the early stage of charge scattering:

_{r,t σ,t}≈4π_r−ψ,t2 ^f′

As shown in Figure 16, the quantities r, σ, ^{f′ f′} , ψ, |1−^K2| |1−^K2| , and m change dramatically at the beginning of charge scattering where r≈_r− r≈_r− . Note that, near the central singularity, r, σ, and ^{f′ f′} have approximate analytic expressions in terms of ξ=_t0−t ξ=_t0−t , where _{t0 t0} is the time coordinate of the singularity curve. So it is natural to guess that, at the early stage of charge scattering, the above quantities may also have approximate analytic expressions of ζ=t−_ts ζ=t−_ts , where _{ts ts} is a certain time value related to the early stage of charge scattering. We plot the evolutions of these quantities at the early stage of charge scattering in Figure 18, from which one can see that r, _{r,t r,t}, and σ may have the following approximate analytic expressions:

r≈_r−

_r,t≈a^ζλ

σ≈blnζ+_σ0

In fact, a logarithmic expression for σ is supported by its behavior near the inner horizon in the Reissner-Nordström black hole case. From Equations (28) and (47), one obtains that, as r approaches the inner horizon r=_r− r=_r− , in the case of t≫x t≫x , σ can be approximated by a logarithmic function of t. As shown in Figure 15c, describing the terms in the equation of motion for σ in charge scattering, at the very early stage ( t≈1 t≈1) of the collision between the scalar fields and the inner horizon, compared to those from other terms, the contributions from the terms related to ϕ and ψ are tiny. Therefore, at this stage, the evolution of σ should not be much different from the corresponding one in the Reissner-Nordström geometry.

As shown in Figure 18c,d, ϕ and ψ can be well fitted by power law functions of ζ. Take ϕ as an example,

ϕ≈c^ζd+_ϕ0

We plot (1−^K2) (1−^K2) in Figure 18d and find that ln(1−^K2) ln(1−^K2)can be well fitted linearly with respect to ζ,

ln(1−^K2)≈fζ+h

Figure 18. (color online) Constraint equations and solutions along the slice (x=0.5,t=t) (x=0.5,t=t) at the early/slow stage of null scattering. (a) and (b): constraint Equations (41) and (42). (c) r,t≈a(t+b)c+d r,t≈a(t+b)c+d , a=(−1.67±0.01)×10−3 a=(−1.67±0.01)×10−3 , b=−0.6622±0.0008 b=−0.6622±0.0008 , c=7.897±0.005 c=7.897±0.005 , d=(−2.610±0.003)×10−4 d=(−2.610±0.003)×10−4 . σ≈aln(t+b)+c σ≈aln(t+b)+c , a=34.39±0.03 a=34.39±0.03 , b=−0.245±0.001 b=−0.245±0.001 , c=3.02±0.04 c=3.02±0.04 . ϕ≈a(t+b)c+d ϕ≈a(t+b)c+d , a=(−9.7±2.3)×10−4 a=(−9.7±2.3)×10−4 , b=−0.71±0.03 b=−0.71±0.03 , c=8.3±0.2 c=8.3±0.2 , d=−0.133340±0.000005 d=−0.133340±0.000005 . (d) ψ≈a(t+b)c+d ψ≈a(t+b)c+d , a=(5.0±0.1)×10−3 a=(5.0±0.1)×10−3 , b=−0.362±0.006 b=−0.362±0.006 , c=5.40±0.02 c=5.40±0.02 , d=(1.19±0.02)×10−3 d=(1.19±0.02)×10−3 . ln(1−K2)≈at+b ln(1−K2)≈at+b , a≈−8.73±0.01 a≈−8.73±0.01 , b≈12.71±0.02 b≈12.71±0.02 . lnm≈aln(t+b)+c lnm≈aln(t+b)+c , a=77.77±0.02 a=77.77±0.02 , b=−0.1994±0.0004 b=−0.1994±0.0004 , c=−15.80±0.03 c=−15.80±0.03 .

Currently we do not have derivations for this linear relation. The good thing is that, in the mass function, (1−^K2) (1−^K2) is a minor factor. Therefore, as verified in Figure 18d, the mass function can be reduced to

m=r21+^{q2 r2}+^e2σ(_r,t2−_r,x2)∼_r−2·^e2σ·_r,t2∼_r−2^e2_σ0·^ζ2b·^{a2 ζ2λ}

where a, b, λ, and _{σ0 σ0}are defined in Equations (103) and (104).

We list the fitting results below:

(1)

_r,t≈a^(t+b)c+d _r,t≈a^(t+b)c+d , a=(−1.71±0.01)×¹⁰⁻³ a=(−1.71±0.01)×¹⁰⁻³ , b=(−5.877±0.008)×¹⁰⁻¹ b=(−5.877±0.008)×¹⁰⁻¹ , c=7.880±0.005 c=7.880±0.005 , d=(−2.639±0.003)×¹⁰⁻⁴ d=(−2.639±0.003)×¹⁰⁻⁴.

(2)

σ≈aln(t+b)+c σ≈aln(t+b)+c , a=34.37±0.03 a=34.37±0.03 , b=(−1.72±0.01)×¹⁰⁻¹ b=(−1.72±0.01)×¹⁰⁻¹ , c=3.06±0.04 c=3.06±0.04.

(3)

ϕ≈a^(t+b)c+d ϕ≈a^(t+b)c+d , a=(−2.7±0.6)×¹⁰⁻³ a=(−2.7±0.6)×¹⁰⁻³ , b=(−7.6±0.3)×¹⁰⁻¹ b=(−7.6±0.3)×¹⁰⁻¹ , c=7.4±0.2 c=7.4±0.2 , d=(−1.3336±0.0001)×¹⁰⁻¹ d=(−1.3336±0.0001)×¹⁰⁻¹.

(4)

ψ≈a^(t+b)c+d ψ≈a^(t+b)c+d , a=(5.1±0.1)×¹⁰⁻³ a=(5.1±0.1)×¹⁰⁻³ , b=(−2.91±0.06)×¹⁰⁻¹ b=(−2.91±0.06)×¹⁰⁻¹ , c=5.38±0.02 c=5.38±0.02 , d=(1.21±0.02)×¹⁰⁻³ d=(1.21±0.02)×¹⁰⁻³.

(5)

ln(1−^K2)≈at+b ln(1−^K2)≈at+b , a≈−8.73±0.01 a≈−8.73±0.01 , b≈12.05±0.02 b≈12.05±0.02.

(6)

lnm≈aln(t+b)+c lnm≈aln(t+b)+c , a=77.69±0.02 a=77.69±0.02 , b=(−1.268±0.004)×¹⁰⁻¹ b=(−1.268±0.004)×¹⁰⁻¹ , c=−15.65±0.03 c=−15.65±0.03.

7. Weak Scalar Charge Scattering

In this section, we consider charge scattering with a weak scalar field. Parameter settings in this section are almost the same as those in the last section with the following exceptions:

(1)

Physical scalar field: _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b , a=0.03 a=0.03 , b=1 b=1 , and _x0=4 _x0=4.

(2)

Grid. Spatial range: x∈[−1212] x∈[−1212] . Grid spacings: Δx=Δt=0.002 Δx=Δt=0.002.

Figure 19. Evolutions for charge scattering with a weak scalar field. (a)–(d): evolutions for r, f′ f′ , and ψ. The time interval between two consecutive slices is 120Δt=0.24 120Δt=0.24 . The central singularity is not approached.

As discussed in the above section, in some spacetime regions where the scalar field is strong, the inner horizon can contract to zero volume, and the central singularity becomes spacelike. However, this does not always necessarily happen. After all, it takes energy for the inner horizon to contract. When the scalar field carries less energy, the inner horizon may only contract to a nonzero value. This is confirmed by our numerical results plotted in Figure 19. The evolution of σ is similar to that in strong scalar field case and is skipped. These results are in agreement with the mathematical proof in Reference [23] and the numerical work in Reference [34]. Since in this case the inner horizon is not totally destructed, one needs to reconsider whether the strong cosmic censorship conjecture is valid here. In addition, in Reference [74], it was argued that, when Hawking radiation is taken into account, this censorship may also be violated. Note that in Reference [75] the interior of a Schwarzschild black hole was also discussed with the backreaction from the Hawking radiation being taken into account.

Figure 20. (color online) Dynamics and evolutions on the slice (x=4,t=t) (x=4,t=t) in weak scattering. (a)–(e): dynamical equations for r, η, σ, ϕ, and ψ. (f) evolutions of r and m.

The dynamics for the quantities r, η, σ, ϕ, and ψ are plotted in Figure 20a–e, respectively. The numerical results show that, at the late stage, the field equations for such quantities become null, in the sense that the temporal and spatial derivatives are almost equal, i.e., _σ,tt≈_{σ,xx σ,tt}≈_σ,xx . Moreover, the derivatives have oscillations. As shown in Figure 20f, the mass function keeps growing even as r approaches a constant value. Further details are skipped.

8. Dark Energy f(R) f(R)Singularity Problem in Charge Scattering

Up to this point, in the numerical simulations of charge scattering that we have implemented in this paper, the scalar degree of freedom ^{f′ f′} asymptotes to zero as the center is approached. Note that, as discussed in Section 5, in neutral scattering in dark energy f(R) f(R) gravity, when the initial velocity or acceleration of ^{f′ f′} is large enough, ^{f′ f′}can go to 1 before the central singularity is approached. Consequently, the Ricci scalar R goes to infinity, and the simulation breaks down. Next we will show that such a problem also happens in charge scattering.

In the simulation, the values of the parameters are the same as those described at the beginning of Section 6, including grid spacings Δx=Δt=0.005 Δx=Δt=0.005 . However, the initial value for ϕ takes the following format: _ϕ(x,t)|t=0=aexp[−^(x−_x0)2]+_{ϕ0 ϕ(x,t)|t=0}=aexp[−^(x−_x0)2]+_ϕ0 , with ^V′(_ϕ0)=0 ^V′(_ϕ0)=0 , a=−0.05 a=−0.05 , and _x0=0.3 _x0=0.3 . We plot the numerical results in Figure 21. As shown in Figure 21b, as the inner horizon is approached, ^{f′ f′}goes to 1, and the Ricci scalar R becomes singular. The simulation breaks down.

Now we explore the causes of this singularity problem. As an example, in Figure 21d,e, we plot the terms in the dynamical equation for ϕ on the slice (x=−2,t=t) (x=−2,t=t) , and find that initially because of the contribution from _{ϕ,xx ϕ,xx} , _{ϕ,t ϕ,t} changes from _{ϕ,t |t=0}=0 _{ϕ,t |t=0}=0 to _ϕ,t>0 _ϕ,t>0at late time. Near the inner horizon, there is

_ϕ,tt≈−2r_{r,t ϕ,t}≈_ϕ,xx≈−2r_{r,x ϕ,x}

Then ϕ is accelerated by gravity. After the inner horizon is met, there is _ϕ,tt≈_ϕ,xx>0 _ϕ,tt≈_ϕ,xx>0 . Eventually, ϕ, ^{f′ f′} , and R go to 0, 1, and +∞ +∞, respectively. The simulation breaks down.

Similar to neutral scattering, the singularity problem can be avoided by adding an ^{R2 R2}term to the Hu-Sawicki model,

f(R)=R−D_R0RR+_R0+α^R2

The parameters take the same values as in the last subsection with the following exceptions:

(1)

Grid spacings: Δx=Δt=0.0025 Δx=Δt=0.0025.

(2)

Physical scalar field ψ: _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b _ψ(x,t)|t=0=a·exp−^(x−_x0)2/b , a=0.05 a=0.05 , b=1 b=1 , and _x0=2 _x0=2.

(3)

Scalar degree of freedom ϕ: _ϕ(x,t)|t=0=a·exp[−^(x−_x0)2]+_{ϕ0 ϕ(x,t)|t=0}=a·exp[−^(x−_x0)2]+_ϕ0 , with ^V′(_ϕ0)=0 ^V′(_ϕ0)=0 , a=−0.05 a=−0.05 , and _x0=−2 _x0=−2.

(4)

f(R) f(R) model Equation (109): D=1.2 D=1.2 , _R0=¹⁰⁻⁵ _R0=¹⁰⁻⁵ , and α=1 α=1.

The numerical results for charge scattering for this modified model are plotted in Figure 22. As shown in Figure 22b, for this model, ^{f′ f′}can cross 1 without difficulty. The simulation can run smoothly.

Figure 21. (color online) A singularity problem in charge scattering for the Hu-Sawicki model Equation (30), f(R)=R−DR0R/(R+R0) f(R)=R−DR0R/(R+R0) . (a) and (b): evolutions of r and f′ f′ . The time interval between two consecutive slices is 20Δt=0.01 20Δt=0.01 . (c) evolutions of r, f′ f′ , and f,t′ f,t′ on the the slice (x=−2,t=t) (x=−2,t=t) . (d) and (e): dynamical equation for ϕ on the slice (x=−2,t=t) (x=−2,t=t) . In (e), as the inner horizon is approached, ϕ,tt≈−2r,tϕ,t ϕ,tt≈−2r,tϕ,t . This equation describes a positive feedback, since −2r,t/r −2r,t/r is positive. As a result, when ϕ,t ϕ,t is positive, ϕ can be accelerated to zero rapidly. Correspondingly, f′ f′ goes to 1 as plotted in (b) and (c), and the Ricci scalar R becomes singular. Then the simulation breaks down as shown in (a).

Figure 22. (color online) Avoidance of the singularity problem in charge scattering for the combined model Equation (109), f(R)=R−DR0R/(R+R0)+αR2 f(R)=R−DR0R/(R+R0)+αR2 . (a) and (b): evolutions of r and f′ f′ . The time interval between two consecutive slices is 60Δt=0.15 60Δt=0.15 . (c) dynamical equation for ϕ on the slice (x=−2,t=t) (x=−2,t=t) . (d) evolution of m on the slice (x=−2,t=t) (x=−2,t=t) .

9. Summary

In this paper, we studied scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in f(R) f(R) gravity numerically. Approximate analytic solutions for different types of collapses were partially obtained. One dark energy f(R) f(R)singularity problem was discussed. We summarize our work on computational and physical issues separately below.

9.1. Computational Issues

(1)

The Jordan frame vs the Einstein frame. The field equations for f(R) f(R) gravity in the Jordan frame are more complex than those in general relativity. Therefore, for ease of computation, we transform f(R) f(R)gravity from the Jordan frame into the Einstein frame, in which the formalism can be formally treated as Einstein gravity coupled to a scalar field.

(2)

dudv dudv vs. (−d^t2+d^x2) (−d^t2+d^x2) in double-null coordinates. In the studies of mass inflation, the dudv dudv format of the Kruskal-like coordinates, d^s2=4^e−2σdudv+^r2d^Ω2 d^s2=4^e−2σdudv+^r2d^Ω2 , is usually used. In the field equations, many terms are mixed derivatives of u and v, e.g., _{r,uv r,uv} . In this paper, we used the (−d^t2+d^x2) (−d^t2+d^x2) format instead, d^s2=^e−2σ(−d^t2+d^x2)+^r2d^Ω2 d^s2=^e−2σ(−d^t2+d^x2)+^r2d^Ω2 , with u=(t−x)/2=const u=(t−x)/2=const and v=(t+x)/2=const v=(t+x)/2=const . In the (t,x) (t,x) line element, one coordinate is timelike, and the rest are spacelike. We are used to this setup. It is more convenient and more intuitive to use this set of coordinates. Moreover, for the (t,x) (t,x) choice, spatial and temporal derivatives are usually separated, e.g., (_r,tt−_r,xx) (_r,tt−_r,xx).

We set the initial conditions close to those in a Reissner-Nordström geometry. With this setup, it is convenient to test the code. Removing the terms related to the scalar fields, we can test our code by comparing the numerical results to the analytic ones in a Reissner-Nordström geometry. Moreover, by comparing numerical results for charge scattering to the dynamics in the Reissner-Nordström geometry, we can obtain intuitions as to how the scalar fields affect the geometry.

(3)

Cauchy horizon: infinite or local regions? As implied by Equation (28), the exact inner horizon r=_r− r=_r− is at the regions where uv uv and (^t2−^x2) (^t2−^x2) are infinite. However, r still can be very close to the inner horizon even when uv uv and (^t2−^x2) (^t2−^x2) take moderate values. Consequently, at regions where uv uv and (^t2−^x2) (^t2−^x2)take some moderate values, the scalar fields and the inner horizon still can have strong interactions, resulting in mass inflation.

9.2. Physical Issues

(1)

Scalar collapse in f(R) f(R) gravity vs scalar collapse in general relativity. In scalar collapse, the scalar degree of freedom ϕ(≡3/2ln^f′/8πG) ϕ(≡3/2ln^f′/8πG) plays a similar role as a physical scalar field in general relativity. Regarding the physical scalar field in the f(R) f(R) case, when _{ϕ,t ϕ,t}is negative (positive), the physical scalar field is suppressed (magnified) by ϕ.

(2)

The inner horizon in a Reissner-Nordström black hole vs the central singularity in a Schwarzschild black hole. These two share some similarities.

For Reissner-Nordström and Schwarzschild black holes, throughout the whole spacetime, the Misner-Sharp mass function is constant. When a scalar field impacts the inner horizon of a Reissner-Nordström black hole, the scalar field can modify the geometry in the vicinity of the inner horizon significantly, especially on _{r,t r,t} . The inner horizon contracts and mass inflation takes place. In neutral scalar collapse toward a Schwarzschild black hole formation, the scalar field can also modify the geometry in the vicinity of the central singularity dramatically, especially on the metric component σ [20]. Then mass inflation also happens.

The Belinskii, Khalatnikov, and Lifshitz (BKL) conjecture is an important result on dynamics in the vicinity of a spacelike singularity [76,77,78,79]. The first statement of this conjecture is that as the singularity is approached, the dynamical terms dominate the spatial terms in the field equations. In other words, the way gravity changes over time is more important than the variation of the gravitational field from one location to the next [79]. We would like to say that, to a large extent, later evolutions in a strong gravitational field largely erase away the initial information on the connections between neighboring points. As discussed in Reference [20] and also in this paper, in double-null coordinates, using the above argument, one can interpret the following behaviors displayed in numerical simulations: near the central singularity of a Schwarzschild black hole and also near the inner horizon of a Reissner-Nordström black hole, there are

_{ψ,x ψ,t}≈_{r,x r,t}<1

_ϕ,tt≈−2_r,tr_ϕ,t

In this paper, it was shown that Equation (110) can explain the causes of mass inflation, while Equation (111) can explain the dark energy f(R) f(R)singularity problem in collapse.

The second and third statements of the BKL conjecture are that i) the metric terms will dominate the matter field terms, while the matter field may not be negligible if it is a scalar field; ii) the dynamics of the metric components and the matter fields is described by the Kasner solution. These two statements were confirmed in simulations of neutral scalar collapse in f(R) f(R) gravity in Reference [20] and in general relativity in Reference [41]. The second statement was also verified in charge scattering in this paper. However, the third statement on Kasner solution may not apply to the dynamics near the inner horizon in charge scattering.

(3)

Compact stars vs black holes in f(R) f(R) gravity. The internal structure of compact stars is usually in an equilibrium state and is static. In f(R) f(R) gravity, inside compact stars, the scalar degree of freedom, ^{f′ f′} , can be coupled to the energy density of the stars and then is not very free to move. However, due to strong gravity, the internal structure of black holes is dynamical. ^{f′ f′} and the matter fields are decoupled. As a result, ^{f′ f′}is more free to move than in the compact stars case. It can keep increasing or decreasing until singularities are met.

(4)

Dark energy f(R) f(R) singularity problem: cosmology (or static compact objects) vs black hole physics. In dark energy f(R) f(R)gravity, the Ricci scalar R can be singular in both cosmology and black hole physics. We consider a homogeneous cosmological model. Using the flat Friedmann-Robertson-Walker metric,

d^s2=−d^t2+^a2(t)d^x2

the equation of motion for ^{f′ f′}is

^f′′¨+3H^f′˙+U(^f′)+8π3T=0

where H is the Hubble parameter. Due to the finiteness of the potential barrier, the force from a perturbation of 8πT/3 8πT/3 may push ^{f′ f′} to 1. Correspondingly, the Ricci scalar goes to singularity [72]. In the black hole case that we discussed in this paper, the equation of motion for ϕ(≡3/2ln^f′/8πG) ϕ(≡3/2ln^f′/8πG) has more complex structure Equation (36),

−_ϕ,tt+_ϕ,xx+2r(−_{r,t ϕ,t}+_{r,x ϕ,x})=^e−2σV′(ϕ)+16κ^T(ψ)

As the central singularity of a Schwarzschild black hole or the inner horizon of a Reissner-Nordström black hole is approached, the above equation can be simplified as Equation (111), and gravity from the black hole, −2_{r,t ϕ,t}/r −2_{r,t ϕ,t}/r, can cause a similar singularity problem as in cosmology or static compact stars.

(5)

The combined f(R) f(R) model and the ^{R2 R2} model: singular or non-singular? At the center of a Schwarzschild black hole and in the very early Universe, the tidal forces are singular, and general relativity fails. Taking into account quantum-gravitational effects, Starobinsky obtained an ^{R2 R2} model, f(R)=R+α^R2 f(R)=R+α^R2 . This model has a non-singular de Sitter solution, which is unstable both to the past and to the future [42,43]. In Section 5, scalar collapse in a Schwarzschild geometry for the combined model (a combination of dark energy model and ^{R2 R2} model) was explored. A new Schwarzschild black hole, including a new central singularity, can be formed. Moreover, under certain initial conditions, ^{f′ f′} and R can be pushed to infinity as the central singularity is approached. In Reference [73], scalar collapse in flat geometry for the ^{R2 R2}model was simulated. Similar results were obtained. Namely, the classical singularity problem, which is present in general relativity, remains in collapse in these models.

(6)

Inside vs. outside black holes: local vs global. Throughout the whole spacetime of stationary Schwarzschild and Reissner-Nordström black holes, the Misner-Sharp mass function is equal to the black hole mass. For a gravitational collapsing system, at asymptotic flat regions, the mass function describes the total mass of the dynamical system. However, in this system, near the central singularity of a Schwarzschild black hole or near the inner horizon of a Reissner-Nordström black hole, the dynamics is local. Then the mass function does not provide global information on the mass of the collapsing system.

In summary, in this paper, we studied scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in f(R) f(R) gravity. For convenience and intuitiveness, in simulating scalar collapses in Schwarzschild and Reissner-Nordström geometries, Kruskal and Kruskal-like coordinates were used, respectively. Approximate analytic solutions for different types of collapses were partially obtained. Causes and avoidance of a dark energy f(R) f(R)singularity problem in collapse were discussed.

Acknowledgments

The authors would like to thank Andrei V. Frolov, José T. Gálvez Ghersi, Ken-ichi Nakao, and Dong-han Yeom for useful discussions. The authors also thank the referees for the helpful comments.

Author Contributions

All authors have contributed to this paper. All authors have read and approved the final version.

Conflicts of Interest

The authors declare on conflict of interest.

A. Equations of Motion for a Physical Scalar Field and ^{f′ f′}in the Einstein Frame

In this appendix, based on the equations of motion for a massless physical scalar field ψ and the scalar degree of freedom ^{f′ f′} in the Jordan frame for f(R) f(R) gravity, we derive the corresponding equations in the Einstein frame. In the transformation from the Jordan frame into the Einstein frame, we use _g˜μν=χ·_{gμν g˜μν}=χ·_gμν and κϕ≡3/2lnχ κϕ≡3/2lnχ, where a tilde denotes that the quantity is in the Einstein frame.

For a scalar field ψ, the first covariant derivatives of ψ are equal in two frames, since they are both equal to the partial derivative [80]:

_∇˜μψ=_∇μψ=_∂μψ

Then for the first contravariant derivative of ψ, there is

^∇μψ=^gμν _∇νψ=χ·^∇˜μψ

For □ψ □ψ , we have [81]

□ψ=1|g|_∂μ|g|^gμν _∂νψ=1|g˜|·^χ−4_∂μ|g˜|·^χ−4·χ·^g˜μν _∂νψ=χ□˜ψ−^g˜μν _∂μψ_∂ν(lnχ)=χ□˜ψ−23κ^g˜μν _∂μϕ_∂νψ

In the Jordan frame, for a massless scalar field ψ, there is

□ψ=0

Combining Equations (117) and (118) yields the equation of motion for ψ in the Einstein frame Equation (22).

In the special case of ψ=χ≡exp(2/3κϕ) ψ=χ≡exp(2/3κϕ), there is

□χ=23κ^χ2□˜ϕ

Combining Equations (9), (10), (17), (24), and (119) gives the equation of motion for ϕ in the Einstein frame Equation (21).

B. Reissner-Nordström Metric in Kruskal-Like Coordinates

The Reissner-Nordström metric in Kruskal-like coordinates was obtained in Reference [52] (also see Reference [53]). For reference concern, we list the derivations here via analog to the Schwarzschild black hole case in the form of Reference [53].

We start from the Reissner-Nordström metric in the conventional form,

d^s2=−▵d^t2+^▵−1d^r2+^r2d^Ω2

where

▵=1−2mr+^{q2 r2}=(r−_r+)(r−_r−)^r2

_r±=m±^m2−^q2

Considering the radial null curves along which d^s2=0 d^s2=0, we have

dtdr=±1−2mr+^{q2 r2}

The solution to the above equation can be expressed as

t=±^r*+const

where

^r*=r+12_k+ln|1−r_r+|+12_k−ln|1−r_r−|

_k±=_r±−_r∓2_r±2

Define

u=t−^r*,hhv=t+^r*

and then the line element Equation (120) becomes

d^s2=▵(−d^t2+d^r*2)+^r2d^Ω2=−12▵(dudv+dvdu)+^r2d^Ω2

B.1. Patch I: _r−r+∞ _r−r+∞

(i) _r+r+∞ _r+r+∞of Patch I

As implied in Equation (125), the inner and outer horizons in the metric Equation (128) are pushed to infinity. To compact them into finite space, we define

^U+=−^e−_k+u,hh^V+=^e_k+v.ddforr>_r+

Rewrite Equation (125) as

^r*=r+12_k+ln|1−r_r+|r_r−−1−12_k+1+_k+|_k−|ln|1−r_r−|

Using Equations (121), (127), (129), and (130), the line element Equation (128) can be converted into the following format

d^s2=−_{r+ r−k+2} ^{r2e−2_k+r r_r−−11+_k+|_k−|}d^U+d^V++^r2d^Ω2

Note that we get used to coordinate systems in which one coordinate is timelike and the rest are spacelike. So defining

T≡^V++^U+2=^{e_k+ r*}sinh(_k+t)

X≡^V+−^U+2=^{e_k+ r*}cosh(_k+t)

one can rewrite the metric Equation (131) as

d^s2=_{r+ r−k+2} ^{r2e−2_k+r r_r−−11+_k+|_k−|}(−d^T2+d^X2)+^r2d^Ω2

Moreover, using Equations (125), (132), and (133), r and t can be expressed as

^T2−^X2=^e2_k+r1−r_r+^{r_r−−1−_k+|_k−|}

TX=tanh(_k+t)

(ii) _r−r_{r+ r−}r_r+of Patch I

Define

^U+=−^e−_k+u,hh^V+=−^e_k+v.ddfor_r−r_r+

Similar to the case of _r+r+∞ _r+r+∞, one can obtain

T=^V++^U+2=^{e_k+ r*}cosh(_k+t)

X=^V+−^U+2=^{e_k+ r*}sinh(_k+t)

TX=coth(_k+t)

The line element is the same as Equation (134), and the expression for r is the same as Equation (135).

B.2. Patch II: 0r_r+ 0r_r+

(i) _r−r_{r+ r−}r_r+of Patch II

Define

^U−=−^e−_k−u,hh^V−=−^e_k−v.ddfor_r−r_r+

In a similar routine as in Patch I, the metric for the spacetime of _r−r_{r+ r−}r_r+can be obtained as

d^s2=_{r+ r−k−2} ^{r2e2|_k−|r 1−r_r+1+|_k−|_k+}(−d^T2+d^X2)+^r2d^Ω2

with

^T2−^X2=^{e−2|_k−|r}r_r−−1^{1−r_r+−|_k−|_k+}

TX=coth(_k−t)

(ii) 0r_r− 0r_r−of Patch II

Define

^U−=−^e−_k−u,hh^V−=^e_k−v.ddfor0r_r−

Then we have

TX=tanh(_k−t)

The line element is the same as Equation (142), and the expression for r is the same as Equation (143).

C. Einstein Tensor and Energy-Momentum Tensor for a Massive Scalar Field

In this appendix, we give specific expressions for the Einstein tensor and the energy-momentum tensor for a massive scalar field. In double-null coordinates expressed by Equations (25) and (26), some components of the Einstein tensor can be expressed as follows:

_Gtt=2^e2σr2r(_{r,t σ,t}+_{r,x σ,x})+r_rxx+12(−^_r,t2+^_r,x2)−12^e−2σ

_Gxx=2^e2σr2−r(_{r,t σ,t}+_{r,x σ,x})−r_rtt+12(−^_r,t2+^_r,x2)−12^e−2σ

_Gθθ=_Gϕϕ=^e2σr−_r,tt+_r,xx−r(−_σ,tt+_σ,xx)

_Guu=−2r(_r,uu+2_{σ,u r,u})

_Gvv=−2r(_r,vv+2_{σ,v r,v})

For a massive scalar field with an energy-momentum tensor

_Tμν=_{ψ,μ ψ,ν}−12^gαβ _{ψ,α ψ,β}+V(ψ)_gμν

there are

_Ttt=−^e2σ12(_ψ,t2+_ψ,x2)+^e−2σV(ψ)

_Txx=^e2σ12(_ψ,t2+_ψ,x2)−^e−2σV(ψ)

_Tθθ=_Tϕϕ=−^e2σ12(−_ψ,t2+_ψ,x2)+^e−2σV(ψ)

_Tuu=_ψ,u2

_Tvv=_ψ,v2

T≡_Tαα=−^e2σ(−_ψ,t2+_ψ,x2)−4V(ψ)

The equations obtained in this appendix can be used to derive the field equations discussed in Section 3.1.

References

1. Burko, L.M., Ori, A., Eds.; Internal Structure of Black Holes and Spacetime Singularities; Institute of Physics Publishing: Bristol, UK, 1998.

2. Brady, P.R. The Internal Structure of Black holes. Prog. Theor. Phys. Suppl. 1999, 136, 29–44.

3. Berger, B.K. Numerical Approaches to Spacetime Singularities. Living Rev. Relativ. 2002, 5, 1.

4. Joshi, P.S. Gravitational Collapse and Spacetime Singularities; Cambridge University Press: Cambridge, UK, 2007.

5. Henneaux, M.; Persson, D.; Spindel, P. Spacelike Singularities and Hidden Symmetries of Gravity. Living Rev. Relativ. 2008, 11, 1.

6. Price, R.H. Nonspherical Perturbations of Relativistic Gravitational Collapse. I. Scalar and Gravitational Perturbations. Phys. Rev. D 1972, 5, 2419–2438.

7. Clifton, T.; Ferreira, P.G.; Padilla, A.; Skordis, C. Modified gravity and cosmology. Phys. Rep. 2012, 513, 1–189.

8. Sotiriou, T.P.; Faraoni, V. f(R) Theories Of Gravity. Rev. Mod. Phys. 2010, 82, 451–497.

9. Felice, A.D.; Tsujikawa, S. f(R) Theories. Living Rev. Relativ. 2010.

10. Nojiri, S.; Odintsov, S. D. Unified cosmic history in modified gravity: From F(R) theory to Lorentz non-invariant models. Phys. Rept. 2011, 505, 59–144.

11. Capozziello, S.; de Laurentis, M. Extended Theories of Gravity. Phys. Rept. 2011, 509, 167–321.

12. Cruz-Dombriz, A.; de la Dobado, A.; Maroto, A.L. Black Holes in f(R) theories. Phys. Rev. D 2009, 80, 124011.

13. Olmo, G.J.; Rubiera-Garcia, D. Palatini f(R) Black Holes in Nonlinear Electrodynamics. Phys. Rev. D 2011, 84, 124059.

14. Nojiri, S.; Odintsov, S.D. Anti-Evaporation of Schwarzschild-de Sitter Black Holes in F(R) gravity. Class. Quantum Gravity 2013, 30, 125003.

15. Sebastiani, L.; Momeni, D.; Myrzakulov, R.; Odintsov, S.D. The instabilities and (anti)-evaporation of Schwarzschild-de Sitter black holes in modified gravity. Phys. Rev. D 2013, 88, 104022.

16. Nojiri, S.; Odintsov, S.D. Instabilities and Anti-Evaporation of Reissner-Nordström Black Holes in modified F(R) gravity. Phys. Lett. B 2014, 735, 376–382.

17. Avelino, P.P.; Hamilton, A.J.S.; Herdeiro, C.A.R. Mass Inflation in Brans-Dicke gravity. Phys. Rev. D 2009, 79, 124045.

18. Borkowska, A.; Rogatko, M.; Moderski, R. Collapse of Charged Scalar Field in Dilaton Gravity. Phys. Rev. D 2011, 83, 084007.

19. Hwang, D.-I.; Lee, B.-H.; Yeom, D.-H. Mass inflation in f(R) gravity: A conjecture on the resolution of the mass inflation singularity. J. Cosmol. Astropart. Phys. 2011, 2011, 006.

20. Guo, J.-Q.; Wang, D.; Frolov, A.V. Spherical collapse in f(R) gravity and the Belinskii-Khalatnikov-Lifshitz conjecture. Phys. Rev. D 2014, 90, 024017.

21. Simpson, M.; Penrose, R. Internal instability in a Reissner-Nordström black hole. Int. J. Theor. Phys. 1973, 7, 183–197.

22. Dafermos, M. Stability and instability of the Cauchy horizon for the spherically symmetric Einstein-Maxwell-scalar field equations. Ann. Math. (N. Y.) 2003, 158, 875–928.

23. Dafermos, M. Black holes without spacelike singularities. Commun. Math. Phys. 2014, 332, 729–757.

24. Ringström, H. Origins and development of the Cauchy problem in general relativity. Class. Quantum Gravity 2015, 32, 124003.

25. Isenberg, J. On Strong Cosmic Censorship. 2015. arXiv:1505.06390[gr-qc].

26. Poisson, E.; Israel, W. Inner-horizon instability and mass inflation in black holes. Phys. Rev. Lett. 1989, 63, 1663–1666.

27. Poisson, E.; Israel, W. Internal structure of black holes. Phys. Rev. D 1990, 41, 1796–1809.

28. Barrabes, C.; Israel, W.; Poisson, E. Collision of light-like shells and mass inflation in rotating black holes. Class. Quantum Gravity 1990, 7, L273–L278.

29. Gnedin, N.Y.; Gnedina, M.L. Instability of the internal structure of a Reissner-Nordström black hole. Sov. Astron. 1992, 36, 296, (Astron. Zh. 1992, 69, 584–592.).

30. Gnedin, M.L.; Gnedin, N.Y. Destruction of the Cauchy horizon in the Reissner-Nordstrom black hole. Class. Quantum Gravity 1993, 10, 1083–1102.

31. Brady, P.R.; Smith, J. Black hole singularities: A numerical approach. Phys. Rev. Lett. 1995, 75, 1256–1259.

32. Burko, L.M. Structure of the black hole’s Cauchy horizon singularity. Phys. Rev. Lett. 1997, 79, 4958–4961.

33. Burko, L.M.; Ori, A. Late-time evolution of nonlinear gravitational collapse. Phys. Rev. D 1997, 56, 7820–7832.

34. Hansen, J.; Khokhlov, A.; Novikov, I. Physics of the interior of a spherical, charged black hole with a scalar field. Phys. Rev. D 2005, 71, 064013.

35. Hod, S.; Piran, T. Mass Inflation in Dynamical Gravitational Collapse of a Charged Scalar Field. Phys. Rev. Lett. 1998, 81, 1554–1557.

36. Oren, Y.; Piran, T. Collapse of charged scalar fields. Phys. Rev. D 2003, 68, 044013.

37. Burko, L.M.; Ori, A. Analytic study of the null singularity inside spherical charged black holes. Phys. Rev. D 1998, 57, R7084–R7088.

38. Burko, L.M. Strength of the null singularity inside black holes. Phys. Rev. D 1999, 60, 104033.

39. Hu, W.; Sawicki, I. Models of f(R) Cosmic Acceleration that Evade Solar-System Tests. Phys. Rev. D 2007, 76, 064004.

40. Cai, R.G.; Cao, L.M.; Hu, Y.P.; Ohta, N. Generalized Misner-Sharp Energy in f(R) Gravity. Phys. Rev. D 2009, 80, 104016.

41. Guo, J.-Q.; Joshi, P.S. Interior dynamics of neutral and charged black holes. 2015. arXiv:1507.01806[gr-qc].

42. Starobinsky, A.A. A new type of isotropic cosmological models without singularity. Phys. Lett. B 1980, 91, 99–102.

43. Vilenkin, A. Classical and quantum cosmology of the Starobinsky inflationary model. Phys. Rev. D 1985, 32, 2511–2521.

44. Nojiri, S.; Odintsov, S.D. The future evolution and finite-time singularities in F(R)-gravity unifying the inflation and cosmic acceleration. Phys. Rev. D 2008, 78, 046006.

45. Bamba, K.; Nojiri, S.; Odintsov, S.D. Future of the universe in modified gravitational theories: Approaching to the finite-time future singularity. J. Cosmol. Astropart. Phys. 2008, 10, 045.

46. Capozziello, S.; de Laurentis, M.; Nojiri, S.; Odintsov, S.D. Classifying and avoiding singularities in the alternative gravity dark energy models. Phys. Rev. D 2009, 79, 124007.

47. Appleby, S.A.; Battye, R.A.; Starobinsky, A.A. Curing singularities in cosmological evolution of F(R) gravity. J. Cosmol. Astropart. Phys. 2010, 06, 005.

48. Elizalde, E.; Nojiri, S.; Odintsov, S.D.; Sebastiani, L.; Zerbini, S. Non-singular exponential gravity: A simple theory for early- and late-time accelerated expansion. Phys. Rev. D 2011, 83, 086006.

49. Bamba, K.; Nojiri, S.; Odintsov, S.D. Time-dependent matter instability and star singularity in F(R) gravity. Phys. Lett. B 2011, 698, 451–456.

50. Poisson, E. A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics; Cambridge University Press: Cambridge, UK, 2004.

51. Frolov, A.V. Is It Really Naked? On Cosmic Censorship in String Theory. Phys. Rev. D 2004, 70, 104023.

52. Graves, J.C.; Brill, D.R. Oscillatory Character of Reissner-Nordström Metric for an Ideal Charged Wormhole. Phys. Rev. 1960, 120, 1507–1513.

53. Reall, H. Lecture Notes on Black Holes. Available online: http://www.damtp.cam.ac.uk/user/hsr1000/black_holes_lectures_2015.pdf (accessed on 29 August 2015).

54. Nunez, A.; Solganik, S. The content of f(R) gravity. 2004. arXiv:hep-th/0403159.

55. Dolgov, A.D.; Kawasaki, M. Can modified gravity explain accelerated cosmic expansion? Phys. Lett. B 2003, 573, 1–4.

56. Amendola, L.; Gannouji, R.; Polarski, D.; Tsujikawa, S. Conditions for the cosmological viability of f(R) dark energy models. Phys. Rev. D 2007, 75, 083504.

57. Guo, J.-Q.; Frolov, A.V. Cosmological dynamics in f(R) gravity. Phys. Rev. D 2013, 88, 124036.

58. Khoury, J.; Weltman, A. Chameleon Fields: Awaiting Surprises for Tests of Gravity in Space. Phys. Rev. Lett. 2004, 93, 171104.

59. Khoury, J.; Weltman, A. Chameleon Cosmology. Phys. Rev. D 2004, 69, 044026.

60. Chiba, T.; Smith, T.L.; Erickcek, A.L. Solar System constraints to general f(R) gravity. Phys. Rev. D 2007, 75, 124014.

61. Tamaki, T.; Tsujikawa, S. Revisiting chameleon gravity - thin-shells and no-shells with appropriate boundary conditions. Phys. Rev. D 2008, 78, 084028.

62. Tsujikawa, S.; Tamaki, T.; Tavakol, R. Chameleon scalar fields in relativistic gravitational backgrounds. J. Cosmol. Astropart. Phys. 2009, 05, 020.

63. Guo, J.-Q. Solar system tests of f(R) gravity. Int. J. Mod. Phys. D 2014, 23, 1450036.

64. Pretorius, F. Numerical Relativity Using a Generalized Harmonic Decomposition. Class. Quantum Gravity 2005, 22, 425–452.

65. Sorkin, E.; Piran, T. Effects of Pair Creation on Charged Gravitational Collapse. Phys. Rev. D 2001, 63, 084006.

66. Golod, S.; Piran, T. Choptuik’s Critical Phenomenon in Einstein-Gauss-Bonnet Gravity. Phys. Rev. D 2012, 85, 104015.

67. Misner, C.W.; Sharp, D.H. Relativistic Equations for Adiabatic, Spherically Symmetric Gravitational Collapse. Phys. Rev. 1964, 136, B571–B576.

68. Hayward, S.A. Gravitational energy in spherical symmetry. Phys. Rev. D 1996, 53, 1938–1949.

69. Baumgarte, T.W.; Shapiro, S.L. Numerical Relativity: Solving Einstein’s Equations on the Computer; Cambridge University Press: Cambridge, UK, 2010.

70. Csizmadia, P.; Racz, I. Gravitational collapse and topology change in spherically symmetric dynamical systems. Class. Quantum Gravity 2010, 27, 015001.

71. Garfinkle, D. Choptuik scaling in null coordinates. Phys. Rev. D 1995, 51, 5558–5561.

72. Frolov, A.V. A Singularity Problem with f(R) Dark Energy. Phys. Rev. Lett. 2008, 101, 061103.

73. Guo, J.-Q.; Joshi, P.S. Spherical collapse for the Starobinsky R2 model. 2015. in preparation.

74. Hong, S.E.; Hwang, D.-I.; Stewart, E.D.; Yeom, D.-H. The causal structure of dynamical charged black holes. Class. Quantum Gravity 2010, 27, 045014.

75. Kawai, H.; Yokokura, Y. Phenomenological Description of the Interior of the Schwarzschild Black Hole. Int. J. Mod. Phys. A 2015, 30, 1550091.

76. Belinskii, V.A.; Kalathnikov, I.M.; Lifshitz, E.M. Oscillatory Approach to a Singular Point in the Relativistic Cosmology. Adv. Phys. 1970, 19, 525–573.

77. Belinskii, V.A.; Khalatnikov, I.M. Effect of scalar and vector fields on the nature of the cosmological singularity. Zh. Eksp. Teor. Fiz. 1972, 63, 1121–1134.

78. Belinskii, V.A. On the cosmological singularity. Int. J. Mod. Phys. D 2014, 23, 1430016.

79. BKL conjecture. Available online: http://www.einstein-online.info/dictionary/bkl-conjecture (accessed on 29 August 2015).

80. Carroll, S. Spacetime and Geometry: An Introduction to General Relativity; Addison-Wesley: San Francisco, U.S.A., 2003.

81. Landau, L.D.; Lifshitz, E.M. The Classical Theory of Fields; Pergamon Press: Oxford, UK, 1971; Volume 2.