Dynamical evolution and stability of quantum corrected

Schwarzschild black holes in semiclassical gravity

Ben Kain
Department of Physics, College of the Holy Cross, Worcester, Massachusetts 01610, USA

The Schwarzschild solution describes a classical static black hole in general relativity. When general
relativity is extended by including semiclassical corrections in the form of a renormalized energy-
momentum tensor, the horizon of the Schwarzschild black hole disappears and is replaced by a
wormhole. We study the stability of this quantum corrected static Schwarzschild solution in semi-
classical gravity by using it as the initial data of a dynamical evolution. We find that the quantum
corrected solution is unstable and that the wormhole can expand or collapse when perturbed. In
vacuum, the wormhole expands, but in the presence of even a small amount of classical matter, the
wormhole collapses, forming a horizon and evolving to an evaporating black hole.
arXiv:2507.00109v1 [gr-qc] 30 Jun 2025

I. INTRODUCTION correction to general relativity, which is the expectation

value of a renormalized energy-momentum tensor oper-
The Schwarzschild solution is the prototypical black ator. We use the Polyakov approximation to compute
hole solution in general relativity. It contains a horizon ⟨Tbµν ⟩ [9–11]. This approximation was used to compute
and a central singularity. When semiclassical corrections the quantum corrected static Schwarzschild solutions [1–
are included—which are arguably the simplest quantum 4] and has been used in dynamical evolutions [12–16].
corrections for extending general relativity—the horizon Moreover, the Polyakov approximation allows for black
disappears and is replaced by a wormhole [1–4]. A curva- hole evaporation [12, 13].
ture singularity is located beyond the wormhole throat, From our simulations, we find that the quantum cor-
but it is no longer centrally located. Similar phenomena rected Schwarzschild black hole is unstable and that the
occur for the quantum corrected Reissner-Nordström [5] wormhole can exhibit both expansion and collapse. In
and Einstein-Yang-Mills [6] black holes. the case of collapse, which occurs when the system is
Semiclassical corrections to general relativity are ex- perturbed through the inclusion of classical matter, a
pected to remove the classical horizon for nonextremal horizon and a central singularity form and the system
static black holes [4, 7]. Static solutions are ubiquitous evolves to an evaporating black hole.
in the study of spherically symmetric gravity. However, In dynamically evolving the quantum corrected
it is important to remember that a static solution may be Schwarzschild black hole, we are dynamically evolving
unstable with respect to time-dependent perturbations. a wormhole. As far as we are aware, very few wormholes
If this is the case, then it is unlikely that the static solu- have been studied dynamically. The Ellis-Bronnikov
tion will occur naturally and it certainly cannot occur as wormhole [17, 18] has been dynamically evolved by var-
the final state in an evolving system. Whenever a static ious groups [19–22] and, recently, the Einstein-Dirac-
solution is found, an immediate and important question Maxwell wormhole [23, 24] has been dynamically evolved
is whether or not it is stable. [25, 26]. We find similarities between the evolution of
A common method for determining stability is to per- the quantum corrected Schwarzschild wormhole and the
turb the static fields with time-dependent perturbations Ellis-Bronnikov wormhole. When relevant, we comment
and then to linearize the system of equations with re- on these similarities.
spect to the perturbations. This method can determine In Sec. II, we briefly review the static quantum cor-
the linear stability of the static solution. Alternatively, rected Schwarzschild black hole. In Sec. III, we describe
the nonlinear stability of the static solution can be de- our dynamical model and explain the numerical methods
termined by dynamically evolving the system with the we use to solve it. In Sec. IV, we present our results and
static solution supplying the initial conditions. If the show that the static solution is unstable. We conclude
static solution is unstable, the nonlinear evolution has in Sec. V. In the Appendix, we present tests of our code.
the advantage of also determining the final state to which Throughout we use units such that c = G = ℏ = 1.
the static solution will evolve. In this work, we study the
stability of the quantum corrected Schwarzschild black
hole by dynamically evolving the static solution. II. STATIC SOLUTIONS
In semiclassical gravity [8], the field equations take the
form In this section, we review static quantum corrected
  Schwarzschild black holes in semiclassical gravity. Our
Gµν = 8π Tµν + ⟨Tbµν ⟩ , (1) main goal is to use these static solutions as initial data for
the dynamical simulations we present in Sec. IV. Com-
where Gµν is the Einstein tensor, Tµν is the classical prehensive studies of these static solutions are given in
energy-momentum tensor, and ⟨Tbµν ⟩ is the semiclassical [1–4], to which we refer the reader for additional details.
 The quantum corrected solutions are found numeri- dynamical simulations and we will comment on our in-
