**1. Introduction**

In recent years, studies of black holes with scalar hair received much interest, both in the context of generalized gravity theories such as Einstein–scalar–Gauß–Bonnet (EsGB) theories [1–22] and in the context of simpler models such as Einstein–Maxwell–scalar (EMs) models [23–40]. In both cases, the way to couple scalar fields to the respective invariants, the Gauß–Bonnet and the Maxwell invariant, is crucial for the resulting types of black hole solutions and their properties. In fact, according to the choice of coupling function, distinct classes arise. In the first class, black holes with scalar hair are present but with no general relativity (GR) black holes. Examples here are dilatonic black holes, first obtained long ago [41,42]. In contrast, the second class allows for both black holes with scalar hair and GR black holes. In the latter case, we should distinguish whether the GR black holes can exhibit tachyonic instability, causing spontaneous scalarization of black holes [1–3,23], or whether the GR black holes will never succumb to such an instability [31].

An example of the latter has been studied in [38,39]. In this set of models, the coupling function *f*(Φ) = 1 + *α*Φ<sup>4</sup> has been chosen, coupling the scalar field *φ* to the Maxwell invariant with coupling constant *α*. Clearly, this coupling function allows for GR black hole solutions, since the resulting coupled set of field equations can always be satisfied with a vanishing scalar field. Thus, Reissner–Nordström (RN) black holes are solutions of this set of models with their usual domain of existence, limited by the set of extremal RN black holes. However, this set of models allows also for black holes with scalar hair. These scalarized solutions can be found in two different sets or branches: the cold branch and the hot branch [38,39], that are labeled according to their horizon temperature. For fixed coupling *α*, their domain of existence can be expressed in terms of their charge to mass ratio *q* = *Q*/*M*. The cold branch resides in the interval [*q*min(*α*), 1] and the hot branch resides in [*q*min(*α*), *q*max(*α*)], while the bald RN branch resides in [0, 1].

On the other hand, hairy black holes with non-Abelian gauge fields have been studied for a long time (see, e.g., [43–45] for reviews). The Einstein–Yang–Mills (EYM) and Einstein–Yang–Mills–Higgs

(EYMH) theories typically allow not only for black holes with non-Abelian hair but also for embedded Abelian solutions, such as RN black holes. This is thus analogous to the EMs models of the second class. Unlike most EMs models though, the EYM and EYMH models also feature globally regular solutions, solitons. In the SO(3)-EYMH case, these solitons correspond to gravitating magnetic monopoles and dyons, which can be endowed with a horizon, generating hairy black holes [46–52]. These non-Abelian solutions do not exist for arbitrary values of the coupling constant or horizon radius but possess a limited domain of existence. As one of the boundaries is approached, an interesting limiting behavior is observed: At a critical radial coordinate *rcr* = *Qcr*, the space time divides into two parts. In the exterior part *r* > *rcr*, the scalar field assumes its vacuum expectation value while the non-Abelian gauge field vanishes except for an embedded Abelian field, yielding the exterior region of an extremal RN black hole with degenerate horizon *rcr*. In the interior part *r* < *rcr*, nontrivial non-Abelian and scalar fields remain. These tend to their respective vacuum values at *rcr*. As this critical solution is approached, both parts of the space time become infinite in extent, since the metric coefficient of the radial coordinate features the double zero of a degenerate horizon.

Here, we will consider the limiting behavior of the EMs hairy black holes on the cold branch as the upper boundary of their domain of existence, *q* = 1, is approached. Since the upper limit indeed has *q* = 1 and a vanishing horizon temperature, one may expect that an extremal RN solution is approached. However, as we will show, this is only part of the truth. In fact, an analogous critical behavior is encountered, as has been known in the case of non-Abelian solutions for long. As the critical solution is approached, the space time splits into two parts, an exterior part *r* > *rcr* corresponding to the exterior region of an extremal RN black hole and an interior part *r* < *rcr* with a finite scalar field that vanishes at *rcr*.

In Section 2, we present the EMs theory and the equations of motion. We then recall the Ansatz for spherically symmetric black hole solutions and the expansions at the horizon and at infinity. The properties of the black hole solutions are recalled in Section 3. Next, we consider the critical behavior for fixed coupling constant *α* and then address the *α*-dependence of the properties of the critical solutions in the same section. In Section 4, we show that the excited solutions exhibit an analogous critical behavior. We conclude in Section 5.

#### **2. EMs Theory**

We consider EMs theory described by the action

$$\mathcal{S} = \int d^4 \mathbf{x} \sqrt{-\mathbf{g}} \left[ R - 2 \partial\_{\mu} \Phi \partial^{\mu} \Phi - f(\Phi) F\_{\mu \nu} F^{\mu \nu} \right],\tag{1}$$

with the Ricci scalar *R*, the real scalar field Φ, the Maxwell field strength tensor *<sup>F</sup>μν* and the coupling function *f*(Φ). We use units so that *c* = *G* = 1.

In this paper, we will assume a quartic dependence:

$$f(\Phi) = 1 + a\Phi^4. \tag{2}$$

We will focus on positive values of the coupling constant *α*, meaning Φ = 0 is the global minimum of the coupling function. This type of coupling is representative of the type IIB (scalarised-disconnected-type following the nomenclature introduced in [31]), since it allows for the existence of both bald (RN) and hairy (scalarised) black holes. However, RN solutions are not unstable with respect to scalar perturbations [38,39].

*Symmetry* **2020**, *12*, 2057

The Einstein–, Maxwell– and scalar-field equations follow from the variational principle and read as follows:

$$R\_{\mu\nu} - \frac{1}{2}\mathfrak{g}\_{\mu\nu}R = T^{\Phi}\_{\mu\nu} + T^{EM}\_{\mu\nu} \tag{3}$$

$$\nabla\_{\mu} (\sqrt{-\overline{\emptyset}} f(\Phi) F^{\mu \nu}) = 0 \,, \tag{4}$$

$$\frac{1}{\sqrt{-g}}\partial\_{\mu}(\sqrt{-g}g^{\mu\nu}\partial\_{\nu}\Phi) = \dot{f}(\Phi)F\_{\mu\nu}F^{\mu\nu},\tag{5}$$

with ˙ *f*(Φ) = *d f*(Φ)/*d*Φ, electromagnetic stress-energy tensor

$$T^{EM}\_{\mu\nu} \equiv 2f(\Phi) \left( F\_{\mu a} F\_{\nu}{}^{a} - \frac{1}{4} g\_{\mu\nu} F^2 \right) \tag{6}$$

and scalar stress-energy tensor

$$T^{\Phi}\_{\mu\nu} \equiv \frac{1}{2} \partial\_{\mu} \Phi \partial\_{\nu} \Phi - \frac{1}{2} g\_{\mu\nu} \frac{1}{2} (\partial\_{\alpha} \Phi)^2 \,. \tag{7}$$

To study static spherically symmetric black holes, we employ the following parametrization of the metric:

$$ds^2 = -N(r)e^{-2\delta}dt^2 + \frac{dr^2}{N(r)} + r^2(d\theta^2 + \sin^2\theta d\phi^2) \,\,,\tag{8}$$

with the metric functions *<sup>N</sup>*(*r*) and *<sup>δ</sup>*(*r*). However, for some expressions, it is useful to define the mass function *m*(*r*) as *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* and the metric function as *g*(*r*) = *<sup>N</sup>*(*r*)*e*<sup>−</sup>2*<sup>δ</sup>*(*r*) = <sup>−</sup>*gtt*(*r*).

To obtain black holes with electric charge and, in the scalarized case, also scalar charge, we parametrize the gauge potential and the scalar field by

$$A\_{\!\!\!\!=} = \left(A\_{\!\!\!=}(r), 0, 0, 0\right), \quad \Phi \quad = \quad \Phi(r) \,. \tag{9}$$

The above Ansatz for a static and spherically symmetric configuration is substituted into the equations of motion (3)–(5). This leads to the following set of coupled ordinary differential equations (ODEs) for the unknown functions of the metric and the fields:

$$m' = \frac{r^2 N \Phi'^2}{2} + \frac{Q^2}{2r^2 f(\Phi)},\\ \delta' + r \Phi'^2 = 0,\\ A'\_t = -\frac{Q \varepsilon^{-\delta}}{f(\Phi) r^2},\tag{10}$$

$$\Phi''(r) + \frac{1+N}{rN} \Phi' - \frac{Q^2}{r^3 N f(\Phi)} \left(\Phi' - \frac{f(\Phi)}{2rf(\Phi)}\right) = 0 \,,\tag{11}$$

where a prime denotes a derivative with respect to the radial coordinate and *Q* is the electric charge of the black holes.

To address the vicinity of the black hole horizon, we perform a power series expansion of the functions of the metric and the fields in *r* − *rH* at the horizon, where we denote the horizon radius by *rH* and the horizon values of the functions by the subscript *H*:

$$m(r) = \frac{r\_H}{2} + m\_1(r - r\_H) + \cdots \quad , \qquad \qquad \delta(r) = \delta\_H - \Phi\_1^2 r\_H (r - r\_H) + \cdots \quad , \tag{12}$$

$$A\_{l}(r) = \Psi\_{H} - \frac{e^{-\delta\_{H}}Q}{r\_{H}^{2}f(\Phi\_{H})}(r - r\_{H}) + \cdots \ , \qquad \qquad \Phi(r) = \Phi\_{H} + \Phi\_{1}(r - r\_{H}) + \cdots \ , \tag{13}$$

with

$$m\_1 = \frac{Q^2}{2r\_H^2 f(\Phi\_H)}, \quad \Phi\_1 = \frac{Q^2 \dot{f}(\Phi\_H)}{2r\_H f(\Phi\_H) \left[Q^2 - r\_H^2 f(\Phi\_H)\right]} \tag{14}$$

*Symmetry* **2020**, *12*, 2057

Some important properties of the horizon that we will use are the horizon area *AH* = <sup>4</sup>*πr*2*H* and the horizon temperature *TH* = *<sup>N</sup>* (*rH*)*e*<sup>−</sup>*<sup>δ</sup>*(*rH*)/4*π* (which is determined by the surface gravity at the horizon in the usual way).

A power series expansion in 1/*r* at infinity of the functions of the metric and the fields allows us to read off the global charges of the black holes:

$$m(r) = M - \frac{Q^2 + Q\_s^2}{2r} + \cdots, \qquad \delta(r) = \frac{Q\_s^2}{2r^2} + \cdots,\tag{15}$$

$$A\_t(r) = -\frac{\mathcal{Q}}{r} + \cdots, \qquad \qquad \qquad \Phi(r) = \frac{\mathcal{Q}\_s}{r} + \frac{MQ\_s}{r^2} + \cdots, \tag{16}$$

with Arnowitt-Deser-Misner mass *M* and scalar charge *Qs*, both of them also defined in the standard way.

#### **3. Limit of Cold Black Holes**

#### *3.1. Branches of Black Holes*

We now briefly recall the properties of static spherically symmetric electrically charged black hole solutions with quartic coupling function (2). The black holes of the RN branch are given by

$$\delta(r) = 0 \text{ , } m(r) = M - \frac{Q^2}{2r} \text{ , } A\_l(r) = -\frac{Q}{r} \text{ , } \Phi(r) = 0 \text{ .} \tag{17}$$

The scalarized black hole solutions are obtained numerically [38]. We solve the field equations subject to the boundary conditions that follow from the above expansions at the horizon and at infinity, with input parameters *α*, *rH* and *Q*. We employ the professional solver COLSYS [53], which is based on a collocation method for boundary-value ODEs and on a damped Newton method of quasi-linearization. Since this solver includes an adaptive mesh selection procedure, it is very suitable for the problem at hand, where high accuracy is needed in a very small interval close to the black hole horizon. Consequently, the grid is successively refined until the required accuracy is reached, typically 10−16.

The solutions are characterized by a set of dimensionless quantities: the charge to mass ratio *q*, the reduced horizon area *aH* and the reduced horizon temperature *tH*, for which

$$q \equiv \frac{\mathcal{Q}}{M}, \qquad a\_H \equiv \frac{A\_H}{16\pi M^2} = \frac{r\_H^2}{4M^2}, \qquad t\_H \equiv 8\pi M T\_H = 2MN'(r\_H)e^{-\delta(r\_H)} \quad . \tag{18}$$

We illustrate the domain of existence of the black holes in Figure 1a, where we show the mass to charge ratio *q* versus the coupling constant *α*. The horizontal black line *q* = 1 represents the set of extremal RN black holes and, thus, the upper boundary for RN black holes. At the same time, that line represents the set of critical scalarized solutions forming the upper boundary of the cold black holes, which reside in the lower green area. The solid green line marks the bifurcation line *q*min(*α*) separating the cold and hot black holes. The hot black holes then extend from this bifurcation line to the dashed red critical line *q*max(*α*), i.e., they fill the whole shaded region.

In Figure 1b, we exhibit the reduced area *aH* and reduced temperature *tH* (inset) versus the charge to mass ratio *q* of the black hole solutions for the particular coupling *α* = 200 [38]. The RN branch is shown in solid black, the cold branch is shown in dotted blue and the hot branch is shown in dashed red. Along the cold branch, the mass to charge ratio *q* decreases, while the reduced area *aH* and temperature *tH* increase. At the minimal value *q*min, the cold branch bifurcates with the hot branch. Along the hot branch, *aH* decreases again, while *tH* increases with increasing *q*.

From the figure, it seems that the cold branch starts from an extremal RN black hole. Clearly, the charge to mass ratio at its endpoint agrees with the ratio *q* = 1 of an extremal RN black hole and the horizon temperature vanishes at its endpoint, *TH* = 0. Looking at the horizon area (Figure 1b) and further properties, one is indeed tempted to conclude that the cold branch emerges from an extremal

RN black hole. However, as we will demonstrate in the following, this is only partially true. In fact, the endpoint of the cold branch is a more intriguing configuration.

**Figure 1.** (**a**) Phase diagram of scalarized black holes: mass to charge ratio *q* vs. coupling constant *α*. (**b**) Reduced horizon area *aH* (inset: reduced temperature *tH*) vs. *q* for *α* = 200: cold branch (dotted blue), hot branch (dashed red) and RN branch (solid black).

#### *3.2. Approach to the Critical Solution for α* = 200

To gain an understanding of the limiting configurations, we now inspect a family of cold black holes with fixed coupling constant *α* = 200, fixed horizon radius *rH* = 2 and increasing charge *Q* in order to reach the critical solution with *q* = 1 at some *Qcr*. If the family of cold black holes simply approached an extremal RN black hole, this extremal black hole would have its degenerate horizon at *rH* = 2, and it would satisfy *rH* = *Q* = *M*.

In Figure 2, we exhibit a set of these cold black hole solutions as they approach the endpoint of the cold branch. Notably, all the solutions exhibited here possess a charge *Q* ≥ 2, i.e., they range between the extremal RN limit of *Q* = 2 and a critical value *Qcr* = 2.0016701 (+*<sup>O</sup>*(<sup>10</sup>−<sup>8</sup>)). Thus, they slightly exceed the extremal RN limit for fixed horizon radius *rH* = 2. The figure shows the metric functions *g*(*r*) (a) and *<sup>N</sup>*(*r*) (b), the electromagnetic function *At*(*r*) (c) and the scalar field function <sup>Φ</sup>(*r*) (d) versus the compactified radial coordinate *x* = 1 − *rHr* in a logarithmic scale. Note that *r* = *rH* then corresponds to *x* = 0 and *r* = ∞ corresponds to *x* = 1.

Whereas the metric and electromagnetic functions appear smooth at first glance, as seen in the inlets of Figures 2a,b and in Figure 2c, the scalar field function (Figure 2d) immediately reveals a critical behavior: As *Qcr* is approached, the scalar field function assumes a finite limiting value <sup>Φ</sup>(*rH*) at the imposed horizon *rH* = 2. However, the function <sup>Φ</sup>(*r*) decreases more and more steeply as it approaches zero, with its boundary value at infinity, <sup>Φ</sup>(∞) = 0. In fact, in the limit *Q* → *Qcr*, it tends to zero already at a critical value of the radial coordinate *r* = *rcr*, which coincides with the value of the critical charge, *rcr* = *Qcr*. In the limit, therefore, the solution features a finite scalar field in the interior *r* < *rcr*, whereas the scalar field vanishes identically in the exterior *r* > *rcr*.

Since the critical exterior solution is a pure electrovacuum solution, starting at *rcr* and possessing electric charge *Qcr*, this suggests that the exterior critical solution is described by an extremal RN black hole. However, instead of carrying charge *Q* = 2, it carries the critical charge *Qcr*. Comparing the *Qcr* numerical solution with a *Qcr* extremal RN black hole shows that this conclusion holds true. In the interior, however, not only is the scalar field function finite but also all the functions differ from this *Qcr* extremal RN black hole, as they must in order to satisfy the imposed boundary conditions at *rH* = 2.

**Figure 2.** Approach to the critical solution for *α* = 200: (**a**) metric function *g*(*r*) = *<sup>N</sup>*(*r*)*e*<sup>−</sup>2*<sup>δ</sup>*(*r*), (**b**) metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* , (**c**) electromagnetic function *At*(*r*) and (**d**) scalar function <sup>Φ</sup>(*r*) vs. the compactified radial coordinate *x* = 1 − *rHr* . The insets highlight the vicinity of the critical radius *r* = *rcr*.

We now inspect the behavior of the functions in the interior in more detail as *Qcr* is approached, starting with the electromagnetic function *At*. As seen in Figure 2c, also in the interior, a limiting critical solution is reached. At the critical radius *rcr*, the electromagnetic function *At* of this critical solution assumes the value *At*(*rcr*) = −1. In fact, in the full interior *r* < *rcr*, it assumes this value *At*(*r* ≤ *rcr*) = −1, as seen in the inset of the figure. Therefore, there is no electric field in the interior region.

To reveal the critical behavior of the metric functions, we need to consider double logarithmic plots, as exhibited in Figure 2a,b. The extremal RN with charge *Qcr* would have a double zero at *rcr*, for both functions *g*(*x*) and *<sup>N</sup>*(*r*). Indeed, we observe a very sharp drop at *rcr* as the limiting solution is approached for both metric functions *g*(*r*) and *<sup>N</sup>*(*r*), confirming our interpretation of the exterior solution. However, in the interior, it becomes apparent that we have not ye<sup>t</sup> fully reached the critical solution but are only very close to it.

In the interior, both functions *g*(*r*) and *<sup>N</sup>*(*r*) differ distinctly. The function *<sup>N</sup>*(*r*) tends to a finite limiting solution in the interior, except at *rH* = 2, where the boundary conditions force it to vanish. In contrast to *<sup>N</sup>*(*r*), the function *g*(*r*) approaches zero in the limit. (With every further digit determined at the critical value *Qcr*, the function *g*(*r*) assumes smaller values in the interior.) Recalling that the function *g*(*r*) has been decomposed into the factors *<sup>N</sup>*(*r*) and exp(−2*<sup>δ</sup>*(*r*)), we conclude that it is the

function *<sup>δ</sup>*(*r*) which causes *g*(*r*) to vanish in the interior in the limit, since *<sup>N</sup>*(*r*) has a finite limit. It is instructive now to look again at the horizon temperature *TH*. Having observed above (e.g., in the inset of Figure 1b) that *TH* → 0 on the cold branch in the limit, one may expect that the reason was that a degenerate horizon arose at *rH* = 2 and therefore the derivative *<sup>N</sup>* (*r*) vanished there. However, now, we see that there is no degenerate horizon at *rH* = 2 in the limit. Instead, *TH* vanishes because of the factor exp(−*<sup>δ</sup>*(*r*)) in Equation (18).

We have thus obtained the following scenario: In the limit *Q* → *Qcr*, the space time splits into an exterior and an interior part. The exterior is described by an extremal RN black hole; the interior has a finite scalar field but no electric field. Since at the critical radius *rcr*, a double zero is approached, as featured by a degenerate horizon, the radial distance *l*(*r*) to and from *rcr* increases as the critical solution is approached. In the limit, both parts of the space time become infinite in extent.

To demonstrate this effect for the interior part, we consider in Figure 3 the radial distance *l*(*r*), defined via

$$l(r) = \int\_{r\_H}^{r} \frac{dr}{\sqrt{N(r)}} \quad . \tag{19}$$

The dependence of the radial coordinate *r*(*l*) on the radial distance is illustrated in Figure 3a for the approach to the critical solution. For *Qcr*, an infinite radial distance is reached at the finite value of the radial coordinate *rcr*; thus, an infinite throat is formed. For *Q* → *Qcr*, the metric function *<sup>N</sup>*(*r*) (Figure 3b), the electromagnetic function *At*(*r*) (Figure 3c) and the scalar function <sup>Φ</sup>(*r*) (Figure 3d) are also shown versus the radial distance *l*(*r*). The metric function *<sup>N</sup>*(*r*) again highlights the formation of an infinite throat in the limit, while the electromagnetic function *At*(*r*) approaches a constant value in the interior region. The scalar function <sup>Φ</sup>(*r*), on the other hand, demonstrates that, when considered as a function of the radial distance *l* instead of the radial coordinate *r*, there is remarkably little dependence on the value of the charge *Q* during the approach *Q* → *Qcr*. An analogous observation was made for the matter functions of the magnetic monopoles during their approach to their respective critical solution [46].

**Figure 3.** *Cont.*

**Figure 3.** Approach to the critical solution for *α* = 200: (**a**) radial coordinate *r*, (**b**) metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* , (**c**) electromagnetic function *At*(*r*) and (**d**) scalar function <sup>Φ</sup>(*r*) vs. the radial distance *l*(*r*) (Equation (19)).

#### *3.3. α-Dependence of the Critical Solution*

We now demonstrate that the critical scenario described above holds for a large range of couplings *α*, showing that the scenario is rather generic. We exhibit the critical solutions in Figure 4, where we show the metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* (Figure 4a) and the scalar function <sup>Φ</sup>(*r*) (Figure 4b) for a set of couplings *α* in the interval [3.2, <sup>200</sup>].

**Figure 4.** Critical solutions for a set of couplings *α*: (**a**) metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* and (**b**) scalar function <sup>Φ</sup>(*r*) vs. the compactified radial coordinate *x* = 1 − *rHr* . Note that the key applies to both figures.

The figure shows that, with decreasing *α*, the critical radius *rcr* increases. We recall that it coincides with the critical charge, *rcr* = *Qcr*. At the same time, the scalar field assumes larger values in the interior. We highlight the *α*-dependence of the critical charge *Qcr* and of the horizon value of the scalar field <sup>Φ</sup>(*rH*) in Figure 5a,b, respectively. We note the steep increase of the critical charge *Qcr* for small *α*, while *Qcr* → *rH* for large *α*. This dependence can be well described by the simple relation

$$\frac{Q\_{cr}}{r\_H} - 1 = \frac{1}{4\sqrt{2}\,\alpha'}\tag{20}$$

as demonstrated in the figure. The horizon value of the scalar field <sup>Φ</sup>(*rH*) satisfies the even simpler relation

$$
\Phi(r\_H) = \frac{1}{\sqrt{\alpha}},
\tag{21}
$$

as seen in the figure as well. In Figure 5c,d, we show that, for the interior critical solution, the derivatives of the metric function *m*(*r*) and the scalar function <sup>Φ</sup>(*r*) at the horizon precisely respect the expansion at the horizon (Equations (12) and (13), respectively).

**Figure 5.** Properties of critical solutions: (**a**) *Qcr*(*α*) (solid red) and limit for *α* → ∞ according to Equation (20) (dashed green), (**b**) <sup>Φ</sup>(*rH*)(*α*) (solid red) and relation from Equation (21) (dashed green), (**c**) <sup>Φ</sup> (*rH*)(*α*) (solid red) and expansion coefficient <sup>Φ</sup>1(*rH*)(*α*) from Equation (14) (dashed green), and (**d**) *m* (*rH*)(*α*) (solid red) and expansion coefficient *<sup>m</sup>*1(*rH*)(*α*) from Equation (14) (dashed green).

#### **4. Excited Solutions**

Until now, we have focused the discussion on the fundamental branch of scalarized black holes, with scalar field functions that possess no nodes. However, in certain cases, the model allows for the existence of excited solutions, meaning solutions with nodes in the scalar field functions. Let us discuss them briefly here.

These solutions present very similar properties to the *n* = 0 solutions, with a hot/cold branch structure, the hot one bifurcating from the cold one. The cold branch also reaches a critical end-point with *q* = 1. However, an important difference is that all of these excited branches are radially unstable, as previously discussed in [38,39]. Also, the number of excited branches depends on the coupling constant *α*, i.e., the larger the value of *α*, the more excited branches exist.

As an example, in Figure 6a, we show again the reduced area *aH* as a function of the reduced charge *q* for *α* = 200. This is the same as Figure 1b, but now, we include the *n* = 1 solutions in solid green. As we can see, these excited solutions could also be distinguished in two different branches: one branch (cold) would extend from *q* = 1 to a minimum *q*min < 1; the second branch (hot) would extend from this *q*min up to a certain *q*max > 1. It turns out that the interval (*q*min, *q*max), where *n* = 1 solutions exist, is smaller than the domain interval of the *n* = 0 solutions. On the other hand, for this particular value of the coupling constant (*α* = 200), only *n* = 1 excited solutions exist. However, sufficiently large values of *α* allow for additional branches with *n* > 1 excited solutions appear.

**Figure 6.** (**a**) Reduced horizon area *aH* vs. *q* for *α* = 200: In addition to the RN branch (solid black), the (*n* = 0) cold branch (finely dotted blue) and the (*n* = 0) hot branch (coarsely dotted red), the *n* = 1 excited solutions (solid green) are shown. (**b**) A similar figure for the scalar field at the horizon Φ*H* vs. reduced temperature *tH*.

In Figure 6b, we show the scalar field at the horizon as a function of the reduced temperature *tH*, also for *α* = 200. As we approach the *q* = 1 limit along the cold branch, the temperature goes to zero but the value of the scalar field remains constant. (This is similar to what happens for the *n* = 0 solutions we have discussed in the previous sections, as shown in this figure.) If we compare solutions with the same *th*, the larger the excitation number, the larger the value of the scalar field at the horizon.

As already said, the critical behavior is also present in these excited solutions. When fixing *rH* = 2, the critical solutions possess a value of the electric charge *Q* = *Qcr*(*n*) > 2 (this value depends on *n*). The limit results in a set of critical solutions with split space time, similar to the *n* = 0 solutions we have been discussing in the previous sections. The exterior part is the extremal RN solution. On the interior part, however, we have a nontrivial scalar field for which the number of nodes can be labeled by an integer number *n*.

In Figure 7, some profiles for *n* = 1 critical solutions with different values of *α* are shown as an example. In particular, we show the critical solutions for *α* = 200 with *Qcr* = 2.02407591 (solid red), for *α* = 400 with *Qcr* = 2.01180482 (dashed green) and for *α* = 800 with *Qcr* = 2.00584391 (dotted blue).

In Figure 7a, we show the metric function *<sup>N</sup>*(*r*) versus *x* = 1 − *rH*/*<sup>r</sup>*. It demonstrates that the behavior of these *n* = 1 solutions is very similar to the one presented in Figure 4a, which corresponds to *n* = 0 solutions. Essentially, the metric function *<sup>N</sup>*(*r*) develops a zero at some critical value *r* = *rcr*. For *r* > *rcr*, the solution is extremal RN, but in the interior region, it is a nontrivial solution with a scalar field. The main difference now, as shown more clearly in the inset, is that the function develops a minimum at a certain point *r* < *rcr*.

