Strong cosmic censorship: The nonlinear story

A satisfactory formulation of the laws of physics entails that the future evolution of a physical system should be determined from appropriate initial conditions. The existence of Cauchy horizons in solutions of the Einstein field equations is therefore problematic, and expected to be an unstable artifact of General Relativity. This is asserted by the Strong Cosmic Censorship Conjecture, which was recently put into question by an analysis of the linearized equations in the exterior of charged black holes in an expanding universe. Here, we numerically evolve the nonlinear Einstein-Maxwell-scalar field equations with a positive cosmological constant, under spherical symmetry, and provide strong evidence that mass inflation indeed does not occur in the near extremal regime. This shows that nonlinear effects might not suffice to save the Strong Cosmic Censorship Conjecture.


I. INTRODUCTION
The Strong Cosmic Censorship (SCC) conjecture embodies the expectation that General Relativity is a deterministic theory. It does so by predicting that Cauchy horizons (CH) -the boundaries of the maximal evolution of initial data via the Einstein field equations -are unstable, and give rise, upon perturbation, to singular boundaries, beyond which the field equations cease to make sense.
Quite surprisingly, recent results [1][2][3][4][5][6][7][8] put into question the validity of SCC in the context of highly charged black holes (BHs) immersed in a spacetime with a positive cosmological constant 1 . More precisely, these results provide evidence for the existence of stable (charged de Sitter) BH configurations containing a CH in their interior, along which the Misner-Sharp mass, a scalar invariant measuring the energy content of symmetry spheres, remains bounded -a no mass inflation scenario. In particular, the CH should then retain enough regularity to allow evolving the spacetime metric across it using the Einstein equations. However, such evolution is not unique, condemning determinism (and the SCC conjecture) to failure. This is particularly disturbing in view of the fact that a positive cosmological constant provides the standard mechanism to model the observed accelerated expansion of our universe. Nonetheless, the results mentioned above are restricted to either the linear setting or to the nonlin- 1 These results appeared almost two decades after the pioneering work in [9,10], which made contradicting claims; see section 2.2 in [6] for a clarification of the implications of these papers regarding SCC.
ear analysis of the geometry of the BH interior starting from an already completely formed event horizon -i.e. the corresponding (horizon) data are put in "by hand". These results provide clear expectations concerning the stability/instability of Cauchy horizons within de Sitter BHs, but are not enough to cast a final verdict on SCC for the following reasons: first, the parameter ranges identified in [4] as potentially problematic for SCC are very narrow, and therefore even small non-linear deviation form these might be enough to save SCC; secondly the (non-linear) results in [1] do not take into account the oscillatory behavior of the scalar field along the event horizon 2 which was identified in [4] for the first time. Moreover, our numerical simulations will allow us to gain further information concerning the behavior of various relevant quantities near the Cauchy horizon, as for instance a quantitative understanding of tidal deformations. To go beyond the previous studies we perform a full nonlinear numerical evolution of both massless and massive minimally coupled self-gravitating scalar fields, in a spacetime with a Maxwell field and a positive cosmological constant. For technical reasons we will restrict ourselves to the spherically symmetric setting, but we should stress that, according to the results in [4], the spherically symmetric mode plays a key role in the linear stability/instability of Cauchy horizons for near extremal Reissner-Nordström de Sitter BHs.
All the nonlinear simulations we will be discussing evolve from characteristic initial data whose outgoing component is located in the BH exterior; in particular, the data along the event horizon arises dynamically in this framework. Our numerical code also allows us to probe the BH interior region and to derive a detailed description of the behavior of fundamental quantities, such as the radius function, scalar field, (Misner-Sharp) mass and curvature. By examining the vicinity of the BH parameters identified as potentially problematic for SCC in [4], we find stable no mass inflation scenarios arising from a full nonlinear evolution (of exterior data). These are, to the best of our knowledge, the first results of this kind. They show, in particular, that nonlinear effects are apparently not strong enough to save SCC in the context of highly charged de Sitter BHs 3 .
Our setup also allows the inclusion of a scalar field mass, and so we take the opportunity to investigate the possibility, raised in [5] and [6], that the curvature might be bounded up to the Cauchy horizon for certain choices of this mass. This scenario is disproved in the cases that we analyze.

II. SETUP
We consider here an evolving, electrically charged spacetime, modelled by the Einstein-Maxwell action with a cosmological constant Λ, minimally coupled to a massive scalar field Φ with mass parameter µ, where F 2 = F αβ F αβ and F αβ is the Maxwell tensor. The equations of motion reduce to where is the Hodge dual. We focus on spherically symmetric spacetimes, written in double null coordinates as where u and v are ingoing and outgoing coordinates, respectively. In this framework, Maxwell's equations decouple and imply that with Q a conserved (electric) charge.

