Hairy black holes in AdS with Robin boundary conditions

We study hairy black holes in Einstein-Maxwell-complex scalar theory in four-dimensional asymptotically global anti-de Sitter (AdS) spacetime when the Robin boundary conditions are imposed on the scalar field. This setup is dual to the double trace deformation of strongly interacting field theory on $R \times S^2$ by charged scalar operators. We identify the instability of the Reissner-Nordstr\"{o}m-AdS (RNAdS) black holes under the Robin boundary conditions and construct backreacted geometries branching at the onset of the instability. Also considering associated horizonless geometries called boson stars, we obtain phase diagrams with fairly rich structure in the grand canonical ensemble depending on the boundary condition parameter or the deformation parameter, where phase transition occurs between thermal AdS, RNAdS, charged boson stars, and hairy black holes.


Introduction
Asymptotically anti-de Sitter (AdS) spacetime offers diverse gravitational dynamics. In contrast to asymptotically flat spacetime, black hole geometry can be considered in the canonical ensemble, where asymptotically global AdS experiences the first order phase transition between horizonless and black hole spacetimes [1,2]. Through the AdS/CFT duality [3][4][5], it is interpreted as the confinement/deconfinement phase transition in strongly coupled Yang-Mills theory. When the gravitational theory has U (1) gauge field and charged scalar field, the spontaneous breaking of the gauge symmetry is discussed as the appearance of the superfluid/superconducting phase [6][7][8].
Aforementioned phenomena are often considered with the Dirichlet boundary conditions imposed on the asymptotic behavior of the scalar field at the AdS boundary. However, general conditions known as the Robin boundary conditions (also called mixed boundary conditions) are allowed [9][10][11][12][13] if the field in AdS has a mass close to the Breitenlohner-Freedman bound [14,15]. When the parameter for the Robin boundary conditions exceeds a critical value and the deviation from the Dirichlet boundary condition becomes sufficiently large, the AdS spacetime becomes unstable [11]. The Robin (or mixed) boundary conditions are related to multitrace deformation in the dual field theory in the AdS/CFT interpretation [16][17][18]. Not only for scalar field considered in these literature, but also the Robin boundary conditions can be imposed for vector field and discussed in the context of introducing dynamical gauge field on the AdS boundary [11,19,20]. Robin boundary conditions have also been considered for metric field so as to promote the boundary metric dynamical [21].
In [22], two of the authors studied the linear mode stability of the four-dimensional Reissner-Nordström AdS (RNAdS) spacetime with global AdS asymptotics for neutral and charged complex scalar field perturbations with Robin boundary conditions. 1 The neutral field shows an instability for the Robin boundary conditions with parameters greater than a critical value. The charged scalar field suffers another type of instability due to the electromagnetic interaction with the black hole, which is known as superradiance [25][26][27][28][29]. 2 With the imposition of the Robin boundary conditions, superradiance and the boundary contribution interplay with each other, potentially enhancing the instability caused by the superradiance depending on the parameters of the scalar field and the background spacetime. It was argued in [22] that the instability would change the RNAdS to charged hairy black hole solutions with a nontrivial scalar field satisfying the Robin boundary conditions, which are a candidate for the final fate of the instability. First studied for neutral scalar, the presence of hairy solutions with the Robin boundary conditions has been known; see [31][32][33] for early works. Motivated by [22], we study charged hairy solutions in four dimensional global AdS spacetime in detail.
In this paper, we study hairy black holes that branch at the onset of instability of the charged scalar field with the Robin boundary conditions on the four-dimensional RNAdS, and obtain results that agree with the expectation of [22] explained above. Following [7,8], hairy black holes have been widely studied in Einstein-Maxwell-complex scalar theory in asymptotically AdS spacetime, in both Poincaré and global AdS spacetimes and in various dimensions. In studies of this sort, the Dirichlet (and Neumann) boundary conditions are often considered. For example, the phase diagram in asymptotically global AdS 4 in the grand canonical ensemble was explored in [34]. 3 In this paper, we conduct a comprehensive study on the phase structures realized under the Robin boundary conditions in the grand canonical ensemble. Within the four dimensional global AdS spacetime, charged scalar solitons (boson stars) and hairy black holes in setups including the same model as ours have been considered in [41]. 4 Our work may be viewed as a generalization of this work, clarifying the full phase structure of such solutions under the Robin boundary conditions. This paper is organized as follows. In section 2, we prepare the setup for constructing boson stars and hairy black holes with the Robin boundary conditions. In particular, we study the onset of instability of the four dimensional RNAdS spacetime with respect to the 1 There is a recent work on the quasinormal mode spectrum of a scalar field with the Robin boundary conditions in Schwarzschild AdS4 spacetime [23]. See also superradiance in BTZ black holes with the Robin boundary conditions [24]. 2 Instability of RNAdS can be associated with the violation of near horizon AdS2 BF bound, but it is a necessary condition. For charged scalar, superradiance occurs regardless [30], so here we simply describe the cause of this charged instability as superradiance. 3 Hairy black holes have been also considered in global AdS5 [35][36][37][38]. See also [39,40]. 4 See also prior works in three dimensions [42,43]. See also a recent study of boson stars of mixed boundary conditions deformation [44,45] motivated by the analysis on the large charge limit in CFT [46]. charged scalar field perturbations with the Robin boundary conditions. In section 3, we show results of the phase diagram for our setup under the Robin boundary conditions. Section 4 concludes the paper. In appendix A, we summarize holographic renormalization for the Robin boundary conditions. In appendix B, we discuss the first law of thermodynamics.
In appendix C, we comment on entropies in microcanonical ensemble.