**Figure 7.** Critical solutions with *n* = 1 for a set of couplings *α*: (**a**) metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* and (**b**) scalar function <sup>Φ</sup>(*r*) vs. the compactified radial coordinate *x* = 1 − *rHr*.

(**b**)

Figure 7b exhibits the scalar field Φ versus *x*, with the inset showing a zoomed region around the critical value of the radial coordinate *rcr*. This demonstrates that the profile of the scalar field is quite different from the one in Figure 4b. Now the scalar field function not only vanishes at *rcr* but also has a node at some intermediate point.

Solutions with more nodes also follow this pattern. As an example, we depict in Figure 8 similar plots for *n* = 2 solutions. These correspond to critical solutions for *α* = 640 with *Qcr* = 2.0362514 (solid red), for *α* = 832 with *Qcr* = 2.0274443 (dashed green) and for *α* = 960 with *Qcr* = 2.02361315 (dotted blue). Figure 8 (left) exhibits how the *N* function develops a minimum a bit below *rcr*, and in Figure 8 (right), we demonstrate that the scalar field function has two nodes below *rcr*.

**Figure 8.** Critical solutions with *n* = 2 for a set of couplings *α*: (**a**) metric function *<sup>N</sup>*(*r*) = 1 − <sup>2</sup>*m*(*r*) *r* and (**b**) scalar function <sup>Φ</sup>(*r*) vs. the compactified radial coordinate *x* = 1 − *rHr*.

