Scalarized Einstein-Maxwell-scalar Black Holes in Anti-de Sitter Spacetime

In this paper, we study spontaneous scalarization of asymptotically anti-de Sitter charged black holes in the Einstein-Maxwell-scalar model with a non-minimal coupling between the scalar and Maxwell fields. In this model, Reissner-Nordstr\"om-AdS (RNAdS) black holes are scalar-free black hole solutions, and may induce scalarized black holes due to the presence of a tachyonic instability of the scalar field near the event horizon. For RNAdS and scalarized black hole solutions, we investigate the domain of existence, perturbative stability against spherical perturbations and phase structure. In a micro-canonical ensemble, scalarized solutions are always thermodynamically preferred over RNAdS black holes. However, the system has much rich phase structure and phase transitions in a canonical ensemble. In particular, we report a RNAdS BH/scalarized BH/RNAdS BH reentrant phase transition, which is composed of a zeroth-order phase transition and a second-order one.

The no-hair theorem states that a black hole can be uniquely determined via three parameters, its mass, electric charge and angular momentum [1][2][3]. Although this theorem holds true in the Einstein-Maxwell field theory, it suffers from challenges due to the existence of hairy black holes possessing extra macroscopic degrees of freedom in other theories. In fact, various black hole solutions, e.g., hairy black holes in the Einstein-Yang-Mills theory [4][5][6], black holes with Skyrme hairs [7,8] and black holes with dilaton hairs [9], have been discovered as counter-examples to the no-hair theorem. For a review, see [10].
Spontaneous scalarization typically occurs in models with non-minimal coupling terms of scalar fields, which can source the scalar fields to destabilize scalar-free black hole solutions and form scalarized hairy black holes. This phenomenon was first studied for neutron stars in scalar-tensor models [11] by coupling scalar fields to the Ricci curvature. It demonstrated that there is a coexistence region for scalar-free and scalarized neutron star solutions, where the scalarized ones can be energetically preferred. Later, it was found that there also exists spontaneous scalarization of black holes in scalar-tensor models, provided that black holes are surrounded by non-conformally invariant matter [12,13].
Recently, the phenomenon of spontaneous scalarization has been studied in the extended Scalar-Tensor-Gauss-Bonnet (eSTGB) gravity [14][15][16][17][18][19]. In particular, asymptotically anti-de Sitter (AdS) scalarized black holes have been studied in a scalar-tensor model with non-minimally coupling the scaler field to the Ricci scalar and the Gauss-Bonnet term, as well as applications to holographic phase transitions [20]. In eSTGB models, the scalar field is non-minimally coupled to the Gauss-Bonnet curvature correction of the gravitational sector, which leads to numerical challenges for solving the evolution equations due to non-linear curvature terms. To better understand dynamical evolution into scalarized black holes, a similar, but technically simpler, type of models, i.e., Einstein-Maxwell-scalar (EMS) models with non-minimal couplings between the scalar and Maxwell fields, have been put forward in [21], where fully non-linear numerical evolutions of spontaneous scalarization were presented. Subsequently, further studies of spontaneous scalarization in the EMS models were discussed in the context of various non-minimal coupling functions [22,23], dyons including magnetic charges [24], axionic-type couplings [25], massive and self-interacting scalar fields [26,27], horizonless reflecting stars [28], stability analysis of scalarized black holes [29][30][31][32][33], higher dimensional scalar-tensor models [34], quasinormal modes of scalarized black holes [35,36], two U(1) fields [37], quasi-topological electromagnetism [38], topology and spacetime structure influences [39] and the Einstein-Born-Infeld-scalar theory [40]. Besides the above asymptotically flat scalarized black holes, spontaneous scalarization was also discussed in the EMS model with a positive cosmological constant [41]. Additionally, spontaneous vectorization of electrically charged black holes was also proposed [42], analytic treatments were applied to study spontaneous scalarization in the EMS models [43][44][45][46], and an infinite family of exact topological charged hairy black hole solutions were constructed in the EMS gravity system [47].
Studying thermodynamics of the EMS models not only provides evidence for endpoints of the dynamical evolution of unstable scalar-free black holes, but also is interesting per se. Understanding the statistical mechanics of black holes has been a subject of intensive study for several decades. In the pioneering work [48][49][50], Hawking and Bekenstein found that black holes possess the temperature and the entropy. However, asymptotically flat black holes are often thermally unstable since they have negative specific heat. To make black holes thermally stable, appropriate boundary conditions need to be imposed, e.g., putting the black holes in AdS space. Asymptotically AdS black holes become thermally stable since the AdS boundary acts as a reflecting wall. The thermodynamic properties of AdS black holes were first studied by Hawking and Page [51], who discovered the Hawking-Page phase transition between Schwarzschild AdS black holes and thermal AdS space. Later, motivated by AdS/CFT correspondence [52][53][54], there has been much interest in studying phase structure and transitions of AdS black holes [55][56][57][58][59][60][61][62][63][64]. In light of this, it is of great interest to study spontaneous scalarization of asymptotically AdS black holes and associated thermodynamic properties in the EMS models with non-minimal couplings of the scalar and electromagnetic fields.
The remainder of this paper is organized as follows. In section II, we introduce the EMS model with a negative cosmological constant and derive the free energy in a canonical ensemble. Section III is devoted to discussing linear perturbations in scalar-free and scalarized black hole solutions. In section IV, we present our main numerical results, including the domain of existence, entropic preference, effective potentials for radial perturbations, and phase structure and transitions in a canonical ensemble. We summarize our results with a brief discussion in section V.