Reissner-Nordström AdS black hole
We consider Einstein-Maxwell-complex scalar theory in four-dimensional asymptotically global AdS spacetime. The action is The gauge coupling constant is written by q. We use units in which Λ = −3 so that the AdS radius can be set to unity. The mass of the scalar is related to the conformal dimension of the scalar operator in the dual field theory as m 2 = ∆(∆ − 3). We set m 2 = −2 in this paper. Then, this equation is solved by ∆ = 1, 2. The equations of motion are We study spherically symmetric static solutions in the spherical AdS boundary. The ansatz can be given by The conformal boundary of the AdS is R × S 2 and located at r = ∞. When f (r) = 1 and χ(r) = 0 (as well as A = φ = 0), the empty AdS is obtained. For horizonless geometries, r = 0 is the center of the AdS. The RNAdS black hole is given by where r h denotes the location of the outermost horizon, satisfying f (r h ) = 0, and Q is the charge of the black hole per solid angle. The total charge is given by Q = 4πQ. We choose the gauge as A t (r h ) = 0, and then we obtain µ = Q/r h , where µ is identified as the chemical potential of the gauge field. For the diagonal metric (2.4), the Hawking temperature and the Bekenstein-Hawking entropy are given by For the RNAdS, the temperature is If µ 2 < 2, the temperature has the minimum T H = T 0 , when Black holes with r h > r 0 are called large black holes, while those with r h < r 0 are small. In the grand canonical ensemble, the first order transition known as the Hawking-Page transition occurs between the RNAdS and AdS when [1,2,47] The horizon radius, or temperature, of this transition can be determined by comparing grand potentials between Euclidean RNAdS (A.62) and thermal AdS geometries. The solution with the lower grand potential is identified to be realized physically. The RNAdS is favored over the thermal AdS in r > r HP , and vice versa. Note that r HP > r 0 . In the grand canonical ensemble, the phase in T > T H is the RNAdS black hole phase. The phase in T < T H corresponds to horizonless AdS geometry, which we refer to as the thermal AdS phase. If µ 2 > 2, the temperature (2.9) becomes zero when This is when the RNAdS black hole becomes extremal. For fixed r h , the range of µ is bounded from above as µ ≤ µ ext . Note that both T HP and T 0 become zero at the borderline value µ 2 = 2. Therefore, for µ 2 > 2, the Hawking-Page transition does not appear in the phase diagram, and the zero temperature geometry is the extremal RNAdS.
To solve the equations of motion, it is convenient to use the z-coordinate defined by z ≡ 1/r. In this coordinate, the AdS boundary is located at z = 0. By the coordinate change, the metric (2.4) can be rewritten as The RNAdS black hole solution (2.6) becomes where z h ≡ 1/r h .
Because of the presence of the horizon, the frequency ω is complex in general. The imaginary part of the frequency is negative Im ω < 0 if the perturbation is stable, and positive Im ω > 0 if instability is induced in the RNAdS background. The border Im ω = 0 is the onset of instability. In the gauge we use, A t (z h ) = 0, both the real and imaginary parts of ω become zero simultaneously at the onset of instability, Re ω = Im ω = 0. 5 This means that, to search the onset of instability of φ, it is sufficient to assume the static perturbation φ(z) and find nontrivial normal modes. At the onset of instability ω = 0, (2.15) is reduced to a static perturbation equation, which depends on three parameters (r h , µ, q) for given m. For the onset of instability, we search normal mode solutions to (2.16) when boundary conditions are imposed at z = 0 and z = z h . On the horizon z = z h , we impose regularity (which used to be the ingoing wave boundary condition if ω = 0, away from the onset of instability). We impose Robin boundary conditions at the AdS boundary z = 0. For m 2 = −2, the asymptotic behavior of φ in z → 0 takes the form where φ 1 and φ 2 are integration constants. Because the scalar mass is in the range −9/4 ≤ m 2 ≤ −5/4, both asymptotic behaviors φ ∼ z and φ ∼ z 2 are normalizable [48]. This means that both coefficients φ 1 and φ 2 can be nonzero for normalizable normal modes. The boundary conditions with φ 1 = 0 and φ 2 = 0 are called Dirichlet, and those with φ 1 = 0 and φ 2 = 0 are Neumann. The case with general values of φ 1 = 0 and φ 2 = 0 is 5 Another gauge is often used that the gauge field vanishes asymptotically while it is nonzero on the horizon, At → 0 (z → 0) and At(z h ) = 0. In that gauge, the perturbation φ = e −iωt φ(z) has a nonzero real part Re ω = 0 at the onset of instability Im ω = 0 [22]. However, this frequency-dependence in the real part can be absorbed by the gauge choice. In this paper, we use a gauge where Re ω = Im ω = 0 at the onset of instability.
called the Robin boundary conditions. The Robin boundary conditions can be specified by a parameter ζ defined by We choose the domain of ζ to be periodic in 0 ≤ ζ < π. The points ζ = 0 and ζ = π/2 correspond to the Dirichlet and Neumann boundary conditions, respectively. Under the Robin boundary conditions, we search the onset of instability for the scalar field perturbation in the four-dimensional parameter space (ζ, r h , µ, q). Technically, for a set of three parameters (r h , µ, q), we integrate the perturbation equation (2.16) from the horizon to the AdS boundary and read off the asymptotic coefficients φ 1 and φ 2 in (2.17), from which ζ can be obtained. This procedure gives a location of the onset of instability in the (ζ, r h , µ, q) parameter space. Iterating this procedure while varying the values for the three parameters (r h , µ, q), we obtain a relation among the four parameters (ζ, r h , µ, q). Thus, for instance, fixing (r h , q), we obtain the onset of instability is given as a curve in (µ, ζ) plane.
In the horizonless limit r h = 0, the perturbation equation (2.16) can be solved analytically. The background is the global AdS f = 1 with a constant gauge field A t = µ. The perturbation equation (2.16) then becomes When the horizon is absent, we impose φ (z)| z=∞ = 0 at the center of the AdS. With this boundary condition and m 2 = −2, (2.19) is solved by which is normalized as φ(z)| z=∞ = 1. Expanding this around z = 0, we find [11,22] cot ζ = − µq tan(πµq/2) . (2.21) For r h = 0, µ and q always show up in a pair µq. The set of the parameters (ζ, µ, q) satisfying the above relation gives a normal mode in the global AdS. While the global AdS is stable against linear perturbations, nontrivial scalar solutions branch from the AdS at the normal modes. For this reason, with a slight abuse of terminology, we also refer to the location of the AdS normal modes as the onset of instability.
In figure 1, we show (a) the location of the AdS charged scalar field normal modes (r h = 0) and (b) the onset of instability of the RNAdS for r h = 0.1, 0.5, 1 at q = 1. In figure 1(b), the value of µ is bounded from above by extremality as µ ≤ µ ext (2.12), which is marked by the vertical red dashed line for each r h . In the same figure, the RNAdS is unstable to the charged scalar field perturbation above each curve, which can be found by studying full quasinormal modes by including nonzero frequencies ω (see also [22]). Correspondingly, also in figure 1(a), the scalar field will be nonzero in the region upper from the curve.
In figure 1(a), we emphasize that the normal modes can be characterized by the number of nodes in the radial direction, which increases as the curve reaches ζ = 0. The solution  Combining with the analysis of quasinormal modes [22], we find that, in (a), the Schwarzschild AdS is unstable above the curve, and correspondingly in (b), it is unstable to the right of each of the blue and orange curves.
without a node is called the fundamental mode, and the solution with nodes are called overtones. Because overtones cost more energy than the fundamental mode, later in the paper, we consider only the backreacted solutions as a fully nonlinear extension of the fundamental mode.
In figure 1(b), the data for r h = 0.1 shows that, when the coupling q is small, the onset of instability terminates at the extremality before reaching the Dirichlet boundary conditions (ζ = 0). For the Dirichlet boundary conditions to be unstable, a larger r h is necessary.
In figure 2, the onset of instability in the Schwarzschild AdS limit (µ = 0) is shown. The value in the horizonless limit (r h = 0) is analytically given by In figure 2(a), the curve has the minimum at r h 0.4807(< r 0 ) with ζ min 0.6728π and approaches ζ → π as r h → ∞. There are hence no overtones for the Schwarzschild AdS. In figures 3,4,5, we show the onset of instability of the RNAdS for the fundamental modes with different ζ at q = 1, √ 2, 2, respectively. The same onset results are shown in the (µ, r h ) and (µ, T H ) planes. We do so because we will discuss the phase structure in the (µ, T H ) plane of the phase diagram later in the paper, and it will be instructive to have the location of the instability both in the (µ, r h ) and (µ, T H ) planes. In each figure, we show the locations of the onset of instability for 8 parameter values ζ/π = 0, 0.1, . . . , 0.7. (Among the 8 color lines, the lightest color is ζ = 0 and the darkest ζ/π = 0.7.) The red dashed line denotes the extremal RNAdS, below which no regular RNAdS exist. The arc by the thin black line (r h = r HP and T H = T HP ) is the Hawking-Page transition of the RNAdS (2.11). Inside the arc, the thermal AdS is thermodynamically favored over the RNAdS. In figures (a), the gray dashed line (r h = r 0 ) separates the small and large black holes (2.10), and small black holes are inside the arc. In figures (b), the gray dashed line (T H = T 0 ) denotes the minimal temperature T 0 , which is realized when r h = r 0 . In figures (a), the RNAdS is unstable below each onset curve. In the grand canonical ensemble, we are interested in the onset of instability outside the arc given by r h = r HP or T H = T HP . Instability can be understood in terms of superradiance [22]. 6 With the imposition of the Robin boundary condition, superradiance and the boundary contribution interplay with each other, potentially enhancing instability caused by superradiance depending on the parameters of the scalar field and the background spacetime. It is demonstrated in figures (a) that, for fixed q and r h , the value of ζ at the onset increases as µ is decreased. That is, the parameter range of µ for instability is wider as ζ is increased (see also figure 8 in [22]).
In figure 3, we can see that the extremal RNAdS are stable if ζ is small and µ is not sufficiently large, while in figures 4 and 5, the extremal RNAdS are unstable to all ζ. The critical value for the instability of the Dirichlet boundary condition ζ = 0 is q = √ 2, that is, on the phase diagrams on the (µ, r h ) plane (see Panel (a) of figures 3, 4, 5), the onset curve for ζ = 0 ends on the (red dashed) curve of the extremal black hole solutions when q < √ 2, while it ends on the r h = 0 axis when q > √ 2. 7 On the phase diagrams on the (µ, T H ) plane, the extremal black hole solutions correspond to the µ ≥ √ 2 part of the T H = 0 axis, and only a part of it is covered by the instability region for ζ = 0 when q < √ 2, while it is wholly covered by the instability region when q > √ 2.

