# Stability analysis on charged black hole with non-linear complex scalar

Zhan-Feng Mai<EMAIL_ADDRESS>Run-Qiu Yang<EMAIL_ADDRESS>Center for
Joint Quantum Studies and Department of Physics, School of Science, Tianjin
University, Yaguan Road 135, Jinnan District, 300350 Tianjin, P. R. China

###### Abstract

It has been shown recently that the charged black hole can be scalarized if
Maxwell field minimally couples with a complex scalar which has nonnegative
nonlinear potential. We firstly prove that such scalarization cannot be a
result of continuous phase transition for general scalar potential.
Furthermore, we numerically find that it is possible that the RN black hole
will be scalarized by a first order phase transition spontaneously and near
extremal RN black hole is not stable in micro-canonical ensemble. In addition,
considering a massless scalar perturbation, we compute the quasi-normal modes
of the scalarized charged black hole and the results do not only imply that
the spontaneously scalarized charged black hole is favored in thermodynamics
but also suggest that it is kinetically stable against scalar perturbation at
linear level. Our numerical results also definitely gives negative answer to
Penrose-Gibbons conjecture and two new versions of Penrose inequality in
charged case are suggested.

##  I Introduction

Recently, the observation of gravitational wave and black hole image provide a
new motivation and vision in black hole physicsAbbott:2016blz ;
Akiyama:2019cqa . Providing a deep understanding in quantum gravity, black
hole physics will always be a long-live project to study, theoretically and
experimentally. One of the most important topics on black hole physics is no-
hair theorem, which claims that black holes are determined only by mass,
charge and angular momentum respectivelyRuffini:1971bza . It was first
concluded by Bekenstein that the static and spherical neutral black hole in
asymptotic flat spacetime cannot be endowed with real scalar field, Proca
field and spin-2 fieldBekenstein:1972ny . Furthermore, the no-hair theorem was
extended by Mayo and Bekestein that the static and spherical black hole cannot
be endowed with a coupling charged scalar field together with a non-negative
self-interacting potential Mayo:1996mv . However, the requirements of above
no-hair theorems are as strong as its conclusion. Therefore, it is not
difficult to find hairy black hole solutions if one breaks the requirements of
these no-hair theorems. In Refs. Feng:2013tza ; Liu:2013gja ; Fan:2015tua ;
Khodadi:2020jij , taking account into a non positive definite potential or
asymptotic AdS spacetime, $D$-dimension scalarized black holes have been
constructed in Einstein theory minimally coupled with neutral scalar field.
Recently, in the extended Maxwell theory non-minimally coupled with a scalar
field, such as Einstein-Maxwell-Dilaton (EMD) theoryFan:2015oca , Einstein-
Maxwell-Scalar theoryKonoplya:2019goy , Einstein-Born-Infeld (EBI)
theoryWang:2020ohb ; Stefanov:2007eq , and Quasi-topological Electromagnetism
theoryMyung:2020ctt , a class of charged black hole with scalar hair has been
found. For a recent review relevant to the no-hair theorem, See Ref.
Herdeiro:2015waa .

The above examples give the simplest case of the black hole with an additional
scalar hair. As we have mentioned, the Bekenstein’s no-scalar-hair theorem was
generalized by Ref. Mayo:1996mv , claiming a very strong conclusion: the non-
extremal static and spherical charged black hole cannot carry charged scalar
hair whether minimally or nonminimally coupled to gravity, and with a regular
positive semidefinite self-interaction potential. However, a numerical charged
black hole solution with Q-hair was found in Einstein-Maxwell gravity
minimally coupled with a non-linear complex scalar, the self interacting
potential taking the following polynomial function Hong:2020miv

$V(|\psi|^{2})=\frac{m^{2}}{2}|\psi|^{2}-\frac{\lambda}{4}|\psi|^{4}+\frac{\beta}{6}|\psi|^{6}\,.$
(1)

It was also pointed out that the detailed form of the potential
$V(|\psi|^{2})$ is not crucial and the Q-hair can exist for large class of
nonlinear potentialHong:2019mcj ; Herdeiro:2020xmb . The reason why Ref.
Mayo:1996mv  obtained a wrong statement is that it omitted a scalar mass term
at an asymptotic infinity. Though the black holes with scalar hair has been
found in many physical models, the discovery of Refs. Hong:2020miv  has a few
of special interesting aspects. Particularly, in this model, the gravity,
Maxwell field and complex scalar field are all minimally coupled with each
other and the potential of scalar field is positive semidefinite. This gives a
possibility to realize the scalar hairy black hole in the Einstein gravity and
asymptotically flat spacetime. Furthermore, as the potential of scalar field
is nonnegative, such model will have stable true vacuum state $\psi=0$ in
Minkowski spacetime.

It needs to note that the Reissner-Nordström black hole (RN black hole) is
still a solution of field equations even if the Q-hair appears in the models
discussed by Refs. Hong:2020miv ; Hong:2019mcj . Given a scalarized charged
black hole solution, it is worth to investigate whether the scalarized charged
black hole is more stable than RN black hole in thermodynamics, i.e. whether
the RN black hole can be spontaneously scalarized by phase transition.
Recently, it has been well studied the thermodynamic self-stability associated
with heat capacity in various ensembles in Caldarelli:1999xj ; Mo:2013sxa ;
Zhang:2018rlv ; Quevedo:2006xk ; Quevedo:2013pba . However, we focus on
investigating the thermodynamic stability of scalarized charged black holes,
compared with the RN black hole in various ensembles. Concretely, in
microcanonical ensemble, given the same ADM mass $M$ and total charge $Q$, the
stability requires the maximum of the black hole entropy $S(M,Q)$ and a phase
is more stable than the other if it has larger entropy. In canonical ensemble,
with fixing Hawking temperature $T$ and $Q$, the black hole owning less
Helmholtz free energy $F(T,Q)$ will be the more stable, while in grand
canonical ensemble, with the identical $T$ and chemical potentials $\mu$, the
black hole which has smaller Gibbs free energy $G(T,\mu)$ indicates it will be
more stable. Furthermore, the black hole entropy $S(M,Q)$ is given by the area
of horizon according to Bekenstein’s entropy formula while the thermodynamic
potentials $F(T,Q)$ and $G(T,\mu)$ can be read off from partition function via
Euclidean path-integral approach developed by Hawking and York, et al.
Gibbons:1976ue ; Brown:1989fa ; Braden:1990hw .

Remarkably, in astrophysics, a real black hole practically is more closer to
grand canonical ensemble allowing charge and energy exchange with other matter
in universe. However, as theoretical research interest, we still take the case
of microcanonical ensemble and canonical ensemble under our consideration to
study whether the RN black hole will be scalarzied under the process of
discontinuous phase transition spontaneously. Furthermore, by means of
analyzing the corresponding entropy or thermodynamic potentials in various
ensembles, we find that there is a possibility that the scalarized charged
black hole is more stable than the RN black hole in thermodynamics. This
implies that, for a large class of nonlinear complex model, it is possible
that the RN black hole will be scalarized spontaneously via a first order
phase transition. Our numerical results also imply that the scalarized black
hole is always more stable than the RN black hole in microcanonical ensemble
and when $M\rightarrow|Q|$. This implies the isolated near extremal RN black
hole is not stable and will spontaneously scalarize.

In addition to the thermodynamic stability, another natural question is
whether the scalarized charge black hole is kinetically stable, i.e. stable
under against a small perturbation at least at linear level. It has been
developed by Vishveshwara and Teukolsky et al.Teukolsky:1973ha ; Press:1973zz
; Vishveshwara:1970cc  that in static axial or spherical symmetric background,
the equation of motion associated with the perturbation field can reduce to
radial equation in frequency domain. Furthermore, the radial equation can be
interpreted as an eigenvalue problem. Specifically, imposing a physical
boundary condition both in the spacial infinity and at the black hole horizon,
a class of complex frequency called black hole quasi-normal modes(QNMs), which
implies dissipation at both event horizon and spatial infinity, can be picked
out. The instability of black hole might be triggered due to the negative
imaginary part of the complex frequency. In general, the QNMs can be solved by
shooting method in numerics. Approximatively, there have been some analytical
approaches to achieve to calculate QNMs as well: the WKB approximation method,
continue fraction method and Monodromy method. However, recently it has been
argued in Ref. Konoplya:2019hlu  that the WKB approximation method cannot
catch the unstable mode within the spectrum of the QNMs. Furthermore, we also
find that the shooting method cannot efficiently give the stable modes through
numerical error analysis. We thus adopt hybrid method to calculate the QNMs
for studying the kinetic stability of scalarized charged black hole. If a
given solution is not kinetically stable, then such solution cannot exist in a
real physical system. For the situations that the scalarized black hole is
more thermodynamically stable than RN black hole, but it does not
automatically guarantee that the scalarized charged black hole is still stable
kinetically. However, our results show the neutral perturbative scalar field
will not trigger the instability of the scalarized charged black hole at
linear level. In addition, for a recent review on perturbation theory on black
hole, stability analysis and quasi-normal modes of black hole, more detail has
been given in Refs. Pani:2013pma ; Konoplya:2011qq ; Berti:2009kk .

This paper will be organized as follow: In Sec. II, we briefly introduce the
model of Einstein-Maxwell theory minimally coupled with a non-linear complex
field. For general non-linear semi-definite potential, we present a proof that
the scalarized charged black hole cannot be a result of continuously
scalarization. In Sec. III, Giving a logarithmic potential, we obtain a class
of numerical scalarized charged black hole solution and in various ensembles,
investigate their thermodynamic stability compared with the RN black hole. We
find that in both microcanonical ensemble and canonical ensemble, it is
possible that the RN black hole can be spontaneously scalarized via first
order phase transition. In Sec. IV, considering a probing neutral scalar
field, we investigate the stability on the scalarized charged black hole by
means of computing the quasi-normal modes using both the shooting method and
the WKB approximation method. This suggests that the probing neutral scalar
field cannot trigger the instability of the sclarized charged black hole. In
Sec. V, we present our conclusion and further discussions.

##  II Model Setup

In this paper, we consider the following Einstein-Maxwell theory minimally
coupled with a non-linear complex scalar field (we set
$G_{N}=c_{s}=\hbar=k=1$, $c_{s}$ denotes the speed of light),

$\displaystyle S=\frac{1}{16\pi}\int\mathrm{d}^{4}x\sqrt{-g}{\cal L}\,,$ (2)
(4) $\displaystyle{\cal
L}=\left(R-F_{\mu\nu}F^{\mu\nu}-(D_{\mu}\psi)^{\dagger}(D^{\mu}\psi)-W(|\psi|^{2})\right)\,,$

where the covariant derivative operator $D_{\mu}=\nabla_{\mu}-iqA_{\mu}$, $q$
denotes charge of the complex scalar field. In addition,
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the strength
tensor of the U(1) electromagnetic field $A_{\mu}$. The non-linear positively
semi-definite potential $W(x)$ is a smooth function and satisfies following
requirements:

$W(0)=0,~{}~{}W^{\prime}(0)=m^{2}>0,~{}~{}W(x)>0~{}\text{if}~{}x>0\,.$ (5)

These conditions insure that scalar field has stable true vacuum $\psi=0$ in
Minkowski spacetime when $F_{\mu\nu}=0$. Performing variation on the action
(2) with respect to the metric $g_{\mu\nu}$, the gauge field $A_{\mu}$ and the
complex scalar field $\psi$ respectively, we obtain the equation of motion

$\displaystyle
R_{\mu\nu}-\frac{R}{2}g_{\mu\nu}=8\pi[T^{(M)}_{\mu\nu}+T^{(\psi)}_{\mu\nu}]\,,$
(6) (7)
$\displaystyle\nabla^{\mu}F_{\mu\nu}=iq\left(\psi^{\dagger}D_{\mu}\psi-\psi(D_{\mu}\psi)^{\dagger}\right)\,,$
(8) (9) $\displaystyle D^{2}\psi-w(|\psi|^{2})\psi=0\,,$ (10)

where we define $w(x)=W^{\prime}(x)$ and the energy momentum tensor associated
with both the electric field and the complex scalar field reads

$\displaystyle
T^{(M)}_{\mu\nu}=\frac{1}{4\pi}\left(F_{\mu\sigma}{F^{\sigma}}_{\nu}-\frac{1}{4}g_{\mu\nu}F_{\rho\sigma}F^{\rho\sigma}\right)\,,$
(11) (14) $\displaystyle
T^{(\psi)}_{\mu\nu}=\frac{1}{16\pi}\left(D_{\mu}\psi(D_{\nu}\psi)^{\dagger}+D_{\nu}\psi(D_{\mu}\psi)^{\dagger}\right.$
$\displaystyle\left.-g_{\mu\nu}\left((D_{\mu}\psi)^{\dagger}(D^{\mu}\psi)+W(|\psi|^{2})\right)\right)\,.$

For the non-linear potential, the requirement (5) implies that

$w(x)=m^{2}+w_{1}x+w_{2}x^{2}+\cdots$ (15)

with some constants $\\{w_{1},w_{2},\cdots\\}$ when $|x|\ll 1$.

To obtain the scalary charged black hole and investigate some related
properties, we adopt the following spherical symmetric line element,

$\mathrm{d}s^{2}=-f(r)\mathrm{e}^{-\chi(r)}\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{f(r)}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\varphi^{2}).$
(16)

Due to the spherical symmetry, the electromagnetic field and the complex
scalar field take the following form,

$A_{\mu}=\phi(r)(\mathrm{d}t)_{\mu},~{}~{}\psi=\psi(r)\,.$ (17)

From the line element Eq. (16), the Hawking temperature $T$ reads

$T=\frac{f^{\prime}(r_{h})\mathrm{e}^{-\chi(r_{h})/2}}{4\pi}\,,$ (18)

where $r_{h}$ denotes the outer event horizon. With the given spherical
anstaz, the equations of motion reduce to

$\begin{split}&\psi^{\prime\prime}+\left(\frac{f^{\prime}}{f}-\frac{\chi^{\prime}}{2}+\frac{2}{r}\right)\psi^{\prime}+\left(\frac{\mathrm{e}^{\chi}\phi^{2}q^{2}}{f^{2}}-\frac{w(\psi^{2})}{f}\right)\psi=0\\\
&\phi^{\prime\prime}+(\chi^{\prime}/2+2/r)\phi^{\prime}-\frac{\psi^{2}q^{2}}{2f}\phi=0\,,\\\
&\chi^{\prime}+\frac{r\mathrm{e}^{\chi}\psi^{2}\phi^{2}q^{2}}{f^{2}}+r\psi^{\prime
2}=0\,,\\\ &f^{\prime}+\left(\frac{1}{r}+\frac{r\psi^{\prime
2}}{2}\right)f+r\mathrm{e}^{\chi}\phi^{\prime 2}\\\
&-\frac{1}{r}+\frac{\mathrm{e}^{\chi}r\psi^{2}\phi^{2}q^{2}}{2f}+\frac{1}{2}rW(|\psi|^{2})=0\,,\end{split}$
(19)

