Maxwell construction and multi-criticality in uncharged generalized quasi-topological black holes

It is well known that black holes behave like thermodynamic systems that radiate a thermal flux of particles [1] and whose equilibrium states are governed by the four laws of black hole mechanics [2], reinterpreted as the four laws of thermodynamics. For asymptotically AdS black holes, the negative cosmological constant can be identified as a thermodynamic pressure that validates the extended first law of thermodynamics [3–6]. In this context black hole mass is interpreted as enthalpy, and thus the mechanics of black holes can be understood in terms of chemical thermodynamics [7, 8]. Over the past decade, the perspective of black hole chemistry has led to the discovery of a number of rich properties, including van der Waals phase transitions [9], re-entrant phase transitions [10, 11], superfluid-like phase transitions [12–14], and triple points [11, 15–17].

A recent interesting discovery in black hole chemistry was that of multi-critical points. These were first seen in Einstein gravity coupled to non-linear electrodynamics [18], but were shortly afterward found to be present in multiply rotating Kerr-AdS black holes [19], and in Lovelock gravity [20]. In the latter case multi-critical behaviour can even occur for asymptotically flat black holes [21]. An Nth order multi-critical point occurs when N distinct phases merge at a single value of pressure and temperature, generalizing the notion of a triple point (with N = 3). Explicit examples of quadruple and quintuple points have been found in all cases examined so far.

In this paper, we demonstrate an alternative method for finding multi-critical points in black hole phase transitions. Previous methods exploited the fact that each extremum of the temperature (regarded as a function of horizon radius) corresponded to a cusp in the Gibbs free energy [18, 20]. $2N-2$ distinct extrema were thus required to obtain N distinct phases, each with its own swallowtail in the Gibbs free energy diagram. The intersection points of corresponding swallowtails will merge if two adjacent inflections of $T(r_+)$ occur at the same temperature. Extending this to all the inflections will then yield an N-tuple critical point where all such intersections merge. However, the number of parameters used in this approach is more than needed—in fact only approximately half of the horizon radii of the extrema are needed. This redundancy results in extremely low efficiency in computing the (finely tuned) parameters need to obtain a multicritical point, as iteratively tuning the various parameters such that all inflections occur at the same temperature becomes very time-consuming. This problem is particularly acute when dealing with more complicated higher-curvature theories. Indeed the situation deteriorates when considering black holes with large numbers of thermodynamic degrees of freedom since considerably high precision is required because of the finely tuned nature of multi-criticality.

The new method overcomes these difficulties. Inspired by Maxwell's equal area rule [22], we provide an equivalent but more compact description of multi-criticality, which can be perfectly adapted to the construction without the redundancy introduced in previous methods. Instead of iteratively manipulating input parameters, our approach directly indicates whether a given set of non-redundant parameters can or cannot have multi-critical points. If the latter holds, we obtain their accurate values in thermodynamic phase space with notably less computation.

To illustrate our method we consider black holes in higher curvature theories of gravity. Specifically we consider generalized quasi-topological (GQT) gravity. The reasons for this are as follows.

For the past 2 decades there has been a revival of interest in higher curvature gravity in the theoretical physics community. Such theories have proven to be significant in a variety of contexts in physics, including string theory, holography, the AdS/CFT correspondence, tests of general relativity (GR), and black hole thermodynamics. The gravitational action becomes renormalizable when supplemented with higher-curvature terms [23], making such theories candidates for a quantum theory of gravity. String-theoretic versions of quantum gravity motivate the possibility of higher dimensional spacetimes, and the addition of higher-curvature corrections allows for broader generalizations of the Einstein–Hilbert action to dimensions larger than four. These higher-curvature theories provide toy models for studying the AdS/CFT correspondence and allow for holographic study of conformal field theories (CFTs).

But the inclusion of higher order curvature terms comes at a price—it can yield equations of motion containing higher derivatives in the metric that give rise to instabilities and negative energy modes [24, 25]. Notably, this inconsistency with GR is absent in some classes of higher-curvature theories [26–29] in which only a massless spin-2 graviton can propagate to infinity. This subclass of higher-curvature gravity theories is considerably more promising than the others and thus warrants further investigation.

GQT gravities, a class of recently proposed higher-derivative theories, satisfies the requirements noted above. Theories in the GQT class characterize generalizations of GR in any dimension and to any order in curvature insofar as they contain non-hairy black hole solutions and second-order-differential equations for the metric in any linearized maximally symmetric background [30–32]. Generally, the bulk part of their action can be written as

$$\begin{align} I_{\mathrm{bulk}} = \frac{1}{16\pi}\int_{\mathcal{V}} \mathrm{d}^dx\sqrt{-g}\left[-2\Lambda+R+\sum_{n = 2} \sum_k \alpha_{n,k}\mathcal{S}_n^{\left(k\right)}\right], \end{align} \tag{ 1.1 }$$

where Λ is the cosmological constant, $\mathcal{S}_n^{(k)}$ are independent densities constructed from different constructions of n Riemann tensors and the metric, $\alpha_{n,k}$ is the associated kth higher-curvature coupling, and the Newtonian constant G is set to 1 for simplicity. Quite remarkably, in the context of gravitational effective field theory, any higher-curvature theory can be mapped into a subset of GQT theories via field redefinition [30]. Furthermore, there exists a subset in the parameter space of higher-curvature couplings where these theories only allow a massless spin-2 graviton propagating at the linearized level.

In general, the field equations of this class contain metric derivatives up to fourth order. However, for a static spherically symmetric (SSS) ansatz (1.2),

$$\begin{align} \mathrm{d}s^2 = -f\left(r\right)\mathrm{d}t^2+\frac{1}{f\left(r\right)}\mathrm{d}r^2+r^2\mathrm{d}\Omega_{d-2, \kappa}^2, \end{align} \tag{ 1.2 }$$

where $\mathrm{d}\Omega_{d-2, \kappa}^2$ describes the $(d-2)$-dimensional transverse surface of constant curvature normalized to $\kappa = +1,0,-1$ denoting spherical, flat and hyperbolic topologies, respectively. The metric function f(r) is fully determined by the vanishing of a total derivative of a second order differential or algebraic expression, with the constant of integration related to the ADM mass. Theories with an algebraic equation of motion for f are identified as quasi-topological (QT) gravities [33–35], and they likely satisfy a Birkhoff theorem [35–38]. However, theories in the QT class only exist in $d\unicode{x2A7E} 5$. Another interesting GQT subclass is that of Lovelock theories [26, 27], which are the most direct generalizations of GR in that the field equations are always second order differential equations for any metric. Similar to the QT class, for the ansatz (1.2), the metric function f is fully characterized by a single algebraic equation that only differs from that of the QT theories by an overall constant. Thus Lovelock gravity corresponds to a subclass of QT theories.