cally. In computing solutions, different forms for the terpretation of this divergence in the next section. For
static spherically symmetric metric have been used, but details on the derivation of (5) or (6), see, for example,
a form we find convenient is [6, 27].
One can show that as x → ∞, the renormalized energy-
ds2 = e2σ(x) (−dt2 + dx2 ) + r2 (x)dΩ2 , (2) momentum tensor drops off sufficiently quickly that the
classical Schwarzschild spacetime, as parametrized by
where x is the radial coordinate, σ(x) and r(x) the ADM mass M , is a solution to (6). This allows
parametrize the metric, r(x) is the areal radius, and us to use the classical Schwarzschild spacetime for the
dΩ ≡ dθ2 + sin2 θ dϕ2 . outer boundary values when solving (6). For classical
The semiclassical field equations are given in (1). In Schwarzschild, with the metric in (2), x is the well-known
vacuum, the classical energy-momentum tensor vanishes, tortoise coordinate. To solve for quantum corrected solu-
Tµν = 0. For the renormalized energy-momentum ten- tions, we proceed as follows. We choose a large value for
sor, ⟨Tbµν ⟩, we use the Polyakov approximation and build r, which marks the outer boundary of the static solution,
it from the exact renormalized energy-momentum ten- and a value for M . At the outer boundary, we have then
sor computed in 1 + 1 dimensions [9–11]. As is com-  
mon, and for simplicity, we will present the renormalized 1 2M  r 
σ = ln 1 − , x = r + 2M ln −1 ,
energy-momentum tensor for N massless scalar fields. In 2 r 2M
the Polyakov approximation, the renormalized energy- (7)
momentum tensor in spherically symmetric 3 + 1 dimen- and, upon taking x derivatives,
sions is taken to be
M r − 2M
3P σ′ = , r′ = , (8)
⟨Tbab ⟩ = 2 ⟨Tbab ⟩(2D) , (3) r2 r
r
which is the classical Schwarzschild solution. With these
where a, b are nonangular components, ⟨Tbµν ⟩(2D) is the outer boundary values and given a value for P , we can
renormalized energy-momentum tensor in 1 + 1 dimen- integrate (6) inward to obtain the quantum corrected
sions, and Schwarzschild solution.
In Fig. 1(a), we show the areal radius, r, as a function
N 2 of x for M = 1 and P = 0.1. Since we find a minimum for
P ≡ ℓ , (4)
12π P r, which is marked by the vertical dotted line, we have
a wormhole structure. The value of r at the minimum
where ℓP is the Planck length. In the units we are using,
is the wormhole throat radius, rth . For comparison, we
ℓP = 1. The choice of the multiplicative factor 3P/r2
show the classical solution, with P = 0, as the dashed
ensures that ⟨Tbµν ⟩ is conserved and that ⟨Tbθθ ⟩ = ⟨Tbϕϕ ⟩ = line. In Fig. 1(b), we show eσ as a function of x. Since
0. eσ is nonzero, there is no horizon at the wormhole throat.
For the metric in (2), the renormalized energy- It can be shown that as x → −∞, there exists a null cur-
momentum tensor works out to vature singularity at a finite proper distance [1, 3, 4]. As
P suggested by Fig. 1(b), eσ continues to be nonzero and
⟨Tbtt ⟩ = − (σ ′ 2 − 2σ ′′ ) the classical horizon has completely disappeared in the
8πr2 (5)
P ′2 quantum corrected solution. We refer to this as a worm-
⟨Tbxx ⟩ = − σ hole, although we note that it may be more appropriate
8πr2
to refer to it as a “wormhole structure” since the region
and ⟨Tbθθ ⟩ = ⟨Tbϕϕ ⟩ = 0, where a prime denotes an x x → −∞ is not asymptotically flat. Quantum corrected
derivative. From the semiclassical field equations in (1), solutions with different values of M and P are qualita-
we can then derive the following second-order equations tively similar to Fig. 1.
for the metric functions: Increasing the mass M increases the wormhole throat
radius. Increasing P does as well, but not by as much
−1
as when increasing M . For example, in Fig. 1(a), rth =
 
1 P P
σ ′′ = − 2r′ σ ′ + σ ′ 2 1− 2 2.081, but for M = 5 and P = 0.1, the wormhole throat
r r r
  −1 (6) radius is 10.024 and for M = 5 and P = 1 it is 10.186.
′′ ′ ′ ′2P P When we present dynamical evolutions of static solutions
r = 2r σ + σ 1− 2 .
r r in Sec. IV, we cannot expect the relationship between the
√ mass and wormhole throat radius that exists for the static
Note that these equations are divergent for r → P . solution to be maintained exactly during the evolution.
This divergence follows from the choice of the multiplica- However, we can assume an approximate relationship,
tive factor in (3). This divergence will not play a role in with increasing (decreasing) mass corresponding to in-
this section because
√ the static solutions we find will not creasing (decreasing) radius. This will help us gain phys-
include r ≤ P . This divergence will play a role in our ical insight into how the radius increases or decreases.

2
 5 III. DYNAMICAL SOLUTIONS

4
We turn now to our main goal, which is to dynami-
3 cally evolve the quantum corrected Schwarzschild black
hole. A dynamical evolution can determine if the static
2 solutions discussed in Sec. II are stable with respect to
perturbations. If we find static solutions that are unsta-
1 ble, as will be the case, then a dynamical evolution can
(a)
determine the final state to which the system will evolve.
0 To perform the dynamical evolution we will develop a
dynamical model and will use the static solutions found
10 -1 in Sec. II as initial data.
For the dynamical model, we choose to work in double
null coordinates and we parametrize the metric as
10 -3
ds2 = −e2σ(u,v) dudv + r2 (u, v)dΩ2 . (9)
-5
10 (b) The metric fields σ and r are the same metric fields in
(2), but of course now are no longer static, and the out-
-50 -40 -30 -20 -10 0 10 going null coordinate u and the ingoing null coordinate
v are defined with respect to the temporal and radial
coordinates in (2) in the usual way,
FIG. 1. The solid lines of both plots display a quantum cor-
rected Schwarzschild black hole with M = 1 and P = 0.1. (a) u = t − x, v = t + x. (10)
The areal radius, r, as a function of the radial coordinate, x.
A minimum occurs at xth = −8.139, indicating a wormhole
throat with radius rth = 2.081. The vertical dotted line in- We choose to use double null coordinates for two im-
dicates the position of the wormhole throat and the dashed portant reasons. First, double null coordinates can
line is the classical Schwarzschild solution, which is included straightforwardly incorporate the renormalized energy-
for comparison. (b) The wormhole throat is horizonless since momentum tensor in the Polyakov approximation [12,
eσ is nonzero. It can be shown that the classical horizon has 13]. Second, if a horizon forms, double null coordinates
completely disappeared in the quantum corrected solution. are well suited for computing the spacetime behind the
horizon. We mention also that double null coordinates
have been used to dynamically evolve wormholes [19, 21].
As mentioned in the Introduction, we will find similar-
It will be useful to define some language for referenc- ities between the evolution of the quantum corrected
ing different regions of the quantum corrected spacetime. Schwarzschild black hole and the evolution of the Ellis-
We’ll refer to the region x > xth as “outside” the worm- Bronnikov wormhole [19–22].
hole. This is the region where we reside. The region When numerically evolving a static solution, dis-
x < xth is reached by passing through the wormhole cretization error inherent in any numerical code acts as a
throat. For convenience, we will refer to this region as small perturbation. Additionally, we will include at times
“inside” the wormhole, though more properly it is on the an explicit perturbation in the form of a pulse of classical
other side of the wormhole. matter. For simplicity, we will use a real massless scalar
field, ϕ, for the pulse, described by the Lagrangian
That semiclassical corrections cause the classical hori-
zon of a Schwarzschild black hole to disappear and to be
1
replaced with an asymmetric wormhole was first found in L = − (∇µ ϕ)(∇µ ϕ), (11)
[1] and subsequently confirmed in [2–4]. Similar results 2
have been found for the classical Reissner-Nordström [5]
and Einstein-Yang-Mills [6] black holes. The question we which
p we then minimally couple to gravity, L →
ask is whether these static solutions, and hence worm- − det(gµν ) L.
holes like that shown in Fig. 1, are stable with respect to
time-dependent perturbations? If they are, then it may
be possible for the wormhole to form naturally. On the
A. Equations
other hand, if the static solutions are unstable, then it
is impossible for the wormhole to form, at least as a fi-
nal state. Further, if the static solutions are unstable, to For the renormalized energy-momentum tensor, we
what final state do they evolve? continue to use the Polyakov approximation in (3). For

