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

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\"om geometries. In simulating scalar collapses in Schwarzschild and Reissner-Nordstr\"om 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\"om black hole is examined. Approximate analytic solutions for different types of collapses are partially obtained. The scalar degree of freedom $\phi$, 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\phi/dt$ is negative (positive), the physical scalar field is suppressed (magnified) by $\phi$, 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^2$ 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.


I. 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 scalartensor theory, high-dimensional theory, and f (R) gravity, etc. For a review of modified gravity theories, see Ref. [7]. For reviews of f (R) theory, see Refs. [8][9][10][11].
Static and spherically symmetric black hole solutions in f (R) gravity were explored in Ref. [12]. Charged Born-Infeld black holes for f (R) theories were studied in Ref. [13]. Instabilities and (anti)-evaporation of Schwarzschild-de Sitter and Reissner-Nordström black holes in modified gravity were discussed in Refs. [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 Ref. [17]. Spherical collapses of a charged scalar field in dilaton gravity and f (R) gravity were explored in Refs. [18] and [19], respectively. Spherical scalar collapse in f (R) gravity was simulated in Ref. [20]. Asymptotic analysis was implemented in the vicinity of the singularity of a formed black hole.

A. 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 Refs. [22,23]. For reviews on the Cauchy problem in general relativity and strong cosmic censorship, see Refs. [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 Ref. [28].
In Refs. [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 Refs. [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 Refs. [35,36], with regular initial data, spherical collapse of a charged scalar field was simulated. An apparent horizon was formed. A null, weak massinflation singularity along the Cauchy horizon and a final, spacelike, central singularity were obtained. Spherical collapses in Brans-Dicke theory, dilaton gravity, and f (R) gravity were investigated in Refs. [17][18][19]. Mass inflation phenomena were also reported.
It is important to connect approximate analytic candidate expressions with numerical results. In Refs. [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: (i) Neutral collapse: neutral scalar collapse toward a black hole formation.
(iii) 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.

B. New results
In this paper, we explore neutral scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in 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) gravity in the Jordan frame was defined in Ref. [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) singularity problem.
We explore scalar collapses in both general relativity and f (R) gravity. For convenience, we transform the dynamical system in 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/2 ln 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) gravity, the physical scalar field ψ is suppressed (magnified) when 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) gravity in this paper.
According to the strength of the scalar field, charge scattering can be classified into five types as follows: (i) 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.
(ii) Type II: null scattering. When the scalar field is intermediate, the inner horizon can contract to a place close to the center or reach the center. For each variable (the metric elements and physical scalar field), the spatial and temporal derivatives are almost equal. In the case of the center being reached, the central singularity is null. This type has two stages: early/slow and late/fast. In the early stage, the inner horizon contracts slowly, and the scalar field also varies slowly. In the late stage, the inner horizon contracts quickly, and the dynamics is similar to that in the spacelike scattering case.
(iii) Type III: critical scattering. This case is on the edge between the above two cases. The central singularity becomes null.
(iv) Type IV: weak scattering. When the scalar field is very weak, the inner horizon contracts but not much. Then the central singularity remains timelike.
(v) 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. 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 Ref. [19] where spherical charged scalar collapse in f (R) gravity was simulated, a singularity problem was reported. When a dark energy f (R) model is used, f can be pushed to 1 easily. Correspondingly, the Ricci scalar R becomes singular. This singularity problem can be avoided when an R 2 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/2lnf / √ 8πG) with 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 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 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 R 2 term to the Lagrangian for gravity [42,43]. When this R 2 term is added to the f (R) function, a singular R is pushed to regions where 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) singularity problem is avoided in the combined model (a combination of a dark energy f (R) model and the R 2 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 can 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 Sec. II, we build the framework for charge scattering, including action for charge scattering, the coordinate system, and the f (R) model. In Sec. III, we set up the numerical formalism for charge scattering. In Secs. IV, V, and VI, scalar collapses in flat, Schwarzschild, and Reissner-Nordström geomet-ries will be explored, respectively. In Sec. VII, we consider weak charge scattering. In Sec. VIII, we discuss the causes and avoidance of the singularity problem. In Sec. IX, the results will be summarized.

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

A. Action
The action for charge scattering in f (R) gravity can be written as follows: with f (R)/(16πG), L ψ , and L F are the Lagrange densities for f (R) gravity, a physical scalar field ψ, and the electric field for a Reissner-Nordström black hole, respectively. f (R) is a certain function of the Ricci scalar R, and G is the Newtonian gravitational constant. 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 (4) The electric field of a Reissner-Nordström black hole can be treated as a static electric field of a point charge of strength q sitting at the origin r = 0. In the Reissner-Nordström metric, the only nonvanishing components of F µν are F tr = −F tr = −q/r 2 . The corresponding energymomentum tensor for the electric field is [50] Although Eq. (5) is obtained in the Reissner-Nordström metric, it is valid in any coordinate system, since as seen by static observers, the electromagnetic field should be purely electric and radial [27,50].
The equivalent of the Einstein equation in f (R) gravity reads where f ≡df /dR and ≡ ∇ α ∇ α . The trace of Eq. (6) describes the dynamics of f , where T is the trace of the stress-energy tensor T µν . Defining a new variable χ by and a potential U (χ) by one can rewrite Eq. (7) as The field equations for f (R) gravity (6) are somewhat different from the more familiar corresponding ones in general relativity. Therefore, for convenience, we transform 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 one obtains the corresponding action of f (R) gravity in the Einstein frame [9] where κ = √ 8πG,g µν = χ · g µν , V (φ) ≡ (χR − f )/(2κ 2 χ 2 ), and a tilde denotes that the quantities are in the Einstein frame. The Einstein field equations arẽ wherẽ is the ordinary energy-momentum tensor for the physical matter fields in terms of g µν in the Jordan frame. With the expression for the energy-momentum tensor for the scalar field ψ in the Jordan frame, shown in Eq. (4), the corresponding expression in the Einstein frame can be written as which gives 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 can be expressed as [see Eq. (5)] Then we have in the Einstein frame, where r JF and r EF 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 r EF . We denote the total energy-momentum tensor for the source fields as The equations of motion for φ and ψ can be derived from the Lagrange equations as Alternatively, Eqs. (21) and (22) can be obtained from the corresponding ones in the Jordan frame. Some details are given in the Appendix. In the Einstein frame, the potential for φ can be written as Then we have

C. Coordinate system
In the studies of mass inflation, the double-null coordinates described by Eq. (25) are usually used, 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 Eq. (26), obtained by defining u = (t − x)/2 and v = (t + x)/2 [51], This set of coordinates is illustrated in Fig. 1. Similar to the Schwarzschild metric, the Reissner-Nordström metric can be expressed in Kruskal-like coordinates [52] (also see Refs. [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 , and the potential U (f ).
In the form of Ref. [53], the Reissner-Nordström metric in Kruskal-like coordinates in the region of r > r − can be written as where r ± (= m ± m 2 − q 2 ) and k ± [= (r ± − r ∓ ]/(2r 2 ± )] are the locations and surface gravities for the outer and inner horizons of a Reissner-Nordström black hole, respectively. r(t, x) is defined implicitly below [53], In this set of coordinates, as implied by Eq. (28), the exact inner horizon is at regions where uv and (t 2 − x 2 ) are infinite. However, it is found that, even when uv and (t 2 − x 2 ) take moderate values, r still can be very close to the inner horizon, e.g., r = (1 + 10 −10 )r − . (See Fig. 2.) Therefore, at such regions, the interaction between the scalar fields and the inner horizon still can be very strong, then we can investigate mass inflation numerically. This formalism has several advantages as follows: (i) In the line element (26), one coordinate is timelike and the rest are spacelike. This is a conventional setup. It is more convenient and more intuitive to use this set of coordinates. For the set of coordinates described by Eq. (25), in the equations of motion, many terms are mixed derivatives of u and v; while for the set of coordinates described by Eq. (26), in the equations of motion, spatial and temporal derivatives are usually separated.
(ii) 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.
(iii) The interactions between scalar fields and the geometry are local effects. In Refs. [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.
For a viable dark energy f (R) model, f has to be positive to avoid ghosts [54], and 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) 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) model, the Hu-Sawicki model, as an example. This model reads [39] f where n is a positive parameter, D 1 and D 2 are dimensionless parameters, R 0 = 8πρ 0 /3, andρ 0 is the average matter density of the current Universe. We consider one of the simplest versions of this model, i.e., n = 1, where D is a dimensionless parameter. In this model, (32) As implied in Eq. (32), to make sure that the de Sitter curvature, for which V (φ) = 0, has a positive value, the parameter D needs to be greater than 1. In this paper, we set D to 1.2 and set R 0 to 10 −5 or 10 −6 . Then, together with Eqs. (23) and (32), these values imply that the radius of the de Sitter horizon is about 1/R 0 ∼

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

A. Field equations
In double-null coordinates (26), using one obtains the equation of motion for r, where r ,t ≡ dr/dt and other quantities are defined analogously. For simplicity, we define η ≡ r 2 and integrate the equation of motion for η, instead [51]. The equation of motion for η can be obtained by rewriting Eq. (33) as In double-null coordinates, the equations of motion for φ (21) and ψ (22) become, respectively, wherẽ The {uu} and {vv} components of the Einstein equations yield the constraint equations

B. Initial conditions
We set the initial data to be time symmetric: Therefore, in this configuration, the constraint equation (41) is satisfied identically. Note that, in this configuration, the values of r ,t and σ ,t at t = 0 are the same as those in the Reissner-Nordström metric case. We set the initial value for ψ as In this paper, we give φ two sets of initial conditions: where φ 0 is de Sitter value, defined by V (φ 0 ) = 0. The initial value for σ is defined to be the same as the corresponding one in the Reissner-Nordström case (27), where r is defined by Eq. (28) with t = 0. We obtain the initial value for r in charge scattering by combining Eqs. (33) and (42), We set r ,x = σ ,x = 0 at the origin (x = 0, t = 0) as in the Reissner-Nordström metric case. The definition domain for the spatial coordinate Then r(x, t)| t=0 can be obtained by integrating Eq. (48) via the fourth-order Runge-Kutta method from x = 0 to x = ±x b , respectively. The initial values of r, σ, f , and ψ are shown in Fig. 11.
In this paper, we employ the finite difference method. The leapfrog integration scheme is implemented, which is a three-level scheme and requires initial data on two different time levels. With the initial data at t = 0, we compute the data at t = ∆t with a second-order Taylor series expansion as done in Ref. [64]. Take the variable ψ as an example, The values of ψ| t=0 and ψ ,t | t=0 are set up as discussed above, and the value of ψ ,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 firstorder time derivatives of r, σ, φ, and ψ at t = 0 described by Eq. (43) ensure that the constraint equation (41) is satisfied. The equation for r ,xx at t = 0 expressed by (48) implies that the constraint equation (42) is satisfied. Computations of r, σ, φ, and ψ at t = ∆t via a secondorder Taylor series expansion, as expressed by Eq. (49) for the case of ψ, satisfy all the equations of motion.
Note that the region of x < 0 is included in the initial conditions. This may not be physical. However, we are mainly interested in the interior dynamics of black holes, and then it is not important where the scalar field originally comes from. This setup makes it convenient for us to compare the results of charge scattering to the known solutions of the Reissner-Nordström geometry. Therefore, we use this setup as a toy model.

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

D. 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.
The equations of motion for φ (36) and ψ (37) are coupled. Newton's iteration method is employed to solve this problem [64]. With the illustration of Fig. 4, the initial conditions provide the data at the levels of "down" and "here". We need to obtain the data on the level of "up". We take the values at the level of "here" to be the initial guess for the level of "up". Then, taking ψ as an example, we update the values at the level of "up" using the following iteration: where G(ψ up ) is the residual of the differential equation for the function ψ up , and J(ψ up ) is the Jacobian defined by We do the iterations for the coupled equations one by one, and run the iteration loops until the desired accuracies are achieved.

E. Tests of numerical code
To make sure that the numerical results are trustworthy, one needs to test the numerical code. We compare the numerical results obtained by the code with the analytic ones for the dynamics in a Reissner-Nordström geometry, and examine the convergence of the constraint equations and dynamical equations in charge scattering.
The dynamics in a Reissner-Nordström geometry is a special case for charge scattering, in which the contributions from the scalar fields are set to zero. This special case has analytic solutions expressed by Eqs. (27) and (28). Therefore, we can test our code by comparing the numerical and analytic results in the Reissner-Nordström geometry. Set m = 1, q = 0.7, and ∆x = ∆t = 10 −4 . We plot the evolutions of r and σ along the slice (x = 3×10 −4 , t = t) in Fig. 5(a). As shown in Fig. 5(a), numerical and analytic results match well at an early stage; while at a late stage where gravity and electric field become strong, the numerical evolutions have a time delay compared to the analytic solutions.
When the numerical results are obtained, we substitute the numerical results into the discretized equations of motion and the constraint equations, and find that they are well satisfied. See Figs. 13 and 15, for example. Moreover, the convergence of the constraint equations (41) and (42) is examined. We assume one constraint equation is nth-order convergent: residual=O(h n ), 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, 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} constraint equation (41) in Fig. 5(b) for the slice (x = x, t = 0.65). 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 where F h is the numerical solution. Then, for step sizes equal to h/2 and h/4, we have Defining The convergence tests for η ≡ r 2 , σ, φ, and ψ are investigated. They are all second-order convergent. As a representative, the results for ψ are plotted in Fig. 5(b) for the slice (x = x, t = 0.65). The values of the parameters in charge scattering in this section are described at the beginning of Sec. VI. We use the spatial range of x ∈ [−10 10] and the grid spacings of h = ∆x = ∆t = 0.02.

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

A. Numerical setup
The numerical setup in neutral scalar collapse in f (R) gravity is discussed in Ref. [20]. The dynamical equations for r, η, σ, φ, and ψ can be obtained by setting the terms related to the electric field in the corresponding equations in Sec. III A to zero: whereT (ψ) = e 2σ (ψ 2 ,t − ψ 2 ,x )/χ. In the equation of motion for σ (54), the term (r ,tt − r ,xx )/r can create big errors near the center x = r = 0. To avoid such a problem, we use the constraint equation (39) alternatively [51]. Defining a new variable g one can rewrite Eq. (39) as the equation of motion for g, 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 Eq. (57).
We set the initial data as The initial values for χ[≡ exp( 2/3κφ)] and ψ(r) are defined as with a = 0.2, b = 0.1, r 1 = r 2 = 4, and U (χ 0 ) = 0. The parameters for the Hu-Sawicki model (30) are set as D = 1.2 and R 0 = 10 −6 . The local Misner-Sharp mass m is defined as [67] (See Ref. [68] for details on various properties of the Misner-Sharp mass/energy in spherical symmetry.) Then on the initial slice (x = x, t = 0), the equations for r, m, and g are [20,51] Set r = m = g = 0 at the origin (x = 0, t = 0). Then the values of r, m, and g on the initial slice (x = x, t = 0) can be obtained by integrating Eqs. (63)-(65) from the center x = 0 to the outer boundary x = x b via the fourth-order Runge-Kutta method. The values of r, σ, φ, and ψ at t = ∆t can be obtained with a second-order Taylor series expansion, as discussed in Sec. III B. The value of g at t = ∆t can be obtained using Eq. (57). The range for the spatial coordinate is defined to be x ∈ [0 20]. At the inner boundary x = 0, r is always set to zero. The terms 2(−r ,t φ ,t + r ,x φ ,x )/r in Eq. (55) and 2(−r ,t ψ ,t + r ,x ψ ,x )/r in Eq. (56) need to be regular at x = r = 0. Since r is always set to zero at the center, so is r ,t . Then we enforce φ and ψ to satisfy φ ,x = ψ ,x = 0 at x = 0. The value of g at x = 0 is obtained via extrapolation. We set up the outer boundary conditions at x = 20 via extrapolation. The temporal and the spatial grid spacings are ∆t = ∆x = 0.005.
The numerical code is second-order convergent and is the one developed in Ref. [20].

B. Black hole formation
The evolutions of r, σ, φ, and ψ are plotted in Figs. 6(a)-6(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 ,ν = e 2σ (−r 2 ,t + r 2 , Using this property, we locate the apparent horizon and plot it in Figs. 6(e) and 6(f). As shown in Fig. 6(a), the central singularity is also approached in the collapse. Therefore, a black hole is formed.
C. 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.

Then we have
Using Eqs. (11) and (73), there is So as shown in Fig. 6(c), when C is positive, χ also approaches zero as ξ and r approach zero. Then the transformation between the Jordan and Einstein frames, g µν , 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.

D. Mass inflation
In the vicinity of the singularity, in the equation of motion for σ (68), because of the contribution from φ, σ(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 Eqs. (62), (71)-(74), the mass function can be written as where K ≡ |r ,x /r ,t |. In the Schwarzschild black hole case, 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. However, the parameter C is not zero. Then the metric quantity σ is modified. [See Eq. (72).] As a result, the delicate balance between r and e 2σ (−r 2 ,t + r 2 ,x ) is broken. Consequently, as implied in Eq. (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 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 Refs. [20,71], and plot two sample sets of results on the slices (x = 1, t = t) and Since 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) gravity in the Jordan frame was defined in Ref. [40]. However, due to the complexity of some integrals, an explicit quasilocal 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 µν , there are e 2σ | JF = χ·e 2σ | EF and r JF = χ −1/2 ·r EF . Considering that near the central singularity [20,41] one can obtain that K JF ≈K EF . Then the mass function in the Jordan frame can be written as In the case of C > 0, due to the factor r − √ 2/3κC , m JF is greater than m EF . The m JF along the slice (x = 1, t = t) is plotted in Fig. 7(d).
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.

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

A. A dark energy f (R) singularity problem
In neutral scattering, for usual initial conditions of the scalar degree of freedom f , f asymptotes to zero as the central singularity is approached, which is similar to what happens in the neutral scalar collapse discussed in Sec. IV. Details are skipped here. On the other hand, when the initial velocity or acceleration of f is large enough, f can become 1 before the central singularity is approached. As implied in Eq. (29), for dark energy 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.
The parameters are set as follows: (i) Schwarzschild geometry: m = 1.
The results in this circumstance are plotted in Fig. 9. We plot the dynamics of φ on the slice (x = −2, t = t) in Fig. 9(d), from which one can see that φ is mainly accelerated by the spatial derivative φ ,xx and the geometrical term −2r ,t φ ,t /r. Eventually, φ and f approach 0 and 1, respectively. Then the Ricci scalar becomes singularity, and the simulation stops, as shown in Fig. 9(a).

B. 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 R 2 term to the dark energy f (R) model [42][43][44][45][46][47][48][49]. In this paper, we add the R 2 term to the Hu-Sawicki model, In this combined model, at high curvature scale, f ≈ 2αR. So a singular R is pushed to regions where 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) model (79) is set to 1. The numerical results for neutral scattering for this modified model are plotted in Fig. 10. As shown in Fig. 10(b), for this model, f can cross 1 without difficulty. What are the limits that f and R can reach? As shown in Fig. 10(d), as the central singularity is approached, there is Since r ,t is negative, the above equation describes a positive feedback system of φ ,tt and φ ,t : |φ ,t | produces more of |φ ,tt |, and in turn |φ ,tt | produces more of |φ ,t |. Since φ ,t is positive, φ can be rapidly accelerated to positive infinity. Then 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 + αR 2 , 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 R 2 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.

VI. 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: (i) Reissner-Nordström geometry: m = 1, and q = 0.7.
A. Evolutions

Outline
In this subsection, we describe the evolutions of r, σ, f , and ψ that are plotted in Fig. 11. Examining the equations of motion (33)-(37) and the numerical results  (i) Metric components: r and σ. They contribute as gravity.
(iii) Electric field and V (φ). As implied in Eq. (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 (φ) is negligible.
Furthermore, these quantities can be separated into two sides: the gravitating side (r, σ, φ, and ψ) and the repulsive side [electric field and 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: (i) 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) in Fig. 11. See Sec. VI B.
(ii) 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) in Fig. 11. See Secs. VI C and VI D.
(iii) 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: Fig. 11. Due to the similarity to that in general relativity discussed in Ref. [41], details on this type of scattering in f (R) gravity are skipped in this paper.
(iv) Type IV: weak scattering. When the scalar field is very weak, the inner horizon does not contract much. Sample case: Fig. 19. See Sec. VII.
(v) 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) in Fig. 11.
In this paper, we will discuss Types I, II, and IV.

Causes of mass inflation and evolutions
The local Misner-Sharp mass for a charge black hole hole is In a Reissner-Nordström geometry, in Kruskal-like coordinates expressed by Eq. (27), near the inner horizon, although σ asymptotes to positive infinity, (r 2 ,t − r 2 ,x ) is much less than e −2σ . Consequently, e 2σ (r 2 ,t − r 2 ,x ) approaches zero. Then as implied in Eq. (81), the mass function m takes a finite value and is equal to the black hole mass. For more details, see Ref. [41].
In charge scattering, the equations of motion for η≡r 2 When the scalar field is weak enough (e.g., −4 < x < −1), the inner horizon does not change much. At the intermediate state (e.g., 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, are not that accurate. We are aware that the inner horizon is actually at infinity, while r still can be very close to r− when x and t take moderate values. and σ are At the beginning, |φ ,t | and |ψ ,t | may be less than |φ ,x | and |ψ ,x |, respectively. However, as r decreases toward the central singularity, gravity becomes stronger. Then |φ ,t | and |ψ ,t | become greater than |φ ,x | and |ψ ,x |, respectively. [See Fig. 15.] As a result, in this case, the repulsive "force", 4π[φ 2 ,t − φ 2 ,x + (ψ 2 ,t − ψ 2 ,x )/χ] + e −2σ q 2 /r 4 , is greater than the corresponding one, e −2σ q 2 /r 4 , in a Reissner-Nordström geometry. This makes σ accelerate faster than in the Reissner-Nordström geometry. Consequently, the repulsive force from 2e −2σ (q 2 /r 2 − 1) for η≡r 2 is much weaker than the corresponding value in the Reissner-Nordström geometry. As a result, near the inner horizon, |r ,t | is much greater than the corresponding one in the Reissner-Nordström case. Then (r 2 ,t − r 2 ,x ) moves from extremely tiny values in the Reissner-Nordström metric case to moderate values, and r crosses the inner horizon r = r − for the given Reissner-Nordström geometry. With Eq. (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 , and ψ are plotted in Figs. 11(a)-11(d).
Because of gravity from the black hole and the contribution from the physical scalar field, for some configurations, f can decrease to zero as the central singularity is approached. [See Eq. (36) and Figs. 11,13,and 15.] Nothing is wrong with this circumstance. However, for certain configurations, f can be pushed to 1, as shown in Fig. 21. Correspondingly, the Ricci scalar R becomes singular. We will discuss the former case in this section and the latter case in Sec. VIII.
The evolution of ψ is plotted in Fig. 11(d). 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.

Locations of horizons
We locate the outer and inner horizons using the following equation, The results are plotted in Figs. 11(e) and 11(f). Due to the absorptions of the energies of φ and ψ, the outer horizon increases from the original value of 1.7 to 3.7. Note that the results for the inner horizon, especially at regions where 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 − even when x and t take moderate values.

B. 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) as an example, we plot the evolutions of r, σ, f , ψ, and m on this slice in Fig. 12. The mass function remains equal to the mass of the original Reissner-Nordström black hole, m 0 = 1, until r is very close to 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 Fig. 13.
The strongness of the scalar field causes several consequences as below. Since φ ,t is negative, the term 2/3κφ ,t ψ ,t in Eq. (37) functions as a friction force for ψ. Consequently, compared to φ, ψ grows slowly. As the singularity is approached, ψ even approaches constant values. [See Figs. 14(c) and 14(d).] As shown in Fig. 13(e), near the singularity, two major sets of terms can be expressed as Alternatively, (iv) Growth of mass function. In spacelike scattering, the equation of motion for φ can be simplified as Consequently, with the results obtained in Sec. IV C, σ has the following asymptotic solution: with φ ≈ C ln ξ. Then similar to the mass function (76) in neutral collapse, with Eq. (81), the mass function in spacelike scattering can be written as Numerical results show that, for the sample slice (x = 1.52, t = t), near the central singularity, the slope of the singularity curve K is about 0.04. As shown in Fig.14(d), we linearly fit the numerical results of m via One can see that the above three sets of results match well.

C. 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 .
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     Figure 13: (color online). Dynamics along the slice (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 2 ,t . (b) η,tt ≈ η 2 ,xx . (c) σ,tt ≈ 4πφ 2 ,t . (d) φ,tt ≈ −2r,tφ,t/r. (e) r,tψ,t≈r,xψ,x and φ,tψ,t ≈ φ,xψ,x. Then we have ψ,t/ψ,x≈r,x/r,t ≈ φ,x/φ,t.
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. Take the slice (x = 0.5, t = t) as a sample slice, we plot the terms in the field equations for r, σ, φ, and ψ in Fig. 15, and the evolutions of r, σ, φ, ψ, m, and |1 − K 2 | in Fig. 16. We investigate the dynamics in the vicinity of the central singularity via mesh refinement and plot the results in Fig. 17. The null scattering has two stages: early/slow and late/fast. As shown in Figs. 15 and 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 Fig. 15, when r is very small, the equations of motion for r (33), σ (35), and φ (36) can be rewritten as 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 Fig. 17. Then we have As shown in Fig. 16(f       (36), and ψ (37) are reduced as follows: r ,tt ≈r ,xx , r 2 ,t ≈r 2 ,x ; σ ,tt ≈ σ ,xx , ψ 2 ,t ≈ψ 2 ,x , r ,tt ≈r ,xx , φ 2 ,t ≈φ 2 ,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) or (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 Fig. 18, the constraint equations (41) and (42) provide some useful information on the connections between some variables at the early stage of charge scattering: As shown in Fig. 16, the quantities r, σ, f , ψ, |1−K 2 |, and m change dramatically at the beginning of charge scattering where r≈r − . Note that, near the central singularity, r, σ, and f have approximate analytic expressions in terms of ξ = t 0 −t, where t 0 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 − t s , where t s 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 Fig. 18, from which one can see that r, r ,t , and σ may have the following approximate analytic expressions: σ≈b ln ζ + σ 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 Eqs. (28) and (47), one obtains that, as r approaches the inner horizon r = r − , in the case of t x, σ can be approximated by a logarithmic function of t. As shown in Fig. 15(c), describing the terms in the equation of motion for σ in charge scattering, at the very early stage (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 Figs. 18(c) and 18(d), φ and ψ can be well fitted by power law functions of ζ. Take φ as an example, We plot (1 − K 2 ) in Fig. 18(d) and find that ln(1 − K 2 ) can be well fitted linearly with respect to ζ, Currently we do not have derivations for this linear relation. The good thing is that, in the mass function, (1 − K 2 ) is a minor factor. Therefore, as verified in Fig. 18(d), the mass function can be reduced to where a, b, λ, and σ 0 are defined in Eqs. (103) and (104). We list the fitting results below: Terms in {tt + xx} const. eq.

VII. 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: (i) Physical scalar field: ψ(x, t)| t=0 = a · exp −(x − x 0 ) 2 /b , a = 0.03, b = 1, and x 0 = 4.
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 Fig. 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 Ref. [23] and the numerical work in Ref. [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 Ref. [74], it was argued that, when Hawking radiation is taken into account, this censorship may also be violated. Note that in Ref. [75] the interior of a Schwarzschild black hole was also discussed with the backreaction from the Hawking radiation being taken into account.
The dynamics for the quantities r, η, σ, φ, and ψ are plotted in Figs. 20(a)-20(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 . Moreover, the derivatives have oscillations. As shown in Fig. 20(f), the mass function keeps growing even as r approaches a constant value. Further details are skipped.

VIII. DARK ENERGY 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 asymptotes to zero as the center is approached. Note that, as discussed in Sec. V, in neutral scattering in dark energy f (R) gravity, when the initial velocity or acceleration of f is large enough, 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 Sec. VI, including grid spacings ∆x = ∆t = 0.005. However, the initial value for φ takes the following format:      φ(x, t)| t=0 = a exp[−(x − x 0 ) 2 ] + φ 0 , with V (φ 0 ) = 0, a = −0.05, and x 0 = 0.3. We plot the numerical results in Fig. 21. As shown in Fig. 21(b), as the inner horizon is approached, 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 Figs. 21(d) and 21(e), we plot the terms in the dynamical equation for φ on the slice (x = −2, t = t), and find that initially because of the contribution from φ ,xx , φ ,t changes from φ ,t | t=0 = 0 to   In (e), as the inner horizon is approached, φ,tt ≈ −2r,tφ,t. This equation describes a positive feedback, since −2r,t/r is positive. As a result, when φ,t is positive, φ can be accelerated to zero rapidly. Correspondingly, 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). φ ,t > 0 at late time. Near the inner horizon, there is Then φ is accelerated by gravity. After the inner horizon is met, there is φ ,tt ≈ φ ,xx > 0. Eventually, φ, 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 R 2 term to the Hu-Sawicki model, The parameters take the same values as in the last subsection with the following exceptions: (i) Grid spacings: ∆x = ∆t = 0.0025.
The numerical results for charge scattering for this modified model are plotted in Fig. 22. As shown in Fig. 22(b), for this model, f can cross 1 without difficulty. The simulation can run smoothly.

IX. SUMMARY
In this paper, we studied scalar collapses in flat, Schwarzschild, and Reissner-Nordström geometries in f (R) gravity numerically. Approximate analytic solutions for different types of collapses were partially obtained. One dark energy f (R) singularity problem was discussed. We summarize our work on computational and physical issues separately below. equilibrium state and is static. In f (R) gravity, inside compact stars, the scalar degree of freedom, 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 and the matter fields are decoupled. As a result, f is more free to move than in the compact stars case. It can keep increasing or decreasing until singularities are met.
(iv) Dark energy f (R) singularity problem: cosmology (or static compact objects) vs black hole physics.
In dark energy 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, the equation of motion for f is where H is the Hubble parameter. Due to the finiteness of the potential barrier, the force from a perturbation of 8πT /3 may push 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/2lnf / √ 8πG) has more complex structure (36), − φ ,tt + φ ,xx + 2 r (−r ,t φ ,t + r ,x φ ,x ) = e −2σ V (φ) + 1 √ 6 κ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 (111), and gravity from the black hole, −2r ,t φ ,t /r, can cause a similar singularity problem as in cosmology or static compact stars.
(v) The combined f (R) model and the R 2 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 R 2 model, f (R) = R + αR 2 . This model has a non-singular de Sitter solution, which is unstable both to the past and to the future [42,43]. In Sec. V, scalar collapse in a Schwarzschild geometry for the combined model (a combination of dark energy model and R 2 model) was explored. A new Schwarzschild black hole, including a new central singularity, can be formed. Moreover, under certain initial conditions, f and R can be pushed to infinity as the central singularity is approached. In Ref. [73], scalar collapse in flat geometry for the R 2 model was simulated. Similar results were obtained. Namely, the classical singularity problem, which is present in general relativity, remains in collapse in these models.
(vi) 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) 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) singularity problem in collapse were discussed.
In this appendix, based on the equations of motion for a massless physical scalar field ψ and the scalar degree of freedom f in the Jordan frame for f (R) gravity, we derive the corresponding equations in the Einstein frame. In the transformation from the Jordan frame into the Einstein frame, we useg µν = χ · g µν and κφ ≡ 3/2 ln χ, where a tilde denotes that the quantity is in the Einstein frame.