where the prime denotes the derivative with respect to $r$. In this paper, we
consider the asymptotically flat space time. The scalar field, gauge field and
metric components should satisfy the following regular boundary conditions
when $r\rightarrow\infty$

$\displaystyle
f=1-\frac{2M}{r}+\cdots,~{}~{}\chi=\frac{\chi_{2}}{r^{2}}+\cdots\,,$ (20) (21)
$\displaystyle\phi=\mu-\frac{Q}{r}+\cdots,~{}~{}|\psi|\leq\mathcal{O}(1/r^{2})\,,$
(22)

where $\mu$ is the chemical potential, $M$ is the ADM mass and $Q$ is the
total charge. With this boundary conditions and noting the fact that
$w(|\psi|^{2})\rightarrow m^{2}$ near the boundary, we find that the first
equation of Eq. (19) reduces into following simple form near the infinity

$\psi^{\prime\prime}+\frac{2\psi^{\prime}}{r}+(q^{2}\mu^{2}-m^{2})\psi=0\,.$
(23)

The solution reads

$\psi(r)=\frac{\psi_{+}}{r}\mathrm{e}^{r\sqrt{m^{2}-q^{2}\mu^{2}}}+\frac{\psi_{-}}{r}\mathrm{e}^{-r\sqrt{m^{2}-q^{2}\mu^{2}}}\,.$
(24)

The boundary conditions Eq. (20) implies following constraints

$m^{2}-q^{2}\mu^{2}>0,~{}~{}\psi_{+}=0\,.$ (25)

These are two necessary conditions on the spontaneous scalarization for
asymptotically flat black holes.

Taking into account a polynomial potential, as mentioned above, a numerical
charged black hole solution with Q-hair has been found in Mayo:1996mv ;
Hong:2020miv . A natural question arises whether this class of scalarized
black hole solution will arise from continuously spontaneous scalarization for
specific non-linear potential. This question is important because if the
answer is yes, the spacetime geometry can transit into scalarized black hole
smoothly, otherwise, the latent heat will be relaxed or absorbed during phase
transition between the RN black hole and scalarized black hole. Moreover, such
latent heat will leave some observable effects if such a phase transition
happened in our universe.

In the following, we shall present a proof that for arbitrary non-linear
potential $W(|\psi|^{2})$ satisfying the requirement (5), the spontaneous
scalarization of RN black hole cannot happen via continuous phase transition.
As shown in Fig. 1, if such continuous phase transition can happen, when we
tune the parameters (total charge, chemical potential, temperature, etc.) of
the black hole, there is a critical point where the strength of scalar field
begins to increase into nonzero continuously. In other words, there must be a
small region near the critical value associated with the black hole parameters
(See also the top panel of Fig. 1) where the complex scalar field is
infinitesimal.

Figure 1: Top: Scalarization appear via continuous phase transition. Bottom:
Scalarization appear via discontinuous phase transition.

Without losing generality, we can assume this infinitesimal scalar field has
following form,

$\psi(r,\theta,\varphi)=\varepsilon\sum_{l,\nu}a_{l\nu}(r)Y_{l\nu}(\theta,\varphi)\,\quad\varepsilon\rightarrow
0\,.$ (26)

where $Y_{l\nu}$ denotes the spherical harmonic function, $\nu=0,\pm 1,\pm
2,\cdots$ is the magnetic quantum number and $l>\nu$ is the azimuthal quantum
number. Taking it into Eq. (6) and neglecting all the non-linear terms of
$\varepsilon$, we find that the metric and gauge field decouple from the
scalar field. Then spacetime geometry and gauge field are given by a RN
solution

$\displaystyle\chi=0,~{}~{}f(r)=\frac{(r-r_{h})(r-\mu^{2}r_{h})}{r^{2}}\,,$
(27) (28) $\displaystyle\phi=\mu(1-r_{h}/r),~{}~{}\mu\in[-1,1]\,.$ (29)

Here we require $\mu\in[-1,1]$ since the $r_{h}$ is defined to be the most
outer horizon. The E.O.M of scalar field reads

$a_{l\nu}^{\prime\prime}+\left(\frac{f^{\prime}}{f}+\frac{2}{r}\right)a^{\prime}+\left(\frac{\phi^{2}q^{2}}{f^{2}}-\frac{m^{2}}{f}-\frac{l(l+1)}{fr^{2}}\right)a_{l\nu}=0\,.$
(30)

This can be rewritten into following form

$\frac{\mathrm{d}}{\mathrm{d}r}\left(r^{2}fa_{l\nu}\frac{\mathrm{d}a_{l\nu}}{\mathrm{d}r}\right)=\left[\left(\frac{\tilde{m}(r)^{2}}{f}-\frac{\phi^{2}q^{2}}{f^{2}}\right)a_{l\nu}^{2}+a_{l\nu}^{\prime
2}\right]r^{2}f\,.$ (31)

where we have redefine the effective mass as
$\tilde{m}(r)^{2}=m(r)^{2}+l(l+1)/r^{2}$ and it is obviously that
$\tilde{m}^{2}\geq m^{2}$. Integrate it from horizon to infinity and we find

$\left.r^{2}fa_{l\nu}\frac{\mathrm{d}a_{l\nu}}{\mathrm{d}r}\right|_{r_{h}}^{\infty}=\int_{r_{h}}^{\infty}\left[\left(\frac{\tilde{m}(r)^{2}}{f}-\frac{\phi^{2}q^{2}}{f^{2}}\right)a_{l\nu}^{2}+a_{l\nu}^{\prime
2}\right]r^{2}f\mathrm{d}r\,.$ (32)

As $a(r)$ is a regular at horizon and decays to zero at infinity, the left
side of Eq. (32) is zero. Thus, we have

$\int_{r_{h}}^{\infty}\left[\left(\frac{\tilde{m}(r)^{2}}{f}-\frac{\phi^{2}q^{2}}{f^{2}}\right)a_{l\nu}^{2}+a_{l\nu}^{\prime
2}\right]r^{2}f\mathrm{d}r=0\,.$ (33)

This means that $\exists r_{0}\in(r_{h},\infty)$ such that

$\frac{\tilde{m}(r)^{2}}{f(r_{0})}-\frac{\phi(r_{0})^{2}q^{2}}{f(r_{0})^{2}}<0\,.$
(34)

We then obtain

$m^{2}\leq\tilde{m}(r)^{2}<\frac{\phi(r_{0})^{2}q^{2}}{f(r_{0})}=\mu^{2}q^{2}\frac{r_{0}-r_{h}}{r_{0}-\mu^{2}r_{h}}<\mu^{2}q^{2}\,.$
(35)

Here we have used the solution (27). Above result and requirement (25) are
contradictory. This shows that, if the complex scalar field appears in a
static black hole, its strength cannot be infinitesimal no matter how we
choose the black hole parameters. Thus, the continuous phase transition from
RN black hole to scalarized charged black cannot occur.

The demonstration above implies that the scalarized charged black hole
solution cannot produce from the continuous spontaneous phase transition.
However, the spontaneous scalarization on charged black hole may arise through
non-continuous phase transition, namely the first-order phase transition. In
this case, the complex scalar field jumps into nonzero from zero when we tune
the parameters of black hole, see the bottom panel of Fig. 1. Recall the
equivalence between the black hole system and the thermal mechanic system, it
is thus intriguing to study the thermodynamic stability of the scalarized
charged black hole, compared with the corresponding charged black hole with no
hair. Motivated by this, interpreted charged black hole with scalar hair as a
class of ensemble, in the following, we focus on investigating the
thermodynamical stability on scalarized charged black hole in various
ensembles: microcanonical ensemble, canonical ensemble and grand canonical
ensemble, and try to check if there is a first order phase transition.

##  III Thermodynamic Instability

In this section, we analyze the thermodynamic stability of scalarized charged
black hole by proposing a specific example. As mentioned, the relation of
Beskein entropy implies that black holes can be investigated as a
thermodynamical system in term of three different ensemble: microcanonical
ensemble, canonical ensemble and grand canonical ensemble. The microcanonical
ensemble implies that black holes are interpreted as an isolated system where
there does not exist any charge and energy exchange of black holes. The
canonical ensemble indicates that black holes only exchange energy with
environment of which the temperature is fixed. The grand canonical ensemble is
similar to canonical ensemble but also admits the black hole to exchange the
particles with environment. In practice, black holes in astrophysics are more
likely to be grand canonical ensemble with exchange of particles and energy.
However, as theoretical investigation on scalaized charged black hole, in
following we still take microcanonical ensemble, canonical ensemble and grand
canonical ensemble into account, studying the thermodynamical stability
associated with the scalarized charged black hole and RN black hole
respectively. Furthermore, we find that in grand canonical ensemble, the RN
black hole is more stable than the scalarized charged black hole in
thermodynamics, corresponding to general expectation. Nevertheless, we also
find that the discontinuous scarization on RN black hole may happen in both
microcanonical ensemble and canonical ensemble respectively. In the following,
we will give more detail discussion.

Before go on discussing our result, let us first make short comment on
different ensembles in black holes. The Euclidean path-integral approach to
black hole thermodynamics originally proposed by Hawking in microcanonical
ensemble Hawking:1976de . Later on, the canonical ensemble was investigated by
York et al. York:1986it ; Whiting:1988qr ; Brown:1994su . It was found that
suitable boundary conditions must be added in canonical ensemble. Then the
York’s approach was generalized into other ensembles such as the charged black
hole in the grand canonical ensemble PhysRevD.42.3376 ; Brown_1990 . It has
been found that the results obtained by using the path-integral approach
depend on the boundary conditions Brown:1994gs ; Hawking:1982dh . The
stability of black holes then also depends on the choice of boundary
conditions and, consequently, on the choice of ensembles Comer_1992 . In fact,
the stability properties of a black hole are drastically influenced by the
boundary conditions that determine ensemble.

In following we proceed to our discussion. We specify the non-linear potential
as a logarithmic function with respect to the scalar field

$W(x):=m^{2}c^{2}\log(1+\frac{x^{2}}{c^{2}}),$ (36)

where $c$ is a constant and $m$ is the effective mass of the scalar field.
Considering the flat directions in gauge-mediated supersymmetric model
deGouvea:1997afu ; Kusenko:1997si , this potential is proposed in Ref.
Hong:2019mcj , in which the supersymmetric breaking has been absorbed into the
rescaling of $\psi$. In addition to satisfying the stable vacuum requirement
Eq.(5) , the shape of potential is asymptotic flat when taking large field
limit $\psi\gg c$. From the equation of motion of $\psi$ in Eq. (19), one will
find that the scalar field becomes massless large field limit $\psi\gg c$.
Moreover, the logarithmic potential can bring better numerical stability as
well.

Since this paper tries to find the black hole solutions with scalar hair,
there should be a horizon locating at position $r=r_{h}$. In static
spherically symmetric case, this implies $f(r_{h})=0$. To set up the numerical
method to solve the equation of motion Eq. (19), we in principle still needs
five additional independent boundary conditions at horizon. Practically, the
regularity of physical fields at $r=r_{h}$ implies that the solution can
approximatively be written as the following Taylor’s series with respect to
$r-r_{h}$,

$\displaystyle
f(r)=f_{1}(r-r_{h})+\cdots,\quad\chi(r)=\chi_{0}+\chi_{1}(r-r_{h})+\cdots\,,$
(37) (39)
$\displaystyle\phi(r)=\phi_{0}+\phi_{1}(r-r_{h})+\phi_{2}(r-r_{h})^{2}+\cdots\,,$
(41) $\displaystyle\psi(r)=\psi_{0}+\psi_{1}(r-r_{h})+\cdots\,.$

Take them into Eq.(19) and we will find

$\displaystyle\psi_{1}=\frac{2\psi_{0}c^{2}m^{2}}{(\psi_{0}^{2}+c^{2})f_{1}}\,,\quad\phi_{0}=0\,.$
(42)

This leaves three independent variables $\\{\psi_{0},\chi_{0},\phi_{1}\\}$ at
horizon and so we have 7 independent parameters in solving Eq. (19)

$\\{r_{h},\psi_{0},\chi_{0},\phi_{1},c,m,q\\}\,,$ (43)

in which $\\{r_{h},\psi_{0},\chi_{0},\phi_{1}\\}$ come from the value of
various fields at the horizon, $q$ is the charge of the scalar field and
$\\{c,m\\}$ are the parameters of the non-linear potential. It is remarkable
that there is a scaling symmetry of the equation of motion

$r\to\lambda r\,,\quad t\to\lambda t\,,\quad q\to\frac{q}{\lambda}\,,\quad
m\to\frac{m}{\lambda},$ (44)

which equivalently rescale the metric $g_{\mu\nu}\to\lambda^{2}g_{\mu\nu}$ and
the electronic field $A_{\mu}\to\lambda A_{\mu}$. Due to such a symmetry, we
fix $m=0.01$ for convenience. In addition, we must impose

$\chi\to 0,\quad\text{when}\quad r\to\infty,$ (45)

which is a requirement about the normalisation of $t$ associated with the
gravitational redshiftHartnoll:2008kx . Practically, given any arbitrary value
of $\chi_{0}$, the equation of motion is invariant when performing the
following scaling transformation

$t\to\mathrm{e}^{-\frac{\chi_{\infty}}{2}}t\,,\quad\chi\to\chi-\chi_{\infty}\,,\quad\phi(r)\to\mathrm{e}^{\frac{\chi_{\infty}}{2}}\phi(r),$
(46)

where $\chi_{\infty}$ denotes the value of the solution $\chi(r)$ in the
asymptotic infinity. In other words, we perform a time rescaling $t\to
e^{-\frac{\chi_{\infty}}{2}}t$ to set $\chi=0$ at the boundary. In order to
numerically solve the equations simply, we in general set $\chi_{0}=0$.
However, such choice will lead to $\chi(\infty)=\chi_{\infty}\neq 0$, we can
thus finally transform the solution to satisfy $\chi(\infty)=0$ by the
transformation (46).

Base on these two symmetries, two physical parameters, the charge of the
complex scalar field $q$ and the parameter related to the scalar potential
$c$, and three parameters as the initial value at the horizon $r_{h}$,
$\phi_{1}$ and $\psi_{0}$ are left. Therefore, The integration of the equation
of motion Eq. (19) from the event horizon to the infinity will give us a map:

$\\{r_{h},\phi_{1},\psi_{0},c,q\\}\mapsto\\{\mu,T,Q,M,\psi_{+}\\}\,.$ (47)

If one chooses five parameters $\\{r_{h},\phi_{1},\psi_{0},c,q\\}$
arbitrarily, the $\psi_{+}$ may be or may not be vanish. From Eq. (25), we
know that only the one satisfying $\psi_{+}=0$ is an admitted solution,
leading to only four of parameters $\\{r_{h},\phi_{1},\psi_{0},c,q\\}$ are
independent.