II. EMS MODEL IN ADS SPACE
In this section, we derive the equations of motion, asymptotic behavior, the Smarr relation and the Helmholtz free energy for asymptotically AdS scalarized black hole solutions in the EMS model. The action of the EMS model with a negative cosmological constant is where we take G = 1 for simplicity throughout this paper. In the action (1), the scalar field φ is minimally coupled to the metric g µν and non-minimally coupled to the gauge field A µ , F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field strength tensor, Λ = −3/L 2 is the cosmological constant with the AdS radius L, and f (φ) is the non-minimal coupling function of the scalar and the gauge fields.

A. Equations of motion
Varying the action (1) with respect to the metric g µν , the scalar field φ and the gauge field A µ , one obtains the equations of motion, The energy-momentum tensor T µν in eqn.
(2) is given by In the following, we focus on the spherically symmetric ansatz for the metric, the electromagnetic field and the scalar field, Plugging the above ansatz into the equations of motion (2) yields where primes denote the derivatives with respect to the radial coordinate r, and the integration constant Q can be interpreted as the electric charge of the black hole solution. For later use, we introduce the Misner-Sharp mass function m (r) by N (r) = 1 − 2m (r) /r + r 2 /L 2 .

B. Asymptotic behavior
To obtain non-trivial hairy black hole solutions of the non-linear ordinary differential equations (5), one should impose appropriate boundary conditions at the event horizon and the spatial infinity. Accordingly, in the vicinity of the event horizon at r = r + , we find that the solutions can be approximated as a power series expansion in terms of (r − r + ), where The two essential parameters, φ 0 and δ 0 , can be determined after matching the asymptotic expansion of the solutions at the spatial infinity, where f (0) = 1 is assumed, M is identified as the ADM mass, φ + can be interpreted as the expectation value of the dual operator in the conformal boundary theory from AdS/CFT correspondence, and Φ is the electrostatic potential with Φ = ∞ r+ dre −δ(r) Q/ r 2 f (φ (r)) . Note that the scalar field usually behaves as φ (r) ∼ φ − + φ+ r 3 in the asymptotic region, and we set φ − = 0 in this paper, corresponding to the absence of the external source in the conformal boundary theory. Consequently, we can use the shooting method to solve the non-linear differential equations (5) for solutions satisfying the asymptotic expansions at the boundaries. It is also noteworthy that there is a scaling symmetry among the physical quantities, which allows us to solve eqn. (5) numerically in terms of redefined dimensionless quantities.