III. NUMERICAL EVOLUTIONS
To numerically evolve the field equations we specify initial conditions along two null segments, u = u i and v = v i . We fix the residual gauge freedom as follows: where r u0 is a constant and r 0 = v i . The profile of the scalar field is set as purely ingoing with the outgoing flux being set to zero, Φ ,u (u, v i ) = 0. See the Supplemental Material for more information about the integration procedure.
To interpret our results it will be convenient to consider the following alternative outgoing null coordinates: • v, an Eddington-Finkelstein type coordinate, defined by integrating along the event horizon (EH), and t, the affine parameter of an outgoing null geodesic, obtained by integrating along a constant u line. In these expressions M stands for the Misner-Sharp mass function, which we also closely monitor during the integration, given by The constant r u0 is thus related to the initial BH mass, Recall that the blow-up of this scalar signals the breakdown of the field equations [17] (compare with [12,13]).
To estimate the curvature we compute the Kretschmann scalar K computed from the field equations (see Supplemental Material; a direct evaluation of this scalar in terms of the metric was found to lead to important round-off error-related problems).
According to the results in Refs. [1,4], concerning the massless case, we expect the curvature to blow up for all non-trivial initial data throughout the entire subextremal parameter range. Although it is a potentially interesting nonlinear effect, we recall that the blow-up of K, per se, is of little significance: it implies neither the breakdown of the field equations [18] nor the destruction of macroscopic observers [19]. Recall that the results in Refs. [5,6] suggest that the introduction of scalar mass could lead, for appropriate choices of BH parameters, to solutions with bounded curvature. As we will see below, our results contradict this expectation.