## **5. Conclusions**

(**a**)

We have considered black holes in EMs models with a quartic coupling function *f*(Φ) = 1 + *α*<sup>Φ</sup>4, that allows for RN black holes as well as cold and hot scalarized black holes [38]. The models are parameterized by the coupling constant *α* and exhibit generic features. In particular, the RN black holes never become unstable to grow scalar hair [38,39]. Still, the cold branch follows the RN branch to a large extent and starts at *q* = *Q*/*M* = 1, the endpoint of the RN branch.

In particular, we have investigated the approach of the cold branch to this critical point *q* = 1 in detail for these EMs models. For that purpose, we have fixed the horizon radius *rH* of the black holes and then increased the electromagnetic charge *Q* until, for some fixed coupling *α*, a critical solution with *q* = 1 has been reached. Whereas an extremal RN black hole with horizon radius *rH* satisfies *q* = 1 with *rH* = *Q* = *M*, the critical solution on the cold branch satisfies *q* = 1 with *rcr* = *Qcr* = *Mcr* and *rH* < *rcr*. Consequently, the critical solution cannot simply be an extremal RN black hole.

Inspection of the functions of the critical solution then revealed that the critical solution splits the space time into two parts: an exterior part which indeed corresponds to an extremal RN black hole solution but with *rcr* = *Qcr* = *Mcr*, and an interior part, which has a finite scalar field that vanishes at *rcr* and a vanishing electromagnetic field. Since the radial metric function develops a double zero at *rcr*, both parts of the space time become infinite in size as the space time splits.

By studying the critical solution for many values of the coupling *α*, we have shown that this observed phenomenon is generic for these EMs models. In fact, the critical charge *Qcr* possesses a rather simple *α*-dependence, with (*Qcr*/*rH* − 1) being inversely proportional to *α*, while the horizon value of the scalar field Φ *H* is inversely proportional the square root of *α*. At the same time, the functions of the critical solution satisfy the horizon expansion at *rH* as well as the expansion at infinity, the latter of course with vanishing scalar charge.