We choosing the parameters as
$r_{h}=1,c=0.1,\phi_{1}=0.8,m=0.01,q=\frac{11}{10}m$ and set
$r=r_{\infty}=1000$ as the infinity cut off, finding that the complex scalar
$\psi$ has been efficiently decay to $0$ at the infinity cutoff $r=1000$, we
firstly numerically solve the equation of motion Eq. (19) by Runge-Kutta
methods under the initial condition Eq. (37) and Eq. (42). As we have
explained, the $\psi_{0}$ now is not free because we need to satisfy
$\psi_{+}=0$. Then we use the standard shooting method to find the smallest
value of $|\psi_{0}|$ to satisfy the constrain Eq. (25). We find a numerical
charged black hole solution with scalar hair when $\psi_{0}\approx 0.1988$.
After Performing the scaling transformation Eq. (46), we show the numerical
solution are shown in Fig. 2 from which we find that
$g_{tt}=f(r)\mathrm{e}^{-\chi(r)}$ is still a monotone increasing function
with respect with $r$, indicating that in scalarized charged black hole, the
gravity is still attraction, as the case of the RN black hole.

In addition, taking $r=1000$ as an infinity cutoff, from Fig. 2 one can see
that the scalar field $\psi(r)$ have efficiently decay to $0$ smoothly at
$r=1000$, implying that the spacetime manifold has efficiently reduce to RN
black hole. In addition, if increase the infinity cutoff, to maintain
effective numerical accuracy, one must increase the working precise, leading
to increase the computational time. In practice, the best cut-off is chosen by
following way: in double float accuracy, we take $\max
r=2^{6},2^{7},2^{8},2^{9}$ and $2^{10}\approx 1000$. We observed that the
differences of $\psi_{0}$ are smaller and smaller. However, if we increase the
cut-off to be $2^{11}$ and more, we found the differences of $\psi_{0}$ are
larger and larger. This implies 1000 is best cut-off. If we using quadruple
float accuracy, the best cut-off is around 2000, however, the computational
time will increase more 10 times. Therefore, for investigating the stability
of scalarized charged black hole effectively and efficiently, it is reasonable
to set $r=1000$ as an infinity cutoff.

Figure 2: Numerical solutions of the metric components
$g_{tt}=f(r)\mathrm{e}^{-\chi(r)}$, $1/g_{rr}=f(r)$, the electric potential
$\phi(r)$ and the scalar field $\psi(r)$ respectively. Here we take parameters
{$r_{h}=1,c=0.1,\phi_{1}=0.8,m=0.01,q=\frac{11}{10}m$}. The $\psi_{0}\approx
0.1988$ is given by shooting method.

###  III.1 Microcanonical Ensemble

####  III.1.1 Thermodynamical Stability for scalarized charged black hole

To compare which one is more stable, RN black hole or scalarized black hole,
we need to specify what ensemble we will consider. We first consider the
microcanonical ensemble describing an isolated system. In the microcanonical
ensemble, the characteristic thermodynamic variables is entropy and a physical
realistic process will be always towards the direction of increasing entropy,
a phase transition thus can happen only if the entropy will be increased.
Therefore, in the case microcanonical ensemble, we consider an isolated black
hole where the total mass $M$ and total charge $Q$ are fixed parameters and
the black hole entropy can be interpreted as $S=S(M,Q)$. Since the black hole
entropy is proportion to the area of black hole $S\sim 4\pi r_{h}^{2}$, which
implies the larger the black hole radius indicates more stable in
microcanonical ensemble, we need to compare horizon radii of scalarized black
hole and RN black hole for the same mass $M$ and total charge $Q$.

Following our illustration above, we numerically analyze the behavior of the
event horizon radius working with $(Q/M,\Delta r_{h}/M)$ plane with fixing
total mass $M$, where $\Delta r_{h}=r_{h}-r_{h_{\text{RN}}}$ denotes
difference between the scalarized charged black hole and the corresponding
charged black hole with the same $M$ and $Q$. Specifically, given a numerical
scalarized charged black hole solution, the ADM mass and the total charge read

$M=-\frac{r}{2}(f(r)-1)|_{r\to r_{\infty}}\,,\quad Q=r^{2}\phi(r)|_{r\to
r_{\infty}}.$ (48)

Then, the event horizon radius of the corresponding RN black hole is

$r_{h_{\text{RN}}}=M+\sqrt{M^{2}-Q^{2}}.$ (49)

Naively, we firstly consider the solution given in Fig. 2. The mass and the
total charge directly read $M\approx 11.53$ and $Q\approx 11.81$, implying
that there does not exist the corresponding RN black hole sharing the same
mass and charge. Furthermore, it also indicates that due to the non-linear
complex scalar field, the mass of the scalarized charged black hole could
smaller than its total charge.

In following we search the parameter region in which the mass of the
scalarized charged black hole is larger than its charge using shooting method
where both scalarized black hole (if exists) and RN black hole are solutions
of Eq. (19). In this paper, we mainly adopt the shooting method in numerics
for investigating the thermodynamical stability associated with the scalarized
charged black hole. We give a plot in $(c,\psi_{0})$ plane (See also the first
plot of Fig. 3), where every point on the curve denotes a numerical solution
of charged black hole with scalar hair. Remarkably, from this plot one can
easily see that $\psi_{0}$ will vanish when $c\ll 1$, implying that the
scalarized charged black hole will reduce to RN black hole. Said another way,
our result is consistent with the case of the probe limit in which the action
Eq. (2) will reduce to Einstein-Maxwell theory after rescaling $\psi\to c\psi$
then taking $c\to 0$.

Recall the Eq. (48), one can pick out the ADM mass $M$ and total charge $Q$ of
every scalarized charged black hole denoted by $\\{c_{i},\psi_{0_{i}}\\}$.
Therefore, we present a plot on $(c,M-Q)$ plane in points and find $M-Q$ exist
a zero point around $c\in[0.00101,0.02]$, (See also the second plot of Fig.
3). Using both the interpolation method and the shooting method, we obtain a
series of numerical scalarized charged black hole solution with $M=Q\approx
0.9676$, denoted by $\\{c\approx 0.01220,\psi_{0}\approx 0.02473\\}$, which is
actually a good seed numerical solution for for investigating the relation
between $\Delta r_{h}$ and $Q$ practically. Furthermore, recall that we are
working with the microcanonical ensemble, the radius of the corresponding RN
black hole $r_{h_{RNi}}$ reads $r_{h_{RNi}}=M+\sqrt{M^{2}-Q_{i}^{2}}$. In
other words. we established a relation between $\Delta
r_{h_{i}}=r_{h_{i}}-r_{h_{RNi}}$ and $Q_{i}$ (See also the last plot in Fig. 3
and the Appendix. B.1 for detail discussion.).

Figure 3: Top: Relationship between $c$ and $\psi_{0}$ based on our numerical
solutions. Middle: A typical example about the relation between $M-Q$ and $c$.
Bottom: The relation between $\Delta r_{h}/M$ and $Q/M$ with fixing
$M=0.9676$.

From $S\sim r_{h}^{2}$ and the plot presented in $(Q/M,\Delta r_{h}/M)$ plane,
one can find that there does exist a small interval $Q/M\in[0.9993,1]$ in
which the entropy of the scalarized charged black hole is larger than the RN
black hole in microcaonical ensemble. To conclude, in microcanonical ensemble,
the scalarized charged black hole with the mass-charge ratio lying on the
interval $0.9993<Q/M<1$ will more stable than the RN black hole. We note the
when $|Q|\sim M$, i.e. the corresponding RN black hole approaches to be
extreme, the scalarized black hole has larger entropy. We have carefully
checked and find that this is also true if the parameters are taken to be
other different values. Thus, our numerical results imply that the near
extreme RN black hole is not stable and will be spontaneously scalarized via a
first order phase transition. Remarkably, when $Q>M$, there still exists
scalarized charged black hole solutions but the corresponding RN black hole
does not.

####  III.1.2 New versions for Penrose-Gibbons conjecture

In fact our numerical result does also definitely give a negative answer to a
long-standing conjecture named Penrose-Gibbons conjecture. It was conjectured
that, for all asymptotically flat black holes, the area of horizon $A_{h}$,
the total mass $M$ and the total charge $Q$ will satisfy following inequality
Mars:2009cj

$M\geq\sqrt{\frac{A_{h}}{16\pi}}+Q^{2}\sqrt{\frac{\pi}{A_{h}}}\,$ (50)

or a weaker version

$\sqrt{\frac{A_{h}}{16\pi}}\leq\frac{1}{2}\left(M+\sqrt{M^{2}-Q^{2}}\right)\,$
(51)

if the weak energy condition is satisfied. The saturation will appear only in
RN black holes. pThese two inequalities are charged generalization of
following Penrose inequality

$M\geq\sqrt{\frac{A_{h}}{16\pi}}$ (52)

Eq. (52) has been proven in general static case by a few of different methods
Mars:2009cj . However, the proofs of charged generalizations (50) and (51) are
still open. In spherical case, they reduce to

$M\geq\frac{r_{h}}{2}+\frac{Q^{2}}{2r_{h}}\,$ (53)

and

$r_{h}\leq\left(M+\sqrt{M^{2}-Q^{2}}\right)\,.$ (54)

and the weak energy condition reduces into the requirement of $T_{00}\geq 0$
outside event horizon. In our model, recall the energy momentum tensor given
in Eq. (11), the $tt$ component of energy momentum tensor in spherical anstanz
Eq. (16) reads

$T_{tt}=\frac{1}{2}q^{2}\phi^{2}\psi^{2}+\frac{1}{2}\mathrm{e}^{-\chi}fV(\phi)+\frac{1}{4}f(r)\phi^{\prime
2}+\frac{1}{2}f^{2}\psi^{2}\,.$ (55)

For a semi-definite non-linear potential $V(\phi)$, it can obviously to
observe that the weak energy condition $T_{00}\geq 0$ is always satisfied. In
addition, The inequalities (50) and (51) are two generalization of Penrose
inequality in charged case. Though the Penrose inequality in static case
(which is called Riemannian-Penrose inequality) has several proofs, the
strength version (50) has not been prove even in static spherical case. A
serval proofs have been obtained by assuming that outside the black hole there
are no charge current sources i.e. in electrovacuum and horizon is connected,
see Refs. Malec:1994sy ; Hayward:1998jj ; Gibbons:1998zr ; Khuri2013 . For
inequality (50), a couterexample was found when the horizon is not connected
Weinstein2005 . Though the proof of inequality (51) has not been obtained yet
in general case, as far as we know, no counterexample of weaker version (51)
was reported. Base on the inequality of arithmetic and geometric means
$\frac{a+b}{2}\geq\sqrt{ab}$, a natural deduction of Eq. (53) is $M\geq Q$,
satisfying in RN black hole. For inequality (54) to make sense it is necessary
that the spacetime should satisfies $M\geq|Q|$.

However, our numerical results offer a counterexample for both inequalities
(53) and (54). Therefore, even in spherically symmetric case which has only
one connected horizon, the inequalities (50) and (51) can still be broken.

We note that inequalities (50) is not the only natural generalization of
original Penrose inequality. Here we offer two new generalizations. One
natural generalization of the Penrose inequality in scalarized charged black
hole is that interpreting $Q$ as the charge enclosed with event horizon
$Q_{h}$, which is identical in RN black hole due to the charge conservation
law. The other natural generalization is that we use chemical potential $\mu$
to replace the role of charge $Q$. Therefore, we have two new versions of the
generalized Penrose inequality,

$M\geq\sqrt{\frac{A_{h}}{16\pi}}+Q_{h}^{2}\sqrt{\frac{\pi}{A_{h}}}\,,\quad
M\geq\sqrt{\frac{A_{h}}{16\pi}}+\mu^{2}r_{h}^{2}\sqrt{\frac{\pi}{A_{h}}}.$
(56)

In spherical case, Eq. (56) will reduce to

$M\geq\frac{r_{h}}{2}+\frac{Q_{h}^{2}}{2r_{h}}\,,\quad
M\geq\frac{r_{h}}{2}+\frac{\mu^{2}r_{h}}{2}\,.$ (57)

Here we only numerically verify these two inequalities Eq. (57) respectively.
Recall the class of scalarized charged black hole solution denoted by
$\\{c,\psi_{0}\\}$, we work with the
$\left(c,M-\left(\frac{r_{h}}{2}+\frac{Q_{h}^{2}}{2r_{h}}\right)\right)$ plane
and $\left(c,M-\left(\frac{r_{h}}{2}+\frac{\mu^{2}r_{h}}{2}\right)\right)$
respectively, where $Q_{h}=r_{h}^{2}\mathrm{e}^{\chi(r_{h})}\phi(r_{h})$. From
the plots given in Fig. 4, we find that
$M>\frac{r_{h}}{2}+\frac{Q_{h}^{2}}{2r_{h}}$ and
$M>\frac{r_{h}}{2}+\frac{\mu_{h}^{2}r_{h}}{2}$ always hold as the growth of
$c$, implying that the two possible generalization of the Penrose inequality
we imposed hold in spherical case. Furthermore, it is an intriguing topic to
proof Eq. (54) in general scalarized charged black hole and we will leave it
as future work.

Figure 4: A numerical check on the inequalities shown in Eq. (57). Here we
take $\\{m=0.01,q=1.1m,r_{h}=1,\phi_{1}=0.8\\}$.

###  III.2 Canonical Ensemble and Grand Canonical Ensemble

In this section, we turn to investigate the thermodynamic stability on
scalarized charged black hole in both canonical ensemble and grand canonical
ensemble, by means of a series of specific numerical solution as an example.
In canonical ensemble, the associated thermodynamic potential is the Helmholtz
Free Energy $F(T,Q)$ with respect to the temperature $T$ and the total charge
$Q$. The thermodynamics of the grand canonical ensemble can be described by
the Gibbs Free Energy $G(T,\mu)$, in which the thermodynamic variables are the
temperature $T$ and the chemical potential $\mu$ respectively. In both these
two ensembles, the real physical process will be towards the direction
increasing $F(T,Q)$ or $G(T,\mu)$. Thus, the RN black hole will spontaneously
scalarize if the hairy black hole has smaller free energy.

Following the procedure developed by Gibbons:1976ue ; Caldarelli:1999xj , one
can read off the free energy from the on-shell Euclidean action, namely
$Z=\mathrm{e}^{-S_{\text{E}_{\text{on-shell}}}}$ and

$F:=-\frac{1}{\beta}\ln Z|_{\text{canonical}}\,,\quad G:=-\frac{1}{\beta}\ln
Z|_{\text{grand canonical}}\,,$ (58)

where $Z$ is the thermodynamic partition function and $\beta=\frac{1}{T}$ is
the inverse temperature. We start with the Euclidean action associated with
(2)

$S_{E}=S_{\text{bulk}}+S_{\text{surf}}+S_{\text{ct}}\,,$ (59)

where