Hairy black holes
Knowing the onset of instability for the charged scalar field perturbation of the RNAdS, we will construct backreacted hairy black hole solutions branching at the onset of instability.
With the ansatz (2.13), the equations of motion (2.2) are reduced to coupled ODEs for We need the asymptotic behavior of the field variables in z = 0 and z = z h or z → ∞. In the AdS boundary z = 0, the asymptotic solutions are given by where (f 3 , χ 0 , φ 1 , φ 2 , a 0 , a 1 ) are six integration coefficients not determined in the asymptotic analysis. We read off them from the asymptotic form of numerical solutions. With these asymptotic behavior, the metric (2.13) in z → 0 naively becomes This can be rescaled to χ 0 = 0 by the scaling symmetry (redefinition of t) as we will see shortly.
In the presence of the black hole horizon, the regular asymptotic solutions near the horizon z = z h = 1/r h are given by where (χ h , φ h , A h ) are integration constants, and the higher order coefficients are determined fully in terms of them. Two degrees of freedom are considered to be correlated to physical parameters, while the remaining one can be fixed by the scaling symmetry discussed below. In the absence of the horizon, the solutions (2.32)-(2.35) are replaced with the following series in z → ∞, There are again three integration constants. Our ansatz, (2.4) and (2.5), has the following scaling symmetry: where c is an arbitrary constant. By this scaling, solutions with χ 0 = 0 can be rescaled to those with canonical boundary metric satisfying χ 0 = 0. This means that, in numerical calculations, we can set the normalization of χ to an arbitrary value convenient for us without loss of generality. We fix χ h = 0 when we compute and then rescale numerical results by (2.40) to satisfy χ| z=0 = 0. From numerical results, we construct thermodynamic quantities. Carrying out the holographic renormalization as will be described in appendix A, we obtain the expressions of the thermodynamic quantities in terms of the asymptotic coefficients given in (2.27)-(2.30). For the Robin boundary conditions, the scalar field is dual to the dimension 1 operator O 1 . After rescaling to χ 0 = 0, the expression of the total energy, charge, and scalar expectation value for the Robin boundary conditions are obtained in (A.52) and (A.56) as (the subscript R is removed here) We also have the temperature T H through (2.7) and entropy S BH ≡ 8πG N S BH = 8π 2 r 2 h through (2.8).
We consider the grand canonical ensemble to discuss the phase structure. The grand potential is given by where µ = a 0 . The grand potential Ω can be evaluated in two different expressions. One is by the combination of thermodynamic quantities as in the RHS of (2.42), and the other is directly by a bulk integral (A.60). These give the same physical quantity. In practice, the latter is less convenient and costly because of the necessity of numerically cancelling the divergent terms in the integrand. Hence, we use Ω given by (2.42) when we evaluate the phase structure. Numerical solutions to (2.23)-(2.26) satisfying the Robin boundary conditions can be obtained simply by integrating the equations of motion. Specifying (φ h , A h , q, r h ), we integrate (2.23)-(2.26) from the horizon z = z h (or AdS center z = ∞) to the boundary z = 0 and read off (f 3 , χ 0 , φ 1 , φ 2 , a 0 , a 1 ) in the asymptotic boundary behavior (2.27)-(2.30). After the rescaling to set χ 0 → 0, we calculate the thermodynamic quantities and ζ (2.18). By these quantities, the grand canonical phase diagram is given as a four-dimensional space (µ, T H , ζ, q). When we present our results, we use data slices in the four dimensional parameter space.
To check numerical results, we can evaluate first-law-like relations generalizing the first law of thermodynamics/black hole mechanics to the case with a nontrivial scalar field. The expressions are discussed in appendix B. For our solutions in the presence of the scalar field satisfying the Robin boundary conditions, we can use (B.4), Note that this contains an atypical variation with respect to cot ζ = φ 1 /φ 2 , which is not a thermodynamic quantity but is a parameter in the model. However, if we compare between numerical solutions where both φ 1 and φ 2 vary while their ratio is not fixed, the first-lawlike equation (2.43) is useful. We find that the above relation is satisfied within numerical errors.