C. Smarr relation
The Smarr relation [65] can be used to test the accuracy of numerical scalarized black hole solutions, since it associates the black hole mass with other physical quantities. The Smarr relation can be derived from computing the Komar integral for a time-like Killing vector K µ = (1, 0, 0, 0) in a manifold M. Integrating the identity ∇ µ (∇ ν K µ ) = K µ R µν over the time constant hypersurface Σ, whose boundary ∂Σ consists of the event horizon r = r + and the spatial infinity r = +∞, one can use Gauss's law to obtain where dS µν is the surface element on ∂Σ, and dS µ is the volume element on Σ, accordingly. Making use of eqns. (3) and (5), we find that the Smarr relation is given by where A H = 4πr 2 + is the horizon area, and T = N (r + ) e −δ(r+) /4π is the Hawking temperature. For a RNAdS black hole with δ (r) = 0, the Smarr relation (11) reduces to where the last term is the P V term in the extended phase space of AdS black holes [66].

D. Free energy
Given a family of scalarized black holes, it is of interest to compute the Helmholtz free energy, which can be used to investigate phase structure and transitions in a canonical ensemble with fixed charge Q and temperature T . The free energy, which is identified as the thermal partition function of black holes, can be derived via constructing the Euclidean path integral. In the semiclassical approximation, the partition function is evaluated by exponentiating the on-shell Euclidean action S E on-shell , where the on-shell action S E on-shell is obtained by substituting the classical solution into the action. However, the onshell action S E on-shell normally diverges in asymptotically AdS spacetime. One then needs holographic renormalization to remove divergences appearing in the asymptotic region. There are several methods to regularize S E on-shell , such as the background-subtraction method [57] and the Kounterterms method [67][68][69]. Here, we adopt the counterterm subtraction method [70,71] to regularize the action by adding a series of boundary terms to the bulk action.
Specifically for the aforementioned bulk action S bulk in eqn. (1), the regularized action S R is supplied with three boundary terms where S GH is the Gibbon-Hawking boundary term to render the variational problem well-defined, S ct includes counterterms to eliminate divergences on asymptotic boundaries, and S surf is used to fix the charge rather than the electrostatic potential when the action is varied [62,72]. The three boundary terms are given by where the integrals are performed on the hypersurface at the spatial infinity, γ is the determinant of the induced metric on the hypersurface, Θ is the trace of the extrinsic curvature, R 3 is the scalar curvature of the induced metric γ, and n µ is the unit vector normal to the hypersurface. Using the equations of motion (2) and the asymptotic expansion (8), we obtain the on-shell Euclidean version of S bulk , S GH , S ct and S surf , where S = πr 2 + is the entropy of the black hole. Consequently, the regularized on-shell Euclidean action for the black hole solution (4) is The Helmholtz free energy F is related to the Euclidean action S E on-shell via which gives

III. PERTURBATIONS AROUND BLACK HOLE SOLUTION
In this section, we investigate linear perturbations around black hole solutions, which can help us understand the stability of the solutions.

A. Scalar perturbation around RNAdS black holes
We first examine a scalar perturbation δφ in a RNAdS black hole background. Note that ifḟ (0) = 0 is imposed, a RNAdS black hole, which is described by is manifestly a solution of the equations of motion (2). In this scalar-free solution background, we can linearize the scalar equation in eqn.
(2) with a scalar perturbation δφ, where µ 2 ef f = −f (0) Q 2 / 2r 4 . In a (3 + 1)-dimensional asymptotically AdS spacetime of AdS radius L, a scalar field can cause a tachyonic instability only if its mass-squared is less than the so-called Breitenlohner-Freedman (BF) bound µ 2 BF = −9/ 4L 2 [73]. For the scalar perturbation δφ, one always has µ 2 ef f > µ 2 BF for large enough r, and hence, asymptotically, the RNAdS black hole is stable against the formation of the scalar field, which guarantees that scalarized black holes induced by the tachyonic instability of the scalar field are asymptotically AdS. However, if µ 2 ef f < µ 2 BF in some region (e.g., near the event horizon), a RNAdS black hole may evolve to a scalarized black hole under a scalar perturbation. Note that iff (0) < 0, one always has µ 2 ef f > µ 2 BF , and hence a tachyonic instability can not occur.
To study how a scalarized black hole solution bifurcates from a scalar-free black hole solution, we calculate zero modes of the scalar perturbation in the scalar-free black hole background. For simplicity, the scalar perturbation δφ is written as the decomposition with spherical harmonics functions, With this decomposition, the scalar equation (21) then reduces to Given the fixed values of l andf (0), requiring that the radial field U l (r) is regular at the event horizon and vanishes at the spatial infinity selects a family of discrete black hole solutions, which can be labelled by a non-negative integer node number n. In this paper, we focus on the l = 0 = n fundamental mode since it gives the smallest q of the black hole solutions [21]. Due to the tachyonic instability, scalarized RNAdS black holes may emerge from these zero modes, which composes bifurcation lines in the domain of existence for the scalarized black holes.

