11institutetext: Department of Physics, Rikkyo University, Nishi-Ikebukuro,
Toshima, Tokyo 171-8501, Japan22institutetext: Astronomical Institute,
Graduate School of Science, Tohoku University, Aoba, Sendai 980-8578,
Japan33institutetext: Niels Bohr International Academy, Niels Bohr Institute,
Blegdamsvej 17, 2100 Copenhagen, Denmark44institutetext: Department of
Physics, Chuo University, Kasuga, Bunkyo, Tokyo 112-8551, Japan

# Hairy black holes in AdS with Robin boundary conditions

Tomohiro Harada 1 Takaaki Ishii 2,3 Takuya Katagiri 4 and Norihiro Tanahashi
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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ö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.

††preprint: RUP-23-6

##  1 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 Hawking:1982dh ; Witten:1998zw . Through the AdS/CFT duality
Maldacena:1997re ; Witten:1998qj ; Gubser:1998bc , 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 Gubser:2008px ;
Hartnoll:2008vx ; Hartnoll:2008kx .

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 if the field in AdS has a
mass close to the Breitenlohner-Freedman bound Breitenlohner:1982bm ;
Breitenlohner:1982jf . 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
Ishibashi:2004wx . The Robin (or mixed) boundary conditions are related to
multitrace deformation in the dual field theory in the AdS/CFT interpretation
Klebanov:1999tb ; Witten:2001ua ; Berkooz:2002ug . 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 Ishibashi:2004wx ; Witten:2003ya ;
Marolf:2006nd . Robin boundary conditions have also been considered for metric
field so as to promote the boundary metric dynamical Compere:2008us .

In Katagiri:2020mvm , 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.111There is a recent work on the quasinormal mode
spectrum of a scalar field with the Robin boundary conditions in Schwarzschild
AdS4 spacetime KKMS . See also superradiance in BTZ black holes with the Robin
boundary conditions Dappiaggi:2017pbe . 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 PhysRevD.7.949 ; Hawking:1999dp ; Uchikata:2011zz ;
Green:2015kur ; Dias:2016pma .222Instability 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 Ishii:2022lwc , so here we
simply describe the cause of this charged instability as superradiance. 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 conjectured in
Katagiri:2020mvm  that there exist 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.

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 support the conjecture of
Katagiri:2020mvm  explained above. Following Hartnoll:2008vx ; Hartnoll:2008kx
, 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 AdS4 in the grand canonical ensemble was
explored in Basu:2016mol .333Hairy black holes have been also considered in
global AdS5 Maeda:2010hf ; Basu:2010uz ; Bhattacharyya:2010yg ; Dias:2011tj .
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 global AdS spacetime, charged scalar solitons (boson stars) and
hairy black holes in settings similar to ours have been considered
Gentle:2011kv .444See also a recent study of boson stars of mixed boundary
conditions deformation Liu:2020uaz ; Guo:2020bqz  motivated by the analysis on
the large charge limit in CFT Hellerman:2015nra . 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 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.

##  2 Setup

###  2.1 Reissner-Nordström AdS black hole

We consider Einstein-Maxwell-complex scalar theory in four-dimensional
asymptotically global AdS spacetime. The action is

$S=\frac{1}{8\pi
G_{N}}\int\mathrm{d}^{4}x\sqrt{-g}\left(\frac{1}{2}\left(R-2\Lambda\right)-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-|D\phi|^{2}-m^{2}|\phi|^{2}\right),$
(1)

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, and
$D_{\mu}\phi=\partial_{\mu}\phi-iqA_{\mu}\phi$. The gauge coupling constant is
written by $q$. We use units in which $\Lambda=-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}=\Delta(\Delta-3)$.
We set $m^{2}=-2$ in this paper. Then, this equation is solved by
$\Delta=1,2$. The equations of motion are

$G_{\mu\nu}+\Lambda
g_{\mu\nu}=T_{\mu\nu},\quad\nabla_{\mu}F^{\mu\nu}=J^{\nu},\quad\left(D_{\mu}D^{\mu}-m^{2}\right)\phi=0,$
(2)

where

$\begin{split}T_{\mu\nu}&=F_{\mu\lambda}F_{\nu}{}^{\lambda}+(D_{\mu}\phi)^{*}D_{\nu}\phi+(D_{\nu}\phi)^{*}D_{\mu}\phi+g_{\mu\nu}\mathcal{L},\\\
\mathcal{L}&=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-|D\phi|^{2}-m^{2}|\phi|^{2},\\\
J^{\mu}&=2q^{2}|\phi|^{2}A^{\mu}+iq(\phi^{*}\partial^{\mu}\phi-\phi\partial^{\mu}\phi^{*}).\end{split}$
(3)

We study spherically symmetric static solutions in the spherical AdS boundary.
The ansatz can be given by

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=-\left(1+r^{2}\right)f(r)e^{-\chi(r)}\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{\left(1+r^{2}\right)f(r)}+r^{2}\mathrm{d}\Omega_{2}^{2},$
(4) $\displaystyle A$ $\displaystyle=A_{t}(r)\mathrm{d}t,\qquad\phi=\phi(r).$
(5)

The conformal boundary of the AdS is $R\times S^{2}$ and located at
$r=\infty$. When $f(r)=1$ and $\chi(r)=0$ (as well as $A=\phi=0$), the empty
AdS is obtained. This is called thermal AdS. For horizonless geometries, $r=0$
is the center of the AdS.

The RNAdS black hole is given by

$\begin{split}\left(1+r^{2}\right)f(r)&=1+r^{2}-\left(1+r_{h}^{2}+\frac{Q^{2}}{2r_{h}^{2}}\right)\frac{r_{h}}{r}+\frac{Q^{2}}{2r^{2}},\\\
A_{t}(r)&=\mu-\frac{Q}{r},\quad\chi(r)=\phi(r)=0,\end{split}$ (6)

where $r_{h}$ denotes the location of the outermost horizon, satisfying
$f(r_{h})=0$, and $Q$ is the charge density of the black hole. We choose the
gauge as $A_{t}(r_{h})=0$, and then we obtain $\mu=Q/r_{h}$, where $\mu$ is
identified as the chemical potential of the gauge field. For the diagonal
metric (4), the Hawking temperature and the Bekenstein-Hawking entropy are
given by

$\displaystyle T_{\mathrm{H}}$
$\displaystyle=\frac{1}{4\pi}(1+r_{h}^{2})f^{\prime}(r_{h})e^{-\chi(r_{h})/2},$
(7) $\displaystyle S_{\mathrm{BH}}$
$\displaystyle=\frac{8\pi^{2}r_{h}^{2}}{8\pi G_{N}}.$ (8)

For the RNAdS, the temperature is

$T_{\mathrm{H}}=\frac{2(1+3r_{h}^{2})-\mu^{2}}{8\pi r_{h}}.$ (9)

If $\mu^{2}<2$, the temperature has the minimum $T_{\mathrm{H}}=T_{0}$, when

$r_{h}=\frac{\sqrt{2-\mu^{2}}}{\sqrt{6}}\equiv r_{0},\quad
T_{0}=\frac{\sqrt{3(2-\mu^{2})}}{2\sqrt{2}\,\pi}=\frac{3r_{0}}{2\pi}.$ (10)

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
thermal AdS when Hawking:1982dh ; Witten:1998zw ; Chamblin:1999hg

$r_{h}=\frac{\sqrt{2-\mu^{2}}}{\sqrt{2}}\equiv r_{\mathrm{HP}},\quad
T_{\mathrm{HP}}=\frac{\sqrt{2-\mu^{2}}}{\sqrt{2}\,\pi}=\frac{r_{\mathrm{HP}}}{\pi}.$
(11)

The horizon radius, or temperature, of this transition can be determined by
comparing grand potentials between the RNAdS and thermal AdS (108). The
solution with the lower grand potential is identified to be realized
physically. The RNAdS is favored over the thermal AdS in $r>r_{\mathrm{HP}}$,
and vice versa. Note that $r_{\mathrm{HP}}>r_{0}$. If $\mu^{2}>2$, the
temperature (9) becomes zero when

$\mu=\sqrt{2(1+3r_{h}^{2})}\equiv\mu_{\mathrm{ext}}.$ (12)