Neutral boson stars and black holes
First, we consider neutral solutions. 8 Here, we focus on the phase transition in the canonical ensemble. 9 8 See also [31][32][33] for neutral scalar hair solutions with the Robin boundary conditions (in the presence of nonlinear scalar potential). 9 In a recent work [49], Hawking-Page transition was discussed for neutral black holes with scalar field in the Dirichlet (Neumann) theory when the scalar source was nonzero. Here, this gravitational setup is studied as the double trace deformation with zero scalar source. As we will explain, using the free energy formula for the Robin boundary conditions, we obtain the phase structure comprehensively. Qualitatively similar to [49], we also observe that the Hawking-Page transition temperature increases when the scalar field is nonzero. Before discussing the black holes, let us recall the basic features of the horizonless solutions (see also [50]). In figure 6, we show the energy and expectation value of neutral horizonless solutions branching at the appearance of the zero normal mode of AdS. Because the scalar field is subject to the Robin boundary conditions, we call the horizonless solutions Robin boson stars. These are a one-parameter family of solutions parametrized by ζ. Scalar hair grows in ζ > ζ c 0.6805π, where the phase transition is of second order. The quantities in the figure approach O 1 → +∞ and E → −∞ in ζ → π. In the following, we will consider two kinds of generalization: black holes by introducing temperature T H , and gauge field by adding (µ, q).
Without the gauge field, the phase structure is specified by two parameters (T H , ζ). In this situation, the free energy we compare for determining the phase structure is nothing but the grand potential (2.42) with µ = 0, Ω| µ=0 = (E − T H S BH ) µ=0 . We compare free energies among thermal AdS, Schwarzschild AdS, Robin boson stars, and black holes with neutral scalar hair, which we call Robin black holes. The free energy for the thermal AdS is zero, and that for the Schwarzschild AdS is given by (A.62) with µ = 0.
In figure 7(a), we show an example of the comparison of free energies among neutral solutions. For r h 1, the two solutions experience the first order phase transition. In the figure, the lines of the r h = 0.9 Robin black holes (blue) and boson stars (black dashed) cross around ζ/π ∼ 0.8. For r h 1, the free energy of the Robin black holes is always lower than that of the Robin boson stars. The free energy for r h = 1.1 (orange) is shown in the figure.
The phase diagram for the neutral solutions is summarized in figure 7(b). The vertical green line at ζ = ζ c is the second order phase transition from thermal AdS to Robin boson stars. The blue line in ζ ≥ ζ t at the border of the Schwarzschild AdS and hairy Robin black holes is the second order phase transition for growing scalar hair, where ζ t 0.6847π. (3.1) Because the source of the scalar field is assumed to be zero, the scalar becomes nonzero spotaneously when the temperature is decreased [51,52]. As ζ increases, the critical tem-  perature for this scalar hair formation rises, and in the limit ζ → π (cot ζ → −∞), the Robin black holes dominate at any high temperatures. The red line in ζ ≥ ζ t marks the first order Hawking-Page transition between Robin black holes and Robin boson stars. The short orange segment in ζ c ≤ ζ ≤ ζ t (see the inset) is the first order phase transition between Schwarzschild AdS and Robin boson stars; for ζ in this region, Robin black holes have the higher free energy than these two, and hence the first order phase transition is between the Schwarzschild AdS and Robin boson stars. The three lines (red, orange, blue) merge at ζ = ζ t and T H 0.3184, (3.2) which corresponds to the triple point at which the Schwarzschild AdS black hole, Robin black hole, and the Robin boson star have the same free energy. The temperature (3.2) at the triple point is slightly higher than the transition temperature T HP for the Schwarzschild AdS and thermal AdS phases (2.11), T HP | µ=0 = 1/π 0.3183. We find that the Hawking-Page transition temperature depends on ζ very mildly. We were not able to pin down the line of the Hawking-Page transition up to ζ → π because of numerical limitations. But, as long as we could confirm, the transition temperature (red line) behaves as which is close to T HP | µ=0 0.3183 and is mostly insensitive to ζ. Thus, for the Hawking-Page transition temperature of neutral geometries, the effect of the Robin boundary conditions on the free energy is minor. This behavior suggests that the free energies of the Robin boson star and the Robin black hole changes by almost the same amount when ζ changes.