However, Lovelock theories seem too restrictive—Einstein gravity is the only possible Lovelock theory in d = 4, and a Lovelock curvature density of order n yields non-trivial dynamics only if $d\gt2n+1$. So the first non-trivial Lovelock theory appears in d = 5, corresponding to Gauss–Bonnet gravity with n = 2, whereas non-trivial GQT gravity theories exist for any order $n\unicode{x2A7E}3$ in $d\unicode{x2A7E} 4$ [30, 31]. In addition, as far as (1.2) is concerned, GQT theories allow multiple inequivalent densities at a given order in $d\unicode{x2A7E} 5$, but QT (Lovelock) theories have one unique density at any order [31]. Thus GQTs constitute a much broader class of higher-curvature theories than have already been specified.

We therefore choose to illustrate our approach in GQT gravity theories. The thermodynamics of these theories have not been explored to the same extent that Lovelock theories have. Previous studies have been carried out in limited contexts, with only a few low order couplings in dimensions not much larger than four [39–44]. Inspired by the multi-critical behaviour found for a broad range of black holes in different contexts [18–21], our interests lie both in demonstrating our method and in understanding multi-criticality in GQT gravity.

Our paper is organized as follows. In section 2 and appendix A we review some properties of GQT black hole solutions. In section 3, we provide an interpretation of multi-criticality and introduce the K-rule, obtained by the reformulation of the Maxwell area rule. We develop in section 4 a method based on the K-rule, carry out a discussion on its feasibility, and construct quadruple and quintuple points based on this approach.

Evaluated on a SSS metric of the form (1.2), the GQT class of curvature order $n\unicode{x2A7E} 2$ has exactly n − 1 inequivalent densities in $d\unicode{x2A7E} 5$ [31]. As required by the integrability of the field equation for f(r), on-shell GQT Lagrangian densities should be total derivatives of the form [31]

$$\begin{align} \mathcal{S}_n^{\left(k\right)} = \sum_{j = 0}^{n}\lambda_{n,j}^{\left(k\right)}\mathcal{S}_{\left(n,j\right)},\quad k = 1,\ldots, k_n\equiv\mathrm{max}\left(1,n-1\right), \quad n\unicode{x2A7E} 1 \end{align} \tag{ 2.1 }$$

where k labels one of the n − 1 inequivalent densities, $\mathcal{S}_{(n,j)}$ is given by

$$\begin{align} \mathcal{S}_{\left(n,j\right)}\equiv\frac{1}{r^{{\,}d-2}}\frac{\mathrm{d}}{\mathrm{d} r}\left[r^{{\,}d-1}\left(\frac{\kappa-f}{r^2}\right)^{n-j}\left(\frac{-f^{^{\prime}}}{2r}\right)^j\right], \end{align} \tag{ 2.2 }$$

and $\lambda_{n,j}^{(k)}$ are some constrained coefficients such that $\mathcal{S}_{n}^{(k)}$ is induced by a real off-shell GQT density. Then the constraints

$$\begin{align} \sum_{j = 0}^n\lambda_{n,j}^{\left(k\right)}\left(2n-dj\right) = 0,\quad \sum_{j = 0}^n\lambda_{n,j}^{\left(k\right)}j\left(j-1\right) = 0. \end{align} \tag{ 2.3 }$$

are required to obtain the most general Lagrangian density [31]. The first constraint ensures that all densities contribute with a power of $r^{{\,}d-1}$; the second ensures that the field equations are of 2nd order when linearized about constant curvature backgrounds. Note that any linear combination of on-shell densities still satisfies (2.1)–(2.3); therefore any density can be decomposed into n − 1 independent densities. In other words, it is sufficient to study one particular choice of $\lambda_{n,j}^{(k)}$.

Incorporating the constraints (2.3), the following choice

$$\begin{align} \lambda_{n,k+1}^{\left(k\right)} & = \left(k-1\right)\left(dk-4n\right),\quad \lambda_{n,k}^{\left(k\right)} = 2d\left(1-k^2\right)-4n\left(1-2k\right),\nonumber\\ \lambda_{n,k-1}^{\left(k\right)} & = k\left[d\left(k+1\right)-4n\right], \end{align} \tag{ 2.4 }$$

of $\lambda_{n,j}^{(k)}$ form a family of GQT densities, with the remaining coefficients identical to 0. This simple choice is employed in our analysis but the results in the remaining part of this section hold generally.

Upon integrating the equations of motion, we obtain [30]

$$\begin{align} \frac{16\pi M}{\left(d-2\right)\Omega_{d-2,\kappa}} = \sum_{n = 0}^{n_{\mathrm{max}}}\alpha_{n,k}\sum_{k = 1}^{k_n}\sum_{j = 0}^{n}\lambda_{n,j}^{\left(k\right)}\mathcal{F}_{\left(n,j\right)}, \end{align} \tag{ 2.5 }$$

where $\Omega_{d-2,\kappa} = \frac{2\pi^{\frac{d-1}{2}}}{\Gamma(\frac{d-1}{2})}$, M is the ADM mass [45–48], and

$$\begin{align} \mathcal{F}_{\left(n,j\right)} & = \frac{\left(-1\right)^{j+1}}{2^{j+1}}r^{d-2+j-2n}\left(\kappa-f{\,}\right)^{n-j-1}\left(f^{^{\prime}}\right)^{j-2}\nonumber\\ & \quad \times\left\{f^{^{\prime}}\left[{\,}j\left(d-1+j-2n\right)\left(\kappa-f{\,}\right)f-\left(j-1\right)r\left(\kappa+\left(n-j-1\right)f{\,}\right)f^{^{\prime}}\right]\right.\nonumber\\ & \quad + \left. j\left(j-1\right)r\left(\kappa-f{\,}\right)ff^{^{\prime\prime}}\right\} \end{align} \tag{ 2.6 }$$

where prime is defined as taking derivatives with respect to r.

For $j = 0, 1$, $\mathcal{F}_{(n,j)}$ becomes an algebraic quantity (having no derivative terms), which implies the density defined by (2.4) with k = 1 is purely QT (Lovelock). To make a consistent notation with common definitions, we separate QT densities from GQT ones (2.4) and define

$$\begin{align} \mathcal{S}^{\mathrm{QT}}_{n}\equiv -\frac{1}{2}\mathcal{S}^{\left(1\right)}_{n} = +\left(2n-d\right)\mathcal{S}_{\left(n,0\right)}-2n\mathcal{S}_{\left(n,1\right)}, \end{align} \tag{ 2.7 }$$

as QT densities. The remainder, with $k\unicode{x2A7E}2$, are genuine GQT densities.

Despite the complexity of (2.5), the thermodynamics of GQT black holes is determined by two simple equations: $f(r_+) = 0$, which defines the outermost black hole horizon at $r = r_+$, and $f^{^{\prime}}(r_+) = 4\pi T$, which defines the temperature of the black hole. These two relations are given by [31]

$$\begin{align} \begin{aligned} M& = \frac{\Omega_{d-2,\kappa}}{16\pi G}\sum_{n = 0}^{n_{\mathrm{max}}}\sum_{k = 1}^{k_n}\sum_{j = 0}^{n}\alpha_{n,k}\lambda_{n,j}^{\left(k\right)}\left(j-1\right)\kappa^{n-j}r_+^{d-2n-1}\left(-2\pi r_+T\right)^j, \end{aligned} \end{align} \tag{ 2.8 }$$