B. Time-dependent perturbation around scalarized black holes
To investigate perturbative stability of the scalarized black hole solution (4), we then consider spherically symmetric and time-dependent linear perturbations. Specifically including the perturbations, the metric ansatz is written as [22] ds 2 = −Ñ (r, t) e −2δ(r,t) dt 2 + dr 2 N (r, t) + r 2 dθ 2 + sin 2 θdϕ 2 , where the time dependence of the perturbations is assumed to be Fourier modes with frequency Ω. Similarly, the ansatzes of the scalar and electromagnetic fields are given bỹ respectively. Solving eqn. (2) with the ansatzes (24) and (25), we can extract a Schrodinger-like equation for the perturbative scalar field φ 1 (r), where Ψ (r) ≡ rφ 1 (r), and the tortoise coordinate r * is defined by dr * /dr ≡ e δ(r) N −1 (r). Here, the effective potential U Ω is given by One can show that the effective potential U Ω vanishes at the event horizon, whereas approaches positive infinity at the spatial infinity. From quantum mechanics, the existence of an unstable mode with Ω 2 < 0 requires the presence of U Ω < 0 in some regions. Nevertheless, a positive definite U Ω ensures that scalarized black hole solutions are stable against the spherically symmetric perturbations. It is noteworthy that the appearance of a negative region of U Ω cannot sufficiently guarantee the presence of an instability [26]. One can utilize other techniques, like the S-deformation method [74], to further discuss the stability of these solutions.

IV. NUMERICAL RESULTS
In this section, we first study various properties of scalarized RNAdS black hole solutions and then investigate their phase structure and transitions in a canonical ensemble. After the non-linear differential equations (5) are expressed in terms of a new dimensionless coordinate they are numerically solved for scalarized black hole solutions using the NDSolve function in Wolfram Mathematica.
In the remainder of this paper, we focus on the coupling function f (φ) = e αφ 2 with α > 0. For this coupling function, one has f (0) = 1 andḟ (0) = 0, which ensures that RNAdS black holes are solutions of the EMS model, andf (0) > 0, which may trigger a tachyonic instability of the scalar field to induce scalarized black hole solutions. For later use, we define reduced quantities, which are dimensionless and invariant under the scaling symmetry (9). To test the accuracy of our numerical method, we use the Smarr relation (11) and find that the numerical error can be maintained around the order of 10 −6 .