$\displaystyle
S_{\text{bulk}}:=-\frac{1}{16\pi}\int\mathrm{d}^{4}x\sqrt{g_{E}}$ (61)
$\displaystyle\left(R-F_{\mu\nu}F^{\mu\nu}-(D_{\mu}\psi)^{\dagger}(D^{\mu}\psi)-W(|\psi|^{2})\right)\,,$
(63) $\displaystyle
S_{\text{surf}}:=-\frac{1}{8\pi}\int\mathrm{d}^{3}x\sqrt{h_{E}}K\,,$ (65)
$\displaystyle S_{\text{ct}}:=\frac{1}{8\pi}\int\mathrm{d}^{3}x\sqrt{h_{\rm
E}}K_{0}-\frac{\alpha}{4\pi}\int\mathrm{d}^{3}x\sqrt{h_{\rm
E}}n_{\nu}F^{\mu\nu}A_{\nu}\,.$

$S_{\text{surf}}$ and $S_{\text{ct}}$ denote the Gibbons-Hawking term and the
counter terms respectively. Explicitly, $h_{E}$ and $K$ denote induce metric
and the intrinsic curvature on arbitrary boundary $\partial M$. In addition,
the first term of $S_{\text{ct}}$, $K_{0}=\frac{2}{r}$, can remove the
infinity arising from the spherical coordinates, while the second term with
$\alpha$ is introduced for removing the boundary term relevant to the
electromagnetic field in the case of canonical ensemble. In other words,
$\alpha=1,0$ denotes the case of canonical ensemble and grand-canonical
ensemble respectively. In following we will present more illustration about
the $\alpha$ term.

In the grand canonical ensemble, performing variation on the Euclidean action
Eq. (59) with respect to the fields gives

$\displaystyle\delta S=$
$\displaystyle-\frac{1}{16\pi}\int_{M}\mathrm{d}^{4}x\sqrt{g_{E}}\left(E_{g}^{\mu\nu}\delta
g_{\mu\nu}+E_{\psi}\delta\psi\right)$ (66) (68)
$\displaystyle-\frac{1}{4\pi}\int_{M}\mathrm{d}^{4}x\sqrt{g_{E}}E_{A}^{\mu}\delta
A_{\mu}$ (70) $\displaystyle+\frac{1}{4\pi}\int_{\partial
M}\mathrm{d}^{3}x\sqrt{h_{E}}n_{\mu}F^{\mu\nu}\delta A_{\nu}\,,$

where $n^{\mu}$ denotes the unit normal vector orthogonal to the boundary
$\partial M$, $E_{g}^{\mu\nu},E_{\psi},E_{A}^{\mu}$ are the E.O.M given in Eq.
(6). For well-posed variation principle, one must impose an appropriate
boundary condition. In general, $\delta A_{\mu}=0$ on $\partial M$ is imposed,
which is appropriate for grand canonical ensemble with fixing chemical
potential $\mu$. However, in the case of canonical ensemble with fixing the
charge $Q$ on $\partial M$ as the boundary condition, the second term of
$S_{\text{ct}}$ must be under consideration, namely $\alpha=1$. Moreover, the
well-posed variation principle requires the boundary condition
$\delta(n_{\mu}F^{\mu\nu})=0$ Hawking:1995ap ; Caldarelli:1999xj .

In the case that the scalar field is zero, the background is a RN black hole,
of which the metric reads

$\displaystyle\mathrm{d}s_{\text{RN}}^{2}=-f_{\text{RN}}(r)\mathrm{d}t^{2}+\frac{\mathrm{d}r^{2}}{f_{\text{RN}}(r)}+r^{2}\mathrm{d}\Omega^{2}\,,$
(71) (73) $\displaystyle
f_{\text{RN}}(r)=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\,.$

Plugging the RN black hole background Eq. (71) into $S_{\text{E}}$ and setting
$r\to\infty$ as the boundary, we find

$S_{\text{E}_{\text{RN}}}=\frac{\beta}{2}(M\pm Q\Phi_{r_{+}})\,,$ (74)

in which $\pm$ denotes the canonical ensemble case and the grand canonical
ensemble case respectively. Then the Helmholtz free energy $F(T,Q)$ and the
Gibbs free energy $G(\mu,T)$ of RN black hole can directly read

$\displaystyle F_{\text{RN}}(T,Q)=\frac{1}{2}\left(M(T,Q)-\mu(T,Q)Q\right)\,,$
(75) (77) $\displaystyle G_{\text{RN}}=\frac{1}{2}\left(M(\mu,T)+\mu
Q(\mu,T)\right)\,,$

in which the chemical potential $\mu=\frac{Q}{r_{+}}$ is interpreted as the
electric potential at the infinity. In the canonical ensemble, the ADM mass of
the RN black hole $M_{\rm RN}$ can be solved by the following relation

$\displaystyle M_{\rm RN}=\frac{1}{2}\left(r_{h_{\rm RN}}+r_{-}\right)\,,\quad
Q_{\rm RN}^{2}=r_{h_{\rm RN}}r_{-}\,,$ (78) (80) $\displaystyle T_{\rm
RN}=\frac{r_{h_{\rm RN}}-r_{-}}{4\pi r_{h_{\rm
RN}}^{2}}\,,\quad\mu=\frac{Q}{r_{h_{\rm RN}}}\,.$

where $r_{-}$ is the inner horizon of the RN black hole. However, the
analytical expression of $M_{\rm RN}(T_{\rm RN},Q_{\rm RN})$ is too
complicated to present since there are three real or complex roots. In
practice, numerically selecting the real root of $M_{\rm RN}(T_{\rm RN},Q_{\rm
RN})$ with $M_{\rm RN}>\frac{r_{h_{\rm RN}}}{2}$, we then read off $F_{\rm
RN}(T_{\rm RN},Q_{\rm RN})$. As to in grand canonical ensemble, the Eq. (78)
gives

$M_{\rm RN}=\frac{1-\mu_{\rm RN}^{4}}{8\pi T_{\rm RN}}\,,\quad Q_{\rm
RN}=\frac{\mu_{\rm RN}(\mu_{\rm RN}^{2}-1)}{4\pi T_{\rm RN}}\,.$ (81)

We thus read off the value of $G_{\rm RN}(\mu,T_{\rm RN})$ from Eq. (75) as
well.

####  III.2.1 The case of probe limit $c\ll 1$

Now we turn to consider the case of the scalarized charged black hole. We
firstly present a proof that, in the probe limit $c\to 0$, the RN black hole
will be more stable than the scalarized charged black hole in both canonical
ensemble and the grand canonical ensemble. Firstly, rescaling $\psi$, namely
$\tilde{\psi}\to c\psi$, gives

$S_{\text{E}}=S_{\text{E}_{\text{EM}}}+\frac{\beta
c^{2}}{16\pi}\int\mathrm{d}^{4}x\sqrt{g_{\text{E}}}(D^{\mu}\tilde{\psi})^{{\dagger}}(D_{\mu}\tilde{\psi})+m^{2}\log(1+|\tilde{\psi}|^{2})\,,$
(82)

where $S_{\text{E}_{\text{EM}}}$ denotes the Euclidean action contributed by
Einstein-Maxwell theory. At the limit $c\to 0$, it is clear that the solution
will be a RN black hole with the metric $g_{\mu\nu}^{\text{RN}}$ and gauge
potential $A_{\mu}^{\text{RN}}$. When $c\neq 0$ but $c\ll q$, one can treat
the contribution of scalar field as a perturbation of order
$\mathcal{O}(c^{2})$. Let us assume that the metric and gauge potential become

$g_{\mu\nu}=g_{\mu\nu}^{\text{RN}}+c^{2}g_{\mu\nu}^{(1)},~{}~{}A_{\mu}=A_{\mu}^{\text{RN}}+c^{2}A_{\mu}^{(1)},~{}~{}~{}c^{2}\ll
1\,.$ (83)

Upon the $\mathcal{O}(c^{2})$ order, we find the Euclidean action contributed
by Einstein-Maxwell theory reads

$\displaystyle S_{\text{E}_{\text{EM}}}=S_{\text{E}_{\rm
RN}}+c^{2}\int_{M}\sqrt{g_{E}}\mathrm{d}^{4}x\left[(G_{\mu\nu}|_{g_{\mu\nu}=g_{\mu\nu}^{\text{RN}}})g_{\mu\nu}^{(1)}\right.$
(84) (86)
$\displaystyle+\left.(\nabla_{\mu}F^{\mu\nu}|_{A_{\mu}=A_{\mu}^{\text{RN}}})A_{\nu}^{(1)}\right]+\mathcal{O}(c^{4})\,,$

where $S_{\text{E}_{\rm RN}}$ is the on-shell action in RN background given by

$G_{\mu\nu}|_{g_{\mu\nu}=g_{\mu\nu}^{\text{RN}}}=\nabla_{\mu}F^{\mu\nu}|_{A_{\mu}=A_{\mu}^{\text{RN}}}=0\,.$
(87)

Therefore, we find the $S_{\text{E}_{\text{EM}}}=S_{\text{E}_{\rm
RN}}+\mathcal{O}(c^{4})$, which implies upon the leading order of $c$, the on-
shell Euclidean action of Eq. (82) can be regard as

$S_{\text{E}_{\text{on-shell}}}=S_{\text{E}_{\rm
RN}}+S_{\text{E}_{\text{scalar}}}\,$ (88)

with

$S_{\text{E}_{\text{scalar}}}=\frac{\beta
c^{2}}{16\pi}\int\mathrm{d}^{4}x\sqrt{g_{\text{E}}}[(\bar{D}^{\mu}\tilde{\psi})^{{\dagger}}(\bar{D}^{\mu}\tilde{\psi})+m^{2}\log(1+|\tilde{\psi}|^{2})]\,,$
(89)

where $\bar{D}_{\mu}$ denotes the covariant derivative under the RN black hole
background (71).

Plugging the RN black hole background Eq. (71) into Eq. (89), we have

$\displaystyle S_{\text{E}_{\text{scalar}}}=\frac{\beta
c^{2}}{4}\int^{\infty}_{r_{h_{\rm RN}}}\mathrm{d}r~{}r^{2}$ (90) (92)
$\displaystyle\left(f_{\text{RN}}\tilde{\psi^{\prime}}^{2}+\frac{q^{2}\phi_{\rm
RN}^{2}\tilde{\psi}^{2}}{f_{\text{RN}}}+m^{2}\log(1+|\tilde{\psi}|^{2})\right),\,$

where the prime means the derivative with respect to $r$. It needs to note
$\phi_{\rm RN}^{2}=-Q^{2}/r^{2}\leq 0$, since $\phi_{\rm RN}$ denotes the
electric potential of the RN black hole in Euclidean spacetime, namely $A^{\rm
RN}_{\mu}=\phi_{\rm RN}(\mathrm{d}\tau)_{\mu}=\frac{iQ}{r}$ (here $\tau$ is
Euclidean time). This leads that the sign in the integration is undefined.
Note that the equation of motion associated with the scalar field in RN
background

$\bar{D}^{\mu}\bar{D}_{\mu}\tilde{\psi}-\frac{m^{2}}{1+|\tilde{\psi|}^{2}}\psi=0,$
(93)

where explicitly given

$f_{\rm RN}\tilde{\psi}^{\prime\prime}+f_{\rm
RN}^{\prime}\tilde{\psi}^{\prime}+\frac{2f_{\rm
RN}}{r}\tilde{\psi}^{\prime}-\frac{m^{2}}{1+|\tilde{\psi}|^{2}}\tilde{\psi}-\frac{q^{2}\phi_{\rm
RN}^{2}}{f}\tilde{\psi}=0$ (94)

Then the Eq. (90) gives

$\displaystyle S_{\text{E}_{\text{scalar}}}=\frac{\beta
c^{2}}{4}\left(\left(r^{2}\tilde{\psi}\tilde{\psi}^{\prime}f_{\rm
RN}\right)|^{\infty}_{r_{h_{\rm RN}}}\right.$ (95) (97)
$\displaystyle+\left.\int^{\infty}_{r_{h_{\rm
RN}}}r^{2}m^{2}\left(\log(1+|\tilde{\psi}|^{2})-\frac{\tilde{\psi}^{2}}{1+|\tilde{\psi}|^{2}}\right)\mathrm{d}r\right)\,.$

Given appropriate boundary condition $\psi(\infty)=0$ and $f_{\rm
RN}|_{r=r_{h_{\rm RN}}}=0$, the first term of Eq. (95) will vanish. To
proceed, we consider the second term of Eq. (95). One can construct a
auxiliary function $t(x)$ as

$t(x)=\log(1+x)-\frac{x}{1+x}\,,\quad\frac{\mathrm{d}t(x)}{\mathrm{d}x}=\frac{x}{(1+x)^{2}}.$
(98)

It is easy to verify that for $x>0$, $\frac{\mathrm{d}t(x)}{\mathrm{d}x}>0$
always hold, indicating $t(x)$ is a monotone increasing function where the
minimum value is lying on $x=0,t(0)=0$. Since $|\psi|^{2}$ is positive for any
$r$, the integrand of the second term of Eq. (95) is always positive.
Therefore, we have proved

$S_{\rm{E}_{\text{scalar}}}>0.$ (99)

Furthermore, according to $F=\frac{1}{\beta}S_{\text{E}},$ and
$G=\frac{1}{\beta}S_{\text{E}}$, we have

$\displaystyle\Delta F=F-F_{RN}=F_{\text{scalar}}>0\,,$ (100)
$\displaystyle\Delta G=G-F_{RN}=G_{\text{scalar}}>0\,.$ (101)

To conclude, base on our demonstration above, the contribution of the complex
scalar field to the $F(T,Q)$ and $G(T,\mu)$ will be always positive. According
to the stability requirement that the smaller free energy indicates the more
stable of the black hole, we thus prove that in probe limit, the RN black hole
is more stable than the scalarized charged black hole in both canonical
ensemble and grand canonical ensemble.

####  III.2.2 The case of general $c$: Canonical Ensemble

When $c$ is not infinitesimal, the higher order terms of $c$ play role and
above proof is broken. In the following, we consider the Euclidean action Eq.
(59) in general. As there is not analytical solution for scalarized black
hole, we can only compute the free energy numerically. In order to do that,
let us first some useful formulas, which can simplify the numerical
computation of free energy. Performing the Wick rotation, the Euclidean line
element gives

$\mathrm{d}s^{2}=f(r)\mathrm{e}^{-\chi(r)}\mathrm{d}\tau^{2}+\frac{\mathrm{d}r^{2}}{f(r)}+r^{2}\mathrm{d}\Omega^{2}$
(102)

and

$A_{\mu}=i\phi(r)(\mathrm{d}\tau)_{\mu},~{}~{}\psi=\psi(r)\,.$ (103)

Note that the trick given in Hartnoll:2008kx , we also find the following
relation between the on shell Lagrangian and the $\theta\theta$ component of
the energy momentum tensor,

$2T^{\theta}{}_{\theta}={\cal L_{\text{on-shell}}}-R.$ (104)

Consider the equation of motion given in Eq. (6), we arrive

${\cal L}_{\rm on-
shell}=-{G^{t}}_{t}-{G^{r}}_{r}=\frac{2}{r^{2}}[(rf)^{\prime}+r\chi^{\prime}f-1]\,,$
(105)

