Higher-order Skyrme hair of black holes

Higher-order derivative terms are considered as replacement for the Skyrme term in an Einstein-Skyrme-like model in order to pinpoint which properties are necessary for a black hole to possess stable static scalar hair. We find two new models able to support stable black hole hair in the limit of the Skyrme term being turned off. They contain 8 and 12 derivatives, respectively, and are roughly the Skyrme-term squared and the so-called BPS-Skyrme-term squared. In the twelfth-order model we find that the lower branches, which are normally unstable, become stable in the limit where the Skyrme term is turned off. We check this claim with a linear stability analysis. Finally, we find for a certain range of the gravitational coupling and horizon radius, that the twelfth-order model contains 4 solutions as opposed to 2. More surprisingly, the lowest part of the would-be unstable branch turns out to be the stable one of the 4 solutions.


Introduction
Black holes pose interesting questions about fundamental physics, such as quantum information and whether quantum theory is truly unitary at all scales. A simpler question is whether a black hole (BH) is characterized just by global charges at infinity or it has the ability to possess hair. That is, can the BH stably couple to fields of the standard model yielding a non-vacuum configuration that surrounds said BH. An attempt at answering this question in the negatory is made by the weak no-hair conjecture: for a scalar field coupled to Einstein gravity, all BHs fall in the Kerr-Newman family of solutions, i.e. being characterized by their mass, charge and spin.
The first stable counterexample was provided in the Einstein-Skyrme model [1][2][3][4][5][6][7][8]. In some sense a counterexample in a specific scalar field theory is not too exciting. However, the Skyrme model [9,10] is believed to be an effective field theory describing QCD at low energies, at least in the limit of a large number of colors [11,12], where baryons are identified with solitons -called the Skyrmions. This means that at low energies, the presence of effective operators induced by QCD would potentially be able to give rise to stable BH hair.
is positive (from f h π); hence no solutions exist. Thus the instabilities present in the construction are immediately seen by the gravitational background. This obstacle lead us to construct the above-mentioned minimal Lagrangians.
Let us briefly mention some further activities in the field of Skyrme-type BH hair. The BH Skyrme hair was extended from the Schwarzschild case to spacetimes which are asymptotically anti-de Sitter (AdS) [32][33][34] and de Sitter (dS) [35]. In particular, Ref. [34] extended the result of Refs. [14,15] to AdS spacetime, i.e. that the sextic BPS-Skyrme term is not able to support stable BH hair.
Gravitating Skyrmions have also been considered in the literature, see e.g. [4,6,7,36,37]. Collective quantization of their zero modes have also been studied [38]. The spherically symmetric case of the gravitating Skyrmions has been extended to a higher charged axially symmetric case [39,40] and it was again quantized [41]. Axially symmetric BH hair was also constructed in Ref. [40]. Spinning gravitating Skyrmions were considered [42,43] and it was further found that the BPS-Skyrme model does not possess spinning Skyrmions [43] -neither in a curved nor in a flat background. More exotic configurations like the gravitating axisymmetric sphalerons in the Einstein-Skyrme model have also been studied [44]. A lower-dimensional example was considered, where exactly solvable gravitating baby-Skyrmions were found in 2 + 1 dimensions [45], see Ref. [46] for a gauged version thereof. In 5 dimensions, a generalization of the Skyrme model to O(5) also possesses solitons and in particular, solitonic hair of BHs [47] with and without spin. An exotic study in this direction involves nonstandard boundary conditions for the metric, identifies 3-space with RP 3 without the point at infinity and contemplates a π 2 -valued Skyrmion that can give rise to a negative gravitational mass, thus antigravity [48]. Neutron stars and in particular the equation of state, which is of crucial importance in that subject, have been studied in the BPS-Skyrme model coupled to gravity [49,50]. Due to the simplicity of the model, the solutions could be compared to the mean field approximation [50], which interestingly showed some deviations. Unfortunately, since neutron stars are expected to possess some amount of spin, which is acquired during a star's core collapse, the BPS-Skyrme model cannot quite be a good approximation due to the instability found in Ref. [43]. In some limit it may give reasonable answers, nevertheless.
Another interesting idea was proposed by Dvali and Gußmann, stating that the baryon number may actually be conserved by BHs and become Skyrmionic hair once swallowed by the BH [51,52]. They also proposed that there may be ways for an observer outside the horizon to measure the number of swallowed baryons, although they rely on certain assumptions.
In this paper, we will focus on what terms are able to sustain stable static hair on a Schwarzschild-type BH. We will only consider the above mentioned minimal Lagrangians as components in the full Lagrangian and leave other possibilities for future work. Specifically, the new models that we consider have a kinetic term, the Skyrme term as well as a 2n-th order derivative term, with n = 4, 5, 6. For brevity, we shall call them the 2 + 4 + 2n models. The idea is then to slowly turn off the Skyrme term, the limit will be denoted the 2 + 2n model. We find that only the 2 + 8 and the 2 + 12 models are stable. In the case of the 2 + 8 model, there is a one-parameter family of models (which we will denote γ ∈ [0, 1]) and only at the endpoint (γ = 0), BH hair solutions cease to exist. The 2 + 8 model (with γ = 1 3 ) can be thought of as the kinetic term and the Skyrme-term squared, while the 2 + 12 model is the kinetic term and the BPS-Skyrme-term squared.
As known in the literature, usually for the Skyrme-type BH hair, there are two branches of solutions; one upper branch (in the value of the profile function at the horizon) and one lower branch. The lower branch in the standard Skyrme model is found to be unstable; it has a higher Arnowitt-Deser-Misner (ADM) mass and this is backed up by a linear perturbation analysis, which shows that the lower branch has one negative eigenvalue in the perturbation spectrum [4,5,8]. In the 2 + 12 model, we find a new behavior, where the ADM mass of the lower branch crosses over that of the upper branch and thus becomes the stable one. To this end, we carry out a linear perturbation analysis to check the eigenvalue spectrum, and indeed the lower branches have regions where they contain only positive eigenvalues. Further studies in this direction, however, is needed in order to determine full nonlinear stability. Finally, in the 2 + 12 model, we find a range of the gravitational coupling where not 2, but 4 solutions exist. That is, the lower branch ceases to be single valued in the gravitational coupling. The surprising result is that the lower part of the lower branch -farthest away from the upper branch -turns out to be the one with the lowest ADM mass.
The paper is organized as follows. In the next section, we will set up the notation of the higher-derivative terms and construct the component Lagrangians that we will use for the model mentioned above. Sec. 3 then couples our generic model Lagrangian with Einstein gravity and in Sec. 4 the explicit Einstein equations, equations of motion and boundary conditions are written down for each model and finally, the numerical results will be presented. Linear stability of the 2 + 4 + 12 model will then be analyzed in Sec. 5. Due to the surprising spectrum of eigenvalues of the perturbations, we check the Skyrme model limit in Sec. 6. The models with BH hair surviving without the presence of the Skyrme term are then studied as functions of the gravitational coupling in Sec. 7. Finally, Sec. 8 concludes with a summary and discussion of our results, as well as an outlook on future work.

Lagrangian components of a higher-order Skyrme model
The class of models we will study in this paper is built out of several component Lagrangians with different numbers of derivatives. In this section, we will set up the framework for the component Lagrangians and their corresponding energy-momentum tensors. The 2n-th order Lagrangians that we consider in this paper are those proposed in Ref. [28] and were termed minimal Lagrangians. They are minimal in the sense that they contain the smallest possible number of derivatives in the i-th direction. There are of course many more possibilities for terms with 2n derivatives, which we will not consider here.

