{'fraction_non_alphanumeric': 0.029909255898366606, 'fraction_numerical': 0.009219600725952813, 'mean_word_length': 4.360116731517509, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical evolution of unstable self-gravitating scalar solitons

Recently, static and spherically symmetric configurations of globally regular self-gravitating scalar solitons were found. These configurations are unstable with respect to radial linear perturbations. In this paper we study the dynamical evolution of such configurations and show that, depending on the sign of the initial perturbation, the solitons either collapse to a Schwarzschild black hole or else ``explode'' into an outward moving domain wall.

which corresponds to an "exploding" configuration where a global phase transition is triggered through the formation of an outward moving domain wall.
This paper is organized as follows. In Section II we derive the system of evolution equations and constraints. Section III describes the initial conditions corresponding to a static soliton with a small perturbation. We discuss the slicing conditions we use for our simulations in Section IV. Section V describes our numerical techniques. In Section VI we show the numerical results for the two types of perturbations mentioned above. We conclude in Section VII. In the Appendix we discuss the hyperbolicity properties of our system of evolution equation.

II. FIELD EQUATIONS
We will consider a model of a scalar field minimally coupled to gravity and with a non-trivial self-interaction potential. The model is described by the Lagrangian (we will use units such that G = c = 1): We choose the following asymmetric scalar-field potential leading to the desired asymptotically flat solutions: with σ, η 1 , η 2 and a constant parameters. For this class of potentials one can easily show that, if η 1 > 2η 2 > 0, φ = a corresponds to a local minimum, φ = a + η 1 to the global minimum and φ = a+η 2 to a local maximum. The key feature of this potential for the asymptotically flat and static solutions to exist is that the local minimum at φ = a is also a zero of V (φ) [1]. The factor σ in front of the potential fixes the scale, so one can always take σ = 1 and just re-scale everything for a different σ afterward. For the simulations discussed here we will take the following values for the parameters: σ = 1, η 1 = 0.5, η 2 = 0.1 and a = 0 (see Figure 1). The field equations following from the Lagrangian (2.1) are the Einstein's field equations and the the Klein-Gordon (KG) equation: 3) The stress-energy tensor for the scalar field is In order to perform a numerical analysis of the problem at hand we shall use a 3+1 approach based on the standard ADM equations [8,9]. Moreover, we shall assume that the shift vanishes. The evolution equations for the 3-metric (γ ij ) and the extrinsic curvature (K ij ) are 6) and the Hamiltonian and momentum constraints are with α the lapse function, D i and R ij the covariant derivative and Ricci tensor associated with γ ij , R := trR ij , K := trK ij , and where the matter sources are defined in terms of the stress-energy tensor as with n µ the unit normal to the spatial hypersurfaces.
We shall focus on the dynamics of a spherically symmetric space-time described by the metric ds 2 = −α 2 dt 2 + Adr 2 + Br 2 dθ 2 + sin 2 θdϕ 2 . (2.13) The spherical symmetry implies that all dynamical functions depend only on r and t.
To write down our evolution equations we start by defining the quantities (2.16) We will work with the mixed components of the extrinsic curvature K A := K r r , K B := K θ θ = K φ φ , and with the matter variables J A := J r , M A := M r r and We also introduce the extra variables and use them instead of D A and K A . The evolution equations for {A, B, D, D B , K, K B } can be obtained directly from the ADM equations. However, in order to have a hyperbolic evolution system (see Appendix) we will remove the terms proportional to ∂ r D B and ∂ r K B from the ADM evolution equations for K and D, respectively, using the constraints.
Spherical coordinates can be problematic at the origin. In order to regularize the coordinate singularity we use the regularization scheme described in Ref. [10] and introduce the auxiliary variable which will be promoted to an independent dynamical quantity and evolved explicitly in time.
Since we have used the momentum constraint to modify the evolution equation for D, for regularizing the equations we need to replace this variable with The final set of dynamical variables is then and their (regularized) evolution equations are These equations can be easily shown to form a strongly hyperbolic system (see Appendix).
In terms of the new variables, the Hamiltonian and momentum constraints take the form Finally, the KG equation can be written as where we have introduced the variables For numerical purposes, the evolution equation for Π above is further transformed into the equivalent form: The last term on the right hand side of this equation includes a first derivative with respect to r 3 . The reason behind this transformation is related to the regularization near the origin of the 1/r 2 factor in the original equation.

III. INITIAL CONDITIONS
The initial conditions used to study the evolution are the static soliton solutions computed in Ref. [1], plus/minus a small Gaussian perturbation in Π. Strictly speaking, pure static initial conditions would not result in non-trivial evolution. However, since the configurations are unstable, the truncation errors are in fact enough to trigger the evolution. Nonetheless, since the numerical errors become smaller when resolution is increased, it is better to include a finite perturbation as "detonator".

A. Static soliton
As mentioned above, we shall consider small perturbations to the the static soliton found in Ref. [1]. The static soliton is obtained by first taking and then solving the set of equations The value φ 0 is used as a shooting parameter such that the scalar-field settles asymptotically to the local minimum V = 0 of the potential, thus guaranteeing asymptotic flatness. The value α 0 , on the other hand, is completely arbitrary. The reason for this is that α appears only in two of the equations (in one of them through the combination D α = ∂ r α/α) and both these equations are in fact linear in α. In practice one takes α 0 = 1, solves the system of equations, and later re-scales the lapse to make sure that its asymptotic value is 1. The numerical code we use to solve for the initial data is fourth order accurate, in contrast with our evolution code which is only second order accurate (see Sect. VI).
For the particular values of the parameters in the potential used here, the static soliton solution is shown in Figures 2 to 4. In all the plots the numerical grid extends to r = 100 and the radial coordinate is shown on a logarithmic scale. The value we find for φ 0 is (3.9) The uncertainty in the eleventh significative figure is estimated by taking the difference from the results of our two finest grids with ∆r = 0.025 and ∆r = 0.0125. It is also important to mention that, having fixed the parameters of the potential, this static solution is unique, as opposed to the case of the corresponding hairy black holes where the family of static solutions is parametrized by the horizon radius [1].

B. Perturbed soliton
In order to perturb the static soliton we add a small Gaussian to the time derivative of the scalar field but still take B = 1, D B = K = 0 (notice that it is now inconsistent to ask for K B = 0 as well). We also keep the conditions ∂ t Π = ∂ t K = 0. The new initial-data set of equations to be solved is 14) where the last equation arises from the fact that the momentum constraint is no longer trivial. Notice also that we now have a second order equation for the lapse α. This is because we no longer have the condition ∂ t K B = 0, but rather the condition ∂ t K = 0. In the static case both conditions are equivalent, so one can solve the simpler first order equation corresponding to ∂ t K B = 0, but in the perturbed case this is no longer true. The data obtained in this way will be consistent with the constraints but will not be static, and will reduce to the static soliton solution when ǫ = 0. For the cases discussed here, we have taken the Gaussian parameters to be ǫ = ±0.002 and s = 10.0. We do not show plots for the initial values for φ, α and A, as they look very similar to those of the static soliton.