This is when the RNAdS black hole becomes extremal. For fixed $r_{h}$, the
range of $\mu$ is bounded from above as $\mu\leq\mu_{\mathrm{ext}}$. Note that
both $T_{\mathrm{HP}}$ and $T_{0}$ become zero at the borderline value
$\mu^{2}=2$. Therefore, for $\mu^{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\equiv 1/r$. In this coordinate, the AdS boundary is located at
$z=0$. By the coordinate change, the metric (4) can be rewritten as

$\mathrm{d}s^{2}=\frac{1}{z^{2}}\left[-\left(1+z^{2}\right)f(z)e^{-\chi(z)}\mathrm{d}t^{2}+\frac{\mathrm{d}z^{2}}{\left(1+z^{2}\right)f(z)}+\mathrm{d}\Omega_{2}^{2}\right].$
(13)

The RNAdS black hole solution (6) becomes

$\left(1+z^{2}\right)f(z)=1+z^{2}-\left(1+z_{h}^{2}+\frac{Q^{2}z_{h}^{4}}{2}\right)\frac{z^{3}}{z_{h}^{3}}+\frac{Q^{2}z^{4}}{2},\quad
A_{t}=\mu\left(1-\frac{z}{z_{h}}\right),$ (14)

where $z_{h}\equiv 1/r_{h}$.

###  2.2 Instability of RNAdS

We consider spherically symmetric scalar field perturbations of the RNAdS,
$\phi(z)e^{-i\omega t}$. The Klein-Gordon equation takes the form

$\phi^{\prime\prime}+\left(\frac{F^{\prime}}{F}-\frac{2}{z}\right)\phi^{\prime}-\left(\frac{m^{2}}{z^{2}F}-\frac{(\omega+qA_{t})^{2}}{F^{2}}\right)\phi=0,$
(15)

where ${}^{\prime}\equiv\partial_{z}$, $F=(1+z^{2})f$, $m^{2}=-2$, and
$f,A_{t}$ are given by the RNAdS background (14). Because of the presence of
the horizon, the frequency $\omega$ is complex in general. The imaginary part
of the frequency is negative $\mathrm{Im}\,\omega<0$ if the perturbation is
stable, and positive $\mathrm{Im}\,\omega>0$ if instability is induced in the
RNAdS background. The border $\mathrm{Im}\,\omega=0$ is the onset of
instability. In the gauge we use, $A_{t}(z_{h})=0$, both the real and
imaginary parts of $\omega$ become zero simultaneously at the onset of
instability, $\mathrm{Re}\,\omega=\mathrm{Im}\,\omega=0$.555Another gauge is
often used that the gauge field vanishes asymptotically while it is nonzero on
the horizon, ${A}_{t}\to 0\ (z\to 0)$ and ${A}_{t}(z_{h})\neq 0$. In that
gauge, the perturbation $\phi=e^{-i\omega t}\phi(z)$ has a nonzero real part
$\mathrm{Re}\,\omega\neq 0$ at the onset of instability
$\mathrm{Im}\,\omega=0$ Katagiri:2020mvm . However, this frequency-dependence
in the real part can be absorbed by the gauge choice. In this paper, we use a
gauge where $\mathrm{Re}\,\omega=\mathrm{Im}\,\omega=0$ at the onset of
instability. This means that, to search the onset of instability of $\phi$, it
is sufficient to assume the static perturbation $\phi(z)$ and find nontrivial
normal modes.

At the onset of instability $\omega=0$, (15) is reduced to a static
perturbation equation,

$\phi^{\prime\prime}+\left(\frac{F^{\prime}}{F}-\frac{2}{z}\right)\phi^{\prime}-\left(\frac{m^{2}}{z^{2}F}-\frac{q^{2}A_{t}^{2}}{F^{2}}\right)\phi=0,$
(16)

which depends on three parameters $(r_{h},\mu,q)$ for given $m$. For the onset
of instability, we search normal mode solutions to (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
$\omega\neq 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 $\phi$ in $z\to 0$ takes the form

$\phi=\phi_{1}z+\phi_{2}z^{2}+\cdots,$ (17)

where $\phi_{1}$ and $\phi_{2}$ are integration constants. Because the scalar
mass is in the range $-9/4\leq m^{2}\leq-5/4$, both asymptotic behaviors
$\phi\sim z$ and $\phi\sim z^{2}$ are normalizable Balasubramanian:1998sn .
This means that both coefficients $\phi_{1}$ and $\phi_{2}$ can be nonzero for
normalizable normal modes. The boundary conditions with $\phi_{1}=0$ and
$\phi_{2}\neq 0$ are called Dirichlet, and those with $\phi_{1}\neq 0$ and
$\phi_{2}=0$ are Neumann. The case with general values of $\phi_{1}\neq 0$ and
$\phi_{2}\neq 0$ is called the Robin boundary conditions. The Robin boundary
conditions can be specified by a parameter $\zeta$ defined by

$\cot\zeta=\frac{\phi_{2}}{\phi_{1}}.$ (18)

We choose the domain of $\zeta$ to be periodic in $0\leq\zeta<\pi$. The points
$\zeta=0$ and $\zeta=\pi/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
$(\zeta,r_{h},\mu,q)$. Technically, for a set of three parameters
$(r_{h},\mu,q)$, we integrate the perturbation equation (16) from the horizon
to the AdS boundary and read off the asymptotic coefficients $\phi_{1}$ and
$\phi_{2}$ in (17), from which $\zeta$ can be obtained. This procedure gives a
location of the onset of instability in the $(\zeta,r_{h},\mu,q)$ parameter
space. Iterating this procedure while varying the values for the three
parameters $(r_{h},\mu,q)$, we obtain a relation among the four parameters
$(\zeta,r_{h},\mu,q)$. Thus, for instance, fixing $(r_{h},q)$, we obtain the
onset of instability is given as a curve in $(\mu,\zeta)$ plane.

In the horizonless limit $r_{h}=0$, the perturbation equation (16) can be
solved analytically. The background is the global AdS $f=1$ with a constant
gauge field $A_{t}=\mu$. The perturbation equation (16) then becomes

$\phi^{\prime\prime}-\frac{2}{z(1+z^{2})}\phi^{\prime}-\left(\frac{m^{2}}{z^{2}(1+z^{2})}-\frac{\mu^{2}q^{2}}{(1+z^{2})^{2}}\right)\phi=0.$
(19)

When the horizon is absent, we impose $\phi^{\prime}(z)|_{z=\infty}=0$ at the
center of the AdS. With this boundary condition and $m^{2}=-2$, (19) is solved
by

$\phi(z)=\frac{z}{\mu q}\sin\left(\mu q\cot^{-1}z\right),$ (20)

which is normalized as $\phi(z)|_{z=\infty}=1$. Expanding this around $z=0$,
we find Katagiri:2020mvm ; Ishibashi:2004wx

$\cot\zeta=-\frac{\mu q}{\tan(\pi\mu q/2)}.$ (21)

For $r_{h}=0$, $\mu$ and $q$ always show up in a pair $\mu q$. The set of the
parameters $(\zeta,\mu,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.

(a)

(b)

Figure 1: (a) The charged scalar field normal modes of global AdS with
constant $A_{t}=\mu$. (b) The onset of instability of the RNAdS for
$r_{h}=0.1,0.5,1$ at $q=1$.

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 1LABEL:sub@fig:onset_q1_zeta, the value
of $\mu$ is bounded from above by extremality as $\mu\leq\mu_{\mathrm{ext}}$
(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 $\omega$ (see also Katagiri:2020mvm ).
Correspondingly, also in figure 1LABEL:sub@fig:onset_rh0, the scalar field
will be nonzero in the region upper from the curve.

In figure 1LABEL:sub@fig:onset_rh0, we emphasize that the normal modes can be
characterized by the number of nodes in the radial direction, which increases
as the curve reaches $\zeta=0$. The solution 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 1LABEL:sub@fig:onset_q1_zeta, 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 ($\zeta=0$). For
the Dirichlet boundary conditions to be unstable, a larger $r_{h}$ is
necessary.

(a)

(b)

Figure 2: The onset of instability of the Schwarzschild AdS ($\mu=0$). The
blue and orange parts are respectively in the large and small black hole
branches of the Schwarzschild AdS. Combining with the analysis of quasinormal
modes Katagiri:2020mvm , 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.

In figure 2, the onset of instability in the Schwarzschild AdS limit $(\mu=0)$
is shown. The value in the horizonless limit $(r_{h}=0)$ is analytically given
by

$\zeta_{c}=\pi-\tan^{-1}(\pi/2)\simeq 0.6805\pi.$ (22)

In figure 2LABEL:sub@fig:onset_sch_rh, the curve has the minimum at
$r_{h}\simeq 0.4807(<r_{0})$ with $\zeta_{\mathrm{min}}\simeq 0.6728\pi$ and
approaches $\zeta\to\pi$ as $r_{h}\to\infty$. There are hence no overtones for
the Schwarzschild AdS.

(a)

(b)

Figure 3: The onset of instability of the RNAdS for $q=1$. Color lines show
the locations of the onset of instability for $\zeta/\pi=0,0.1,\dots,0.7$ from
bottom to top (lighter to darker). For each value of $\zeta$, instability
occurs in the region below the corresponding curve. Panel (a): phase diagram
on the $(\mu,r_{H})$ plane. The red dashed line corresponds to the extremal
solutions, below which no black hole solutions exist. The curves for
$r_{h}=r_{\text{HP}}$ (black solid) and $r_{h}=r_{0}$ (gray dashed) denote the
Hawking-Page transition and the transition between the small/large black
holes. Panel (b): phase diagram on the $(\mu,T_{\mathrm{H}})$ plane. The curve
for $T_{\mathrm{H}}=T_{0}$ (gray dashed) denote the minimum horizon
temperature, which corresponds to $r_{H}=r_{0}$ and no black hole solutions
with $T<T_{\mathrm{H}}$ exist. Note that $r_{h}=0$ on the $(\mu,r_{h})$ plane
is mapped to $T_{\mathrm{H}}\to\infty$ on the $(\mu,T_{\mathrm{H}})$ plane.
The red dashed line at $T_{\mathrm{H}}=0$ corresponds to the extremal
solutions.

(a)

(b)

Figure 4: Same as figure 3 but for $q=\sqrt{2}$.

(a)

(b)

Figure 5: Same as figure 3 but for $q=2$.

In figures 3,4,5, we show the onset of instability of the RNAdS for the
fundamental modes with different $\zeta$ at $q=1,\sqrt{2},2$, respectively.
The same onset results are shown in the $(\mu,r_{h})$ and
$(\mu,T_{\mathrm{H}})$ planes. We do so because we will discuss the phase
structure in the $(\mu,T_{\mathrm{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 $(\mu,r_{h})$ and $(\mu,T_{\mathrm{H}})$ planes. In each figure,
we show the locations of the the onset of instability for 8 parameter values
$\zeta/\pi=0,0.1,\dots,0.7$. (Among the 8 color lines, the lightest color is
$\zeta=0$ and the darkest $\zeta/\pi=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_{\mathrm{HP}}$ and $T_{\mathrm{H}}=T_{\mathrm{HP}}$) is the
Hawking-Page transition of the RNAdS (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 (10), and small
black holes are inside the arc. In figures (b), the gray dashed line
($T_{\mathrm{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_{\mathrm{HP}}$ or
$T_{\mathrm{H}}=T_{\mathrm{HP}}$.

Instability can be understood in terms of superradiance Katagiri:2020mvm
.666For charged scalar fields in non-rotating charged black hole spacetimes,
superradiance is an amplified scattering of the scalar field around the black
hole if its frequency $\tilde{\omega}$ satisfies $\tilde{\omega}<q\mu$
PhysRevD.7.949 ; Hawking:1999dp ; Dias:2016pma . In the context of linear
modal stability analysis, the frequency $\tilde{\omega}$ is replaced by a real
part of the quasinormal mode frequency, ${\rm Re}[\tilde{\omega}_{\rm QNM}]$,
which is a function of $\zeta$ under the Robin boundary condition. The
condition for instability is then ${\rm Re}[\tilde{\omega}_{\rm QNM}]<q\mu$;
therefore, instability caused by superradiance is controlled by $\zeta$ (see
the details in Katagiri:2020mvm ).  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
$\zeta$ at the onset increases as $\mu$ is decreased. That is, the parameter
range of $\mu$ for instability is wider as $\zeta$ is increased (see also
figure 8 in Katagiri:2020mvm ).

In figure 3, we can see that the extremal RNAdS are stable if $\zeta$ is small
and $\mu$ is not sufficiently large, while in figures 4 and 5, the extremal
RNAdS are unstable to all $\zeta$. The critical value for the instability of
the Dirichlet boundary condition $\zeta=0$ is $q=\sqrt{2}$, that is, on the
phase diagrams on the $(\mu,r_{h})$ plane (see Panel (a) of figures 3, 4, 5),
the onset curve for $\zeta=0$ ends on the (red dashed) curve of the extremal
black hole solutions when $q<\sqrt{2}$, while it ends on the $r_{h}=0$ axis
when $q>\sqrt{2}$.777See Basu:2016mol  for an earlier discussion on the phase
diagram for the Dirichlet boundary condition. Note that the Hawking-Page
transition occurs at $\mu=1$ in this reference due to the normalization
different from ours. See also Dias:2016pma , which studied massless scalar but
observations given there can be easily generalized to massive cases. On the
phase diagrams on the $(\mu,T_{\mathrm{H}})$ plane, the extremal black hole
solutions correspond to the $\mu\geq\sqrt{2}$ part of the $T_{\mathrm{H}}=0$
axis, and only a part of it is covered by the instability region for $\zeta=0$
when $q<\sqrt{2}$, while it is wholly covered by the instability region when
$q>\sqrt{2}$.

###  2.3 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 (13), the equations of motion (2)
are reduced to coupled ODEs for $f(z),\chi(z),\phi(z),A_{t}(z)$ as

$\displaystyle F^{\prime}-\left(\frac{3}{z}+z\phi^{\prime
2}\right)F-e^{\chi}\left(\frac{z^{3}A_{t}^{\prime
2}}{2}+\frac{zq^{2}A_{t}^{2}\phi^{2}}{F}\right)+z+\frac{3}{z}+\frac{2\phi^{2}}{z}$
$\displaystyle=0,$ (23) $\displaystyle\chi^{\prime}-2z\phi^{\prime
2}-\frac{2ze^{\chi}q^{2}A_{t}^{2}\phi^{2}}{F^{2}}$ $\displaystyle=0,$ (24)
$\displaystyle\phi^{\prime\prime}+\left(\frac{F^{\prime}}{F}-\frac{\chi^{\prime}}{2}-\frac{2}{z}\right)+\left(\frac{2}{z^{2}F}+\frac{e^{\chi}q^{2}A_{t}^{2}}{z^{2}F}\right)\phi$
$\displaystyle=0,$ (25) $\displaystyle
A_{t}^{\prime\prime}+\frac{\chi^{\prime}}{2}A_{t}^{\prime}-\frac{2q^{2}\phi^{2}}{z^{2}F}A_{t}$
$\displaystyle=0,$ (26)

where $F=(1+z^{2})f$.

We need the asymptotic behavior of the field variables in $z=0$ and $z=z_{h}$
or $z\to\infty$. In the AdS boundary $z=0$, the asymptotic solutions are given
by

$\displaystyle f(z)$
$\displaystyle=1+\phi_{1}^{2}z^{2}+f_{3}z^{3}+\left(2\phi_{1}^{4}+2\phi_{2}^{2}+\frac{a_{1}^{2}e^{\chi_{0}}}{2}\right)z^{4}+\cdots,$
(27) $\displaystyle\chi(z)$
$\displaystyle=\chi_{0}+\phi_{1}^{2}z^{2}+\frac{8}{3}\phi_{1}\phi_{2}z^{3}+\left(\frac{3}{2}\phi_{1}^{4}+2\phi_{2}^{2}-q^{2}a_{0}^{2}\phi_{1}^{2}e^{\chi_{0}}\right)z^{4}+\cdots,$
(28) $\displaystyle\phi(z)$
$\displaystyle=\phi_{1}z+\phi_{2}z^{2}+\frac{1}{2}\phi_{1}\left(\phi_{1}^{2}-q^{2}a_{0}^{2}e^{\chi_{0}}\right)z^{3}+\cdots,$
(29) $\displaystyle A_{t}(z)$
$\displaystyle=a_{0}+a_{1}z+q^{2}a_{0}\phi_{1}^{2}z^{2}+\frac{1}{6}\phi_{1}\left(4q^{2}a_{0}\phi_{2}+(2q^{2}-1)a_{1}\phi_{1}\right)z^{3}+\cdots,$
(30)

where $(f_{3},\chi_{0},\phi_{1},\phi_{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 (13) in $z\to 0$ naively becomes

$\mathrm{d}s^{2}|_{z\to
0}=\frac{1}{z^{2}}\left(-e^{-\chi_{0}}\mathrm{d}t^{2}+\mathrm{d}z^{2}+\mathrm{d}\Omega_{2}^{2}\right).$
(31)

This can be rescaled to $\chi_{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

$\displaystyle f(z)$
$\displaystyle=-\frac{3+z_{h}^{2}+2\phi_{h}^{2}-e^{\chi_{h}}A_{h}^{2}z_{h}^{3}}{z_{h}(1+z_{h}^{2})}(z-z_{h})+\cdots,$
(32) $\displaystyle\chi(z)$ $\displaystyle=\chi_{h}+\cdots,$ (33)
$\displaystyle\phi(z)$ $\displaystyle=\phi_{h}+\cdots,$ (34) $\displaystyle
A_{t}(z)$ $\displaystyle=A_{h}(z-z_{h})+\cdots,$ (35)

where $(\chi_{h},\phi_{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 (32)–(35) are replaced with the
following series in $z\to\infty$,

$\displaystyle f(z)$ $\displaystyle=1+O(z^{-2}),$ (36) $\displaystyle\chi(z)$
$\displaystyle=\chi_{h}+O(z^{-2}),$ (37) $\displaystyle\phi(z)$
$\displaystyle=\phi_{h}+O(z^{-2}),$ (38) $\displaystyle A_{t}(z)$
$\displaystyle=A_{h}+O(z^{-2}).$ (39)

There are again three integration constants.

Our ansatz, (4) and (5), has the following scaling symmetry:

$t\to e^{-c/2}t,\quad\chi\to\chi-c,\quad A_{t}\to e^{c/2}A_{t},$ (40)

where $c$ is an arbitrary constant. By this scaling, solutions with
$\chi_{0}\neq 0$ can be rescaled to those with canonical boundary metric
satisfying $\chi_{0}=0$. This means that, in numerical calculations, we can
set the normalization of $\chi$ to an arbitrary value convenient for us
without loss of generality. We fix $\chi_{h}=0$ when we compute and then
rescale numerical results by (40) to satisfy $\chi|_{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 (27)–(30). For the Robin boundary conditions, the scalar
field is dual to the dimension 1 operator $\mathcal{O}_{1}$. After rescaling
to $\chi_{0}=0$, the expression of the total energy, charge, and scalar
expectation value for the Robin boundary conditions are obtained in (98) and
(102) as (the subscript $R$ is removed here)

$\begin{split}\mathcal{E}&=4\pi(-f_{3}+3\phi_{1}^{2}\cot\zeta)=4\pi(-f_{3}+3\phi_{1}\phi_{2}),\\\
\mathcal{Q}&=-4\pi
a_{1},\qquad\langle\mathcal{O}_{1}\rangle=4\pi\sqrt{2}\,\phi_{1}.\end{split}$
(41)

We also have the temperature $T_{\mathrm{H}}$ through (7) and entropy
$\mathcal{S}_{\mathrm{BH}}\equiv 8\pi G_{N}S_{\mathrm{BH}}=8\pi^{2}r_{h}^{2}$
through (8).

We consider the grand canonical ensemble to discuss the phase structure. The
grand potential is given by

$\Omega=\mathcal{E}-T_{\mathrm{H}}\mathcal{S}_{\mathrm{BH}}-\mu\mathcal{Q},$
(42)

where $\mu=a_{0}$. The grand potential $\Omega$ can be evaluated in two
different expressions. One is by the combination of thermodynamic quantities
as in the RHS of (42), and the other is directly by a bulk integral (106).
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 $\Omega$ given by (42) when we
evaluate the phase structure.

Numerical solutions to (23)–(26) satisfying the Robin boundary conditions can
be obtained simply by integrating the equations of motion. Specifying
$(\phi_{h},A_{h},q,r_{h})$, we integrate (23)–(26) from the horizon $z=z_{h}$
(or AdS center $z=\infty)$ to the boundary $z=0$ and read off
$(f_{3},\chi_{0},\phi_{1},\phi_{2},a_{0},a_{1})$ in the asymptotic boundary
behavior (27)–(30). After the rescaling to set $\chi_{0}\to 0$, we calculate
the thermodynamic quantities and $\zeta$ (18). To check numerical results, we
can evaluate the first law of thermodynamics/black hole mechanics for these
thermodynamic quantities. The expressions of the first law are given in
appendix B. In particular, the formula for the Robin boundary conditions is
obtained in (112) as

$\mathrm{d}\mathcal{E}=T_{\mathrm{H}}\mathrm{d}\mathcal{S}_{\mathrm{BH}}+\mu\mathrm{d}\mathcal{Q}+\frac{1}{8\pi}\langle\mathcal{O}_{1}\rangle^{2}\mathrm{d}(\cot\zeta).$
(43)

We find that this first law is satisfied within numerical errors. The grand
canonical phase diagram is four dimensional $(\mu,T_{\mathrm{H}},\zeta,q)$.
When we present our results, we use data slices in the four dimensional
parameter space.

##  3 Results

###  3.1 Neutral boson stars and black holes

First, we consider neutral solutions.888See also Hertog:2004dr ; Hertog:2004ns
for neutral scalar hair solutions with the Robin boundary conditions (in the
presence of nonlinear scalar potential). Here, we focus on the phase
transition in the canonical ensemble.999In a recent work Basu:2021bpr ,
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 Basu:2021bpr , we also observe that the Hawking-Page
transition temperature increases when the scalar field is nonzero.

(a)

(b)

Figure 6: Neutral Robin boson stars.

Before discussing the black holes, let us recall the basic features of the
horizonless solutions (see also Bizon:2020yqs ). 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
$\zeta$. Scalar hair grows in $\zeta>\zeta_{c}\simeq 0.6805\pi$, where the
phase transition is of second order. The quantities in the figure approach
$\langle\mathcal{O}_{1}\rangle\to+\infty$ and $\mathcal{E}\to-\infty$ in
$\zeta\to\pi$. In the following, we will consider two kinds of generalization:
black holes by introducing temperature $T_{\mathrm{H}}$, and gauge field by
adding $(\mu,q)$.

Without the gauge field, the phase structure is specified by two parameters
$(T_{\mathrm{H}},\zeta)$. In this situation, the free energy we compare for
determining the phase structure is nothing but the grand potential (42) with
$\mu=0$,
$\Omega|_{\mu=0}=\left(\mathcal{E}-T_{\mathrm{H}}\mathcal{S}_{\mathrm{BH}}\right)_{\mu=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 (108) with $\mu=0$.

(a)

(b)

Figure 7: (a) Comparison of free energies between Robin black holes with
$r_{h}=0.9$ (top, blue) and 1.1 (bottom, orange) and boson stars (black
dashed). (b) The phase diagram for neutral solutions.

In figure 7LABEL:sub@fig:neutral_bh_F, we show an example of the comparison of
free energies among neutral solutions. For $r_{h}\lesssim 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 $\zeta/\pi\sim 0.8$. For $r_{h}\gtrsim 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
7LABEL:sub@fig:neutral_bh_phase. The vertical green line at $\zeta=\zeta_{c}$
is the second order phase transition from thermal AdS to Robin boson stars.
The blue line in $\zeta\geq\zeta_{t}$ at the border of the Schwarzschild AdS
and hairy Robin black holes is the second order phase transition for growing
scalar hair, where

$\zeta_{t}\simeq 0.6847\pi.$ (44)

Because the source of the scalar field is assumed to be zero, the scalar
becomes nonzero spotaneously when the temperature is decreased Faulkner:2010fh
; Faulkner:2010gj . As $\zeta$ increases, the critical temperature for this
scalar hair formation rises, and in the limit $\zeta\to\pi$
($\cot\zeta\to-\infty$), the Robin black holes dominate at any high
temperatures. The red line in $\zeta\geq\zeta_{t}$ marks the first order
Hawking-Page transition between Robin black holes and Robin boson stars. The
short orange segment in $\zeta_{c}\leq\zeta\leq\zeta_{t}$ (see the inset) is
the first order phase transition between Schwarzschild AdS and Robin boson
stars; for $\zeta$ 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 $\zeta=\zeta_{t}$ and

$T_{\mathrm{H}}\simeq 0.3184,$ (45)

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 (45) at the triple point is slightly higher than the
transition temperature $T_{\mathrm{HP}}$ for the Schwarzschild AdS and thermal
AdS (11), $T_{\mathrm{HP}}|_{\mu=0}=1/\pi\simeq 0.3183$.

We find that the Hawking-Page transition temperature depends on $\zeta$ very
mildly. We were not able to pin down the line of the Hawking-Page transition
up to $\zeta\to\pi$ because of numerical limitations. But, as long as we could
confirm, the transition temperature (red line) behaves as

$T_{\mathrm{H}}\simeq 0.03(\zeta-\zeta_{t})/\pi+0.3184,$ (46)

which is close to $T_{\mathrm{HP}}|_{\mu=0}\simeq 0.3183$ and is mostly
insensitive to $\zeta$. 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
$\zeta$ changes.

###  3.2 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
Gentle:2011kv  in the same setup as ours, but here we discuss the solutions in
the phase space parametrized by $(\mu,\zeta,q)$.

(a)

(b)

Figure 8: $\langle\mathcal{O}_{1}\rangle$ of charged Robin boson stars for
$\zeta>\zeta_{c}$, $\zeta=\zeta_{c}$, $\zeta<\zeta_{c}$ at $q=1$ and $2$.

In figure 8, the expectation value $\langle\mathcal{O}_{1}\rangle$ is compared
for three cases with $\zeta>\zeta_{c}$, $\zeta=\zeta_{c}$, and
$\zeta<\zeta_{c}$ for $q=1,2$. Recall that AdS at $\mu=0$ is unstable for
$\zeta\geq\zeta_{c}\simeq 0.6805\pi$ for forming neutral boson stars. This
implies that, for $\zeta>\zeta_{c}$, charged Robin boson stars are connected
to neutral Robin boson stars (with $\mu=0$) by turning on finite $\mu$.
Meanwhile, for $\zeta<\zeta_{c}$, they branch at the appearance of the zero
normal mode of AdS with finite $\mu$. For example, when $\zeta/\pi=0.6$, the
value of $\mu$ at the branching point of the condensed solution in figures
8LABEL:sub@fig:hairy_star_O1_q1 and 8LABEL:sub@fig:hairy_star_O1_q2 (i.e. the
limit of $\langle O_{1}\rangle\to 0$) corresponds to $\mu$ in the $r_{h}=0$
limit in figures 3LABEL:sub@fig:onset_q1_rh and 5LABEL:sub@fig:onset_q2_rh,
respectively. The boundary between these two families of the solutions is
$\zeta=\zeta_{c}$. In addition, in figure 8LABEL:sub@fig:hairy_star_O1_q1,
these charged Robin boson stars have the maximal $\mu$ above which solutions
do not exist. In the inset, the data region near the maximal $\mu$ for
$\zeta/\pi=0.6$ is enlarged. While it might be visually unclear even in the
inset, the region near the largest $\mu$ has a spiral structure, corresponding
to the attractor solutions discussed in Gentle:2011kv . In figure
8LABEL:sub@fig:hairy_star_O1_q2, the expectation value can be arbitrarily
large. This corresponds to solutions allowing the planar limit discussed in
Gentle:2011kv .

The boundary between these two distinct behaviors depends both on $\zeta$ and
$q$. The tendency is that the spiral structure disappears (moves to infinity
on the $(\mu,\langle O_{1}\rangle)$ plane) as $q$ and $\zeta$ are increased.
Not only $\langle O_{1}\rangle$ 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 $\langle O_{1}\rangle$) 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 $\langle
O_{1}\rangle$ can become arbitrarily large.

The grand potential of charged Robin boson stars always satisfy $\Omega<0$,
where thermal AdS has $\Omega=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
7LABEL:sub@fig:neutral_bh_F).

###  3.3 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 $(\mu,T_{\mathrm{H}},\zeta,q)$.

(a)

(b)

(c)

(d)

(e)

Figure 9: Phase diagram for $q=1$ and $\zeta/\pi=0.7,0.682,0.6,0.239,0.2$. The
dashed and solid black curves denote the Hawking-Page transition temperature
$T_{\mathrm{HP}}$ (Eq. 10) between the RN AdS black holes and the thermal AdS.
When the grand potential of these two solutions are bigger than that of the
charged Robin boson star, the corresponding part of this curve is irrelevant
and does not represent a physical phase boundary, but we added it with a
dashed line for reader’s convenience.

In figure 9, phase diagrams for $q=1$ are shown for different $\zeta$. 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 (11), although it is not
physically dominant because it is superseded by the charged Robin boson star
phase.

Starting from a large value of $\zeta$, we browse notable features in the
phase structure by decreasing $\zeta$.

  * •

$\zeta>\zeta_{t}\simeq 0.6847\pi$: In figure 9LABEL:sub@fig:hairy_bh_phase7_q1
(see Eq. (44) for the definition of $\zeta_{t}$), neutral solutions ($\mu=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 $\mu$.

  * •

$\zeta_{t}>\zeta>\zeta_{c}\simeq 0.6805\pi$: In figure
9LABEL:sub@fig:hairy_bh_phase682_q1, the phase structure for this parameter
region is shown for $\zeta/\pi=0.682$. When $\zeta$ is decreased to
$\zeta_{t}$, the two phase transition lines (blue and red) first meet at
$\mu=0$. As shown in figure 7LABEL:sub@fig:neutral_bh_phase, $\zeta=\zeta_{t}$
is bigger than $\zeta=\zeta_{c}$ where the thermal AdS phase shows up. This
means that, in $\zeta<\zeta_{t}$, the phase transition from the RNAdS to
charged Robin boson stars (orange line) appears.

  * •

$\zeta_{c}>\zeta\gtrsim 0.24\pi$: In $\zeta<\zeta_{c}$, the thermal AdS phase
can be present as $\mu$ 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 9LABEL:sub@fig:hairy_bh_phase6_q1.
The vertical green line is the second order phase transition between thermal
AdS and charged Robin boson stars.

  * •

In figures
9LABEL:sub@fig:hairy_bh_phase7_q1–9LABEL:sub@fig:hairy_bh_phase6_q1, the
Hawking-Page transition (red line) will approach $T_{\mathrm{H}}\to 0$ as
$\mu$ 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_{\mathrm{H}}\to 0$ for a wide parameter range (in $\zeta/\pi\gtrsim 0.24$).
We also expect that the Hawking-Page transition should go to
$T_{\mathrm{HP}}\to 0$ before boson star solutions disappear at the upper
limit in $\mu$ for Robin boson stars with small $q$ (discussed in section
3.2).

  * •

$\zeta\simeq 0.24\pi$: When $\zeta$ 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 $\zeta$ near $\zeta/\pi\simeq 0.24$. Figure
9LABEL:sub@fig:hairy_bh_phase239_q1 is the phase diagram for
$\zeta/\pi=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 $\mu$ where the extremal RNAdS survives in the phase diagram
at zero temperature. The charged Robin boson stars branch from the extremal
RNAdS.

  * •

$\zeta\lesssim 0.24\pi$: The charged Robin boson star phase then disappears
when $\zeta$ is decreased further. In figure
9LABEL:sub@fig:hairy_bh_phase2_q1, the phase diagram at $\zeta/\pi=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 $\zeta$ dependence of the phase
structure can be different.

  * •

For $q=\sqrt{2}$, the phase structures of figures
9LABEL:sub@fig:hairy_bh_phase7_q1, 9LABEL:sub@fig:hairy_bh_phase682_q1, and
9LABEL:sub@fig:hairy_bh_phase6_q1 are observed for $\zeta>0$, but those of
figures 9LABEL:sub@fig:hairy_bh_phase239_q1 and
9LABEL:sub@fig:hairy_bh_phase2_q1 are absent because no stable extremal RNAdS
exist even for the Dirichlet boundary condition $\zeta=0$. Instead, at
$\zeta=0$, a phase structure not shown here appears (see figure 7(a) in
Basu:2016mol ). 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_{\mathrm{H}}=0$ axis Basu:2016mol .

  * •

For $q>\sqrt{2}$, the scalar hair grows at finite temperatures before
extremality is reached, because all the extremal solutions with
$T_{\mathrm{H}}=0$ are unstable toward scalar hair formation when
$q>\sqrt{2}$, as explained in section 2.2. Therefore, the phase structures
depicted in figures 9LABEL:sub@fig:hairy_bh_phase239_q1 and
9LABEL:sub@fig:hairy_bh_phase2_q1 are absent in $q>\sqrt{2}$.

There is qualitative difference for the phase structures of the Robin boundary
conditions from those for the Dirichlet boundary condition (see Basu:2016mol )
that some of the phase structures in figure 9 are absent in the same system
under the Dirichlet boundary condition. The structures of figures
9LABEL:sub@fig:hairy_bh_phase7_q1 and 9LABEL:sub@fig:hairy_bh_phase682_q1 do
not exist for the Dirichlet boundary condition, because thermal AdS phase
should appear in small $\mu$ region when $\zeta<\zeta_{c}$ and particularly in
the Dirichlet case ($\zeta=0$). Adding to that, the presence of neutral Robin
black hole phase (with $\mu=0$) for $\zeta>\zeta_{t}$ observed in figure
9LABEL:sub@fig:hairy_bh_phase7_q1 is another unique feature of the Robin
boundary condition. The structure of figure 9LABEL:sub@fig:hairy_bh_phase6_q1
is observed for the Dirichlet boundary condition with a gauge coupling
$q>\sqrt{2}$ (in our normalization) Basu:2016mol . For the Robin boundary
condition, however, this phase structure can be seen even for small $q$ if
$\zeta$ is sufficiently large. The structure of figure
9LABEL:sub@fig:hairy_bh_phase239_q1 is not seen for the Dirichlet boundary
condition because the charged boson star phase disappears at the same time as
$T_{\mathrm{HP}}\to 0$ at $q=\sqrt{2}$ (see Basu:2016mol ). The structure of
figure 9LABEL:sub@fig:hairy_bh_phase2_q1 is typical in $q<\sqrt{2}$.

##  4 Conclusion

We considered charged boson stars and black holes in four-dimensional
Einstein-Maxwell-complex 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\times S^{2}$ with a dimension 1 charged scalar operator. The current setup
has the four-dimensional parameter space $(T_{\mathrm{H}},\mu,q,\zeta)$, 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 Ishibashi:2004wx ; Witten:2003ya ; Marolf:2006nd ; Compere:2008us . 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.
Ahn:2022azl ; Jeong:2023las ; Nozawa:2022upa ). 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.

###### Acknowledgements.

The authors thank Li Li, Keiju Murata and Matthew Roberts for useful
discussions, and also thank the JGRG webinar series where this work was
initiated. The work of T.H. was supported in part by JSPS KAKENHI Grant
Numbers 19H01895, 20H05853, and 19K0387617. The work of T.I. was supported in
part by JSPS KAKENHI Grant Number 19K03871. The work of T.K. was supported in
part by JSPS KAKENHI Grant Number 17H06360. T.K. thanks for support by VILLUM
FONDEN (grant no. 37766), by the Danish Research Foundation, and under the
European Union’s H2020 ERC Advanced Grant “Black holes: gravitational engines
of discovery” grant agreement no. Gravitas–101052587. The work of N.T. was
supported in part by JSPS KAKENHI Grant Numbers 18K03623 and 21H05189.

##  Appendix 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) Papadimitriou:2007sj . We follow the calculations in
Caldarelli:2016nni , 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 (27)–(30) take the form

$\displaystyle f(r)$
$\displaystyle=1+\frac{\phi_{1}^{2}}{r^{2}}+\frac{f_{3}}{r^{3}}+\left(2\phi_{1}^{4}+2\phi_{2}^{2}+\frac{a_{1}^{2}e^{\chi_{0}}}{2}\right)\frac{1}{r^{4}}+\cdots,$
(47) $\displaystyle\chi(r)$
$\displaystyle=\chi_{0}+\frac{\phi_{1}^{2}}{r^{2}}+\frac{8}{3}\frac{\phi_{1}\phi_{2}}{r^{3}}+\left(\frac{3}{2}\phi_{1}^{4}+2\phi_{2}^{2}-q^{2}a_{0}^{2}\phi_{1}^{2}e^{\chi_{0}}\right)\frac{1}{r^{4}}+\cdots,$
(48) $\displaystyle\phi(r)$
$\displaystyle=\frac{\phi_{1}}{r}+\frac{\phi_{2}}{r^{2}}+\frac{\phi_{1}}{2}\left(\phi_{1}^{2}-q^{2}a_{0}^{2}e^{\chi_{0}}\right)\frac{1}{r^{3}}+\cdots,$
(49) $\displaystyle A_{t}(r)$
$\displaystyle=a_{0}+\frac{a_{1}}{r}+\frac{q^{2}a_{0}\phi_{1}^{2}}{r^{2}}+\frac{1}{6}\phi_{1}\left(4q^{2}a_{0}\phi_{2}+(2q^{2}-1)a_{1}\phi_{1}\right)\frac{1}{r^{3}}+\cdots.$
(50)

In the following, we assume that the scaling (40) has been applied so that
$\chi_{0}=0$.

The action is regularized by introducing a cutoff surface at $r=r_{\Lambda}$.
Let $M$ denote the regularized spacetime manifold defined in $r\leq
r_{\Lambda}$ and $\partial M$ the cutoff surface at $r=r_{\Lambda}$. The bulk
action (1) accompanied by the Gibbons-Hawking term can be regularized as

$S_{\mathrm{reg}}=\frac{1}{8\pi
G_{N}}\int_{M}\mathrm{d}^{4}x\sqrt{-g}\left(\frac{1}{2}\left(R-2\Lambda\right)-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-|D\phi|^{2}-m^{2}|\phi|^{2}\right)+S_{\mathrm{GH}},$
(51)

where

$S_{\mathrm{GH}}=\frac{1}{8\pi G_{N}}\int_{\partial
M}\mathrm{d}^{3}x\sqrt{-\gamma}\,K$ (52)

with $K\equiv K_{ij}\gamma^{ij}$ being the trace of the extrinsic curvature
$K_{ij}$ with respect to the induced metric $\gamma_{ij}$ on $\partial M$
($i,j$ run over the three-dimensional coordinates on $\partial M$). The
extrinsic curvature is given by

$K_{ij}=\frac{1}{2}\delta_{i}^{\mu}\delta_{j}^{\nu}\left(\nabla_{\mu}n_{\nu}+\nabla_{\nu}n_{\mu}\right),$
(53)

where $n^{\mu}$ is an outward unit normal $g_{\mu\nu}n^{\mu}n^{\nu}=1$.
However, the “bare” action (51) diverges when the cutoff is simply removed by
taking the limit of $r_{\Lambda}\to\infty$.

This divergence can be cancelled by counterterms $S_{\mathrm{ct}}$. Including
$S_{\mathrm{ct}}$ formally, we can define a subtracted action that is finite
in the limit $r_{\Lambda}\to\infty$ as

$S_{\mathrm{sub}}=S_{\mathrm{reg}}+S_{\mathrm{ct}}.$ (54)

Then, removing the cutoff gives a renormalized action,

$S_{\mathrm{ren}}=\lim_{r_{\Lambda}\to\infty}S_{\mathrm{sub}}.$ (55)

The form of $S_{\mathrm{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 $\zeta=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 Balasubramanian:1999re ;
deHaro:2000vlm ; Caldarelli:2016nni

$S_{\mathrm{ct}}=-\frac{1}{8\pi G_{N}}\int_{\partial
M}\mathrm{d}^{3}x\sqrt{-\gamma}\left[2+\frac{R_{\gamma}}{2}+\phi^{2}\right],$
(56)

where $R_{\gamma}$ is the Ricci scalar for $\gamma_{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_{\mathrm{ren}}^{D}$ denote the renormalized action for the Dirichlet
theory.

The expectation values of field theory operators can be obtained through
variation as

$\delta
S_{\mathrm{ren}}^{D}=\int\mathrm{d}^{3}x\sqrt{-h}\left(\frac{1}{2}\langle
T^{ij}\rangle\delta h_{ij}+\langle J^{i}\rangle\delta\Psi_{i}+\langle
O_{2}\rangle\delta\Phi_{D}\right),$ (57)

where $\Phi_{D}=\sqrt{2}\phi_{1}$ is the source of the scalar operator,
$\Psi_{i}$ denotes that of the gauge field, and $h_{ij}$ are the metric
components of the boundary $R\times S^{2}$. For the gauge field, we turn on
the chemical potential $\Psi_{t}=\mu=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

$(T_{\gamma})_{ij}=-\frac{2}{\sqrt{-\gamma}}\frac{\delta
S_{\mathrm{sub}}^{D}}{\delta\gamma^{ij}}=\frac{1}{8\pi
G_{N}}\left(-K_{ij}+K\gamma_{ij}-2\gamma_{ij}+(G_{\gamma})_{ij}-\phi^{2}\gamma_{ij}\right),$
(58)

where $(G_{\gamma})_{ij}=(R_{\gamma})_{ij}-\frac{1}{2}R_{\gamma}\gamma_{ij}$
is the Einstein tensor for the induced metric. This scales as
$(T_{\gamma})_{ij}\sim 1/r_{\Lambda}$ because $\gamma^{ij}\sim
1/r_{\Lambda}^{2}$ and $\sqrt{-\gamma}\sim r_{\Lambda}^{3}$. Hence, by
switching from $\gamma_{ij}$ to $h_{ij}$, the expectation value of the
boundary stress energy tensor (57) reads

$\langle
T_{ij}\rangle=\lim_{r_{\Lambda}\to\infty}r_{\Lambda}(T_{\gamma})_{ij}.$ (59)

Explicitly, the components are given by

$\displaystyle 8\pi G_{N}\langle T_{tt}\rangle$
$\displaystyle=-f_{3}+2\phi_{1}\phi_{2},$ (60) $\displaystyle 8\pi
G_{N}\langle T_{\theta\theta}\rangle$
$\displaystyle=-\frac{f_{3}}{2}+2\phi_{1}\phi_{2},$ (61) $\displaystyle 8\pi
G_{N}\langle T_{\psi\psi}\rangle$
$\displaystyle=\sin^{2}\theta\left(-\frac{f_{3}}{2}+2\phi_{1}\phi_{2}\right),$
(62)

where $(\theta,\psi)$ denote the coordinates on $S^{2}$ introduced as
$\mathrm{d}\Omega_{2}^{2}=\mathrm{d}\theta^{2}+\sin^{2}\\!\theta\mathrm{d}\psi^{2}$.
From the stress energy tensor, the total energy (also called the total mass)
of the Dirichlet theory is expressed as

$\mathcal{E}_{D}\equiv 8\pi G_{N}\int\mathrm{d}\Omega_{2}\langle
T_{tt}\rangle=4\pi(-f_{3}+2\phi_{1}\phi_{2}),$ (63)

where $8\pi G_{N}$ is included in the definition of the LHS.

The expectation values for the matter fields are

$8\pi G_{N}\langle J^{t}\rangle=-a_{1},\quad 8\pi G_{N}\langle
O_{2}\rangle=\sqrt{2}\phi_{2}.$ (64)

These quantities are the densities per solid angle. The total charge is given
by

$\mathcal{Q}\equiv 4\pi\cdot 8\pi G_{N}\langle J^{t}\rangle=-4\pi a_{1}.$ (65)

Similarly, the scalar expectation value integrated over the sphere is

$\langle\mathcal{O}_{2}\rangle\equiv 4\pi\cdot 8\pi G_{N}\langle
O_{2}\rangle=4\pi\sqrt{2}\,\phi_{2}.$ (66)

The trace of the stress energy tensor satisfies

$\langle T^{i}{}_{i}\rangle=\frac{2\phi_{1}\phi_{2}}{8\pi
G_{N}}=\Phi_{D}\langle O_{2}\rangle.$ (67)

If both $\phi_{1}$ and $\phi_{2}$ are nonzero, the theory that gives the
variation (57) can be interpreted as Dirichlet theory in the presence of a
nonzero source $\Phi_{D}$. The nonzero trace (67) then indicates that the
conformal symmetry is explicitly broken by the source. When the source is
absent $\Phi_{D}=0$, i.e. $\zeta=0$, the expression of the energy (63) reduces
to $\mathcal{E}_{D}|_{\phi_{1}=0}=-4\pi f_{3}$.

We can also calculate the finite Euclidean on-shell action when the
counterterms are added.101010We thank Li Li for discussions on this
calculation. Using the equations of motion, we obtain (see section 3.4 in
Hartnoll:2008kx )

$\frac{1}{2}\left(R-2\Lambda\right)+\mathcal{L}=-\frac{1}{2}\left(G^{t}{}_{t}+G^{r}{}_{r}\right)=-\frac{1}{\sqrt{-g}}\left(\frac{(1+r^{2})f\sqrt{-g}}{r}\right)^{\prime}+\frac{1}{r^{2}},$
(68)

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 (1) is
simplified to

$\displaystyle S_{\mathrm{bulk}}$ $\displaystyle=\frac{4\pi}{8\pi
G_{N}}\int\mathrm{d}t\mathrm{d}r\left[\left(-r(1+r^{2})fe^{-\frac{\chi}{2}}\right)^{\prime}+e^{-\frac{\chi}{2}}\right]$
$\displaystyle=\frac{4\pi}{8\pi
G_{N}}\int\mathrm{d}t\left[-r_{\Lambda}^{3}-\left(1+\frac{\phi_{1}^{2}}{2}\right)r_{\Lambda}-f_{3}+\frac{4\phi_{1}\phi_{2}}{3}+O\left(\frac{1}{r_{\Lambda}}\right)+\int^{r_{\Lambda}}_{r_{h}}\mathrm{d}r\,e^{-\frac{\chi}{2}}\right],$
(69)

where (47)–(50) were used. This diverges for $r_{\Lambda}\to\infty$. The
divergence can be cancelled by adding the counterterms as well as the Gibbons-
Hawking term,

$S_{\mathrm{GH}}+S_{\mathrm{ct}}=\frac{4\pi}{8\pi
G_{N}}\int\mathrm{d}t\left[r_{\Lambda}^{3}+\frac{\phi_{1}^{2}}{2}r_{\Lambda}+\frac{f_{3}}{2}+\frac{2\phi_{1}\phi_{2}}{3}+O\left(\frac{1}{r_{\Lambda}}\right)\right].$
(70)

Combining (69) and (70), we obtain the finite Lorentzian renormalized on-shell
action,

$S_{L}=S_{\mathrm{ren}}^{D}=\lim_{r_{\Lambda}\to\infty}\left(S_{\mathrm{bulk}}+S_{\mathrm{GH}}+S_{\mathrm{ct}}\right).$
(71)

The Euclidean on-shell action $S_{E}$ can be obtained by replacing
$\int\mathrm{d}t\to-\int_{0}^{1/T_{\mathrm{H}}}d\tau$ where $\tau$ denotes the
Euclidean time. It is related to the grand potential for the Dirichlet theory
as $\Omega_{D}\equiv 8\pi G_{N}T_{\mathrm{H}}S_{E}$. The expression of the
grand potential in terms of the bulk integral is hence given by

$\Omega_{D}=4\pi\int_{r_{h}}^{\infty}\mathrm{d}r(1-e^{-\frac{\chi}{2}})+4\pi\left(\frac{f_{3}}{2}-2\phi_{1}\phi_{2}+r_{h}\right),$
(72)

where we used $r_{\Lambda}=\int_{r_{h}}^{r_{\Lambda}}\mathrm{d}r+r_{h}$ to
rewrite the cutoff dependence for numerical evaluation of the $r$-integral.

#### Neumann theory

For $\zeta\neq 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 $\zeta=\pi/2$ is the Neumann theory. It
turns out that the source of the scalar operator is identified as
$\Phi_{N}=-\sqrt{2}\phi_{2}$, and the expectation value of the scalar operator
is $8\pi G_{N}\langle O_{1}\rangle=\sqrt{2}\phi_{1}$. The renormalized action
is modified from the Dirichlet theory as follows.

The Neumann theory is the Legendre transform of the Dirichlet theory
Klebanov:1999tb ,

$S_{\mathrm{ren}}^{N}=S_{\mathrm{ren}}^{D}+S_{\mathrm{LT}},$ (73)

where $N$ denotes the Neumann theory, and

$S_{\mathrm{LT}}=-\frac{2}{8\pi
G_{N}}\int\mathrm{d}^{3}x\sqrt{-h}\,\phi_{1}\phi_{2}.$ (74)

The variation with respect to the scalar field gives

$\delta_{\phi}S_{\mathrm{LT}}=-\frac{2}{{8\pi
G_{N}}}\int\mathrm{d}^{3}x\sqrt{-h}\left(\phi_{2}\delta\phi_{1}+\phi_{1}\delta\phi_{2}\right)=\int\mathrm{d}^{3}x\sqrt{-h}\left(-\langle
O_{2}\rangle\delta\Phi_{D}+\langle O_{1}\rangle\delta\Phi_{N}\right).$ (75)

The variation of the renormalized Neumann action hence takes the form

$\delta
S_{\mathrm{ren}}^{N}=\int\mathrm{d}^{3}x\sqrt{-h}\left(\frac{1}{2}\langle
T^{ij}\rangle\delta h_{ij}+\langle J^{i}\rangle\delta\Psi_{i}+\langle
O_{1}\rangle\delta\Phi_{N}\right).$ (76)

In the above equation, $\langle T^{ij}\rangle\delta h_{ij}$ contains the
contribution from the variation of (74) by $h_{ij}$, which shifts (60)–(62).
The stress energy tensor for the Neumann theory is thus given by

$\displaystyle 8\pi G_{N}\langle T_{tt}\rangle$
$\displaystyle=-f_{3}+4\phi_{1}\phi_{2},$ (77) $\displaystyle 8\pi
G_{N}\langle T_{\theta\theta}\rangle$ $\displaystyle=-\frac{f_{3}}{2},$ (78)
$\displaystyle 8\pi G_{N}\langle T_{\psi\psi}\rangle$
$\displaystyle=-\sin^{2}\theta\,\frac{f_{3}}{2}.$ (79)

The trace of the stress energy tensor is

$\langle T^{i}{}_{i}\rangle=-\frac{4\phi_{1}\phi_{2}}{8\pi
G_{N}}=2\Phi_{N}\langle O_{1}\rangle.$ (80)

Correspondingly, the total energy is

$\mathcal{E}_{N}=4\pi(-f_{3}+4\phi_{1}\phi_{2})=\mathcal{E}_{D}+\mathcal{E}_{\mathrm{LT}},$
(81)

and so is the grand potential, $\Omega_{N}=\Omega_{D}+\Omega_{\mathrm{LT}}$,
where $\mathcal{E}_{\mathrm{LT}}=\Omega_{\mathrm{LT}}=8\pi\phi_{1}\phi_{2}$.
When the source is absent $\Phi_{N}=0$, i.e. $\zeta=\pi/2$, the energy is
given by $\mathcal{E}_{N}|_{\phi_{2}=0}=-4f_{3}$.

#### Double trace deformation

For $\zeta\neq 0$ nor $\pi/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

$\Phi_{R}=-\sqrt{2}\left(\phi_{2}+\alpha\phi_{1}\right),$ (82)

where $\alpha$ is a real parameter. The undeformed Neumann theory corresponds
to $\alpha=0$. For this source, we need an additional finite boundary term,

$S_{\mathrm{Dtr}}=-\frac{\alpha}{8\pi
G_{N}}\int\mathrm{d}^{3}x\sqrt{-h}\,\phi_{1}^{2}.$ (83)

This term corresponds to the relevant double trace deformation of the dual
field theory. The renormalized action is modified to

$S_{\mathrm{ren}}^{R}=S_{\mathrm{ren}}^{N}+S_{\mathrm{Dtr}}=S_{\mathrm{ren}}^{D}+S_{\mathrm{LT}}+S_{\mathrm{Dtr}}.$
(84)

The renormalized action equipped with the finite term $S_{\mathrm{Dtr}}$ gives
the correct variation with respect to the source $\Phi_{R}$. The scalar field
variation of (84) is

$\delta_{\phi}S_{\mathrm{ren}}^{R}=\frac{1}{8\pi
G_{N}}\int\mathrm{d}^{3}x\sqrt{-h}\left(-2\phi_{1}\delta\phi_{2}-2\alpha\phi_{1}\delta\phi_{1}\right)=\int\mathrm{d}^{3}x\sqrt{-h}\langle
O_{1}\rangle\delta\Phi_{R}.$ (85)

The full variation of $S_{\mathrm{ren}}^{R}$ takes the form

$\delta
S_{\mathrm{ren}}^{R}=\int\mathrm{d}^{3}x\sqrt{-h}\left(\frac{1}{2}\langle
T^{ij}\rangle\delta h_{ij}+\langle J^{i}\rangle\delta\Psi_{i}+\langle
O_{1}\rangle\delta\Phi_{R}\right).$ (86)

The above stress energy tensor $\langle T^{ij}\rangle$ contains finite
contribution from the variation of $S_{\mathrm{Dtr}}$ with respect to
$h_{ij}$, shifting the expressions of the Neumann theory (77)–(79). The
expectation values in (86) are given by

$\displaystyle 8\pi G_{N}\langle T_{tt}\rangle$
$\displaystyle=-f_{3}+4\phi_{1}\phi_{2}+\alpha\phi_{1}^{2},$ (87)
$\displaystyle 8\pi G_{N}\langle T_{\theta\theta}\rangle$
$\displaystyle=-\frac{f_{3}}{2}-\alpha\phi_{1}^{2},$ (88) $\displaystyle 8\pi
G_{N}\langle T_{\psi\psi}\rangle$
$\displaystyle=\sin^{2}\theta\left(-\frac{f_{3}}{2}-\alpha\phi_{1}^{2}\right),$
(89) $\displaystyle 8\pi G_{N}\langle J^{t}\rangle$ $\displaystyle=-a_{1},$
(90) $\displaystyle 8\pi G_{N}\langle O_{1}\rangle$
$\displaystyle=\sqrt{2}\phi_{1}.$ (91)

In our setup, we consider the Robin boundary conditions (18) as the double
trace deformation with vanishing source $\Phi_{R}=0$. From (18), we choose
$\alpha=-\cot\zeta$, and the condition for the source is reduced to

$\Phi_{R}=-\sqrt{2}\left(\phi_{2}-\phi_{1}\cot\zeta\right)=0.$ (92)

Under this condition, the components of the stress energy tensor (87)–(89)
become

$\displaystyle 8\pi G_{N}\langle T_{tt}\rangle$
$\displaystyle=-f_{3}+4\phi_{1}\phi_{2}-\phi_{1}^{2}\cot\zeta,$ (93)
$\displaystyle 8\pi G_{N}\langle T_{\theta\theta}\rangle$
$\displaystyle=-\frac{f_{3}}{2}+\phi_{1}^{2}\cot\zeta,$ (94) $\displaystyle
8\pi G_{N}\langle T_{\psi\psi}\rangle$
$\displaystyle=\sin^{2}\theta\left(-\frac{f_{3}}{2}+\phi_{1}^{2}\cot\zeta\right).$
(95)

When $\zeta=\pi/2$, these expressions reduce to those for the Neumann boundary
conditions (77)–(79). The total energy is given by

$\mathcal{E}_{R}\equiv 8\pi G_{N}\int\mathrm{d}\Omega_{2}\langle
T_{tt}\rangle=4\pi(-f_{3}+4\phi_{1}\phi_{2}-\phi_{1}^{2}\cot\zeta).$ (96)

This can be decomposed into individual contributions as

$\mathcal{E}_{R}=\mathcal{E}_{N}+\mathcal{E}_{\mathrm{Dtr}}=\mathcal{E}_{D}+\mathcal{E}_{\mathrm{LT}}+\mathcal{E}_{\mathrm{Dtr}},$
(97)

where $\mathcal{E}_{\mathrm{Dtr}}=-4\pi\phi_{1}^{2}\cot\zeta$ is the
contribution of $S_{\mathrm{Dtr}}$. Among these,
$\mathcal{E}_{\mathrm{LT}}+\mathcal{E}_{\mathrm{Dtr}}$ is interpreted as the
energy stored on the AdS boundary.111111See also allwright2016robin  for the
relationship between the Robin boundary condition and the modification of the
potential energy near a boundary, which may suggest that the parameter $\zeta$
for the Robin boundary condition controls the amount of the energy stored in
the near-boundary region. Note that
$\mathcal{E}_{\mathrm{LT}}+\mathcal{E}_{\mathrm{Dtr}}=0$ when $\zeta=\pi/2$,
while it is not when $\zeta\neq\pi/2$ (and $\zeta\neq 0$, of course). The
total charge and scalar expectation value are given by

$\mathcal{Q}\equiv-4\pi a_{1},\quad\langle\mathcal{O}_{1}\rangle\equiv
4\pi\sqrt{2}\,\phi_{1}.$ (98)

Using $\cot\zeta=\phi_{2}/\phi_{1}$, we can rewrite (93)–(95) as

$\displaystyle 8\pi G_{N}\langle T_{tt}\rangle$
$\displaystyle=-f_{3}+3\phi_{1}\phi_{2},$ (99) $\displaystyle 8\pi
G_{N}\langle T_{\theta\theta}\rangle$
$\displaystyle=-\frac{f_{3}}{2}+\phi_{1}\phi_{2},$ (100) $\displaystyle 8\pi
G_{N}\langle T_{\psi\psi}\rangle$
$\displaystyle=\sin^{2}\theta\left(-\frac{f_{3}}{2}+\phi_{1}\phi_{2}\right).$
(101)

The total energy is expressed as

$\mathcal{E}_{R}=4\pi(-f_{3}+3\phi_{1}\phi_{2}).$ (102)

The trace of the energy momentum tensor can be written in the form

$8\pi G_{N}\langle
T^{i}{}_{i}\rangle=-\phi_{1}\phi_{2}=-\cot\zeta\,\phi_{1}^{2}=-\frac{\cot\zeta}{2}(8\pi
G_{N})^{2}\langle O_{1}\rangle^{2}.$ (103)

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

$\Omega_{R}=\Omega_{N}+\Omega_{\mathrm{Dtr}}=\Omega_{D}+\Omega_{\mathrm{LT}}+\Omega_{\mathrm{Dtr}},$
(104)

where

$\Omega_{\mathrm{Dtr}}=\mathcal{E}_{\mathrm{Dtr}}=-4\pi\cot\zeta\,\phi_{1}^{2}=-4\pi\phi_{1}\phi_{2}.$
(105)

The expression of the grand potential in terms of the bulk integral is shifted
from (72) as

$\Omega_{R}=4\pi\int_{r_{h}}^{\infty}\mathrm{d}r(1-e^{-\frac{\chi}{2}})+4\pi\left(\frac{f_{3}}{2}-\phi_{1}\phi_{2}+r_{h}\right).$
(106)

#### RNAdS

For the RNAdS black holes (Eqs. (4)–(6)), we have (the label of $D,N,R$ is
removed because the scalar field is zero)

$\mathcal{E}=4\pi r_{h}\left(1+r_{h}^{2}+\frac{\mu^{2}}{2}\right).$ (107)

The grand potential is

$\Omega=2\pi r_{h}\left(1-r_{h}^{2}-\frac{\mu^{2}}{2}\right).$ (108)

In thermal AdS, $r_{h}=0$, we obtain $\mathcal{E}=\Omega=0$. The Hawking-Page
transition between the RNAdS and thermal AdS (11) occurs when the black hole
reaches $\Omega=0$. The grand potential of the RNAdS (108) is $\Omega>0$ for
$r_{h}<r_{\mathrm{HP}}$ and $\Omega<0$ for $r_{h}>r_{\mathrm{HP}}$ (11). For
the RNAdS, we can analytically check that (42) is satisfied, where
$\mathcal{Q}=4\pi\mu r_{h}$ for the RNAdS.

##  Appendix B First law of thermodynamics

To check numerical results, we wish to evaluate the first law of
thermodynamics/black hole mechanics. If we regard solutions with nonzero
$\phi_{1}$ and $\phi_{2}$ as the Dirichlet theory with explicit scalar source,
the first law is given by121212On 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
Ishii:2018ucz , this was discussed for the holographic superconductor model
same as this paper except in the probe limit with the planar AdS boundary.
Recently in Li:2020spf , this scalar source contribution to the first law was
derived by using Wald’s formalism Wald:1993nt ; Iyer:1994ys . See also
Lu:2013ura ; Liu:2013gja ; Lu:2014maa  for ealier discussions.

$\mathrm{d}\mathcal{E}_{D}=T_{\mathrm{H}}\mathrm{d}\mathcal{S}_{\mathrm{BH}}+\mu\mathrm{d}\mathcal{Q}-\langle\mathcal{O}_{2}\rangle\mathrm{d}\Phi_{D}.$
(109)

By the Legendre transform (81), this becomes the first law for the Neumann
theory as

$\mathrm{d}\mathcal{E}_{N}=T_{\mathrm{H}}\mathrm{d}\mathcal{S}_{\mathrm{BH}}+\mu\mathrm{d}\mathcal{Q}-\langle\mathcal{O}_{1}\rangle\mathrm{d}\Phi_{N}.$
(110)

Adding the double trace deformation (83), we can rewrite this as the first law
for the double trace deformed theory as

$\mathrm{d}\mathcal{E}_{R}=T_{\mathrm{H}}\mathrm{d}\mathcal{S}_{\mathrm{BH}}+\mu\mathrm{d}\mathcal{Q}-\langle\mathcal{O}_{1}\rangle\mathrm{d}\Phi_{R}-\frac{1}{8\pi}\langle\mathcal{O}_{1}\rangle^{2}\mathrm{d}\alpha,$
(111)

where we treat $\alpha$ as a variable parameter and therefore
$\mathrm{d}\Phi_{R}=-\sqrt{2}(\mathrm{d}\phi_{2}+\alpha\mathrm{d}\phi_{1}+\phi_{1}\mathrm{d}\alpha)$.
The coefficient of the last term $1/(8\pi)=1/(2\cdot 4\pi)$ is due to the
normalization of $\langle\mathcal{O}_{1}\rangle$ (98). When we impose the
Robin boundary conditions, i.e. $\alpha=-\cot\zeta$ and $\Phi_{R}=0$ (92),
this equation reduces to

$\mathrm{d}\mathcal{E}_{R}=T_{\mathrm{H}}\mathrm{d}\mathcal{S}_{\mathrm{BH}}+\mu\mathrm{d}\mathcal{Q}+\frac{1}{8\pi}\langle\mathcal{O}_{1}\rangle^{2}\mathrm{d}(\cot\zeta).$
(112)

If we consider the family of solutions for fixed $\zeta$, the last term drops,
but we keep this term so that we can compare solutions where $\zeta$ varies in
general. We can use any of the above equations to check numerical results
because these are rewriting of the same relation.

##  Appendix C Comparison of entropy in microcanonical ensemble

Figure 10: Comparison of entropies in the microcanonical ensemble for
$\mathcal{E}=10$, $\zeta/\pi=0.6$, and $q=1$. The end point of of the blue
line is at $(\mathcal{S}_{\mathrm{BH}},\mathcal{Q})=(0,9.955)$, and it
corresponds to a charged Robin boson star.

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) $\mathcal{E}$ and charge $\mathcal{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 $(\mathcal{E},\mathcal{Q})$ (see also Dirichlet boundary
condition Maeda:2010hf ; Basu:2010uz ).

In figure 10, we show the entropies of the two kinds of the solutions in the
$(\mathcal{Q},\mathcal{S}_{\mathrm{BH}})$ plane for $\mathcal{E}=10$,
$\zeta/\pi=0.6$ and $q=1$. The black curve is the entropy of the RNAdS with
$\mathcal{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 $(\mathcal{E},\mathcal{Q},\zeta,q)$, the latter has the higher
entropy than the former. We also examined other values of the parameters
$(\mathcal{E},\mathcal{Q},\zeta,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 Maeda:2010hf ; Basu:2010uz ). 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}\to 0$ with diverging temperature
$T_{\mathrm{H}}\to\infty$. The profile of the field variables
$(f,\chi,\phi,A_{t})$ approaches that of a charged Robin boson star.

## References

  * (1) S.W. Hawking and D.N. Page, _Thermodynamics of Black Holes in anti-De Sitter Space_ , _Commun. Math. Phys._ 87 (1983) 577.
  * (2) E. Witten, _Anti-de Sitter space, thermal phase transition, and confinement in gauge theories_ , _Adv. Theor. Math. Phys._ 2 (1998) 505 [hep-th/9803131].
  * (3) J.M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231 [hep-th/9711200].
  * (4) E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253 [hep-th/9802150].
  * (5) S.S. Gubser, I.R. Klebanov and A.M. Polyakov, _Gauge theory correlators from noncritical string theory_ , _Phys. Lett. B_ 428 (1998) 105 [hep-th/9802109].
  * (6) S.S. Gubser, _Breaking an Abelian gauge symmetry near a black hole horizon_ , _Phys. Rev. D_ 78 (2008) 065034 [0801.2977].
  * (7) S.A. Hartnoll, C.P. Herzog and G.T. Horowitz, _Building a Holographic Superconductor_ , _Phys. Rev. Lett._ 101 (2008) 031601 [0803.3295].
  * (8) S.A. Hartnoll, C.P. Herzog and G.T. Horowitz, _Holographic Superconductors_ , _JHEP_ 12 (2008) 015 [0810.1563].
  * (9) P. Breitenlohner and D.Z. Freedman, _Positive Energy in anti-De Sitter Backgrounds and Gauged Extended Supergravity_ , _Phys. Lett. B_ 115 (1982) 197.
  * (10) P. Breitenlohner and D.Z. Freedman, _Stability in Gauged Extended Supergravity_ , _Annals Phys._ 144 (1982) 249.
  * (11) A. Ishibashi and R.M. Wald, _Dynamics in nonglobally hyperbolic static space-times. 3. Anti-de Sitter space-time_ , _Class. Quant. Grav._ 21 (2004) 2981 [hep-th/0402184].
  * (12) I.R. Klebanov and E. Witten, _AdS / CFT correspondence and symmetry breaking_ , _Nucl. Phys. B_ 556 (1999) 89 [hep-th/9905104].
  * (13) E. Witten, _Multitrace operators, boundary conditions, and AdS / CFT correspondence_ , hep-th/0112258.
  * (14) M. Berkooz, A. Sever and A. Shomer, _’Double trace’ deformations, boundary conditions and space-time singularities_ , _JHEP_ 05 (2002) 034 [hep-th/0112264].
  * (15) E. Witten, _SL(2,Z) action on three-dimensional conformal field theories with Abelian symmetry_ , in _From Fields to Strings: Circumnavigating Theoretical Physics: A Conference in Tribute to Ian Kogan_ , pp. 1173–1200, 7, 2003 [hep-th/0307041].
  * (16) D. Marolf and S.F. Ross, _Boundary Conditions and New Dualities: Vector Fields in AdS/CFT_ , _JHEP_ 11 (2006) 085 [hep-th/0606113].
  * (17) G. Compere and D. Marolf, _Setting the boundary free in AdS/CFT_ , _Class. Quant. Grav._ 25 (2008) 195014 [0805.1902].
  * (18) T. Katagiri and T. Harada, _Stability of small charged anti-de Sitter black holes in the Robin boundary_ , _Class. Quant. Grav._ 38 (2021) 135026 [2006.10301].
  * (19) S. Kinoshita, T. Kozuka, K. Murata and K. Sugawara, _to appear_ , .
  * (20) C. Dappiaggi, H.R.C. Ferreira and C.A.R. Herdeiro, _Superradiance in the BTZ black hole with Robin boundary conditions_ , _Phys. Lett. B_ 778 (2018) 146 [1710.08039].
  * (21) J.D. Bekenstein, _Extraction of energy and charge from a black hole_ , _Phys. Rev. D_ 7 (1973) 949.
  * (22) S.W. Hawking and H.S. Reall, _Charged and rotating AdS black holes and their CFT duals_ , _Phys. Rev. D_ 61 (2000) 024014 [hep-th/9908109].
  * (23) N. Uchikata and S. Yoshida, _Quasinormal modes of a massless charged scalar field on a small Reissner-Nordström-anti-de Sitter black hole_ , _Phys. Rev. D_ 83 (2011) 064020 [1109.6737].
  * (24) S.R. Green, S. Hollands, A. Ishibashi and R.M. Wald, _Superradiant instabilities of asymptotically anti-de Sitter black holes_ , _Class. Quant. Grav._ 33 (2016) 125022 [1512.02644].
  * (25) O.J.C. Dias and R. Masachs, _Hairy black holes and the endpoint of AdS 4 charged superradiance_, _JHEP_ 02 (2017) 128 [1610.03496].
  * (26) T. Ishii, Y. Kaku and K. Murata, _Energy extraction from AdS black holes via superradiance_ , _JHEP_ 10 (2022) 024 [2207.03123].
  * (27) P. Basu, C. Krishnan and P.N. Bala Subramanian, _Phases of Global AdS Black Holes_ , _JHEP_ 06 (2016) 139 [1602.07211].
  * (28) K. Maeda, S. Fujii and J.-i. Koga, _The final fate of instability of Reissner-Nordström-anti-de Sitter black holes by charged complex scalar fields_ , _Phys. Rev. D_ 81 (2010) 124020 [1003.2689].
  * (29) P. Basu, J. Bhattacharya, S. Bhattacharyya, R. Loganayagam, S. Minwalla and V. Umesh, _Small Hairy Black Holes in Global AdS Spacetime_ , _JHEP_ 10 (2010) 045 [1003.3232].
  * (30) S. Bhattacharyya, S. Minwalla and K. Papadodimas, _Small Hairy Black Holes in $AdS_{5}xS^{5}$_, _JHEP_ 11 (2011) 035 [1005.1287].
  * (31) O.J.C. Dias, P. Figueras, S. Minwalla, P. Mitra, R. Monteiro and J.E. Santos, _Hairy black holes and solitons in global $AdS_{5}$_, _JHEP_ 08 (2012) 117 [1112.4447].
  * (32) S.A. Gentle, M. Rangamani and B. Withers, _A Soliton Menagerie in AdS_ , _JHEP_ 05 (2012) 106 [1112.3979].
  * (33) H.-S. Liu, H. Lu and Y. Pang, _Revisiting the AdS Boson Stars: the Mass-Charge Relations_ , _Phys. Rev. D_ 102 (2020) 126008 [2007.15017].
  * (34) S.-F. Guo, H.-S. Liu, H. Lü and Y. Pang, _Large-charge limit of AdS boson stars with mixed boundary conditions_ , _JHEP_ 04 (2021) 220 [2101.00017].
  * (35) S. Hellerman, D. Orlando, S. Reffert and M. Watanabe, _On the CFT Operator Spectrum at Large Global Charge_ , _JHEP_ 12 (2015) 071 [1505.01537].
  * (36) A. Chamblin, R. Emparan, C.V. Johnson and R.C. Myers, _Holography, thermodynamics and fluctuations of charged AdS black holes_ , _Phys. Rev. D_ 60 (1999) 104026 [hep-th/9904197].
  * (37) V. Balasubramanian, P. Kraus and A.E. Lawrence, _Bulk versus boundary dynamics in anti-de Sitter space-time_ , _Phys. Rev. D_ 59 (1999) 046003 [hep-th/9805171].
  * (38) T. Hertog and K. Maeda, _Black holes with scalar hair and asymptotics in N = 8 supergravity_ , _JHEP_ 07 (2004) 051 [hep-th/0404261].
  * (39) T. Hertog and G.T. Horowitz, _Designer gravity and field theory effective potentials_ , _Phys. Rev. Lett._ 94 (2005) 221301 [hep-th/0412169].
  * (40) P. Basu, P. Kumar and C. Mosomane, _Incomplete RG: Hawking-Page transition, C-theorem and relevant scalar deformations of global AdS_ , 2112.07750.
  * (41) P. Bizoń, D. Hunik-Kostyra and M. Maliborski, _AdS Robin solitons and their stability_ , _Class. Quant. Grav._ 37 (2020) 105010 [2001.03980].
  * (42) T. Faulkner, G.T. Horowitz and M.M. Roberts, _New stability results for Einstein scalar gravity_ , _Class. Quant. Grav._ 27 (2010) 205007 [1006.2387].
  * (43) T. Faulkner, G.T. Horowitz and M.M. Roberts, _Holographic quantum criticality from multi-trace deformations_ , _JHEP_ 04 (2011) 051 [1008.1581].
  * (44) Y. Ahn, M. Baggioli, K.-B. Huh, H.-S. Jeong, K.-Y. Kim and Y.-W. Sun, _Holography and magnetohydrodynamics with dynamical gauge fields_ , _JHEP_ 02 (2023) 012 [2211.01760].
  * (45) H.-S. Jeong, M. Baggioli, K.-Y. Kim and Y.-W. Sun, _Collective dynamics and the Anderson-Higgs mechanism in a bona fide holographic superconductor_ , _JHEP_ 03 (2023) 206 [2302.02364].
  * (46) M. Nozawa and T. Torii, _New family of C metrics in N=2 gauged supergravity_ , _Phys. Rev. D_ 107 (2023) 064064 [2211.06517].
  * (47) I. Papadimitriou, _Multi-Trace Deformations in AdS/CFT: Exploring the Vacuum Structure of the Deformed CFT_ , _JHEP_ 05 (2007) 075 [hep-th/0703152].
  * (48) M.M. Caldarelli, A. Christodoulou, I. Papadimitriou and K. Skenderis, _Phases of planar AdS black holes with axionic charge_ , _JHEP_ 04 (2017) 001 [1612.07214].
  * (49) V. Balasubramanian and P. Kraus, _A Stress tensor for Anti-de Sitter gravity_ , _Commun. Math. Phys._ 208 (1999) 413 [hep-th/9902121].
  * (50) S. de Haro, S.N. Solodukhin and K. Skenderis, _Holographic reconstruction of space-time and renormalization in the AdS / CFT correspondence_ , _Commun. Math. Phys._ 217 (2001) 595 [hep-th/0002230].
  * (51) G. Allwright and D.M. Jacobs, _Robin boundary conditions are generic in quantum mechanics_ , 1610.09581.
  * (52) T. Ishii and K. Murata, _Floquet superconductor in holography_ , _Phys. Rev. D_ 98 (2018) 126005 [1804.06785].
  * (53) L. Li, _On Thermodynamics of AdS Black Holes with Scalar Hair_ , _Phys. Lett. B_ 815 (2021) 136123 [2008.05597].
  * (54) R.M. Wald, _Black hole entropy is the Noether charge_ , _Phys. Rev. D_ 48 (1993) R3427 [gr-qc/9307038].
  * (55) V. Iyer and R.M. Wald, _Some properties of Noether charge and a proposal for dynamical black hole entropy_ , _Phys. Rev. D_ 50 (1994) 846 [gr-qc/9403028].
  * (56) H. Lü, Y. Pang and C.N. Pope, _AdS Dyonic Black Hole and its Thermodynamics_ , _JHEP_ 11 (2013) 033 [1307.6243].
  * (57) H.-S. Liu and H. Lü, _Scalar Charges in Asymptotic AdS Geometries_ , _Phys. Lett. B_ 730 (2014) 267 [1401.0010].
  * (58) H. Lu, C.N. Pope and Q. Wen, _Thermodynamics of AdS Black Holes in Einstein-Scalar Gravity_ , _JHEP_ 03 (2015) 165 [1408.1514].