A. Scalarized black holes
Here, we present the numerical results, e.g., the domain of existence, entropic preference and effective potentials, for scalarized black hole solutions, which are dynamically induced from RNAdS black holes. Without loss of generality, we focus onQ = Q/L = 0.1 in this subsection. For a RNAdS black hole, a tachyonic instability of the scalar field occurs if µ 2 ef f < µ 2 BF somewhere in the spacetime. Since the minimum value of µ 2 ef f occurs at the event horizon r = r + , we only need to check µ 2 ef f < µ 2 BF at r = r + . The region in the (α, q) parameter space of RNAdS black holes where µ 2 ef f < µ 2 BF at r = r + is plotted in the upper left panel of Fig. 1. The distribution of values of (µ 2 ef f − µ 2 BF )Q 2 | r=r+ is also displayed, which shows that the tachyonic instability region becomes larger as α increases, and the scalar field suffers from a strong tachyonic instability when black holes are near-extremal. The bifurcation line is composed of the l = 0 = n zero modes of eqn. (23) for the scalar perturbation in RNAdS black holes, and represented by blue dashed lines in Fig. 1. When the tachyonic instability is strong enough, the scalar perturbation can lead to scalarized black holes with non-trivial scalar fields above the bifurcation line. In the upper right panel of Fig. 1, we display the domain of existence for scalarized black holes, which is exhibited by a light blue region. The domain of existence is bounded by the bifurcation and critical lines, and resembles that of RN scalarized black holes [21]. On the critical line, the mass and the charge of scalarized solutions remain finite, whereas its horizon radius vanishes. On the other hand, the mass-to-charge ratio q of RNAdS black holes reaches the maximum in the extremal limit, which is shown by a horizontal dashed gray line. Moreover, there is a certain region bounded by the extremal and bifurcation lines, where scalarized and RNAdS black holes coexist.
The reduced area a H as a function of reduced charge q is plotted for RNAdS and scalarized black holes in the lower left panel of FIG. 1, which demonstrates that scalarized black holes emerge from RNAdS black holes at the bifurcation points, marked by B, and eventually terminate on the critical line with zero a H . For a multiphase system in a micro-canonical ensemble with conserved energy, the phase of maximum entropy is globally stable and will be present at equilibrium. Therefore, in the scalarized and RNAdS black holes coexisting region, our numerical results show that scalarized solutions are entropically preferred over RNAdS black hole solutions, and hence are the globally stable phase in the micro-canonical ensemble.
To study the stability of scalarized solutions, the effective potentials U Ω of scalarized solutions with α = 5 are plotted for several values of q in the lower right panel of FIG. 1, where solid and dashed colored lines correspond to potentials with and without negative regions, respectively. The scalarized solutions have positive effective potentials for a large enough value of q, and thus are free of radial instabilities. However, as q decreases towards the bifurcation line, there appear negative regions in the effective potentials, which means that radial instabilities cannot be excluded near the bifurcation line.