7)
C 2 ≡ − 2 + 1 2 , (2.8) as well as the quantities 14) Here, r = 1, . . . , 6, n µ ≡ ∂ µ n, and the modulo function in the first index, p + 1|r (meaning p + 1 mod r), simply ensures that the index µ r+1 is just µ 1 . The coefficients c 2n ≥ 0 in the Lagrangians in Eqs. (2.1)-(2.6) all have to be positive semi-definite, whereas the coefficients a can take any real values. In fact, we will show shortly that the coefficients a are irrelevant for the Lagrangian formulation of the theory. The scalar fields n(x) = (n 0 , n 1 , n 2 , n 3 ) are related to the chiral Lagrangian field U (x) ∈ SU(2) as with τ a being the Pauli matrices and a = 1, 2, 3 is summed over. The field U transforms as U → g L U g † R with g L,R ∈ SU(2) L,R and hence the symmetry group is SU(2) L ×SU(2) R , which is spontaneously broken to the diagonal SU(2) L+R (or it is also simultaneously broken explicitly by a pion mass term, which however we will not consider in this paper).
The first three Lagrangians (Eqs. (2.1)-(2.3)) are well known. The first, L 2 , is the standard kinetic term, the second, L 4 , is the Skyrme term and the third, L 6 , is the BPS-Skyrme term [18]. The physical interpretation of the three terms is that the kinetic or Dirichlet term accounts for the kinetic energy, the Skyrme term measures the curvature on the O(4) target space and finally, the BPS-Skyrme term acts like a perfect fluid term [53].
We will now show that I 4,5,6 vanish identically. In order to do this, we will utilize the eigenvalues of the four-dimensional strain tensor, defined by using the left invariant form for which the invariants (2.16) simply read Here, λ 1,2,3 are the eigenvalues of the strain tensor and V is the corresponding diagonalization matrix. Inserting the above relation into I 4,5,6 of Eqs. (2.13)-(2.15) and using the definitions in Eqs. (2.7)-(2.12), it follows that they vanish identically (2.20) Using this, we can simplify the Lagrangians L 4,5,6 to Since we are interested in black holes, we need the stress-energy tensor corresponding to the above Lagrangian densities where we used the chain rule to express the stress-energy tensor in terms of for which we get where we have used δ r δg µν = r r µν , (2.32) and defined 1 µν ≡ n µ · n ν (2.33) 2 µν ≡ g µ 2 ν 1 (n µ · n ν 1 )(n µ 2 · n ν ) (2.34) r µν ≡ g µ 2 ν 1 (n µ · n ν 1 )(n µr · n ν ) r−1 p=2 g µ p+1 νp n µp · n νp , for r > 2. (2.35) Finally, we can write down the stress-energy tensor T (2n) µν for each component Lagrangian, L 2n , as We will now specialize to the case with spherical symmetry, which is relevant for the 1-Skyrmion sector. Thus we can assume the hedgehog Ansatz for the Skyrme field with a profile function f (r) satisfying the boundary condition f (r → ∞) → 0, which in terms of the four-vector n reads n = sin f (r) sin θ cos φ, sin f (r) sin θ sin φ, sin f (r) cos θ, cos f (r) . (2.43) In this paper, we will choose a metric tensor of the form which is compatible with a spherically symmetric Schwarzschild black hole.
We can now plug in the hedgehog Ansatz and metric to the component Lagrangians and corresponding stress-energy tensors. First we note that the invariants have an astonishingly simple form for the hedgehog where f r ≡ ∂ r f is the radial derivative of f . We will use this notation for derivatives throughout the paper. We will first evaluate the building blocks C s with the hedgehog and Schwarzschild metric (2.51) and the identities I s are readily checked to vanish It is now easy to write down the Lagrangian densities It is convenient to note that the invariant derived with respect to the inverse metric, has a very simple form once we plug in the hedgehog Ansatz and Schwarzschild metric. Due to spherical symmetry and diagonal metric, they remain diagonal and their components read Let us now evaluate all the once-derived building blocks with respect to the inverse metric, namely, C (s) µν .
and due to spherical symmetry, we have that T φφ = sin 2 (θ)T θθ .

Black hole Skyrme hair
Once we have chosen the Lagrangian, L, out of our set of Lagrangians and fixed the constants, we are ready to solve the Einstein equation fleshed out as 1 for the metric in Eq. (2.44). After taking suitable linear combinations, we can write them as which with two boundary conditions determine N, C and hence the metric. Finally, we also need the equation of motion coming from varying the matter Lagrangian The black hole horizon is defined as C(r h ) = 0, where r h is the horizon radius. The other boundary conditions we need to impose, is that the profile function goes to zero at spatial infinity, f (∞) = 0, and the correct value of the first derivative at the horizon radius, f r (r h ). The latter can be derived by taking the r → r h limit of the first Einstein equation and the equation of motion, yielding where lim r→r h C(r) = 0. Since the metric function C accompanies the radial derivative squared, the limit is relatively simple. Indeed, by looking at the time-time component of the energy-momentum tensors (2.69)-(2.84), only the kinetic term, the Skyrme term and the eighth-order term have a nonvanishing contribution to the right-hand side of Eq. (3.8).
We can rewrite Eq. (3.9) by using the fact that f r is necessarily accompanied by at least a factor of C. Therefore, the first term can only give a nonzero contribution when the radial derivative hits a factor of C. 2 Thus we can write where we have used Eq. (3.6).
In our calculations, we will use a variable simply related to C, namely i.e. m and it has the physical interpretation as a gravitational mass function; more precisely, lim r→∞ m = m(∞) = m ADM is the ADM mass. It is a good observable to distinguish stable and unstable branches of solutions by (without turning to a linear stability analysis).
As the inverse metric that accompanies the double derivative of the profile function, f rr , is g rr = C, the equation of motion is singular exactly at the horizon. For this reason, we will start all calculations a tiny step r δ from the horizon: i.e. at r = r h + r δ and extrapolate the values of the fields linearly using the derivatives of f and µ.
It is well known that the topological charge or topological degree is only a full integer when there is no BH horizon. That is, part of the topological charge is swallowed by the BH. In this paper we only consider the spherically symmetric hedgehog, for which the topological charge outside the horizon is less than unity. In particular, we have where f h ≡ f (r h ) is the value of the profile function at the horizon. B < 1 for f h < π, which is the case for all the BH hair solutions we find in this paper. Note that since the first derivative of B with respect to f h is zero at f h = π, the charge is close to one for a range of f h π.

The 2 + + 2n model
This model is based on the standard black hole Skyrme hair with an added higher-order term, which is higher order than the Skyrme term. The purpose of this model is to construct a stable modified black hole hair, and then slowly turn off the Skyrme term to see if the hair persists.

