Stability of Charged Global AdS$_4$ Spacetimes

We study linear and nonlinear stability of asymptotically AdS$_4$ solutions in Einstein-Maxwell-scalar theory. After summarizing the set of static solutions we first examine thermodynamical stability in the grand canonical ensemble and the phase transitions that occur among them. In the second part of the paper we focus on nonlinear stability in the microcanonical ensemble by evolving radial perturbations numerically. We find hints of an instability corner for vanishingly small perturbations of the same kind as the ones present in the uncharged case. Collapses are avoided, instead, if the charge and mass of the perturbations come to close the line of solitons. Finally we examine the soliton solutions. The linear spectrum of normal modes is not resonant and instability turns on at extrema of the mass curve. Linear stability extends to nonlinear stability up to some threshold for the amplitude of the perturbation. Beyond that, the soliton is destroyed and collapses to a hairy black hole. The relative width of this stability band scales down with the charge Q, and does not survive the blow up limit to a planar geometry.


Introduction
Stability of Anti de Sitter vacua has received much attention in the last years since the numerical experiment in [1] . There it was found that nonlinear evolution of some family of arbitrarily small scalar field perturbations inevitably end up in the collapse and formation of a black hole. This is the so called instability corner [2]. Perturbatively, the problem was also examined [3] [4] in the context of purely gravitational perturbations, and the importance of two ingredients was signalled: the presence of a fully resonant spectrum for the linearized perturbations, and the existence of periodic solutions that acted as centers of some stability islands in the space of initial conditions. Back to the scalar field case, this suggestion received further backup from other contributions [5] [6]. In this case, the periodic solutions were named oscillons, and indeed, the two previous observations became consistent in that the spectrum of linearized perturbations around an oscillon turns out not to be fully resonant. After some years of analytic and numerical work, it has become clear that there is a wealth of situations that one can encounter. One may choose to change either the dynamics (the action) or the kinematics (the boundary conditions). Generically, when departing from the easiest case of perturbations on pure AdS, the resonant property is lost [7] [8]. An odd case is that of AdS 3 , where the mass gap (of conical singularities) kills the instability corner despite the fact that the spectrum is fully resonant [9] [10]. In this paper we study both linear and nonlinear stability of asymptotically AdS vacua in Einstein-Maxwell theory interacting with a charged massless scalar. There's little we can add to motivate focussing on this theory. The space of static solutions is rich and involves regular solitons, hairy and Reissner-Nordström (RN) black holes. The first two are interpreted in the context of AdS/CFT as holographic duals of quantum states with spontaneously broken global U (1) symmetry. We will be working in AdS 4 in global coordinates. Our findings build up in parallel with the achievements of [11] [12] where the landscape of static solutions was unraveled for the case of AdS 5 . The thermodynamical (linear) stability issue concerning the competition of these three solutions depends upon the ensemble that is considered. We will construct the grand potential to find the correct vacua in the grand canonical ensemble. For the nonlinear stability we must carry out numerical evolution analysis, searching for endpoints of the evolution in one of the above possible static forms. Performing these simulations with Dirichlet boundary conditions at the boundary, is tantamount to studying the thermalization of the dual quantum system in the microcanonical ensemble. The Lagrangian that governs the dynamics is given by 1 for an Anti de Sitter space with curvature radius l. Here D µ φ = (∂ µ − ieA µ )φ and there is a U (1) gauge symmetry. The coupling e is a free parameter and we can measure all lengths in units of l, hence setting l = 1. Analytic expressions for generic d are provided in the appendices, but the numerical analysis will examine only the d = 3 case. RN and hairy black holes are dual to normal and superfluid phases of the dual field theory respectively. In this last case, it is the scalar field that condenses, breaking spontaneously the U (1) global symmetry. The dual state is called an s-wave superfluid [13]. A similar setup was considered in [14], albeit in the probe limit, while our analysis takes full account of the backreaction of the fields on the geometry. In this context, it is known that black holes come typically in pairs, one of which is small and thermodynamically unstable. Moreover, there are solitonic solutions which are completely regular in the bulk. Whenever these are the ground states, they should correspond in the dual field theory to zero entropy condensates. Finally, we also find extremal black holes, both among RN and hairy solutions. The plan of the paper is the following. Section 2 contains a brief summary of the set of solutions that can be obtained as a function of the electromagnetic coupling e. It is written in the same style as in [11] and [12] for AdS 5 , highlighting differences and similarities. In section 3 we address the issue of the connection of such solutions with superfluid states. This requires identifying the thermodynamically dominant phase in the grand canonical ensemble. We compute the renormalized grand potential and find out both first and second order phase transitions. nonlinear stability is the topic of concern in section 4. It requires numerically evolving arbitrarily small initial perturbations over AdS 4 . The central question in this section is whether similar conclusions as those obtained in [1] can be extrapolated to the present situation. The answer points in the affirmative, namely, the vacuum exhibits a corner of nonlinear instability even at finite charge. On the other hand, the islands of stability get amplified, probably as a consequence of the electrostatic repulsion. This matches with expectations, since the AdS 4 linear scalar perturbations are still fully resonant in this theory. In section 5 we study the stability of solitonic ground states. In a band above them, oscillating solutions never decay. We found the upper bound of this protected from collapse region and examine the way it scales upon blowing up into a black brane geometry.
Our findings indicate that, in such limit, this region does not survive. We also determine the scalar field normal mode spectrum above a solitonic background, finding that it has a dispersive character. In the light of findings in [4] this implies that the weakly turbulent instability present for the AdS 4 vacuum is absent in these geometries.

Static solutions
In this section we describe the space of static solutions as a function of the coupling e. We shall start by providing some background material, and relegate to the appendices all the cumbersome expressions.