B. Phase structure in a canonical ensemble
In this subsection, we consider phase structure and transitions of scalarized and RNAdS black holes in a canonical ensemble maintained at a given temperature of T and a given charge of Q. In a canonical ensemble, the globally ef f − µ 2 BF Q 2 evaluated at the event horizon r = r+ as a function of q and α for the scalar perturbation in RNAdS black holes. We only display the tachyonically unstable region, where µ 2 ef f < µ 2 BF , in the α-q plane. The closer RNAdS black holes are to the extremal limit, the more unstable the scalar field becomes. The blue dashed line represents the bifurcation line, where tachyonic instabilities are strong enough to induce scalarized black holes. Upper Right: Domain of existence for scalarized RNAdS black holes in the α-q plane, which is highlighted by the light blue region and bounded by the critical and bifurcation lines. The critical line is depicted by a red solid line, on which the reduced horizon area aH vanishes. The horizontal dashed gray line denotes extremal RNAdS black holes, above which RNAdS black hole solutions do not exist. Lower Left: Reduced horizon area aH against q for RNAdS and scalarized black holes. The scalarized black hole solutions are always entropically preferred, which means that they are globally stable in a micro-canonical ensemble. Lower Right: Effective potentials of scalarized black holes with α = 5 for several values of q. Solid red lines denote positive definite effective potentials between the event horizon and the spatial infinity, while dashed red lines represent those possessing negative regions. When q is large enough, the scalarized black hole solutions are stable against radial perturbations. stable phase of a multiphase system, which exists at equilibrium, has the lowest possible Helmholtz free energy F , which can be computed via eqn. (19). The rich phase structure of black holes usually comes from expressing the horizon radius r + as a function of temperature T . If the function r + (T ) is multivalued, there will be more than one black hole phase, corresponding to different branches of r + (T ).
To illustrate phase structure and transitions, we plot the reduced horizon radiusr + and the free energyF as functions of reduced temperatureT for scalarized and RNAdS black holes with three representative values ofQ in Fig. 2, where we have α = 5. In the left column of Fig. 2 with a smallQ, the upper panel shows that three branches of the RNAdS black hole solution coexist in some range ofT , and are dubbed as large, intermediate and small RNAdS BHs, respectively, based on their values of horizon radius. At a high (low) enough temperature, only the large (small) RNAdS BH phase exists. On the other hand, there is only one phase for the scalarized black hole solution, which bifurcates from the RNAdS black hole solution at the bifurcation point B, and does not exist at a low temperature. The reduced free energyF is plotted againstT for these four phases in the lower panel, which shows that the scalarized black hole can not be the globally stable phase since there always exists some RNAdS black hole phase of a lower free energy at a givenT . In the coexisting region of the RNAdS black hole phases, a first-order phase transition between large and small RNAdS BHs occurs at some point, where their free energies intersect each other. There is a band of temperatures where three branches of RNAdS black hole solutions coexist, and a first-order phase transition occurs between the large RNAdS BH phase (i.e., the branch with the largest horizon radius) and the small RNAdS BH phase (i.e., the branch with the smallest horizon radius). Scalarized black holes emerge from the bifurcation point. Nevertheless, they are not globally preferred since they always have a higher free energy than RNAdS black holes. Center column: RNAdS black hole solutions have only one branch, whose free energy is smaller than that of scalarized black holes. There is no phase transition. Right column: RNAdS black hole solutions have only one branch, whereas scalarized black hole solutions have two branches of different sizes. AsT increases, the globally stable phase jumps from RNAdS black holes to scalarized black holes (the branch with a larger horizon radius), corresponding to a zeroth-order phase transition atTmin. Further increasingT , there would be a second-order phase transition returning to RNAdS black holes at the bifurcation point B. Here, we observe a RNAdS BH/scalarized BH/RNAdS BH reentrant phase transition.
The coexisting region of the RNAdS black hole phases shrinks asQ increases until reaching a critical point, where a second-order phase transition occurs between large and small RNAdS BHs. Beyond the critical point, large RNAdS BH is indistinguishable from small RNAdS BH, hence RNAdS black hole solutions have a single phase. In the center column of Fig. 2, we depictr + (T ) andF (T ) for the case with a value ofQ greater than the critical value. The upper panel shows that both RNAdS and scalarized black holes have a single phase. Moreover, as displayed in the lower panel, the RNAdS black hole always has a smaller free energy than the scalarized black hole, and hence is the globally stable phase. Consequently, there is no phase transition in this case.
However whenQ is large enough, phase structure of scalarized and RNAdS black holes becomes much richer. For example,r + (T ) andF (T ) are plotted for scalarized and RNAdS black holes with a large enoughQ in the right column of Fig. 2. It shows in the upper panel that the RNAdS black hole solution possesses a single phase, whereas the scalarized solution can have two phases at some givenT , namely large scalarized BH (i.e., the one with a larger horizon radius) and small scalarized BH (i.e., the one with a smaller horizon radius). In fact, the scalarized black hole solution has a minimum temperatureT min , and large scalarized BH coexists with small scalarized BH betweeñ T =T min and the bifurcation point B, where large scalarized BH and the RNAdS black hole merge. The lower panel exhibits the free energy as a function ofT for the three phases. If we start increasing the temperature fromT = 0, the system follows the blue line of the RNAdS black hole untilT =T min , where the free energy has a discontinuity at its global minimum. Further increasingT , the inset shows that the system jumps to the lower green line of large scalarized BH, which corresponds to a zeroth-order phase transition between scalarized and RNAdS black holes. As T continues to increase, the system follows the lower green line until it joins the blue line at the bifurcation point B, which corresponds to a second-order phase transition between scalarized and RNAdS black holes. Note that since scalarized and RNAdS black holes have the same entropy at the bifurcation point, a phase transition between them at the bifurcation point is second-order. In short, a RNAdS BH/scalarized BH/RNAdS BH reentrant phase transition Here, we take α = 5. The phase diagram exhibits the globally stable phases, which have the lowest free energy, and the phase transitions between them. The light yellow and blue regions correspond to RNAdS and scalarized black holes, respectively. A first-order phase transition line (the purple line) separates the large RNAdS BH phase, which is above the line, and the small RNAdS BH phase, which is below the line. The first-order phase transition line terminates at the critical point, labelled by C. The scalarized black hole phase is delimited by a zeroth-order phase transition line (the red line) and a second-order one (the blue dashed line), which coincides with the bifurcation line.
is observed asT increases.
In addition, it is interesting to consider the local thermodynamic stability of black hole phases against thermal fluctuations. In a canonical ensemble, the quantity of particular interest is the specific heat at constant charge, Since the entropy is proportional to the size of the black hole, a positive specific heat means that black holes radiate less when they are smaller. Thus, the thermal stability of a phase follows from C Q > 0 (or equivalently, ∂r + /∂T > 0). From the upper row of Fig. 2, we notice that large and small RNAdS BHs in the left column, RNAdS black hole in the center column, and large scalarized BH and RNAdS black hole in the right column all have ∂r + /∂T > 0. In consequence, the globally stable phases of scalarized and RNAdS black holes possess a positive C Q , and are thermally stable.
To better illustrate the globally stable phases of lowest free energy and the associated phase transitions, we display the phase diagram of scalarized and RNAdS black holes in theQ-T plane in Fig. 3, where α = 5. There is a firstorder phase transition line (the purple line) separating large and small RNAdS BHs for smallQ, which terminates at the critical point. This first-order phase transition is quite similar to the liquid/gas phase transition. WhenQ is large enough, large scalarized BH (the light blue region) appears, and is bounded by zeroth-order (the red line) and second-order (the blue dashed line) phase transition lines.