3
the metric in (9), we have [12–16] which leads to
P r
1 + 4e−2σ r,u r,v ,


⟨Tbuu ⟩ = 2
(σ,uu − σ,u ) m= (18)
4πr2 2
P 2 which gives the total mass inside a sphere of radius
⟨Tbvv ⟩ = (σ,vv − σ,v ) (12)
4πr2 r(u, v).
P The evolution equations and the mass function are in-
⟨Tbuv ⟩ = − σ,uv variant under transformations of the form
4πr2
and ⟨Tbθθ ⟩ = ⟨Tbϕϕ ⟩ = 0, where we use the notation σ,µ = u → ũ = ũ(u)
∂µ σ and similarly for other variables. Details on how v → ṽ = ṽ(v) (19)
these components are computed can be found in [27, 28]. 1 1
The only classical matter in the system is the scalar σ → σ̃ = σ − ln (∂u ũ) − ln (∂v ṽ)
2 2
field described by the Lagrangian in (11). From this La-
grangian and the metric in (9), the equation of motion and r and ϕ unchanged. This is a coordinate gauge trans-
for the scalar field is the evolution equation formation and σ is a coordinate gauge field. In principle,
we can use this gauge transformation to set σ to any-
1 thing we would like. One of the numerical methods we
ϕ,uv = − (r,u ϕ,v + r,v ϕ,u ) (13)
r use for solving the system of equations, which we out-
and the classical energy momentum tensor is line below, will make heavy use of this gauge transforma-
tion. For this numerical method, the uu component of
Tuu = ϕ2,u the renormalized energy-momentum tensor in (12) trans-
Tvv = ϕ2,v forms nontrivially under the gauge transformation. As
(14) a consequence, the top constraint equation in (15) also
Tuv = 0
transforms nontrivially. We do not use this constraint
Tθθ = 2r2 e−2σ ϕ,u ϕ,v equation in solving for our dynamical solutions, but we
do use it for testing our code. Additional details are given
and Tϕϕ = Tθθ sin2 θ.
in the Appendix.
Inserting the renormalized energy-momentum tensor
into the semiclassical field equations in (1), we find two
constraint equations, B. Initial data
 
P 2
r,uu = 2σ,u r,u − 4πr Tuu + (σ,uu − σ,u )
4πr2 The computational domain is a two-dimensional grid of
  (15) (u, v) values in the ranges ui ≤ u ≤ uf and vi ≤ v ≤ vf .
P 2
r,vv = 2σ,v r,v − 4πr Tvv + (σ ,vv − σ ,v ) , We choose to set ui = 0 and vi = 0, but will continue
4πr2 to write ui and vi for completeness. u = ui and v = vi
and two evolution equations, are the initial hypersurfaces, for which we must supply
initial data.
1 
σ,uv = 2 4r,u r,v + e2σ − 8π 2r2 Tuv + e2σ Tθθ


The initial data will be a static quantum corrected
4r Schwarzschild black hole, as described in Sec. II, some-
 −1
P times augmented with a pulse of scalar field. The static
× 1− 2 (16)
r solutions are functions of x and contain a null curvature
1
 
P
 singularity at x → −∞. This singularity will not be
r,uv = − 4r,u r,v + e2σ − 16πr2 Tuv − σ,uv . included in the initial data, since we will not include
4r 4πr2
x → −∞, which in our computational domain corre-
Note
√ that the σ evolution equation is divergent for r → sponds to u → ∞. Moreover, since the singularity is null,
P . As with the analogous divergence in the static equa- no additional boundary condition would be necessary on
tions, this divergence follows from the choice of the mul- the null initial hypersurface.
tiplicative factor in (3). This divergence will play a role We are at liberty to choose the value of x at the ori-
in our simulations and we make the not uncommon inter- gin of the computational domain, x0 ≡ x(ui , vi ). With
pretation that this divergence corresponds to the central this choice, the value of x at any point on either initial
singularity, which has been shifted from r = 0 by√semi- hypersurface is given by
classical effects [14–16, 28]. We therefore require P to
1
be small compared to any length scale in the system, such x(u, vi ) = x0 − (u − ui )
as the radius of a black hole horizon. 2 (20)
The Misner-Sharp mass function, m(u, v), is defined 1
x(ui , v) = x0 + (v − vi ),
by 2
2m which we use to determine, in the absence of a scalar field
g µν r,µ r,ν = 1 − , (17) pulse, the values of σ and r on the initial hypersurfaces.
r

4
 We will always place the wormhole throat at the origin the row is equal to zero. Since σ is a coordinate gauge
of the computational domain, so that x0 = xth . This is field, this is perfectly consistent. Eilon and Ori refer to
convenient, since then the u = ui initial hypersurface is this as σ gauge. For our system, the maximum value
outside the wormhole and the v = vi initial hypersurface always occurs at the edge of the computational domain
is inside the wormhole. at v = vf . Making this gauge choice has the effect of
For the pulse of scalar field, we have in mind classical increasing the number of rows near the event horizon.
matter that we might try to fire into the wormhole. We In σ gauge, the u coordinates for the grid points are dif-
will therefore place the pulse on the u = ui hypersur- ferent than the u coordinates in the original gauge, where
face, which is outside the wormhole, and ϕ will be zero in the original gauge the values of σ(u, vi ) are equal to
everywhere along the v = vi hypersurface. A look at the those from the static solution. Our code maintains a
equations in Sec. III A shows that they depend on deriva- uniform grid of u coordinates in σ gauge. The u coordi-
tives of ϕ and not on ϕ itself. As such, we will define the nates in the original gauge, for the same grid points, are
pulse in terms of ϕ,v , nonuniform. In this way, many u coordinates are used
  near horizons. Indeed, for a uniform step size in σ gauge