$$\begin{align} \begin{aligned} 0 & = \sum_{n = 0}^{n_{\mathrm{max}}}\sum_{k = 1}^{k_n}\sum_{j = 0}^{n}\alpha_{n,k}\lambda_{n,j}^{\left(k\right)}\left(d-2n+j-1\right)\kappa^{n-j}\left(-2\pi r_+T\right)^jr_+^{d-2n-2}, \end{aligned} \end{align} \tag{ 2.9 }$$

where the couplings of the lowest two orders are set to be $\alpha_{0,1}\equiv -(d-1)(d-2)/(2\ell^2)$, and $\alpha_{1,1}\equiv 1/2$ for consistency with Einstein gravity. The parameter $\ell$ here is the AdS length, and so $\alpha_{0,1}$ is identical to the cosmological constant Λ. In the context of black hole chemistry, all couplings except for $\alpha_{1,1}$ are identified as thermodynamic variables, with

$$\begin{align} \Phi_{n}^{\left(k\right)}\left(r_+,T\right)\equiv\frac{\partial M\left(r_+,T,\left\{\alpha_{p,l}\right\}\right)}{\partial \alpha_{n,k}} \end{align} \tag{ 2.10 }$$

as their corresponding conjugate potentials.

As a well-known consequence in GR, the first law of thermodynamics and the Smarr relation hold as well in the extended phase space:

$$\begin{align} \mathrm{d}M & = T\mathrm{d}S+V\mathrm{d}P+\sum_{n\unicode{x2A7E}2}\Phi_{n}^{\left(k\right)}\mathrm{d}\alpha_{n,k}, \end{align} \tag{ 2.11 }$$

$$\begin{align} \left(d-3\right)M & = \left(d-2\right)TS-2PV+2\sum_{n\unicode{x2A7E}2}\left(n-1\right)\Phi_{n}^{\left(k\right)}\alpha_{n,k}, \end{align} \tag{ 2.12 }$$

where the pressure P is a redefinition of $\alpha_{0,1}$ and V is its corresponding conjugate volume

$$\begin{align} P\equiv\frac{\left(d-1\right)\left(d-2\right)}{16\pi\ell^2},\quad V\equiv\frac{\Omega_{d-2,\kappa}r_+^{d-1}}{d-1}\; . \end{align} \tag{ 2.13 }$$

The thermodynamic quantity S is the Wald entropy [49], which reads [31]

$$\begin{align} S = \frac{\Omega_{d-2,\kappa}}{8G}\sum_{n = 0}^{n_{\mathrm{max}}}\sum_{k = 1}^{k_n}\sum_{j = 0}^{n}\alpha_{n,k}\lambda_{n,j}^{\left(k\right)}j\kappa^{n-j}\left(-2\pi r_+T\right)^{j-1}r_+^{d-2n}. \end{align} \tag{ 2.14 }$$

This quantity is not always positive, and in such situations it has been common to simply discard solutions for which this is the case. However ambiguities exist in the definition of the black hole entropy. For example, adding to the Lagrangian a term proportional to the induced metric on the horizon will, without having an effect on the other properties of the solution, shift the entropy by an arbitrary constant. One example is that of adding an Euler density to the action [44]. We shall therefore retain solutions with S < 0 in considerations, appropriately indicating in our figures where this occurs.

Henceforth we shall consider the Gibbs free energy $G = M-TS$ for investigating phase transitions. The global minimum of G yields the most stable thermodynamic phase at any given temperature.

Before continuing, we note that physical theories should only propagate one type of massless spin-2 graviton on constant curvature backgrounds. This in turn implies the effective Newtonian constant must have the same sign as the one in GR, which means that for the class of metrics (1.2) having asymptotically AdS solutions of the form

$$\begin{align} f\left(r\right) = \kappa+f_{\infty}\frac{r^2}{\ell^2}+\frac{m}{h^{^{\prime}}\left(f_{\infty}\right)r^{d-3}}+\cdots, \end{align} \tag{ 2.15 }$$

we shall only consider black holes with $f_{\infty}\gt0$, $h^{^{\prime}}(f_{\infty})\lt0$ and $\gamma^2\gt0$ ³ as satisfying the requisite physical criteria. We discuss these issues in appendix A.

We seek to obtain the conditions under which three or more phases merge at a particular temperature and pressure. The Gibbs free energy provides a diagnostic for this. Its global minimum as a function of the temperature T determines the thermodynamically stable state of the system for a given fixed choice of the other thermodynamic parameters. The presence of swallowtails in the Gibbs free energy indicates multiple phases, with first order phase transitions between two distinct phases taking place at the intersection point of the swallowtail. There must be N − 1 swallowtails in order to have N distinct phases. Whenever the intersection points of j different swallowtails coincide, then there is a jth order multicritical point, where $j \unicode{x2A7D} N$.

Previous methods for finding multiple phases and N-tuple critical points exploited the fact that, regarding temperature as a function of horizon radius, each extremum of $T(r_+)$ corresponds to a cusp in the Gibbs free energy [18, 20]. Hence N distinct phases require $2N-2$ distinct extrema. If two adjacent inflections of $T(r_+)$ occur at the same temperature, then the intersection points of corresponding swallowtails will merge. If this takes place for all the inflections, then all such intersection points will merge, corresponding to an N-tuple critical point. These critical points can be found by finely tuning the other thermodynamic parameters.

Here we demonstrate an alternate method that is considerably more efficient. We start with a brief review of the Maxwell construction [22]. It is well-known that the multiplicity of the Gibbs free energy $G(P,T)$ corresponds to the non-monotonic behaviour of the pressure $P(V,T)$. As illustrated in figure 1, a full oscillation AaBbC in the pressure at a fixed temperature $T^{{\,}*}$ leads to a swallowtail on the Gibbs phase diagram. With some abuse of notation, integrating $\mathrm{d}G$ along the loop $\mathrm{A}\rightarrow \mathrm{b} \rightarrow \mathrm{a}\rightarrow \mathrm{C}$ yields

$$\begin{align} \begin{aligned} 0 = \oint \mathrm{d}G = \oint\frac{\partial G}{\partial P}\bigg|_{T^{{\,}*}}\mathrm{d}P = \int V\mathrm{d}P = PV\bigg|_{V_{\mathrm{A}}}^{V_{\mathrm{C}}}-\int P\mathrm{d}V = \int_{V_{\mathrm{A}}}^{V_{\mathrm{C}}} \left(P^*-P\right)\mathrm{d}V. \end{aligned} \end{align} \tag{ 3.1 }$$

The second equality holds because the temperature is fixed, the fourth one comes from integration by parts, and the last expression follows from $P(V_{\mathrm{A}}) = P(V_{\mathrm{C}}) = P^*$, where $P^*$ characterizes the swallowtail intersection point in the Gibbs phase plot. The geometric interpretation of (3.1) is obvious: $P^*$ corresponds to a pressure that partitions the oscillatory parts of the P − V diagram into equal areas.

