{'fraction_non_alphanumeric': 0.04268308456424451, 'fraction_numerical': 0.016840417000801924, 'mean_word_length': 4.040156453715776, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the instability of charged wormholes supported by a ghost scalar field

In previous work, we analyzed the linear and nonlinear stability of static, spherically symmetric wormhole solutions to Einstein's field equations coupled to a massless ghost scalar field. Our analysis revealed that all these solutions are unstable with respect to linear and nonlinear spherically symmetric perturbations and showed that the perturbation causes the wormholes to either decay to a Schwarzschild black hole or undergo a rapid expansion. Here, we consider charged generalization of the previous models by adding to the gravitational and ghost scalar field an electromagnetic one. We first derive the most general static, spherically symmetric wormholes in this theory and show that they give rise to a four-parameter family of solutions. This family can be naturally divided into subcritical, critical and supercritical solutions depending on the sign of the sum of the asymptotic masses. Then, we analyze the linear stability of these solutions. We prove that all subcritical and all critical solutions possess one exponentially in time growing mode. It follows that all subcritical and critical wormholes are linearly unstable. In the supercritical case we provide numerical evidence for the existence of a similar unstable mode.


I. INTRODUCTION
Wormhole spacetimes in Einstein's gravitational theory have received considerable attention in the literature. Presumably, this is due to their interesting topological and causal properties which open the door to spectacular phenomena such as interstellar travel and time machines [1,2,3]. More recently, wormholes have also been proposed as black hole mimickers [4] (see also [5]). However, there are several problems which make these solutions somehow unattractive from a physical point of view and pushes them to the verge of science fiction. First of all, wormholes need to be supported by exotic matter if they are to be considered as asymptotically flat, globally hyperbolic spacetime solutions of Einstein's field equations [1,6]. This means that they require an energy-stress tensor which violates the (averaged) null energy condition, a phenomena which has not been observed for classical matter fields in the laboratory. Furthermore, to our knowledge, there are no static, asymptotically flat wormhole models which have been shown to be linearly stable with respect to arbitrary linear fluctuations of the metric and matter fields in the given model. This is to be contrasted with static, asymptotically flat black holes with a regular horizon in vacuum or electrovacuum spacetimes, which are known to be stable with respect to linear perturbations [7,8,9,10,11,12]. Therefore, even if some exotic form of matter could be found in the Universe, it is not clear whether or not it could be used to form static wormholes.
In previous work [13,14] we analyzed the question of wormhole stability for a very simple matter model which consists of a massless ghost scalar field, that is, a massless scalar field whose kinetic energy has a reversed sign. We found that all static and spherically symmetric wormholes in this theory are unstable with respect to linear and nonlinear perturbations. Each of these wormholes possesses a single unstable mode which causes the wormhole to collapse to a black hole or to undergo a rapid expansion. Furthermore, the time scale associated to the linear instability is of the order of the areal radius of the wormhole's throat divided by the speed of light.
The purpose of this article is to analyze whether or not one could stabilize these wormholes. One possible mechanism for stabilization is to consider stationary, rotating generalization of such wormhole solutions and hope that they become stable if the angular momentum is large enough. Slowly rotating wormholes have been constructed in [15] by studying linear perturbations of the static solutions. However, such slowly rotating solutions cannot be expected to be stable since the unperturbed solutions are unstable. Wormhole solutions of the full nonlinear field equations which represent rotating generalizations of the static ones have been considered in [16,17], but a detailed stability analysis of such solutions is not expected to be simple. For this reason we consider, here, a different possible stabilization mechanism which is the addition of an electromagnetic charge. Therefore, we add a Maxwell field to the ghost field and consider static, spherically symmetric wormhole solutions in this theory. It turns out that the resulting field equations can be integrated exactly and give rise to a four-parameter family of wormhole solutions. This is described in section II. Some of these solutions have been found in [17]. The phase space of solutions can naturally be divided into subcritical, critical and supercritical wormholes depending on whether the sum of their asymptotic Arnowitt-Deser-Misner (ADM) masses is negative, zero or positive, respectively.
Next, we analyze the stability of such wormholes with respect to linear, spherically symmetric fluctuations. In section III we derive the perturbation equations and cast them into a constrained wave system for two gauge-invariant quantities. We then discuss different ways of decoupling this system. In section IV we first prove that all subcritical and critical wormholes are linearly unstable. We do this using the same techniques as in our previous work, namely the theory of Schrödinger operators on the real line. Then, we focus our attention to supercritical wormholes which turn out to be more interesting. In this case, we are not able to reduce the perturbation equations to a single wave equation with regular potential. Therefore, we analyze the numerical stability of these wormholes by numerical integration of the constrained wave system. The results indicate that these wormholes are linearly unstable as well, with the perturbations exhibiting an exponential growth of the form e βτ for large values of proper time τ at the throat, where β > 0. Interestingly, we find that if the charge is large enough, this exponential growth is accompanied with an oscillating factor of the form cos(ωτ − δ) for some frequency ω and phase δ. Our results also indicate that the growth rate β decreases monotonically when the charge increases, and for the range of parameters used in our numerical simulations we observe that β can be decreased by a factor of more than 100 when compared to the uncharged case. Finally, we show by a numerical matching algorithm that these asymptotic solutions correspond to eigenfunctions of the spatial perturbation operator corresponding to a complex eigenvalue −(β + iω) 2 . A summary of our results and conclusions are given in section V. Technical properties of the static wormhole solutions which are needed for the stability analysis are stated and proved in an appendix.

II. STATIC, SPHERICALLY SYMMETRIC CHARGED WORMHOLES
We consider a spherically symmetric gravitational field which is coupled to a massless ghost scalar field Φ and to an electromagnetic field F . Therefore, choosing suitable local coordinates t, x, ϑ, ϕ, the metric has the form where the functions d = d(t, x), a = a(t, x) and c = c(t, x) depend only on the time coordinate t and the spatial coordinate x. We assume that Φ and F are also spherically symmetric which means that Φ = Φ(t, x) and that F has the form with two functions α = α(t, x) and β = β(t, x). We are interested in traversable wormhole geometries which consist of a throat connecting two asymptotically flat ends at x → +∞ and x → −∞, respectively. This means that the areal radius r = e c is strictly positive and proportional to |x| for large |x| and that the 2-manifold (M ,g) = (R 2 , −e 2d dt 2 + e 2a dx 2 ) is regular and asymptotically flat at x → ±∞.
The equations of motion are where R µν and ∇ µ denote, respectively, the Ricci tensor and the covariant derivative associated with the spacetime metric g µν . In terms of Newton's constant G, the coupling constants κ 0 and κ are given by κ 0 = 2G > 0 and κ = −8πG < 0 which is negative due to the fact that Φ describes a ghost scalar field. For the spherically symmetric ansatz (1,2) Maxwell's equations (4) imply that α = Q e e a+d−2c and β = Q m with Q e and Q m two constants representing, respectively, the electric and magnetic charge. Setting Q := Q 2 e + Q 2 m the remaining Eqs. (3,5) yield the evolution equations which are subject to the constraints Here, the subscript t and x refer to the derivatives with respect to t and x, respectively. For a static configuration, the scalar field Φ and the metric coefficients d, a and c are independent of t. In this case the field equations can be integrated analytically. For zero charge the corresponding solutions have been obtained in [18,19]. Here, we generalize their solutions for the charged case. In the static case, the field equations simplify to Adopting a gauge where a = −d, the first two equations imply that [e 2(c+d) ] xx = 2 which has the general solution e 2(c+d) = x 2 + 2α 1 x + α 0 with two constants α 0 and α 1 . By a suitable translation of the coordinate x it is always possible to obtain α 1 = 0. Furthermore, since we are interested in wormhole geometries with the properties described below equation (1) we need e 2(c+d) > 0 for all x ∈ R. Therefore, we have e 2(c+d) = x 2 + b 2 with some strictly positive constant b > 0. Equation (13) then gives Φ = Φ 1 arctan(x/b) + Φ 0 with two integration constants Φ 0 and Φ 1 . Since only the gradient of Φ appears in the equations we set the parameter Φ 0 to zero in what follows. Next, setting ǫ := κ 0 Q 2 /2b 2 and y := arctan(x/b) ∈ (−π/2, π/2), Eq. (11) gives The unique local solution with initial conditions d| y=0 = γ 0 ∈ R and d y | y=0 = γ 1 ∈ R is We distinguish between the subcritical case where Λ > 0, the critical case where Λ = 0 and d = γ 0 − log(1 − γ 1 y) and the supercritical case where Λ = iµ for some real, strictly positive number µ in which case d = γ 0 − log cos(µy) − γ 1 sin(µy) µ . In order to obtain a global wormhole solution, we need the expression inside the square brackets to be strictly positive for all y ∈ [−π/2, π/2]. This is the case if and only if Λ |γ 1 | < 1 in the subcritical case, π 2 |γ 1 | < 1 in the critical case, µ < 1 and Finally, Eq. (14) yields the relation between the parameters Λ and Φ 1 . Summarizing, we obtain the solutions where y = arctan(x/b) and e 2d = e 2γ0 cosh(Λy) − γ 1 sinh(Λy) Λ −2 with Λ = γ 2 1 − κ 0 (Q 2 e + Q 2 m )e 2γ0 /(2b 2 ). The parameters b, Φ 1 , Q e , Q m , γ 0 and γ 1 are subject to the two constraints (17) and (18). Notice that the constant with Ω a nonvanishing constant leaves the solution unchanged. In particular, we can rescale the coordinates such that either lim x→+∞ d = 0 or lim x→−∞ d = 0 which shows that the spacetime described by the metric (21) has indeed two asymptotically flat ends at x → +∞ and x → −∞, respectively. Therefore, we obtain a four-parameter family of wormhole solutions characterized by the scale invariant quantities B := be −γ0 > 0, γ 1 ∈ R, Q e ∈ R and Q m ∈ R, which are subject to the restrictions (17). In the particular case Q e = Q m = 0 this family reduces to the bi-parametric solution obtained in [18,19]. As shown in [13,14] these uncharged wormholes are unstable with respect to linear and nonlinear perturbations.
Let us analyze the physical properties of the wormhole solutions. First, the wormhole throat is given by the global minimum of the areal radius r = √ In Lemma 1 in the appendix we show that r has a unique minimum which is determined by the unique root of the function c x given in Eq. (22) below. Next, we compute the Misner-Sharp mass function [20]. For the spherically symmetric spacetime metric given by equation (1) it is defined by Specialized to the static family of solutions described by Eqs. (19,20,21) it yields where c x is the derivative of the logarithm of the areal radius, The ADM masses of the two asymptotically flat ends can be computed by considering the asymptotic values m ±∞ := lim x→±∞ m(x) of the mass function, which yields In particular, we have the relations where we have set ν := κ 0 (Q 2 e + Q 2 m )/(2B 2 ). From the first relation and Eqs. (19,20,21) we see that the wormholes are reflection-symmetric about their throat if and only if γ 1 = 0. Since Λ = γ 2 1 − ν 2 the asymmetry parameter γ 1 and the dimensionless charge ν determine the asymptotic masses m +∞ and m −∞ , the total electromagnetic charge Q := Q 2 e + Q 2 m and the areal radius of the throat up to the scale factor B. From the second and third relations we see that in the subcritical case (Λ > 0) the two masses have opposite signs and that their sum is negative. In the critical case (Λ = 0) the sum of the two masses is zero, and the masses are different from zero unless γ 1 = Q e = Q m = 0 in which case m +∞ = m −∞ = 0. In the supercritical case (Λ = iµ, µ > 0), the sum of the two masses is positive, with opposite signs if γ 2 1 > µ 2 tan 2 (µπ/2) and equal signs if γ 2 1 < µ 2 tan 2 (µπ/2) while one of the masses is zero and the other is positive if γ 1 = ±µ tan(µπ/2).

III. DERIVATION OF THE PULSATION AND MASTER EQUATIONS
In this section we derive the relevant equations for analyzing the linear stability of the four-parameter family of static wormhole solutions discussed in the previous section. For this, we consider small perturbations of the form where Φ is the background solution, and where denotes the variation of Φ. The same applies to the other fields d, a and c. A general method for analyzing such perturbations has been developed in [21]. Since in spherical symmetry there are no gravitational nor electromagnetic dynamical degrees of freedom, one obtains a single master equation for the linearized scalar field δΦ. However, as we have discussed in detail in [13] for the uncharged case, the resulting master equation turns out to be singular at the throat which leads to difficulties when studying the stability by standard methods based on Schrödinger operators. For this reason, we will not use the method described in [21] and instead base our treatment on a gauge-invariant approach which leads to a constrained wave system which is everywhere regular.

A. Gauge-invariant quantities
With respect to an infinitesimal coordinate transformation δt → δt + ξ t , δx → δx + ξ x on the 2-manifoldM generated by a vector field (ξ t , ξ x ), we have Since Φ x = 0 everywhere we may construct the following two gauge-invariant fields, which reduce to δa and δc, respectively, in the gauge δΦ = 0.
B. Constrained wave system for A and C Here we derive a constrained wave system for the gauge-invariant quantities A and C. In order to do so, we first rescale the coordinate x such that b = 1. In terms of the coordinate y = arctan(x) which satisfies ∂ y = e d−a+2c ∂ x , the background Eqs. (11)(12)(13) yield where ǫ 2 = κ 0 Q 2 /2. Next, we consider the constraint equations (9,10). Their linearization yields With the help of the background equations (27) one may rewrite this as which shows that the expression inside the square bracket must be equal to a constant σ. In terms of the gaugeinvariant quantities A and C defined above, the resulting first integral is The interpretation of the constant σ is the following. With respect to an infinitesimal variation of the constants Bγ 1 and Λ (keeping the charges Q e and Q m fixed), the family of static solutions (21) yields the linearized solution where the function F is given by Introducing the expressions (33) in (32) gives σ = −δ(Bγ 1 ). Therefore, σ describes variations of the static family of wormhole solutions with respect to the constant Bγ 1 . Since any solution to the linearized equations may be written as the sum of such a variation plus a solution with δ(Bγ 1 ) = 0 we may assume that σ = 0 in the following. Next, we linearize the evolution Eqs. (6,7,8). For simplicity, we choose the gauge such that δΦ = 0, in which case A = δa and C = δc. Linearization of Eq. (8) yields δd − δa + 2δc = h(t) for some function h(t) which we may set to zero by a redefinition of δt. Using δd = δa − 2δc and the first integral (32) with σ = 0 in the linearization of the evolution Eqs. (6,7) a lengthy calculation yields the constrained wave system where we have defined The constrained wave system (35) describes the dynamics of the two gauge-invariant linearized fields A and C. The linear stability properties of the wormholes are determined by the Cauchy evolution of this system. Two difficulties with analyzing the properties of the solutions are the fact that we are confronted with a coupled system of two equations (as opposed to a single, scalar equation) and the presence of the constraint C = 0. In the next subsection we start by deriving a decoupled equation for the constraint field C based on a factorization method. This shows that it is sufficient to enforce the constraint and its time-derivative at an initial time. As a byproduct of our factorization method, we also obtain a master equation for a quantity v 1 defined below, from which u 1 and u 2 can be reconstructed. However, this equation turns out to be singular at the throat, and as mentioned before, this means one has to be careful with the stability analysis. For this reason, we derive in the following subsection a different master equation which, in the critical and subcritical cases is everywhere regular and allows to prove that such wormholes are linearly unstable. In the supercritical case, however, both master equations turn out to be singular, and so the constrained wave system (35) has to be analyzed directly. If the potential V would be symmetric, or if it could be brought into symmetric form by a linear transformation of u, one could analyze the system by spectral analysis of the formally self-adjoint operator H = −e −2c ∂ y e −2c ∂ y + V . The transformation w 1 = u 1 , w 2 = −u 1 + u 2 brings the system into the form with D = 3c 2 y − 2e 2(c+d) + 4ǫ 2 e 2d and E = 4c y d y − e 2(c+d) + ǫ 2 e 2d + 4Λ 2 . If Λ 2 > 0 this can be symmetrized by a trivial rescaling of w; for Λ 2 < 0, however,V 12 andV 21 have different signs. In fact, we will show in the next section that in the latter case the operator H may have complex eigenvalues, and so it cannot be written as a symmetric operator. Here, we try to factorize it in the form H = AB, with the two first-order operators where K is a 2 × 2 matrix which has to satisfy the Riccati matrix equation If K solves (37), the factorization H = AB allows us to rewrite the wave problem u tt + Hu = 0 into first-order form, In particular, it follows that v satisfies the dual wave problem where BA = −∂ 2 +V with the transformed potentialV = −V + 2K 2 . A solution to (37) can be obtained as follows. We demand that the second component of v t = −Bu is proportional to the constraint variable C defined in Eq. (35). This implies that K must be of the form with two unkown functions f and g. (Notice that in this case ∂ t v 2 = e −2c C.) Introducing this ansatz into Eq. (37) yields the unique solution for f and g. The transformed potential is then Therefore, the constraint variable v 2 satisfies a decoupled wave equation with potentialV 22 = e −4c (c yy − c 2 y ), and it is consistent to enforce the constraint C = 0. Setting v 2 = 0 the first-order system (38,39) then reduces to In particular, v 1 satisfies the following master equation, As mentioned above, the resulting potentialV 11 is singular at the wormhole throat, where c y = 0. As discussed in detail in [13] this enforces an unphysical boundary condition at the throat if one tries to define −∂ 2 +V 11 as a selfadjoint operator. Namely, it requires v 1 to approach zero sufficiently rapidly as y converges to the throat's location. However, there is no reason for enforcing such a strong condition on v 1 from a physical point of view. As we will see in the next section, physically permissible perturbations even allow v 1 to diverge at the throat. In particular, this implies that the operator −∂ 2 +V 11 is not symmetric when defined on the space of physically permissible states. As we show next, it is possible to obtain a different master equation which in the subcritical and critical cases is everywhere regular and yields a self-adjoint operator without enforcing unphysical boundary conditions.

D. Master equation II
The derivation of the new master equation is based on the observation that the constrained wave system (35) possesses the particular solution which is obtained from (33) after setting δ(Bγ 1 ) = 0 and δΛ = −Λ. A related solution is obtained by multiplication with t, corresponding to the following solution of the first-order system (43,44,45) Since Ψ 0 is a time-independent solution of the master equation (46) we may rewrite the latter in first-order form The corresponding dual wave equation is with the transformed potential where we have used the identity ∂yΨ0 Ψ0 + f = G H c y . In terms of the new variable χ the first-order system yields which allows to obtain the gauge-invariant perturbation quantities u 1 , u 2 (and c y v 1 ) from χ after a time integration. In contrast to the first master equation, the new master equation (50) is regular at the throat. In fact, the potential W is everywhere regular as long as the function F does not have any zeroes. This turns out to be the case for the critical and subcritical cases, see Lemma 2 in the appendix.

IV. LINEAR STABILITY ANALYSIS
In this section we discuss the linear stability of the wormhole solutions in the subcritical, critical and supercritical cases. In the first two cases, we show instability by proving that the master equation (50) which is regular in those cases possesses precisely one exponentially in time growing mode. In the supercritical case, both master equations are singular, and we analyze the stability by different means.

A. The subcritical case
The results from the previous section allow us to describe the stability problem in the subcritical case by the regular master equation (50) on the Hilbert space X = L 2 (R, e −2d dx) which admits the zero mode Since the function c y /e c is uniformly bounded on −∞ < x < +∞ and possesses exactly one zero, and since the function F is strictly positive and satisfies F/x → Λ 2 (1 + Λ 2 ) −1 π/2 for x → ±∞, this mode belongs to X and represents a bound state of the Schrödinger operator −∂ 2 + W . Because it has one node, it follows from the nodal theorem 1 that it is the first excited bound state and so the operator −∂ 2 + W possesses precisely one negative eigenvalue −β 2 < 0 with eigenfunction χ β , corresponding to an exponentially growing mode of (50) which is of the form Since the coefficients in the Eqs. (52,53,54) are everywhere regular this gives rise to a unique unstable mode for each subcritical wormhole. We conclude that all such wormholes are unstable with respect to linear, spherically symmetric perturbations.

B. The critical case
Linear perturbations in the critical case are also described by the regular master equation (50). However, in contrast to the previous case, the function is not normalizable, so the above argument based on the nodal theorem does not directly apply. Instead, we construct a family f n , n = 1, 2, 3, ..., of static solutions to (50) which are defined on the interval −n ≤ x < +∞ and satisfy the following two properties for large enough n: (i) f n (−n) = 0, (ii) (f n ) x (−n) = 0, (iii) f n has exactly one zero on the interval x > n. It then follows from the results in Ref. [23] that there is a unique bound state with negative energy, as in the subcritical case. Therefore, all critical wormholes are linearly unstable as well.
The family f n is defined as follows. Let x throat be the value of x at the throat, where c y = 0, and let n > −x throat . Then, where k is a constant to be determined. It is simple to verify that f n satisfies the relation χ 0 ∂f n = 1 + f n ∂χ 0 and the master equation (50) on the two open intervals (−n, x throat ), (x throat , +∞), and that f n (−n) = 0 and (f n ) x (−n) < 0. In order to analyze the behavior of f n near x = x throat we first introduce the new coordinate in terms of which we have ∂ = ∂ ρ . Next, since χ 0 satisfies (−∂ 2 + W )χ 0 = 0 and vanishes at ρ = 0, we can write it in the form χ 0 (ρ) = αρ[1 + ρ 2 q(ρ)], where q is a smooth function on R such that q(0) = 0, and α > 0. Also, 1 + ρ 2 q(ρ) has to be strictly positive since otherwise χ 0 would have more than one zero. In terms of this, we find which shows that f n is continuous at ρ = 0. Furthermore, we see that (f n ) ρ is also continuous at ρ = 0 provided we choose k such that Therefore, f n can be extended on the whole interval −n < x < +∞. Finally, we notice that f n is negative on the interval −n < x ≤ x throat and positive for large enough x > x throat . Hence, f n has at least one zero at some x * > x throat . On the other hand, because χ 0 (x * )∂f n (x * ) = 1 + f n (x * )∂χ 0 (x * ) = 1 and χ 0 (x * ) > 0 it follows that ∂f n (x * ) > 0 which means that this zero is unique. In the reflection-symmetric case where γ 1 = 0 we have C. The supercritical case: Numerical integration of the constraint wave problem For the supercritical case, F has always two zeroes and both master equations are singular. Therefore, the previous arguments cannot be used to establish the linear instability of the wormholes in this case. Short of an analytic proof, we shall analyze the stability of supercritical wormholes by numerical means. We start in this subsection with a numerical integration scheme of the constraint wave problem (35). We model this scheme on the recent work in Ref. [24], where hyperboloidal time slices with a compactified space coordinate are used. One of the main advantages of this scheme is that the Cauchy evolution is performed on a compactified domain where the boundaries of the domain correspond to future null infinity. In this way, one avoids the problem of introducing an artificial timelike boundary with absorbing boundary conditions. This is particularly attractive for our stability problem since we do not want the time evolution of our solution to be contaminated by artificial boundary conditions.
The starting point of the numerical scheme is to rewrite the wave equation in Eq. (35) in geometric form whereg = −e 2d dt 2 + e −2d dx 2 is the wormhole 2-metric,∇ is the associated connection, andṼ = e −2d V is the rescaled potential. An observation of later importance is the fact that the rescaled potential decays at least as O(x −2 ) when x → ±∞. Next, we introduce a new time coordinate τ := t − h(x), where the height function h, which will be specified later, is such that the τ = const slices are everywhere spacelike and asymptote to outgoing null geodesics as x → ±∞. This means that the function e 4d h 2 x is strictly less than one and converges to one as x → ±∞. The τ = const slices, as embedded in the four-dimensional, spherically symmetric spacetime with 2-metricg and areal radius r have mean curvature where the function J is given by J := e 3d h x / 1 − e 4d h 2 x . Here, we choose the height function such that J = c 0 x(x 2 + 3b 2 )/(x 2 + b 2 ) with some positive constant c 0 . This implies that the mean curvature of the τ = const slices which converges to c 0 > 0 as x → ±∞, and the function e 2d h x also satisfies the properties described above Eq. (59).
With respect to the new coordinates (τ, z) the 2-metric assumes the form Since ΩJ is everywhere regular and positive near z = ±1, the conformal 2-metricĝ is regular for all z ∈ [−1, 1]. Since the two-dimensional wave operator is conformally covariant, the wave equation (58) This system is to be integrated on the compact interval z ∈ [−1, 1]. The characteristic speeds are Therefore, λ − < 0 < λ + on z ∈ (−1, 1). At the left boundary, λ − = 0 while at the right boundary, λ + = 0 which implies that both boundaries are outflow. This is of course expected from the fact that the boundaries z = ±1 represent future null infinity. In terms of the compactified coordinates, the constraint C = 0 reads A simple way of solving the constraints is to specify initial data for u 2 , Π 2 and D 2 which are compactly supported away from the throat, and to use the constraint C = 0 and its time derivative in order to determine the remaining fields u 1 , Π 1 and D 1 , keeping in mind that The procedure we used to study the solution is as follows.
1. We set up initial data by specifying an initial pulse for u 2 and solving the constraint, as described above.
3. We measure the amplitude of u 1 at the location of the throat. 4. We fit the resulting value of u 1 at the throat using the ansatz A cos(ωτ − δ)e βτ , where τ is proper time at the throat, and determine in this way the growth rate β of the perturbation and the frequency of oscillation ω in case there is one.
We use a finite differences approximation method with a method of lines for the evolution of the perturbation (61,62,63). When the dimensionless charge ν is not too close to its limit value given by the constraint (17), secondorder accurate stencils show convergence and confident results. In the case ν approaches its limit value, the results obtained from the second-order accurate stencils fail to show convergence, and we use eighth-order accurate stencils in order to reduce the errors and work in the convergence regime.
In Fig. 1 we present our results for the reflection-symmetric case γ 1 = 0. For small values of ν the perturbation grows exponentially and no oscillating mode shows up. However, there is a threshold of ν around ν * = 0.55 above which the exponential growth is modulated by a harmonic component. Also, we observe that the growth rate decreases to less than 1% of the rate in the uncharged case as ν increases from zero to 0.96. Above ν = 0.96 the functions in the evolution equations become stiff and our numerical approach breaks down. In order to validate our results, we compare them with the results of the matching method described below.
An interesting question that arises from the plots in Fig. 1 is whether or not the growth rate β may be zero for some value of ν lying between 0.965 and 1, implying the existence of charged wormholes which are stable with respect to linear radial perturbations. A different possibility is that β stays positive for all 0 ≤ ν < 1 and converges to zero for ν → 1, meaning that all supercritical wormholes are linearly unstable, but the growth rate can be made arbitrarily small by adding a sufficient amount of charge. The third possibility is that β is bounded away from zero for all 0 ≤ ν < 1 in which case all supercritical wormholes are linearly unstable as well. The answer to this question requires an analytic understanding of the behavior of β as a function of ν near one and lies beyond the scope of this paper. In Fig. 2 we present results for the asymmetric case for different values of γ 1 = 0. The behavior we find is similar to the massless case: for each γ 1 there is a threshold value for ν above which the perturbation shows an oscillatory behavior while growing exponentially.  (17) for the supercritical case. Similarly to the reflection-symmetric case there is a threshold between purely exponential growth and exponential oscillating growth.
the functions u 1 (t, ·) and u 2 (t, ·) being regular for all times t ≥ 0. Defining the quantity Ψ := e −c c y v 1 = e c ∂ t u 2 and using master equation I, we may rewrite the system (43,44,45) in the form where we notice the fact that c yy is everywhere positive, see Lemma 1 in the appendix. Since all coefficients in these equations are regular, it follows that Ψ is regular if and only if u 1 and u 2 are. Therefore, we look for eigenmodes for which Ψ is regular and has time-dependency e (β+iω)t . As a consequence of master equation I, the quantity Ψ satisfies the wave equation We are looking for solutions of this equation with time-dependency e (β+iω)t which are everywhere regular and vanish in the asymptotic regime y → ±π/2. This leads to the following eigenvalue problem, where Γ = β + iω. This equation has a regular singular point at the throat c y = 0. In the reflection-symmetric case where γ 1 = 0 one finds 2c yy c y = 2 y 1 + 2 3 (1 + µ 2 )y 2 + O(y 4 ) , In order to understand the asymptotic behavior of the solutions of Eq. (70) when y → ±π/2 approaches the asymptotic regime, it is easier to go back to the consideration of master equation I, which in terms of the coordinate ρ defined in (57) reads where we have assumed again a time-dependency e Γt for v 1 . In the reflection-symmetric case, we find According to standard theorems 2 , there exist local solutions near ρ = +∞ such that where h is a C 2 -function satisfying h(ρ) → 0, h ρ (ρ) → 0 when ρ → +∞. Therefore, we obtain the asymptotic solution for y → π/2. Our numerical method consists in integrating the local solutions (71,73) and to match their Wronski determinant at some intermediate point 0 < y 1 < π/2 by fine-tuning the complex eigenvalue Γ. The result is shown in Fig. 1, and agrees very well (within 5% of relative error) with the growth rate and oscillation frequency found from the Cauchy evolution.

V. CONCLUSIONS
In this article we analyze the stability of static, spherically symmetric general relativistic wormhole solutions sourced by a massless ghost scalar field and an electromagnetic field. To this purpose we first construct the complete spectrum of such solutions. Among the solutions we distinguish between three types depending on the values of the charge and mass parameters: subcritical, critical and supercritical. We show that in the first two cases all solutions are unstable with respect to linear spherically symmetric perturbations. This is done by reducing the perturbation equation to a Sturm-Liouville problem and using standard tools in Schrödinger operator theory.
In the supercritical case we are not able to reduce the perturbation equations to a scalar equation which is everywhere regular. We instead obtain a constrained wave system for two gauge-invariant quantities. We study this system numerically as a Cauchy evolution problem based on a domain compactification scheme and observe the growth of the perturbations at the wormhole throat. The analysis reveals: i) for small values of the dimensionless charge parameter ν the growth of the perturbation is exponential in proper time, ii) there is a threshold ν * in the value of the charge above which the exponential growth is modulated by an oscillatory component, iii) the growth rate decreases when the value of the charge increases, iv) as the value of the charge approaches its limiting value the functions in the evolution equations become stiff and our numerical approach breaks down. However, we find consistent results up to values of ν = 0.96. Within this regime we can decrease the growth rate by a factor of more than 100 compared to the uncharged case. To validate our numerical results we also analyzed the eigenvalues of the spatial operator in the perturbation equation. Based on a matching method we find a real eigenvalue for ν < ν * . Above that value the eigenvalue acquires a nontrivial imaginary part, explaining the oscillatory behavior found in the Cauchy evolution. The eigenvalues fit the growth rate and the oscillation frequency obtained with the Cauchy evolution within an accuracy of 5%. The fact that we find complex eigenvalues provides an explanation for the inability to reduce the perturbation equations to a Sturm-Liouville problem, for which all eigenvalues are necessarily real.
The motivation for this work was to study a possible mechanism for stabilizing the static, spherically symmetric wormholes solutions supported by a ghost scalar field. We have analyzed the stability behavior when such solutions are charged up by an electromagnetic field. The fact that this mechanism does not seem to be able to stabilize the wormholes raises some doubts about the success of other similar mechanisms, like adding angular momentum. d yy = ǫ 2 e 2d , c yy = 1 + x 2 − ǫ 2 e 2d , (2d y + c y )c y = x 2 − ǫ 2 e 2d − Λ 2 , which follow from Eqs. (11,12,14) with e 2(c+d) = 1 + x 2 , ǫ 2 = κ 0 Q 2 /2 and the relation (18). This together with c y = x − d y also implies the equations which turn out to be useful. The next result is related to the global behavior of the areal radius r = e c .