V. DISCUSSION AND CONCLUSIONS
In this paper, we investigated spontaneous scalarization of asymptotically AdS charged black holes in the EMS model, and studied phase structure of scalarized and RNAdS black holes in a canonical ensemble. We focused on a non-minimal coupling function f (φ) = e αφ 2 , which leads to spontaneous scalarization due to the tachyonic instability of the scalar field near the event horizon. In practice, scalarized black holes can be induced from RNAdS black holes on the bifurcation line, which consists of zero modes of the scalar perturbation in RNAdS black holes. In the α-q plane with a fixedQ, the domain of existence for scalarized RNAdS black holes is bounded by the bifurcation and critical lines, which resembles that of scalarized RN black holes very closely [21]. In a micro-canonical ensemble, we found that scalarized black hole solutions are always entropically preferred over RNAdS black holes, and hence the globally stable phase.
On the other hand, the system has much richer phase structure in a canonical ensemble. After the Helmholtz free energy of the EMS model was computed, we obtained the phase structure of scalarized and RNAdS black holes. In the small Q regime, scalarized black holes never globally minimize the free energy, and the corresponding phase diagram resembles that of the liquid/gas system closely. Nevertheless in the large Q regime, scalarized black holes can be the globally stable phase in some parameter region. As the temperature increases at a given charge, the system undergoes a RNAdS BH/scalarized BH/RNAdS BH reentrant phase transition, which consists of zeroth-order and second-order phase transitions.
The phenomenon of reentrant phase transition was first observed in a nicotine/water mixture, and later discovered in the context of black hole thermodynamics, e.g., Born-Infeld-AdS black holes [62,75], higher dimensional singly spinning Kerr-AdS black holes [76], AdS black holes in Lovelock gravity [77], AdS black holes in dRGT massive gravity [78], hairy AdS black holes [79]. For these black holes, reentrant phase transitions were found to include zeroth-order and first-order phase transitions. In this paper, we present an example of a reentrant phase transition for black holes, which is composed of zeroth-order and second-order phase transitions. Furthermore, the second-order phase transition between scalarized and RNAdS black holes is of great interest since this implies that our results may provide an interesting model of holographic superconductors. We leave this for future work.
AdS black holes can also been studied in the context of extended phase space thermodynamics, where the cosmological constant is interpreted as a thermodynamic pressure P ≡ 6/L 2 [60,80]. In terms of P , the reduced quantities are expressed asT = T 6/P ,F = F P/6,r + = r + P/6,Q = Q P/6,M = M P/6.
Note that M and F are identified as the enthalpy and the Gibbs free energy, respectively, in extended phase space. With eqn. (31), our results can be generalized to extended phase space thermodynamics.