Construction
The ansatz for the metric of AdS d+1 in global coordinates is the following The standard Schwarzschild radial coordinate r is related to x as r = tan(x). Concerning the gauge field, isotropy is maintained by selecting A = A t (t, x)dt. This ansatz leaves as residual gauge symmetry under which the equations of motion must be covariant. This motivates defining the following The ansatz (2.1) breaks general covariance leaving a residual time-reparametrization symmetry This can be employed to either set δ b or δ o to zero, fixing the time coordinate t to be the proper time at either x b = π/2 or x o , respectively. Equations of motion for the time dependent ansatz can be found in appendix A, eqs. (A.8)-(A.11) and (A.16) . Equations for static solutions follow after setting Π = −ie δ A t φ/f ,Π = 0 and C = A t . We provide them here for completeness The scalar field φ can be taken real by an appropriate gauge choice. Solving (2.6) requires, as usual, specifying boundary conditions at both ends of the radial domain. At the boundary, which is holographically related to the UV regime of the field theory, the following asymptotics follow directly from the equations of motion, setting ρ = x − π 2 , The holographic dictionary identifies the leading term of the scalar field asymptotic series expansion, φ b , with the source of a marginal field theory operator, O, and the subleading term φ b,3 with its vacuum expectation value, O . Likewise, from the asymptotic expansion of the gauge field we can read off both the chemical potential, µ, and charge density, Q ∝ J t , of the dual field theory state. We refer to the infrared end of the radial coordinate as x o , with either x o = 0 for solitonic solutions, or x o = x h > 0 for the horizon of a black hole. The infrared series expansions read where it should be noted that regularity of solutions at the origin demands f o = 1 and forces every odd term of (2.11)-(2.14) to be zero. In the black hole case, the existence of a horizon and regularity of the gauge field one-form implies Each static solution of the equations of motion is completely characterized by its infrared series expansion which, in turn, is totally fixed in terms of a finite number of parameters. In the absence of a horizon, these parameters can be the values of the scalar and gauge fields at x o = 0, (φ 0 , A 0 ). A soliton geometry is dual to a field theory state with spontaneously, and not explicitly, broken symmetry. This further demands that the source φ b vanishes, which provides a nonlinear relation between (φ 0 , A 0 ) that determines completely A 0 in terms of φ 0 . A family of soliton solutions is therefore uniparametric. The same reasoning goes through to the black hole case. Now, each black hole solution would be totally determined by the triplet (x h , φ h , A h,1 ) and, again, the condition φ b = 0 would link φ h and A h,1 at the given x h . In this way, a family of black hole solutions is bi-parametric. The backreaction of the gauge field on the geometry is controlled by the coupling e and the probe limit corresponds to e → ∞.

Summary of phases
Following closely the discussion contained in [11] [12], we parameterise the space of solutions at a given coupling, e, by the charge Q and mass M . Three different regimes appear separated by two threshold values, e t and e sr , that signal the appearance of two distinct instabilities. The lower one, e t , marks the threshold for a near-horizon tachyonic instability that affects RN black holes. It is triggered by the fact that the gauge field of an extremal RN black hole makes a negative contribution to the effective scalar field mass, lowering it below the Breitenlohner-Freedman bound of the AdS 2 factor of the near-horizon geometry. It is of the same kind as the one found in the context of holographic superconductors [13] [14]. As discussed in appendix B, it is only relevant for large enough black holes, a fact that explains why it was found in the context of planar geometries. On the other hand, for small black holes, there is a superradiant instability at work beyond e sr [11] [12]. In the following In this table, e t stands for the tachyonic instability threshold, and e sr for the superradiant one. 3 Within each interval of possible values for the coupling we have encountered a similar phenomenology.
In the remaining part of this section we summarize the phase diagram for the static solutions of the Einstein-Maxwell-scalar action for values of the coupling lying on each interval and discuss the soliton solutions in those regimes.
• e < e t Here, the coupling is below the tachyonic instability threshold. The only static solutions present are RN black holes and solitons which exist for values of the charge smaller than Q crit (e), a sort of Chandrasekar limit on regular selfgravitating solutions. In figure  1 (left) the green curve corresponds to solitons and the blue curve to the extremal RN black holes. Regular RN black hole solutions exist in the blue shaded region above this curve.
• e t < e < e sr In this regime, in addition to RN black holes and regular solitons, there exist large and small hairy black holes in a band about extremality for Q > Q 0 (e) (red shaded region of figure 1, right). The red curve denotes extremal hairy black hole solutions and the black line is the instability curve for the RN solution. This function has the limiting behaviours lim e→et Q 0 (e) = ∞ and lim e→esr Q 0 (e) = 0. Solitons, in green on figure 1, appear in two branches. One is connected with the vacuum and has bounded charge Q < Q crit (e). The other one has instead a lower bound Q > Q dis (e). Both bounds merge Q dis (e) → Q crit (e) in the upper limit e → e sr . When Q > Q dis (e) the lowest mass hairy BH is no longer extremal and has non vanishing (possibly diverging) temperature.
• e sr < e Now hairy black holes exist for all values of Q in a band about extremal RN. Concerning solitonic solutions, the two branches seen in figure 1 merge to form a single branch that extends unbounded for all values of Q ≥ 0. To the best accuracy of our numerics, the soliton line is the limit of vanishingly small hairy black holes with diverging temperature. In this limit, the temperature becomes ill defined for the soliton. The situation is very similar to the one in AdS 5 [12]. However, in that case, a second critical charge was found, beyond which hairy black holes appeared below the soliton line. We have not seen any trace of this behaviour in our numerical scan over the space of static solutions. However we cannot discard a fine structure of the form observed in AdS 5 beyond the reach of our numerical accuracy. For that, a perturbative analytical study along the lines used in [12] would be needed.