C. Mass
To find the mass of the configurations we notice that, since the scalar field decays very rapidly with r, in the asymptotic region we will have essentially the Schwarzschild metric. If we re-parametrize the radial metric A as then the mass of the configuration can be obtained from In practice, there is no need to go very far, since in the region where the scalar field is zero this expression converges very rapidly to a limiting value. The mass function m(r) for the static configuration is shown in Figure 5. Its limiting value for large r turns out to be where again the error bar is estimated by comparing results from runs with ∆r = 0.025 and ∆r = 0.0125. For the static configuration one can show, using the Hamiltonian constraint, that the mass given by (3.17) can also be computed as an integral of the energy density ρ associated with the scalar field: However, this is only true for static initial data, as extrinsic curvature terms will enter this expression for nonstatic configurations. Using the integral expression we find the following value for the mass of the static soliton (we have only integrated up to the boundary of the grid, but since ρ goes to zero rapidly, the value of the integral stops changing to the last decimal figure much before we reach the boundary) The error in this value is considerably larger, which is not surprising since the numerical integral was done using the second order trapezium rule. Still, it agrees with the result above to 7 significant figures. For the perturbed configurations, the value of M computed with (3.17) turns out to be extremely close to the value for the static soliton, differing by about 0.1%: Notice that the mass is the same for both perturbations, which is not surprising as ρ depends only on Π 2 .