General setup
The model is defined as where n = 3, 4, 5, 6 specifies the higher-order term added to the normal black hole Skyrme hair, R is the Ricci scalar and G is Newton's constant.
We will now switch to dimensionless units where a sets the (inverse) length scale and has mass dimension 1. Hence, we can write the Lagrangian as where we have defined i.e. the effective (dimensionless) gravitational coupling, while M 0 will turn out to set the mass scale for the soliton hair, as we will show shortly. We will use the freedom of choosing the scales to remove c 2 and c 2n from the dynamical equations: (4.6) The Hence, we can now write the rescaled Lagrangian as where we have defined which can be verified to be a dimensionless parameter of the model. The dynamical equations now only depend on the effective gravitational coupling α, and the coupling β.
Finally, we will show that the soliton hair in the rescaled units has a mass given by M 0 times a dimensionless integral (number): and as a check, we can see that the Einstein equations are now dimensionless 10) and the suitable linear combinations in Eqs. (3.5)-(3.6) read for n = 3, 4, 5, 6. It will be useful in the next subsections, to have dimensionless expressions for the boundary conditions. Let us start with Eq. (3.8): whereas the matter equation of motion is slightly more involved. We will use the expression (3.10) and write it in dimensionless units as For the numerical calculations, we will use the dimensionless mass function, µ, given via The boundary conditions at the black hole horizon can now be written as the following expansion where f h ≡ f (ρ h ) is the shooting parameter and provides nonperturbative information about the soliton hair. µ ρ (ρ h ) is extracted from Eq. (4.13) by using that C ρ = 2µ i.e. the derivative of the profile function of the Skyrme hair at the horizon, is isolated in Eq. (4.14) and C ρ is eliminated by insertion of Eq. (4.13).
The last observable that we will use here, is the Hawking temperature is almost a local quantity; it can be found by just knowing the value of the profile at the horizon, f h . This is not so with N (ρ h ); this is a nonlocal quantity and it must be obtained by integrating N from infinity down to the horizon radius, ρ h , using Eq. (4.12), where we have used the boundary condition N (∞) = 1. The Hawking temperature defined above is dimensionless; to convert to the dimensionful temperature one should multiply by a yielding aT H . In the following subsections, we will study each case of n = 3, 4, 5, 6 in turn.

The 2 + 4 + 6 model
This model was already studied in Ref. [15], and the result was that the sixth-order BPS-Skyrme term cannot stabilize solitonic black hole hair. For completeness, we will review the results here and provide a few new figures. This will help facilitate a comparison between the characteristics of this model and the other ones. Let us first complete the Einstein equations (4.11)-(4.12), while the equation of motion for the profile function is Finally, the boundary conditions (4.13)-(4.14) can be written as We are now ready to present the numerical results for this model. Fig. 1 shows the following quantities: the profile function at the horizon (the shooting parameter) f h , the derivative of the profile function at the horizon f ρ (ρ h ), the ADM mass µ(∞) and finally, the Hawking temperature T H . They are all plotted as functions of the horizon radius, ρ h . These quantities are used in many figures in this paper and we will henceforth refer to them as the standard quantities. The stable branches are indicated with black solid lines and the red dashed lines display the unstable branches. The two branches meet in all the figures at the bifurcation point, which determines the largest possible black hole size that can support the black hole hair for the given parameters. The numbers in the figures  indicate the values of β for the different branches. From Fig. 1(c) it is easy to verify that the unstable branches all have larger ADM mass, for all values of the horizon radius, ρ h , than their corresponding stable branches do. A feature of the 2 + 4 + 6 model, which is not present in the standard Skyrme model (the 2 + 4 model), is that the unstable branches do not continue all the way back up to the flat-space limit (ρ h → 0) solution [14,15]. We can see two things that happen right about where the unstable branches cease to exist; Fig. 1(b) shows that (minus) the derivative of the profile function at the horizon, −f ρ (ρ h ), increases drastically and Fig. 1(d) shows that the temperature drops drastically, right before the branches cease to exist.
The most important point of Fig. 1 is that as β is turned off (sent to zero), even the stable branches tend to not exist. In order to show it more explicitly, we provide the same figures, but plotted as functions of β in Fig. 2. It is clear from Fig. 2   tends to zero, only parts of the branches with very small horizon radii, ρ h remain. In the limit of β → 0, no branch remains and the Skyrme hair ceases to exist. Fig. 2 shows what happens to the derivative of the profile function at the horizon as β is turned off. A bit counter-intuitive, perhaps, the derivative actually tends to zero. The ADM mass gets smaller as β → 0, which signals that the amount of black hole hair decreases, see Fig. 2(c). It is not completely clear-cut what happens to the Hawking temperature in the limit of β → 0, Fig. 2(d). It is clear that only the branches with very small horizon radii, ρ h , persist for small β. It seems that the bifurcation point moves slightly up as the horizon radius is decreased. In any case, for small ρ h , the unstable branches possess temperature curves that go drastically fast towards zero, until they cease to exist. The Figs. 1 and 2 provide solid evidence for the fact that the BH hair does not exist in the β → 0 limit.

The 2 + 4 + 8 model
This model contains the first term in increasing order of derivatives after the sextic term, also known as the BPS-Skyrme term -that is, an eighth-order derivative term. The minimal construction of Ref. [28] limits the terms such that there are no powers of a derivative in the i-th direction higher than four. This leaves us with two independent terms, as we will see and the coefficients in Ref. [28] are called c 8|4,4 and c 8|4,2,2 due to their composition in terms of eigenvalues, see Eq. (2.19). We can intuitively think of the two terms as the quadratic part of the Skyrme term squared and the cross terms, respectively, see Ref. [28] for details. The term with coefficient c 8|4,2,2 also has the interpretation as being the BPS-Skyrme term multiplied by the kinetic (Dirichlet) term. Furthermore, for c 8|4,2,2 = 2c 8|4,4 the two terms combine to be the Skyrme term squared.
The rescaled Lagrangian density (4.7) divides the term L 8 by a factor of c 8 . Since we have two independent coefficients, let us define After the rescaling, c 8 drops out and γ ∈ [0, 1] interpolates between the two terms. Here we will be interested in the following three possibilities. The two terms separately is one way to disentangle their effect, therefore we will consider γ = 0 and γ = 1. Furthermore, we will consider γ = 1 3 because it corresponds to the Skyrme term squared. The γ = 0 term is composed by one eigenvalue of the strain tensor in Eq. (2.18) to the fourth power, multiplied by the two other (nonzero) eigenvalues squared, i.e. λ 4 1 λ 2 2 λ 2 3 and then the product is symmetrized, yielding The γ = 1 term, on the other hand, is composed by only two eigenvalues of the strain tensor in Eq. (2.18), both to the fourth power, i.e. λ 4 1 λ 4 2 and then the product is again symmetrized over all 3 eigenvalues, yielding Finally, the γ = 1 3 term, which also corresponds to the Skyrme term squared, reads These three cases should be enough to determine which parts of the eighth-order term are essential for stabilizing the BH hair.
We are now ready to complete the Einstein equations (4.11)-(4.12),