Charged boson stars
To proceed with the reduced number of parameters, we discuss charged but horizonless solutions with the Robin boundary conditions, which we call charged Robin boson stars.
Features of these solutions have been explored in [41] in the same setup as ours, but here we discuss the solutions in the phase space parametrized by (µ, ζ, q).
In figure 8, the expectation value O 1 is compared for three cases with ζ > ζ c , ζ = ζ c , and ζ < ζ c for q = 1, 2. Recall that AdS at µ = 0 is unstable for ζ ≥ ζ c 0.6805π for forming neutral boson stars. This implies that, for ζ > ζ c , charged Robin boson stars are connected to neutral Robin boson stars (with µ = 0) by turning on finite µ. Meanwhile, for ζ < ζ c , they branch at the appearance of the zero normal mode of AdS with finite µ. For example, when ζ/π = 0.6, the value of µ at the branching point of the condensed solution in figures 8(a) and 8(b) (i.e. the limit of O 1 → 0) corresponds to µ in the r h = 0 limit in figures 3(a) and 5(a), respectively. The boundary between these two families of the solutions is ζ = ζ c . In addition, in figure 8(a), these charged Robin boson stars have the maximal µ above which solutions do not exist. In the inset, the data region near the maximal µ for ζ/π = 0.6 is enlarged. While it might be visually unclear even in the inset, the region near the largest µ has a spiral structure, corresponding to the attractor solutions discussed in [41]. In figure 8(b), the expectation value can be arbitrarily large. This corresponds to solutions allowing the planar limit discussed in [41].
The boundary between these two distinct behaviors depends both on ζ and q. The tendency is that the spiral structure disappears (moves to infinity on the (µ, O 1 ) plane) as q and ζ are increased. Not only O 1 but also the energy E shows a qualitatively similar behavior. This tendency can be qualitatively understood as an outcome of the balance between the gravitational attraction, scalar field pressure and the electric repulsion. When q is small, the electric repulsion is weak and then there is a critical mass (and O 1 ) for a boson star beyond which the boson star cannot exist. When q is large, the electric repulsion becomes strong enough to sustain the boson star against the gravitational collapse, and correspondingly the mass and O 1 can become arbitrarily large.
The grand potential of charged Robin boson stars always satisfy Ω < 0, where thermal AdS has Ω = 0. Therefore, when the charged Robin boson stars exist, they are always preferred over the thermal AdS. This feature is the same as that in the neutral case, in which the boson stars have the smaller free energy than the thermal AdS (see section 3.1 and figure 7(a)).