Zoom In Zoom Out Reset image size

It is useful to define the function $K(V,V_\mathrm{i})$ and its derivative $K^{^{\prime}}(V,V_\mathrm{i})$ as follows:

$$\begin{align} K\left(V,V_\mathrm{i}\right)\equiv \int_{V_\mathrm{i}}^{V} \left(P^*-P\right)\mathrm{d}V, \quad K^{^{\prime}}\left(V,V_\mathrm{i}\right)\equiv\frac{\partial K}{\partial V} = P^*-P\left(V,T^{{\,}*}\right). \end{align} \tag{ 3.2 }$$

It is obvious that $K(V_{\mathrm{A}},V_{\mathrm{A}}) = 0$ and that the last expression of (3.1) can be rewritten as

$$\begin{align} K\left(V_{\mathrm{C}},V_{\mathrm{A}}\right) = 0. \end{align} \tag{ 3.3 }$$

However for any two points in thermodynamic phase space whose volumes satisfy (3.3), the relation (3.3) alone does not imply that their difference in free energy is zero. It is also necessary to ensure that $P(V_{\mathrm{C}}) = P(V_{\mathrm{A}})$ so that (3.1) holds. In a plot of G vs. P, this requirement is equivalent to the condition that A and C are the same point. This in turn implies that

$$\begin{align} K^{^{\prime}}\left(V_{\mathrm{A}},V_{\mathrm{A}}\right) = K^{^{\prime}}\left(V_{\mathrm{C}},V_{\mathrm{A}}\right) = 0 \end{align} \tag{ 3.4 }$$

where the first condition ensures that $P(V_{\mathrm{A}},T^{{\,}*}) = P^*$. Geometrically, (3.3) and (3.4) ensure that A and C are the same point in the Gibbs energy diagram, and the continuity of K between A and C guarantees that this point is on some closed loop. A true self-intersection point (or double point) therefore emerges.

We note that if the second derivative Kʹʹ

$$\begin{align} K^{^{\prime\prime}} \equiv\frac{\partial^2 K}{\partial V^2} = -\frac{\partial P}{\partial V} \end{align} \tag{ 3.5 }$$

vanishes at some point then K will no longer be an extremum there. This is illustrated for point C in the rightmost diagram of figure 1 by the red curve. The pressure will then be an extremum at this point (as shown in the leftmost diagram in figure 1), and the corresponding part of the curve in the free energy diagram will get reflected through $P^*$, as shown by the red curve in the middle diagram in figure 1. By convention, we still regard this as a double point.

These considerations can be easily generalized to any N-tuple point. We say that an N-tuple point exists at $(P^*,T^{{\,}*})$ if and only if the K-rule is satisfied: namely that the function $K(V,V_0)$ has N real zero points $\{V_n\}$ for some fixed V₀ and Kʹ vanishes for all those roots, namely

$$\begin{align} K\left(V_n,V_0\right) = 0, \quad K^{^{\prime}}\left(V_n,V_0\right) = 0, \quad n = 1,2,\ldots,N \end{align} \tag{ 3.6 }$$

are satisfied by exactly N different values of $\{V_n\}$, including V₀ itself. Since the above argument about multicriticality is quite general, we would expect these discussions apply to any thermodynamic system for any conjugate pair of thermodynamic quantities, such as the temperature and the entropy.

We shall now construct multiple phases and N-tuple points for GQT black holes based on the K-rule introduced in section 3. The procedure is simple.

• 1.  
  Write the function K as

  $$\begin{align} K\left(r_+,r_0\right) = \int_{r_0}^{r_+} \left[P^*-P\left(\tilde{r}_+,T^{{\,}*},\left\{\alpha_{n,k}\right\}\right)\right]\frac{\mathrm{d} V\left(\tilde{r}_+\right)}{\mathrm{d} \tilde{r}_+}\mathrm{d}\tilde{r}_+ , \end{align} \tag{ 4.1 }$$

  where P is given from (2.9) (which can be regarded as the equation of state) with $\ell^2$ replaced by (2.13), V is identified as the thermodynamic volume defined in (2.13), and $\{\alpha_{n,k}\}$ is the set of undetermined couplings.
• 2.  
  Apply the K-rule to N positive distinct values of $r_+$ (where r₀ is taken to be any one of these values), then solve for $P^*, T^{{\,}*}, \{\alpha_{n,k}\}$ from the $2N-1$ independent equations ⁴ (3.6). This implies that a minimum number of $2N-3$ non-zero higher-curvature couplings are required.
• 3.  
  Check if the solution $P^*, T^{{\,}*}, \{\alpha_{n,k}\}$ provides a real N-tuple point in a sense that exactly N roots solve (3.6) as desired. If not, change the choice of as many horizon radii as needed until a real N-tuple point occurs.

We pause to make a few supplementary comments regarding the feasibility of the method. First of all, K is constructed from P and V instead of T and S because we want K to be a simple function such that the equations in step (2) are solvable: the definition (4.1) fulfills this requirement since K is in fact a polynomial in $r_+$. For convenience, K is defined as a function of the radius rather than the volume. It should be pointed out that the density $\mathcal{S}_{n\unicode{x2A7E} 3}^{(2)}$ is quasi-topological and becomes trivial in $d = 2n$. Therefore, for even dimensions, we exclude $\alpha_{d/2,2}$ from our considerations. We shall also restrict ourselves to spherical black holes with κ = 1 for simplicity. Since the pressure becomes a polynomial in $r_+$ (with the temperature T considered as a non-dynamical parameter), then Descartes' rule of signs can be applied, which relates the largest number of oscillations in the region $r_+\gt0$ to the number of sign changes in the sequence of a polynomial's coefficients. Thus we discuss the feasibility of our method through studying the possibility for the occurrence of N − 1 oscillations in P by manipulating signs of couplings in the next paragraph.

The feasibility of step (2) can be seen by induction. As indicated by (2.9) or table 1, switching on an arbitrary genuine GQT coupling ($k\unicode{x2A7E} 2$) always introduces three independent terms proportional to different powers of $r_+$ in the expression for the pressure in addition to

$$\begin{align} P_0 = \frac{\left(d-2\right) T}{4 r_+}-\frac{\left(d-3\right) \left(d-2\right) \kappa}{16 \pi r_+^2}, \end{align} \tag{ 4.2 }$$

which is the expression with all higher-curvature couplings set to zero. Since P₀ has one sign change, then turning on any particular coupling can introduce another sign change. Respecting the fact that temperature and the coupling are free, it is possible to make P have a full oscillation, which means that (3.6) has real solutions and a double point can be obtained. Similarly, for critical points involving more phases, we can obtain additional oscillations by including two additional couplings per oscillation, as long as they switch on at least two monomials in $r_+$ that differ from those already present. Hence not any choice of $2N-3$ higher-order densities yields an equilibrium state with four phases or more. For example, as bolded in table 1, switching on $\{\alpha_{7,2},\alpha_{8,4},\alpha_{9,6},\alpha_{10,8},\alpha_{11,10}\}$ and keeping other couplings zero only introduces $1/r_+^{11}$, $1/r_+^{12}$, $1/r_+^{13}$ into the pressure. Together with P₀, a total of five monomials are present in the expression for the pressure, which is only enough to construct a triple point. In this example we must therefore avoid turning on more than three couplings that contribute to the same three powers of $r_+$.