and the equation of motion for the profile function reads
Finally, we need the boundary conditions (4.13)-(4.14), which can be written as We will start with the case of γ = 0. The numerical results for this case are shown in Fig. 3 and Fig. 4, where the standard quantities are shown as functions of the horizon radius, ρ h , and as functions of the Skyrme-term coefficient, β, respectively.
Qualitatively, everything is very similar to the 2 + 4 + 6 model. That is, the branches shrink as β tends to zero; in particular, the bifurcation point, which represents the largest possible black hole possessing BH hair, moves to smaller and smaller values of the horizon radius, ρ h , as β decreases. Fig. 3(c) clearly shows that the unstable branches have higher ADM masses than their respective stable ones. Like in the case of the 2 + 4 + 6 model, also in this model, the unstable branches end before reaching the flat space limit at f h = π for ρ h → 0. Just before the unstable branches end, (minus) the derivative of the profile function at the horizon, −f ρ (ρ h ) increases drastically, see Fig. 3(b), and the temperature drops drastically, see Fig. 3(d). Fig. 4 shows the same standard quantities as functions of β. Like in the 2 + 4 + 6 model, it is clear that all branches cease to exist in the limit of β → 0. In particular, we can see from Fig. 4(b) and (c) that (minus) the derivative of the profile function, −f ρ (ρ h ), goes towards zero and ADM mass decreases, for β → 0. A difference with respect to the 2 + 4 + 6 model, is that the Hawking temperature really drops drastically to zero for β → 0, see Fig. 4(d). More precisely, the stable branch with ρ h = 0.01 drops below the temperatures of the other branches, whereas in the 2 + 4 + 6 model, the stable branch only drops to about 2 and then turns into an unstable branch which decreases to values of about 10 −1 and then terminates. The reason why we are interested in knowing the  limiting behavior of these observables as functions of β in the limit of β → 0, is to better understand the cause of the collapse of the BH hair.
We will now turn to the case of the 2 + 4 + 8 model with γ = 1. The behavior of the BH hair in this model is vastly different from the 2 + 4 + 6 model and the 2 + 4 + 8 model with γ = 0. In particular, it is the first known Skyrme-type model other than the standard Skyrme model, which possesses stable BH hair (without the Skyrme term to stabilize it). This can easily be seen from Fig. 5(a) as the branches do not collapse to small horizon radii, ρ h , as β is turned off. In fact, a peculiarity of the model, is that the BH hair can support larger BHs without the Skyrme term than with it. This can be seen in Fig. 5(a) as the branches are extending to larger ρ h for smaller β. The largest branch corresponds to β = 0. We can confirm for all β, that the unstable branches have larger ADM masses than their corresponding stable ones, see Fig. 5   all figures at the bifurcation point, which also corresponds to the largest possible BH that can support the BH hair for the given parameters. We can also see that the derivatives of the profile function are steeper for the unstable branches than for their corresponding stable ones, see Fig. 5(b). It is observed from Fig. 5(d) that the temperature dependence on β is quite mild and that the unstable branches have lower temperatures than their corresponding stable ones. For completeness, we will present the standard quantities as functions of β as well, see Fig. 6, in order to show clearly the β → 0 limit. The detailed information in the figure is not so interesting, except for the fact that all stable branches have well-defined β → 0 limits. The unstable branches survive the β → 0 limit for ρ h 0.2. We can confirm that all unstable branches have higher ADM masses than their corresponding stable ones, see Fig. 6  lower temperatures than the corresponding stable ones, see Figs. 6(b) and (d), respectively. Finally, we will consider the last case of γ = 1 3 which corresponds to the Skyrme term squared. The standard quantities are shown in Fig. 7 as functions of the horizon radius, ρ h . Qualitatively, the behavior of this model is very similar to that of γ = 1. Both cases yield a stable BH hair in the limit of vanishing β; that is, the Skyrme term is not needed for stabilization of the solitonic hair. One difference with respect to the γ = 1 case, is that the unstable branches reach down to smaller BH sizes (smaller ρ h ). The branch with β = 0 -which means without the Skyrme term -is the largest branch, see Fig. 7(a) and also the branch with the smallest ADM mass, see Fig. 7  in the β → 0 limit, we will not show the corresponding figure with the standard quantities as functions of β. As can be seen from Fig. 7, the β dependence is mild.

The 2 + 4 + 10 model
The model we consider here is made of the standard Skyrme model with an added higherorder derivative term that consists of the Skyrme term multiplied by the BPS-Skyrme term. This model is simpler than the 2 + 4 + 8 model as it only has one overall coefficient c 10 which we scale away in the units, see Eq. (4.7).   (both stable and unstable), and decreases to β = 0 for the outermost branch. In (b) β = 1 starts as the lowest branch and crosses over to the higher branch as ρ h is increased, for the stable branch. For the unstable branch, β = 1 corresponds to the lowest branch. In (c) the topmost branches (both stable and unstable) correspond to β = 1 and the branches move downward as β → 0. The stable and unstable β = 0 branches are shown in orange and dark green colors, respectively.
Completing the Einstein equations (4.11)-(4.12), we have while the equation of motion for the profile function reads The boundary conditions (4.13)-(4.14) can be written as We can see from Figs. 8 and 9 that the 2 + 4 + 10 model is qualitatively similar to the 2 + 4 + 6 model and the 2 + 4 + 8 model with γ = 0. Indeed, it does not support BH hair when the Skyrme term is turned off, i.e. in the β → 0 limit. As usual for a model that cannot support BH hair (without the Skyrme term), we see in Fig. 9(a) that only the smallest BH sizes remain for small β and no branches are left in the β → 0 limit. From Figs. 8(b) and 9(b), we can see that the derivative of the profile function at the horizon, f ρ (ρ h ), tends to zero as β does. We confirm from Figs. 8(c) and 9(c) that the unstable branches indeed have larger ADM mass than their corresponding stable ones. The Hawking temperature tends asymptotically to zero where the unstable branches end, see Figs. 8(d) and 9(d). When the temperature tends to zero for the unstable branches, we can see from Figs. 8(b) and 9(b) that the derivative of the profile function at the horizon, f ρ (ρ h ), tends to increase sharply until the solution ceases to exist.
Finally, the Figs. 8 and 9 provide solid evidence for the fact that in the 2 + 4 + 10 model, the BH hair does not exist in the β → 0 limit.

