{'fraction_non_alphanumeric': 0.04191832039114179, 'fraction_numerical': 0.008513085993672707, 'mean_word_length': 4.367888846005403, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Simplified Mathematical Model for the Formation of Null Singularities Inside Black Holes I - Basic Formulation and a Conjecture

Einstein's equations are known to lead to the formation of black holes and spacetime singularities. This appears to be a manifestation of the mathematical phenomenon of finite-time blowup: a formation of singularities from regular initial data. We present a simple hyperbolic system of two semi-linear equations inspired by the Einstein equations. We explore a class of solutions to this system which are analogous to static black-hole models. These solutions exhibit a black-hole structure with a finite-time blowup on a characteristic line mimicking the null inner horizon of spinning or charged black holes. We conjecture that this behavior - namely black-hole formation with blow-up on a characteristic line - is a generic feature of our semi-linear system. Our simple system may provide insight into the formation of null singularities inside spinning or charged black holes in the full system of Einstein equations.


Introduction
This paper examines a simple system of two equations inspired by the Einstein equations. The main purpose is to gain insight into the onset of null singularities inside spinning or charged black holes (BHs).
To understand the background and motivation for our toy model it will be worthwhile to review the development of our present conception of the null singularity inside BHS. The r = 0 curvature singularity of the Schwarzschild geometry has been regarded for many years as a prototype for the spacetime singularity expected to be present inside BHs. However, the Reissner-Nordstrom (RN) solution, describing a spherically symmetric charged BH, lacks a spacelike r = 0 singularity. Instead it admits an inner horizon (IH)-a perfectly smooth null hypersurface which constitutes a Cauchy horizon (CH) for partial Cauchy surfaces outside the BH. The (analytically extended) RN solution admits an r = 0 singularity too, but this singularity is timelike rather than spacelike, and it is located beyond the IH (hence outside the Cauchy development). A similar situation is found in the Kerr solution, describing a stationary spinning BH: A perfectly smooth IH, which again functions as a CH; and the spacetime singularity is timelike, located beyond this null hypersurface. In both the RN and Kerr solutions, the regular IH is known to be unstable to small perturbations, and this instability leads to the formation of a curvature singularity instead of a smooth IH. Thus, in order to explore the structure of the singularities inside realistic spinning BHs, one must understand the process of singularity formation due to the instability of the IH.
Of the three BH solutions mentioned above-Schwarzschild, RN, and Kerr-the one which is mostly relevant to realistic spinning BHs is obviously the Kerr solution. Nevertheless there is a remarkable similarity between the internal structures of spinning and charged BHs, which allows one to use spherical charged BHs as a useful toy model for the more realistic (but much more complicated) spinning BHs.
The IH of the RN solution is the locus of infinite blue-shift, as was already pointed out by Penrose [1]. Infalling perturbations of various kinds are infinitely blue-shifted there, which leads to instability of the IH [2]. As a consequence the latter becomes the locus of a curvature singularity, to which we shall often refer as the IH singularity. In order to explore this phenomenon, Hiscock [3] modeled the blue-shifted perturbations by a null fluid-a stream of massless particles. He analyzed the geometry inside a charged BH perturbed by a single such stream, an ingoing null fluid, using the charged Vaidya solution [4]. He found that the IH becomes a nonscalar null curvature singularity. Later Poisson and Israel [5] explored the system of a charged BH perturbed by two fluxes, namely both ingoing and outgoing null fluids. They concluded that in this case too the IH becomes a null curvature singularity. This time, however, the singularity is a scalarcurvature one because the mass-function -a scalar quadratic in derivatives of the area coordinate-diverges, a phenomenon known as mass-inflation. The detailed structure of this mass-inflation singularity was later analyzed [6] within a simplified model (in which the outgoing flux is replaced by a discrete null shell). This study showed that the metric tensor (when expressed in appropriate coordinates) has a continuous and non-singular limit at the singularity. Yet derivatives of the metric functions diverge at the IH, yielding a curvature singularity. The continuity of the metric has crucial physical consequences: It implies that the singularity is weak [7], namely an extended object will only experience a finite (and possibly very small) tidal deformation on approaching the IH singularity.
Subsequently more detailed numerical and analytical studies of the mass inflation phenomenon were performed, in which the perturbations were modeled by null fluids or by a self-gravitating scalar field [8,9,10,11,12]. These studies confirmed the conclusions of the earlier analyses ( [3], [6]). In addition, numerical analyses revealed that, at least in the case of scalar field perturbations, a spacelike singularity forms in the asymptotically-late advanced time. More recently Dafermos [13] proved for a characteristic initial value problem for the spherically symmetric Einstein-Maxwell-Scalar Field equations that for an open set of initial data on the event horizon (EH), the future boundary of the maximal domain of development becomes a null surface along which the curvature blows up. Dafermos proved that the metric can be continuously extended beyond the IH, namely, the singularity is weak.
The situation inside a spinning BH is similar in many aspects to that of a spherical charged BH. Here, again, the inner horizon is the locus of unbounded blue shift, suggesting that the regular IH of the Kerr geometry will become a curvature singularity when perturbed. A thorough perturbation analysis [14] [15] showed that indeed a scalar-curvature singularity forms at the early portion of the IH, which is again null and weak. This picture of the spinning IH singularity was later confirmed by an independent perturbative analysis by Brady et al. [16]. The existence of a class of solutions to the vacuum Einstein equations which admit null, weak, scalar-curvature singularities was also demonstrated in exact non-perturbative analyses [17,18] (though these exact analyses, unlike the earlier perturbative analyses, did not demonstrate the actual occurrence of a null weak singularity inside BHs).
If a cosmological constant Λ > 0 is present, the spacetime is no longer asymptotically flat and a cosmological horizon replaces the future null infinity. The spherical charged BH and the stationary spinning BH are then described by the Reissner-Nordström-de Sitter (RNDS) and Kerr-de Sitter solutions, respectively. In both cases there are three horizons, namely cosmological, event, and inner horizons. The surface gravity of these horizons depend on the parameters of the solutions, namely the cosmological con-stant, mass, and the charge or angular momentum. We denote these surface gravities by κ co (cosmological horizon), κ ev (EH) and κ in (IH). If κ in > κ co , there is an infinite blue-shift at the IH, suggesting an instability of the latter. This instability was first investigated by Mellor and Moss [19] in the case of spherical charged BHs, and by Chambers and moss [20] for spinning BHs, using linear perturbations in both cases. In the case of a spherical charged BH the non-linear instability with respect to ingoing null fluid was investigated by Brady and Poisson [21]. Brady, Núñez, and Sinha [22] investigated a model in which both ingoing and outgoing null fluids are present, and found that the mass function diverges provided that κ in > 2κ co . In the range 2κ co > κ in > κ co the mass function is finite, yet the Kretchman curvature scalar R αβγδ R αβγδ diverges at the inner horizon. Later, Brady, Moss, and Myers [23] considered also the contribution of the radiation that is scattered by the curvature in the vicinity of the EH. They found that when this scattering is taken into account, the necessary condition for stability of the IH (namely bounded curvature) is that both κ co and κ ev are greater than κ in . None of the stationary (asymptotically-de Sitter) electro-vacuum black holes satisfy this condition. Chambers [24] studied a simplified massinflation model with a continuous ingoing null fluid and a discrete outgoing shell and confirmed the earlier results [22]. He also found that for all values of κ in the metric functions are continuous and non-singular at the IH, even though the mass function diverges. Namely, the mass-inflation singularity is weak in the Λ > 0 case as well.
The combination of all the above-mentioned investigations strongly suggests (though a mathematical proof is still lacking) that the vacuum (or electro-vacuum) Einstein equations admit a generic class of solutions in which a null weak singularity forms inside a spinning (or charged) black hole. In what follows we shall assume that this is indeed the case. Now, the Einstein equation in 3+1 dimensions (and with the lack of any symmetry) is a rather complicated non-linear dynamical system. The following question therefore naturally arises: Is it possible to extract from the Einstein equations a smaller and simpler dynamical system, which is capable of producing black hole-like configurations with generic null weak singularities inside them? If such a simpler system is found, perhaps it could be viewed as the "active ingredient" of the Einstein equations (as far as the formation of black holes and null singularities is concerned). This may provide insight into the mathematical process of the formation of null singularities. The construction and exploration of such a simple system of equations is our main goal in this paper.
The system of Einstein equations combines both evolution and constraint equations. It appears likely, though, that the property of producing generic null singularities is admitted by the sub-system of evolution equations. We shall therefore extract our simplified toy-system from the evolution equations and simply ignore the constraint equations. Now, when the constraint equations are discarded, one obtains dynamical behavior even in sphericallysymmetric situations (it is the constraint section which "freezes" the dynamics in spherical symmetry). Consequently, we shall extract our toy-system from the evolution section of the electro-vacuum Einstein equations in spherical symmetry. We shall thus proceed in Sec. 2 as follows: We start from the electrovacuum Einstein equations in spherical symmetry. We also add a cosmological constant, for reasons explained below. Then we discard the constraint equations, and re-formulate the evolution equations in a simple form free of first-order derivatives. This yields a semi-linear hyperbolic system of two equations for the two unknowns which we denote R(u, v) and S(u, v), where u, v are two null coordinate and R, S are constructed from the metric functions (specifically R is the square of the area coordinate, and s will be specified below). The new system involves a "generating function" h(R) which in the above construction emerges in a very specific form [see Eq. (14)]. However, from the mathematical point of view it appears likely that the global properties of the solutions such as BH and singularity formation will not be sensitive to the detailed functional form of h(R), but only to certain global and/or asymptotic features of this function. For this reason we extend our view point and explore this semi-linear system with a rather general function h(R). This generalizes our investigation and simplifies it at the same time.
Our strategy of considering a general function h(R) also has a side benefit: As it turns out, certain two-dimensional general-relativistic dilatonic models can be re-formulated such that their evolution sector is described by our semi-linear system, with a certain function h(R). This includes the model by Callan et al. [25] and its charged generalization [26,27,28]. We describe this at the end of Sec. 2.
In Sec. 3 we describe some basic mathematical properties of our semilinear system, including conserved fluxes, a generalized mass function, and gauge freedom. The latter means that the semi-linear system is invariant under coordinate transformations of the form u → u ′ (u), v → v ′ (v). Then in Sec. 4 we construct, for any h(R), a class of exact solutions with vanishing fluxes, to which we shall refer as the flux-free solutions. This is a one-parameter family of solutions (for given h(R)), a generalization of the RNDS solution to arbitrary h(R). We then observe that for functions h(R) admitting three roots (or more), the corresponding flux-free solution describes a RNDS-like styatic BH, with three horizons, namely three null lines of constant R: an event horizon located at a line u = const, and cosmological and inner horizons, both located at v = const. (In the "Eddington-like" coordinates, in which the flux-free solution is first derived, the three horizons are located at infinite value of the relevant null coordinate, but this is later fixed by a coordinate transformation as described below.) The three horizons intersect at a single point P, representing the timelike infinity for the external region between the cosmological and event horizons. The function s diverges at P, but this divergence does not represent a spacetime singularity: Instead it reflects the fact that the proper-time interval between P and any point to its past is infinite. The function R is many-valued at P.
The singularity structure of the flux-free solution is studied in Sec. 5. The solution becomes singular at the horizons (s diverges), but this is merely a coordinate singularity. To regularize the solution we transform u, v to new, "Kruskal-like", coordinates. (Specifically we "Kruskalize" u with respect to the event horizon and v with respect to the cosmological horizon, such that the initial data for the BH formation are regular.) In these new coordinates the solution extends smoothly into the BH, and provides a description of the internal geometry up to the inner horizon.
The asymptotic form of the functions R, s near the IH is the primary objective of this paper. In the flux-free solution (expressed in Kruskal coordinates and extended into the BH as described above), R admits a constant finite value along the IH but s diverges there (for a generic h(R)). This divergence, too, does not indicate a true singularity, because it can be removed by "Kruskalizing" v with respect to the inner horizon. With such a coordinate transformation, the variables R, s become perfectly regular (in fact analytic) in the IH neighborhood. 1 In fact, this divergence of s reflects the infinite blue-shift (or red-shift in some cases) which takes place at the IH, just as in the standard RN and RNDS geometries. It should be noted that all invariant quantities involving the variables R, s and their derivatives are regular at the IH in the flux-free solution.
Consider now the initial-value problem for our semi-linear system. The initial hypersurface is taken to be a spacelike hypersurface which intersects both the event and the cosmological horizons. In the first stage we assume that the initial data agree with those of the flux-free solution. Then the evolving solution will be just the flux-free solution, with an IH of the form described above. We assume that the initial data are everywhere regular, which means that the flux-free solution is obtained not in the Eddingtonlike coordinates, but in other coordinates which are regular at the event and cosmological horizons-e.g. the above mentioned Kruskal-like coordinates.
The major challenge is now to understand how will the functions R, s be affected if the initial data are modified such that they no longer agree with those corresponding to the flux-free solution. What features of the BH and the IH will survive the perturbation, and which features will be modified? We do not have a full answer to this question, but we do have a conjecture that we present in Sec. 6, based on several compelling indications. These include the linear perturbation of the flux-free solution, the "generalized Vaidya solution" (valid for arbitrary h(R); see the appendix), and also some specific examples of h(R) for which the general solution may be constructed. Our conjecture may be stated very briefly as follows: First, the global blackhole structure is unchanged (this is manifested by the persistent divergence of s at a point P where the three horizons meet); Second, R remains finite (though no longer constant) along the IH; Third, the divergence of s on approaching the IH persists and preserves its leading asymptotic form; and after "re-Kruskalization" s becomes finite along the IH, just as in the unperturbed flux-free solution. However, one important difference occurs due to the deviation from the flux-free initial data: Although the variables R, s are continuous (after "re-Kruskalization") at the IH, they are no longer smooth. Certain invariant quantities involving the derivatives of R now diverge on the IH (this holds provided that the "surface gravity" of the IH is sufficienly large).
The h(R) function corresponding to the electro-vacuum solutions in four dimensions without a cosmological constant, namely Eq. (14) with Λ = 0, has only two roots. The corresponding BH solution has two horizons, the event and inner horizons. The cosmological horizon disappears when Λ vanishes, and instead there is a future null infinity. Our semi-linear system is useful in this case too, but the initial-value problem described by this system is conceptually more complicated in this case. To understand the reason, consider a black-hole solution with a cosmological horizon, and consider an initial spacelike hypersurface Σ which intersects both the event and cosmological horizons. We can pick a compact portion Σ 0 of Σ which still intersects the event and cosmological horizons. Then the early portion of the IH is included in the closure of D + (Σ 0 ), where D + denotes the future domain of dependence. On the other hand, in the analogous asymptotically-flat case the initial hypersurface Σ must extend to spacelike infinity (or alternatively to future null infinity) in order to have any portion of the IH being included in the closure of D + (Σ). This means that the behavior of R, s near the IH will depend on the asymptotic behavior of the initial data as the initial hypersurface approaches spacelike infinity. No such complication occurs in the case of a BH with a cosmological horizon: Here it is sufficient to require that the initial data are sufficiently regular on Σ 0 , and the issue of their large-R asymptotic behavior does not arise. For this reason, we shall restrict our attention in this paper to functions h(R) with at least three roots. Note that in this case the divergence of s at P and at the IH is, from the PDE point of view, a manifestation of the finite-time blow-up phenomenon, caused by the non-linearity of the hyperbolic system. (The standard General-Relativistic point of view is somewhat different, however, because the proper-time distance of P is infinite, due to the divergence of s, so this divergence is not a spacetime singularity.) As was mentioned above, we view our semi-linear system as a toy model for the much more complicated system of Einstein equations in four dimensions. Obviously not all properties of the Einstein equations are mimicked by our semi-linear system. The properties we expect our toy system to display are (i) the very formation of the BH (expressed in our system by the finite-time blow-up of s at the point P), (ii) the no-hair properties of the BH-namely the decay of external perturbations, and (iii) the generic formation of a null, scalar-curvature, weak singularity on the IH. We do not expect our toy system to properly address the spacelike singularity (which may intersect the IH at later times). Also this simple system is incapable of describing the oscillatory character of the null IH singularity inside a generically-perturbed spinning BH [29]. Nevertheless our toy system correctly demonstrates the basic properties of the IH singularity-in particular its weakness.
We conclude in Sec. 7 by outlining some directions for future research.

The Field equations 2.1 Maxwell-Einstein equations in spherical symmetry
We start by considering the Maxwell-Einstein equations in a spherically symmetric spacetime. We write the metric in double-null coordinates as where dΩ 2 = dθ 2 + sin 2 θdφ 2 . The Maxwell equations are easily solved, yielding with all other components vanishing. Here Q is a free parameter, to be interpreted as the charge. The electromagnetic energy-momentum tensor is then substituted in the Einstein equations with a cosmological constant Λ, This yields a system of two evolution equations, and two constraint equations, 2 The latter two equations (unlike the two evolution equations) are in fact ordinary equations along the lines v = const or u = const, respectively.

Constructing our semi-linear system
We first re-formulate the field equations such that no first-order derivatives appear in the evolution equations. To this end we introduce two new variables R and s instead of r and f : With the new variables, the evolution equations take the convenient form The constraint equations become In the next stage we simply omit the constraint equations (9) and keep the evolution equations (7,8) as our dynamical system. Recall that the evolution equations form a closed hyperbolic system, which uniquely determines the evolution of (properly-formulated) initial data. By this we achieve several goals: First, a non-constrained dynamical system is conceptually simpler to analyze than a constrained one; Second, this allows us to explore and test our hypothesis that the phenomenon of generic null-singularity formation inside four-dimensional spinning BHs is essentially a property of the evolution sector of the Einstein equations. 3 Finally, omitting the constraint equations retain the dynamics to the problem (it is the constraint sector which is responsible to properties like e.g. the Birkhoff theorem). This provides us with an effectively two-dimensional toy system aimed at mimicking dynamical properties of the Einstein equations in four dimensions.
Next we recognize that the expression in the parentheses in Eq. (8) is the derivative of the expression in the parentheses in Eq. (7) with respect to R. We can therefore write the two equations as where at this stage F (R) denotes the specific function and hereafter a prime denotes a derivative of a function of one variable with respect to this variable. These equations, which are semi-linear nonhomogeneous wave equations, constitute the core of our model. For later convenience we introduce another function, h(R), defined by its derivative: Note that h(R) is defined up to an integration constant. The semi-linear system now reads In the specific case (11) we have Note that this is 2r times the standard RNDS function where here the integration constant was expressed as −4m, m being the ADM mass of the corresponding static RNDS solution.
The final stage in constructing our toy system is to abandon the specific function h(R) of Eq. (14) and instead to explore the semi-linear system (13) for a general function h(R). Again we define F ≡ −dh/dR, or This generalization is advantageous for several reasons. First, if indeed the system (13) leads to generic null singularities, it is plausible that this property will not be sensitive to the specific functional form (14). Rather, we expect the qualitative properties of our dynamical system to depend only on certain qualitative features of h(R). Note also that our primary goal is to provide a simple toy model aimed at mimicking certain dynamical features of e.g. the vacuum Einstein equations in four dimensions, and from this perspective the spherically-symmetric electro-vacuum system of the previous subsection should itself be regarded as a toy model; hence there is no reason to firmly stick to the specific function (14). Second, this extension of our view-point will allow us to seek simple examples of functions h(R) for which the general solution of the system (13) may be constructed. Such solvable examples would provide valuable insight into the dynamical properties of this system.

Application to two-dimensional black holes
In addition to its role as a toy model for singularity formation, our semilinear system is also directly applicable to certain dilatonic models of twodimensional BHs. In the model developed by Callan et al [25] there is a dilaton φ(u, v), a cosmological constant of the two-dimensional model λ, and the metric is The classical, matter-free, Einstein equations then yield two evolution equations and two constraint equations. Transforming to the new variables R = e −2φ and S = 2(ρ − φ), the constraint equations reduce to Eq. (9), and the evolution equations take the form (13), this time with the generating function This dilatonic two-dimensional model was later generalized to include a Maxwell field as well as charged matter fields [26,27,28]. Here, again, with the same substitution R = e −2φ , S = 2(ρ − φ) the classical matter-free Einstein equations are reduced to Eqs. (9,13) with the generating function where Q is a parameter proportional to the Maxwell field's charge.

Basic mathematical properties
In this section we introduce some basic features of our semi-linear system (10).

The gauge freedom
Our semi-linear system (10) is invariant under a family of gauge transformations. These are coordinate transformations which preserve the double-null form of the metric: u →ũ(u), v →ṽ(v). The variable R is invariant under this coordinate transformation, but s changes. Since e s ∝ g uv , it transforms like a covariant tensor of rank two, and one finds: The various quantities made of R and s may be classified according to the way they transform under a gauge transformation. The simplest are the scalars, namely quantities which are unchanged. Obviously R is a scalar. Apart from R itself, there is only one scalar made of R and s and their firstorder derivatives: e −s R ,u R ,v . Another useful, non-scalar, quantity is e −s R ,w , where hereafter w stands for either u or v. This quantity is invariant to a transformation of w, but not to transformation of the other null coordinate.

The conserved fluxes
Consider the quantities Differentiation Φ with respect to u, one observes that identically vanishes by virtue of the field equations (10). In a similar manner one finds that the derivative of Ψ(u) with respect to v vanishes. Namely, We shall refer to Φ(v) and Ψ(u) as the two conserved fluxes (or simply fluxes). It is sometimes useful to express these fluxes as One can easily verify that in a gauge transformation u →ũ(u), v →ṽ(v) the two fluxes transform asΦ (namely like components of a covariant second-rank tensor). Note that Φ(v) and Ψ(u) are uniquely determined by the initial data for R and s (this is most easily seen when the characteristic initial-value formulation is used [30]).
Application to the spherically-symmetric charged case: In the four-dimensional spherically-symmetric case, an important problem is that of the RN solution perturbed by two fluxes of null fluids, namely ingoing and outgoing fluxes. In this case, the dust contribution to the energy-momentum tensor is [5] where r is the area coordinate, u and v are two null coordinates, and L in , L out denote the ingoing and outgoing dust fluxes, respectively. From the Einstein equations for T vv and T uu with the line element (1) one finds that These quantities are directly related to the conserved fluxes Φ(v), Ψ(u) discussed above. In fact one can easily show, using Eqs. (6,21), that We may therefore regard the quantities Φ, Ψ as the generalization of the spherically-symmetric null-fluid fluxes to arbitrary h(R).
It is important to recall that in the spherically symmetric case (11) the semi-linear hyperbolic system (10) is mathematically equivalent to the massinflation model [5] with two arbitrary fluxes L in and L out .

The generalized mass function
Consider the scalar quantity where h 0 (R) is a certain member of the one-parameter family (16) (namely one associated with a certain choice of the integration constant). This is a generalization of the mass parameter to dynamical cases (see below). Note that for a given F (R) the mass function is defined up to an additive constant (associated with different choices of h 0 (R)). One can easily show, using the field equation (13) for R, that the derivatives of M satisfy Also, differentiating the last equation with respect to v and recalling Eq. (24), one observes that M satisfies the simple field equation From Eq. (32) we see that when the fluxes Φ(v) and Ψ(u) vanish, the mass function becomes a fixed parameter. This is the situation in the "fluxfree solution" described below (section 4). Also when one flux vanishes (e.g. Ψ), the mass function only depends on one null coordinate (v in this case), which is the situation in the generalized Vaidya solution (discussed in the Appendix).
Application to the spherically-symmetric charged case: As an illustration we consider here the mass function for spherically symmetric charged black holes. In this case one naturally defines h 0 by omitting the constant in Eq. (14), namely In terms of the original variables r, f the mass function then reads In the case of Λ = 0, the last expression is just four times the function m(u, v) defined by Poisson and Israel [5]. 4

The flux-free solution
In this section we investigate a class of solutions which is the generalization of the static RNDS family to general h(R). These are the solutions in which both Ψ and Φ vanish. We first construct these solutions in Eddington-like coordinates and then transform to Kruskal-like coordinates. We then explore the singularities of these solutions.

Construction in Eddington-like coordinates
In the case considered here, Eqs. (24,25) read The first integral of these two equations is The equation for R ,v is invariant to a transformation of v, but a transformation of u affects the function c u (u). However, the signs of c u is preserved in such a gauge transformation, because we require the new null coordinateũ to be future-directed, just like the original u. The situation with the equation for R ,u is exactly the same (with the obvious interchange of u and v). Thus, with the aid of a gauge transformation we can bring both functions c u (u) and c v (v) to ±1, with the signs corresponding to those of the original functions.
Consider first the case where R is increasing with v and decreasing with u, namely c v > 0 and c u < 0 (this is typically the situation outside a BH, though no further than the cosmological horizon -namely region I in Fig.  1). Then the gauge transformation described above leads to This implies that both R and s are functions of a single variable, and these two functions satisfy Then the field equation (13) for R ,uv reads yielding the first order ODE Thus, R(x) is given by its inverse function 5 Then s(x) is given by s = ln(h(R)) .
The second field equation, namely Eq. (13) for s ,uv , is automatically satisfied, as one can easily verify. Note that this solution is static, in the sense that it only depends on the spatial variable x = v − u. Note also from Eqs.
(43) or (45) that h must be positive in this case. Next let us consider the case where R is decreasing with both v and u, which is typically the situation inside a BH (region II in Fig. 1). Then instead of Eq. (39) we now get Correspondingly we now define and both R and s are functions of x only, satisfying this time. Substituting again in the field equations (13), one finds that Eqs.
(43) and Eq. (44) for x(R) still hold, but there is a sign change in the expression for s, namely s = ln(−h(R)). Note that h is negative in this case. The general expression for s which holds in both cases is obviously Thus, the flux-free solution comes in two versions: The "external-type" version, which depends on the spatial variable v − u, and the "internaltype" version, which depends on the temporal variable v + u. The first version occurs in regions where h > 0, and the second occurs when h < 0 (typically inside a BH). Equations (44) and (49) hold in both cases. 5 The full integral of Eq. (42) obviously involves two integration constants. The first one is already embodied in the definition of h in Eq. (43). The second one is an arbitrary constant to be added to the right-hand side of Eq. (44). This constant has no physical meaning, however, as it may be absorbed by a gauge transformation e.g. v → v + const, which merely shifts x by a constant.
In its both versions, the flux-free solution is in fact a one-parameter family of solutions (for given F (R)), due to the arbitrary integration constant in the definition of h in Eq. (16). The situation is somewhat confusing because to a specific flux-free solution a function h(R) is associated in two different ways: (i) through the field equation (13) (to be satisfied by the solution), and (ii) through its explicit construction via Eq. (44). It must be noticed that h(R) appears in two conceptually different ways in these two occasions: In the field equation it appears as an equivalence class (because the field equations only depend on dh/dR), and in the construction procedure it appears as a one-parameter family of distinct functions. To clarify this notational confusion we reserve h(R) of Eq. (16) to denote the equivalence class, and introduce the notation where h 0 (R) is a representative of the equivalence class h(R) (already mentioned above), and c is an arbitrary constant which distinguishes the various members of this equivalence class. Thus, we rewrite the flux-free solution as (with x = v ± u as above), namely The parameter c may be interpreted as one associated with the system's mass. Indeed, the mass function (31) reads for flux-free solutions In the spherically-symmetric four-dimensional case, the flux-free solution is just the RNDS family of solutions, whose ADM mass is related c through m = −c/4. We therefore write H in this case as

Horizons
The solutions constructed above become pathological at any value R = R 0 for which H vanishes. From Eq. (51) s diverges there to −∞. Also Eq.
(52) implies (assuming finite F (R 0 )) that x diverges at R = R 0 , meaning that either u or v is unbounded there. This phenomenon is analogous to the coordinate singularity at the horizon of the Schwarzschild solution, when the metric is expressed in double-null Eddington coordinates. In our case, too, this coordinate singularity may be overcome by transforming to new, Kruskal-like, null coordinates, as shown in the next subsection. To this end, however, we must first analyze the asymptotic behavior of s and x on approaching the horizon in the original Eddington-like gauge. We define (for each horizon): and assume K = 0, therefore in the neighborhood of the horizon. Note that for K > 0, dx/dR = 1/H is negative for R < R 0 and positive for R > R 0 , hence x → −∞ at both sides of the horizon. Similarly, for K < 0 at both sides x → +∞. Therefore Kx → −∞ for both K > 0 and K < 0, and at both sides of the horizon.
Thus, both s and x diverge logarithmically in R − R 0 . The occurrence of sign flips at the horizon in some of the above expressions complicates the analysis. To help clarifying this confusion we define, for each horizon, the quantity X (unlike x) is continuous and monotonous across the horizon, and it vanishes at the horizon itself. The comparison of Eqs. (58) and the definition of X immediately yields dR/dX = ±1 at the horizon, and a closer look at the signs reveals that at both sides Since dX/dx = KX, R(X) satisfies the same differential equation at both sides: This, combined with Eq. (60), implies that X(R) is analytic across the horizon (provided that H(R) itself is analytic in a neighborhood of R = R 0 , which we assume). Note that H/X is analytic too, and it gets the nonvanishing value at the horizon. We shall primarily be interested in functions H(R) admitting (at least) three simple roots R i (i = 1, 2, 3), ordered R 3 < R 2 < R 1 , such that H is positive at R 2 < R < R 1 and negative at R 3 < R < R 2 , as shown in Fig. 2. An archetype is the function h(R) corresponding to the sphericallysymmetric electro-vacuum solutions, Eq. (54), which (for sufficiently small Q and Λ, and sufficiently large m) admits three roots. These roots correspond to the cosmological horizon (R 1 ), the event horizon (R 2 ), and the inner horizon (R 3 ), as shown in Fig. 1. The horizons divide the spacetime into three regions which we denote I,II,III, as shown in Fig. 1. We shall primarily be concerned here with the regions I and II (region III will not concern us here, except at the very neighborhood of the cosmological horizon). Note that H is positive in region I and negative in regions II and III.
As was shown above, in the Eddington gauge s diverges on the three horizons. The divergence at the IH does not pose any difficulty-in fact investigating this divergence and its physical implications is one of our primary goals. However, the divergence of s on the event and cosmological horizons does pose an undesired feature: We would like to explore a situation in which a (flux-free) solution of the type described above, which includes the three horizons, emerges from regular initial data prescribed on some compact spacelike initial hypersurface Σ 0 . Furthermore we want the horizons' intersection point to be included in D + (Σ 0 ). To this end, Σ 0 must intersect both the event and cosmological horizons (but not the inner horizon). The divergence of s (and also u or v) on these two horizons renders the Eddington gauge inappropriate for such a regular initial-value set-up. We shall therefore proceed now to transform the Eddington coordinates into Kruskal-like coordinates with respect to the event and cosmological horizons.

Transforming to Kruskal-like coordinates
The construction of the Kruskal-like coordinates in our case is similar to the standard procedure in e.g. the Schwarzschild spacetime-except that here we "Kruskalize" u with respect to the EH and v with respect to the cosmological horizon.
Let us define on each horizon R = R i : (namely, it is the |K| value associated with the i'th horizon.) 7 Consider first the EH, R = R 2 . Here K > 0, hence x diverges to −∞. In both regions I and II the Eddington coordinate v is regular along the EH but u diverges (see Fig. 1). In region I H > 0, hence x = v − u, and the divergence of x means that u → +∞. On the other hand in region II H < 0, hence x = v + u, and the divergence of x now implies that u → −∞. Correspondingly we define in region I, and in region II. Then U is continuous across the EH, and is monotonouslyincreasing (namely future-directed) everywhere; it is negative at region I, positive at region II, and vanishes at the EH. This transformation cures the divergence of s, as we show below.
Next we consider the cosmological horizon, R = R 1 . Here K < 0, hence x diverges to +∞, meaning that v → +∞ on approaching the horizon from region I, whereas u is regular. 8 We thus define (in region I) Again, V is a future-directed null coordinate which takes negative values in region I and vanishes at the cosmological horizon. Since V is continuous across the EH, it takes negative values in region II as well.
The variable R is invariant under the coordinate transformation (u → U, v → V ). Therefore, R is formally given as a function of U and V through We now proceed to show the regularity of R and S, and their smoothness as functions of U and V , at both the event and cosmological horizons. Considering the EH first, we denote by X 2 the function X(R) (defined above) associated with the EH. Noting that at the EH K > 0 and hence k 2 = K, we write One finds (treating carefully the flipping signs) that where hereafter k ij ≡ k i /k j > 0 for any i, j = 1...3. Since V is strictly negative along the EH, and R is an analytic function of X 2 (as establish above for a general horizon), we conclude that R(U, V ) is analytic in the neighborhood of the EH. To analyze the variable S, we note that As was established in the previous subsection, H/X 2 is a regular function of R (or X 2 ) which takes the value K = k 2 > 0 at the EH. Therefore S(U, V ) too is a regular (in fact analytic) function of U and V . The cosmological horizon is treated in an analogous manner, except that here we do not need to explicitly analyze the various functions in region III: Instead we simply extend the relevant functions analytically from region I into region III. Note that here K < 0 and therefor k 1 = −K. Correspondingly we find (for region I) This yields Again we see that R(U, V ) is analytic in the neighborhood of the cosmological horizon, because R(X 1 ) is analytic (and U is strictly negative). Recalling that we obtain Again, −h/X 1 is a regular function of R (or X 1 ) which takes the value −K = k 1 > 0 at the cosmological horizon. Therefore we conclude again that both R(U, V ) and S(U, V ) are analytic functions in the neighborhood of the cosmological horizon. Then in region III R(U, V ) and S(U, V ) are defined to be the analytic extension of the corresponding functions from region I across the cosmological horizon.
In fact it is straightforward to show that R(U, V ) and S(U, V ) are regular at U, V < 0 not only in the neighborhood of the event and cosmological horizons, but also in the entire range R 2 ≥ R ≥ R 1 and its neighborhood (provided that H(R) itself is regular throughout). Furthermore, for any smooth spacelike initial hypersurface Σ 0 which intersects the event and cosmological horizons, the initial data for R and S (corresponding to the flux-free solution in the U, V coordinates) are regular throughout. 9

Singularities in the flux-free solution
In this section we analyze the singularities that appear in the flux-free solution using the expressions that were derived in the previous section. As before, we shall assume that H(R) has three simple roots at R 1 , R 2 and R 3 , with H positive at R 2 < R < R 1 and negative at R 3 < R < R 2 (additional roots at R > R 1 and/or R < R 3 are allowed; see last subsection). We divide the discussion into two types of singularities -the vertex singularity and the IH singularity.

The vertex singularity
In the Kruskal-type coordinates defined in the previous section (U, V ) = (0, 0) is the intersection point of the three horizons (see Fig. 1). At this point R is many-valued (for example it admits a fixed value R = R i along the i'th horizon). Furthermore, S → +∞ as U, V → 0. This divergence proceeds in a slightly different manner along various paths towards the vertex (0, 0). For example, along curves of constant R = R i , S ∝ log |U | + log |V |, as can be seen from Eq. (68). Also, Eqs. (71) and (74) imply that along the event horizon, and along the cosmological horizon. Since k 12 > 0, S → +∞ along both horizons on approaching U = V = 0.

The singularity at the inner horizon
Next we analyze the asymptotic behavior of the flux-free solution on approaching the IH from region II, namely the limit V → 0 − in the range U > 0. To this end we use the same Kruskal-like coordinates (U, V ) introduced above, which allow for regular flux-free initial data on Σ 0 . From Eq. (61) applied to the event horizon we deduce that R decreases throughout region II, hence R ≥ R 3 . We denote by X 3 the function X(R) (defined in the previous section) associated with the inner horizon. Noting that at the IH K < 0 and hence k 3 = −K, we write for region II (recalling Substitution in Eq. (68) (which is valid in region II as well) yields From Eq. (62) we obtain along the IH H(R)/X 3 = −k 3 . Therefore, for any curve in region II which approaches the IH (at some U > 0), the term in squared brackets approaches a finite value as V → 0. Along such a curve e S behaves as We find that S → ∞ if k 1 > k 3 , S → −∞ if k 1 < k 3 and S is regular in the case of k 1 = k 3 . The singularity at the IH is locally gauge-removable; Namely, it may be removed by transforming from V to a new Kruskal-type coordinate defined with respect to the IH. One then obtains in the (U, V 3 ) coordinates, in full analogy with Eqs. (74) and (76), which admits the regular limit at the IH. However, this transformation will destroy the regularity of S on the cosmological horizon. In particular it will spoil the regularity of the initial data on Σ 0 as the latter intersects the cosmological horizon. Hence the divergence of S at the IH (for k 1 = k 3 ) is globally non-removable. From the perspective of the finite-time blow-up phenomenon which concerns us here, it is mandatory to introduce regular initial data, and the divergence of S on the IH is inevitable. We now introduce three additional quantities and explore their asymptotic behavior at the IH. These quantities will serve as gauge-invariant indicators for certain aspects of regularity or irregularity at the IH.
i) e −S R ,V R ,U : This quantity is interesting because (unlike S) it is a scalar, namely invariant under a gauge transformation. We can therefore calculate it directly from the original Eddington-like gauge through Eq. (51): which actually vanishes on the IH.
(ii) e −S R ,V : This quantity is not a scalar, yet it is invariant under a transformation of V . In a gauge transformation u →ũ, v →ṽ it is simply multiplied by dũ/du. We restrict attention here to transformations u →ũ which are regular in the interior of region II (i.e. at U > 0). We then conclude that although the actual value of e −s R ,v is modified, the divergence of this quantity at the IH is invariant to a gauge transformation. In the Eddington gauge Eq. (51) implies e −s R ,v = −1. In the Kruskal-like gauge we obtain e −S R ,V = −dU/du = −k 2 U .
(iii) S ,U : This quantity too is invariant under a transformation v →ṽ, because the term −ln(dṽ/dv) in Eq. (20) is independent of U . Therefore its regularity at the IH is gauge invariant. We calculate this quantity in the Eddington gauge: which is presumably regular. The regularity of S ,U (or, more generally, of its gauge-transformed counterparts ,ũ ) is a useful indicator for whether the divergence of S at V = 0 is locally gauge-removable or not. To show this, we choose in region II a line U = U 0 > 0 which intersects the IH. We then define the new coordinateṽ to be the affine parameter along this line, setting e.g.ṽ = 0 at the IH. Theñ s = const ≡s 0 everywhere along this line, and in particular atṽ → 0. Now, at any U > 0 we haves Therefore regularity ofs ,ũ implies regularity ofs at the IH, and vice versa. 10 On the other hand, the divergence of e −S R ,V R ,U and e −S R ,V indicates (qualitatively speaking) to what extent the gauge-invariant variable R is smooth at the IH. Recall that after the divergence of S has been removed (locally) by a gauge transformation V → V 3 , the solution (R, S 3 ) is continuous at the IH, so we are left with the issue of the next-level regularity, namely C 1 smoothness at the IH. The quantities e −S R ,V R ,U and e −S R ,V both provide information on this issue, in a gauge-invariant manner. 11 In the flux-free solution all the above indicators (i-iii) are regular (even though S itself diverges). As we shall discuss in the next section, this situation presumably changes for indicators (i) and (ii) in the more general flux-carrying solutions.
The discussion above was based on the Kruskal-like coordinates U, V . It should be emphasized, though, that the singularity structure is in fact gauge-invariant, in the following sense: For any new coordinatesŨ (U ),Ṽ (V ) related to U, V in a regular manner, 12 the structure of both the vertex and the IH singularities will be the same as in the original coordinates. In other words, any choice of null coordinatesŨ ,Ṽ for which the flux-free initial data are regular on Σ 0 and its immediate neighborhood, will lead to the same structure of singularity. (Though, obviously the vertex singularity will shift from (0,0) to some other pointŨ 0 ,Ṽ 0 .)

The relation between the PDE and general-relativistic notions of singularity
The issue which primarily concerns us in this paper is that of finite-time blowup in a non-linear hyperbolic system of PDEs mimicking the Einstein equations. To this end we throughout adopt the standard terminology of 10 Obviously in the flux-free solutions the local gauge removal of the divergence of S can be demonstrated explicitly, see Eq. (81). Yet the indicator S,U is useful for the wider class of flux-carrying solutions. 11 We consider here both quantities (i,ii), even though both essentially probe the smoothness of R at the IH, because each of them has its own advantages and disadvantages. The quantity (i) is advantageous because it is a scalar. On the other hand in certain cases e −S R,V diverges, indicating the non-smoothness of R, and yet e −S R,V R,U vanishes. This happens because (qualitatively speaking) although the V -derivative of R diverges, its Uderivative vanishes. This situation typically happens in Vaidya-like solutions where only ingoing flux is present. In this sense the indicator e −S R,V is more robust.
12 by this we also mean thatŨ (U ) is monotonic with non-vanishing derivative, and the same forṼ (V ).
PDEs in considering the onset of irregularities ("singularities") in the evolving solutions. This terminology is quite different from that used in General relativity (GR) for discussing spacetime singularities. In this subsection we discuss how the various irregularity phenomena admitted by our PDE toy system are interpreted from the GR point of view.
We have seen before that the three-roots flux-free solutions typically develop two types of singular phenomena: the vertex singularity and the IH singularity. It should thus be emphasized that both phenomena do not imply a General-relativistic spacetime singularity, as we now discuss.
The divergence of S at the vertex U = V = 0 fails to mark a General-Relativistic singularity because (in a typical application of our system to a general-relativistic spacetime) no divergence of curvature is involved. (Note that S is a gauge-dependent quantity, and gauge-invariant quantities like e −S R ,V R ,U do not diverge at the vertex.) The divergence of the metric function g uv ∝ e S simply implies that the proper-time distances between the vertex and points in region I are infinite-namely, the vertex becomes future timelike infinity.
Similarly, the divergence of S at the IH fails to mark a spacetime singularity because it may be locally removed by a gauge transformation V → V 3 .
Smoothness indicators like e −S R ,V R ,U or e −S R ,V (and others not mentioned here) all indicate that subsequent to such a coordinate transformation the manifold is perfectly regular (in fact analytic) at the IH.
Despite the failure of the vertex and IH singularities to mark genuine GR singularities, both phenomena actually bear crucial implications to the structure and features of the GR black-hole spacetime, as we now discuss.
The vertex singularity functions as timelike infinity, as we already mentioned. In fact it defines the event and cosmological horizons-both are null lines which intersect the vertex singularity at their future. 13 Furthermore, from the standard hyperbolic-PDE point of view the occurrence of a singularity at the vertex implies the presence of a Cauchy horizon in the (U, V ) manifold at V = 0. This property carries over to the GR terminology-the inner horizon at V = 0 is indeed a Cauchy horizon of the spacetime. Thus, the very presence of a Cauchy horizon in e.g. the RNDS family of black holes may be attributed to the divergence of S at the vertex.
Next let us consider the General-relativistic implications of the IH singularity. As was mentioned above the divergence of S at V = 0 may be gauge-removed in a local manner only: Curing the divergent S at the IH will inevitably lead to blow-up of S at the cosmological horizon. The mismatch of S between the inner and cosmological horizons is an inherent, gauge-invariant, global feature of the typical three-root flux-free solutions. In the general-relativistic language this is translated into infinite blue-shift which takes place at the IH of the RNDS black holes (as well as their rotating counterparts). 14 This divergent blue-shift, in turn, leads to a null curvature singularity at the IH when the BH is generically perturbed. But this relativistic inner-horizon singularity is represented in our PDE terminology by the divergence of quantities like e −S R ,V R ,U and e −S R ,V (which, as was mentioned above, only diverge when fluxes are added).

Three roots versus four roots
The analysis and discussion so far was independent of the behavior of H(R) at R > R 1 . The extensions of the spacetime diagram in Figs 1 or 3 to the right of the cosmological horizon will depend the properties of H(R > R 1 ). We may conceive of several options: (i) H(R) admits an additional root at R 0 > R 1 ; (ii) No additional root exists at R > R 1 ; either (iia) H(R) diverges (or admits an irregularity) at a certain finite R s > R 1 ; or (iib) H(R) is regular-and negative-for any finite R > R 1 (the asymptotic behavior of this function as R → ∞ will not concern us here).
Consider first the case (i). In this case the dynamics in region III is essentially the same as in region II. In particular, the line H 3 in Fig. 3 (the IH) which borders region II has a counterpart at the future boundary of region III-namely the null line located at U = 0, V > 0, which we denote H 0 . The behavior of R and S on approaching H 0 just reflects that at H 3 , with the obvious interchange of the horizons' indices 3 → 0 and 1 → 2, 2 → 1. In particular R gets the fixed value R 0 , whereas S diverges. The sign of divergence is −∞ (implying infinite blue-shift) for k 0 > k 2 and +∞ (infinite red-shift) for k 0 < k 2 . Essentially the spacetime diagram of Fig. 3 (with the H 0 singularity added) would exhibit a reflection symmetry with respect to a vertical line which passes through the vertex at U = V = 0. The singularity will have a symmetric "V-shape" in this case, namely two null singularities which emerge from a common vertex singularity.
Note that additional roots at R > R 0 will not affect the spacetime diagram because we terminate the solution at the H 0 singularity, namely at R = R 0 . For the same reason, additional roots at R < R 3 will not affect the spacetime diagram of Fig. 3 (regardless of whether a root R 0 > R 1 exists or not).
Turn now to discuss the case (ii), namely H(R) has no roots at R > R 1 . In this case region III will terminate at a spacelike boundary line, corresponding to R = R s in case (iia) or R → ∞ in case (iib). This spacelike line, which intersects the vertex U = V = 0, will constitute a singularity of our PDE system. 15 Recall, however, that in the GR terminology the spacelike boundary will not necessarily be a spacetime singularity; Instead it may well mark a regular boundary at timelike infinity.
The function H(R) of Eq. (54) corresponds to the RNDS black-hole solutions. In this case there are only three roots, corresponding to the cosmological, event, and inner horizons. This function belongs to the class (iib), and the future spacelike boundary of region III marks the future timelike infinity of the external de Sitter universe.

The conjecture about the singularity formation
The semi-linear system (13) is the main objective of this paper and the subsequent paper [30]. In the previous sections a class of special solutions was constructed which demonstrate a finite time blow-up. We saw various features of the singularity which evolved from regular initial data. The solutions that we constructed were flux-free, namely "static". They all correspond to flux-free (Ψ = Φ = 0) initial data. This leaves us with the following open question: How will these solutions change if we perturb them by modifying the initial data, such that (Ψ, Φ = 0)? We are particularly interested in the finite-time blowup phenomenon and the structure of the singularity exhibited by the flux-free solutions: To what extent it is stable to the perturbations?
To be more specific, and to set the notation for the discussion below, let us consider a specific generating function h(R) and a specific three-root (or more) flux-free solution H(R) = h 0 (R) + c. The flux-free initial data are prescribed on an initial hypersurface Σ 0 which intersects the cosmological and event horizons (where R = R 1 and R = R 2 , respectively). The initial data need not correspond to the Kruskal-like coordinates U, V , but could be in any other coordinatesŨ ,Ṽ which are regular functions of respectively U and V -namely, we only require the flux-free initial data to be regular through Σ 0 . The evolving solution then develops a vertex singularity at P ≡ (Ũ 0 ,Ṽ 0 ), and a null IH singularity atṼ =Ṽ 0 . Now we add small (but finite), regular, perturbations to the initial functions on Σ 0 such that (Ψ, Φ = 0). Which properties of the flux-free singularity will survive the perturbations and which will be modified?
We conjecture that the basic structure of a vertex singularity and a null IH singularity is rather robust and will survive the perturbation (see Fig.  3). However, some of the features of the IH singularity will be modified. In particular, the character of the flux-free IH singularity of being fully locally-removable (by means of a local gauge transformation) will no longer hold when fluxes are added. A more detailed description of the robust and non-robust properties of the flux-free singularity is given below.
Our conjecture is based on several sources of evidence: (i) A linear perturbation analysis [31]. Note that at the linear level the two fluxes Ψ and Φ do not interact. The effect of each flux may be separately explored by means of the exact Vaidya-like solutions (see Appendix A), and then these two perturbations are simply superposed. This drastically simplifies the analysis, as no PDE needs be solved. In addition, the perturbation scheme suggests that the effects of non-linear coupling between the two fluxes are negligible compared to the basic effect of the blue-shifted incoming flux. (This insight-namely the dominance of the linear perturbation at the IH over the non-linear effects-has emerged in previous studies of the IH singularity inside perturbed spherical charged BHs [6,11], and, more generally, of four-dimensional general-relativistic null spacetime singularities [14,17,18]; and it applies to our toy semi-linear system as well.) (ii) The general exact solution of our semi-linear system for the case h = cos(R) [31].
This case is exactly solvable, and it demonstrates all the elements of the conjecture-except that it satisfies k 3 = k 1 hence it does not test issues related to the blue-shift phenomenon. Nevertheless it provides support to our conjecture concerning all properties of the vertex singularity, as well as some of the properties of the IH singularity-in particular its null character. (iii) Mathematical study of the general, flux-carrying solutions in the case of a saw-tooth function h(R) [30]. In this case h(R) is made of three linear sections, patched together to form a continuous but non-smooth function. This type of h(R) confirms all elements of the conjecture below-for both cases k 3 = k 1 and k 3 = k 1 .
In what follows we present in detail our conjecture about the singularity structure when the flux-free solution undergoes a generic small perturbation. We divide the discussion into the vertex singularity and the IH singularity.

The perturbed vertex singularity
The perturbed solution exhibits a finite-time blow-up and develops a vertex singularity at a slightly-modified point 16 . where • s diverges logarithmically to +∞ on approaching P ′ (e.g. from the past); • The function R is many-valued.
We stress again that this divergence of s does not represent a spacetime singularity -it merely reflects the fact that the proper time interval between P ′ and any point to its past is infinite. Therefore P ′ functions as timelike infinity in the BH spacetime. The presence of the singular point P ′ defines the horizons, and in this sense the very formation of the black hole. More specifically, the cosmological and event horizons are now defied to be the two null lines which intersect P ′ at their future [namely, (Ũ <Ũ ′ 0 ,Ṽ =Ṽ ′ 0 ) and (Ũ =Ũ ′ 0 ,Ṽ <Ṽ ′ 0 ), respectively]. These two null lines are no longer lines of constant R. Nevertheless we conjecture that • along each of these lines R admits a well-defined value R ′ i (i = 1...2) on approaching P ′ ; • Furthermore these two values are both roots of the same function H ′ (R) = h 0 (R) + c ′ , with c ′ close to the original c (hence R ′ 1,2 are close to the two original roots R 1,2 ).
From the GR point of view, this conjectured structure of the perturbed vertex singularity is a manifestation of the no-hair principle to our toy system; Namely all initial perturbations (the non-vanishing fluxes) should disperse away and leave a black-hole spacetime which at late time is well described by the corresponding static solution.

The perturbed IH singularity
First we conjecture that, when a flux-free solution with an inner horizon is weakly-perturbaed, no spacelike singularity will form in the neighborhood of P ′ . Namely, in the neighborhood of P ′ the perturbed solution will be regular throughoutṼ <Ṽ ′ 0 ,Ũ >Ũ ′ 0 . However the solution will generically develop a null singularity at the IH, namely atṼ =Ṽ ′ 0 (forŨ >Ũ ′ 0 ). This singularity will be similar in many respects to its flux-free counterpart: • R has a well-defined finite limit R IH (Ũ ) atṼ →Ṽ ′ 0 ; • Furthermore asŨ →Ũ ′ 0 (the P ′ limit) R IH (Ũ ) approaches a limiting value R ′ 3 which too is a root of H ′ (R), close to the original root R 3 .
• Just as in the unperturbed case, the function s generically diverges on the IH. The sign of divergence is determined by the relation between k ′ 1 and k ′ 3 , where k ′ i is |dH ′ /dR| at the root R = R ′ i : Namely s → −∞ for k ′ 3 > k ′ 1 and s → +∞ for k ′ 3 < k ′ 1 . 18 These two cases correspond to infinite blue shift or infinite red shift, respectively.
• In the special case k ′ 1 = k ′ 3 we expect s to remain regular atṼ →Ṽ ′ 0 . In this case the IH becomes a perfectly regular Cauchy horizon of our nonlinear hyperbolic system.
• Returning to the generic case k ′ 1 = k ′ 3 : Although s diverges, this divergence is locally removable by a gauge transformatioñ 17 In fact, the same property is shared by R ′ 3 , namely the limiting value of R on approaching P ′ along the perturbed IH, as discussed in the next subsection. 18 If the original flux-free solution satisfies k3 < k1 or k3 > k1, then we may assume that the perturbation is sufficiently small-and correspondingly c ′ is sufficiently close to c-such that the same inequality is satisfied by the perturbaed quantities k ′ 1,3 . Obviously in such a case, in the above condition for the sign of the divergent s we may replace k ′ 1,3 by k1,3.
But this transformation leads to a divergence of s along the cosmological horizon. Namely the divergence of s at the IH is globally nonremovable.
• The indicator ∂s/∂Ũ is regular at the IH limitṼ →Ṽ ′ 0 . As discussed in the previous section, this indicates that indeed the divergence of s is locally removable.
All the above features merely reflected the similarity between the perturbed and unperturbed (flux-free) IH singularities. But there also is an important difference between the perturbed and unperturbed singularities in the case of infinite blue shift. This difference is manifested by the following properties which all hold for generically-perturbaed solutions with k ′ 3 > 2k ′ 1 19 (to which we may refer as "strongly-divergent blue shift"): • The quantity e −s (∂R/∂Ṽ ) diverges at the IH limitṼ →Ṽ ′ 0 .
• As a direct consequence of the above the mass function, which was fixed (M=-c) in the flux-free case, diverges at the IH of the genericallyperturbed solution.
The meaning of these last indicators is simple: The divergence of s at the IH may be locally removed by a gauge transformation as mentioned above. This transformation renders the solution continuous at the IH. Yet generically the flux-carrying solution fails to be smooth at the IH. This lack of smoothness is manifested by the divergence of certain gauge-invariant quantities, so it cannot be cured by any gauge transformation. In the case of stronglydivergent blue shift (k ′ 3 > 2k ′ 1 ) the solution fails to be C 1 . As it turns out, in the case k ′ 1 < k ′ 3 < 2k ′ 1 the solution is C 1 but fails to be C 2 at the IH. In the case k ′ 3 < k ′ 1 (infinite red-shift) the solution will be C 2 , yet generically smoothness will fail at a certain derivative orders (unless 2k ′ 1 /k ′ 3 is an integer). 20 In both blue-shift cases (k ′ 3 > 2k ′ 1 and k ′ 1 < k ′ 3 < 2k ′ 1 ) the generallyperturbed solution fails to be C 2 at the IH. In General-Relativistic terms this means that the Reimann tensor will diverge there (e.g. as measured by polynomial scalars or by its p.p. components)-even if the mass function remains bounded. On the other hand, in the red-shift case (k ′ 3 < 2k ′ 1 ) the 19 This condition was derived in the perturbed RNDS case in Ref. [22]. 20 The general rule is the following: The conjectured smoothness level of R as a function ofṼ3 is the same as that of the function |Ṽ3| 2k ′ 1 /k ′

.
Riemann tensor will remain bounded, though generically its derivatives will blow up at a certain order.
Our primary concern here is the case of generic perturbations, in which both ingoing and outgoing fluxes are present. However, it is also interesting to consider the case a BH perturbed by an ingoing flux while the outgoing flux vanishes. This case is exactly solvable by means of the Vaidya-like solution (see Appendix). In this case scalar quantities as M and e −s (∂R/∂Ṽ )(∂R/∂Ũ ) remain bounded. However e −s (∂R/∂Ṽ ) generically diverges at the IH if k ′ 3 > 2k ′ 1 . Although this quantity is not a scalar, its divergence is a gaugeinvariant phenomenon as discussed in the previous section. Again this divergence indicates an unbounded curvature (though only in the p.p. sense).
As was mentioned above, the divergence of s at the IH is conjectured to be removable by means of a local gauge transformation, leading to a continuous (though not smooth) solution at the IH. From the GR point of view this means that the metric tensor has a continuous (though not smooth) non-singular limit at the IH. This means that the IH singularity is weak (namely physically non-destructive), as discussed at the Introduction.
We summarize some of the conjectured properties of the genericallyperturbed singularity in the table below.

The quantity
Type of divergence Condition for the divergence no null fluid only influx two null fluids +∞ +∞ +∞

Future research directions
The problem which motivated this paper concerned the internal structure of generically-perturbed spinning vacuum black holes. One would like to prove (or disprove) the generic formation of a null singularity inside the BH, and its weakness. We find it likely that our simple two-dimensional semi-linear system presented here properly mimics these properties of General-Relativistic spinning black holes. From this perspective the task of mathematically exploring the internal structure of realistic spinning black holes may be divided into the following stages: (i) Proving (or disproving) the conjecture presented in the previous section concerning the generic behavior of our semi-linear system; (ii) Extending the analysis to the asymptotically-flat case, namely to functions H(R) with two roots and with large-R asymptotic behavior similar to the Λ = 0 case of Eq. (54); and (iii) Further extending the analysis to the real problem of a generically-perturbed spinning black hole.
Although stages (i) and (ii) are not easy, stage (iii) will probably pose a much harder challenge, as it involves the transition from two to four dimensions (and with a much larger number of unknown functions). Nevertheless it may be hoped that the insight gained from our simple toy system in stages (i,ii) may make this challenge a bit easier.
Differentiation of Eq. (82) with respect to u, substitution of the evolution equation for R and then a second differentiation with respect to v, one finds that the second evolution equation (13) is satisfied as well.
Thus, the original system (13) of evolution equations has been reduced to a single ordinary differential equation (83) along lines of constant v. Solving this ODE (with a v-dependent initial condition R 0 (v) at a certain initial u value u = u 0 ) yields the function R(u, v). Subsequently s is obtained by Eq. (82). 21 Substituting this solution in Eq. (31) one obtains The outflux is then given by Eq. (25) or (32), (and, recall, Φ = 0).