Soliton branches in more detail
Let us describe the soliton solutions in more detail. They are horizonless, fully backreacted, charged solutions, sourced by a normalizable scalar field profile φ s (x), where equilibrium is attained by an exact compensation between gravitational an electric forces. Due to the normalizability condition, on the field theory side, a soliton represents a macroscopic Bose-Einstein condensate that spontaneously breaks the U (1) symmetry. They come in oneparameter families that can be parametrized by the central value φ 0 ≡ φ s (x = 0). Depending on the value of the coupling e, soliton families display different aspects 4 • For e < e t , there exists a single soliton branch. It is continuously connected to the AdS 4 vacuum, in the sense that, for φ 0 → 0, it reduces to global AdS 4 . Therefore, for φ 0 1, this solution family admits a perturbative construction, and can be described as a nonlinearly dressed ω = 0 AdS scalar normal mode. Besides this fact, the trademark of this branch is the existence of a critical value φ 0 = φ c,1 for which both the soliton mass M (φ 0 ) and charge Q(φ 0 ) attain a maximum. When φ 0 > φ c,1 , M and Q spiral around a limiting value that is reached for φ 0 → ∞.
• For e t < e < e sr , there exist two different soliton branches. The first one, connected with the vacuum, was already present in the e < e t case. The second one, disconnected from the vacuum, is not amenable to a perturbative construction. In this branch, solutions exist for φ 0 larger than some critical value φ c,2a for which the conserved charges (M, Q) diverge. They decrease with φ 0 > φ c,2a until they reach a minimum value at some φ 0 = φ c,2b > φ c,2a . In parallel with the first soliton branch, for φ 0 > φ c,2b , M and Q show damped oscillations around a limiting value that is attained in the φ 0 → ∞ limit. Representative plots of the behavior just described are provided in figure 3 (left).
• For e ≥ e sr , the two soliton branches described in the previous item fuse into a single solution family that is vacuum connected (see figure 3, right). Again, there exists a critical φ c,3 such that M and Q seem to diverge in the limit φ 0 → φ c,3 .

The blow up limit
In figure 4 we plot f s and φ s for representative solutions of the unbounded branch at e = 2.
It is clearly appreciated that, as the conserved charges associated with the solution become larger, the field gradients become more localized near the boundary region. As emphasized in [15], the fact that solitons come in one-parameter families with unbounded charge, allows to take a blow up limit that maps onto a solution with planar geometry. The procedure starts by looking at the near-boundary expansion Introducing new coordinates and redefinitions So far, this is only a reparametrisation of our initial solution. However, if we now take the singular limit λ → ∞ and, simultaneously, move along the soliton branch in such a way that the vector of rescaled quantities χ(λ) ≡ (m(λ),q(λ),μ(λ),φ 3 (λ)) remains finite, we obtain a planar geometry 5 characterized by the hatted quantities. For example, identifying the parameter λ with µ, the dimensionless ratios must tend to constant quantities when µ → ∞ for the planar limiting geometry to have finite energy, charge and vev densities, respectively. In figure 5 the ratios (2.25) are plotted for the soliton branch at e = 5. The planar geometry after the blow up limit is taken clearly matches the extremal hairy black brane geometry studied in [16], at the given e. 6 Incidentally, this observation explains why there does not exist a second soliton branch when e < e t . In that case, in planar AdS, there is no near-horizon tachyonic instability that can trigger hair condensation, so there is no limiting extremal solution which this hypothetical branch could map onto.

Grand Canonical ensemble
The previous section contains a classification of static solutions without examining the issue about their stability. This is not an unambiguous question, as to compare different solutions one must first specify the ensemble. Whereas the plots we have shown are more akin to studying the system in the microcanonical ensemble, physics is more likely related to the grand canonical ensemble, where temperature and chemical potential are feasible knobs. A very generic feature to take into account is that black hole solutions in global AdS (both RN and hairy) split in two kinds: small and large, having negative and positive specific heat respectively. It is the later ones that should be brought into correspondence with stable superfluid quantum states in the dual field theory. In addition to these thermal solutions, there is the one dual to a thermal gas represented by pure Euclidean AdS 4 with a compactified time coordinate. Finally, we have the regular solitons. These static solutions are the building blocks of a rich landscape of first and second order phase transitions. Along this section we are going to report results for e = 3. For any coupling e > e t the situation is qualitatively the same. Below, instead, the transition involves only RN black holes and solitons, since, as evident from figure 1, there are no hairy solutions in this regime of coupling. In this case, our system differs from the one analyzed in [15] [17], where the scalar has a tachyonic mass to start with.

Small and large hairy black holes
The thermodynamics of AdS-RN solutions, both in the canonical and grand canonical ensemble has been examined with great care in the past [18] [19]. Here we are going to study the thermodynamical behavior of hairy solutions. The Hawking temperature is given by and the entropy is where A h is the area of the horizon. In figure 6 (left) we plot the vev of the dual operator O for fixed chemical potential µ = 1.5 and coupling e = 3 as a function of the temperature. The observed behaviour is typical when the gravitational solution corresponds to a small black hole: condensation appears for temperatures greater than some critical value. We will show that this condensed phase is not physically relevant because it has greater free energy than the RN and soliton solutions. This phenomenon, dubbed retrograde condensation in the literature, has been observed in different contexts [20] [21]. In figure 6 (right) we show the temperature as function of the horizon position for small hairy black holes, where their negative specific heat is manifest.  As an example of a small-large black hole transition we rise the value of the chemical potential up to µ = 3. Figure 7 (left) shows the behavior of the vev as a function of the temperature. Points on the upper (lower) segment are small (large) hairy black hole solutions. Upon lowering the temperature the system undergoes a second order phase transition along the lower (stable) branch from a normal to a superconducting state. This can be appreciated in the typical mean field theory behavior of the condensate near the critical temperature, since it goes to zero as O ∝ (T c − T ) 1/2 . The denomination large/small comes, as usual, from the double valuedness of the temperature as a function of the horizon radius x h . This is seen on the right plot of figure 7. Further increasing the chemical potential, we observe the behaviour shown in figure 8 for µ = 15. The same second order phase transition happens here, but the small black hole branch has now disappeared.