2 v − v1 equal to ∆ũ = 0.01, the step size in the original gauge can
ϕ,v (ui , v) = A sin π for v1 < v < v2 (21)
v2 − v1 become ∆u ∼ 10−10 or even smaller. The v coordinates
of the grid points will always be uniform.
and zero everywhere else on the initial hypersurfaces, The relationship between the u coordinates in σ gauge
where A is a constant. This form is commonly used be- and in the original gauge follows from the bottom equa-
cause both ϕ,v and ϕ,vv are zero at v = v1 and v = v2 . tion in (19). We have implemented a couple of differ-
When including a pulse, we choose to keep σ(ui , v) ent ways for computing the u coordinates in the original
unchanged from the static solution, which we are at lib- gauge. One way is to integrate (19), giving
erty to do since σ is a coordinate gauge field. σ and ϕ Z Z
are then determined on the initial hypersurfaces and r is e2σ(u,vi ) du = e2σ̃(ũ,vi ) dũ, (22)
determined on the v = vi hypersurface, where it is un-
changed from the static solution. It remains to determine where ũ and σ̃ are the values in σ gauge. The right-hand
r on the u = ui hypersurface. We can do this by solving side is computed as our code is running using the trape-
the bottom constraint equation in (15). In this equation, zoidal rule, which is second-order accurate. The left-hand
σ,v and σ,vv are unchanged from the static solution as side, for a range of u values, is computed beforehand for
are r(ui , vi ) and r,v (ui , vi ). We can solve the constraint the initial data. Given the value of the right-hand side,
equation numerically by integrating outward from v = vi we can determine the value of u that gives the left-hand
along the u = ui hypersurface. side using a standard interpolation method (e.g. a spline).
Another way is to write the integral as
Z
C. Numerical methods u = e2(σ̃−σ) dũ (23)

To dynamically evolve the system, we solve numeri- and then to write down a formal solution using the trape-
cally the evolution equations in (13) and (16) using a zoidal rule. We then search for the value of u which solves
standard second-order predictor-corrector scheme (for a the formal solution using the Newton-Raphson method.
description of the scheme, see [29]). Our code solves for We find that both methods work well. The results pre-
field values at all grid points on a u = constant “row,” sented in this paper make use of the first method.
starting at v = vi +∆v, where ∆v is the step size between With the u coordinates in the original gauge we can
grid points, and ending at v = vf , before moving to the compute the remaining fields on the v = vi initial hyper-
next row. surface for σ gauge. For us, this is just r (since ϕ = 0
This method exhibits second-order convergence and is on the v = vi hypersurface). r is gauge invariant, so the
sufficient for determining the stability of the static solu- value of r we use in σ gauge is the value of r for the
tions. As we will see, apparent horizons form and we will corresponding u coordinate in the original gauge.
be interested in accurately computing an apparent hori-
zon out to large values of v. Doing so requires improved
numerical accuracy. We use Eilon and Ori’s adaptive IV. RESULTS
gauge method, which we find to be highly efficient [29].
In Fig. 2(a), we show the dynamical evolution of a
static solution with M = 1, P = 0.1, and no scalar field
1. Adaptive gauge method pulse. The perturbation in Fig. 2(a) is from discretiza-
tion error alone. The thin gray lines are contour lines
The adaptive gauge method makes use of the coordi- for the areal radius, r. The thick black and blue lines
nate gauge freedom of the system. In each row, the value are apparent horizons, defined by r,u = 0 and r,v = 0,
of σ(u, vi ) is chosen such that the maximal value of σ on respectively.

5
 3

3
2.

2.

2.
2
2.
100 100 2
100
2.

2
2.
75 75 75

50 50 50

1
2.
1
2.

3
3

10

15
3

10

15

20
10

15

20

2.
2.

3
5
2.

3
5
3
5
1

25 25 25

2
2.

2.
2
2.

20
2
2.

25

1
25

25
2.
30

30

30
1
2.
35

35

35
1

40

40

40
2.

0 0 0
0 25 50 75 100 0 25 50 75 100 0 25 50 75 100

FIG. 2. Dynamical evolution of a static quantum corrected Schwarzschild black hole with M = 1 and P = 0.1. The gray lines
are contours for the areal radius, r. The thick black line is an apparent horizon defined by r,u = 0 and the thick blue line
is an apparent horizon defined by r,v = 0. There is no scalar field pulse included and the perturbation is from discretization
error alone. We find that the wormhole is expanding. Each plot uses the same initial data, but evolves the system using a
uniform grid with grid spacing ∆u = ∆v = 1/N , where (a) N = 100, (b) 200, and (c) 400. As N increases, the strength of
the perturbation decreases and the static solution holds its configuration longer. This is the expected behavior for an unstable
static solution.

We recall that the wormhole throat of the static solu- we find that the wormhole throat takes longer before ex-
tion is located at the origin, that the region to the right of panding. In other words, if we decrease the strength of
the origin is outside the wormhole and is where we reside, the perturbation, the static solution is able to hold its
and the region above the origin is inside the wormhole. configuration longer. This is precisely the expected be-
The region above contains relatively few contour lines. havior for an unstable static solution. Indeed, the same
We can understand why by looking at Fig. 1(a), where behavior was seen for the Ellis-Bronnikov wormhole [19],
we see that the solid blue curve passes through a rela- which is known to be unstable [30]. We conclude that
tively small range of r values as we move left from the this static quantum corrected Schwarzschild black hole
wormhole throat. solution is unstable.
Moving away from the lower left corner of Fig. 2(a), We now introduce a scalar field pulse as an explicit
we find apparent horizons along the wormhole throat. perturbation. As previously mentioned, we have in mind
This is expected, since from Fig. 1(a) we can see that that we are firing into the wormhole some classical mat-
the wormhole throat in the static solution is defined by ter. We make use of a simple measure for the strength
∂x r = 0. From (10), the wormhole throat in the static of the pulse, defined as follows. We compute the mass of
solution is then also defined by r,u = r,v = 0. the system along the u = ui hypersurface, m(ui , v), using
At around u ≈ 50, the two apparent horizons sepa- (18). At large v, the mass approaches a constant value
rate. It is at this point that the system begins evolving which, in general, is larger than the ADM mass M used
away from the static solution. In between the apparent for the static solution (and is equal to M in the absence
horizons, timelike and null directions necessarily find in- of a pulse). We use the percent increase of the mass at
creasing values of the areal radius: the wormhole throat is large v, with respect to M , as a measure of the strength
expanding. Given that the system is evolving away from of the pulse.
the static solution, it is not surprising that we find ex- In Fig. 3, we show results for a static solution with
pansion, since the only energy-momentum in the system M = 1, P = 0.1, and a pulse with parameters A =
is from the renormalized energy-momentum tensor. For 0.0035, v1 = 10, and v2 = 30. This pulse increases the
comparison, the Ellis-Bronnikov wormhole can also ex- mass by roughly 1%. As in Fig. 2, the gray lines are
hibit expansion and the analogous diagram looks similar contour lines for r and the thick black and blues lines
to Fig. 2(a) (cf. figure 7(a) in [21]). are apparent horizons. Starting in the lower left corner,
Since the only perturbation is from discretization error, we find similar behavior when compared to Fig. 2, in
we can decrease the size of the perturbation by decreas- that the apparent horizons lie along the wormhole throat
ing the spacing between grid points. Figure 2(a) is made and the static solution is holding its configuration. At
with a uniform grid with grid spacing ∆u = ∆v = 1/N around u ≈ 17, the evolution moves away from the static
and N = 100. Figures 2(b) and 2(c) are evolutions with configuration. In between the apparent horizons, timelike
the same initial data as Fig. 2(a), but with N = 200 and and null directions now find decreasing values of the areal
400 respectively. As the discretization error decreases, radius: the wormhole throat is collapsing.