IV. SLICING CONDITIONS
We have assumed so far a vanishing shift. We now need to specify a slicing condition, that is, a way to find the value of the lapse function α during the evolution. A common choice is the so-called polar slicing, which in the case of vanishing shift is equivalent to the choice of areal (or radial) coordinates throughout the evolution. This consists of imposing B(t, r) = 1, which in turn implies K B (t, r) = 0, and hence K = K A . From the evolution equation for K B one then finds an ordinary differential equation in r for the lapse that can be solved at each time step. The main drawback of this approach is that the radial-polar slicing gauge does not penetrate BH horizons, since inside a horizon it is impossible to keep the areas of spheres fixed without a non-trivial shift vector.
Since we will be interested in looking for BH horizons during the evolution, we have decided to use horizon penetrating coordinates instead of the above gauge. In fact, we will use two different types of slicing conditions, namely, harmonic slicing and maximal slicing (see below), each adapted to the physical situation that results from a given type of perturbation. The adequate choice, of course, is only known a posteriori, so in practice we have performed short trial runs in each case with har-monic slicing in order to gain some insight about the dynamic behavior of the system.
Harmonic slicing is a well known condition that relates the lapse function to the spatial volume elements: α = f (x i ) γ 1/2 , with γ the determinant of the spatial metric and f (x i ) an arbitrary time-independent function. This slicing condition is in fact equivalent to the requirement that the time coordinate t satisfies the wave equation ✷t = 0, i.e. t is a harmonic function. The harmonic condition can also be written as an evolution equation for the lapse in the form which is a particular case of the more general Bona-Masso family of slicing conditions [11]. It can be shown that harmonic slicing avoids so-called "focusing singularities" [12,13] (those for which the spatial volume elements vanish at a bounded rate): the lapse collapses to zero at the same rate as the volume elements. The singularity avoidance of harmonic slicing, however, is only marginal in the sense that the slices come arbitrarily close to the singularity after a finite time. For this reason harmonic slicing is usually not used for black hole evolutions where one wishes to avoid the singularity.
For black hole evolutions, choices different than harmonic are usually better. In particular, maximal slicing, defined by K = ∂ t K = 0, has been a standard workhorse for evolving black holes because of its strong singularity avoidance properties (see e.g. [14]). Maximal slicing leads to the following elliptic equation for the lapse In our case, we have chosen initial data that has precisely the property that K = ∂ t K = 0 (even in the perturbed case). Moreover, as we will show below, one particular perturbation of the scalar soliton leads to a BH formation, so maximal slicing is a natural choice. In spherical symmetry, the maximal slicing condition reduces to the following second order ordinary differential equation for the lapse function α This equation is solved by imposing the following boundary conditions with r b the position of the outer boundary. The first condition is just the regularity condition at the origin, and the second is a Robin outer boundary condition that guarantees that the lapse behaves asymptotically as 1 + k/r (with k constant). Now, although maximal slicing is well adapted for the case when the scalar field collapses to a BH, for the "exploding" configuration (i.e. an outward moving domain wall), the harmonic slicing condition is much better suited (4.1). This is because, as shown in Sec. VI below, in that case the spacetime inside the outward moving wall behaves in a manner similar to an anti-de-Sitter spacetime, as the scalar field moves toward the true minimum of the potential. This behavior will produce a big-crunch type singularity in a finite proper time, to which maximal slicing responds by making the lapse collapse extremely rapidly everywhere, thus completely freezing the evolution [15,16]. Harmonic slicing, on the other hand, makes the lapse collapse much slower, allowing a much longer evolution. As mentioned above, such slicing allows the hypersurfaces to move arbitrarily close to the singularity.

V. NUMERICAL METHODOLOGY
For the time integration in our code we use an iterative Crank-Nicholson scheme with 3 iterations (see e.g. [17]). Derivatives are represented by second order centered finite differences on the radial grid. The numerical evolution is therefore expected to be second order accurate in both ∆r and ∆t. We also typically take ∆t = ∆r/4 in order to be sure that we satisfy the Courant-Friedrichs-Levy stability condition.

A. Boundary conditions
From the evolution equations (2.22), (2.23) and (2.35) it is clear that the metric functions A and B and the scalar field φ evolve only through source terms and can be updated point-wise all the way to the boundary of the numerical grid. The evolution equations for other variables, however, have spatial derivatives in the right hand side so we require a boundary condition. For these variables we use an outgoing wave (Sommerfeld) boundary condition. That is, we assume that near the boundary all dynamical variables behave as spherical waves: where v is the wave speed. In practice, we do not use the boundary condition (5.1) as above, but rather we use it in differential form We also assume that the boundary is sufficiently far from the origin for the wave speed to be essentially unit (α and A are close to 1). So we actually use We use finite differences for the above equation (consistent to second order in both space and time) and apply it to {λ, U, D B , K, K B , ξ, Π}. When using harmonic slicing, we apply the condition also to D α and evolve α pointwise all the way to the boundary. It is important to mention that the boundary conditions just described are in fact not compatible with the constraints. We therefore expect a small amount of constraint violation to be introduced at the boundaries, which will later propagate inward. We have checked that this has essentially no impact on the results presented here by moving the outer boundary to different locations and comparing the results. Constraint preserving boundary conditions are of course possible to implement, and are certainly desirable, but they are not necessary for the results discussed here.