The 2 + 4 + 12 model
The last model that we will consider in this section, is the 2 + 4 + 12 model, where the higher-order derivative term that is added to the standard Skyrme model is made of the baryon charge density to the fourth power -or equivalently the BPS-Skyrme term squared. This model also only has one overall coefficient of the highest-order term, c 12 , which we scaled away in Eq. (4.7).  Completing the Einstein equations (4.11)-(4.12), we get   whereas the equation of motion for the profile function is Finally, we need the boundary conditions (4.13)-(4.14), which can be written as  We are now ready to present the results of the numerical calculations for the 2 + 4 + 12 model; they are shown in Figs. 10 and 11. This model is unique among the models we have considered in this paper. Indeed it has existing branches of solutions in the β → 0 limit, but the unforeseen twist is that the lower branches start to go below the upper branches  in the ADM mass when β becomes smaller than about 0.3, see Fig. 10(c). This is evidence for a kind of phase transition, where the lower branch becomes the stable one. We note that the general tendency is that as β becomes smaller, the ADM mass decreases. Before the phase transition takes place (β > β crit ), the lower branches have higher ADM mass for the same value of horizon radius as their corresponding upper branches. There is a range in β where the ADM mass of the lower branch intersects that of the upper branch. Then finally, as mentioned just above, in the β → 0 limit, the lower branch has a lower ADM mass than the corresponding upper branch, for all values of ρ h that the branch exists.
As common for many of the models with higher higher-order derivative terms than fourth order, the unstable branches terminate on the way back from the bifurcation point at a finite horizon radius (ρ h > 0), see e.g. Fig. 10(a). This is related to the line in the (f h , ρ h ) diagram where Ξ vanishes, see e.g. Eq. (4.42). At that point, the first derivative of the profile function at the horizon, f ρ (ρ h ), would tend to minus infinity, which in turn would prevent a smooth soliton solution [15]. The implication here, of the unstable branch ending at a finite horizon radius, is that the upper branch (in f h ) will be the stable solution for small enough ρ h and then become unstable at a horizon radius where the lower branch starts, see Fig. 10(a). This means that if a BH possesses Skyrme hair in this model and has a size ρ h < ρ crit h , the upper branch will dictate the value of the profile function. If we now throw in some stuff, say a piano, the BH grows and at the point where ρ h > ρ crit h , the Skyrme hair will become unstable and decay to the stable solution, i.e. the lower branch. In order to check our claims, we will perform a stability analysis of this model in the next section.
In Fig. 10(b) we can see another peculiar feature happening after our claimed phase transition. Namely, (minus) the first derivative of the profile function at the horizon −f ρ (ρ h ), for the lower branch with β = 0 tends towards zero, whereas before the transition β > β crit , they tend to become very large.
The Hawking temperature for β > β crit has the usual feature for a stable model, e.g. the 2 + 4 + 8 model with γ = 1 or γ = 1 3 ; that is, the temperature for the upper branches all rise as the horizon radius tends to zero and decreases for larger BHs until the bifurcation point of the particular branch is reached, see Fig. 10(d). The lower branches for β > β crit tend backwards in the (T H , ρ h ) diagram (see Fig. 10(d)) until they suddenly drop and terminate. What happens when the phase transition, β < β crit is reached, is that the lower branches no longer tend to vanishing temperature and in the β → 0 limit, the Hawking temperature is at the same level as the upper branches.
We will now plot all the standard quantities as functions of the Skyrme-term coefficient, β, instead of the horizon radius, ρ h , see Fig. 11. The main point to notice is that for ρ h 0.1, the upper branches do not have a bifurcation point at any finite β and thus have a well-defined β → 0 limit. The last branch with a bifurcation point in Fig. 11(a) is the ρ h = 0.1 branch and said point is around β = 0.25. This indicates that β crit is probably between 0.2 and 0. 25  All the lower branches depicted in Fig. 11 are connected to the upper branches with β > β crit and therefore do not exhibit the peculiar behavior seen in Fig. 10. The derivative of the profile function at the horizon goes up slightly, but remains of order one in the β → 0 limit, see Fig. 11(b) and the Hawking temperature remains basically constant in the limit, see Fig. 11(d). This applies to the upper branches. Since the lower branches possess somewhat more complicated behavior around and after the phase transition, they are not depicted in Fig. 11.
We can see from Fig. 10(a) that the bifurcation point moves to smaller BH radii as β decreases from 1 to about 0.04 and then it turns around and continues to increase as β goes to zero. The bifurcation point, which corresponds to the largest BH size possessing hair, is plotted in Fig. 12.
We will now make a more detailed plot of β branches as functions of ρ h and in order to avoid clutter, we will make the first series of solutions with β ∈ [0.04, 0.44], see Fig. 13, and the second series of solutions with β ∈ [0, 0.04], see Fig. 14. In Fig. 13 the bifurcation point is decreasing in ρ h , as β decreases, and in Fig. 14 the bifurcation point is increasing in ρ h . The physics in Fig. 13 is very versatile; three different situations arise. The classification is based on the ADM masses, see Fig. 13(c). For β ≥ 0.44 the ADM mass of the lower branch is always above that of the upper branch. This makes lower branches of solutions unstable and everything is as usual. In a range of β ∈ (0.44, 0.12) (approximately), the lower branch has a lower ADM mass than that of the upper branch for a finite range in BH sizes ρ h ∈ [ρ h , ρ birfurcation ]. At ρ h = ρ h the ADM masses of the lower and upper branches cross and the lower branch becomes the unstable branch for ρ h < ρ h . Finally, for β ∈ [0, 0.12] the lower branch possesses a lower ADM mass than that of the upper branch for its entire domain of existence. The last branch, i.e. the smallest value of β for which the ADM mass of the lower branch crosses that of the upper branch is β = 0.14. For β = 0.12, the branch terminates just about where the ADM mass of the lower branch would cross over that of the upper branch. For β < 0.12, the lower branches behave drastically different than those with β ≥ 0.12. In the case of the ADM mass, instead of bending down, going up and then continuing back with a constant ADM mass as ρ h decreases (for β ≥ 0.12), it just goes downwards in an arc-like shape until the branch terminates, see Fig. 13(c).
This same transition point -β ≥ 0.12 -also has an impact on f h , see Fig. 13(a). Indeed, the lower branches keep a relatively high value (∼ 2) of the profile function at the horizon. For β < 0.12, approximately, the lower branches experience a paradigm shift and they continue downwards in an arc-like shape until they terminate. The derivative of the profile function at the horizon, −f ρ (ρ h ), also experiences the same paradigm shift around β ∼ 0.12, see Fig. 13(b). For β ≥ 0.12, the derivatives go upwards until the lower branches terminate whereas for β < 0.12, the derivatives go below that of the upper branches and only shortly before they terminate, they start to go upwards. Finally, this "paradigm shift" can also be seen in the Hawking temperature, see Fig. 13(d). For β ≥ 0.12 the end point of the lower branches moves downwards in ρ h in a monotonic fashion as β is decreased, whereas for β < 0.12 the behavior is different. Instead, the temperature is raised a bit and the branches move in a higher level, see Fig. 13(d). Fig. 13 is the most interesting one as it contains the "paradigm shifting" behavior and possible a phase transition. Fig. 14, however, is the remaining details for β ∈ [0, 0.04] which -as we have mentioned above -is separated from Fig. 13 in order to avoid clutter (overlapping curves) due to the fact that the bifurcation point turns around and starts moving up again, see Fig. 12.
Indeed the behavior of the curves in Fig. 14 is practically monotonic and the most interesting fact is that the β → 0 limit exists. The catch in this model is that the lower branch possesses a lower ADM mass than its corresponding upper branch. In order to verify the stability and be able to claim which of them is the stable branch of solutions, we will perform a linear stability analysis in the next section.