In addition, we have discussed the case of the excited solutions, with nodes in the scalar field functions. These solutions posses a similar (hot/cold) branch structure, with a critical limit that splits the space time just like in the nodeless case. The excited solutions can be labeled by an integer number *n*, counting the number of nodes of the scalar in the interior part of the space time.

It is interesting that the critical phenomenon which has previously been encountered for non-Abelian magnetic monopoles and their associated hairy black holes as well as for similar non-Abelian solutions arises also in these EMs models. Note a somewhat similar observation in the recent investigation of strong gravity effects of charged *Q*-clouds and inflating black holes [54] (see also [55] in this issue). For non-Abelian monopoles, the interior solution has nontrivial non-Abelian gauge fields and Higgs field, whereas the exterior solution is simply an embedded RN black hole with magnetic charge. Not surprisingly, for solutions with both electric and magnetic charge, this phenomenon persists. Even pure non-Abelian solutions (without a Higgs field) were shown to exhibit a somewhat analogous phenomenon, which emerges in the limit of infinite node number (see e.g., [43]). It would thus be interesting to find general criteria to identify models that allow for this type of phenomenon where the space time splits into two infinitely extended parts, with the interior part containing nontrivial fields that vanish identically in the exterior part, where a simple extremal black hole solution emerges.

**Author Contributions:** Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing—original draft preparation, writing—review and editing, visualization, supervision, project administration, funding acquisition: S.K., J.K. and J.L.B.-S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by: DFG Research Training Group 1620 *Models of Gravity*; DFG Project No. BL 1553; COST Action CA15117; COST Action CA16104; FCT project PTDC/FISOUT/28407/2017; FCT project PTDC/FIS-AST/3041/2020

**Acknowledgments:** J.L.B.-S., S.K. and J.K. gratefully acknowledge support by the DFG Research Training Group 1620 *Models of Gravity*. J.L.B.-S. would like to acknowledge support from the DFG Project No. BL 1553 and the FCT projects PTDC/FISOUT/28407/2017 and PTDC/FIS-AST/3041/2020. We also would like to acknowledge networking support by the COST Actions CA15117 and CA16104.

**Conflicts of Interest:** The authors declare no conflict of interest.