Grand potential
Linear stability of the previous solutions amounts to the minimisation of certain thermodynamical potential that depends on the ensemble. In the grand canonical ensemble this is the grand potential Ω(T, µ), which the AdS/CFT correspondence identifies with the on-shell renormalized Euclidean action Ω = T S on−shell (see appendix C for calculations). There is a critical temperature, T c , below which the soliton is the dominant solution. We observe that the entropy is discontinuous at the meeting point between the brown and blue curve showing the appearance of a first order phase transition at T c .
In figure 9 we plot the grand potential (C.7) and the entropy (3.2) respectively as functions of the temperature. In both figures the blue curve refers to the RN solution and the orange curve to the small hairy black holes for µ = 1.5 and e = 3.  The brown dashed line is the free energy of the soliton. 7 We observe that a first order phase transition occurs, whereby increasing T the solitons switch into RN black holes. The small hairy black holes are never the ground state at this value of µ. Figure 10 shows the same functions for the µ = 3, e = 3 solutions. The preferred branch is that of large black holes (lower orange segment) instead of small black ones (upper orange line) which have negative specific heat (see figure 7). Again, the derivative of the entropy is discontinuous at a critical temperature T c 1 and we have a second order phase transition denoting the normalsuperconductor transition of the dual QFT. We detect a second first order phase transition temperature T c 2 , below which solitons are preferred. Finally, figure 11 corresponds to µ = 15. The phase space is the same, as well as the nature of the phase transitions, the only difference being the disappearance of the small black hole solutions.

AdS nonlinear stability
The above sections have relied on a combination of analytical and numerical arguments. The construction of static solutions is performed by solving ODE's and setting up a shooting procedure. It involves fixing, for example, the radius of the desired solution and the value of the scalar field φ o at such radius. Then the value of A t (x o ) is also varied until one obtains a solution with vanishing source φ b = 0 and non vanishing vev φ b,3 = 0. The end result is scrutinized to find the actual value of the mass M and the charge Q of the obtained stationary solution. If more than one solution is available, the free energy is to be invoked in order to select the correct ground state. In a sense, this strategy relies on a certain amount of guesswork. Prior to the construction of hairy black holes in [14], the space of known static vacua consisted of either pure AdS or RN black holes. Later on, solitons where first inferred, and then constructed [11] from a limit whereby the hairy black hole's horizon is shrunk to zero.
In this section we will use a complementary approach. A numerical simulation code for time evolution is the closest thing one can have to a real experiment. In this spirit, the approach starts from the other end: one devises a certain initial radial profile for the bulk fields, with a given total mass M and charge Q, and lets it evolve under a scheme that preserves these values. 8 If the evolution settles down to a certain stationary state, it must necessarily be one in the list above. And if two of them are available with the same values of M and Q, the evolution will select the ground state in the microcanonical ensemble. 9 By evolving an initial condition below the blue curve in figure 4 this would have shown that black holes with abelian hair exist in global AdS, had this paper been written prior to 2008. Moreover, imagine there was any other exotic type of black hole that nobody has constructed yet using static methods (ODE) and, moreover, suppose it has larger entropy. The collapse simulation would smell its existence and the fields decay to that solution after exploring large portions of phase space. We must admit we have not found any new such solution using this, admittedly expensive, method. In figure 12 we have plotted some snapshots of a typical collapse to a hairy black hole. For the values of charge, mass, and coupling e used in this simulation there is no RN black hole available. Nevertheless at t = 14.8 the scalar develops a spike at a point where the metric approaches an apparent horizon, signalled by a zero of the function f (blue curve). Even if the zero value is never reached, the dynamics close to this point becomes extremely slowed down in terms of the time at the boundary x = π/2. At later times, the outer oscillations of the scalar field start piling up on top of the first spike, and the metric function f tries to reach zero at higher values of the coordinate x (see inset). A very high precision and up to 2 17 grid points are needed to push this numerical evolution safely, and resolve the region close to the collapse with enough accuracy, in particular monitoring the constancy of M and Q values throughout the process. The exponentially decaying ringdown ends up in a static solution where the outside hair profile resembles the ones in figure 4 for the soliton solutions. In ref. [4] a strong claim was made, that the geometries exhibiting a corner of instability, were those for which the spectrum of linear perturbations is fully resonant. An important case that shows how sensitive the dynamics is to the resonance property can be found in the case of a scalar field on asymptotically flat space enclosed in a cavity. An instability corner is obtained for Dirichlet boundary conditions [7], where the spectrum is resonant, but not instead in the case of Neumann boundary conditions [6], where it is only asymptotically resonant. 10 The upshot of these analysis seems to indicate that, not only an UV cascade is important, but also a non dispersive spectrum, such that the initial wave packet remains focused at all times.
To test this picture further, one would be interested in families of situations that depart smoothly from it and, in this sense, there are two ways to achieve this. One can keep the original lagrangian setup and, say, consider perturbations around equilibrium solutions other than pure AdS. In generic cases, the spectrum of perturbations of the scalar laplacian on this new background will not be resonant and the collapse will not exhibit a corner of instability. This is the proposed mechanism for the stability of broad initial data, as studied in [6]. They can be understood as perturbations of a single oscillon solutions which are known to be stable. A different strategy is to deform the theory by a continuous parameter. Examples of this are Gauss Bonnet [8], or the inclusion of a hard wall [24]. A remarkable example that stands out is AdS 3 , which despite being resonant, has a mass gap [9] [10]. The present setup offers another interesting example where the theory is deformed by the presence of additional degrees of freedom, while still preserving the resonant spectrum. At the linear level, scalar perturbations obey the same equations of motion as in Einstein-scalar theory. 11 A perturbation of the (now charged) scalar field of amplitude sources a gauge field A µ at the same order, which will backreact on the scalar equation of motion at next order 2 since it couples through the covariant derivatives. Hence we have another work bench to test whether nonlinear perturbations exhibit a corner of instability. In ref. [25] collapse of a complex scalar was considered and no significant difference was obtained for the phenomenology of charged vs. neutral configurations, apart from a small decrease in collapse time. The results in that paper can be recovered in the limit of vanishing coupling e → 0 of our work. However, at finite e we find opposite results as we will now see.