6
 2.082
2.2

2.1
2.08
40

1.5
1.85
2.078
2

2
2.076

r
0 5 10 15 20 25
1.8
30

1.6

(a)
20 1.4
0 50 100 150 200
v
1
2.

3
2.

2.13
3

10
5

(b)
10

15
2.12

20

25
0 2.11
0 25 50 75

r
2.1

FIG. 3. Dynamical evolution of a static quantum corrected 2.09

Schwarzschild black hole with M = 1, P = 0.1, and a scalar
field pulse as an explicit perturbation. The gray lines are 2.08
contours for r and the thick black and blue lines are apparent 0 50 100 150 200
horizons, just as in Fig. 2. The thick red line marks the central v
singularity and the region beyond the red line is not part of
the spacetime. The inclusion of classical matter in the form FIG. 4. (a) The same apparent horizon (defined by r,v = 0)
of a scalar field pulse causes the wormhole to collapse and to shown as the blue curve in Fig. 3, but computed out to large
form event horizons and a central singularity. v. (b) The apparent horizon computed using the same ini-
tial data as Fig. 3, but with evolution equations with P = 0.
This is an inconsistent evolution and should not be taken too
seriously, but it shows that with P = 0 the apparent horizon
The contour lines and the apparent horizons show the
asymptotically approaches a constant value indicating the sys-
formation of event horizons, one at around u ≈ 17.5 and tem evolves to a static black hole. This strongly suggests that
another at around v ≈ 34. We therefore find that, as the the radius of the apparent horizon in (a) is decreasing at large
wormhole collapses, black holes form over the mouths of v because of black hole evaporation.
the wormhole.
√ Eventually the wormhole throat collapses
to rth ≈ P , at which point we reach the central sin-
gularity, which we have indicated with the red line. The The first reason is that wormholes, which violate the null
region beyond the red line is not part of the spacetime. energy condition, can have apparent horizons with de-
For comparison, the Ellis-Bronnikov wormhole also ex- creasing radii when collapsing. Such behavior occurs for
hibits collapse and the analogous diagrams look similar the Ellis-Bronnikov wormhole [20]. The second reason is
to Fig. 3 (cf. figures 2(a) and 5(a) in [21]). that the renormalized energy-momentum tensor allows
To study the region outside the wormhole, we focus for black hole evaporation, during which the radius of
on the apparent horizon defined by r,v = 0, which is the the apparent horizon can decrease [12, 13].
thick blue line in Fig. 3. From (18), the mass inside an In the absence of the renormalized energy-momentum
apparent horizon is always equal to r/2, where r is the tensor, our expectation is that the outside of a col-
radius of the horizon. Since the apparent horizon aligns lapsing wormhole will eventually settle down to a static
with the event horizon at large v, the mass inside the Schwarzschild black hole. In this case, the radius of the
event horizon at large v is equal to r/2. apparent horizon would asymptotically approach a con-
To compute the apparent horizon out to large values stant value, which would be the radius of the black hole.
of v, we use the adaptive gauge method, as described in This is precisely what was found with the Ellis-Bronnikov
Sec. III C 1. In Fig. 4(a), we show the radius of the appar- wormhole [20]. With the inclusion of the renormalized
ent horizon as a function of v. In the inset, we zoom into energy-momentum tensor, this cannot happen because
the small v region. We can see that the apparent horizon the black hole will evaporate and therefore cannot settle
starts out at r = rth = 2.081. It holds its configuration down to Schwarzschild.
until around v ≈ 15 and then starts decreasing. It con- To gain some insight, we consider the following incon-
tinues to decrease out to large v. Typically, the radius sistent evolution: We use the same initial data as used
of an apparent horizon cannot decrease. There are two in Figs. 3 and 4(a), which are based on a static solu-
reasons at play here, which allow the radius to decrease. tion with P = 0.1 and an initial pulse computed using