IV. INITIAL CONDITIONS
The physical problem is then fully determined upon specifying Q, Λ, µ, M 0 , A, v c and w. Since our purpose here is to determine whether the linearized predictions of Refs. [4,5] hold in the full nonlinear regime, we focus on M 0 = 1, Λ = 0.06 and use the following configurations: A: Q = 0.9000, µ = 0, corresponding to Q = 0.890Q max . In this case, the results in Ref. [4] (lower left panel of Fig. 3) predict mass inflation. B: Q = 1.0068, µ = 0, corresponding to Q = 0.996Q max . Linearized studies provide evidence in favor of a no mass inflation scenario [4]. C: Q = 1.0068, µ = 1.0. The results of Ref. [4] together with those of Ref.
[5] -see Fig. 2 -also provide evidence in favor of a no mass inflation scenario. Here we are considering the superposition of both neutral massless scalar perturbations [4] and charged massive scalar perturbations [5] as being the most predictive of the full non-linear evolution. If we just take into account massive scalar perturbations, then the results in Ref. [5] (see Fig. 2]) and [6] (page 22) suggest that curvature might also be bounded.
To test the dependence of our results on initial data, we use the following initial profiles for the scalar field:  We have evolved the relevant system of equations using the DoNuTS code, described briefly in the Supplemental Material. It is based on the formulation presented in Refs. [20][21][22], but the integration technique makes it spectrally accurate in the v-direction and, correspondingly, runs with trivial memory requirements and orders of magnitude faster than previously reported codes. It is important to start by noticing that, as widely expected [1,11], our numerics show that all solutions contain a non-empty CH in their BH interior. This can be attested by monitoring the radius function -shown in Fig. 1 -along null lines u = u 1 , for u 1 > u EH , where u = u EH is the event horizon. In fact, for u 1 larger but close to u EH , the radius converges, in v, to a nonvanishing constant. It is also interesting to note that, for some initial configurations and large enough u 1 , the radius does converge to zero, signaling (in that region) a singularity beyond which the metric cannot be extended [23].
As is well known the behavior of the scalar field along the event horizon is of great significance for the structure of the BH interior region. The first noteworthy feature of our results (see Supplemental Material for details) is that, as expected, the field decays exponentially (inv). More surprisingly, we also clearly observe an oscillatory profile; this might seem odd at first, since it is in contrast with what happens for Λ = 0 and with the expectation created by the study of sufficiently sub-extremal BHs with 0 < M 2 Λ 1 [24]. However, it turns out that such behavior should be expected from the linearized analysis of Refs. [4,5], where it is shown that, for a configuration resembling our configuration B, there are two modes which dominate the response: a non-oscillatory "near extremal" (NE) mode with characteristic frequency ω NE ∼ −0.081i, and a "photon sphere" (PS) mode with ω PS ∼ 0.096 − 0.095i (these numbers are given in the units and time-coordinate of Ref. [4]). Here we find very good agreement with the PS mode (when translated to ourv coordinate) which is oscillatory in nature. Similar agreement can be found for the remaining configurations A and C. We also recall that, according to the results in [4], in the M 2 Λ 1 case, the dominant mode is a non-oscillatory "de Sitter" mode, in agreement with [24].
Our main results (concerning mass and curvature) are summarized in Figs. 2-4. Fig. 2 shows the evolution of the mass function and the Kretschmann scalar for configurations A: in these cases mass inflation occurs, and, consequently, the curvature invariant K diverges. Note that an observer crossing one such region will be subjected to physical deformations which are not necessarily infinite (see discussion below). Nonetheless, because there is mass inflation, the singularity is strong enough to deserve the classification of terminal boundary, since it corresponds to a locus where the field equations cease to make sense. These conclusions are consistent with the linear results in [4] and the nonlinear results in [1].
The main novel result of this work concerns Fig. 3: there exist configurations for which no mass inflation occurs. For configurations B, in the BH interior, after a small "accretion" stage the mass settles to a constant value. Moreover, as recently predicted [1,4], the CH remains a curvature singularity, since the curvature scalar K diverges. However, the lack of mass inflation makes the singularity so "mild" that, in principle, one should be able to continue the evolution of the space-time metric across it, by solving the Einstein field equations! A somewhat unexpected feature (of configuration B) is the oscillatory way in which the curvature scalar diverges. In hindsight, such behavior could be expected from the previously discussed oscillatory behavior of the scalar field along the event horizon. Note that in a no mass inflation situation it is the blow up of Φ ,v /r ,v that dominates the behavior of K. This should be contrasted with what happens when mass inflation occurs: then it is the monotone divergence of the mass that controls the Kretschmann; this last fact provides an explanation for the non-oscillatory behavior observed for configuration A.
Concerning massive scalars, the results presented in Fig. 4 identify configuration C as another no mass inflation configuration. Once again we find that the corresponding CH is a "weak" curvature singularity. In fact,  the presence of scalar mass seems to have no attenuation effect on the growth of K, in contrast with what might be expected from linear analysis [5] and [6].
We finish this section with some further remarks concerning the blow up of curvature. In configuration A, our results indicate that the Kretschmann scalar blows up as t −2 (possibly modulated by logarithmic terms), where t is the affine parameter defined in (9) with the Cauchy horizon located at t = 0. This might suggest that the curvature blows up as t −1 , but, as noted in [19,25], there are curvature components that may blow up even faster. In fact, the quantities that determine the blow-up of the Kretschmann scalar are the square of the (Misner-Sharp) mass M and the square of the gradient of Φ, which is dominated by (Φ ,v /r ,v ) 2 . However, all curvature components are controlled by M and (Φ ,v /r ,v ) 2 (the origin of the last term can be traced to the energy-momentum tensor). From the behavior of the Kretschmann scalar we can then conclude that the components of the curvature blow up at most as t −2 ; we also expect inverse logarithmic powers [12,13,19] that are hard to detect numerically. Although divergent, these curvature components should yield (with the help of the logarithmic terms) a finite "tidal deformation" when integrated twice with respect to t, in agreement with the picture in [19].
From the equation (see Ref. [26]) we conclude that no mass inflation is essentially equivalent to the integrability of (Φ ,v /r ,v ) 2 , with respect to t 4 . Moreover, when both occur we immediately see that the curvature can only give rise to finite "tidal deformations". This reasoning is verified by our results concerning configuration B, for which the mass is bounded, Φ ,v /r ,v grows slower than t −1/2 and the Kretschmann scalar and the curvature components blow up at most as t −1 .

VI. DISCUSSION
The main motivation for this study was to understand whether nonlinear effects could trigger mass inflation, even when the linearized analysis suggests otherwise [4]. We found that nonlinear effects are not strong enough to change the picture: in fact, the nonlinear results are in full agreement with the linearized predictions. The linearized analysis of Ref. [4] suggests that no mass inflation should occur for BH charge above a threshold Q * 0.95. Within a nonlinear evolution, the precise linearized results are difficult to reproduce (for instance, the final 4 In fact, there is a mathematical equivalence under the reasonable, in principle generic, assumption that the quantity 1− 2M r + Q 2 r 2 − Λ 3 r 2 does not vanish at the Cauchy horizon. spacetime parameters depend on the initial parameters and on the size of the initial data). However, for small scalar amplitudes, our results are indeed consistent with the previous threshold. The (numerical) no mass inflation solutions presented here are the first solutions of this kind arising from the full nonlinear evolution of exterior data. They contain a Cauchy horizon in their BH interior region that can be seen as ("weakly") singular, due to the divergence of curvature invariants. However, these divergent tidal forces are not necessarily strong enough to lead to a divergent tidal deformation and the consequent unequivocal destruction of all macroscopic objects [19]. Even more problematic, the lack of mass inflation indicates that these Cauchy horizons should maintain enough regularity as to allow the field equations to determine (classically), in a highly non-unique way, the evolution of the metric to their future. This corresponds to a potential severe violation of SCC.
Our results concern spherically symmetric spacetimes. The picture is unlikely to change even with asymmetric initial conditions [4]. Thus, from the conceptual point of view [5], our results show that SCC is not enforced by the field equations. In the meantime, interesting suggestions to remedy SCC, in the presence of a positive cosmological constant, have been put forward: these include enlarging the allowed set of initial data by weakening their regularity [27], or restricting the scope of SCC to the uncharged BH setting [28]. It thus seems plausible that the astrophysical interpretation of SCC remains valid, once other fields and realistic BH charges are considered.

