Hairy black holes and the endpoint of AdS$_4$ charged superradiance

We construct hairy black hole solutions that merge with the anti-de Sitter (AdS$_4$) Reissner-Nordstr\"om black hole at the onset of superradiance. These hairy black holes have, for a given mass and charge, higher entropy than the corresponding AdS$_4$-Reissner-Nordstr\"om black hole. Therefore, they are natural candidates for the endpoint of the charged superradiant instability. On the other hand, hairy black holes never dominate the canonical and grand-canonical ensembles. The zero-horizon radius of the hairy black holes is a soliton (i.e. a boson star under a gauge transformation). We construct our solutions perturbatively, for small mass and charge, so that the properties of hairy black holes can be used to testify and compare with the endpoint of initial value simulations. We further discuss the near-horizon scalar condensation instability which is also present in global AdS$_4$-Reissner-Nordstr\"om black holes. We highlight the different nature of the near-horizon and superradiant instabilities and that hairy black holes ultimately exist because of the non-linear instability of AdS.

One of the most challenging open questions of general relativity concerns finding the time evolution and endpoint of the superradiant instability in the Kerr-AdS black hole (or in systems that mimic this system). When a wave with frequency ω and azimuthal angular momentum m φ scatters a Kerr-AdS black hole with angular velocity Ω H , its amplitude can grow if ω < m φ Ω H [1][2][3] (see [4,5] for recent reviews on superradiance). This follows from a linear perturbation analysis [6][7][8]. In global AdS, the amplified wave reflects at the AdS boundary and keeps being amplified driving the system unstable. Finding the properties of the instability evolution and its endpoint requires a non-trivial (3 + 1)-dimensional numerical simulation not yet available.
However, valuable insights on the possible evolution can be obtained without solving the actual initial value problem. Indeed, it was found that the zero mode (i.e. onset) of the instability signals the existence of a new branch of black hole solutions with scalar hair [9][10][11][12] or even gravitational hair [13]. Per se the existence of these hairy black holes is interesting for two main reasons: 1) they immediately rule out any attempts of proving no-hair theorems (under non-restrictive assumptions), and 2) they have a single Killing vector field that is also the horizon generator, thus showing that the assumptions of Hawking's rigidity theorem can be evaded [11,13,14]. Of relevance for the initial value problem, it turns out that, given a mass and angular momentum, these hairy black holes always have higher entropy than the Kerr-AdS black hole. This property and the fact that they are connected to the onset of superradiance suggest that these black holes could be the endpoint of rotational superradiance. However, this is not the case because all of those constructed to date have Ω H L > 1 (where L is the AdS radius) and are thus still superradiantly unstable [11,13,15], this time to superradiant modes with higher m φ . This, together with other observations that can be found in the original articles, led the authors of [11,13,16] to conjecture that the superradiant instability might evolve into one of two possibilities: 1) the endpoint of the superradiant instability is a singular solution reached in finite time which implies cosmic censorship violation, or 2) the time evolution develops higher and higher m φ structure (with increasing entropy) and the system never settles down. In the latter case quantum gravity effects would become relevant at some point and cosmic censorship would also be violated, at least in spirit. In either case, we would have an explicit example of cosmic censorship violation in four dimensions and in a system that can be formed by gravitational collapse. This would show that the violation found in [17,18] can also occur in a four-dimensional geometry.
To downsize the level of difficulty of the problem while still keeping some essential ingredients and associated questions, we can turn-off the rotation of the black hole and consider a static, but electrically charged, black hole. Indeed the global AdS Reissner-Nordström black hole still has a generalized (effective) ergoregion where negative energy states can be excited [15,[19][20][21][22] and the initial value problem is now a simpler (1 + 1)-dimensional problem. Such a black hole with chemical potential µ can then be unstable to superradiance when a scalar field with frequency ω and charge q scatters it if ω < µq [15,[20][21][22]. Similar to the rotating geometry case, the onset of charged superradiance is associated to novel hairy black holes with a charged scalar field floating above the horizon: electric repulsion balances the scalar condensate against gravitational collapse. In global AdS 5 , these charged hairy black holes were constructed perturbatively and then numerically in [23][24][25][26][27].
For a given energy and charge, the hairy solutions of [23][24][25]27] have higher entropy than the AdS 5 Reissner-Nordström black hole. Therefore, unless the phase diagram of the theory contains a (unlikely) solution with even higher entropy, these hairy black holes should be the endpoint of charged superradiance. Unlike the rotating case, this is a robust statement since q is fixed once we choose the theory and thus the hairy black holes are not unstable to superradiance (in the rotating case m φ is a quantized integer that is not fixed by the theory so hairy black holes associated with a given m φ are always unstable to higher m φ 's). This is where the charged superradiant system becomes a less interesting toy-model to discuss the endpoint of rotating superradiance. 1 Nevertheless, moving beyond the initial motivation, there are good reasons to investigate the charged system. Firstly, the properties of the evolution of the unstable charged system should be known. This was addressed by a recent numerical time evolution simulation which confirmed that, in AdS 4 , an unstable Reissner-Nordström black hole indeed evolves towards a charged hairy black hole [28].
A second reason is related to the non-linear, weakly turbulent, instability of global AdS. If we perturb global AdS with a scalar field or with a wave packet of gravitons, the system evolves towards the formation of a black hole, for arbitrarily small initial data, as several studies suggest [29][30][31][32][33][34][35][36]. If the spherically symmetric initial data is made of a neutral scalar field the final black hole is an AdS-Schwarzschild black hole. However, if the scalar field is charged then the endpoint of the nonlinear instability should be a hairy black hole since it has higher entropy than the AdS-Reissner-Nordström with same charge and energy. This was confirmed to be the case by numerical simulations in global AdS 4 and AdS 5 [31,37,38].
Since non-linear simulations are computationally costly, the numerical evolution studies of [28,31,37,38] were necessarily limited in the range of parameters and initial data. Thus, it would be desirable to construct the hairy black hole solutions of Einstein-Maxwell theory in global AdS 4 directly as solutions of the elliptic boundary value problem. This would confirm (or not) that the qualitative features of hairy black holes of AdS 5 extend to AdS 4 . It would also give quantitative thermodynamic quantities that can be used to compare with the endpoint of future numerical simulations similar to [28,31,37,38]. This is a main aim of our manuscript.
A third reason to be interested on the charged AdS system, even in the absence of rotation, is related to a third instability of AdS known as the near-horizon scalar condensation instability [39]. The superradiant instability is present only in global AdS black holes and requires (for charged superradiance) that both the scalar field and the geometry are charged. On the other hand, the near-horizon scalar condensation instability (and associated hairy black holes) was first found in planar AdS and, in the AdS/CFT duality context, triggered the research field of "holographic superconductors" and AdS/condensed matter correspondence [39][40][41]. This instability is present in any AdS geometry that has an extremal (zero temperature) configuration even if the spacetime and the scalar field are neutral [41,42]. It occurs whenever the scalar field satisfies the asymptotic AdS d Breitenlohner-Freedman (BF) bound [43] of the full theory but violates the AdS 2 BF bound (associated with the near-horizon geometry of the extremal black hole of the theory). In particular, the near-horizon instability is present in the global AdS-Reissner-Nordström black hole, where it co-exists with the superradiant instability [25,42]. These two instabilities are entangled but in certain limits their distinct nature emerges: large (small) global AdS-Reissner-Nordström black holes, in AdS radius units, only have the near-horizon (superradiant) instability. In Section III, we will discuss the sharp distinction of these two instabilities in AdS 4 .
Since the superradiant instability is present in small AdS-Reissner-Nordström black holes, the associated hairy black holes can be constructed perturbatively around global AdS as a small expansion in the horizon radius and in the scalar condensate amplitude [23][24][25]. (On the other hand, large hairy black holes that emerge from the nearhorizon instability do not have such a perturbative expansion). A few observations immediately suggest the existence of a perturbative expansion in the superradiant case. As mentioned above, a complex scalar field with frequency ω and charge q scattering a charged black hole with chemical potential µ is superradiantly amplified when [20, 21] 2 Re ω < q µ. (1) For definiteness, consider a Reissner-Nordström black hole in global AdS 4 (RN-AdS 4 ), with cosmological radius L, horizon radius r + and chemical potential µ. Its gauge potential is thus A = µ (1 − r + /r) dt. Alternatively, we can parametrize this solution with r + and the temperature T in which case µ = √ 2 1 + 3r 2 + L 2 − 4πr + T . These black holes are regular all the way to extremality (where their temperature vanishes), and the chemical potential is where in the last approximation we considered the small horizon radius limit of the condition. Consider now a scalar perturbation in this RN-AdS 4 background that we Fourier decompose as φ(t, r) = e −iωt φ(r). If we impose reflecting asymptotic boundary conditions that preserve the energy and charge of the system, RN-AdS 4 behaves as a confined box whereby waves can be dissipated only through the horizon. Therefore, only certain frequencies can fit inside the AdS 4 box: the frequency spectrum is quantized. The inverse of the imaginary part of the frequency gives the dissipation timescale (if Im ω < 0) or an instability growth rate (if Im ω > 0). In the limit where r + → 0, one has Im ω → 0 and the frequency spectrum reduces to the normal modes of global AdS 4 , where p is the radial overtone of the perturbation that gives the number of radial nodes of the associated wavefunction (in this manuscript we will be interested on the scalar solutions with lowest energy which correspond to p = 0). Altogether, conditions (1)−(3) imply that incident scalar waves on a small near-extremal RN-AdS 4 black hole drive it unstable to superradiance if where we introduced the dimensionless scalar field charge e ≡ qL. Thus, for p = 0 small near-extremal black holes are unstable for e 3 √ 2 . By continuity, they should also be unstable when we move away from extremality up to a point where the instability switches-off (this point is the onset of superradiance if we do the reversed path, from large to low temperatures).
In the absence of a horizon, a normal mode of global AdS 4 , with frequency ω and amplitude ε, solves the Klein-Gordon equation. Going beyond the linear order, this scalar field sources the Einstein equation at order O ε 2 and higher and back-reacts in the background geometry. If we impose reflecting boundary conditions on the boundary of global AdS 4 , it is conceivable that the solution can be back-reacted (while being smooth everywhere) to any order in a perturbative expansion in the scalar amplitude ε. In the back-reaction process the normal mode frequency should itself receive corrections at each order. This possibility turns out to be correct and the back-reaction of a normal mode of global AdS is known as a boson star. Indeed they were explicitly constructed in AdS 5 in [23][24][25][26] and studied with much detail in [26]. Similarly, these solutions exist in AdS 4 and will be constructed in Section IV. We have the freedom to do a U (1) gauge transformation that eliminates the phase of the scalar field at the expense of changing the gauge potential. So the 'boson star' can be equivalently written as a 'soliton', i.e. a time-independent horizonless solution that is regular everywhere. In most of our manuscript we will adopt this perspective.
When a theory has a solitonic solution and a black hole solution that is unstable, we should consider, as a rule of thumb, the possibility that there might exist a third solution that describes a hairy black hole constructed by placing the small black hole at the core of the soliton 3 . This was first observed in the Einstein−Yang-Mills theory: a black hole can be added to the Bartnik-Mckinnon soliton [48] of the theory leading to a coloured black hole with non-abelian Yang-Mills hair [49][50][51]. This rule of thumb also materializes in the gravitational Abelian Higgs model (i.e. Einstein-Maxwell theory with a complex scalar field) that we are considering [23][24][25]. The theory has a small RN-AdS 4 black hole that we can put on top of the soliton that we construct in Section IV.
Quite remarkably, the idea encoded in the above rule of thumb can be realized in a simple 'non-interacting thermodynamic model' that captures the leading order thermodynamics of the system without solving the equations of motion of the theory 4 . The model simply assumes that at leading order the energy and charge of the hairy black hole is just the sum of the energies and charges of its RN-AdS 4 and solitonic components. The distribution of charges among the two constituents is determined by the condition that the entropy must be maximized. This is equivalent to require that the chemical potential of the RN-AdS 4 and solitonic component match, i.e. that the system is in thermodynamic equilibrium. We apply this model in Section V and will find that it gets the correct leading order thermodynamics of the system.
Moving beyond the leading order inspection of the system, the hairy black holes will be constructed perturbatively in a double expansion in the dimensionless horizon radius r + /L and scalar condensate amplitude ε in Section VI. The equations of motion of the gravitational Abelian Higgs theory are solved using a standard matching asymptotic expansion procedure (see e.g. [11,12,[23][24][25]) whereby we consider a near region r+ L ≤ r L 1 and a far region r r + which overlap in the region r+ L r L 1 if we consider small black holes, r+ L 1. The resulting thermodynamic quantities of the hairy black hole − energy, charge, chemical potential, temperature, entropy, and the Helmoltz and Gibbs free energies − are then computed in (87) and obey the first law of thermodynamics. Confirming the rationale beyond their construction, hairy black holes merge with the RN-AdS 4 family at the onset of superradiance and their zero-horizon radius limit is a soliton. The reader not interested on the construction details can immediately jump to these expressions that describe uniquely the hairy black holes and to their physical discussion in Section VII.
The hairy black holes and solitons of the theory add new competitions in the phase diagram of thermal phases of the gravitational Abelian Higgs model. The discussion of these competitions depends on the thermal ensemble − microcanonical, canonical or grand-canonical − that we consider. It turns out that, as summarized in the phase diagrams of Fig. 1, in the microcanonical ensemble hairy black holes are always the dominant thermal phase. However, large RN-AdS 4 are still the preferred thermal phase both in the canonical ensemble (see phase diagram of Fig. 3) and grand-canonical ensemble (see phase diagram of Fig. 6).
A final comment on the near-horizon scalar condensation instability and its associated hairy black holes is in order. In Section III B we will find the conditions under which a scalar field in an extremal RN-AdS 4 is unstable to the near-horizon instability. This happens when the dimensionless scalar charge e is above the bound (39). In the limit where the dimensionless horizon radius R + → ∞, we recover the known result for planar RN-AdS 4 black holes, (for a massless field) [39][40][41]. Global RN-AdS 4 black holes are unstable for larger values of the scalar field charge. However, in the opposite limit, R + → 0, we get e 2 ≥ 1 8R 2 + + O(1), which indicates that the nearhorizon scalar condensation instability is suppressed for small (global) RN-AdS 4 black holes. This is to be contrasted with the superradiant instability and highlights the different nature of the two instabilities. In particular, this also means that hairy black holes branching-off from the RN-AdS 4 family do not admit a small perturbative expansion in the dimensionless horizon radius R + (and ε). A numerical non-linear construction, that we leave for future work, would be required.
Summarizing our conclusions, the phase diagram of AdS 4 -Einstein-Maxwell theory with a complex scalar field depends on the window of the dimensionless scalar field charge e that we consider. There are three cases: • e 2 < 3 2 : in this case RN-AdS 4 is stable both against superradiance and the near-horizon scalar condensation instability. Thus, there are no hairy black hole solutions. We do have a 1-parameter soliton family. Based on the AdS 5 results of [23][24][25] and specially of [26], very likely this soliton has an intermediate Chandrasekhar limit at a critical charge and mass and then an intricate discontinuous branch structure. Finding wether this is the case for large charges also in AdS 4 requires a numerical study beyond the perturbative analysis done here.
• 3 2 ≤ e 2 ≤ 9 2 : the near horizon instability takes place near extremality for large RN-AdS 4 black holes, i.e. only above a critical mass and charge. Consequently, the associated large hairy black hole solutions cannot be constructed perturbatively and would require a full non-linear numerical construction. Since they do not have a perturbative construction their zero-horizon radius limit is likely singular and certainly not the soliton (this is the case in AdS 5 [23][24][25]). The soliton, on the other hand has similar properties to those of the previous case.
• e 2 ≥ 9 2 : Small RN-AdS 4 black holes are unstable to the superradiant instability and the associated hairy black holes can be constructed perturbatively for small mass and charge (see Section VI). They are a 2-parameter family of solutions that in a phase diagram span an area bounded by the onset of superradiance (where they merge with RN-AdS 4 ) and, at their zero-horizon radius limit, the soliton (see Fig. 1). Eventually, for larger values of the charges (that are not captured by our perturbative analysis) the zero-horizon radius limit of the hairy black hole might no longer be the soliton (this is the case in AdS 5 [25]) but answering this requires a numerical analysis. These hairy black holes always have higher entropy than the RN-AdS 4 black hole in the window of energy and charge where the two phases co-exist. Interestingly, in the phase diagram of solutions of the theory, the hairy black holes connect the onset of the (linear) superradiant instability with the key player − the soliton − of the non-linear weakly turbulent instability of AdS 4 . Indeed, at the simplest but fundamental level, this non-linear instability follows from the appearance of irremovable secular resonances when two or more solitons are taken as initial data and collide [29,30,53].
A perturbative construction of the hairy solutions for the gravitational Abelian Higgs model in AdS 5 was done in [23][24][25] and agree extremely well (for small charges) with the full nonlinear solutions constructed in [25,26]. Our perturbative construction finds hairy solutions for the AdS 4 theory. Qualitatively, we do not find differences between the AdS 5 and AdS 4 phase diagrams. Therefore, [25,26] and our results suggest that (at least for small black holes) the qualitative phase diagram of gravitational Abelian Higgs theory is independent of the dimension. In addition, our quantitative results can be used to compare against and testify the hairy solutions that are the endpoint of initial value simulations like the ones reported in [28,31,37,38]. Hairy black hole solutions of gravitational Abelian Higgs model in AdS 4 for a scalar field mass of m 2 = −4/L 2 were recently constructed numerically in [54] and discussed in the grandcanonical ensemble. Their qualitative results agree with our conclusions for the grand-canonical ensemble (namely, hairy black holes are always subdominant with respect to RN-AdS 4 ) indicating that the qualitative conclusions, at least in the grand-canonical ensemble, might be similar for different scalar masses. The reader interested on superradiance and asymptotically AdS 5 ×S 5 hairy black holes that are solutions of supergravity can find an exhaustive discussion in [24,27].
The plan of the manuscript is as follows. In Section II we introduce the theory, its equations of motion and the boundary conditions for the elliptic problem, and describe how the thermodynamic quantities can be computed. Section III briefly reviews the RN-AdS 4 solution and discusses in detail the conditions where they are unstable to nearhorizon scalar condensation and/or superradiant instabilities. The soliton (boson star) of the theory is constructed perturbatively in Section IV. The leading order thermodynamics of hairy black holes is discussed in Section V using the non-interacting thermodynamic model. In Section VI, hairy black hole solutions are constructed perturbatively using a matched asymptotic expansion analysis. The thermodynamical and physical properties of the hairy solutions are analysed in detail in the three thermal ensembles in Section VII. For the benefit of a reader who wants to reproduce our results, the technical appendices A, B, and C give the expressions for the field coefficients that appear in the perturbative expansions of the superradiant instability of Section III D, of the soliton construction of Section IV and of the hairy black hole of Section VI.

A. Field ansatz, equations of motion and boundary conditions
We consider 4-dimensional Einstein-Maxwell gravity with negative cosmological constant Λ = − 3 L 2 , minimally coupled to a charged complex scalar field φ with charge q. The associated action is where R is the Ricci scalar, F = dA with A being a gauge potential and D µ φ = ∇ µ φ − iqA µ φ is the associated gauge covariant derivative. We consider a potential V (|φ|) = m 2 φφ * , with m being the mass of the scalar field. For definiteness, onwards we set m = 0 but our nonlinear construction can be extended to a massive scalar field or to a more generic potential. We fix Newton's constant as G 4 = 1.
We are interested on solitonic and black hole solutions of (5) that are static, spherically symmetric and asymptotically global AdS 4 . We can use reparametrization of the time t and radial coordinates r, r →r(r) and t →t = t+H(t, r) to fix the gauge to be such that the radius of a round S 2 is r and there is no cross term dtdr (this is often called the radial or Schwarzschild gauge) 5 . A field ansatz that accommodates the desired symmetries is where dΩ 2 (2) describes the line element of a unit radius S 2 (parametrized by {θ, ψ}, say). Once a static solution is found we can always use a U (1) gauge transformation, ϕ → ϕ + q χ , A t → A t + ∇ t χ to rewrite a scalar field φ = |φ|e iϕ in a gauge where the scalar field oscillates with a frequency ω, φ = |φ|e −iωt . However, since the energy-momentum tensor of the scalar field only depends on φφ † and ∂φ(∂φ) † , the gravitational and Maxwell fields are always invariant under the action of the Killing vector field ∂ t . In this gauge, the equations of motion require that A r+ = ω/q, if the solution has a horizon at r = r + . In particular, we shall adopt this gauge in sections III B and III C.
At this point we have fixed all the diffeomorphic and U (1) gauge freedom to simplify our field ansatz. However, our system still has two scaling symmetries that we can use for further simplifications 6 . The first scaling symmetry is {t, r, θ, ψ} → {λ 1 t, λ 1 r, θ, ψ}, {f, g, A, φ} → {f, g, A, φ}, {q, L, r + } → q λ 1 , λ 1 L, λr + , which leaves the equations of motion invariant and rescales the line element and the gauge field 1-form as ds 2 → λ 2 1 ds 2 and Adt → λ 1 Adt. We use this scaling symmetry to work with dimensionless coordinates and measure our thermodynamic quantities in units of L, the natural scale of AdS (this amounts to set L ≡ 1), The second scaling symmetry is A Taylor expansion of the equations of motion at the asymptotic boundary yields f c 0 (R 2 + 1) + · · · and g −1 R 2 + 1 + · · · . We use the scaling symmetry (9) to require that g approaches 1/f as r → ∞, i.e. to fix c 0 ≡ 1. This choice can also be motivated by the following observation (see e.g. [41]). In the context of the AdS 4 /CFT 3 duality, AdS 4 black holes are dual to thermal states on the boundary CFT. The choice g = 1 f as r → ∞ fixes the normalization of the time coordinate with respect to the gravitational redshift normalization of the radial coordinate to be such that the Hawking temperature of the black hole in the bulk matches the temperature of the dual CFT on the boundary.
Variation of (5) yields the equations of motion for the fields f, g, A and φ: We can use (10c) to get an algebraic relation for g in terms of the other fields and their first order derivatives: We insert this algebraic relation into (10a), (10b) and (10d) to obtain a coupled system of three second order ODE's. Onwards, these are the three equations of motion that we will solve to find the three fields {f, A, φ} (and thus also g).
The most notable solution of (5) is global AdS 4 spacetime: We are interested in solutions that asymptote to this background.
We can rewrite our ansatz and equations of motion in Fefferman-Graham coordinates (defined such that g zz = 1/z 2 and g zb = 0, where z is the radial distance with boundary at z = 0, and x b = {t, θ, ψ} are the boundary coordinates) [55][56][57]. Any asymptotically AdS 4 spacetime has the following Taylor expansion of the metric around the holographic boundary z = 0 [58,59]: ab (x) (12) and where g (0) (x) and g (3) (x) are the two integration "constants" of the expansion 7 ; the first dots include only even powers of z (smaller than 3) and depend only on g (0) (thus being the same for any solution that asymptotes to global AdS 4 ) while the second dots depend on the two independent terms g (0) , g (3) . Within the AdS/CFT duality we are (typically) interested in Dirichlet boundary conditions (BCs) that do not deform the conformal metric g (0) : the dual holographic CFT is formulated on this fixed background. This is the static Einstein Universe R t × S 2 , ds 2 ∂ = −dT 2 + dΩ 2 (2) . Given this Dirichlet BC our task is then to solve the equations of motion in the bulk, subject to regular BCs at the horizon or radial origin, to read g (3) . Using the holographic renormalization formalism, this gives us the expectation value of the holographic stress tensor T ab (x) which specifies and describes the boundary CFT [58][59][60]. This Dirichlet BC also ensures that there is no dissipation of energy at the asymptotic boundary (see Appendix A of [8]). For this reason these Dirichlet BCs can be denoted as 'reflecting boundary conditions' 8 .
We still need to discuss the asymptotic BCs for the Maxwell and massless scalar fields. A Taylor expansion of the equations of motion at the conformal boundary yields, where µ and ρ are the chemical potential and charge density, respectively, and ε and σ are the two arbitrary integration constants describing the decay of the massless scalar field. In d = 4 the unitarity bound for the scalar field mass is m 2 unit L 2 = −5/4. We have chosen to work with a massless scalar field. Hence it is above the unitarity bound. It follows that only the mode ε/r ∆+ with faster fall-off is normalisable (since it has finite canonical energy) and σ/r ∆− is a non-normalisable mode [43,61]. We thus choose the BC σ = 0. According to the AdS/CFT dictionary, the boundary CFT is then not sourced and the scalar field amplitude ε is proportional to the expectation value O of the boundary operator O that has dimension ∆ + = 3 [62].
To summarize, at the asymptotic boundary (R → ∞) we will impose Dirichlet BCs such that our solutions at large R behave as with the dots representing terms that are a function of m 0 , µ, ρ, ε.
Consider now the inner boundary of our solution. This is a fictitious boundary since it is a coordinate singularity at the edge of our integration domain, but is at a finite proper distance from other points [14]. In the present study, this can be the origin of the radial coordinate, R = 0, if we look into a soliton (boson star). Or, it can be a horizon at R = R + if the solution is a black hole, a condition that is imposed by demanding that f (R + ) = 0. In both cases to find the physical BCs we can Euclideanise the metric with a Wick rotation and require regularity of the metric and matter fields in Cartesian coordinates, as explained in detail in the review [14] For a soliton, a Taylor expansion of the equations of motion at the origin yields, where {f 0 , a 0 , φ 0 } are arbitrary integration constants and all higher order terms in (15) are even powers of R with coefficients fixed in terms of {f 0 , a 0 , φ 0 }. Here we have already imposed smoothness of the fields at the origin by imposing a defining BC (in the nomenclature of [14]). Namely, we set to zero three integration constants (one for each field f , g, and φ) which are associated to terms that diverge as R → 0. Following the nomenclature of [14], if we wish, we can now use this defining BC together with the equations of motion to impose a derived BC to find the 7 Our solutions are static and spherically symmetric so there is no x dependence. 8 We can however have other reflecting BCs that preserve the charges but not the static Einstein Universe conformal boundary. 9 For the gauge field, this refers to the gauge invariant field strength tensor F = dA, rather than the gauge potential A soliton. For example, it follows from (15) that a good derived BC is a Neumann BC for the three fields at R = 0 since the derivative of all of them vanishes.
On the other hand, if the solution is a black hole a Taylor expansion around its horizon − defined as the locus after imposing defining BCs that set three integration constants to zero (one for each function) to have smoothness as R → R + . We are left with the remaining three arbitrary integration constants {f 0 , a 0 , φ 0 }, and all higher order terms in (16) are fixed as a function of {f 0 , a 0 , φ 0 }. At this stage we can ask how many parameters we need to describe the soliton and hairy black hole that we want to construct. The theory (5) is fixed once we choose the mass and charge of the scalar field. At the UV asymptotic boundary we have a total of five free parameters {c 0 , m 0 , µ, ρ, ε} to which we need to add the cosmological radius L and the horizon radius R + (if present). As explained previously, we can use the two scaling symmetries (7) and (9) to fix L and c 0 . We are therefore left with four UV free parameters {m 0 , µ, ρ, ε} and the horizon radius R + . On the other hand, in the IR we have a total of three free parameters {f 0 , a 0 , φ 0 } (both in the soliton and black hole case). Therefore, any black hole of the theory is described by 4 U V + 1 H − 3 IR = 2 parameters. The AdS 4 −Reissner-Nördstrom black hole is indeed a 2-parameter family of black holes. More importantly, we also find that the hairy black holes we want to construct are described by two parameters. These are the mass and the electric charge of the black hole. Alternatively (but equivalently) we can take these to be the horizon radius R + and the amplitude of the scalar condensate ε. On the other hand, any soliton of the theory is described by 4 U V − 3 IR = 1 parameter. This can be the mass of the soliton (which fixes its charge) or, alternatively, the amplitude of the scalar condensate ε.

B. Conserved and thermodynamic quantities
Once we have found our fields {f, g, A, φ} we can read the gauge invariant thermodynamic quantities that describe hairy black holes or solitons.
To find the energy we can use either holographic renormalization [58][59][60] or the Ashtekar-Das formalism [63], [64]. We describe briefly the latter formalism. Energy is a conserved asymptotic quantity associated to the Killing vector field ξ a ∂ a = ∂ t . Consider a conformal transformation of the metric (6) defined as :ĝ ab = Ωg ab . The appropriate conformal factor to use here is Ω = 1 R . Then consider the Weyl tensor of the conformal metric:Ĉ abcd . Consider also the conformal hypersurface defined by Σ| Ω=0 and the normal vector to it: n a = ∂ a Ω. We define the tensorŝ K abcd = Ω 3−dĈ abcd | Ω=0 andε ac =K abcd n b n d . Consider, on the conformal hypersurface Σ| Ω=0 , the timelike surfaces of constant t,Σ| Ω=0,t=const with normal vector t a = ∂ a t. This surface has dimension 2 and induced metric h ab . Then the conserved energy is defined as: which, in terms of the expansion parameters defined in (14) yields, To define the conserved electric charge, we use Gauss' law. Consider a 2-sphere S 2 at constant time t and R in the limit R → ∞. The electrical charge is where F is the Hodge dual of the field strength. As a function of the asymptotic expansion parameters defined in (14) this gives, The value of the temperature for a soliton is undefined, and its entropy is zero as it does not have a horizon. Our black holes have a Killing horizon at R = R + generated by the Killing vector field K = ∂ t . Their temperature is thus given by its surface gravity over 2π which yields, The entropy of the black hole is a quarter of horizon area. In terms of the expansion parameters introduced in (16), the temperature and entropy are then The chemical potential µ is defined in (14). It is the difference between the value of the electromagnetic field at infinity and at the horizon 10 .
To discuss the competition between thermal phases in the canonical and in the grand-canonical ensembles we have to compute the associated thermodynamic potentials, namely the Helmholtz free energy F = E − T S and the Gibbs free energy, G = M − T S − µQ, respectively.
Solutions of (10) must satisfy the first law of thermodynamics: When the scalar field vanishes, global AdS 4 −Reissner-Nördstrom (RN-AdS 4 ) black hole is a known solution of (10). This is a 2-parameter family of solutions, that we can take to be the horizon radius R + and the chemical potential µ. It is described by ansatz (6) with Alternatively, we can parametrize RN-AdS 4 with R + and the temperature T in which case the chemical potential is given by It follows that regular RN-AdS 4 black holes exist if their chemical potential satisfies When µ = 0 one recovers the AdS 4 -Schwarzschild solution while the upper bound of (27) occurs for the extremal configuration (T = 0) where the entropy is finite but the temperature vanishes.
The mass, charge, entropy and temperature expressed as a function of R + and µ are In Section V, we will find the leading order thermodynamics of hairy black holes assuming that they can be described as a non-interacting mixture of a soliton and a small RN-AdS 4 black hole. To complete that analysis, we will only consider the leading order (in the R + expansion) thermodynamics of the RN-AdS 4 black hole. For future reference this is: When perturbed by an infinitesimal scalar field perturbation, global AdS 4 -RN black holes can develop two types of linear instabilities that have different physical nature. One is the near-horizon scalar condensation instability and the other is the superradiant instability. The former is already present in planar RN-AdS 4 black holes (where it was first found) while the later is only present in global AdS RN black holes. The zero mode of these instabilities is closely connected to the existence of hairy black hole solutions. Therefore, in the next two subsections we discuss the main properties of these two instabilities with some detail.

B. Near horizon scalar condensation instability of large black holes
Consider a scalar field with mass m and charge q in an asymptotically AdS 4 background. A Taylor expansion around the asymptotic boundary yields the two independent solutions 11 which reduces to (13) when m = 0. Breitenlöhner and Freedman (BF) found that a scalar field in AdS 4 is normalizable, i.e. it has finite energy, as long has its mass obeys the AdS 4 BF bound, When this is the case we can have asymptotically AdS 4 solutions that are stable in the UV region. But this analysis is blind to the behaviour of the scalar field in the interior of the spacetime. In particular, the properties of the scalar field in the IR region of a particular spacetime might drive the system unstable. This is indeed the case: as first found in the context of holographic superconductors [39][40][41], planar AdS black holes can be unstable to a near horizon scalar condensation mechanism [39]. The zero mode of this instability then signals a second order phase transition to planar hairy black holes that are dual to holographic superconductors [40,41]. It was later understood that the near horizon scalar condensation instability is quite universal and present for any asymptotically AdS background that has a black hole with an extremal (zero temperature) configuration [41,42] (this includes even some systems with neutral black holes and uncharged scalar fields). The criterion for this instability is the following. Start with a scalar field whose mass obeys the UV BF bound (31). When it violates the AdS 2 BF bound, the system, although asymptotically stable, is unstable in its near-horizon region. We might then expect the system to evolve into a hairy black hole with the same UV asymptotics but different IR behaviour.
We can work out in more detail the consequences of this instability for our current system, following the original analysis in [42] and [25]. Global AdS 4 -RN black holes described by (6) and (25) and let λ → 0. Setting, this procedure yields the solution of (10), So, the near-horizon limit of extremal RN-AdS 4 is the direct product of AdS 2 with a sphere S 2 . Applying this nearhorizon limit to the Klein-Gordon equation that describes a perturbation δφ(t, r) = e −iωt R(r) of a scalar field (with mass m and charge q) about the extremal RN-AdS 4 black hole yields which, not surprisingly, is the Klein-Gordon equation for a scalar field around AdS 2 with an electromagnetic potential A τ = αρ. A Taylor expansion of (36) yields which dictates the effective mass of the scalar field from the perspective of a near-horizon observer, where in the second line we re-introduced the dimensionless quantities R + = r + /L and e = qL. A scalar field with mass (38) in AdS 2 is unstable if it violates the AdS 2 BF bound (32). It follows that extremal RN-AdS 4 black holes are unstable whenever the charge of the scalar field obeys Taking the large R + expansion, this yields In the strict limit R + → ∞, we recover the known result for planar RN-AdS 4 black holes, e 2 ≥ 3 2 + m 2 L 2 [39][40][41]. Global RN-AdS 4 black holes are unstable for larger values of the scalar field charge.
In the opposite R + → 0 limit, an expansion of (39) yields which indicates that the near-horizon scalar condensation instability is suppressed for small (global) RN-AdS 4 black holes. This is to be contrasted with the superradiant instability (discussed in the next subsection) which is present for small black holes but not large black holes. This property highlights the different nature of the two instabilities. The near-horizon instability criterion that we have been discussing applies strictly to extremal black holes. Continuity indicates that this instability should extend to near-extremal black holes. Thus large, near extremal RN-AdS 4 black holes are unstable to condensation of the scalar field when condition (39) is satisfied.
We have chosen to work with a massless (m = 0) scalar field, so the near horizon instability is present in our system for large black holes with e > 3 2 .

C. Normal modes of AdS4
Consider a massive charged scalar field perturbation φ(t, r) in global AdS 4 with a constant gauge field A = µdt, which is governed by the Klein-Gordon equation. The background is time independent so we can do a Fourier decomposition of the perturbation in time such that a particular mode is described by φ(T, R) = e −iωT ψ(R) (ω is the frequency of the mode). If the scalar field mass m is above the BF bound (31), we can impose the boundary condition σ = 0 in (30) and this yields the normalizable solution where ∆ + is defined in (30), ε is an arbitrary amplitude and 2 F 1 is the hypergeometric function. To have a regular solution at the origin, the frequency spectrum must then be quantized as 12 where the non-negative integer n is the radial overtone that gives the number of radial zeros of ψ(R). Equations (42) and (43) give, respectively, the eigenvectors and eigenvalues of the normal modes of a scalar field. Note however that we can use the U (1) gauge transformation of the system to choose the gauge potential to be such that, e.g. µ = (∆ + + 2n) /e and the normal mode frequency (i.e. for a given n) vanishes, ω n = 0. In our current study, we are interested in the case m = 0, which implies that ∆ + = 3, and onwards we choose to work with the solution that has the lowest frequency (and thus energy), namely we set n = 0. Then, (43) reads

D. Superradiant instability of small black holes
In subsection III B we have seen that large RN-AdS 4 black holes are unstable to the near-horizon scalar condensation instability. In the present subsection we show that, in the opposite limit, small RN-AdS 4 black holes are unstable to the superradiant instability. For that we solve the Klein-Gordon equation for a massless scalar field perturbation with charge e ≥ 3 √ 2 (see discussion associated with (4)), around the RN-AdS 4 in a perturbative expansion in the horizon radius R + . We will find that the imaginary part of the perturbation frequency ω is positive, and thus the system is unstable, when µ ≥ 3 e . To impose smoothness of the perturbations at the horizon boundary it is convenient to work in ingoing Eddington- The Klein-Gordon equation cannot be solved analytically at each order. To circumvent this limitation, we consider small black holes, R + 1, and use the matched asymptotic expansion method to solve the resulting Klein-Gordon equation at each order in R + . Namely, we consider two different regions: the far region, R R + , and the near region, R 1 which have an overlapping region R + R 1 since we are assuming R + 1. In the near and far regions, certain contributions in the equation of motion are sub-dominant and can be discarded. We are left with ODEs, with a source term that depends on the previous order solutions and their derivatives, that can now be solved analytically. At each order, the near and far region solutions are then matched in the overlapping region.
Both in the far (superscript out ) and in the near (superscript in ) regions, we write the perturbative expansion in R + of φ(R) and ω as: and solve the Klein-Gordon equation order by order. The choice of ω 0 is justified because the perturbation mode of RN-AdS 4 with lowest frequency reduces to the lowest normal mode of AdS (44) in the limit R + → 0. At any order k ≥ 1 we have five parameters to determine: two integration constants from φ out k (R), say α k andα k , two integration constants from φ in k (R), say β k andβ k , and the frequency coefficient ω k . For our purposes, it is enough to present the results up to order k = 2: this is the order at which the frequency acquires an imaginary contribution.

Far region analysis:
In the far region, R R + , the Klein-Gordon equation for φ out (R) around the RN-AdS 4 black hole (25) effectively describes linearized perturbations around global AdS. Indeed at each order in the R + expansion it reads, whereD µ =∇ m u − ieA µ is the gauge covariant derivative of the global AdS backgroundḡ with a gauge potential A = µdt, and the RHS is a source term that depends on the lower order fields φ out j<k and their derivatives. As justified when discussing (13), at any expansion order, in the far region we impose the Dirichelet boundary condition: At leading order, k = 0 the far region solution is The boundary condition (49) requires we setα 0 = 0. Moreover, without loss of generality we can set α 0 ≡ 1, and (49) and (50) then fix the amplitude of the scalar field to be ε = −ie − 1 2 iπeµ . For higher orders, k ≥ 1, the boundary condition (49) always fixes the two integration constants α k andα k but leaves the frequency coefficient ω k undetermined. It is fixed by matching the far region with the near region. The final far region solutions φ out k for k = 0, 1, 2, after imposing the boundary condition (49) and fixing the frequency coefficients, are listed in (A1) of Appendix A.

Near region analysis:
In the near region R 1, it is convenient to work with the rescaled variable y = R R+ . The reason for this rescaling will become fully clear in the discussion associated to (67) and (68) that we avoid repeating here. In short, in the near region the perturbation problem simplifies since, for R + 1, the electric field is weak and the perturbation system reduces to a small perturbation around the neutral solution. After this rescaling y = R R+ , at each order, the Klein-Gordon equation for φ in (y) around the AdS 4 black hole − see (68) − reads whereD µ is the leading order part of the gauge covariant derivative of the rescaled coordinate in the RN-AdS 4 black hole 13 , and the RHS is a source term that depends on the lower order fields φ in j<k and their derivatives. At each order, we solve (51) subject to the boundary condition that the solution is regular at the horizon in Eddington-Finkelstein coordinates. At leading order, k = 0, the near region solution is and smoothness requires that we setβ 0 = 0. At any order k ≥ 1, the boundary condition at the horizon always fixes one of the integration constants, sayβ k , but the other is left undetermined until we do the matching of the near and far regions. The final near region solutions φ in k for k = 0, 1, 2, after imposing the boundary condition and fixing the integration constants, are listed in (A2) of Appendix A.
Matching the near and far region solutions: At this stage, at each order k ≥ 1 we have two undetermined parameters, namely the near region integration constant β k and the frequency coefficient ω k . They are fixed by requiring that the large radius expansion of the near solution φ in matches the small radius expansion of the far solution φ out . In these expansions we keep only terms that will not receive contributions from higher orders.
Let us illustrate this matching at the leading order, k = 0. The small radius limit of φ out 0 and the large radius limit of φ in 0 yield Matching requires that we set β 0 = 1.
As another example consider the matching at order k = 1. The small radius expansion of φ out 1 and the large radius expansion of φ in 1 are: At order O(R 0 + ) the two solutions agree as a result of the previous order matching. Moving to next-to-leading order, the divergent R + /R term in the far region has no counterpart in the near region expansion so we fix ω 1 to eliminate it. Matching the remaining terms, namely the contributions proportional to R + R 0 and R + log R then fixes uniquely the integration constant β 1 . A similar matching procedure can be done for higher orders.
By the end of the day, we find that the frequency coefficients in (47) up to order k = 2 are: The property in (55) to be highlighted is that at order k = 2 the frequency acquires an imaginary part proportional to µe − 3. This is negative whenever µ ≤ 3 e , which signals exponential damping, but it becomes positive for µ > 3 e . The latter case describes the superradiant instability of RN-AdS 4 black holes. The reader interested on the properties of RN-AdS 4 superradiant modes beyond the perturbative analysis done here can find them in [22,65].
The above analysis shows that global RN-AdS 4 black holes are unstable to superradiance in the small horizon radius regime, R + 1. On the other hand, large radius global RN-AdS 4 black holes are expected to be stable to superradiance. Indeed, in the most extreme case of a large radius black hole whereby R + → ∞, the studies of [39,40] indicate that planar RN-AdS 4 black holes are unstable only to the near-horizon instability described in Section III B. Starting from this planar limit, it seems natural to expect that, as R + decreases, the superradiant instability will kick in at (and below) a critical horizon radius (in addition, the near-horizon instability becomes supressed in the limit R + → 0 as seen in (41)). It would be interesting to do a numerical linear perturbation analysis that spans the full parameter space of RN-AdS 4 black holes to find the instability properties of RN-AdS 4 . However, it should be noticed that such an analysis will not be able to disentangle the nature of the two instabilities except in the two limiting cases R + → 0 and R + → ∞ that are covered by the analytical analysis done in this manuscript.

IV. SMALL SOLITONS (BOSON STARS)
In subsection III C we revisited the normal mode frequency spectrum of global AdS 4 . These frequencies are quantized according to (43). A far-reaching property of this spectrum is that the imaginary part of the frequency vanishes: the associated eingenmodes are not dissipative. Moving beyond linear order in perturbation theory in the amplitude of the scalar field, we might then consider the back-reaction of these scalar normal modes on the gravitational and electromagnetic field. Since the leading order is not radiative we might suspect that the normal mode can be backreacted to all orders and yield a boson star, i.e. a horizonless solution with static electromagnetic and gravitational fields and a complex scalar field with time dependence e −iωt . We should further expect that the frequency ω is given at leading order by (43) but gets corrected as we climb the expansion ladder. We can further explore the U(1) gauge transformation freedom to rewrite this solution in a gauge where the frequency vanishes and the scalar field is real. This static solution is equivalent to the boson star but usually called a soliton. As explained in section II, this soliton is a 1-parameter family of solutions that we can take to be the asymptotic amplitude ε of the scalar field as described in (14).
In this section, we consider the massless scalar field case m = 0, and will confirm these ideas. Namely, we construct the ground state soliton − i.e. the soliton with lowest energy that at leading order is described by the lowest normal mode (44) − up to fourth O(ε 4 ). A U (1) gauge transformation allows to set the frequency of the solution to zero; then at leading order we have A = 3 e dT . To get this soliton, we solve perturbatively the three second order ODEs for {f, A, φ}, that we obtain after replacing the algebraic equation (11) for g into (10), subject to the boundary conditions (14) and (15) at the asymptotic boundary and at the origin, respectively. The pair of integration constants associated to each ODE is fixed by the boundary conditions. We build up the perturbative solution in a power series expansion of the fields in ε around global AdS 4 , When ε = 0 we recover global AdS 4 as expected. The O(ε) scalar field and its derivatives source the first non-trivial contribution to the electromagnetic and gravitational fields at order O(ε 2 ). The structure of the equations of motion is then such that the odd (even) powers of ε never contribute to the electro-gravitational (scalar) fields.
The equations of motion can be solved analytically order by order for arbitrary order but the expressions for the solutions quickly become quite long. In Appendix B we display the final (i.e. after imposing the boundary conditions and thus fixing all integration constants) coefficient functions f 2n (R), A 2n (R) and φ 2n+1 (R) up to order n = 2. This order is enough to extract the relevant physical conclusions of our study.
Using these field expansions we can now compute the gauge invariant thermodynamic quantities discussed in Section II B. Namely, the mass M sol , charge Q sol , chemical potential µ sol , Helmoltz free energy F sol and Gibbs free energy G sol of the ground state soliton up to order O(ε 6 ) are: As a non-trivial check of our computation, we confirm that these quantities satisfy the first law of thermodynamics for solitons (24) up to the required order O(ε 6 ). We postpone the physical discussion of our soliton results (57) to Section VII.

V. NON-INTERACTING THERMODYNAMIC MODEL FOR HAIRY BLACK HOLES
In the last Section we have constructed the ground state AdS 4 soliton of the Einstein-Scalar-Maxwell theory (5). On the other hand, this theory also has the RN-AdS 4 black hole as a solution. In Section (III D) we have seen that the latter is unstable to superradiance when scattered by a scalar field. It is natural to expect that this superradiant instability drives the system to a hairy black hole that has a charged scalar condensate floating above the horizon, with the electromagnetic repulsion balancing the gravitational collapse of the scalar condensate into the horizon. For this to be true we need that: 1) a hairy black solution of (5) does indeed exist, and 2) that for a given energy and charge, the entropy of the hairy black hole is higher than the entropy of the unstable RN-AdS 4 black hole so that the system can evolve from the latter to the former solution while preserving the second law of thermodynamics.
In this section we confirm the above expectations, i.e. we establish the existence of the hairy black hole solutions and their leading order thermodynamic properties using a very simple non-interacting thermodynamic model introduced in [23][24][25]. In the next section we will confirm that this thermodynamic model yields the correct physical properties of the sytem by solving directly the equations of motion (10)-(11) as a boundary value problem. Our analysis thus complements a recent numerical time evolution simulation which, for some solution parameters and initial values, showed that RN-AdS 4 are superradiantly unstable and decay into hairy black holes [28].
The leading order thermodynamics of the hairy black hole can be constructed using the so-called non-interacting thermodynamic model of [23,24] (see also [8,11,13,25] ). One interprets the hairy BH as a non-interacting mix of a soliton and a RN-AdS 4 BH. One first observes that at leading order a charged boson star (soliton) is just a normal mode of the ambient background (whose frequency ω is then corrected as we walk along the perturbation expansion ladder). This is a 1-parameter family of solutions with E = ω q Q + O(Q 2 ), where q is the charge of the scalar condensate. This soliton further obeys the first law, dE = ω q dQ. We can now place a small RN-AdS 4 BH on top of this soliton to get a 2-parameter family of hairy BHs. At leading order, we will assume a non-interacting mixture whereby we take the mass M and the charge Q of the resulting hairy BH to be just the sum of the masses and charges of the BH (M RN and Q RN ) and soliton (M sol and Q sol ) components: 14 The soliton carries no entropy, so the entropy S of the hairy BH simply reads The hairy BH can partition its charge Q and mass M between the RN-AdS 4 BH and soliton components. On physical grounds one expects this distribution to be such that, for fixed mass M and charge Q, the entropy S is maximised, dS = dS RN = 0, while respecting the first law of thermodynamics, For fixed total mass and charge, dM = 0 = dQ, and thus dM RN = −dM sol and dQ RN = −dQ sol . It follows from the first law for the soliton that the LHS in (60) yields the chemical potential of the soliton, dM RN dQ RN = µ sol . On the other hand, the first law for the RN-AdS 4 BH implies that ∂ M RN S RN = 1/T H and ∂ Q RN S RN = −µ RN /T H where T RN = T H is the temperature of the RN-AdS 4 BH (and thus of the hairy BH) and µ RN is its chemical potential. Therefore the RHS of (60) is simply − The simple non-interacting thermodynamic model and associated maximization of entropy therefore requires That is, the hairy BH inherits the chemical potential of its solitonic component and the system is in chemical and thermal equilibrium. The mass and charge of the soliton are related by M sol = µQ sol , while the leading order thermodynamics of the RN-AdS 4 black hole is given by (29). Using these properties together with the non-interacting relation (58) and equilibrium conditions (61), we can express the leading order thermodynamic quantities of the hairy black hole (and of its two components) in terms of its mass M , charge Q and chemical potential µ: The domain of existence of the hairy black hole can be inferred from this analysis. In one extremum, the soliton component is absent and all the mass and charge of the hairy BH is carried by the RN-AdS 4 component. This describes the hairy black hole that merges with the RN-AdS 4 black hole at the zero-mode of its superradiant instability. On the opposite extremum configuration, the RN-AdS 4 component is absent and the soliton component carries all the mass and charge of the solution. This is the zero-radius or zero-entropy limit of the hairy black hole. It follows that the hairy black hole mass must be within these two boundaries: Equating the two extrema configurations, yields an interval of existence for the chemical potential: where in the last relation we used that for the ground state solution the leading order potential is µ = ω/e = 3/e + O(R + ). Note that condition (64) is the same we found in the superradiant analysis leading to (4). The above leading order thermodynamic analysis must be considered with a few grains of salt. Indeed, first note that a theory can have hairy black holes that do not have a zero-radius limit, i.e. a solitonic limit. Second, there is no reason why the non-interacting mixture assumption (58) should hold, even at leading order.
In the next section we will construct the hairy black hole. In Section VI C we compute the thermodynamic quantities of the hairy black hole and we check that our intuition was correct, and that, to leading order, the hairy black hole thermodynamic quantities are indeed given by (62).

VI. SMALL HAIRY BLACK HOLES
In this section we construct perturbatively the hairy black holes whose leading order thermodynamics was discussed in the previous section. For that we solve the elliptic equations of motion (10)- (11), subject to the boundary conditions (14) and (16), to find a perturbative analytical expression for the hairy black hole.
A. Setting up the perturbation problem As explained in Section II, the hairy black holes of our theory are a two-parameter family of solutions that we can take to be the asymptotic scalar amplitude ε and the horizon radius R + . Thus, their perturbative construction requires that we do a double expansion of the fields in powers of ε and R + . In order to be able to solve analytically the equations of motion (10)- (11) we must resort to a matched asymptotic expansion, similar to those done in a similar context in [11,12,[23][24][25]. Much like in Section III D, we divide the outer domain of communications of our black hole into two regions; a near-region where r + ≤ r L and a far-region where r r + . Restricting the analysis to small black holes that have r + /L 1, the two regions have an overlapping zone, r + r L. In this overlapping region, we can match/relate the set of independent parameters that are generated by solving the perturbative equations of motion in each of the two regions.
The chemical potential of the solution should itself have a double expansion in powers of ε and R + , Indeed, recall that the soliton is the back-reaction of a normal mode of AdS to higher orders. At leading oder the chemical potential of the soliton is related, via a gauge transformation, to an AdS normal mode frequency. We saw that this is corrected at higher orders. Thus, we must permit similar corrections when the horizon is present. We shall construct the hairy black hole family whose zero-radius limit is the ground state soliton of Section IV (so with lowest energy for a given charge). 15 In the far region, R R + , the hairy black hole can be seen as a small perturbation in R + and ε around global AdS. We use the superscript out when referring to far region fields which have the double expansion: This expansion already anticipates that odd (even) powers of the scalar condensate do not correct the fields f, A (φ). At each order {n, k} the perturbed equations of motion can be solved analytically (with the help of Mathematica) to find the far fields up to a total of 6 integration constants (recall that the equations of motion are a system of 3 second order ODEs). In addition, the system also depends on the chemical potential corrections µ 2n,k . The requirement that the far-region fields obey the asymptotic boundary conditions (14) typically fixes 4 of the integration constants (namely, those associated to φ out 2n+1,k , one associated to f out 2n,k and another to A out 2n,k ). We are left with two integration constants and µ 2n,k that will be determined by the requirement that the far fields match the near-region fields in the overlapping region. To prepare the system for the matching we need to take the small radius R expansion of f (out) , A (out) , φ out (with the expansion coefficients available at the order under consideration). We find that these are singular as R → 0, diverging with a power of R+ R . This indicates that the far-region analysis breaks down at R ∼ R + . This justifies why the far-region analysis is valid only for R R + . Also, it follows that in the far-region we can safely do a Taylor expansion in R 1 and ε 1 since the large hierarchy of scales between the solution parameters and the distance guarantees that they do not compete.
Consider now the near-region, R + ≤ R 1. This time the Taylor expansions in R 1 and ε 1 should proceed with some caution since the small parameters can now be of similar order as the radius R + . This is closely connected with the fact that the far-region solution breaks down when R/R + ∼ O(1). This suggests that to proceed with the near-region analysis we should define a new radial, y, and time, τ , coordinates as 15 A similar construction could be done for the excited hairy black holes.
The near region now corresponds to 1 ≤ y R −1 + . If we further require that R + 1 one sees that the near region corresponds to y ≥ 1 R + (and y ε). In particular, we can now safely do Taylor expansions in R 1 and ε 1 since the radial coordinate y and the black hole parameters have a large hierarchy of scales 16 . To have further physical insight it is also instructive to rewrite the RN-AdS 4 solution (25) in the new coordinate system (67) The explicit factor of R + 1 in A τ (y) indicates that in the near region the electric field is weak and the system can be seen as a small perturbation around the neutral solution. The same should hold when we add a small scalar condensate to the system. The near fields of the hairy black hole thus have the double expansion, At each order {n, k}, we can solve the perturbed equations of motion analytically to find the near fields up to a total of 6 near region integration constants. Imposing the horizon boundary conditions (16) typically fixes 3 of the integration constants (namely, one associated to φ in 2n+1,k , one to f in 2n,k and another to A in 2n,k ). The 3 leftover near field integration constants are left undetermined until we match the far and near region fields in the overlapping region. For the matching, we first need to restore the original coordinates {T, R} and take the large radius limit of f in (R), A in (R), φ in (R) (with the expansion coefficients available at the given order). We find that these diverge as a power of R, which shows that the near region analysis breaks down at R ∼ 1. This explains why the near region analysis is valid only for R 1. In the following subsection, we give some low order detailed examples of the matching asymptotic expansion implementation. We are necessarily constrained on this exposition to lowest orders since, as the expansion order grows, the analytical expressions for the far and near fields become quite large. For the benefit of the reader that wants to reproduce our results in full, we give the far and near field functions in Appendices C 1 and C 2, respectively. These are the final fields, i.e. after fixing all the integration constants of the system using the boundary conditions (14) and (16) and the matching of far and near fields.
We leave to subsection (VI C) in the main text, the physically relevant outcome of our matching asymptotic construction, namely the thermodynamic quantities that describe uniquely the hairy black holes. At lowest order, O(ε 0 ) in the scalar amplitude expansion, the scalar field is naturally absent and the far field coefficients {f out 0,k (R), A out 0,k (R)} can be read directly from an R + 1 expansion of the RN-AdS 4 solution (25). A similar Taylor expansion of (68) yields the near field coefficients {f in 0,k (y), A in 0,k (y)}. The O(ε 1 ) correction turns-on the scalar field φ without back-reacting yet in the gravitational background: it describes a small perturbation of the scalar field around the RN-AdS 4 black hole. That is, the non-trivial equations of motion (10)- (11) reduce to the Klein-Gordon equation without a source 17 .
Next, we illustrate how the matching asymptotic expansion determines φ 1,k . A similar procedure yields the fields φ 2n+1,k at the odd orders O(ε 2n+1 ).
• Far region, R R + : The definition of ε introduced in the boundary condition (14) determines completely the solution at order O(R 0 + ): 18 and the equation of motion is obeyed only for the choice µ 0,0 = 3/e (which corresponds, via a gauge transformation, to the lowest normal mode frequency of AdS 4 ). We now expand the equation to O(R 1 + ). We get a Klein-Gordon equation with an inhomogeneous term sourced by a derivative of φ out 1,0 (R) which can be solved analytically. The two integration constants are directly determined by the asymptotic boundary condition (14) (see also footnote 18). This yields The value of the chemical potential correction µ 0,1 will be determined only in the matching of the far region with the near region.
As explained above, we skip here the details of the computation of the field coefficient φ out 1,2 (R). However, it is perhaps useful to highlight that it is typically a good idea to complete the matching procedure, that fixes constants not determined by the boundary conditions, before proceeding to the next order. This keeps the size of the expressions that carry to next order smaller.

The homogeneous Klein Gordon equation in the near region at
We need to impose the boundary conditions (16) at the horizon y = 1. Namely regularity of the scalar field at the horizon requires that we eliminate the divergent logarithmic term by choosing C 2 = 0. Hence, the regular near region solution is At order O(R 1 + ), φ in 1,1 (y) still solves the homogeneous Klein-Gordon equation 19 . Therefore, its regular solution is • Matching in the overlapping region R + R 1: In order to do the matching we take the small R limit of the far region solution and the large R limit of the near solution. Notice that at this order only terms up to R 0 R 1 + and R 1 R 0 + are correctly accounted for (terms higher than this will receive corrections from the next order). The small R expansion of the far region is: 17 At higher orders in n, the equation of motion for φ still has the form of a Klein-Gordon equation but with an inhomogeneous term sourced by the lower order fields and their derivatives. 18 To be precise, ε is an integration constant that parametrizes our solution and, without loss of generality, we choose not to allow corrections to its definition (14) at any order. 19 A would be source term is proportional to the derivatives of the previous order field φ in 1,0 , but the latter is just a constant, see (73).
On the other hand, the large R expansion of the near region solution is: The divergent term in the far region expansion has no counterpart in the near region and needs to be eliminated by a judicious choice of µ 0,1 . Then, the matching of the two leftover terms fixes the integration constants C 1 and C 1 that were not fixed by the boundary conditions. Altogether we find, With this we have completed the solution up to O(ε 1 , R 1 + ). We can now proceed to order O(ε 1 , R 2 + ) following a similar procedure. The final far and near field coefficients {f out 0,2 , A out 0,2 , φ out 1,2 } and {f in 0,2 , A in 0,2 , φ in 1,2 }, after imposing the boundary conditions and the matching are listed in Appendices C 1 and C 2, respectively.
2. Matching asymptotic expansion at O ε 2 , R k + At order O ε 2 , the order O ε 1 scalar field back-reacts on the metric and gauge potential and the hairy black hole starts differentiating from the RN-AdS 4 solution. On the other hand the Klein Gordon equation is trivially satisfied. Next, we illustrate how the matching asymptotic expansion determines f 2,k and A 2,k . The procedure applies similarly for all the even orders ε 2n .
• Far region, R R + : Starting at order O ε 2 , R 0 + the far field equation of motion can be solved analytically to get f out 2,0 and A out 2,0 . Each one of these fields has two integration constants. Imposing the asymptotic Dirichlet boundary condition (14) we fix one of the constants in each pair. We are left with the integration constants C f and C A that will be fixed in the matching step. At this stage the fields read: • Near region, R + ≤ R 1. Also at order O ε 2 , R 0 + , we solve for the near fields f in 2,0 and A in 2,0 up to a pair of integration constants. Requiring that these fields vanish linearly at the horizon y = 1, as required by the horizon boundary condition (16), we fix two of the constants. The near region fields are then still a function of two arbitrary constants K f and K A : • Matching in the overlapping region R + R 1: To match the far and near fields, we take the small R expansion of the far fields (78): and the large R expansion of the near fields (80): A in Note that we keep only terms that will not be corrected by higher order contributions. Matching the two expansions fixes the four integration constants, that were not determined by the boundary conditions, as At this stage the chemical potential correction µ 2,0 is still left undetermined. It is fixed at order O ε 3 , R 0 + , when the scalar field φ 3,0 is found.

C. Thermodynamic quantities
The final results for the field coefficients introduced in (66) and (69) are given in Appendices C 1 (for the far region) and C 2 (for the near field). These are the final fields, i.e. after fixing all the integration constants of the system using the boundary conditions (14) and (16) and the matching of the far and near fields.
In the presentation of our results, we take 1 and R + 1 and we assume that O( 2 ) ∼ O(R + ). The latter assumption implies that terms with the same (n + k) contribute equally to the perturbative expansion, i.e.
We did the consistent perturbative expansion analysis up to (n + k) = 2, which is sufficient for our purposes. Extending the perturbative analysis to higher orders is possible, but increasingly cumbersome.
The fields presented in Appendix C allow to compute the relevant physical properties of hairy black holes. Namely, as described in II B, we can compute the thermodynamic quantities of the hairy solutions. The dimensionless mass M , electric charge Q, chemical potential µ, entropy S, temperature T are: and the Helmoltz free energy F and Gibbs free energy G are: In the next section we discuss physical checks to these quantities and withdraw physical conclusions.

A. Checking the hairy solutions and their interpretation
The thermodynamic quantities (87) for the hairy black hole must pass a few checks that already give valuable information about their physical interpretation.
Firstly, when we set the horizon radius to zero, R + → 0, we do get the thermodynamical quantities (57) of the soliton. Recall that the latter was obtained without any matching asymptotic expansion. This also confirms that the soliton is the zero horizon radius limit of the hairy black hole. We have assumed this was the case in our construction and the fact that the computation can be done is consistent with it.
Secondly, when we set the scalar condensate amplitude to zero, ε → 0, we get the thermodynamical quantities (28) of the RN-AdS 4 black hole, with a very specific chemical potential namely, That is, although the RN-AdS 4 black hole is a 2-parameter family of solutions, the ε → 0 limit of the hairy black hole gives a particular 1-parameter sub-family of RN-AdS 4 parametrized by R + and chemical potential fixed by (88). This is precisely the RN-AdS 4 1-parameter family at the onset of the superradiant instability studied in Section III D. To see this is indeed the case start by recalling that in Section III D we studied scalar field perturbations with Fourier time dependence φ(T, R) ∼ e −iωT φ(R) and found that the frequency is quantized as described in (47) and (55). In particular, the frequency has a positive imaginary part for ω > eµ signalling the superradiant instability. The onset of the superradiant instability occurs when the imaginary part of the frequency vanishes and we can work in the U (1) gauge where the real part of the frequency vanishes. If we impose these superradiant onset conditions on (55) and solve with respect to the chemical potential in a series expansion in R + up to order O R 2 + we precisely get (88). This is an important check to our computations and, together with the information from the last paragraph, confirms that hairy black holes are indeed a 2-parameter family of solutions that spans an area in a phase diagram that has the onset curve of the RN-AdS 4 superradiance and the soliton curve as boundaries. At leading order, these are the boundaries (63) predicted by the simple thermodynamical model of Section V and that will be very clear in the left panel of Fig. 1 (red and black curves).
The first law of thermodynamics provides the third check. Indeed, we can explicitly verify that the thermodynamic quantities (87) do obey the first law (23) up to order (n + k) = 2. We emphasize this is a fundamental and non-trivial check of our results.
To summarize, we have constructed analytically the hairy black hole within perturbation theory up to order (n+k) = 2 (and we could extend this construction to higher order). This was done solving directly the perturbed equations of motion (10)-(11) of the theory as a boundary value problem.
We can confirm that the simple non-interacting thermodynamic model of Section V − which does not use the equations of motion − yields the correct leading order thermodynamics, i.e. the leading terms of (87). First note that at leading order it follows from (87) that µ = 3 e . Assuming, as justified above, that O( 2 ) ∼ O(R + ), the leading order contributions of the expansion (87) allows to express analytically R + and ε in terms of M and Q, which we insert in the expressions for the other thermodynamic quantities to find that at leading order in M and Q one has: These quantities do match the result of the non-interacting model (62). This confirms that the non-interacting thermodynamic model, in spite of its crude simplicity, is quite robust and does indeed capture the fundamental leading order properties of the hairy black hole system at very low cost. This explicit confirmation adds to those done in similar hairy black hole systems in [8,11,13,[23][24][25].
A more detailed discussion of the regime of validity of our perturbation theory is also in order. By construction, it should be valid only for ε 1 and R + 1. On the other hand, the scalar charge of the hair must obey e 3 √ 2 ∼ 2.12 to have superradiant hairy black holes: see discussion that leads to (4). However, since the black holes are a small expansion around AdS 4 , e should not be too large to avoid large back-reactions. This suggests that it is appropriate to require at most 0.1, R + 0.1 and 3 √ 2 e 3, say. Inserting these bounds into the thermodynamic formulas (28), (57), and (87) we get approximate upper bounds for the physical charges. For example, for e = 2.5 we should look into masses and charges that, in AdS units, are themselves below 10 −1 . However the reader interested on making a direct comparison between our perturbative results and exact numerical constructions that emerge from a nonlinear code or from the endpoint of the superradiant time evolution might require more precise statements. We will have the opportunity to give precise criteria in the following discussion: see footnotes 21 and 25.
We have found that the gravitational Higgs model with action (5) admits several solutions. Namely, global AdS 4 , the RN-AdS 4 black hole, the ground state soliton, the ground state hairy black hole and an infinite tower of excited solitons and hairy black holes 20 . For a given electric charge, the latter excited solutions always have larger energy than the ground state partners so we do not discuss them further.
Naturally, these thermal phases of the theory compete with each other. This competition can be framed in one of the three possible ensembles: the microcanonical, canonical or grand-canonical ensembles. We discuss the phase diagram in these three ensembles in the next three subsections. Recall that the thermodynamic quantities of the RN-AdS 4 black hole, (ground state) soliton and hairy black hole are revisited in (28), (57), and (87), respectively [66,67].
Our results are better illustrated with some plots. Recall that the gravitational Higgs theory with action (5) is only fully defined once the particular value for the scalar field charge q = e/L is specified. For definiteness we will choose the particular value of e = 2.5 in our plots.

B. Phase diagram in the microcanonical ensemble
In the microcanonical ensemble the energy M and electric charge Q of the system are fixed and the the entropy S is the relevant thermodynamic potential to discuss the competition between the several thermal phases: the phase with higher entropy is the favoured one (global AdS 4 and the soliton have vanishing entropy).
In Fig. 1 we display the phase diagram M vs Q of the microcanonical ensemble for e = 2.5 (our perturbative expressions are valid for small M, Q). To make the diagram more clear, in the vertical axis we actually plot the difference between the mass of the solution and the mass of the extremal RN-AdS 4 with the same electric charge, ∆M ≡ M − M ext . So RN-AdS 4 black holes exist in regions I and II where ∆M ≥ 0. The black line with negative slope describes the soliton. Hairy black holes exist in regions II and III. Namely they fill the area limited by the soliton curve (where R + → 0) all the way up to the magenta line with positive slope (with ε = 0). The latter also describes RN-AdS 4 black holes at the onset of superradiance. So RN-AdS 4 black holes below (above) the magenta line are unstable (stable) to superradiance.
In region II there is no uniqueness since RN-AdS 4 and hairy black holes coexist with the same M and Q. To find the preferred phase in the microcanonical ensemble we must compare their entropy S. We find that for a given electric charge Q the hairy black hole always has higher entropy and is thus the favoured solution, whenever they coexist. This is illustrated with a particular example in the right panel of Fig. 1  Our findings also permit robust conclusions about the endpoint of the superradiant instability of the RN-AdS 4 black holes. Typically, numerical simulations are done at fixed energy and charge. Consider starting with a RN-AdS 4 that is within the curve AB and thus unstable to superradiance. By the second law of thermodynamics, the entropy of the system can only increase. Assuming that there is no other black hole solution in the spectrum of the theory besides the ones discussed so far, the system must evolve towards the hairy black hole that has the same M and Q but higher entropy. The latter can be pinpointed in (87) or in Fig. 1. By construction, this hairy black hole is stable to superradiance and we have no arguments suggesting that it is unstable to any other mechanism (see however footnote 1). Of course our findings say nothing about how the system actually evolves in time but the endpoint of the numerical simulation that gives this information can be tested against (87). 22 21 This provides a further check of our computations. Indeed, note that the first law requires that at a second order phase transition the slope dS/dM is the same for the two branches since T is the same at the bifurcation point (and dQ = 0 in the right panel of Fig. 1). This is clearly the case in our plot. However, if we start departing from the regime of validity of our perturbation analysis we increasingly find that the merger is not perfect and the slopes of the two branches no longer match. This is the best criterion to identify the regime of validity of (87). 22 Also recall that for 3/2 < e < 3/ √ 2 there are hairy black holes that emerge from the near-horizon scalar condensation instability that are not captured by our "superradiant" solutions (87). But they can be found as the endpoint of a time evolution simulation or by solving numerically the elliptic system of equations (10)- (11).

C. Phase diagram in the canonical ensemble
In the canonical ensemble the temperature and electric charge of the system are held fixed. The relevant thermodynamic potential in this ensemble is the Helmholtz free energy F = E − T S and the preferred phase is the one with lowest free energy.
Before discussing the hairy phase and thermal competition it is convenient to review the thermal properties of the RN-AdS 4 solution. These were discussed in detail in the seminal works of [66] and [67]. The left panel of Fig. 2 displays the temperature of RN-AdS 4 as a function of the horizon radius for four different values of the charge. This figure reveals two key properties (the upper, brown, curve describes the Schwarzschild-AdS 4 solution and is present only for reference). First, notice that there is a critical charge , Q crit /L = 1 6 √ 2 , associated to the dashed green curve. For charges above this (lower black curve) there is only one RN-AdS 4 black hole for a given {Q, T } pair. However, for Q > Q crit (second, blue, curve from the top) we see that there is a window of T where {Q, T } do not uniquely describe the RN-AdS 4 solution since we can have three different RN-AdS 4 branches that can be denoted as the 'small', 'intermediate' and 'large' black holes depending on their horizon radius (which is proportional to the square root of their entropy). 23 A second property observed in the left panel of Fig. 2 is that T → 0 as we decrease R + towards the extremal configuration (which is absent when Q = 0). Therefore, small horizon radius corresponds to small temperatures. 23 For reference, when Q = 0 (top brown curve in left panel of Fig. 2) we have not three but two black hole branches because T → ∞ as R + → 0. These are commonly designated by the 'small' and 'large' Schwarzschild-AdS 4 black holes. For Q ≤ Q crit , the first branch, that corresponds to smaller values of the horizon radius, goes from extremality to a maximum found at R + =  However, T = 0 is reached at a finite minimal radius that increases with Q. We will come back to this observation later when discussing the validity regime of our perturbation theory.
We can now introduce these RN-AdS 4 solutions in a phase diagram for the canonical ensemble and populate it with the novel hairy solutions. The canonical phase diagram consists of plotting the free energy F as a function of T for a fixed Q. 24 In Fig. 3 we show this phase diagram F (T ) for the three relevant cases: Q = 0, 0 < Q < Q crit and Q > Q crit (left, middle and right panels). For the phases without hair we keep the same colour code as in the left panel of Fig. 2. Moreover we necessarily do the plots for specific values of Q and e but they are qualitatively similar for all other values where the perturbation results (87) are valid. The left panel with Q = 0 describes the Schwarzschild (Schw-AdS 4 ) black hole and is given here again for a familiar reference [68]. We have the 'large' branch with negative specific heat (left) and the 'small' branch with positive specific heat (right) . The free energies of large Schw-AdS 4 are always lower than that of the corresponding small Schw-AdS 4 with the same T . Therefore, large Schw-AdS 4 are favoured over small Schw-AdS 4 black holes. In the phase diagram of the left panel of Fig. 2 the large and small Schw-AdS 4 branches meet at the regular cusp (where the specific heat vanishes). This is not the whole story since there is a third phase − thermal AdS 4 − which is just the Euclidean solution of global AdS 4 with an arbitrary period (and thus T ) chosen for the Euclidean time circle. For a temperature below the Hawking-Page T HP (defined as F (T HP ) = 0), thermal AdS 4 has lower free energy than both large and small Schw-AdS 4 black holes. However, at temperatures above T HP , large Schw-AdS 4 black holes are the preferred phase. This is the familiar Hawking-Page (HP) first-order phase transition [68]. In the AdS/CFT context, this is interpreted as a confinement/deconfinement transition in the dual conformal field theory [69].
Consider now the case 0 < Q < Q crit , whose phase diagram (for Q/L = 0.05) is plotted in the middle panel in Fig.  3. As explained above we have now three RN-AdS 4 branches. The large branch is the natural extension to Q = 0 of the large Schw-AdS 4 branch. However, the small branch now extends all the way down to the extremal, T = 0, solution (for Q = 0 it extends instead to T → ∞). These two branches are connected by a third, 'intermediate' branch. When they coexist small RN-AdS 4 black holes are the dominant branch below a critical temperature (determined by the intersection of the small and large branches in the middle panel). Above this critical temperature, large RN-AdS 4 black holes have lower free energy than the small and intermediate branches and thus dominate the canonical horizon radius. The intervals for the temperatures of the three RN-AdS 4 branches are: The critical charge Q crit is the case where the maximum and the minimum of T (R + ) coincide (inflection point). 24 Ideally we would show the 3-dimensional plot F (T, Q) but it is not very clear. The plots F (T ) for fixed Q (or F (Q) for fixed T ) illustrate better the relevant conclusions.
ensemble. Above T HP (defined such that F (T HP ) = 0) large RN-AdS 4 are the preferred phase, but below it there is a first order Hawking-Page phase transition and thermal AdS 4 is the favoured phase for all T < T HP . In particular, as T decreases to T = 0 and small RN-AdS 4 becomes the only branch, thermal AdS 4 is still the preferred phase. When Q = 0 the theory also has hairy black hole solutions. In the perturbative regime where our results are valid 25 , the hairy family bifurcates from the branch of small RN-AdS 4 at the onset of superradiance and extends to larger values of the temperature all the way up to the soliton limit (R + → 0) that is reached at a finite value of T . 26 However, the free energy of the hairy black hole (and soliton) is always larger that the free energies of small and/or large RN-AdS 4 and of thermal AdS 4 . Therefore, hairy black holes and the soliton are never the preferred phase in the canonical ensemble. This is thus in sharp contrast with the situation in the microcanonical ensemble.
Finally we have to discuss the case Q > Q crit (right panel in Fig. 3). Recall from the discussion associated with the left panel of Fig. 2 that in this case there is a single branch of RN-AdS 4 black holes. Its free energy is plotted in the right panel in Fig. 3. As before, the Hawking-Page phase transition is present at T = T HP with thermal AdS 4 (RN-AdS 4 ) being the preferred phase at T < T HP (T > T HP ). For Q > Q crit there are also hairy black hole solutions. We choose not to show them in the right plot of Fig. 3 because we do not expect our perturbative results (87) to be valid for such large charge. To understand the reason we go back to the left panel of Fig. 2. For Q > Q crit , RN-AdS 4 has horizon radius R + 0.2. But our perturbative results for the hairy black holes that merge with the RN-AdS 4 solution are valid for R + 1 so we should not expect our perturbative hairy results to be accurate for this case, at least around the merger point.
Given that the canonical ensemble keeps the temperature T and electric charge Q fixed, it is also relevant to summarize the conclusions above for the thermal phases of the theory in a phase diagram T vs Q. This effectively describes a projection of the 3-dimensional plot F (T, Q). This is done in the right panel of Fig. 2, and the reader can find its discussion in the caption of the figure.
So far we have discussed the preferred global thermodynamic phases of the canonical phase diagram. Given that in this ensemble we keep an intensive thermodynamic quantity − the temperature − fixed we can also discuss the local thermodynamic stability of the solutions. Local thermodynamic stability in the canonical ensemble requires that the specific heat at constant charge is non-negative. To get the second relation we used the first law, dF = −SdT + µdQ. We find that small hairy black holes have C Q < 0. Thus, they are locally thermodynamically unstable. This is illustrated in Fig. 3 (for Q = 0.05 and e = 2.5): the sign of C Q = −T ∂ 2 T F Q is negative 27 . Essentially, a thermal fluctuation that increases the temperature of the small hairy black hole leads to a decrease of its entropy (horizon size).

D. Phase diagram in the grand-canonical ensemble
In the grand-canonical ensemble the system is kept at fixed temperature T and fixed chemical potential µ. The dominant thermal phase is the one that minimizes the Gibbs free energy, G = M − T S − µQ.
The discussion of the hairy solutions in this ensemble requires that we revisit first the properties of the RN-AdS 4 black hole [66,67]. Taking the RN-AdS 4 temperature (28) and solving it with respect to R + we find 25 Much like in the microcanonical ensemble (see footnote 21) the first law of thermodynamics (with dQ = 0) requires that the slope dF/dT is the same for the two branches at the bifurcation point (since they have the same S). This is definitely the case in our plot. However, if we start departing from the regime of validity of our perturbation analysis we increasingly find that the merger is not perfect and the slopes of the two branches no longer match. This is the best criteria to identify the regime of validity of (87) for the canonical ensemble. 26 We can identify the Euclidean time circle of the soliton with any period so that it can have an arbitrary temperature. However, the free energy of this 'thermal' soliton coincides with its mass, hence it is always positive and the soliton is always dominated by thermal AdS 4 . 27 The local thermodynamic stability of the RN-AdS 4 branches can also be inferred from the concavity/convexity of the free energy. This stability was already discussed in detail in [66] and [67] and we do not repeat it here. Recent discussions of local thermodynamic instability can be found in [70][71][72] and references therein. which indicates that for µ < √ 2 there are two branches of RN-AdS 4 solutions: the 'small' and 'large' black holes described by R + − and R + + , respectively. For µ > √ 2 only the large branch is present. This is best illustrated in the plots and caption of Figs. 4 where we display T (R + ) for several fixed µ's (left panel) and µ(R + ) for several fixed T 's (right panel). For later use, notice that small horizon radius corresponds to large temperatures.
In the grand-canonical ensemble, the appropriate phase diagram to describe the thermal phases of the theory is G(µ, T ). To make the presentation clear, in Fig. 5 we plot the G of the RN-AdS 4 black hole as a function of µ for several fixed temperatures. At extremality (left panel), the RN-AdS 4 black hole dominates over thermal AdS 4 . Above extremality, large black holes always have lower G than small RN-AdS 4 black holes whenever they co-exist (G of small black holes is always non-negative). For 0 < T < T c2 = 1 π ∼ 0.32, large black holes have G > 0 for µ < µ HP but G < 0 for µ > µ HP (see middle and right panels). So thermal-AdS 4 is preferred over RN-AdS 4 for µ < µ HP , and a Hawking-Page phase transition occurs at µ = µ HP , with large black holes becoming the dominant phase. However, for T > T c2 = 1 π (not shown in this figure; case represented in the left panel of Fig. 6), there is no Hawking-Page phase transition since large RN-AdS 4 black holes always have G ≤ 0 and dominate the ensemble.
We can now add the hairy thermal phases to the grand-canonical phase diagram. Our perturbative analysis describes hairy black holes that bifurcate, at the onset of superradiance, from small R + RN-AdS 4 black holes. It follows from Figs. 4 that our perturbative computation is valid for large temperatures T T c = √ 3 2π and chemical potential µ < √ 2 ∼ 1.41. Therefore, as an illustrative example, in the left panel of Fig. 6 we fix the temperature to be T = 5 > T c2 > T c and the scalar field charge to be e = 2.5, and we show all the thermal phases (small and large RN-AdS 4 , thermal AdS 4 , thermal soliton and hairy black hole) in the phase diagram G vs µ. In the main plot we find the large RN-AdS 4 black hole branch which exists for any µ and − since we are at T > T c2 − always has G < 0. We also display the small RN-AdS 4 that has G ≥ 0 and exists for µ ≤ √ 2. The magenta point in this small RN-AdS 4 branch signals the zero-mode of superradiance, with small RN-AdS 4 to the right of this point being unstable. The inset plot then zooms the phase diagram to capture the region where the small RN-AdS 4 black holes are unstable. It also shows the hairy black hole branch (red curve) that bifurcates from the small RN-AdS 4 at the onset of superradiance. These black holes have lower G, and are thus preferred, than small RN-AdS 4 black holes, whenever the two phases coexist. For sufficiently large µ, the hairy black holes also have lower G than thermal AdS 4 (which, recall, has G = 0). However, hairy black holes always have higher G than large RN-AdS 4 black holes, when The magenta point signals the onset of superradiance. The inset plot, is a zoom around the superradiant onset. It now also includes the hairy black hole (red curve) that branches-off from the small RN-AdS4 black hole at the onset of superradiance, as well as the thermal soliton (black curve) that starts at µ = 0, G = 0 and extends to negative values of G. The dashed horizontal line describes the Giggs free energy of thermal AdS4. The dominant thermal phase has the lowest Gibbs free energy. Right Panel: Phase diagram β = 1 T vs µ of the grand-canonical ensemble (shown for e = 2.5). This is discussed in the text.
they co-exist. Thermal solitons (black curve) have G < 0 28 . These qualitative properties extend to all other values of T (and e) where our perturbative analysis is valid. Therefore, we conclude that in the perturbative regime where our hairy results hold and the solutions co-exist, large RN-AdS 4 black holes are always the preferred thermal phase over the hairy solutions in the grand-canonical ensemble.

Near region
Here we present the coefficients of the near region expansion (47). The integration constants have already been determined by the boundary conditions and by the matching conditions, as described in Section III D.
Due to the length of the expression, we only show the large y expansion of φ in 2 (y) needed to determine ω 2 . However, note that the next order contributions to φ in 2 (y) and φ out 2 (R) are required to fix a remaining integration constant.