Linear stability of the 2 + 4 + 12 model
In all models, except the 2 + 4 + 12 model, the unstable branches had lower values of the profile function at the horizon, f h , and a higher ADM mass. Hence, we have not made a detailed stability analysis for those cases. The 2 + 4 + 12 model is different. In this model, the there is a critical value of β for which the upper (in f h ) branches have a higher ADM mass and hence are unstable (or at best metastable). For β < β crit this is the case and the lower (in f h ) branches start to have a lower ADM mass for a finite range in ρ h than their corresponding upper ones. In the limit of β → 0, the lower branch retains a lower ADM mass than the upper branch throughout its entire domain of existence.
Our claim is that there is a phase transition and the upper branch becomes unstable (or metastable) for a range in BH sizes (horizon radii). In order to back up this claim, we will here make a linear stability analysis of the model in question.
Starting with the Lagrangian of the 2 + 4 + 12 model and turning on time dependence of the profile function f (ρ) → f (ρ, τ ), we can write the matter part of the action as where we have defined the functions and the dimensionless time coordinate τ ≡ at, as well as δ ≡ log N, (5.5) which will be convenient shortly.
Turning on time dependence of the metric functions δ(ρ) → δ(ρ, τ ) and C(ρ) → C(ρ, τ ) does not alter the metric part of the usual two combinations of the Einstein equations but it gives and additional equation from the τ ρ component of the Einstein equations Finally, we also need the full time-dependent equation of motion for the profile function where u f ≡ ∂u ∂f , u ρ ≡ ∂u ∂ρ is the partial derivative of u (the total derivative of u is du dρ = u ρ + u f f ρ ), and similarly for the other functions.
Armed with the full time-dependent system of equations, it is now easy to write down the linearized perturbations. Hence, let us define C(ρ, τ ) ≡C(ρ) + C (ρ, τ ), (5.12) wheref ,N andC are the solutions of the soliton background while f , N and C are fluctuation fields about the background soliton. The linearization greatly simplifies especially the high powers of the time derivatives since and hence any power larger than one will eliminate the term. Following Refs. [4,5,8], we start by determining the perturbation of the radial metric function, C , by integrating Eq. (5.8), whereū means that u(f ) is evaluated with the background field: u(f ), and q(ρ) is an integration constant. We will now prove that the integration constant vanishes. In order to do so, we will first combine the two Einstein equations (5.6) and (5.7) to get Next, we will linearize Eq. (5.7): 17) and insert this as well as the equation of motion (5.9) evaluated on the background soliton (i.e. ∂ τ = 0) into Eq. (5.16), yielding which when integrated yields −ρC = 4α ū + 2wCf 2 ρ Cf ρ f + e −δq (τ ). (5.19) Comparing the above equation with Eq. (5.14), we can conclude that the integration constant cannot be time dependent. Using now that the appropriate boundary conditions for the fluctuations are that all fluctuations vanish at spatial infinity, we can finally conclude that q =q = 0. Finally, we need the linearized equation of motion for the perturbation of the profile function, f . After some massage, we can write the full perturbation in terms of that of the profile function as where the potential is where we have used the equation of motion for the background field,f , to eliminatev f . We will now set for which Eq. (5.20) is a regular Sturm-Liouville problem with a nontrivial (non-constant) weight or density function The right-hand side of the above equation is the weight function of the Sturm-Liouville problem and ω 2 ∈ R is the eigenvalue and by Sturm-Liouville theory it has to be real. Physically, if ω 2 < 0 then ω is imaginary and the perturbation mode signals an instability at the linear level. Full nonlinear stability requires the absence of linear instabilities, although that is generally not sufficient for claiming nonlinear stability. In this paper, we will content ourselves with the above linear stability analysis.
In the case of the standard Skyrme model, which corresponds to the case of w = 0, the stability is considerably simpler and the Sturm-Liouville problem can be transformed into a free eigenvalue problem with a complicated potential by setting ξ = ζ/ √ū and using tortoise coordinates defined by dx = dρ/(e δ C) [8]. The latter transformation is able to simplify the problem considerably because the weight function, r(x), and the kernel function, p(x), are equal to each other. In our case with w = 0, the kernel and the weight functions are different and we have not been able to find a suitable transformation to simplify the problem. We will therefore solve the Sturm-Liouville problem (5.23) directly using a numerical finite difference method.
The reason why the kernel and weight functions are different in our higher-order model as compared to the standard Skyrme model is due to the linearization of a higher-order derivative operator on a static background. That is, the linearization of f 4 ρ yields 4f 3 ρ f ρ , whereas the linearization of f 4 τ vanishes because of the static background. The higher-order term still gives a contribution to the time derivative of the fluctuation, but it comes from a cross term, f 2 ρ f 2 τ and hence does not give the same factor as the radial derivative of the fluctuation does.
We are now ready to present the numerical results for the lowest eigenvalue, ω 2 1 , of the perturbation of the profile function, f , shown in Fig. 15. If we start with the upper branches, they start off at a small horizon radius, ρ h , with a positive lowest eigenvalue (ω 2 1 > 0), which decreases as the horizon radius is increased. We would expect the eigenvalue to go towards zero and at the bifurcation point meet a negative eigenvalue coming from below zero. This does not happen in this model. If we start with the small values of β ∈ [0, 0.4], the behavior is as follows. Instead the lowest eigenvalue of the upper branch goes upwards and meets that of the lower branch at the bifurcation point. The lower branch hence does not have a negative mode and the lowest eigenvalue of the fluctuation for the lower branch is higher than that of the upper branch. Now, when β is increased to above about 0.5, the lower branch has a hybrid behavior. Instead of having a negative lowest eigenvalue of the fluctuation, it emanates from bifurcation point over that of the upper branch -stays there for a finite range in horizon radius -and then turns negative shortly before the branch terminates at a smaller horizon radius, see Fig. 15. There are only quite few data points with a negative eigenvalue of the fluctuation, but enough to conclude that the lower branch develops an unstable mode for large enough β 0.5.
The fact that the lower branch (in f h ) does not have an unstable mode over its entire domain of existence is indeed a surprise. Since the analysis carried out here is limited to a linear stability analysis, this does not rule out a nonlinear instability. Physically, the absence of an unstable mode in the fluctuation spectrum could be interpreted as the lower branch being a local minimum or vice versa. In order to determine which local minimum is the global minimum of solution space, we would still turn to the ADM mass. Our explanation at this point is simply that the highly nonlinear nature of our model has created a situation where the two branches both become local minima. Except for very particular points in parameter space, only one of them is a global minimum. With the above discussion in mind, we will now make a phase diagram of the 2 + 4 + 12 model, see Fig. 16. The expected features are the radius of the bifurcation point (black upper curve) and the point where the lower branch terminates (red lower curve). The crossover from where the lower branch has a lower ADM mass happens at β 0.5 and the curve decreases in β as the horizon radius, ρ h , gets smaller (orange dashed curve).
Since it is puzzling that the lower branches do not possess the expected unstable modes, we will consider the limit of turning off the twelfth-order term (but keeping the Skyrme term) in the 2 + 4 + 12 model in the next section.
6 Taking the Skyrme model limit of the 2 + 4 + 12 model In this section we will take the limit of the 2 + 4 + 12 model becoming the standard Skyrme model; i.e. turning off the twelfth-order derivative term. In order to do so, we need to rescale the Lagrangian such that the free parameter is the coefficient of the twelfth-order term and not the Skyrme term.
The model is defined as where the Cs are defined in Eqs. (2.7)-(2.9). Similarly to Sec. 4, we will now switch to dimensionless units, ρ ≡ ar, for which we get where the effective (dimensionless) gravitational coupling is still given by Eq. (4.4) and M 0 is the mass scale of the soliton. Instead of letting c 4 be the free parameter, we will fix the coefficient of L 4 and let c 12 be the free parameter. This is done by which both have mass dimension 1 as they must. The rescaled Lagrangian is now written as where the new parameter, η is defined as (6.6) In the limit of η → 0, this model becomes the Skyrme model. The soliton mass then reads The dimensionless Einstein equations read 8) and the equation of motion becomes  model (coupled to gravity). First we note that varying η in the range [0.1, 1] has little effect on both the upper and lower branches, see Fig. 17(a). The approach from the branch with η ∼ 1 to the Skyrme model is almost logarithmic in η. The η = 0 (i.e. the standard Skyrme model limit) branches of solutions are distinct from the nonzero η ones by the fact that the lower branch of the η = 0 returns smoothly to f ∼ π in the limit of vanishing horizon radius, ρ h → 0, see the dark green dashed curve in Fig. 17(a). Once a tiny nonzero η is turned on, the lower branch does not return all the way when ρ h → 0, but terminates at a finite horizon radius. The derivative of the profile function at the horizon also behaves differently for the lower η = 0 branch as compared to the nonzero η branches, see Fig. 17(b). The lower η = 0 branch remains almost constant in the ρ h → 0 limit, whereas the lower nonzero η branches raise up (in negative values) about two orders of magnitude before they terminate at a finite horizon radius. A final distinction between the lower η = 0 branch and the nonzero ones is seen in Fig. 17(d), where the Hawking temperature goes smoothly back from the bifurcation point in the ρ h → 0 limit. The lower nonzero η branches on the other hand, drop very drastically and suddenly in temperature where they terminate at a finite horizon radius.
The ADM masses for all the lower branches in this model are larger than all those of the upper branches, which is the expected (from the standard Skyrme model scenario) behavior, see Fig. 17(c). So far, except for the distinct behavior of the η = 0 branch, everything seems in line with the results for the 2 + 4 + 6 model studied in Refs. [14,15]. The expectation from the standard Skyrme model is that the lower branches have larger ADM masses than the upper branches, and they also possess a negative lowest eigenvalue in the linear fluctuation spectrum -signaling an instability at the linear level. Now we will consider the lowest eigenvalue for this model, see Fig. 18. Although all the ADM masses for the lower branches are higher, they do not all have a negative eigenvalue in their linear fluctuation spectrum over the entire range of horizon radii. If we start with the η = 0 branch, everything is as expected. The upper branch has a positive eigenvalue which drops suddenly near the bifurcation point (see Fig. 18(a)) and the lower branch has one single negative mode all the way in ρ h and it only increases near the bifurcation point to meet with that of the upper branch. If we turn on a very tiny η in the range [10 −5 , 10 −4 ], the trend continues as just described. Then for η around 10 −3 a transition occurs and the eigenvalue for the upper branch no longer drops to zero near the bifurcation point. Instead the lower branch now possesses only (linearly) stable modes for a finite range in ρ h from a finite value larger than where the branch terminates, up to the bifurcation point. For η = 10 −3 , the upper branch still has the largest eigenvalue compared to the lower branch, but that quickly changes and for η 0.01, there is a finite range in ρ h where the eigenvalue of the lower branch is larger than that of the upper branch. After this range, the eigenvalue of the lower branch drops suddenly to negative values where it remains until the branch terminates at a finite horizon radius, see Fig. 18(b).