Lemma 1
The function c has a unique local minimum and c yy is strictly positive on the interval −∞ < x < +∞.
Proof. First, we notice from Eq. (A2) that at points where c y = 0 we must have c yy = 1 + Λ 2 > 0 which shows that local extrema of c are necessarily local minima. Since c → +∞ when x → ±∞ it follows that c has a unique local minimum at some point x = x throat . Next, we prove that c yy is everywhere positive. For this, we first notice that in the uncharged case ǫ = 0, it follows from Eq. (A1) that c yy = 1 + x 2 ≥ 1. Next, assume that ε > 0 and that c yy has a zero at some point x * . Since c yy = 1 + Λ 2 > 0 at x throat , x * lies either to the left or to the right of x throat . Suppose that x * < x throat . Then, from Eq. (A2), it follows that c yyy = 2d yy c y < 0 at x = x * , where we have used d yy = ǫ 2 e 2d > 0 and c y < 0 for points to the left of the throat. This means that c yy crosses the x-axis from above at all its zeroes lying to the left of x throat . However, since c yy > 0 at x throat , this means that c yy cannot have zeroes to the left of x throat . A similar argument shows that c yy cannot have zeroes to the right of x throat either.
Next, we analyze the properties of the function F = 1 + xy − yc y 1 + Λ 2 = 1 + (Λ 2 x + d y )y 1 + Λ 2 defined in (34). Using Eqs. (A2) it is not difficult to show that F satisfies the relation F y = d y F + (1 + xy) Λ 2 1 + Λ 2 c y . (A3) For the linear stability analysis, it is important to know whether or not F has zeroes.
Lemma 2 In the subcritical and critical cases, the function F is strictly positive on the interval −∞ < x < +∞. In the supercritical case, F has precisely two zeroes.
Proof. In the critical case, F = (1 − γ 1 y) −1 which is everywhere positive. In the other cases, since d y is bounded and Λ 2 = 0, we see from its definition that F satisfies F → +∞ (F → −∞) as |x| → +∞ in the subcritical (supercritical) case. On the other hand, at x throat where c y = 0, we have F = 1 + xy > 0. Finally, let x * be a zero of F which lies to the left of x throat . From Eq. (A3) it follows that (1 + Λ 2 )F y = (1 + xy)Λ 2 c y at x * which is negative (positive) in