Table 1. The table shows a general pattern that the pressure follows when some genuine GQT densities ($k\unicode{x2A7E} 2$) are turned on. The top element of each column indicates the power of $1/r_+$, and each row contains all possible couplings with a constant sum of subscripts. The table tells what powers of $1/r_+$ in the expression of pressure are influenced by which couplings. For example, if $\alpha_{5,2}$ is the only non-zero coupling, then the pressure will contain three monomials of the form $1/r_+^7, 1/r_+^8, 1/r_+^9$ in additional to P₀; if $\alpha_{7,4}$ is turned on as well, then two extra monomials $1/r_+^{10},1/r_+^{11}$ will be introduced and both of $\alpha_{5,2}$ and $\alpha_{7,4}$ contribute to coefficient of the term $1/r_+^9$.

3

4

5

$\alpha_{3,2}$

$\alpha_{3,2}$

$\alpha_{3,2}$

6

7

$\alpha_{4,2}$

$\alpha_{4,2}$

$\alpha_{4,2}$

8

9

$\alpha_{4,3}$

$\alpha_{4,3}$

$\alpha_{4,3}$

$\alpha_{5,2}$

$\alpha_{5,2}$

$\alpha_{5,2}$

10

11

$\alpha_{5,3}$

$\alpha_{5,3}$

$\alpha_{5,3}$

$\alpha_{6,2}$

$\alpha_{6,2}$

$\alpha_{6,2}$

12

13

$\alpha_{5,4}$

$\alpha_{5,4}$

$\alpha_{5,4}$

$\alpha_{6,3}$

$\alpha_{6,3}$

$\alpha_{6,3}$

$\boldsymbol{\alpha_{7,2}}$

$\boldsymbol{\alpha_{7,2}}$

$\boldsymbol{\alpha_{7,2}}$

14

15

$\alpha_{6,4}$

$\alpha_{6,4}$

$\alpha_{6,4}$

$\alpha_{7,3}$

$\alpha_{7,3}$

$\alpha_{7,3}$

$\alpha_{8,2}$

$\alpha_{8,2}$

$\alpha_{8,2}$

16

17

$\alpha_{6,5}$

$\alpha_{6,5}$

$\alpha_{6,5}$

$\alpha_{7,4}$

$\alpha_{7,4}$

$\alpha_{7,4}$

$\alpha_{8,3}$

$\alpha_{8,3}$

$\alpha_{8,3}$

$\alpha_{9,2}$

$\alpha_{9,2}$

$\alpha_{9,2}$

18

19

$\alpha_{7,5}$

$\alpha_{7,5}$

$\alpha_{7,5}$

$\boldsymbol{\alpha_{8,4}}$

$\boldsymbol{\alpha_{8,4}}$

$\boldsymbol{\alpha_{8,4}}$

$\alpha_{9,3}$

$\alpha_{9,3}$

$\alpha_{9,3}$

$\alpha_{10,2}$

$\alpha_{10,2}$

$\alpha_{10,2}$

⋯

$\alpha_{7,6}$

$\alpha_{7,6}$

$\alpha_{7,6}$

$\alpha_{8,5}$

$\alpha_{8,5}$

$\alpha_{8,5}$

$\alpha_{9,4}$

$\alpha_{9,4}$

$\alpha_{9,4}$

$\alpha_{10,3}$

$\alpha_{10,3}$

$\alpha_{10,3}$

$\alpha_{11,2}$

⋯

$\alpha_{8,6}$

$\alpha_{8,6}$

$\alpha_{8,6}$

$\alpha_{9,5}$

$\alpha_{9,5}$

$\alpha_{9,5}$

$\alpha_{10,4}$

$\alpha_{10,4}$

$\alpha_{10,4}$

$\alpha_{11,3}$

$\alpha_{11,3}$

⋯

$\alpha_{8,7}$

$\alpha_{8,7}$

$\alpha_{8,7}$

$\boldsymbol{\alpha_{9,6}}$

$\boldsymbol{\alpha_{9,6}}$

$\boldsymbol{\alpha_{9,6}}$

$\alpha_{10,5}$

$\alpha_{10,5}$

$\alpha_{10,5}$

$\alpha_{11,4}$

$\alpha_{11,4}$

$\alpha_{11,4}$

⋯

$\alpha_{9,7}$

$\alpha_{9,7}$

$\alpha_{9,7}$

$\alpha_{10,6}$

$\alpha_{10,6}$

$\alpha_{10,6}$

$\alpha_{11,5}$

$\alpha_{11,5}$

$\alpha_{11,5}$

$\alpha_{12,4}$

⋯

$\alpha_{9,8}$

$\alpha_{9,8}$

$\alpha_{9,8}$

${\alpha_{10,7}}$

${\alpha_{10,7}}$

${\alpha_{10,7}}$

$\alpha_{11,6}$

$\alpha_{11,6}$

$\alpha_{11,6}$

$\alpha_{12,5}$

$\alpha_{12,5}$

⋯

$\boldsymbol{\alpha_{10,8}}$

$\boldsymbol{\alpha_{10,8}}$

$\boldsymbol{\alpha_{10,8}}$

${\alpha_{11,7}}$

${\alpha_{11,7}}$

${\alpha_{11,7}}$

$\alpha_{12,6}$

$\alpha_{12,6}$

$\alpha_{12,6}$

⋯

$\alpha_{10,9}$

$\alpha_{10,9}$

$\alpha_{10,9}$

$\alpha_{11,8}$

$\alpha_{11,8}$

$\alpha_{11,8}$

${\alpha_{12,7}}$

${\alpha_{12,7}}$

${\alpha_{12,7}}$

${\alpha_{13,6}}$

⋯

$\alpha_{11,9}$

$\alpha_{11,9}$

$\alpha_{11,9}$

$\alpha_{12,8}$

$\alpha_{12,8}$

$\alpha_{12,8}$

${\alpha_{13,7}}$

${\alpha_{13,7}}$

⋯

$\boldsymbol{\alpha_{11,10}}$

$\boldsymbol{\alpha_{11,10}}$

$\boldsymbol{\alpha_{11,10}}$

$\alpha_{12,9}$

$\alpha_{12,9}$

$\alpha_{12,9}$

$\alpha_{13,8}$

$\alpha_{13,8}$

$\alpha_{13,8}$

⋯

$\alpha_{12,10}$

$\alpha_{12,10}$

$\alpha_{12,10}$

$\alpha_{13,9}$

$\alpha_{13,9}$

$\alpha_{13,9}$

$\alpha_{14,9}$

⋯

$\alpha_{12,11}$

$\alpha_{12,11}$

$\alpha_{12,11}$

$\alpha_{13,10}$