Charged black holes
Finally, we consider black holes with nontrivial charged scalar field with the Robin boundary conditions. We call these hairy Robin black holes. The phase space depends on the all four parameters (µ, T H , ζ, q).
In figure 9, phase diagrams for q = 1 are shown for different ζ. In each figure, the blue line on the border between the RNAdS and hairy Robin black holes denotes the second-order phase transition below which the scalar hair forms. The red line is the first order Hawking-Page transition between hairy Robin black holes and charged Robin boson stars. The orange segment denotes the first order phase transition between the RNAdS and charged Robin boson stars. The black dashed line is plotted for reference of the Hawking-Page transition between the thermal AdS and RNAdS (2.11), although it is not physically dominant because it is superseded by the charged Robin boson star phase.
Starting from a large value of ζ, we browse notable features in the phase structure by decreasing ζ.
• ζ > ζ t 0.6847π: In figure 9(a) (see Eq. (3.1) for the definition of ζ t ), neutral solutions (µ = 0) can have nontrivial scalar hair. Thermal AdS does not appear because its free energy is always higher than Robin boson stars when the latter exist as solutions. Hence, the phase diagram contains three phases: zero scalar RNAdS, hairy Robin black holes, and charged Robin boson stars. By decreasing the temperature, the RNAdS spontaneously grows the scalar hair, and then the hairy Robin black hole transitions to the charged Robin boson stars. This feature is common to all µ.
• ζ t > ζ > ζ c 0.6805π: In figure 9(b), the phase structure for this parameter region is shown for ζ/π = 0.682. When ζ is decreased to ζ t , the two phase transition lines (blue and red) first meet at µ = 0. As shown in figure 7(b), ζ = ζ t is bigger than ζ = ζ c where the thermal AdS phase shows up. This means that, in ζ < ζ t , the phase transition from the RNAdS to charged Robin boson stars (orange line) appears.
• ζ c > ζ 0.24π: In ζ < ζ c , the thermal AdS phase can be present as µ is increased from 0 until the charged Robin boson stars branch from thermal AdS as discussed in figure 8. The phase diagram in this parameter region is shown in figure 9(c). The vertical green line is the second order phase transition between thermal AdS and charged Robin boson stars.
• In figures 9(a)-9(c), the Hawking-Page transition (red line) will approach T H → 0 as µ is increased. We were not able to compute until this limit due to tough numerics, but we can see that the transition line will go down towards T H → 0 for a wide  parameter range (in ζ/π 0.24). We also expect that the Hawking-Page transition should go to T HP → 0 before boson star solutions disappear at the upper limit in µ for Robin boson stars with small q (discussed in section 3.2).
• ζ 0.24π: When ζ is decreased further, the Hawking-Page transition between hairy Robin black holes and charged Robin boson stars reaches zero temperature and disappears. For q = 1, this occurs in a small parameter window of ζ near ζ/π 0.24. Figure 9(d) is the phase diagram for ζ/π = 0.239. This has four phases, but the charged Robin boson stars and hairy Robin black holes are separated by the RNAdS, and correspondingly there is a small gap of µ where the extremal RNAdS survives in the phase diagram at zero temperature. The hairy Robin black holes branch from the extremal RNAdS.
• ζ 0.24π: The charged Robin boson star phase then disappears when ζ is decreased further. In figure 9(e), the phase diagram at ζ/π = 0.2 is shown. While charged Robin boson stars also exist as solutions in this parameter region, their grand potential is always bigger than that of hairy Robin black holes, and hence they do not show up in the grand canonical phase diagram.
When the coupling q is increased, the ζ dependence of the phase structure can be different.
• For q = √ 2, the phase structures of figures 9(a), 9(b), and 9(c) are observed for ζ > 0, but those of figures 9(d) and 9(e) are absent because no stable extremal RNAdS exist even for the Dirichlet boundary condition ζ = 0. Instead, at ζ = 0, a phase structure not shown here appears (see figure 7(a) in [34]). It contains three phases, where thermal AdS and hairy Robin BH are separated by the RN AdS BH reaching sufficiently low temperature. There, the phase of the charged boson stars also disappears because the onset is exactly on the T H = 0 axis [34].
• For q > √ 2, the scalar hair grows at finite temperatures before extremality is reached, because all the extremal solutions with T H = 0 are unstable toward scalar hair formation when q > √ 2, as explained in section 2.2. Therefore, the phase structures depicted in figures 9(d) and 9(e) are absent in q > √ 2.
There is qualitative difference for the phase structures of the Robin boundary conditions from those for the Dirichlet boundary condition (see [34]) that some of the phase structures in figure 9 are absent in the same system under the Dirichlet boundary condition. The structures of figures 9(a) and 9(b) do not exist for the Dirichlet boundary condition, because thermal AdS phase should appear in small µ region when ζ < ζ c and particularly in the Dirichlet case (ζ = 0). Adding to that, the presence of neutral Robin black hole phase (with µ = 0) for ζ > ζ t observed in figure 9(a) is another unique feature of the Robin boundary condition. The structure of figure 9(c) is observed for the Dirichlet boundary condition with a gauge coupling q > √ 2 (in our normalization) [34]. For the Robin boundary condition, however, this phase structure can be seen even for small q if ζ is sufficiently large. The structure of figure 9(d) is not seen for the Dirichlet boundary condition because the charged boson star phase disappears at the same time as T HP → 0 at q = √ 2 (see [34]). The structure of figure 9(e) is typical in q < √ 2.