B. Apparent Horizon
Throughout the evolution we look for apparent horizons as an indicator of the formation of a black hole. Of course, an apparent horizon will only coincide with an event horizon if the spacetime reaches a stationary state at late times. But if an apparent horizon is present, an event horizon is guaranteed to exist outside of it as long as the cosmic censor conjecture and the null energy condition hold [18,19]. Apparent horizons have the advantage that their definition is local in the sense that they can be located within a given spatial hypersurface. Event horizons, on the other hand, are defined globally and can therefore only be located once the whole evolution of the spacetime is known (or at least, once the evolution is known up to the point where the spacetime is essentially stationary). In our case, the null energy condition is independent of V (φ), and will therefore hold since for a null vector ℓ µ we will have T µν ℓ µ ℓ ν = (ℓ µ ∇ µ φ) 2 ≥ 0. Moreover, as it turns out, when an apparent horizon forms in our simulations the scalar field is in the region where the potential is positive, so an event horizon is guaranteed to exist out of it.
An apparent horizon (AH) is defined as the outermost marginally trapped surface [18], that is, a closed twodimensional surface for which the expansion of the outgoing null geodesics is zero. In terms of 3+1 quantities, the expansion H of a congruence of null rays moving in the outward normal direction to a given surface takes the form [20] with s i the unit outward pointing normal vector to the surface. An AH is then the outer-most closed surface such that H = 0 everywhere on the surface. In spherical symmetry, this expression can be shown to reduce to the simple form In the code we evaluate this expression over the whole numerical grid looking for places where it becomes negative (notice that for Minkowski we simply have H = 2/r > 0). If H ever changes sign, the outermost place where this happens is identified as the AH.
Once an apparent horizon has been located at r = r AH , its area will be given by: where B AH is the metric function B evaluated at r = r AH . We can also associate a mass to this AH through the formula

VI. RESULTS FROM NUMERICAL SIMULATIONS
As already mentioned, the perturbations corresponding to ǫ = 0.002 and ǫ = −0.002 result in very different evolutions, so we will consider each case in turn.

A. Collapse to a black hole
We will first consider the perturbed configuration corresponding to ǫ = −0.002. In this case we will use maximal slicing to determine the lapse. We have done numerical simulations using four different grid resolutions, ∆r = (0.1, 0.05, 0.025, 0.0125). In all cases the numerical grid extended to r = 100. Figures 6 and 7 show the evolution of the metric functions A and B. In all the plots solid lines correspond to initial and final configurations, and dotted lines to intermediate stages (the separation in time between these lines is ∆t = 25). At first, there are small oscillations around the initial value, but later on the radial metric starts to grow in a form characteristic of the slice stretching associated with BH spacetimes. Figure 8 shows the behavior of the lapse function. Again, we see the characteristic collapse of the lapse indicative of the approach to a singularity.
The corresponding evolution of the scalar field can be seen in Figure 9. Notice how the scalar field moves toward the local minimum at φ = 0 everywhere. At late times, the evolution of the inner regions is frozen due to the collapse of the lapse there. By that time, however, the scalar field has values below 0.1 everywhere, which implies that we are in the region where the potential is positive.
The slice stretching and collapse of the lapse are indicative of the presence of a BH, but they are not enough to prove that one is present, so we have looked for the appearance of an AH. Figure 10 shows the mass associated with the AH's found during the simulation. From the figure one can see that an AH first appears at t ∼ 115. The mass associated with this horizon grows slightly at first as scalar field falls into the BH, but later settles. For consistency, this mass should always remain below the initial ADM mass of the spacetime. From the figure we see that for low resolution the numerical error causes the horizon mass to eventually become larger than the ADM mass. At higher resolutions, however, the horizon mass remains below the ADM mass throughout the entire evolution (see inset of Fig. 10). The small difference between the horizon mass and the ADM mass indicates that a small amount of scalar field has been radiated away.
A crucial test of the validity of our code is the behavior of the constraints. Analytically, the constraints must be identically zero, but numerical errors mean that for our simulations the constraints have in fact finite values. Still, those values should converge to zero as resolution is increased. Figure 11 shows the logarithm of the root As expected, the Hamiltonian constraint is converging to zero to second order in ∆r.