Initial conditions
In this section we will show the results of uniparametric families of collapses that scan across the (Q, M ) plane. The first protocol will involve a set of initial conditions, parameterized with some amplitude . Consider the following family of gaussian initial data that initially fall from the boundary Φ = Φ 1 + iΦ 2 , and Π = Π 1 + iΠ 2 with Φ 2 = Π 1 = 0 and 10 See also [23] 11 We thank Oscar Dias for this observation.
Here the angle β is fixed, and will decrease monotonically towards zero. The cases β = 0, π/2 corresponds to the uncharged initial conditions studied in [1], albeit with an initial pulse that starts infalling from the boundary [26], a fact that is inspired from the physics of a quench. For β = 0, π/2 we are dealing with a shell of charged scalar and gauge field collapsing together.  We want to stress that it is by no means easy to engineer initial conditions that come close to the line of soliton solutions. In particular, within the family of gaussians spelled out in (4.1) and (4.2), by letting β sweep from 0 to π/2, the lines incline up to some point, for some β 0 , where the initial conditions approaches maximally the soliton line and then turn back towards the vertical. These values are, for example, β 0 = 82 • for σ = 0.1 and β 0 = 75 • for σ = 0.2. In principle one can engineer initial conditions that come closer to the soliton line by starting from the other end: namely, by perturbing a soliton, and this will be the subject of the next section. Most remarkable is the fact that is seems impossible to even write initial data whose charge, Q, and mass, M , give a point below the soliton line (white region). This seems to point out that soliton solutions extremize certain positive definite functional that can be derived from the action. It would be very interesting to elaborate on this point further.
Time for collapse is one of the important observables in the game. We can see in plot 13 the case of sharp pulses with σ = 0.1 and β = 0 • , 48 • and 82 • . The horizontal axis represents the initial mass, M , which, for fixed σ, grows with 2 . On the right, the times for collapse for each family are plotted. When moving along an individual series from right to left, the plateaux reflect the number of oscillations that the system undergoes before the final collapse is reached. The fact that electromagnetic repulsion counteracts gravitational attraction is probably behind the fact that oscillations exist at higher values of M when the charge is larger. We see that the behaviour points towards the existence of a corner of instability at the origin even in the charged situation, i.e., no sign of a threshold for stability is appreciated. This seems to confirm the expectation coming from the resonant character of the linearized approximation. From the right figure 13 we draw the important conclusion that charged configurations take longer time to collapse. This is in sharp contrast with the case without gauge field where charged initial conditions were collapsing sooner than neutral ones of the same mass (see figure 6 in [25]).
In the present situation, the charge of the pulse also adds to the defocusing. Still, what figure 13 says is that, for σ = 0.1, this is not enough to erase the instability corner, even for the most charged gaussians that one can device (the right most blue magenta diagonal).
Both plots for collapse time scale with 1/M ∼ 1/ 2 in the limit → 0. However, a closer look at the evolution of the scalar field reveals that, for β = 89 • , the initial gaussian develops subpulses. A reflection of this can be observed in figure 14 which plots the minimum of the metric function min x f (t, x) as a function of time. The evolution still exhibits a quasiperiodic structure where the action of the weak turbulent cascade is apparent in that the minima become sharper and deeper until finally collapse takes over. The 1/ 2 scaling is apparent in the plots of figure 15 for the maxima of the scalar curvature at the origin. The orange curves correspond to the orange dots in figure 13 and the same is true for the brown curves and dots.
After the works in [2,27,28], it has become clear that both focusing and defocusing dynamics (i.e. direct and inverse cascade) seem to be in action and in a delicate equilibrium. For very sharp initial data, small σ, focusing wins. In figure 16, times for collapse with initial width σ = 0.2 are plotted. For vanishing charge β = 0, π/2 ⇒ Q = 0, it scales indefinitely as expected as 1/ 2 as for sharper pulses. However, as soon as some charge is added, we start seeing a divergence at finite values of t collapse . This is exactly the same effect encountered in [5] but here it appears for smaller values of σ than in that case.

Soliton stability
In recent years, configurations analogous to the solitons considered in this paper, known as boson stars, have occupied a prominent role in the study of the AdS nonlinear stability problem. A boson star is a stationary complex scalar field configuration with non zero charge Q that backreacts non trivially on the metric. However, in this case, the bulk U(1) symmetry is global and, therefore, boson stars carry no gauge field A. These solutions provide extended configurations that, once perturbed, help to shed light on how relevant an exactly resonant spectrum of linearized scalar field perturbations is for a weakly turbulent instability to be present on the system. This issue was raised in [4], where the authors showed that just an asymptotically resonant spectrum is not enough to trigger a turbulent cascade along the lines of [1]. Later, it was recognized that boson stars are not endowed with an exactly resonant spectrum and, furthermore, they were shown to be nonlinearly stable [5] [6]. It remains to be seen if the correlation between these two facts survives the present situation. We expect a similar mechanism to be here behind the behaviour of the σ = 0.2 vertical lines in figure  16 which approach the soliton (green) line at constant charge Q. One of the main aims of this section is to provide evidence for this.