7
P = 0.1. We then evolve the system, but we set P = 0 and v2 = 20. This range of v values puts the initial pulse
in the evolution equations. Since the evolution equations close to the interesting region where the wormhole throat
and the initial data are inconsistent with one another, we is collapsing, which ends up leading to somewhat chaotic
should not take the results too seriously. Nevertheless, values for energy-momentum tensor in this region. Re-
as shown in Fig. 4(b), the radius of the apparent horizon sults are simplified for a pulse with v1 = 5, v2 = 10,
asymptotically approaches a constant value at large v. and amplitude A = 0.0005. The contour diagram for r
This strongly suggests that it is black hole evaporation, is shown in Fig. 6(a). We can see that it is very simi-
as caused by the renormalized energy-momentum tensor, lar to Fig. 3, except that it takes longer for the worm-
which accounts for the decreasing radius at large v in hole to begin collapsing because the pulse is weaker. In
Fig. 4(a). Figs. 6(b) and 6(c), we show the uu and vv components
The collapse shown in Fig. 3 occurs for a scalar field of the energy-momentum tensor (the apparent horizons
pulse that increases the mass by approximately 1%. We from Fig. 6(a) are shown in yellow). A horizon forms
continue to find that the system collapses as we lower outside the wormhole at u ≈ 29. Outside this horizon
the amplitude of the pulse down to A = 5 × 10−5 , which (u > 29) and at large v, which approaches future null
corresponds to a 0.001% increase in the mass. In general, infinity, the uu component in Fig. 6(b) is positive, in-
we find that a relatively small amount of classical matter dicating outgoing energy-momentum, and the vv com-
is needed to trigger collapse. ponent in Fig. 6(c) is negative, also indicating outgoing
As mentioned in Sec. II, increasing the mass M in- energy-momentum. This net effect of outgoing energy-
creases the wormhole throat radius of the static solution. momentum is consistent with decreasing mass and worm-
As the wormhole throat radius increases, it requires more hole collapse.
computational resources to dynamically evolve the sys- For a wormhole to exist, the energy-momentum tensor
tem because the size of the computational grid must in- around the wormhole throat must violate the null energy
crease. We have dynamically evolved static solutions up condition [31]. In double null coordinates, the null energy
to M = 5 and found that, in the absence of a scalar condition is violated if Tuu +⟨Tbuu ⟩ < 0 or Tvv +⟨Tbvv ⟩ < 0.
field pulse, they are unstable and the wormhole expands. If the null energy condition is satisfied and the wormhole
For M = 5 and P = 0.1, using a uniform grid with is collapsing, we should expect focusing of null geodesics
∆u = ∆v = 1/100, the apparent horizons separate and and the formation of caustics, as follows from the Ray-
the wormhole begins expanding at roughly u ≈ 300. chaudhuri equation and the focusing theorem [32]. This
Aside from this, the resulting diagram looks similar to strongly suggests the formation of a singularity. Indeed,
Fig. 2. We also continue to find for static solutions with from Figs. 6(b) and 6(c) we can see that the null energy
larger masses that a relatively small amount of classical condition flips from being violated to being satisfied right
matter triggers collapse and that the resulting diagrams about where the singularity forms at u ≈ 29 and v ≈ 45.
look similar to Fig. 3.
We can gain some insight as to why the wormhole in
Fig. 2 expands and why the wormhole in Fig. 3 collapses V. CONCLUSION
by looking at components of the energy-momentum ten-
sor. For the evolution shown in Fig. 2(a), contour di- The classical static spherically symmetric vacuum so-
agrams for the uu and vv components of the energy- lution is the Schwarzschild black hole. When extended
momentum tensor are shown in Figs. 5(a) and 5(b) to include semiclassical corrections in the form of a
(the apparent horizons from Fig. 2(a) are shown in yel- renormalized energy-momentum tensor, the horizon dis-
low). The uu component describes outgoing energy- appears and is replaced by a wormhole [1–4]. Since
momentum and the vv component describes ingoing the renormalized energy-momentum tensor can describe
energy-momentum. As can be seen in Figs. 5(a) and black hole evaporation [12, 13], it is perhaps not surpris-
5(b), both are negative. We therefore show in Fig. 5(c) ing that one does not find a static black hole solution,
their difference, ⟨Tbvv ⟩ − ⟨Tbuu ⟩. Since this difference since an evaporating black hole is not static.
is non-negative, we have a net flow of ingoing energy- We have studied the stability of the quantum corrected
momentum. This is consistent with the mass increasing, static vacuum solution by using it as the initial data of
which is easily seen to be the case along the apparent a dynamical evolution. We have shown that the static
horizons in Fig. 2(a), since they move along increasing solution is unstable and that the wormhole will expand
radii. From the relationship between mass and the worm- or collapse.
hole throat radius we previously reviewed for static so- In the absence of classical matter, the wormhole ex-
lutions, the net flow of ingoing energy-momentum is also pands since the only energy-momentum in the system is
consistent with the wormhole expanding. from the renormalized energy-momentum tensor. On the
We could show analogous plots for the collapsing other hand, if there is even a small amount of classical
wormhole in Fig. 3. However, it is easier to see what matter present, our results indicate that the wormhole
is happening if we consider a slightly different evolution. collapses, that it forms a horizon and a central singu-
The evolution shown in Fig. 3 makes use of an initial larity, and that it evolves to an evaporating black hole.
pulse that is nonzero for v1 < v < v2 , where v1 = 10 Since only a small amount of classical matter is necessary

8
 100 100 100

2e-06
1e-09

5e-06
0 5
e-
.5
-2

75 75 75

2e-06
-2.5e-05

1e-07
5
-1.5e-0
05

-06
e-

-2 .8e-05
-5e
.8

5
-2

-1 .5e-0
-5 -05
06
06

07

08
50 50 50

-2
7

.5e
5e-0

e-
e-

e-

e-
06

-1

-1

-1
- 07
- 1e -1e-

1e-08
-08
-1e
25 25 25
25 50 75 100 25 50 75 100 25 50 75 100

FIG. 5. Contour diagrams for (a) the uu component and (b) the vv component of the energy-momentum tensor for the
dynamical evolution shown in Fig. 2(a). The yellow curves are the apparent horizons shown in Fig. 2(a). The difference
between the energy-momentum tensor components is displayed in (c). The region on the left side in (c) is very close to zero,
with the precise values being difficult to compute numerically.

50 50 50

0.0001
05

5
e-0
2.1

1.5
1.85

6e-

40
2

76
-2.7

-2.
-3.3e-05

40 40

-3.2e
30
e-05

-05
-2.88

20 30 30
05
2.1

5e-06
2.3

e- 06
76 e-
3

.
5

-1
5

-2
e-0
10

10
-06
15

76

20 20
-5e
-2.

1e-
07
0
0 25 50 25 50 75 25 50

FIG. 6. Dynamical evolution of a static quantum corrected Schwarzschild black hole with M = 1, P = 0.1, and a scalar field
pulse with v1 = 5, v2 = 10, and A = 0.0005. (a) Contour diagram for the areal radius, r. The black and blue lines are apparent
horizons r,u = 0 and r,v = 0, respectively, and the red line is the central singularity. (b) Contour diagram for the uu component
of the energy-momentum tensor. The apparent horizons are included for convenience as the yellow lines. (c) Contour diagram
for the vv component of the energy-momentum tensor. Some of the contours are unlabeled: Starting at the upper left and
moving to the right, the contour values are -2.76e-5, -2.88e-5, -3.2e-5, -2.88e-5, -2.76e-5, -1e-6, and +1e-4.