Conclusion
We considered charged boson stars and black holes in four-dimensional Einstein-Maxwellcomplex scalar theory with the Robin boundary conditions for the charged scalar field in asymptotically global AdS spacetime. This setup is dual to the double trace deformation of three-dimensional dual field theory on R × S 2 with a dimension 1 charged scalar operator. The current setup has the four-dimensional parameter space (T H , µ, q, ζ), and the consideration of the Robin boundary conditions offers the most general solutions in the four-dimensional Einstein-Maxwell-complex scalar theory. The phase structure and phase transition are studied in the grand canonical ensemble. There are four phases characterized by the presence and absence of the black hole horizon and nontrivial scalar hair. There is an interplay between two kinds of instability on the formation of a charged scalar hair, the one caused by the Robin boundary conditions and the other by the chemical potential or the black hole charge. These introduce the richer phase structure compared with the case of the Dirichlet boundary condition, as explained in section 3. We considered the Robin boundary conditions for scalar field in this paper. This type of boundary conditions can be also imposed on vector and metric fields [11,[19][20][21]. It will be interesting to consider phases of gravitational solutions where the Robin boundary conditions are imposed on such different kinds of field. Rather recently, the Robin boundary conditions are utilized for studies in various context including the holography and also the supergravity (see e.g. [53][54][55]). Our study would provide useful information to clarify various properties in these cases, such as the thermodynamical phase structures and also the dynamical (in)stabilities.

A Holographic renormalization
We carry out holographic renormalization in the asymptotically global AdS spacetime with the Robin boundary conditions (also called the mixed boundary conditions) [56]. We follow the calculations in [57], the application of which to complex scalar theory in global AdS is straightforward.
We use the r-coordinate in calculation. The asymptotic solutions near the AdS boundary (2.27)-(2.30) take the form In the following, we assume that the scaling (2.40) has been applied so that χ 0 = 0. The action is regularized by introducing a cutoff surface at r = r Λ . Let M denote the regularized spacetime manifold defined in r ≤ r Λ and ∂M the cutoff surface at r = r Λ . The bulk action (2.1) accompanied by the Gibbons-Hawking term can be regularized as with K ≡ K ij γ ij being the trace of the extrinsic curvature K ij with respect to the induced metric γ ij on ∂M (i, j run over the three-dimensional coordinates on ∂M ). The extrinsic curvature is given by where n µ is an outward unit normal g µν n µ n ν = 1. However, the "bare" action (A.5) diverges when the cutoff is simply removed by taking the limit of r Λ → ∞. This divergence can be cancelled by counterterms S ct . Including S ct formally, we can define a subtracted action that is finite in the limit r Λ → ∞ as Then, removing the cutoff gives a renormalized action, The form of S ct depends on the boundary conditions at the AdS boundary. We will discuss the cases of the Dirichlet theory, Neumann theory, and double trace deformation in turn.
Dirichlet theory When ζ = 0, our Einstein-Maxwell complex scalar system is treated as the Dirichlet theory that has a dimension 2 operator O 2 in the dual field theory on the AdS boundary. This is also called the standard quantization. The counterterms for this case can be given by [57][58][59] where R γ is the Ricci scalar for γ ij , and we ignored derivative terms of the scalar field that do not contribute in our spherically symmetric static solutions. With these counterterms, let S D ren denote the renormalized action for the Dirichlet theory.
The expectation values of field theory operators can be obtained through variation as where Φ D = √ 2φ 1 is the source of the scalar operator, Ψ i denotes that of the gauge field, and h ij are the metric components of the boundary R × S 2 . For the gauge field, we turn on the chemical potential Ψ t = µ = a 0 .
The boundary stress energy tensor can be practically calculated as follows. From the subtracted action, the stress energy tensor on the cutoff surface can be obtained as where (G γ ) ij = (R γ ) ij − 1 2 R γ γ ij is the Einstein tensor for the induced metric. This scales as (T γ ) ij ∼ 1/r Λ because γ ij ∼ 1/r 2 Λ and √ −γ ∼ r 3 Λ . Hence, by switching from γ ij to h ij , the expectation value of the boundary stress energy tensor (A.11) reads Explicitly, the components are given by 16) where (θ, ψ) denote the coordinates on S 2 introduced as dΩ 2 2 = dθ 2 + sin 2 θdψ 2 . From the stress energy tensor, the total energy (also called the total mass) of the Dirichlet theory is expressed as where 8πG N is included in the definition of the LHS.
The expectation values for the matter fields are These quantities are the densities per solid angle. The total charge is given by Similarly, the scalar expectation value integrated over the sphere is The trace of the stress energy tensor satisfies If both φ 1 and φ 2 are nonzero, the theory that gives the variation (A.11) can be interpreted as Dirichlet theory in the presence of a nonzero source Φ D . The nonzero trace (A. 21) then indicates that the conformal symmetry is explicitly broken by the source. When the source is absent Φ D = 0, i.e. ζ = 0, the expression of the energy (A.17) reduces to We can also calculate the finite Euclidean on-shell action when the counterterms are added. 10 Using the equations of motion, we obtain (see section 3.4 in [8]) where the last term, which is not a total derivative, is due to the spherical topology of the global AdS. By this relation, the bulk action (2.1) is simplified to Combining (A.23) and (A.24), we obtain the finite Lorentzian renormalized on-shell action, The Euclidean on-shell action S E can be obtained by replacing dt → − 1/T H 0 dτ where τ denotes the Euclidean time. It is related to the grand potential for the Dirichlet theory as Ω D ≡ 8πG N T H S E . The expression of the grand potential in terms of the bulk integral is hence given by where we used r Λ = r Λ r h dr + r h to rewrite the cutoff dependence for numerical evaluation of the r-integral.
Neumann theory For ζ = 0, the bulk theory is considered to be dual to the boundary field theory with a dimension 1 scalar operator O 1 . This is known as the alternative quantization. The case of ζ = π/2 is the Neumann theory. It turns out that the source of the scalar operator is identified as Φ N = − √ 2φ 2 , and the expectation value of the scalar operator is 8πG N O 1 = √ 2φ 1 . The renormalized action is modified from the Dirichlet theory as follows.
The Neumann theory is the Legendre transform of the Dirichlet theory [16], where N denotes the Neumann theory, and The variation with respect to the scalar field gives (A.29) The variation of the renormalized Neumann action hence takes the form In the above equation, T ij δh ij contains the contribution from the variation of (A.28) by h ij , which shifts (A.14)-(A.16). The stress energy tensor for the Neumann theory is thus given by The trace of the stress energy tensor is Correspondingly, the total energy is and so is the grand potential, Ω N = Ω D + Ω LT , where E LT = Ω LT = 8πφ 1 φ 2 . When the source is absent Φ N = 0, i.e. ζ = π/2, the energy is given by Double trace deformation For ζ = 0 nor π/2, the theory is interpreted as double trace deformation of the Neumann theory. For this, we need to include additional finite boundary terms in order for consistent variation with respect to the source in the deformed theory. We give the source in the form where α is a real parameter. The undeformed Neumann theory corresponds to α = 0. For this source, we need an additional finite boundary term, This term corresponds to the relevant double trace deformation of the dual field theory. The renormalized action is modified to The renormalized action equipped with the finite term S Dtr gives the correct variation with respect to the source Φ R . The scalar field variation of (A.38) is The full variation of S R ren takes the form The above stress energy tensor T ij contains finite contribution from the variation of S Dtr with respect to h ij , shifting the expressions of the Neumann theory (A.31)-(A.33). The expectation values in (A.40) are given by In our setup, we consider the Robin boundary conditions (2.18) as the double trace deformation with vanishing source Φ R = 0. From (2.18), we choose α = − cot ζ, and the condition for the source is reduced to Under this condition, the components of the stress energy tensor (A.41)-(A.43) become When ζ = π/2, these expressions reduce to those for the Neumann boundary conditions (A.31)-(A.33). The total energy is given by This can be decomposed into individual contributions as where E Dtr = −4πφ 2 1 cot ζ is the contribution of S Dtr . Among these, E LT +E Dtr is interpreted as the energy stored on the AdS boundary. 11 Note that E LT +E Dtr = 0 when ζ = π/2, while it is not when ζ = π/2 (and ζ = 0, of course). The total charge and scalar expectation value are given by Using cot ζ = φ 2 /φ 1 , we can rewrite (A.47)-(A.49) as The total energy is expressed as The trace of the energy momentum tensor can be written in the form This implies the spontaneous breaking of the conformal symmetry in the double trace deformed theory when the scalar operator acquires an expectation value. The grand potential of the double trace deformed theory is also shifted from the Dirichlet and Neumann theories by a finite term as Ω R = Ω N + Ω Dtr = Ω D + Ω LT + Ω Dtr , (A. 58) where Ω Dtr = E Dtr = −4π cot ζ φ 2 1 = −4πφ 1 φ 2 . (A.59) The expression of the grand potential in terms of the bulk integral is shifted from (A. 26) as RNAdS For the RNAdS black holes (Eqs. (2.4)-(2.6)), we have (the label of D, N, R is removed because the scalar field is zero) The grand potential is In thermal AdS, r h = 0, we obtain E = Ω = 0. The Hawking-Page transition between the RNAdS and thermal AdS phases (2.11) occurs when the black hole reaches Ω = 0. The grand potential of the RNAdS (A.62) is Ω > 0 for r h < r HP and Ω < 0 for r h > r HP (2.11  Adding the double trace deformation (A.37), we can rewrite this for the double trace deformed theory. In this step, we can treat also α (defined by Eq. (A.36)) as an independent variable. By doing this, we can compare solutions with different values of α. We obtain 12 On general grounds, this first law in the presence of a nonzero scalar source follows from the fact that the grand potential is the generating function for responses of sources. In [61], this was discussed for the holographic superconductor model same as this paper except in the probe limit with the planar AdS boundary. Recently in [62], this scalar source contribution to the first law was derived by using Wald's formalism [63,64]. See also [65][66][67] for earlier discussions. However, if the last term is taken into account, we can use (B.4) as a relation useful to compare solutions where ζ varies in general. We can use any of the above equations to check numerical results because these are rewriting of the same relation.

C Comparison of entropy in microcanonical ensemble
In the main text, we have seen the phase structures in the grand canonical ensemble. We can also consider the microcanonical ensemble where the total energy (mass) E and charge Q are treated as independent variables. In this ensemble, we can argue the fate of an unstable RNAdS black hole by comparing the entropies between solutions with and without scalar at the same (E, Q) (see also Dirichlet boundary condition [35,36]).
In figure 10, we show the entropies of the two kinds of the solutions in the (Q, S BH ) plane for E = 10, ζ/π = 0.6 and q = 1. The black curve is the entropy of the RNAdS with E = 10. The extremal RNAdS is marked by the red dot, and the onset of instability for the branching of the hairy black holes is shown by the blue dot. When the RNAdS and hairy Robin black hole both exist at the same parameters (E, Q, ζ, q), the latter has the higher entropy than the former. We also examined other values of the parameters (E, Q, ζ, q) and found that hairy black holes have higher entropy than RNAdS when solutions overlap (see also the same comparison in the Dirichlet boundary condition [35,36]). This implies that an unstable RNAdS can dynamically evolve into a hairy black hole in the microcanonical ensemble when it is perturbed and nonlinear time evolution is considered. In figure 10, the zero entropy limit of the hairy Robin black hole is the zero size limit r h → 0 with diverging temperature T H → ∞. The profile of the field variables (f, χ, φ, A t ) approaches that of a charged Robin boson star.