where $G_{\mu\nu}$ is the Einstein tensor. Recall the Euclidean on-shell
action of RN black hole Eq. (74), we obtain the Euclidean on-shell action of
Eq. (59),

$S_{\text{E}_{\text{on-shell}}}=\frac{\beta
M}{2}+\frac{\beta}{2}\int^{\infty}_{r_{h}}\mathrm{d}r[(rf\mathrm{e}^{-\chi/2})^{\prime}-\mathrm{e}^{-\chi/2}]+\alpha\mu
Q\,.$ (106)

After performing integration by parts, we respectively obtain the Gibbs free
energy in grand canonical ensemble

$G(\mu,T)=\frac{r_{h}}{2}-\frac{M}{2}+\frac{1}{2}\int_{r_{h}}^{\infty}(1-\mathrm{e}^{-\chi/2})\mathrm{d}r\,,$
(107)

and the Helmholtz free energy in canonical ensemble

$F(Q,T)=\frac{r_{h}}{2}-\frac{M}{2}+\frac{1}{2}\int_{r_{h}}^{\infty}(1-\mathrm{e}^{-\chi/2})\mathrm{d}r+\mu
Q\,.$ (108)

Given a scalarized charged black hole $\\{c_{i},\psi_{0_{i}}\\}$, the Hawking
temperature and the chemical potential read

$T=\frac{f^{\prime}(r_{h})\mathrm{e}^{-\chi(r_{h})/2}}{4\pi}\,,\quad\mu=[\phi(r)+r\phi^{\prime}(r)]|_{r\to\infty}.$
(109)

In the case of canonical ensemble, we pick out the temperature $T$, total
charge $Q$ and evaluate the helmholtz free energy $F(T,Q)$ of every scalarized
charged black hole solution $\\{c_{i},\psi_{0_{i}}\\}$ base on Eq. (48), Eq.
(109) and Eq. (108). Using shooting method, we work with dimensionless
$(Q/T,\Delta F/T)$ with fixing temperature $T=T_{0}=0.02858$ and $(T/Q,\Delta
F/Q)$ plane with fixing $Q=Q_{0}=1.318$ without losing generality. (See also
the middle and the right plot in Fig. 5 and Appendix. B.2 for detail
discussion).

Figure 5: Top: A typical example about relationship between $\Delta F$ and
$c$. Here we take $\\{r_{h}=1,m=0.01,q=11m,\phi_{1}=0.8\\}$. Middle: The
relation between $Q/T$ and $\Delta F/T$ by fixing $T=0.02858$. Bottom: The
relation between $T/Q$ and $\Delta F/Q$ with fixing $Q=1.318$.

From the plot in $(Q/T,\Delta F/T)$ and $(T/Q,\Delta F/Q)$ plane we present in
Fig. 5, one can find that there also exist small intervals
$Q/T\in[46.15,49.30]$ and $T/Q\in[0.02167,0.02318]$ in which the $F(T,Q)$ of
sclarized charged black hole is smaller than the corresponding RN black hole
in canonical ensemble. To summarize, this result indicates that in the case of
canonical ensemble, the scalarized charged black hole is more stable than the
corresponding RN black hole in the region $Q/T\in[46.15,49.30]$ and
$T/Q\in[0.02167,0.02318]$. Remarkably, when $Q/T>49.30$ and $T/Q>0.02318$
region there does not exist the corresponding RN black hole.

####  III.2.3 The case of general $c$: Grand Canonical Ensemble

Finally, we turn to consider the grand canonical ensemble case where the
thermodynamic variables is temperature $T$ and chemical potential $\mu$.
Different from in the case of microcanonical ensemble and canonical ensemble,
from Eq. (81) one will find that given a numerical scalarized charged black
hole solution with temperature $T_{i}$ and chemical chemistry $\mu_{i}$, there
is always a corresponding RN black hole sharing the same $T_{i}$ and
$\mu_{i}$. With this in mind, we pick out $\Delta G(T,\mu)=G(T,\mu)-G_{\rm
RN}(T,\mu)$ and present a plot in $(c,\Delta G)$ plane (See also the plots in
Fig. 6) from which one can observe that $\Delta G$ is always positive for any
$c_{i}$.

Figure 6: Top: A typical example about the relationship between $\Delta G$ and
$c$. Here we take $c\in[0.001,0.4]$. Bottom: Zooming in the left plot in the
small $c$ region $c\in[0.001,0.004]$. Here we take
$\\{r_{h}=1,m=0.01,q=11m,\phi_{1}=0.8\\}$.

The absence of appropriate seed solution indicates that the Gibbs free energy
$G(T,\mu)$ of an arbitrary scalarized charged black hole solution might be
always larger than the corresponding RN black hole in grand canonical
ensemble. In following we take a small $c=0.00101$ and larger $c=0.1$ as
example to study the behavior of $\Delta G$ related to the fixing temperature
$T$ and chemical potential $\mu$. In this case, one can read off the $G_{\rm
RN}(T,\mu)$ of the corresponding RN black hole with temperature $T_{i}$ and
chemical potential $\mu_{i}$ from Eq. (75) and Eq. (78). We then work with a
dimensionless $(\mu/T,\Delta G/T)$ plane by fixing $T=0.02864$ and
dimensionless $(T/\mu,\Delta G/\mu)$ plane with fixing $\mu=0.800$
respectively, without losing generality as we mentioned above. We also offer
two plots as well (See also the plots presented in Fig. 7 and Appendix. B.3
for detail).

Figure 7: In the case of $c=0.00101$, Top: A typical example about the
relation between $\mu/T$ and $\Delta G/T$ with fixing $T=0.02864$. Bottom: The
other typical example about the relation between $T/\mu$ and $\Delta G/\mu)$
with fixing $\mu=0.800$.

These results show that the $G(T,\mu)$ of the scalarized charged black hole is
always larger than the corresponding RN black hole in grand canonical ensemble
when $c$ is small. Therefore it can be concluded that the RN black hole is
more stable than scalarized charged black hole in grand canonical ensemble,
which is consistent with our proof in probe limit.

As to Working with dimensionless $(\mu/T,\Delta G/T)$ plane with fixing
$T=0.02721$ and dimensionless $(T/\mu,\Delta G/\mu)$ plane with fixing
$\mu=0.8114$ respectively, we present a plot (See the plots of Fig. 8)

Figure 8: In the case of $c=0.1$, Top: A typical figure about the relation
between of $\mu$ and $\Delta G$ in dimensionless $(\mu/T,\Delta G/T)$ plane
with fixing $T=0.02721$. Bottom: Another typical example about the relation
between of $T$ and $\Delta G$ in dimensionless $(T/\mu,\Delta G/\mu)$ plane
with fixing $\mu=0.8114$.

and observe that $\Delta G$ is also positive for any $T$ and $\mu$ from the
plots, it can be concluded that for a larger $c$, the RN black hole is also
more stable than the scalarized charged black hole in grand canonical
ensemble. We have check carefully for other different parameters and find the
same conclusion. Hence, we summarize that in thermodynamics, the RN black hole
is more stable than the scalarized charge black hole in grand canonical
ensemble.

##  IV Kinetic Stability

In this section, we turn to study the kinetic stability of the scalarized
charged black hole solution, i.e. the stability against a small perturbation.
At first, based on perturbation theory of black hole, we setup a general model
to obtain one-dimension, Schrödinger like radial function in frequency domain.
We then investigate the validity of the shooting method numerically through
analyzing numerical error, pointing out that the shooting method can catch the
unstable modes efficiently, rather than stable modes. Based on the argument in
Konoplya:2019hlu , we adopt the WKB approximation method to calculate the
stable mode of perturbative field since the damping mode can be given by means
of the WKB approximation method.

Practically, considering a massless real perturbative scalar field, we obtain
the linearized equation of motion associated with the scalar perturbation. Due
to the static and spherical black hole background, the equation of motion of
the scalar perturbation can be reduced to 1 dimension in frequency
domainChandrasekhar:1975zza . Imposing an appropriate physical boundary
condition, ingoing wave near the horizon and outgoing wave at the spatial
infinity, resonance state of the scalar perturbation arise, picking out a
class of complex frequency $\omega=\omega_{\rm R}+i\omega_{\rm I}$ so called
black hole quasi-normal modes (QNMs). Moreover, the imaginary part of
frequency $\omega_{\rm I}$ indicates energy dissipation at both the horizon
and the spactial infinity. If $\omega_{\rm I}>0$, the scalar perturbation will
grow exponentially, leading to the instability of the black hole background at
least at linear level. If $\omega_{\rm I}<0$, the scalar field will not
trigger on the unstable of spacetime for exponentially damp, and finally
dissipative out rapidly.

Firstly, we begin with reducing the master equation of the scalar perturbation
in radial equation in frequency domain. Given a spherical line element Eq.
(16), the equation of motion of the perturbative scalar field gives

$\square\psi_{2}=0.$ (110)

where $\square=\nabla^{\mu}\nabla_{\mu}$ is the d’Alembert operator in
spherical line element background, Eq. (16). By adopting the separation
variables method, we take the following anstaz of $\psi_{2}$ under the
spherical background,

$\psi_{2}=\mathrm{e}^{-i\omega t}R(r)Y_{lm}(\theta,\varphi),$ (111)

where $l$ and $m$ are the azimuthal quantum number and magnetic quantum number
respectively, satisfying $l>m,m=0,1,2,...$. The radial equation reads

$\Delta\frac{d}{dr}\Delta\frac{dR}{dr}+UR=0,$ (112)

where

$\Delta=r^{2}f(r)\mathrm{e}^{-\frac{\chi}{2}}\,,\quad
U=\frac{\omega^{2}r^{2}}{f\mathrm{e}^{-\frac{\chi}{2}}}-l(l+1).$ (113)

By introducing the tortoise coordinate and a new radial equation $\tilde{R}$
as

$\frac{dy}{dr}=\frac{r^{2}}{\Delta}\,,\quad R=r\tilde{R}\,,$ (114)

the radial equation Eq.(112) can be written in the following standard wave
function form

$\frac{d^{2}\tilde{R}}{dy^{2}}+\tilde{U}\tilde{R}=0\,,$ (115)

where

$\displaystyle\tilde{U}$ $\displaystyle=$ $\displaystyle\omega^{2}-V$ (116)
(117) $\displaystyle=$
$\displaystyle\omega^{2}-\frac{l(l+1)}{r(y)^{2}}f(y)\mathrm{e}^{-\chi(y)}$
(118) (120)
$\displaystyle-\frac{1}{2r(y)}\frac{d}{dr}\left(f(y)^{2}\mathrm{e}^{-\chi(y)}\right)\,.$

In Eq. (116), we denotes $f(y),\chi(y)$ and $\chi(y)$ as a function with
respect to $y$. From the definition of the tortoise coordinate, one can read
$y\to-\infty$ corresponding to $r\to r_{h}$, while $y\to\infty$ corresponding
to $r\to\infty$. To single out a series of QNMs, we impose the following
boundary condition,