$\alpha_{13,10}$

$\alpha_{13,10}$

$\alpha_{14,9}$

$\alpha_{14,9}$

⋯

$\vdots$

$\vdots$

$\vdots$

$\vdots$

$\vdots$

$\vdots$

$\ddots$

The number $2N-3$ should be considered as the minimum number of couplings required by our method. This may not be the smallest number of couplings needed for the emergence of an N-tuple point in general, since it is possible to have another sign change internally between the three monomials corresponding to the $\alpha_{n,k}$. The explicit form of new terms that are activated by switching on $\alpha_{n,k}$ reads

$$\begin{align} \begin{aligned} \left(-2\right)^{k-3}\alpha_{n,k} & \Bigg\{\left(d + k - 2 n\right)\left(-1 + k\right) \left(d k - 4 n\right) \frac{\pi^kT^{k+1}\kappa^{n-k-1}}{r_+^{2n-k-1}}\nonumber\\ & \quad + \left(d+k-2n-1\right) \left[d \left(-1 + k^2\right) + \left(2 - 4 k\right) n\right]\frac{\pi^{k-1}T^{k}\kappa^{n-k}}{r_+^{2n-k}}\nonumber\\ & \quad +\left(d + k - 2n-2\right)k\left(d + d k - 4 n\right) \frac{\pi^{k-2}T^{k-1}\kappa^{n-k+1}}{4r_+^{2n-k+1}}\Bigg\}. \end{aligned} \end{align} \tag{ 4.3 }$$

We claim that there can be at most one sign change between these three terms. If the sign switches twice, the first two coefficients must have a negative product, namely,

$$\begin{equation*} \left(d + k - 2 n\right)\left(-1 + k\right) \left(d k - 4 n\right)\left(d+k-2n-1\right) \left[d \left(-1 + k^2\right) + \left(2 - 4 k\right) n\right]<0. \end{equation*}$$

Since we restrict ourselves to genuine GQT densities ($k\unicode{x2A7E}2$) in $d\unicode{x2A7E} 5$ only, the above can be simplified to

$$\begin{align} \left(d + k - 2 n\right) \left(d+k-2n-1\right) \left(d k - 4 n\right) \left[d \left(-1 + k^2\right) + \left(2 - 4 k\right) n\right]<0. \end{align} \tag{ 4.4 }$$

Note that the dimension d cannot be any of $2n-k,2n-k+1,2n-k+2$; otherwise one of the terms would vanish and it would be no longer possible to have two sign changes between two terms. Therefore the first two factors in (4.4) must be positive, which leads to

$$\begin{equation*} \left(d k - 4 n\right) \left[d \left(-1 + k^2\right) + \left(2 - 4 k\right) n\right]<0, \end{equation*}$$

finally giving

$$\begin{align} \frac{dk}{4} < n < \frac{d\left(k^2-1\right)}{4k-2}, \end{align} \tag{ 4.5 }$$

where $k\unicode{x2A7E} 2$ is applied to determine the directions of inequalities. Meanwhile, we require the last two terms produce a sign change as well. Through a similar analysis, we arrive at

$$\begin{equation*} \left[d \left(-1 + k^2\right) + \left(2 - 4 k\right) n\right] \left(d + d k - 4 n\right)<0, \end{equation*}$$

which yields

$$\begin{align} \frac{d\left(k^2-1\right)}{4k-2}< n < \frac{d\left(k+1\right)}{4}. \end{align} \tag{ 4.6 }$$

Our claim is thus proved, since (4.6) contradicts (4.5). Even if there is an extra internal sign change, a half oscillation in the pressure does not necessarily occur, since these variables are not only integers but also constrained relative to each other in a complicated way. However, notice that the range of n in (4.6) is proportional to d, implying that we have a rather large parameter space for n and k. Hence we can still expect that there exists some choices of $d,k,n$ that can give rise to an N-tuple point with less than $2N-3$ couplings. As a consequence, we would also expect N − 1 to be a lower bound on the number of couplings needed for the occurrence of an N-tuple point. Because of these considerations, step (3) is added to guarantee that exactly N phases are obtained.

From the previous discussion, we can see that to figure out a general Gibbs phase rule is challenging. For neutral multi-rotating black holes [19], because the phase structure is invariant under the exchange of any two angular momenta, turning on any two additional couplings always creates a new phase and vice versa. However in the GQT scenario, this symmetry is broken between any two couplings, and theories with distinct values of k (even if they have the same n) differ a lot in their phase structures. A multi-critical point may not occur even if infinitely many couplings are turned on. This implies that the problem is not as simple as it is in Lovelock gravity where k = 1 and each density (2.7) similarly contributes to the thermodynamics [20]. We leave this question for future investigation.

In order to obtain multiple phases and multi-critical points, the physical constraints and ensuring P > 0 everywhere must also be considered. Since it is difficult to find a case with positive γ² everywhere, we only impose $\gamma^2 \gt 0$ in a neighbourhood of each critical point. In practice, we keep manipulating $r_+$'s until a critical point satisfying all constraints occurs. Under these considerations, we explicitly obtain a quadruple point (figures 2 and 3) and a quintuple point (figures 4 and 5) for two different spherical GQT black holes. Note that to see the merging of multiple swallowtails requires high precision in computations due to the finely tuned nature of multi-critical points. Both multi-critical points have negative Gibbs free energies, implying stable phase transitions. For the quintuple point, an extra coupling $\alpha_{3,2}$ is fixed to be 1 before running the procedure in order to make physical cases easier to emerge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared with the previous methods, the advantages of our procedure are quite significant. First, the application of the Maxwell construction turns the problem into algebra which enables our method to produce critical points with arbitrarily high precision such that very tiny phase structures can be discovered easily. More importantly, the method is more efficient in the sense that it does not need any fine-tuning procedure such as those previously employed [18–21] where N-tuple points were obtained by manipulating thermodynamic variables so that a common point of inflection occurred between multiple maxima and minima of the temperature as a function of $r_+$. However, due to the additional physical constraint $\gamma^2\gt0$ induced from the non-algebraic nature of equation of motions in genuine GQTGs, finding a physical multicritical point with a large N is still time-consuming.

We have exploited Maxwell's equal area law to find an interphase equilibrium for black holes with multiple phases, reformulating it into what we call the K-rule. Utilizing the K-rule, we developed a novel approach for constructing N-tuple points in the phase space of black holes.

We applied our results to GQT theories with $2N-3$ genuine couplings. Our analysis suggests that the minimum number of couplings required for the formation of an N-tuple point is likely confined between N − 1 and $2N-3$. We presented quadruple and quintuple points to illustrate the effectiveness of the method.

The discovery that black holes can have multiple phases and multicritical behaviour has interesting implications for quantum gravity. The distinct phases correspond to different thermodynamically stable states having different entropy, analogous to the way in which steam and ice are respectively high-entropy and low-entropy states of water. Systems having multicritical behaviour, such as colloids and polymers [50, 51], interact not only through short-range hard-sphere (or hard-cylinder) effects but also through additional soft or long-range interactions comprising multiple length sales [50]. It is reasonable to infer that the fundamental degrees of freedom of a black hole likewise have such complicated interactions, commensurate with recent indications that they are molecular in character [52, 53].