to trigger collapse, it appears unlikely for the expand- evolution deserves further study.
ing wormhole to form naturally. Instead, we expect an
evaporating black hole to be the astrophysically relevant There may be some challenges in studying this region.
system. We previously mentioned that only a small range of the
areal radius, r, is probed for a relatively large range
In our study of the collapsing dynamical solution, we of the outgoing null coordinate, u. If we would like to
purposely focused on the region outside the wormhole, evolve further into this region, this will require increased
which is where we reside. However, there is also the computational resources. However, it may be possible
region reached by passing through the wormhole. This is to compress this region using a coordinate gauge trans-
properly the region on the other side of the wormhole, but formation. Additionally, the further we move into this
for convenience we have referred to this region as inside region, the smaller the metric component eσ becomes, as
the wormhole. This is the upper left region of Fig. 3 can be seen from Fig. 1(b), which may cause numerical
and it is in this region that the static solution contains challenges. In terms of the apparent horizon that is in-
a null curvature singularity at u → ∞, whose dynamical side the wormhole (the black curve in Fig. 3 defined by

9
 To present a test of convergence, we focus on the re-
sults shown in Fig. 3, which do not make use of the adap-
Cr 10-5 tive gauge method. Using the same initial data, we have
computed the dynamical solution using three different
10-7 uniform grids, defined by grid spacings ∆u = ∆v = 1/N
with N = 100, 200, and 400. Using these results, we
(a) compute the convergence function
10-9
0 25 50 75 X
v CfN1 ,N2 = |fiN1 − fiN2 |, (A1)
i
10-3
where fiN is the value of field f computed at grid point
10 -5
i using grid spacing N . In (A1), fiN1 and fiN2 must be
Cr

evaluated at the same grid points. In Fig. 7(a), we show

10-7 results for f = r where the sum is performed over rows
(b) with constant u. The dashed black curves are for Cr100,200
10-9 and the dashed yellow curves are for 4Cr200,400 . The lower
0 10 20 30 40 set of curves is for u = 10 and the upper set of curves
u is for u = 30. Since Cr200,400 is multiplied by 4 and the
curves overlap, we have second-order convergence [33].
FIG. 7. The convergence function in (A1) is plotted for f = This continues to be the case for the upper curves as
r and with grid spacing ∆u = ∆v = 1/N and N = 100, they approach the singularity, which is why the upper
200, and 400 where the sum is performed over (a) rows and curves veer upward. Figure 7(b) is an analogous plot,
(b) columns. These curves indicate our code is second-order
but summing over columns with constant v. The lower
convergent. See the explanation given in the Appendix for
details. set of curves is for v = 25 and the upper set of curves
is for v = 60. We find similar results indicating second-
order convergence when using f = σ and f = ϕ.
r,u = 0), we expect that it can be studied similarly to Our code primarily uses the evolution equations in (13)
how we studied the apparent horizon that is outside the and (16). We have therefore confirmed that our dynam-
wormhole. We would need to write code that solves the ical results satisfy both constraint equations in (15), as
evolution equations along “columns” instead of rows so required. When not using the adaptive gauge method,
that we could use the adaptive gauge method applied to we use the constraint equations as written in (15). When
the v coordinate up until the r,u = 0 apparent horizon is using the adaptive gauge method, the upper constraint
computed. equation in (15) takes on a different form. For the coor-
Another potential direction of study is to dynamically dinate transformation u → ũ = ũ(u), defined by
evolve the quantum corrected Reissner-Nordström [5] or
Einstein-Yang-Mills [6] black holes. We expect nonex- ∂u ũ = e2(σ−σ̃) , (A2)
tremal black holes of these types to also be unstable, but
this remains to be confirmed. The quantum corrected the renormalized energy-momentum tensor component
Reissner-Nordström black hole has an extremal solution ⟨Tbuu ⟩ in (12) undergoes a nontrivial transformation, so
and it would be interesting to determine how this static that the constraint equation becomes
solution evolves.  
P 2
r,ũũ = 2σ̃,ũ r,ũ −4πr Tuu + (σ,ũũ + σ,ũ − 2σ,ũ σ̃,ũ ) ,
4πr2
ACKNOWLEDGMENTS (A3)
where Tuu = ϕ2,ũ and where ũ and σ̃ are the values in σ
gauge.
I am thankful to the referee, whose suggestions led to
For an example of our code satisfying the constraint
an improved paper.
equations, we focus on the results in Fig. 4(a), which are
computed using the adaptive gauge method. In comput-
ing these results, we end the computation just after the
Appendix A: Code tests apparent horizon is computed. As such, the computation
does not reach the singularity. We solve the constraint
As a first basic test of our code, we have confirmed equation in (A3) for r using second-order Runge-Kutta
that if we use the static classical Schwarzschild solution and second-order finite differencing for the derivatives.
for initial data, our code reproduces the static solution Note that all derivatives are in terms of ũ, which have
dynamically. This is the case when using the adaptive uniform spacing, making the finite differencing straight-
gauge method and when not using it. forward. At each grid point, we compute the difference

10
 between the result obtained from the constraint equa-
10 -5
tion and the result from the dynamical evolution. We
then compute the root-mean-square (rms) value,
10 -6

10 -7

10 -8
rms(rcon − rdyn ), (A4)
0 50 100 150 200

FIG. 8. Equation (A4), where rcon is computed from the

constraint equation in (A3) and rdyn is from the dynamical
evolution. From top to bottom, the curves are computed us- along columns of constant v. The result is shown in Fig. 8
ing grid spacing ∆ũ = ∆v = 1/N and N = 50, 100, and 200. for grid spacing ∆ũ = ∆v = 1/N and N = 50, 100, and
That the curves drop by a factor of 4 when the grid spac- 200. That the curves drop by a factor of 4 when the
ing drops by a factor of 2 indicates second-order convergence. grid spacing drops by a factor of 2 indicates second-order
That the results are small indicates that the constraint equa- convergence. That the results are small indicates that
tion in (A3) is satisfied by the dynamical evolution. the constraint equation is satisfied.

[1] A. Fabbri, S. Farese, J. Navarro-Salas, G. J. Olmo, arXiv:hep-th/9405007.