$\tilde{R}\varpropto\left\\{\begin{aligned} &\mathrm{e}^{-i\omega
y}\quad\quad&&y\to-\infty\\\ &\mathrm{e}^{i\omega
y}\quad\quad&&y\to\infty\,,\end{aligned}\right.$ (121)

indicating that the probing scalar is pure going wave at the horizon and pure
outgoing wave at the spatial infinity. In general, in order to investigate
both the stable and unstable mode of perturbative field, one can solve QNMs
using numerical method, popularly the shooting method. There are also serval
effective approximative methods to calculate QNMs, for example, the WKB
approximation method. It has been argued that the WKB approximation can only
catch the stable modes with $\omega_{\rm I}<0$ even if the spectrum contains
unstable modes Konoplya:2019hlu . We will show that, on the contrary, the
shooting method can only catch the unstable modes of $\omega_{\rm I}>0$ but
cannot catch the stable modes. Thus, the combination of shooting method and
the WKB approximation method can offer us a complete analyses on QNMs.

###  IV.1 Shooting Method and Analysis on Numerical Error

In this section, we analyze numerical error of shooting method and illustrate
its validity on the calculation of the unstable modes. Let us first briefly
explain how to use shooting method to find QNMs.

To setup the numerical procedure, we firstly interpret Eq. (114) and Eq. (115)
as a set of ODE (one first-order differential equation as well as one second
order differential equation), one thus need three boundary condition. In
addition to Eq.(121), we need the last boundary condition associated to the
radial coordinate $r$ pand tortoise coordinate $y$. we thus asymptotically
expand near horizon. Specifically, in near horizon region, Eq. (114) gives

$\frac{\mathrm{d}y}{\mathrm{d}r}=\frac{1}{f_{1}(r-r_{h})}\,,$ (122)

where $f_{1}$ is given by Eq. (37) and the equation of motion Eq. (19).
Interpreted the radial coordinate as a analytical function with respect to the
tortoise coordinate $y$. The Eq. (122) then can be solved that

$r=r_{h}+\mathrm{e}^{f_{1}y}\,,\quad y\to-\infty\,.$ (123)

Therefore, we obtain a numerical solvable boundary-valued question with two
differential equations

$\left\\{\begin{aligned}
&\frac{d^{2}\tilde{R}}{dy^{2}}+\tilde{U}\tilde{R}=0\,,\cr&\frac{dr}{dy}=f(r(y))\,,\end{aligned}\right.$
(124)

where $\tilde{U}$ has been given in Eq. (116). To obtain a numerical
eigenfunction of QNMs satisfying the boundary condition Eq. (121), working
with complex numerics, we primarily numerically perform integration on
Eq.(124) with Eq. (121) and

$\tilde{R}\varpropto\mathrm{e}^{-i\omega y}\,,\quad y\to\infty.$ (125)

In general, the numerical integration will give two branch solution in
$y\to\infty$ region,

$\tilde{R}\varpropto\mathrm{e}^{i\omega y}+{B}\mathrm{e}^{-i\omega y}\,,\quad
y\to\infty.$ (126)

Furthermore, we have

${B}\varpropto\frac{\mathrm{e}^{i\omega
y}}{2\omega}\left(R\omega+iR^{\prime}(y)\right)\,,\quad y\to\infty.$ (127)

The QNMs will then give

$B\left(\omega_{\text{QNMs}}\right)=0\,.$ (128)

In practice, one can locate the QNMs by drawing a density figure related to
$\frac{1}{B}$ in $\left(\omega_{\text{R}},\omega_{\text{I}}\right)$ plane.
Since the QNMs will give a vanish $R_{-}$, the location of QNMs in density
figure can thus be observed as a bright spot in a density figure. Through
observe the approximative region of the bright spot related to
$\omega_{\text{R}}$ and $\omega_{\text{I}}$, denoted as
$\omega_{\text{Initial}}=\omega_{\text{R}_{\text{Initial}}}+i\omega_{\text{I}_{\text{Initial}}}$.
One can solve the QNMs by standard shooting method.

Numerically, the error includes three parts in general. One comes from the
finite difference when we solve the differential equation numerically, which
can be suppressed by using higher order methods or smaller step-size. The
second one comes from the fact that we have to set two finite cut-off at the
horizon and infinite boundary. The third part comes from the float-point error
of computer. We will show in following that, if imaginary part of QNMs is
positive, errors of the second and third parts lead to the computational
complexity will increase exponentially if we improve the desired accuracy of
QNMs.

Instead of considering the boundary condition (121), we firstly consider a
general boundary-valued eigenvalue problem with following boundary condition
for the radial equation Eq. (115)

$\tilde{R}\varpropto\left\\{\begin{aligned} &\mathrm{e}^{-i\omega
y}+A\mathrm{e}^{i\omega y}\quad\quad&&y\to-\infty\\\ &\mathrm{e}^{i\omega
y}+B\mathrm{e}^{-i\omega y}\quad\quad&&y\to\infty\,,\end{aligned}\right.$
(129)

where $A,B$ are two complex constants. For a pair of given specific $A$ and
$B$, the radial equation Eq. (115) can numerically give a complex frequency
$\omega=\omega_{\text{R}}+i\omega_{\text{I}}$. Therefore, we can interpret a
series of complex frequency as an analytical function with respect to $A$ and
$B$, denoted as $\omega(A,B)$. Theoretically, it is obvious that the QNMs
arise from vanishing both $A$ and $B$, $\omega_{\text{QNM}}=\omega(0,0)$,
matching the boundary condition Eq. (121). However, in the shooting method in
numerics, $A,B$ cannot exactly vanish due to floating-point error.
Practically, in order to solve QNMs in numerics, the real boundary condition
in the numerical computations is

$\tilde{R}\varpropto\left\\{\begin{aligned} &\mathrm{e}^{-i\omega
y}+\epsilon\quad\quad&&y\to-\infty\\\ &\mathrm{e}^{i\omega
y}+\epsilon\quad\quad&&y\to\infty\,,\end{aligned}\right.$ (130)

where $\epsilon$ denote the floating-point error in numerics. For a given
complex frequency, comparing with Eq. (129) and Eq. (130), one can find that
it is equivalent to set a pair of nonzero $A$ and $B$ and

$A\sim\epsilon\mathrm{e}^{\omega_{\text{I}}y_{-}}\,,\quad
B\sim\epsilon\mathrm{e}^{-\omega_{\text{I}}y_{+}}\,,$ (131)

where $y_{\pm}$ are numerical cutoff at both negative infinity and positive
infinity respectively. Thus, for a computer which has fixed machine precision,
the floating-point error restricts our ability to find the QNMs in arbitrary
precision. To improve the accuracy of finding QNMs, we have to increase the
machine precision of computer. In general, $\omega(A,B)$ can be expanded as
following approximation upon the 1st order of $A$ and $B$

$\omega(A,B)=\omega_{\text{QNM}}+\frac{\partial\omega}{\partial
A}A+\frac{\partial\omega}{\partial B}B+\cdots\,.$ (132)

We thus can obtain following estimation on the error of finding QNMs

$\omega_{\text{err}}:=|\omega(A,B)-\omega_{\text{QNM}}|\sim\max\left\\{\left|\frac{\partial\omega}{\partial
A}\right|\left|A\right|,\left|\frac{\partial\omega}{\partial
B}\right|\left|B\right|\right\\}.$ (133)

denotes the numerical error associated with the QNMs. As it is reasonable to
assume that $\omega(A,B)$ is the analytical function of $A$ and $B$, so
$\partial\omega/\partial A$ and $\partial\omega/\partial B$ are both finite.
Thus,it can be concluded that the error is controlled by $A$ and $B$, i.e.
$\epsilon\mathrm{e}^{\omega_{\text{I}}y_{-}}$ and
$\epsilon\mathrm{e}^{-\omega_{\text{I}}y_{+}}$.

For a series of unstable mode, $\omega_{\text{I}}>0$, triggering on
exponential growth as time evolution, from Eq.(131) one can find that an
effective working precision with given $\epsilon$ both $A,B$ will
exponentially decay as large cutoff, ensuing numerical accuracy that
$\omega_{\text{err}}$ is convergence in finite computation time. Therefore,
the error caused by float-point is suppressed exponentially and we only need
to care about the errors caused by finite difference and cut-off. In this
case, the shooting method can catch the unstable mode efficiently.

However, it does not work for stable modes $\omega_{\text{I}}<0$, in which
$\psi_{2}$ will exponential damp as time evolution. In this case, we will
suffer a contradiction that if we suppose $\omega_{\text{err}}$ is convergence
to guarantee accuracy, the required working precision, related to computing
time, will grow exponentially as the growth of cutoffs $y_{\pm}$, while in
order to obtain a precision-guarantee QNMs, we need to set $y_{\pm}$ as large
as possible to reduce cutoff error. Said another specific way, if one improve
ten times of the cutoff, $y_{\pm}\to 10y_{\pm}$, the required working
precision has to be improved up to $\mathrm{e}^{10\omega_{\text{I}}}$ to
ensure computation precision. Consequently, the numerical approach cannot
catch the stable modes efficiently.

To illustrate our analysis more clearly, we take the negative Pöschl-Teller
potential as an example in Appendix. A, in which shows how the shooting method
can catch unstable modes efficiently. Based on our illustration above, we
adopt the numerical approach to calculate the unstable mode
$\omega_{\text{I}}>0$ while the WKB approximation method to calculate the
stable mode $\omega_{\text{I}}<0$.

To end this section, we adopt the numerical solution, denoted as
$\\{c=0.1,\psi_{0}\approx 0.1988\\}$ with
$\\{r_{h}=1,\phi_{1}=0.8,m=\frac{1}{100},q=\frac{11}{10}m\\}$, to investigate
whether there exist unstable mode within QNMs. Following the above procedure,
we show a 3D-plot associated with $1/B$ with $(\omega_{\rm R},\omega_{\rm I})$
plane in Fig. 9. One can easily observe that there does not exist any peak of
singularity which indicates that the existence of unstable mode within QNMs.
Furthermore, we adopt some another numerical solutions
$\\{c_{i},\psi_{0_{i}}\\}$ to investigate whether there exist unstable mode
and cannot observe any other unstable mode within QNMs. Before we give a
conclusion that the scalarized charged black hole is stable against a neutral
scalar field at linear level, in the following section we also adopt the WKB
approximation method to calculate the stable mode in QNMs.

Figure 9: 3D-plot associated with $\frac{1}{B}$ with $(\omega_{\rm
R},\omega_{\rm I})$ plane , taking the numerical solution
$\\{c=0.1,\psi_{0}\approx 0.1988\\}$ with
$\\{r_{h}=1,\phi_{1}=0.8,m=\frac{1}{100},q=\frac{11}{10}m\\}$ as an example.

###  IV.2 The WKB Approximation Method

At first, we give brief introduction on calculating the QNMs by the WKB
approximation method. The WKB approximation method was firstly used to
calculate the QNMs of Schwardzchild black hole by Schutz and Will Shusz . It
then was developed to the 3rd order by Iyer and Will Iyer:1986np ; Iyer:1986nq
and further 6th order by R.A. Konoplya Konoplya:2003ii . Recently it has been
extended to the 12nd order by Matyjasek and Opala Matyjasek:2017psv . In this
paper, we adopt the the 3rd order WKB approximation method to analyze the QNMs
of the scalarized charged black hole. We point out again the we have used
numerical method in the last subsection verified that there is no unstable
QNMs.

Recall the radial function Eq. (115), the effective potential $V(y)$ gives

$V(y)=\frac{l(l+1)}{r(y)^{2}}f(y)\mathrm{e}^{-\chi(y)}+\frac{1}{2r(y)}\frac{d}{dr}\left(f(y)^{2}\mathrm{e}^{-\chi(y)}\right)$
(134)

We firstly consider the numerical scalarized charged black hole solution
$\\{c=0.1,\psi_{0}=0.1987\\}$ as an example, found in Sec. III. Substituting
the numerical solution in Eq. (134), we present a plot of $V(y)$ taking
$l=0,1,2$ as example in Fig. 10. Furthermore, one can numerically read off the
local maximal value $V_{0}(y_{0})$ from Eq. (134) (See also Table. 1)

| $l=0$ | $l=1$ | $l=2$
---|---|---|---
$V_{0}$ | $0.04641$ | 0.2240 | $0.5794$
$y_{0}$ | $-4.485$ | $-4.355$ | $-4.327$

.

Table 1: The local maximal $V_{0}$ of the effective potential $V(y)$ under the
numerical scalarized charged black hole solution
$\\{c=0.1,\psi_{0}=0.1978\\}$.  Figure 10:  The effective potential under the
tortoise coordinate $V(y)$.

We then proceed to compute the QNMs using the WKB approximation method up to
the 3rd order. To maintain accuracy, we respectively take $l=0$ up to $n=2$ as
well as $l=1,2$ up to $n=5$ as example. Recall the formula of the 3rd order
WKB Approximation method, the QNMs have been given explicitly in Ref.
Iyer:1986np ; Iyer:1986nq  as follow,

$\omega^{2}=(V_{0}+\sqrt{-2V^{\prime\prime}_{0}}\Lambda_{2})-i(n+\frac{1}{2})\frac{1}{\sqrt{-2V_{0}^{\prime\prime}}}(1+\Lambda_{3}),$
(135)

where

$\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{-2V_{0}^{\prime\prime}}}\left(\frac{1}{8}\left(\frac{V^{(4)}_{0}}{V_{0}^{\prime\prime}}\right)\left(\alpha^{2}+\frac{1}{4}\right)\right.$
$\displaystyle\left.-\frac{1}{288}\left(\frac{V^{(3)}_{0}}{V_{0}^{\prime\prime}}^{2}(60\alpha^{2}+7)\right)\right),$
(137) $\displaystyle\Lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{-2V_{0}^{\prime\prime}}\left(\frac{5}{6912}\left(\frac{V^{(3)}_{0}}{V_{0}^{\prime\prime}}\right)^{4}(188\alpha^{2}+77)\right.$
(142)
$\displaystyle-\frac{1}{384}\left(\frac{V^{(3)2}_{0}V^{(4)}_{0}}{V^{\prime\prime
3}_{0}}\right)(100\alpha^{2}+51)$
$\displaystyle+\frac{1}{2304}\left(\frac{V^{(4)_{0}}}{V^{\prime\prime}_{0}}\right)^{2}(68\alpha^{2}+67)$
$\displaystyle+\frac{1}{288}\left(\frac{V^{\prime\prime\prime}_{0}V^{(5)}_{0}}{V^{\prime\prime
2}_{0}}\right)(28\alpha^{2}+19)$
$\displaystyle-\left.\frac{1}{288}\left(\frac{V^{(6)}_{0}}{V^{\prime\prime}_{0}}\right)(4\alpha^{2}+5)\right)\,,$

where $V_{0}^{(n)}$ denotes the $n$-th derivative of $V(y)$ with respect to
$y$ located on $y_{0}$ and $\alpha=(n+\frac{1}{2})$. In Eq. (IV.2),
$\Lambda_{2}$ and $\Lambda_{3}$ denote the second order and the third order
approximation associated with the WKB method, respectively Plugging the
specific numerics into Eq. (IV.2), we then show the results of $n\leq l+1$ in
which $l=0,1,2$ respectively for maintaining precision (See also the Table.
2). Since the WKB method is sufficient in high-lying mode but may not be
sufficient enough in low-lying mode, we do not consider the high overtone
case. In Table. 2, we show the numerical results of QNMs up to 3rd and the
associated relative errors $\Delta$, defined as
$\Delta=|\omega-\omega_{2}|/|\omega|$ where $\omega_{2}$ denotes the numerical
QNMs given by the second order WKB method. One can observe that the imaginary
part of the complex frequency $\omega_{\rm I}$ is always negative, indicating
that the numerical solution $\\{c_{0},\psi_{0_{0}}\\}$ is stable against the
scalar perturbation.

$n$ | $l=0$ | $l=1$ | $l=2$
---|---|---|---
$0$ |  | $0.146-0.125i$
---
$~{}24.8\%$
| $0.450-0.112i$
---
$2.03\%$
| $0.748-0.112i$
---
$0.488\%$
$1$ |  | $0.104-0.409i$
---
$~{}22\%$
| $0.419-0.350i$
---
$~{}6.24\%$
| $0.728-0.340i$
---
$1.90\%$
$2$ |  |  | $0.371-0.601i$
---
$~{}10\%$
| $0.692-0.575i$
---
$4.23\%$
$3$ |  |  |  | $0.646-0.816i$
---
7.07 %
Table 2: The QNMs of the scalarized charged black hole with
$\\{c=0.1,\psi_{0}=0.1988\\}$. In each cell we show the numerical results of
QNMs and associated relative error $\Delta=|\omega-\omega_{2}|/|\omega|$.

Recall a series numerical scalarized charged black hole solution, denoted by
$\\{c_{i},\psi_{0_{i}}\\}$ in Sec. III.1, for without lost generality, we
calculate $\omega_{\rm I}$ associated with $c_{i}$ by the WKB approximation
method as well. We find that $\omega_{\rm I}$ are always negative, and do not
change obviously as the growth of $c_{i}$, implying that the QNMs of
scalarized charged black hole is less affected by the amplitude of the non-
linear potential. Furthermore, we also have checked other different parameters
and found the similar results. This suggests that the scalarized charged black
hole should be kinetically stable under a neutral perturbation.111Strictly
speaking, the QNMs do not form a complete bases in mathematics, so above
analysis does not cover all possible perturbations in mathematics.

##  V Conclusion

In this paper, we consider the Einstein-Maxwell theory minimally coupled with
a non-linear complex field. Considering an appropriate boundary condition for
scalarized black hole in asymptotic flat spacetime. We at first briefly list
our main result:

  * •

For general non-linear semi-definite potential, we prove that the
scalarization cannot result from a continuous phase transition.

  * •

Treating the scalarized black hole as a thermodynamical system, we observe
that the discontinuous scalarization on RN black hole will not happen in grand
canonical ensemble but will happen in both microcanonical ensemble and
canonical ensemble.

  * •

We also find that neutral scalar perturbation will not trigger kinetic
instability associated with the scalarized charged black hole by means of
analysing the QNMs using numerical method and the WKB approximation method.

  * •

As a by-product, we use numerical results to give negative answer to Penrose-
Gibbons conjecture and suggest two new versions of Penrose inequality in
charged case.

In detail, motivated by the numerical solution as a counterexample for the no-
hair theorem given in Hong:2020miv , we investigate the thermodynamic
stability of the scalarized charged black hole, compared with the RN black
hole in various ensembles. It needs to note that it is more suitable to choose
grand canonical ensemble for a black hole in astrophysics. However, in this
paper, we still take microcanonical ensemble and canonical ensemble into
account as theoretical research interests.

In microcanonical ensemble, ADM mass $M$ and total charge $Q$ are fixed and
the phase transition will happen towards the direction which increases the
entropy, namely the radius of the event horizon. Giving specific $M$ and $Q$
and working with $(Q/M,\Delta r_{h}/M)$ plane, we find that it is possible the
scalarized charged black hole is more stable than the RN black hole in
thermodynamics. Particularly, our numerical results imply that the scalarized
black hole is more stable than RN black hole when temperature is low enough
and the near extremal RN black hole will always transit into scalarized black
hole via a first order phase transition.

As to in the canonical ensemble and grand canonical ensemble, the stability in
thermodynamics requires the minimal of the Helmholtz free energy $F(T,Q)$ with
fixing Hawking temperature $T$ and total charge $Q$, and the Gibbs free energy
$G(T,\mu)$ with fixing $T$ and chemical potential $\mu$ respectively. Taking
probe limit, we firstly present a proof that the RN black hole is more
thermodynamically stable in both canonical ensemble and grand canonical
ensemble. Following the similar procedure in microcanonical ensemble, giving
specific $T$ and $Q$ and working with $(T/Q,\Delta F/Q)$ and $(Q/T,\Delta
F/T)$ plane, we also find it is possible that the scalarized charged black
hole is more stable in thermodynamics in canonical ensemble and so the RN
black hole may spontaneously scalarize via a first order phase transition in
canonical ensemble. However, in grand canonical ensemble, we find that RN
black hole always have smaller Gibbs free energy and so is more stable than
scalarized black hole, which implies that the RN black hole will not
spontaneously scalarize in grand canonical ensemble.

Finally, we study the kinetic stability of the scalarized charged black hole
against scalar perturbation. Due to the static and spherical spacetime
background, we firstly reduce the master equation of the perturbative scalar
to $1$ dimension in frequency domain. Given pure ingoing wave condition at the
horizon as well as pure outgoing wave condition at the spatial infinity, a
series of complex frequency $\omega=\omega_{\rm R}+i\omega_{\rm I}$ is picked
out, namely the quasi-normal modes (QNMs) of black hole. Through numerical
error analysis, we claim that the shooting method in numerics can efficiently
catch the unstable mode $\omega_{\rm I}>0$, rather than the stable mode
$\omega_{\rm I}<0$. Therefore, we calculate the unstable mode within QNMs by
the shooting method in numerics and conclude that there does not exist any
unstable mode within the QNMs associated with the scalarized charged black
hole. In addition, since the WKB approximation method can effectively catch
the stable modes, rather than the unstable mode, within the QNMs, we adopt the
3rd order WKB approximation method to compute the QNMs of the scalarized
charged black hole and find that the imaginary part of the QNMs $\omega_{\rm
I}$ is always negative. Therefore, we concluded that the scalarized charged
black hole should be kinetically stable under a neutral perturbation. As
further discussion on future work, it is worth to investigate whether other
perturbation, for instance, vector perturbation or tensor perturbation will
trigger on kinetic instability. We thus will keep focusing on this topic in
our future work.

As a by-product of our numerical construction of scalarized black hole, we
also definitely give negative answer to a long-standing conjecture named
Penrose-Gibbons conjecture. Particularly, we find that the the total charge
can be larger than the ADM mass but the temperature is still positive in
scalarized black hole. Based on our numerical results, we propose two new
generalizations of Penrose inequality in charged case and numerically verify
their correctness in our model. Moreover, we will give more discussion in our
future works.

###### Acknowledgements.

We are grateful to Hong Lü, De-Cheng Zou and A. Zhidenko for useful
discussions and Wen-Di Tan, Shi-Fa Guo and Ze Li for proofreading.

##  Appendix A Efficiency of Numerical Method for Searching Unstable Modes:
Negative Pöschl-Teller Potential

In this section, we shall illustrate that the numerical method is efficient
for searching unstable QNMs Negative Pöschl-Teller Potential as example.
Following the procedure given in Ref. Ferrari:1984zz , the "QNMs"(the bound
frequency) of the negative PT potential can be analytically solved and we then
turn to verify it by the shooting method in numerics in this section. At
first, we consider the equation of motion associated with a one-dimension wave
under a Pöschl-Teller potential well (PT potential) in frequency domain can be
written as

$\frac{\mathrm{d}^{2}\Psi}{\mathrm{d}x^{2}}+(\omega^{2}-V_{\text{PT}}(x))\psi(x)=0\,,$
(143)

where

$V_{\text{PT}}=\frac{V_{0}}{\cosh^{2}\alpha x}.$ (144)

In general, the case with $V_{0}>0$, in which the unstable mode
$\Im(\omega)>0$ does not exist, has been under consideration in Ferrari:1984zz
. To stimulate the unstable mode, in this paper we consider the so-called
negative PT potential case where $V_{0}<0$, $V_{0}=V(0)$ and
$\alpha=\left.-\frac{1}{2V_{0}}\left(\frac{\mathrm{d}^{2}U}{\mathrm{d}x^{2}}\right)\right|_{x=0}$
with shape as potential well(See also Fig. 11).

Figure 11: The Pöschl-Teller potential well $V(x)=V_{0}/\cosh^{2}\alpha x$.

The "QNMs" then are given as

$\omega=-i\alpha\left(n+\frac{1}{2}\right)\pm
i\alpha\sqrt{\frac{1}{4}-\frac{V_{0}}{\alpha^{2}}}\,,\quad V_{0}<0\,,$ (145)

where $n=0,1,2,\cdots$. One can find that all $\omega$ lie on the imaginary
axis, where $\omega_{\text{I}}<0$ stands for a stable mode and
$\omega_{\text{I}}>0$ denotes an unstable mode. For convenience, we
specifically taking parameters that $\alpha=1,V_{0}=-\frac{3}{4}$, a series of
frequency is shown in Fig. 12. It is obvious to find that when $n\geq 0$,
there only exist a single unstable mode lying on $\omega=0+i\frac{1}{2}$. We
show our result in Fig. 12.

$n$ | $\omega_{\text{I}}$
---|---
$0$ | $1/2$
$1$ | $-1/2$
$2$ | $-3/2$

Figure 12: The "QNMs" of negative PT potential

In following, we verify the consistency of our illustration in Sec. IV.1 by
the shooting method. In the case of the negative PT potential, the
superfluousness of tortoise coordinate indicates the absence of the
differential equation associated with the definition of tortoise coordinate.
Therefore, we establish a boundary valued problem including a second-order
differential equation Eq. (143) and boundary condition related to QNMs

$\Psi\varpropto\left\\{\begin{aligned} &\mathrm{e}^{-i\omega
x}\quad\quad&&x\to-\infty\\\ &\mathrm{e}^{i\omega
x}\quad\quad&&x\to\infty\,.\end{aligned}\right.$ (146)

In practise, one can numerically integrate Eq. (143) with one-side boundary
condition in $x\to-\infty$, then in $x\to\infty$ region the numerical
integration will generally give

$\Psi\varpropto\Psi_{+}\mathrm{e}^{i\omega x}+\Psi_{-}\mathrm{e}^{-i\omega
x}\,,\quad x\to\infty\,.$ (147)

Furthermore, we have

$\Psi_{-}=\frac{\mathrm{e}^{i\omega
x}}{2\omega}\left(\omega\Psi+i\Psi^{\prime}\right)\,,\quad x\to\infty\,.$
(148)

Comparing the boundary condition associated with the QNMs and interpreting
$\Psi_{-}$ as an analytical function with respect to $\omega$, one can easily
find that vanishing $\Psi_{-}$ will give QNMs, namely

$\Psi_{-}(\omega_{\rm QNMs})=0\,.$ (149)

To calculate the QNMs of negative PT potential, we firstly draw a density plot
of $\frac{1}{\Psi_{-}}$ in $(\omega_{\rm R},\omega_{\rm I})$. As our
illustration in Sec. IV.1, the location of the QNMs will become a bright spot
as a singularity of our density plot related to $1/\Psi_{-}$ since
$\omega_{\rm QNMs}$ indicate $\Psi_{-}=0$. From the density plot one can read
an approximate value of $\omega$ as an initial value, namely
$\omega=\omega_{\rm R_{\rm Initial}}+i\omega_{\rm I_{\rm Initial}}$, then the
QNMs can be solved by standard shooting method

$\\{\Psi_{-}=0\\}\xrightarrow[\omega_{\text{Initial}}=\omega_{\text{R}_{\text{Initial}}}+i\omega_{\text{I}_{\text{Initial}}}]{\text{Shooting}}\\{\omega_{\text{QNMs}}\\}.$
(150)

Choosing $V_{0}=\frac{3}{4}$ and $\alpha=1$. We show a density plot associated
to $1/\Psi_{-}$ in Fig. 13, working in $(\omega_{\text{R}},\omega_{\rm I})$
with $\omega_{\rm R}\in[-0.5,0.5]$ as well as $\omega_{\rm I}\in[0,1]$.

One can easily observe that there exist a bright spot near $\omega_{\rm R}=0$
and $\omega_{\rm I}\approx 0.5$. Using shooting method, we obtain a numerical
result of QNMs that

$\omega_{\rm QNMs}\approx 5.189\times 10^{-10}+i(0.5+5\times 10^{-5})\,,$
(151)

indicating that within the allowable region in numerical error, the numerical
result of the unstable mode is consistent with the value given by theoretical
calculation. Therefore, it is concluded that the shooting method in numerics
can efficiently catch the unstable mode within the QNMs.

Figure 13: The density plot associated with $1/\Psi_{-}$ of negative PT
potential in numerics.

##  Appendix B Shooting Method for Thermodynamical Stability of Scalarized
charged black hole

In this section, we shall present our explicit numerical results of shooting
method for studying the thermodynamical stability associated with the
scalarized charged black hole.

###  B.1 Microcanonical Ensemble

By fixing the independent parameters as
$\\{r_{h}=1,\phi_{1}=0.8,m=0.01,q=\frac{11}{10}m,r_{\infty}=1000\\}$, we set
the $\\{{c_{0}},{\psi_{0_{0}}}\\}$, shown in Fig. 2, as seed solution, then
find next solution $\\{{c_{\pm 1}},\psi_{0_{\pm 1}}\\}$ taking
$\\{{c_{0}}=0.1,{\psi_{0_{0}}}\approx 0.1988\\}$ as the initial value in the
shooting method. For maintaining numerical stability we set the step size
$c_{\pm i}=c_{\pm(i+1)}\pm\frac{1}{1000}$. Starting with
$\\{c_{0}=0.1,\psi_{0}\approx 0.1988\\}$, we then obtain a class of numerical
scalarized charged black hole represented by $\\{c_{\pm i},\psi_{0_{\pm
i}}\\}$ where $i=1,2,3,\cdots$. Given specific $M_{i}$ and $Q_{i}$, the map
Eq. (47) indicates one can leave three parameters free and fix two under the
constrains $\\{M=M_{i},Q=Q_{i},\psi_{+}=0\\}$ (The last one is from the
constrain Eq. (25) for spontaneous scalarization). Said another way, we can
obtain a class of numerical scalarized charged black hole solution with
specific $M_{i}$ and $Q_{i}$ by the shooting method. Concretely, fixing
$\\{c\approx 0.01220,q=1.1m\\}$, we further obtain a series of numerical
charged black hole solutions with scalar hair, denoted by three parameters
$\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}$, for given
$\\{M=M_{i}=0.9676,Q=Q_{i},\psi_{+}=0\\}$. Following the trick we set up
above, the parameters of previous numerical solution
$\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}$ will be used as the initial
value for shooting the next solution
$\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}$ under
$\\{M=M_{i}=0.9676,Q=Q_{i},\psi_{+}=0\\}$. It can be written as

$\displaystyle\\{M=M_{i}=0.9676,Q=Q_{i},\psi_{+}=0\\}\xrightarrow[\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}]{\text{shooting}}$
(152) $\displaystyle\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}\,.$ (153)