ACKNOWLEDGMENTS
We are indebted to Kyriakos Destounis, Aron Jansen, and Peter Hintz for many and very useful conversations and comments. R. L. gratefully acknowledges the hospitality of the group at CENTRA, where this work was started. We are grateful to the Yukawa Institute for Theoretical Physics at Kyoto University, for hospitality while this work was completed during the YITP-T-18-05 on "Dynamics in Strong Gravity Universe." R. L.
and are subjected to the following constraints These equations must be satisfied by the initial data. Then, by virtue of the Bianchi identities, they will be satisfied in the whole computational domain provided that the dynamical equations are accurately satisfied.
To integrate these equations, we start by transforming them into a system of ODEs. Our procedure is as follows. Let h(u, v) be any evolved quantity r(u, v), σ(u, v) and Φ(u, v). Defining f (v) = ∂ u h(u, v), all dynamical equations, for fixed u, have the form where denotes the derivative with respect to v. These equations can be solved by introducing the integrating factor Multiplying Eq. (A6) by λ(v), we get which are ODEs in u for all values of v. Given initial conditions in the two null segments u = u i , h(u i , v) ∀v we can integrate the equations in a rectangular region u i < u < u f and v i < v < v f . For our three functions in Eqs. (A1), (A2) and (A3), p(v) and g(v) are the following: We integrate these equations using the Double Null Through Spectral methods (DoNuTS) code written in Julia [29]. To integrate the system within DoNuTS, all functions are expanded in a Chebyshev basis in the v direction (where all v derivatives and integrations can be readily performed), and the remaining ODEs in the u direction are integrated using an adaptive step integrator through the DifferentialEquations.jl Julia package [30].

Adaptive gauge
When using the initial gauge, r ,u becomes extremely large around the apparent horizon for large v. Therefore, in order to explore the near-horizon region at late times, it is convenient to use an adaptive gauge in u during the numerical evolution.
Since the change u →ũ(u) together with σ → σ − 1 2 log dũ du leaves the equations invariant, we can change the gauge in u along the integration by choosing appropriately the initial condition σ ,u (u, v i ) at each value of u.
To explore the near-horizon geometry, we can choose an Eddington-like gauge for u, i.e., a gauge that brings the event horizon to u → ∞. A good way to do so, as described in [31], is to set σ(u, v f ) = log (2r ,v (u, v f ))+C, where C can be any constant. In the DoNuTS code, this is achieved by picking the initial condition for σ ,u (u, v i ) (A7) The term 3 2 log 2 is chosen so that σ ,u (u i , v i ) is small when σ(u i , v f ) ≈ − 1 2 log 2. With this condition, σ(u, v f ) is damped towards the desired value log (2r ,v (u, v f )) − 3 2 log 2 along the evolution in u. Additionally, in order to satisfy the constraint equation (A4), we must introduce an additional ODE for the initial condition r , with r ,u (u i , v i ) = r u0 obtained using the expression for the Misner-Sharp mass in the main text. By solving this ODE along with all the others, we get the initial conditions at v = v i at each value of u along the integration.

Code tests
As a test of our numerical implementation we have reproduced the late-time decay of an asymptotically flat configuration with a massless scalar field. For this, it was crucial to employ the gauge conditions (A7), (A8). We also compute the "local power" of the scalar field decay, defined as −vΦ ,v /Φ. These are shown in Fig. 5 and are consistent with expected results.
To further test the code, we have analyzed its convergence properties. We thus evaluate the quantity for a given function F N obtained with resolution N at a fixed u coordinate, and where the maximum is evaluated for all values of v. Here, the index m refers to a reference solution obtained using a large number m of grid points while n denotes test solutions using a coarser resolution, n < m.
In Fig. 6 we show the convergence properties of the Kretschmann scalar for configurations B1. The plots show exponential convergence up to N ≈ 40.
Finally, since we use a free-evolution scheme, we have checked that the constraint equations (A4) and (A5) remain satisfied throughout our evolution. We show typical plots for the corresponding constraint violation in Fig. 7.
where M is the Misner-Sharp mass.
Appendix C: The scalar field along the EH The scalar field along the event horizon is shown in Fig. 8 (contrast with the asymptotically flat example of Fig. 5). The late-time behavior is described by ringing exponential falloff of the signal, well described by the lowest quasinormal modes of the spacetime [4].