The 2 + 2n model
In this section, we will study the dependence of the existing models with β = 0 on the effective gravitational coupling (4.4). These models are thus new Skyrme-type models that possess BH hair without having the Skyrme term component in the Lagrangian density. The 2 + 2n model is the β → 0 limit of the 2 + 4 + 2n model of Sec. 4. Only two models survive in the limit, namely the 2 + 4 + 8 model with γ = 0 and the 2 + 4 + 12 model. The γ = 0 condition in the 2 + 4 + 8 model corresponds to having a nonvanishing c 8|4,4 term in the Lagrangian density. Here we will consider the 2 + 8 model with γ = 1 and γ = 1 3 as well as the 2 + 12 model. The γ = 1 3 case of the 2 + 8 model corresponds to the higher-derivative term being the Skyrme term squared, whereas the γ = 1 case is the three curvatures of the Skyrme term individually squared without the corresponding cross terms. After rescaling, and in the case of the 2 + 8 model, fixing γ, the effective gravitational coupling is the only parameter of the model. As the β = 0 branches for fixed α = 0.01 have been studied already in Sec. 4, here we will consider only the standard quantities as functions of α, for a few different BH sizes (different horizon radii). The equations of motion and the boundary conditions for the 2 + 2n models with n = 4 and n = 6 are simply given in Sec. 4 with β = 0. Thus we will not repeat them here.
In Fig. 19, the standard quantities are shown for the γ = 1 case of the 2 + 8 model as functions of the gravitational coupling α. The different branches correspond to different horizon radii and as expected, the branch with the smallest horizon radius has the largest value of the profile function, f h . As the horizon radius is increased, the branches move downwards in f h , Fig. 19(a), as expected, since the branches will move downwards in f h to meet their respective bifurcation point. The upper branches are called stable here, as they everywhere have lower ADM mass than the corresponding lower branches, see Fig. 19(c). We have not carried out a linear perturbation analysis for this case, as we expect the lower branches to be unstable, if not linearly, then at best metastable. Although the ADM mass increases roughly linearly with the gravitational coupling, α (see Fig. 19(c)), the Hawking temperature remains almost constant as α is varied, for the stable branches. We made an extensive search for the unstable or lower branch for ρ h = 0.1, but were unable to find any solutions -both for large and small values of f h . Finally, let us mention that we have found that the unstable branches for ρ h = 0.4, 0.5 continue all the way to α → 0.
For completeness, we have made the same plots for the γ = 1 3 case in Fig. 20. Because these plots are quite similar to those of Fig.19, let us just mention the differences. Indeed the quantitative behavior is the same, but we have been able to find an unstable branch with horizon radius ρ h = 0.1. Both the unstable branches ρ h = 0.1 and ρ h = 0.2 end at a finite horizon radius, while those with ρ h = 0.3, 0.4 continue back in the limit of α → 0. Again the ADM masses suggest that the lower branches are unstable, see Fig. 20(c).  Finally, we will consider the last model, i.e. the 2 + 12 model, which only has the gravitational coupling, α, as a parameter (after rescaling of the length and energy units). The result is shown in Fig. 21. This model possesses a more complicated branch structure than the two flavors of the 2 + 8 model. The upper branches behave as expected; they start from above in f h with small horizon radius and move downwards as ρ h is increased. They have little dependence on α, expect near their bifurcation point. The lower branches, however are far more complicated. The lower branch with ρ h = 0.9 (brown dashed line) is the only one depicted which is single valued in α. The lower branches with ρ h = 0.7 (green dashed line) and ρ h = 0.5 (yellow dashed line) are still continuous, but not single valued in α. Surprisingly, the lower part of the ρ h = 0.5 lower branch exists up till quite large α 0.04, before it sharply turns back, see Fig. 21 Fig. 21(c). As the horizon radius is decreased, the lower branches become discontinuous. Indeed, the ρ h = 0.3 lower branch has an upper part connected to the upper branch and a disconnected lower part, see Fig. 21(a). Finally, the ρ h = 0.1 lower branch only exists very close to the bifurcation point, i.e. for quite large values of α. We did not find a disconnected lower part for this value of the horizon radius.
As mentioned above, the ADM masses interestingly show that the lower parts of the lower branches have lower ADM masses than their corresponding upper branches, which would make them the stable branches, Fig. 21(c). Note, however, that these lower parts do not exist all the way up to the bifurcation point, so for values of α close the bifurcation point, the upper branch would be the stable one. We have not carried out a linear perturbation analysis for this case, see comments in the next section.
It is interesting to see what happens to the derivative of the profile function at the  horizon for the lower parts of the lower branches, which according to the ADM masses are the stable ones, see Fig. 21(b). Let us consider the ρ h = 0.5 lower branch (yellow dashed line). The unstable part of the lower branch has a higher derivative at the horizon and about where this branch becomes the stable one, the derivative drops below that of the upper branch, see Fig. 21(b). The lower parts of the lower branches, which are the stable ones for those values of α, do in fact all have smaller derivatives of the profile function at the horizon.
Finally, we will show the Hawking temperature for this model. This calculation turned out to be the hardest one in the paper and we needed to use over 1 billion lattice points for our integrator to get convergence for the temperature. The interesting result from this calculation is that the lower parts (in f h ) of the lower branches, which are the stable ones for the corresponding values of α, have higher Hawking temperature compared with the upper parts of the same branches, see Fig. 21(d).