Future work would involve applying the K-rule to other kinds of black hole holes, particularly those whose horizon structures are not spherically symmetric. These include accelerating black holes in non-linear electrodynamics [54], multiply rotating black holes [21], and various black hole solutions in supergravity theories [55].

No new data were created or analysed in this study.

The equation of motion for f(r) with non-zero black hole mass admits asymptotically AdS solutions of the form

$$\begin{align} f\left(r\right) = \kappa+f_{\infty}\frac{r^2}{\ell^2}+\frac{m}{h^{^{\prime}}\left(f_{\infty}\right)r^{d-3}}+\cdots, \end{align} \tag{ A.1 }$$

with omitted terms decaying faster than $1/r^{d-3}$. The parameter m is related to the ADM mass M via

$$\begin{equation*} M = \frac{\left(d-2\right)\Omega_{d-2,\kappa}}{16\pi}m. \end{equation*}$$

The constant $f_{\infty}$ satisfies the polynomial equation

$$\begin{align} h\left(f_{\infty}\right)\equiv 1-f_{\infty}+\sum_{n\unicode{x2A7E} 2}^{n_{\mathrm{max}}}\sum_{k = 1}^{k_p}\frac{\left(-1\right)^{n}2\alpha_{n,k}}{\ell^{2n-2}} \frac{d-2n}{d-2}{f_{\infty}^{{\,}n}} = 0, \end{align} \tag{ A.2 }$$

whose first derivative with respect to $f_{\infty}$ is

$$\begin{align} h^{^{\prime}}\left(f_{\infty}\right) = -1+\sum_{n\unicode{x2A7E} 2}^{n_{\mathrm{max}}}\sum_{k = 1}^{k_p}\frac{\left(-1\right)^n2\alpha_{n,k}}{\ell^{2n-2}} \frac{d-2n}{d-2}nf^{{\,}n-1}_{\infty}. \end{align} \tag{ A.3 }$$

These conditions ensure that the graviton is not a ghost in a constant curvature background. In some theories $h^{^{\prime}}(f_{\infty})$ vanishes identically [56, 57], in which cases the asymptotic expansion (2.15) is no longer valid. We shall not consider this situation here, leaving it for future investigation.

For the QT class, a positive value of $f_{\infty}$ that solves (A.2) is sufficient to impose asymptotically AdS boundary conditions, but is insufficient for a theory with genuine GQT densities. This is because of the non-algebraic nature of equation (2.5). Since this is a differential equation, the expansion (2.15) only represents one particular solution for f(r). To gain a complete version of the asymptotic solution, we must find the homogenous part f_h , which solves the differential equation

$$\begin{align} f^{^{\prime\prime}}_h-\frac{4}{r}f^{^{\prime}}_h-\gamma^2r^{d-3}f_h = 0. \end{align} \tag{ A.4 }$$

Although the structure of (A.4) is general for any genuine GQT theories, the expression for the coefficient γ² here depends on the choice of higher-curvature couplings. For example, the expression

$$\begin{align} \gamma^2 = \frac{4\left[ h^{^{\prime}}\left(f_{\infty}\right)\right]^2 f_{\infty}^2}{3 \left(d-1\right) m \left[ h\left(f_{\infty}\right)+f_{\infty}-1\right]\ell^2} \end{align} \tag{ A.5 }$$

holds for theories that have $\{\alpha_{n,k = 2}\}$ as the only non-zero couplings, but becomes invalid if any coupling with $k\unicode{x2A7E}3$ is turned on.

Hence for GQT theories, in addition to the sign of $f_{\infty}$, we require that f_h have only decaying modes faster than $1/r^{d-3}$ to preserve the correct asymptotic condition at the boundary. Solving (A.4) at far region yields

$$\begin{equation*} f_h = \begin{cases} r^{5/2} \left[A I_{\frac{5}{\left(d - 1\right)}}\left( 2\gamma \frac{r^{\left(d-1\right)/2}}{d - 1} \right) + B K_ {\frac{5}{\left(d - 1\right)}}\left( 2\gamma \frac{r^{\left(d-1\right)/2}}{d - 1} \right) \right]& \text{if}~ \gamma^2>0\\ \\ r^{5/2} \left[C J_{\frac{5}{\left(d - 1\right)}}\left( 2\gamma \frac{r^{\left(d-1\right)/2}}{d - 1} \right) + D Y_ {\frac{5}{\left(d - 1\right)}}\left( 2\gamma \frac{r^{\left(d-1\right)/2}}{d - 1} \right) \right]& \text{if} ~\gamma^2<0, \end{cases} \end{equation*}$$

where A, B, C and D are constants fixed by boundary conditions and $I_k(x)$, $K_k(x)$, $J_k(x)$ and $Y_k(x)$ are Bessel functions of order k. In the large-r limit they become

$$\begin{equation*} f_h^{{\,}\left(+\right)}\sim Ar^{\frac{11-d}{4}}\exp\left(\frac{2|\gamma| r^{\frac{d-1}{2}}}{d-1}\right)+Br^{\frac{11-d}{4}}\exp\left(-\frac{2|\gamma| r^{\frac{d-1}{2}}}{d-1}\right), \end{equation*}$$

$$\begin{equation*} f_h^{{\,}\left(-\right)}\sim C r^{\frac{11-d}{4}} \cos{\left(\frac{2|\gamma|r^{\frac{d-1}{2}}}{d-1}+\frac{\left(3d-13\right)\pi}{4\left(d-1\right)}\right)}+D r^{\frac{11-d}{4}} \sin{\left(\frac{2|\gamma|r^{\frac{d-1}{2}}}{d-1}+\frac{\left(3d-13\right)\pi}{4\left(d-1\right)}\right)}, \end{equation*}$$

where the superscripts $(\pm)$ indicate the sign of γ². In order that the homogenous part is subdominant at large r, we require $\gamma^2 \gt 0$ and A = 0.

Apart from the correct asymptote, physical theories should only propagate one type of massless spin-2 graviton on constant curvature backgrounds. This implies the effective Newtonian constant must have the same sign as the one in GR, which means the third term in (2.15) should become negative for positive mass, that is $h^{^{\prime}}(f_{\infty})\lt0$ [42].

Summarizing, we only consider black holes with $f_{\infty}\gt0$, $h^{^{\prime}}(f_{\infty})\lt0$ and $\gamma^2\gt0$ as satisfying the requisite physical criteria.

Here we compare our method with methods previously employed in understanding multicritical behaviour. For this illustration, we consider a quadruple point in the Einstein-power-Maxwell theory comparing previous methods used to find it [18] with our proposed approach.

Power-Maxwell theory is a general form of non-linear electrodynamics minimally coupled to D = 4 Einstein gravity [58]. The action is

$$\begin{align} S = \int \mathrm{d}^4x \sqrt{-g}\left(R-2 \Lambda -\sum_{i = 1}^{N} \alpha _i \left(F^2\right)^i \right), \end{align} \tag{ B.1 }$$