Linear stability properties
One of the aims of this subsection is to show that whenever the soliton mass M (φ 0 ) attains an extremum the solutions become linearly unstable. In order to check this explicitly, in we are going to study linearized radial perturbations of the solitonic solutions with an harmonic time dependence of the form cos ωt. Before delving into the details it is useful to notice that, as we are considering a time-reversal invariant and, therefore, nondissipative problem, ω 2 is going to be purely real. In this way, an exponentially growing mode that signals an instability appears whenever ω 2 < 0. We start by fixing our perturbations to be of the form 12 with real φ 1 , φ 2 . Since, due to spherical symmetry, the metric carries no degrees of freedom in our setup, the perturbations defined by (5.1)-(5.4) are not independent. In fact, the reason for having chosen this particular form for the scalar field perturbation is that it allows to solve for δ 1 and f 1 in terms of φ 1 , φ 2 and A 1 , by making use of the momentum and Maxwell constraints (A.12),(A.14) linearized in . Specifically, we get that where C δ , C f are integrating functions that must be fixed by the correct choice of boundary conditions. For harmonic perturbations, set which forces C δ , C f to be zero. Then, we can obtain the equations of motion forφ where M 1 , M 0 a,b are matrix-valued functions that depend exclusively on the background solution. In order to solve (5.9), we have to choose appropriate boundary conditions for Z. At x = 0 we demand regularity. As for the background soliton, this forces Z to be even at x = 0. At x = π/2, boundary conditions arê The first two conditions come from imposing normalizability on the scalar field perturbation (5.4). The last condition demands a more thoughtful explanation. Let us consider the most general near boundary expansion forÂ 1 , In this case, it can be shown that we haveδ 1 (π/2) ∝Â 1,1 and, in consequence, if we want to maintain our gauge choice for the time coordinate, we must setÂ 1,1 = 0. This is tantamount to demanding that the frequency ω is the one measured by a boundary observer. 13 On the other hand, nothing prevents us from allowing thatÂ 1,0 = 0 i.e., perturbations that don't keep fixed the soliton chemical potential. However, examining the explicit form of equation (5.9) we discover that any solution is invariant under the changeφ 1 →φ 1 ,φ 2 →φ 2 + α andÂ 1 →Â 1 + β, provided that αω 2 + eβ = 0. Therefore, we can employ this residual symmetry 14 to setÂ 1,0 = 0 with no loss of generality, fixing it completely along the way. 15 Before discussing how equation (5.9) was solved numerically, let us make a last general comment. The boundary conditions (5.10)-(5.12) and the relations (5.5)-(5.6) imply that the perturbations here considered don't change the charge and the mass of the soliton at linear order. This observation allows for a better understanding of the relation between the soliton linear stability properties and the fact that the mass curve, M (φ 0 ), encounters an extremum at φ 0 = φ 0,c . 16 First, let us mention that, whenever M (φ 0,c ) = 0, we also find that Q (φ 0,c ) = 0. Therefore, around φ 0,c , two infinitesimally close solitons, parameterized respectively by φ 0,c and φ 0,c + ∆φ, have the same mass and charge, up to O(∆φ 2 ) corrections. This implies that there must be a time-independent linear radial perturbation that connects these static configurations and, in consequence, equation (5.9) admits a solution with ω 2 1 (φ 0,c ) = 0, i.e. a zero mode in the soliton spectrum. For ω 2 1 (φ 0 ) at least a C 2 function around φ 0,c , we find which becomes negative one side of the mass curve extremum, signalling an instability. On the remaining part of this section we are going to proceed by solving (5.9) numerically, confirming this expectation. For this task we employed Tchebychev pseudospectral collocation method. Inserting expansionsφ into equation (5.9) and evaluating on a collocation grid {t k , k = 1...N } an algebraic generalized eigenvalue problem is to be solved which gives the numerical values of the first soliton normal modes. 17 As for the results, first, the spectrum thus found is not resonant. In the spirit of [4], this should entail the absence of a turbulent cascade in the fully nonlinear regime. In figure 17, we plot the first eight normal frequencies for the soliton branch at coupling e = 5. They consistently reduce to their global AdS values as φ 0 → 0. A remarkable feature is the mode splitting that occurs for k ≥ 3. 18 This is common to every perturbative soliton branch we have analyzed. If the spectrum was exactly resonant, we would have that, for k ≥ 2, We plot the right hand side of (5.19) for each k in figure 17 (blue-dashed). It is clearly seen that the equality is not satisfied away from φ 0 = 0. A similar exercise can be performed only between the lower or upper splitted eigenfrequencies, choosing as reference the difference between any two consecutive ones, with identical conclusion.
Let us consider the intermediate region 3/2 ≤ e 2 ≤ 9/2. In figure 18 (left), we plot ω 2 1 (φ 0 ), ω 2 2 (φ 0 ) for the vacuum connected soliton branch at e = 2, together with the rescaled Q(φ 0 ) curve. We clearly see that a zero mode develops precisely at the point where Q(φ 0 ) reaches its first maximum and that, past this point, the solutions become linearly unstable. 19 The same phenomenon can be clearly appreciated on the vacuum disconnected soliton 17 We discard the values that don't converge when N is increased. A convergence test for the pseudospectral code is provided in appendix D. 18 We remind the reader that a mode splitting was previously found in the perturbative computation of [6]. 19 With a resolution of δφ 0 = 10 −3 , we have determined that the maximum lies at φ 0,c = 0.904. Our pseudospectral code produces the values ω  branch after Q(φ 0 ) attains its first minimum ( figure 18, right). The general pattern we find is that, when Q(φ 0 ) hits a new extremum, a new normal mode crosses zero. 20 As an additional comment, note that, regarding the vacuum disconnected branch, and in the Q 1 regime, ω 2 1 (Q) is a decreasing function of Q that stays finite in the Q → ∞ limit ( figure 19, left). Instead, in order for the phase of the harmonic perturbation ωt to remain finite in the blow up limit (2.19), the frequency should scale up as ω ∼ µ ∼ Q 1 2 . We conclude that harmonic linear perturbations die off when the blow up limit is taken. This is consistent with the fact that the soliton branch maps onto a T = 0 hairy black brane, for which linearized perturbations correspond to quasinormal rather than normal modes. The discussion goes through in parallel to the regime e > e sr (see figure 19 (right) for ω 2 1 (Q) at e = 5).
After having discussed the linear stability properties of the solitons, in the next section we move on to the study of their nonlinear stability.