Starting with the seed solution
$\\{r_{h_{0}}=1,\phi_{1_{0}}=0.8,\psi_{0_{0}}\approx 0.02473\\}$ together with
fixing parameters $\\{c=0.01220,\psi_{0}=0.02473\\}$, we eventually obtain a
series of numerical scalarized charged black hole solution with
$M=0.9676,Q=Q_{i}$, where we set the step size
$Q_{i+1}=Q_{i}-\frac{1}{1000},Q_{0}=0.9676$ for maintaining accuracy.

###  B.2 Canonical Ensemble

To search a good seed solution, we firstly work with $(c,\Delta F)$ plane.
From the plot presented in the first panel of Fig. 5, we find that there does
exist a zero point of $\Delta F$ in interval $[0.020,0.023]$. Adopting the
interpolation as well as the shooting method, we search another numerical seed
solution $\\{c\approx 0.02148,\psi_{0}\approx 0.04350\\}$ with $T=0.02858$ and
$Q=1.318$ , sharing the same $F(T,Q)$, $T$ and $Q$ with the corresponding RN
black hole.With the seed solution and the map Eq. (47), for given specific
$T_{i}$ and $Q_{i}$, we can obtain a series of numerical scalarized charged
black hole solution under the constrains $\\{T=T_{i},Q=Q_{i},\psi_{+}=0\\}$ by
the shooting method. Working with $(T,\Delta F)$ plane and $(Q,\Delta F)$
plane respectively and choosing the fixed parameters
$\\{c=0.02148,q=1.1m,m=0.01\\}$, we start with the seed solution
$\\{r_{h_{0}}=1,\phi_{1_{0}}=0.8,\psi_{0_{0}}\approx 0.04350\\}$ and
respectively obtain two series of numerical solutions
$\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}$ with $T=T_{i}=0.02858,Q=Q_{i}$ and
$Q=Q_{i}=1.318,T=T_{i}$ by shooting method,

$\displaystyle\\{T=T_{i}=0.02858,Q=Q_{i},\psi_{+}=0\\}\xrightarrow[\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}]{\text{shooting}}$
(155) $\displaystyle\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}\,,$ (158)
$\displaystyle\\{T=T_{i},Q=Q_{i}=1.318,\psi_{+}=0\\}\xrightarrow[\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}]{\text{shooting}}$
$\displaystyle\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}\,.$

For holding accuracy, we set the step size that $Q_{i}=Q_{i-1}-\frac{1}{100}$
from $Q_{0}=1.318$ and $T_{\pm i}=T_{\pm(i-1)}\pm\frac{1}{1000}$ from
$T_{0}=0.02858$ respectively. In the case of canonical ensemble, one can read
off $F(T_{i},Q_{i})$ of the corresponding RN black hole with temperature
$T_{i}$ and total charge $Q_{i}$ from Eq. (75) and Eq. (78). Therefore, we
numerically build a relation between $\Delta
F(T_{i},Q_{i})=F(T_{i},Q_{i})-F_{\rm RN}(T_{i},Q_{i})$ and $T$ or $Q$.

###  B.3 Grand Canonical Ensemble