with $F^2 \equiv F_{\mu \nu}F^{\mu \nu}$ and $F_{\mu \nu}\equiv \nabla_{\mu} A_\nu - \nabla _\nu A_\nu$, where R is the Ricci scalar. The α_i are dimensional coupling constants ( $[\alpha_i] = L^{2(i-1)}$) and A_µ is the U(1) Maxwell field. We recover Einstein–Maxwell theory when $\alpha_1 = 1$ and $\alpha_i (i\gt1) = 0$.

For the ansatz (1.2) it is straightforward to find a solution of the form $f(r) = 1+\sum_{i = 1}^{K} c_i r^{-i}+\frac{r^2}{l^2}$ [18]. The corresponding thermodynamic variables are

$$\begin{align} T & = \frac{1}{4\pi r_+}\Bigg[1+8\pi P r_+^2-\frac{Q^2}{r_+^2}-\frac{5b_5 Q}{2r_+^6}-\frac{3b_9 Q}{r_+^{10}} -\frac{13b_{13}Q}{4r_+^{14}}-\frac{17b_{17}Q}{5r_+^{18}}-\frac{7b_{21}Q}{5r_+^{22}}+\mathcal{O}\left(\frac{1}{r_+^{26}}\right)\Bigg], \end{align} \tag{ B.2 }$$

$$\begin{align} S & = \pi r_+^2, \quad V = \frac{4}{3}\pi r_+^3, \quad P = \frac{3}{8\pi l^2}, \end{align} \tag{ B.3 }$$

where Q is the charge parameter, and each of the other parameters represent the same physical quantities as before. The field equations imply

$$\begin{align} c_1 = - 2M \qquad c_i = \frac{4Q}{i+2}b_{i-1}, \qquad \text{for} \qquad i>1 \end{align} \tag{ B.4 }$$

where

$$\begin{align} b_1 & = Q \qquad b_5 = \frac{4}{5}Q^3\alpha_2 \qquad b_9 = \frac{4}{3}Q^5\left(4\alpha_2^2-\alpha_3\right) \end{align} \tag{ B.5 }$$

$$\begin{align} b_{13} & = \frac{32}{13}Q^7\left(24 \alpha_2^3-12 \alpha_3\alpha_2+\alpha_4\right) \end{align} \tag{ B.6 }$$

$$\begin{align} b_{17} & = \frac{80}{17}Q^9\left(1 76\alpha_2^4-132 \alpha_2^2\alpha_3+16\alpha_4\alpha_2+9\alpha_3^2-\alpha_5\right) \end{align} \tag{ B.7 }$$

$$\begin{align} b_{21} & = \frac{64}{7}Q^{11}\left(1456\alpha_2^5+234\alpha_3^2\alpha_2+208\alpha_4\alpha_2^2-24\alpha_4\alpha_3 - 1456\alpha_2^3\alpha_3-20\alpha_5\alpha_2+\alpha_6 \right) \end{align} \tag{ B.8 }$$

and we can set $\alpha_1 = 1$ without loss of generality.

The Gibbs free energy is again $G = M-TS$ but with the black hole mass now modified to be

$$\begin{align} \begin{aligned} M = \frac{4}{3} \pi P r_+^3+\frac{r_+}{2}+\frac{Q^2}{2 r_+}+\frac{b_5 Q}{4 r_+^5}+\frac{b_9 Q}{6 r_+^9} +\frac{b_{13} Q}{8 r_+^{13}}+\frac{b_{17} Q}{10 r_+^{17}}+\frac{b_{21}Q}{12r_+^{21}}+\mathcal{O}\left(\frac{1}{r_+^{25}}\right). \end{aligned} \end{align} \tag{ B.9 }$$

The approach ⁵ stated in [18] requires presetting $2N-1$ intersections $r_+^{(i)}$ to construct an N-tuple point by solving $2N-1$ equations

$$\begin{align} T\left(r_+^{\left(1\right)},P^*,\left\{b_n\right\}\right) = \cdots = T\left(r_+^{\left(2N-1\right)},P^*,\left\{b_n\right\}\right) = T^{{\,}*} \end{align} \tag{ B.10 }$$

for $T^{{\,}*}$, $P^*$, Q, b_n . This yields a numerical expression for the temperature (B.2) each time we change the values of intersections. Hence an N-tuple point (likely) emerges eventually by adjusting the $r_+^{(i)}$ such that N − 1 swallowtails of the free enthalpy appears at a single point. The feasibility of the method is lucid by applying a similar analysis in the section 4 to (B.2). As a consequence of the rule of signs, $2N-4$ is the minimum number of couplings mandatory to construct an N-tuple point.

For an illustration, we apply this approach to a quadruple point, and so four non-zero b_n 's are required. For simplicity, we choose them to be $\{b_5,b_9,b_{13},b_{17}\}$. In order to work out a desired set of b_n 's, we pick intersections $\{1,1.04,1.1,1.4,1.9,2.9,4\}$ (set A) to start with. However, this choice of $r_+^{(i)}$ is not optimal in the sense that only two swallowtails appear in the Gibbs free energy, as shown in figure 6, which is far from a quadruple point merger. After a few iterations, we arrive at a better situation displayed by the figure 7, where all three swallowtails appear. These do not get intersect a single point until all radii are tuned to $\{1,1.0424,1.1,1.4243,1.9,2.93,4\}$ (set B) as illustrated in figure 8. Radii with higher precision can be obtained if one iterates the approach for many times, but a considerably longer time is inevitably needed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now consider how things change using our method. We first define the function K as

$$\begin{align} K\left(r_+,r_0\right) = \int_{r_0}^{r_+}4\pi\left[P^*-P\left(\tilde{r}_+,T^{{\,}*},\left\{b_n\right\}\right)\right] \tilde{r}_+^2 \mathrm{d} \tilde{r} \end{align} \tag{ B.11 }$$

where (B.2) and $V = 4\pi r_+^3/3$ are applied to (4.1). Next pick $\{1,1.1,1.9,4\}$ out of the set A to be four different values of $r_+$ satisfying (3.6). After one implementation of the procedure discussed in section 4, we find

$$\begin{align} r_+^{\left(1\right)} & = 1, \quad r_+^{\left(2\right)} = 1.042\,427\,218, \quad r_+^{\left(3\right)} = 1.1,\nonumber\\ r_+^{\left(4\right)} & = 1.424\,291\,548, \quad r_+^{\left(5\right)} = 1.9, \quad r_+^{\left(6\right)} = 2.930\,025\,838, \quad r_+^{\left(7\right)} = 4,\end{align} \tag{ B.12 }$$

which reproduces the set B with considerably higher precision and notably less effort. From this perspective, our approach can be interpreted as a fine-tuning process that automatically identifies all possible proper positions of the N − 1 remaining radii in between the preset N radii such that an N-tuple point emerges. All those radii collectively constitute a workable set of $r_+^{(i)}$ which constructs an N-tuple point in the preceding method.