Nonlinear stability properties
We consider now the effect that a localized scalar field perturbation has on the soliton. We will stick to the same family of initial conditions that were used to perturb the AdS 4 vacuum (4.1)(4.2) but in a purely real setup. Concretely, our initial condition will be φ s + φ with φ(0, x) = 0 and As this configuration has zero charge (see (A.21)), the family of perturbed solitons that we use to start with spans a vertical line in the (Q, M ) plane above the unperturbed soliton solution φ s , like the magenta vertical sets in figure 16 (left). The difference now is that this set of initial conditions explores down to the bottom green line as → 0. Despite the fact that placing a perturbation like (5.20) on top of the AdS vacuum or on top of a soliton leads to very different initial conditions, the phenomenology we discover is remarkably similar. Namely, the magenta lines found in figure 16 (right) are qualitatively reproduced here. Indeed, for high enough prompt collapse is observed. Below some threshold mass M c we have a delayed collapse and a number of oscillations are completed before the system finally undergoes gravitational collapse. This number, and with it the final time for collapse, diverges rapidly at some value of the mass above the soliton curve. Hence we don't see any trace of a nonlinear instability corner centered at the soliton solution (instead of the AdS vacuum). This is presumably again a symptom of the nonresonant character of the spectrum of soliton perturbations. The possible survival of this oscillating region in the blow up limit is a relevant question and, indeed, was one of the main motivations that started the present work; the answer is negative. More precisely, we are interested in establishing whether its width M c − M s has a finite size relative to M s M c − M s ∼ M s ∼ Q 2/3 s . In this case, for M s 1, entering the oscillating regime requires to fine-tune the initial perturbation (5.20) in such a way that its relative contribution to the system energy goes to zero. The behaviour just described is fairly natural given that, in the blow up limit, the soliton branch we are perturbing reduces to an extremal hairy black brane with an AdS 4 near-horizon geometry [16]. It is not unreasonable that, above this background, any perturbation localized in the near-boundary region, regardless of its amplitude, leads to direct collapse to a T = 0 black brane, in the same vein as it happens for zero charge [31]. It remains to be seen if this is also the pattern in other matter models. Our educated guess is that, for any theory that displays an unbounded soliton branch, the blown up extension of the oscillatory regime is finite whenever the theory can support a gapped spectrum of scalar fluctuations in the planar limit. Theories of this kind are not unknown. Consider, for instance, the Improved Holographic QCD models which have recently been analyzed at the dynamical level in [32] [33]. There, the planar geometry dual to the field theory ground state is sourced by a nontrivial scalar field profile that generates a naked singularity in the infrared. 21 When considered in the fully nonlinear regime, weak perturbations localized near the boundary may be noncollapsing and forever oscillating, since the singularity would "repel" them from the infrared, so as they never reach their Schwarzschild radius. This behavior has been explicitly seen in more crude models of gapped field theories in planar AdS, such as a scalar field in a hard wall geometry [24] or the AdS-soliton [34]. It would be interesting to classify, in generic terms, which Einstein-Maxwell-scalar theories support a soliton branch with a gapped planar limit and check if this is correlated with a nonvanishing scaling of the oscillatory regime width.
We have also determined the boundary of the prompt collapse region for the vacuum connected soliton branch at coupling e = 2, upon scalar fluctuations of the form (5.20) with σ = 0.05 (figure 20, left). In accordance with the fact that soliton solutions become linearly unstable when the mass curve reaches its maximum, here we find that the width of the oscillation region shrinks to zero in a linear fashion.

Conclusions
The original motivation of this paper was to examine the status of the nonlinear instability problem in the Einstein-Maxwell-scalar theory in d + 1 = 4 dimensions. The work demanded first a thorough unraveling of the landscape of static solutions which showed almost the same features as in one dimension higher [11] [12]. The issue of stability can be examined at different levels. Studying the thermodynamical linear stability demands selecting first a certain ensemble. Working in the grand canonical ensemble we have pinned down the relevant phase transitions. In the second part of the paper we have studied the nonlinear stability in the microcanonical ensemble. First of all, we have considered perturbations of pure AdS 4 with families of initial conditions of different charges, generalizing the uncharged ones considered in [1]. These initial conditions have an amplitude and a width σ. We have found that for thin initial pulses σ = 0.1, the same type of instability sets in at times of order −2 . This fact confirms the expectation that places the origin of the mechanism in the fully resonant character of the linear spectrum, something that is also true in this case. Wider pulses, σ ≥ 0.2, exhibit a divergence in the time for collapse below some critical amplitude. This is also in parallel with the effect detected for wide radial perturbations of AdS 4 in [5]. The role of the oscillons in that situation is probably taken up by the solitons in this setup. To confirm this picture, we have examined the linear and nonlinear stability of fluctuations placed on top of the soliton solutions. The spectrum of linear normal modes is not resonant. The linear stability, signalled by an imaginary eigenfrequency, is seen to appear, as expected, coinciding with extrema of the mass curve. For higher amplitude perturbations we need to resort to a full fledged simulation of the evolution of the system. In general, linear stability extends to nonlinear stability up to some threshold for the amplitude of the perturbation. Beyond that, the soliton is destroyed and collapses to a hairy black hole. The protection region where oscillations do not decay does not scale properly in the limit of large mass and charge to survive the blow up planar limit. This seems to point to the necessity of having a mass gap, and not just a mass scale, to find such oscillatory behaviours. While this paper was completed in its last part, namely sections 4 and 5, reference [17] appeared in the arXive. It overlaps substantially with the content of our section 3, and we find agreement with their results whenever they can be compared (their scalar field has negative mass).

A Ansatz and equations of motion
Let us give the equations for a situation which is slightly more general than the one considered in the main text, and allow for a scalar potential Sitter. Also, the scalar field is complex and we have D µ φ ≡ (∂ µ − ieA µ )φ. The equations of motion are where the energy-momentum and charge currents are By employing the ansatz described in section 2, the Klein-Gordon equation can be casted into the first order forṁ while, from the Einstein equations, we obtain, There is one additional equation coming from the (t, x) component of Einstein equations that yields the momentum constrainṫ Concerning Maxwell's equations, define C = A t (t, x), then from (A.4) we can derive the two following equations which can be easily shown to be compatible. The first one can be recasted as follows and the second is the Maxwell constraint. With the condition that A t be bounded at the origin we can integrate (A.15) to find The electromagnetic current is given by the following expression Take (A. 16) and notice that the following expression holds 19) or, after dividing by (2 − d) and multiplying by e δ(π/2) , because Q has to be invariant under δ → δ + c, In other words The Reissner-Nordström black hole solution can be recovered by setting φ = Φ = Π = δ = 0. Specializing to the d = 3 case, the stationary solution to (??) and (A.15) is This expression gives the mass and charge of the black hole as The charge is always proportional to the derivative of the gauge field on the boundary, i.e. the charge density, which agrees with (A.20). More generally, the mass of a black hole in the time dependent background (2.1) can be written

B Instabilities
Take the RN metric We will show that extremal RN black holes suffer from two potential instabilities.