B. Explosion
Next we consider the perturbed configuration corresponding to ǫ = +0.002. In this case we will use harmonic slicing to determine the lapse. We will use the same grid resolutions as before, but will extend the numerical grid further out to r = 500. Figures 12 and 13 show the evolution of the metric functions A and B for this case (lines are now separated in time by ∆t = 75). It is evident that the dynamical evolution is completely different to the case described in the previous section. In the first place, there is no indication of the slice stretching effect. Moreover, in the evolution of the radial metric A it is clear that there is a wall moving outward. The wall moves essentially at a uniform speed, even if this is not evident in the log plot (this speed is approximately 1 in our units, which coincides with the speed of light in the outer regions). Inside this wall the radial metric is collapsing to zero. The angular metric is also collapsing to zero in this region, but not as rapidly. The evolution of the lapse function is shown in Fig. 14. Again, the presence of an outward moving wall is evident. Inside the wall the lapse is collapsing to zero, indicating that we are approaching a singularity. However, in this case no apparent horizon was located for the duration of the run. The evolution of the scalar field can be seen in Fig. 15. In contrast to the results of the previous section, in this case the scalar field is moving toward the true minimum of the potential at φ = 0.5. Since this minimum corresponds to a negative value of the potential, the interior of the wall resembles an anti-de-Sitter spacetime, except for the fact that the scalar field is not uniform. Still, one would expect the formation of a big-crunch type singularity in this region in a finite proper time [15,16]. However, because of the singularity avoiding properties of harmonic slicing, this singularity would only be reached after an infinite coordinate time.
Finally, Figure 16 shows the evolution of the logarithm of the root mean square of the Hamiltonian constraint for the same four resolutions. Second order convergence is FIG. 14: Same as Figure 12 for the lapse function α. again evident.
The simulation seems to indicate that the perturbation has triggered a global phase transition from the false vacuum in the exterior to the true vacuum in the interior. This phase transition propagates through a domain wall that always moves outward, which implies that the spacetime will never reach a stationary state.

VII. DISCUSSION
We have considered the dynamical evolution of static self-gravitating scalar solitons. These soliton configurations arise for self-interaction potentials V (φ) that have a local minimum for which V = 0, and a global minimum such that V < 0. This allows for static, asymptotically flat solutions for which the scalar field interpolates between the two minima. The static configurations were shown to be unstable [1], and we have studied their re- sponse to two different types of perturbations differing in the sign of the perturbing term.
In one case the configuration is found to undergo gravitational collapse, as the scalar field moves toward the local minimum of the potential everywhere. A small amount of scalar field is radiated away, and the spacetime finally settles to a stationary Schwarzschild BH.
For the other type of perturbation the scalar field "explodes", triggering a global phase transition that propagates through an outward moving domain wall. This wall separates an inner region where the scalar field moves toward the true minimum of the potential (the true vacuum), and an outer region with the scalar field in the local minimum. As the true minimum has V < 0, the inner region behaves in a manner similar to anti-de-Sitter spacetime, and should produce a big-crunch type singularity in a finite proper time. The spacetime in this case does not reach a stationary state, as the domain wall always keeps moving outwards. A question arises as to whether the singularity that is forming in this case is naked, since there is no evidence for an apparent horizon [15,16]. We believe that this is unlikely, as nothing special seems to happen during the evolution. More likely an event horizon does exist, but this is no BH in the standard sense since, first, no trapped surfaces form inside it, and second, a stationary state is never reached and the BH eventually swallows the whole spacetime.

APPENDIX
Recently, the need for having a hyperbolic system of evolution equations has been stressed by several authors as a necessary condition for the well posedness of the Cauchy initial data problem (and also for the stability of the numerical evolution). In order to comply with this requirement we have re-written the ADM equations together with the harmonic slicing condition as a first order system of the form ∂ t w + Å r ∂ r w = S , and S source terms that include no derivatives (the variables {α, A, B} evolve only through source terms and need not be considered for the hyperbolicity analysis).
The eigenvalues of the matrix Å r turn out to be Clearly, all eigenvalues are real. Also, the eigenvector matrix has determinant equal to B 2 /2A 3 > 0, showing the linear independence of the eigenvectors. The system of equations is therefore strongly hyperbolic.
The above analysis applies to the case of harmonic slicing. For maximal slicing the argument is similar, except that in that case K is no longer a dynamical variable (K = 0), and the lapse is obtained from an elliptic equation. The reduced system of equations for {D B , K B , U, λ} is still strongly hyperbolic, and the full system now has both an elliptic and a hyperbolic sector.