We firstly consider the smaller $c=0.00101$ case. Recall the approach we
adopted previously, fixing parameter $\\{c=0.00101,q=1.1m\\}$ and starting
with the seed solution $\\{r_{h_{0}}=1,\phi_{1_{0}}=0.8,\psi_{0_{0}}\approx
0.002047\\}$ with $\\{T_{0}=0.02864,\mu=0.8000\\}$, we obtain two series of
numerical solution $\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}$ under the
constrains $\\{T=T_{i},\mu=\mu_{i},\psi_{+}=0\\}$ by shooting method,

$\displaystyle\\{T=T_{i}=0.02864,\mu=\mu_{i},\psi_{+}=0\\}\xrightarrow[\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}]{\text{shooting}}$
(160) $\displaystyle\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}\,,$ (163)
$\displaystyle\\{T=T_{i},\mu=\mu_{i}=0.8000,\psi_{+}=0\\}\xrightarrow[\\{r_{h_{i-1}},\phi_{1_{i-1}},\psi_{0_{i-1}}\\}]{\text{shooting}}$
$\displaystyle\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}\,.$

Here we set the step size that $T_{i}=T_{i-1}+\frac{1}{100}$,
$T_{-i}=T_{-(i-1)}-\frac{1}{100}$ while $\mu_{\pm
i}=\mu_{\pm(i-1)}\pm\frac{1}{100}$ for maintaining numerical precision and
stability.

As to the case of $c=0.1$, we briefly illustrate our results. Following the
procedure we used above and starting with the seed solution
$\\{r_{h_{0}}=1,\phi_{1_{0}}=0.8,\psi_{0_{0}}=0.1987\\}$ with temperature
$T_{0}=0.02721$ and chemical potential $\mu_{0}=0.8114$, we obtain another two
series of numerical solutions $\\{r_{h_{i}},\phi_{1_{i}},\psi_{0_{i}}\\}$
under the constrains $\\{T=T_{i},\mu=\mu_{i},\psi_{+}=0\\}$, where the step
size is set as $\mu_{\pm i}=\mu_{\pm(i-1)}\pm\frac{1}{100}$ and
$T_{i}=T_{i-1}+\frac{1}{100}$, $T_{-i}=T_{-i-1}-\frac{1}{500}$. For every
$T_{i}$ and $\mu_{i}$, we pick out $\Delta G=G(T,Q)-G_{\rm RN}(T,Q)$ from Eq.
(75) and Eq. (78).

## References

  * (1) B. P. Abbott et al. [LIGO Scientific and Virgo], _Observation of Gravitational Waves from a Binary Black Hole Merger_ , Phys. Rev. Lett. 116, no.6, 061102 (2016), [arXiv:1602.03837 [gr-qc]].
  * (2) K. Akiyama et al. [Event Horizon Telescope], _First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole_ , Astrophys. J. 875, no.1, L1 (2019), [arXiv:1906.11238 [astro-ph.GA]].
  * (3) R. Ruffini and J. A. Wheeler, _Introducing the black hole_ Phys. Today 24, no.1, 30 (1971)
  * (4) J. D. Bekenstein, _Transcendence of the law of baryon-number conservation in black hole physics,_ Phys. Rev. Lett. 28, 452-455 (1972)
  * (5) A. E. Mayo and J. D. Bekenstein, _No hair for spherical black holes: Charged and nonminimally coupled scalar field with selfinteraction_ , Phys. Rev. D 54, 5059-5069 (1996) [arXiv:gr-qc/9602057 [gr-qc]].
  * (6) X. H. Feng, H. Lu and Q. Wen, _Scalar Hairy Black Holes in General Dimensions_ , Phys. Rev. D 89, no.4, 044014 (2014) [arXiv:1312.5374 [hep-th]].
  * (7) H. S. Liu and H. Lü, _Scalar Charges in Asymptotic AdS Geometries_ , Phys. Lett. B 730, 267-270 (2014) [arXiv:1401.0010 [hep-th]].
  * (8) Z. Y. Fan and H. Lu, _Static and Dynamic Hairy Planar Black Holes_ , Phys. Rev. D 92, no.6, 064008 (2015) [arXiv:1505.03557 [hep-th]].
  * (9) M. Khodadi, A. Allahyari, S. Vagnozzi and D. F. Mota, _Black holes with scalar hair in light of the Event Horizon Telescope_ , JCAP 09, 026 (2020) [arXiv:2005.05992 [gr-qc]].
  * (10) Z. Y. Fan and H. Lu, _Charged Black Holes with Scalar Hair_ , JHEP 09, 060 (2015) [arXiv:1507.04369 [hep-th]].
  * (11) R. A. Konoplya and A. Zhidenko, _Analytical representation for metrics of scalarized Einstein-Maxwell black holes and their shadows_ , Phys. Rev. D 100, no.4, 044015 (2019) [arXiv:1907.05551 [gr-qc]].
  * (12) I. Z. Stefanov, S. S. Yazadjiev and M. D. Todorov, _Phases of 4D scalar-tensor black holes coupled to Born-Infeld nonlinear electrodynamics_ , Mod. Phys. Lett. A 23, 2915-2931 (2008) [arXiv:0708.4141 [gr-qc]].
  * (13) P. Wang, H. Wu and H. Yang, _Scalarized Einstein-Born-Infeld-scalar Black Holes_ , [arXiv:2012.01066 [gr-qc]].
  * (14) Y. S. Myung and D. C. Zou, _Scalarized black holes in the Einstein-Scalar-Maxwell theory with quasi-topological term_ , [arXiv:2011.09665 [gr-qc]].
  * (15) C. A. R. Herdeiro and E. Radu, _Asymptotically flat black holes with scalar hair: a review_ , Int. J. Mod. Phys. D 24, no.09, 1542014 (2015) [arXiv:1504.08209 [gr-qc]].
  * (16) J. P. Hong, M. Suzuki and M. Yamada, _Spherically Symmetric Scalar Hair for Charged Black Holes_ , Phys. Rev. Lett. 125, no.11, 111104 (2020) [arXiv:2004.03148 [gr-qc]].
  * (17) J. P. Hong, M. Suzuki and M. Yamada, _Charged black holes in non-linear Q-clouds with O(3) symmetry_ , Phys. Lett. B 803, 135324 (2020) [arXiv:1907.04982 [gr-qc]].
  * (18) C. A. R. Herdeiro and E. Radu, _Spherical electro-vacuum black holes with resonant, scalar $Q$-hair_, Eur. Phys. J. C 80, no.5, 390 (2020) [arXiv:2004.00336 [gr-qc]].
  * (19) M. M. Caldarelli, G. Cognola and D. Klemm, _Thermodynamics of Kerr-Newman-AdS black holes and conformal field theories_ , Class. Quant. Grav. 17, 399-420 (2000) [arXiv:hep-th/9908022 [hep-th]].
  * (20) J. X. Mo, X. X. Zeng, G. Q. Li, X. Jiang and W. B. Liu, _A unified phase transition picture of the charged topological black hole in Horava-Lifshitz gravity_ , JHEP 10, 056 (2013) [arXiv:1404.2497 [gr-qc]].
  * (21) M. Zhang, _Thermodynamical stability of $f(R)$-AdS black holes in grand canonical ensemble_, Gen. Rel. Grav. 51, no.1, 13 (2019) [arXiv:1812.04220 [gr-qc]].
  * (22) H. Quevedo, _Geometrothermodynamics_ , J. Math. Phys. 48, 013506 (2007) [arXiv:physics/0604164 [physics]].
  * (23) H. Quevedo, M. N. Quevedo, A. Sanchez and S. Taj, _On the ensemble dependence in black hole geometrothermodynamics_ , Phys. Scripta 8, 084007 (2014), [arXiv:1304.3954 [gr-qc]].
  * (24) G. W. Gibbons and S. W. Hawking, _Action Integrals and Partition Functions in Quantum Gravity,_ ” _Phys. Rev. D_ 15, 2752-2756 (1977)
  * (25) J. D. Brown, G. L. Comer, E. A. Martinez, J. Melmed, B. F. Whiting and J. W. York, Jr., _Thermodynamic Ensembles and Gravitation_ , Class. Quant. Grav. 7, 1433-1444 (1990)
  * (26) H. W. Braden, J. D. Brown, B. F. Whiting and J. W. York, Jr., _Charged black hole in a grand canonical ensemble_ , Phys. Rev. D 42, 3376-3385 (1990)
  * (27) S. A. Teukolsky, _Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations_ , Astrophys. J. 185, 635-647 (1973)
  * (28) W. H. Press and S. A. Teukolsky, _Perturbations of a Rotating Black Hole. II. Dynamical Stability of the Kerr Metric_ , Astrophys. J. 185, 649-674 (1973)
  * (29) C. V. Vishveshwara, _Stability of the schwarzschild metric_ , Phys. Rev. D 1, 2870-2879 (1970)
  * (30) R. A. Konoplya, A. Zhidenko and A. F. Zinhailo, _Higher order WKB formula for quasinormal modes and grey-body factors: recipes for quick and accurate calculations_ , Class. Quant. Grav. 36, 155002 (2019) [arXiv:1904.10333 [gr-qc]].
  * (31) P. Pani, _Advanced Methods in Black-Hole Perturbation Theory_ , Int. J. Mod. Phys. A 28, 1340018 (2013), [arXiv:1305.6759 [gr-qc]].
  * (32) R. A. Konoplya and A. Zhidenko, _Quasinormal modes of black holes: From astrophysics to string theory_ , Rev. Mod. Phys. 83, 793-836 (2011) [arXiv:1102.4014 [gr-qc]].
  * (33) E. Berti, V. Cardoso and A. O. Starinets, _Quasinormal modes of black holes and black branes_ , Class. Quant. Grav. 26, 163001 (2009) [arXiv:0905.2975 [gr-qc]].
  * (34) S. A. Hartnoll, C. P. Herzog and G. T. Horowitz, _Holographic Superconductors_ , JHEP 12, 015 (2008) [arXiv:0810.1563 [hep-th]].
  * (35) S. W. Hawking, _Black Holes and Thermodynamics,_ Phys. Rev. D 13, 191-197 (1976),
  * (36) J. W. York, Jr., _Black hole thermodynamics and the Euclidean Einstein action,_ Phys. Rev. D 33, 2092-2099 (1986),
  * (37) B. F. Whiting and J. W. York, Jr., _Action Principle and Partition Function for the Gravitational Field in Black Hole Topologies,_ Phys. Rev. Lett. 61, 1336 (1988),
  * (38) J. D. Brown and J. W. York, Jr., _The Path integral formulation of gravitational thermodynamics,_ [arXiv:gr-qc/9405024 [gr-qc]].
  * (39) H. W. Braden, J. D. Brown, B. F. Whiting and J. W. York _Charged black hole in a grand canonical ensemble,_ Phys. Rev. Lett. 42, 3376 (1990),
  * (40) J. D. Brown, G. L. Comer E. A. Martinez, J. Melmed, B. F. Whiting and J. W. York, _Thermodynamic ensembles and gravitation,_ Class. Quantum Grav. 7, 1433 (1990).,
  * (41) J. D. Brown, J. Creighton and R. B. Mann, _Temperature, energy and heat capacity of asymptotically anti-de Sitter black holes,_ Phys. Rev. D 50, 6394-6403 (1994), [arXiv:gr-qc/9405007 [gr-qc]].
  * (42) S. W. Hawking and D. N. Page, _“Thermodynamics of Black Holes in anti-De Sitter Space,_ Commun. Math. Phys. 87, 577 (1983),
  * (43) G. L. Comer, _Ensemble dependence of the stability of thermal black holes,_ Class. Quantum Grav. 9, 947 (1992).,
  * (44) A. Kusenko and M. E. Shaposhnikov, _Supersymmetric Q balls as dark matter_ , _Phys. Lett. B_ 418 (1998) 46–54, [hep-ph/9709492].
  * (45) A. de Gouvea, T. Moroi and H. Murayama, _Cosmology of supersymmetric models with low-energy gauge mediation_ , _Phys. Rev. D_ 56 (1997) 1281–1299, [hep-ph/9701244].
  * (46) M. Mars, _Present status of the Penrose inequality,_ , _Class. Quant. Grav. 26, 193001 (2009)_, [arXiv:0906.5566 [gr-qc]].
  * (47) S. A. Hayward, “Inequalities relating area, energy, surface gravity and charge of black holes,” Phys. Rev. Lett. 81, 4557-4559 (1998) [arXiv:gr-qc/9807003 [gr-qc]].
  * (48) E. Malec and N. O. Murchadha, “Trapped surfaces and the Penrose inequality in spherically symmetric geometries,” Phys. Rev. D 49, 6931-6934 (1994) [arXiv:gr-qc/9401014 [gr-qc]].
  * (49) G. W. Gibbons, “Some comments on gravitational entropy and the inverse mean curvature flow,” Class. Quant. Grav. 16, 1677-1687 (1999) [arXiv:hep-th/9809167 [hep-th]].
  * (50) Marcus A. Khuri, “A Penrose-like inequality with charge,” Gen. Relativ. Gravit 45, no. 11, 2341-2361 (2013)
  * (51) Gilbert Weinstein and Sumio Yamada, “On a Penrose Inequality with Charge,” Communications in Mathematical Physics 257, 703-723 (2005)
  * (52) S. W. Hawking and S. F. Ross, _Duality between electric and magnetic black holes_ , Phys. Rev. D 52, 5865-5876 (1995) [arXiv:hep-th/9504019 [hep-th]].
  * (53) S. Chandrasekhar and S. L. Detweiler, _The quasi-normal modes of the Schwarzschild black hole_ ,” Proc. Roy. Soc. Lond. A 344 (1975), 441-452
  * (54) H. Furuhashi and Y. Nambu, _Instability of massive scalar fields in Kerr-Newman space-time_ ,” Prog. Theor. Phys. 112, 983-995 (2004) [arXiv:gr-qc/0402037 [gr-qc]].
  * (55) B. F. Schutz and C. M. Will, _Black hole normal modes - A semianalytic approach_ , Astrophysical Journal, Part 2, vol. 291, April 15, 1985, p. L33-L36
  * (56) S. Iyer and C. M. Will, _Black Hole Normal Modes: A WKB Approach. 1. Foundations and Application of a Higher Order WKB Analysis of Potential Barrier Scattering_ ,” Phys. Rev. D 35, 3621 (1987)
  * (57) S. Iyer, _BLACK HOLE NORMAL MODES: A WKB APPROACH. 2. SCHWARZSCHILD BLACK HOLES_ ,” Phys. Rev. D 35, 3632 (1987)
  * (58) R. A. Konoplya, _Quasinormal behavior of the d-dimensional Schwarzschild black hole and higher order WKB approach_ , Phys. Rev. D 68, 024018 (2003) [arXiv:gr-qc/0303052 [gr-qc]].
  * (59) J. Matyjasek and M. Opala, _Quasinormal modes of black holes. The improved semianalytic approach_ , Phys. Rev. D 96, no.2, 024011 (2017), [arXiv:1704.00361 [gr-qc]].
  * (60) V. Ferrari and B. Mashhoon, _New approach to the quasinormal modes of a black hole_ , Phys. Rev. D 30, 295-304 (1984)