Discussion and conclusion
In this paper, we investigated a number of Skyrme-like models with terms containing six to twelve derivatives. The higher-derivative terms considered here are not the most generic ones, but the so-called minimal terms constructed in Ref. [28]. The main motivation was to get a better understanding of the criteria for when a Schwarzschild-type BH can support scalar hair. Indeed in Refs. [14,15], it was shown that although the Skyrme term can support BH hair, the sextic BPS-Skyrme term cannot. In this paper we have checked 3 further models, i.e. a model with a kinetic term and a 2n-th order term, n = 4, 5, 6. The 2 + 8 model is a one-parameter family of models and it turned out that it can support BH hair as long as the model does not purely consist of the BPS-Skyrme term times the standard kinetic term. One of the possibilities that are stable, is the Skyrme-term squared. The 2 + 10 model has turned out not to be able to sustain BH hair. Finally, the 2 + 12 model is basically the kinetic term and the BPS-Skyrme term squared and surprisingly it does possess a stable BH hair. The BH hair comes in two branches, one upper branch (in the profile function at the horizon, f h ) which is typically stable, and one lower branch (in f h ) which is typically unstable (see below, however).
A feature already seen in the generalized Skyrme model coupled to Einstein gravity, is that the unstable branches, for sizable coupling to the BPS-Skyrme term, end at a finite BH horizon radius simply because the temperature approaches zero. This can be viewed as the BH approaching an extremal BH state or more pragmatically as the derivative of the field profile at the BH horizon blowing up. This feature has turned out to be shared by the other models that do not possess BH hair in the limit where the Skyrme term is turned off; in particular, those are the 2 + 4 + 6, the 2 + 4 + 8 (γ = 0), and 2 + 4 + 10 models.
For the standard Skyrme model coupled to Einstein gravity, the upper and the lower branches of solutions correspond to the stable and unstable solutions. This was checked in Refs. [4,5] with a linear perturbation analysis which showed that the lower branch contains a single negative frequency of the perturbation modes of the linearized system. This system contains implicitly both the linear perturbations of the metric as well as of the Skyrme field radial profile. The metric perturbations are then eliminated, yielding a single field master equation. Evidence for the instability of the lower branches was already seen by calculating their respective ADM masses, for which the lower branches always possessed the higher ADM masses compared to the upper (stable) branches. However, in the 2 + 4 + 12 model the situation turned out to be somewhat more intricate. Indeed, when the Skyrme term is slowly turned off, the lower branch switched to become the one with a lower ADM mass. For an intermediate range of the Skyrme-term coefficient, the lower branch possesses a lower ADM mass compared to the upper branch, for a finite range from the bifurcation point down to a critical radius where the ADM masses cross over (details have been shown in Fig. 13). Finally, in the limit of a vanishing Skyrme-term coefficient, the lower branches of solutions remain the ones with the lowest ADM masses. Those lower branches, however, still terminate at a finite horizon radius. This behavior is mapped out in the phase diagram of Fig. 16. To this end, we have carried out a linear perturbation analysis of the BH hair system analogous to that of Refs. [4,5]. It turned out, however, that the problem at hand is more complicated due to the higher nonlinearity of the problem, which causes the kernel and the weight function of the resulting Sturm-Liouville problem to differ. Consequently, the master equation for the perturbations of the 2 + 4 + 12 model cannot be written as a Schrödinger equation and the full Sturm-Liouville problem needed to be solved. The result of this analysis was consistent with the naive conclusion from the ADM masses, namely, the lower branches become the stable ones in the limit of the Skyrme term being turned off in the 2 + 4 + 12 model. This result is surprising. As a double check we have tried turning off the twelfth-order term and thus returning to the standard Skyrme model, and indeed, the lowest eigenvalue of the perturbations returned to the standard behavior. That is, the upper branch has only positive eigenvalues and the lower branch has a single negative eigenvalue. Once both the Skyrme term and the twelfth-order term are turned on with sizable (order one) coefficients, the eigenvalue possesses a hybrid behavior; for small horizon radii the lowest eigenvalue is negative, but it turns positive at a finite radius and crosses over that of the upper branch until they meet at the bifurcation point.
Since the models under study in this paper are highly nonlinear problems, the linearized perturbation analysis may not suffice. Indeed, as a future investigation, full nonlinear stability should be considered seriously in order to understand the situation of the 2+4+12 model in the limit of a small or vanishing Skyrme-term coefficient.
A peculiar observation about the limit of the Skyrme-term coefficient being turned off in the 2 + 4 + 12 model, is that the explanation for the lower branches terminating at a finite horizon radius until now was that the Hawking temperature would go to zero, or equivalently the derivative of the profile function at the horizon would blow up. In this case, however, the lower branches still terminate at a finite radius, but the reason switches from being the derivative of the profile function blowing up (i.e. equal to the Hawking temperature dropping to zero) to the derivative going to zero. That also has as consequence that the solutions cease to exist, but the Hawking temperature remains finite.
To complicate the situation with the termination of the lower branches, in the 2 + 4 + 8 model (for any value of β, the Skyrme-term coefficient) the lower branches terminate at a finite horizon radius, but with a finite Hawking temperature and an apparently finite first derivative of the profile function at the horizon. Further studies are needed to conclude the fate of the small horizon-radius limit for these models.
Finally, in the 2+12 model which has BH hair stabilized by a twelfth-order term and no Skyrme term, a certain range of the gravitational coupling exists for which there are not 2 but 4 solutions with the same horizon radius, ρ h . According to their ADM masses, the lower part of the lower branch is the stable one. The upper branch may possess metastability and the upper part of the lower branch may be either unstable or metastable.
Throughout this paper, we have turned off a potential in order to keep the analysis as clean as possible. Although we have not checked, we think that most, if not all, results will be qualitatively similar if a potential would be turned on; in particular a standard pion mass term. The absence of BH hair in the 2 + 4 + 6 model without a mass term or other potential is thus confirmed; Ref. [14] carried out all their numerical calculations with the pion mass term present. Ref. [14] also gives a physical explanation for the lack of BH hair in the model with only the BPS-Skyrme term, which claims that the pressure becomes negative at the horizon due to the potential. This explanation cannot cover the case considered here, where we have excluded the potential altogether. It may be that the zero pressure from the BPS-Skyrme term is not sufficient for preventing a collapse of the hair. Further studies are needed for a conclusion on this issue.
Finally, a question tightly related to the above mentioned issue, is whether the stressenergy tensor of the model (excluding the kinetic term) can be rewritten in the form of a perfect fluid and whether this may be a criteria for whether a BH can possess stable hair or not. In Ref. [54], it was shown that the BPS-Skyrme model can be rewritten as a non-barotropic perfect fluid using the Eulerian formulation of a relativistic fluid. The fluid element velocity was identified with the baryon charge current. Not all models considered here can be written in terms of just the baryon charge current or baryon charge density. Although some can, others cannot and some are hybrids of the latter two options.
It is an open question whether such potentials can give rise to non-spherical BH Skyrmions, other than axially symmetric ones.
An obvious generalization of this study is to consider more general terms with the same number of derivatives; i.e. non-minimal Lagrangians. The first class of more general Lagrangians will contain only a second-order equation of motion (albeit nonlinear in derivative terms), but with more than four derivatives in the same direction, for the sixth, eighth, tenth and twelfth order terms. A further generalization is to include more than one derivative acting on the same field, yielding a higher-order equation of motion. We avoided such a complication here and in Ref. [28], because in general it could give rise to an Ostrogradsky instability. A first interesting question would be, if the sixth-order Lagrangian was relaxed to contain more general terms than the BPS-Skyrme term, would it be able to sustain stable BH hair. We will leave this question for future work.