B.0.1 Tachyonic instability
This is sourced by the lowest scalar mode close to the horizon. The near horizon limit is obtained by expanding and taking the limit λ → 0. In the extremal case f (r h ) = f (r h ) = 0 and, therefore, the limit has the form of AdS 2 × S d−1 and also One may take the near horizon limit in the equation of motion of the scalar field and extract the mass from the coefficient of φ We observe explicitly that the mass gap is set by the chemical potential. The scalar field will be unstable if its mass violates the BF bound of AdS 2 , hence if m 2 < −1/4. This implies e > e t = 1 4 The values of M and Q at extremality are given by Evaluating f (r h ) with these values results in the following expression , which diverges in the limit r h → 0. The tachyonic instability is thus suppressed for infinitesimally small RN back holes. Instead, for large r h , e t (r h ) admits the following expansion In summary, for e 2 t < 3/2 in d = 3 and for e 2 t < 3 in d = 4, there is no tachyonic instability. As this threshold is surpassed, large black holes are unstable first. Small black holes do not suffer from this instability.

B.0.2 Superradiant instability
If eµ > ∆ 0 , with ∆ 0 being the lowest eigenvalue of the scalar linearized equation, RN black holes becomes superradiant in the limit r h → 0, meaning that this mode scatters off with a reflection coefficient |R| > 1. The final state should be a hairy black hole with greater entropy than the RN one, which dominates the micro-canonical ensemble. For a massless scalar, ∆ 0 = d. Substituting the extremal value gives, for r h 1, which constitutes the relevant lower bound on the coupling, since µ < µ ext . The precise superradiance threshold is Hence all extremal RN black holes above e 2 = 9/2 are unstable.
• in d = 4, This is the result of [12]. For e 2 > 32/2 all extremal RN black holes are superradiantly unstable.
It should be noted that the condition eµ > ∆ 0 is obtained in a perturbative expansion in r h and, in consequence, it can't be extrapolated to sufficiently large RN black holes.

C Renormalized action
The central quantity in the grand canonical ensemble is the grand potential Ω(T, µ). The AdS/CFT correspondence identifies it with the on-shell regularized Euclidean action. To compute it we continue to Euclidean signature and compactify time with a period 1 T . Writing Using the equations of motion one can show that 2L = −(G t t + G x x ) whencê S bulk = 4π x b x h dx sec 2 (x)e −δ(x) − e −δ(x) tan(x) sec 2 (x)f (x) | x=x b (C. 3) where the (diverging) result has been regularized at some boundary value x b = π/2 − . To have a well defined variational problem we must add the Gibbons-Hawking term whereḡ stands for the induced metric at x b , K = g µν ∇ µν n ν is the extrinsic curvature and n µ = cos(x) f (x)δ µ x is the outward pointing unit normal vector to the boundarŷ S GH = −2π tan(x)e −δ(x) tan(x)f (x) + f (x) −2 tan(x)δ (x) + 6 sec 2 (x) − 2 x=x b . (C.5) The sumŜ bulk +Ŝ GH diverges when x b → π/2. This can be handled by the following countertermŜ The grand canonical thermodynamic potential Ω is obtained from the limit → 0 We checked the Smarr relations for this expression which, in particular, imply that, indeed, Ω = M − T S − µQ, as demanded by thermodynamic consistency. If we wish to compute the thermodynamic properties in the canonical ensemble (fixed Q) instead the relevant quantity to consider would be the free energy F = M − T S.

D.1 Time-evolution code
It is a virtue of the coordinate system (2.1) that the Einstein and Maxwell equations appear as constraints (A.10), (A.11) and (A.16) that can be solved at each instant of time. The evolution of the system is then driven by the scalar field equations (A.8) and (A.9). Starting from given nonequilibrium initial data, such as (4.1) and (4.2), we have solved these equations numerically by resorting to a fourth-order accurate finite-difference evolution code.
Time evolution is performed by an explicit Runge-Kutta method. In order to deal with the high-frequency noise generated due to the finiteness of the discretization grid, we implement standard Kreiss-Oliger dissipation. By setting δ(π/2) = 0 we obtain stable evolutions with a constant Courant factor λ with no need of local mesh refinement in time. On the other hand, spatial derivatives are discretized by employing a centered finite-difference stencil, while integrations are handled by a specifically designed routine, based on local polynomial interpolation. To deal with boundary conditions and numerical stability both at x = 0, π/2 requires some detailed procedures that can be found in [35].
The major difficulty in the present setup stems from the fact that, upon evolution, the scalar profile develops very spiky features that demand a high resolution. To resolve these sharp features, which are apparent in figure 12, we used global mesh refinement in space, eventually reaching 2 17 + 1 grid points to discretize the interval x ∈ [0, π/2]. This has required a parallel implementation that employs the MPI infrastructure to run the code on the SVG cluster at the CESGA facility (www.cesga.es). Optimal results have been obtained for ∼ 30 nodes running in parallel. Smarter solutions involving local space mesh refinement are left for the future.
The quality control parameters employed to activate the refinement process are both the norm of the momentum contraint (A.12), as well as the relative mass loss at each time step. As a matter of fact, only at late times in the simulation are such mentioned fine resolutions required. The code stops at a time t f when the minimum value of A(t, x) reaches below an user defined cutoff A c . 22 This is the time that is meant in the right figures 13 and 16. Of course, mathematically speaking, the apparent horizon will only form in the infinite future lim t f →∞ min[A(t f , x)] = 0 in the chosen coordinate gauge.

D.2.2 Soliton eigenfrequencies pseudospectral code
For the computation of the normal modes on top of a soliton background, we have resorted to the pseudospectral method described in the main text. As mentioned, the output of this procedure are the first N soliton normal modes. For pseudospectral methods, we expect exponential convergence, since we are approximating the analytical scalar field eigenmodes. In order to determine the convergence properties of the method, we show in figure 22 the quantity ∆ω 2 k (N ), defined as ∆ω 2 k (N ) = ω 2 k (N + 1) − ω 2 k (N ) . (D.5) The fact that ∆ω 2 k (N ) → 0 exponentially as N → ∞ implies that the sequence {ω 2 k (N ), N = N 0 , N 0 + 1, ...} converges as anticipated.