and H. Sanchis-Alepuz, Semiclassical zero-temperature [13] S. Ayal and T. Piran, Spherical collapse of a massless
corrections to Schwarzschild spacetime and holography, scalar field with semiclassical corrections, Phys. Rev. D
Phys. Rev. D 73, 104023 (2006), arXiv:hep-th/0512167. 56, 4768 (1997), arXiv:gr-qc/9704027.
[2] P.-M. Ho and Y. Matsuo, Static Black Holes With Back [14] E. Sorkin and T. Piran, Formation and evaporation of
Reaction From Vacuum Energy, Class. Quant. Grav. 35, charged black holes, Phys. Rev. D 63, 124024 (2001),
065012 (2018), arXiv:1703.08662 [hep-th]. arXiv:gr-qc/0103090.
[3] P.-M. Ho and Y. Matsuo, Static Black Hole and Vacuum [15] S. E. Hong, D.-i. Hwang, E. D. Stewart, and D.-h. Yeom,
Energy: Thin Shell and Incompressible Fluid, JHEP 03, The Causal structure of dynamical charged black holes,
096, arXiv:1710.10390 [hep-th]. Class. Quant. Grav. 27, 045014 (2010), arXiv:0808.1709
[4] J. Arrechea, C. Barceló, R. Carballo-Rubio, and L. J. [gr-qc].
Garay, Schwarzschild geometry counterpart in semi- [16] D.-i. Hwang and D.-h. Yeom, Internal structure of
classical gravity, Phys. Rev. D 101, 064059 (2020), charged black holes, Phys. Rev. D 84, 064020 (2011),
arXiv:1911.03213 [gr-qc]. arXiv:1010.2585 [gr-qc].
[5] J. Arrechea, C. Barceló, R. Carballo-Rubio, and L. J. [17] H. G. Ellis, Ether flow through a drainhole - a particle
Garay, Reissner–nordström geometry counterpart in model in general relativity, J. Math. Phys. 14, 104 (1973).
semiclassical gravity, Class. Quant. Grav. 38, 115014 [18] K. A. Bronnikov, Scalar-tensor theory and scalar charge,
(2021), arXiv:2102.03544 [gr-qc]. Acta Phys. Polon. B 4, 251 (1973).
[6] B. Kain, Quantum corrected Einstein-Yang-Mills black [19] H.-a. Shinkai and S. A. Hayward, Fate of the first
holes in semiclassical gravity, Phys. Rev. D 111, 044033 traversible wormhole: Black hole collapse or inflationary
(2025), arXiv:2502.08278 [gr-qc]. expansion, Phys. Rev. D 66, 044005 (2002), arXiv:gr-
[7] C. Berthiere, D. Sarkar, and S. N. Solodukhin, The fate qc/0205041.
of black hole horizons in semiclassical gravity, Phys. Lett. [20] J. A. González, F. S. Guzmán, and O. Sarbach, Insta-
B 786, 21 (2018), arXiv:1712.09914 [hep-th]. bility of wormholes supported by a ghost scalar field.
[8] N. D. Birrell and P. C. W. Davies, Quantum Fields in II. Nonlinear evolution, Class. Quant. Grav. 26, 015011
Curved Space (Cambridge Univ. Press, Cambrdige, UK, (2009), arXiv:0806.1370 [gr-qc].
1984). [21] A. Doroshkevich, J. Hansen, I. Novikov, and A. Shatskiy,
[9] P. C. W. Davies, S. A. Fulling, and W. G. Unruh, En- Passage of radiation through wormholes, Int. J. Mod.
ergy Momentum Tensor Near an Evaporating Black Hole, Phys. D 18, 1665 (2009), arXiv:0812.0702 [gr-qc].
Phys. Rev. D 13, 2720 (1976). [22] K. Calhoun, B. Fay, and B. Kain, Matter traveling
[10] P. Davies and S. A. Fulling, Quantum vacuum energy in through a wormhole, Phys. Rev. D 106, 104054 (2022),
two dimensional space-times, in Proc R Soc London Ser arXiv:2210.04905 [gr-qc].
A, Vol. 354 (1977) pp. 59–77, 1676th ed. [23] J. L. Blázquez-Salcedo, C. Knoll, and E. Radu,
[11] A. M. Polyakov, Quantum Geometry of Bosonic Strings, Traversable wormholes in Einstein-Dirac-Maxwell theory,
Phys. Lett. B 103, 207 (1981). Phys. Rev. Lett. 126, 101102 (2021), arXiv:2010.07317
[12] R. Parentani and T. Piran, The Internal geometry of an [gr-qc].
evaporating black hole, Phys. Rev. Lett. 73, 2805 (1994), [24] R. A. Konoplya and A. Zhidenko, Traversable Worm-

11
 holes in General Relativity, Phys. Rev. Lett. 128, 091104 hole spacetimes, Phys. Rev. D 93, 024016 (2016),
(2022), arXiv:2106.05034 [gr-qc]. arXiv:1510.05273 [gr-qc].
[25] B. Kain, Probing the Connection between Entangled Par- [30] J. A. González, F. S. Guzmán, and O. Sarbach, Insta-
ticles and Wormholes in General Relativity, Phys. Rev. bility of wormholes supported by a ghost scalar field. I.
Lett. 131, 101001 (2023), arXiv:2309.03314 [hep-th]. Linear stability analysis, Class. Quant. Grav. 26, 015010
[26] B. Kain, Are Einstein-Dirac-Maxwell wormholes (2009), arXiv:0806.0608 [gr-qc].
traversable?, Phys. Rev. D 108, 044019 (2023), [31] M. S. Morris and K. S. Thorne, Wormholes in space-time
arXiv:2305.11217 [gr-qc]. and their use for interstellar travel: A tool for teaching
[27] C. Barcelo, R. Carballo, and L. J. Garay, Two for- general relativity, Am. J. Phys. 56, 395 (1988).
malisms, one renormalized stress-energy tensor, Phys. [32] S. W. Hawking and G. F. R. Ellis, The Large Scale Struc-
Rev. D 85, 084001 (2012), arXiv:1112.0489 [gr-qc]. ture of Space-Time, Cambridge Monographs on Math-
[28] A. Fabbri and J. Navarro-Salas, Modeling black hole evap- ematical Physics (Cambridge University Press, Cam-
oration (Imperial College Press, London and World Sci- bridge, UK, 1975).
entific Publishing, Singapore, 2005). [33] M. Alcubierre, Introduction to 3+1 numerical relativity
[29] E. Eilon and A. Ori, Adaptive gauge method for (Oxford University Press, Oxford, UK, 2008).
long-time double-null simulations of spherical black-

12