# Black holes with electroweak hair – the detailed derivation

Romain Gervalle, Mikhail S Volkov

<figure><img src="image_1.png" /></figure>

## To cite this version:

Romain Gervalle, Mikhail S Volkov. Black holes with electroweak hair – the detailed derivation. Phys.Rev.D, 2025, 112 (6), pp.064018. 10.1103/gzhg-62fz. hal-05066063

HAL Id: hal-05066063

https://hal.science/hal-05066063v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Black holes with electroweak hair – the detailed derivation

Romain Gervalle¹, * and Mikhail S. Volkov¹, ², †

¹Institut Denis Poisson, UMR - CNRS 7013,
Université de Tours, Parc de Grandmont, 37200 Tours, France

²CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland

We present a very detailed derivation of solutions describing hairy black holes within the gravity-coupled Weinberg-Salam theory, which were previously reported in Phys.Rev.Lett. 133 (2024) 171402. These black holes support a strong magnetic field that polarizes the electroweak vacuum and creates a condensate of massive fields carrying superconducting currents along the black hole horizon. The currents, in turn, generate a “corona” of magnetic vortex segments attached to the horizon at both ends. The condensate and corona together constitute the black hole hair. The extremal solutions approach, in the far field, the magnetic Reissner-Nordström configuration, with a total mass that is lower than the total charge, $M < |Q|$, due to the negative Zeeman energy of the condensate. This makes the removal of the hair energetically unfavorable. The maximally hairy black holes exhibit masses comparable to terrestrial values, with approximately 11% of their total mass stored in the hair. Given that these solutions arise within a well-tested theoretical framework, they are likely to have physical relevance.

<footnote>* romain.gervalle@univ-tours.fr</footnote>

<footnote>† mvolkov@univ-tours.fr</footnote>
CONTENTS

I. INTRODUCTION 4

II. EINSTEIN-WEINBERG-SALAM THEORY 7

A. Field equations 9

B. ADM mass 9

C. Electromagnetic and Z fields, currents, and magnetic charge 10

1. Nambu definition 10

2. 't Hooft definition 11

3. Current densities and magnetic fields 11

4. Integrated currents 12

5. Magnetic charge 13

III. ABELIAN SOLUTIONS 13

A. Reissner-Nordström solution 14

Bertotti-Robinson solution 16

B. Reissner-Nordström-de Sitter solution 16

Extremal limit 17

IV. PERTURBATIONS AROUND THE RN BLACK HOLE 18

V. SPHERICALLY SYMMETRIC HAIRY SOLUTIONS FOR $n = \pm 2$ 21

A. The charge 23

B. The ADM mass 24

C. Extremal hairy black hole versus Cho-Maison monopole 26

D. Stability of the extremal hairy solution 27

E. The excited solutions 28

VI. PERTURBATIVE DESCRIPTION OF HAIRY BLACK HOLES WITH $|n| > 2$ 28

A. Electroweak equations in the unitary gauge 29

B. Condensate energy 30

C. Energy minimization 34

D. Structure of vortices 36

E. Axial symmetry 37

F. Special cases 37

1. $|n| = 2$ 37

2. $|n| = 4$ 37

3. $|n| = 3$ 38

G. Conclusion of the perturbative analysis 38

VII. AXIAL SYMMETRY 38

A. Static and axially symmetric gravitational fields 38

B. The integral formula for the ADM mass 40

C. Axially symmetric electroweak fields 41

D. Removing the string singularity 42

E. Fixing the gauge 43
F. Boundary conditions and numerical procedure 44

G. Virial relations 45

VIII. HAIRY SOLUTIONS WITH $n = \pm 4$ 46

A. Choice of the radial gauge 46

B. Structure of hairy solutions for $n = \pm 4$ 47

IX. PROPERTIES OF HAIRY BLACK HOLES 52

A. The horizon oblateness 54

B. The charge 54

C. The ADM mass 55

D. Quadrupole moments 56

E. Condensate, vector fields and currents 58

1. 't Hooft description 59

2. Nambu description 62

X. EXTREMAL SOLUTIONS 63

A. Phase I – spherical horizon, $|Q| < Q_*$ 64

1. Structure of vortices forming the corona 64

2. Vortex fluxes and the similarity with W-strings 66

3. The origin of vortices within the symmetric phase 67

4. Agreement between the different descriptions 69

5. Summary of the 't Hooft description 71

6. The magnetic ring, ADM mass, and the weakness of gravity 71

7. Solutions with $|Q| \ll Q_*$ versus electroweak monopoles 73

8. Increasing the charge 74

B. Phase II – oblate horizon, $|Q| > Q_*$ 75

C. Maximally hairy black holes 75

D. Summary 76

XI. CONCLUDING REMARKS 77

Appendix A. Energy-momentum tensor 79

Appendix B. Axially-symmetric form of the $n = \pm 2$ solutions 80

A. Non-extremal solutions 81

B. Extremal solutions 82

C. Fixing the radial gauge – summary 84

D. Choice of $r_H$ and diffeomorphisms 84

Appendix C. Anti-screening property of the condensate 85

1. Magnetic field along $z$-axis. 86

2. Magnetic field on a sphere. 86

Appendix D. Estimates of the magnetic field strength 88

Appendix E. Horizon limit of extremal black holes 89

A. The first group of perturbations 90

B. The second group of perturbations 90
C. The third group of perturbations 92

Appendix F. Electroweak vortices 94

1. Z-strings 94

2. W-strings 94

3. B-strings 95

References 95

# I. INTRODUCTION

In 1971, Ruffini and Wheeler formulated the no-hair conjecture according to which the only parameters of an isolated black hole should be its mass, angular momentum, and electric/magnetic charge [1]. This conjecture was supported by the uniqueness theorems showing that all *electrovacuum* black holes are described by the Kerr-Newman metrics [2–4] (see [5] for a review). Besides, the conjecture was confirmed by the no-hair theorems proven for a number of field theory models containing other (scalar, vector, etc.) physical fields [6–13].

The first broadly recognized counterexample to the conjecture was found in 1989 [14]: black holes supporting, outside the horizon, a non-trivial Yang-Mills field not characterized by any charge, hence “a hair” (see [15] for details). This finding was confirmed [16, 17] and then generalized in many ways, leading to a plethora of black holes supporting a Skyrme hair [18, 19], a Higgs hair [20, 21], various types of stringy hair [22–25], hairy black holes without spherical symmetry [26, 27], and many others (see [28] for a review).

More recent findings include black holes with scalar hair obtained by adjusting the scalar potential [29–31] or by considering non-minimally coupled scalars [32–35]. Besides, there exist interesting solutions describing black holes supporting “spinning clouds” made of ultra-light bosons [36–38], “spontaneously scalarized” black holes [39], black holes with a massive graviton hair [40, 41], etc.

Nowadays, solutions describing hairy black holes have become so numerous that even a brief review of this subject would be difficult. This naturally raises the question of the relevance of these solutions for describing something that may really exist in Nature. Unfortunately, it seems one cannot be too optimistic in this respect. All known solutions describing hairy black holes were obtained either within overly simplified models, or within exotic models relying on the existence of yet undiscovered particles and fields, or within alternative theories of gravity whose applicability remains to be confirmed.

It is not completely excluded that gravitons may have a non-zero mass or that there exists a dilaton field. It is equally possible that the Dark Matter contains ultralight bosons or that the Dark Energy can be described by a scalar-tensor theory. However, none of these possibilities have been confirmed yet. Therefore, the physical relevance of solutions describing hairy black holes with massive gravitons, or with a dilaton field, etc., is not immediately obvious. The same can be said about most of the other known solutions.

Some physical relevance may perhaps be attributed to the black holes supporting a pure Yang-Mills field, because the latter arises in QCD. In that case, however, quantum corrections are large and should totally destroy the classical contribution. Therefore, it does not make much sense to describe hairy black holes by classical solutions in this case, although such solutions may probably be used for some effective description [42].
What could be physically relevant solutions describing hairy black holes? To address this question, we explore in this work the coupling between two well-tested physical theories: the Standard Model of fundamental interactions and General Relativity. The former contains the QCD sector, where quantum effects are dominant, but in the electroweak sector, the quantum corrections are not too large. Therefore, it makes sense to study extended classical configurations in the Einstein-Weinberg-Salam (EWS) theory. This theory contains the electrovacuum sector with the Kerr-Newman black holes, but does it allow for some other black holes as well ?

In fact, allowing for a net magnetic charge, more general black hole solutions do exist in the EWS theory. However, as we shall see, magnetically charged hairy electroweak black holes are in general not spherically symmetric, hence they are difficult to construct. Spherical symmetry is possible in the unphysical limit of the theory when the weak mixing angle vanishes and the U(1) hypercharge field decouples. The corresponding hairy solutions were constructed long ago [21]. For the physical choice of the weak mixing angle, $\sin^2 \theta_W \approx 0.22$, spherical symmetry is possible only for one value of the magnetic charge $|P| = 1/e$, where $e$ is the electron charge. The corresponding solutions were obtained only recently in Ref. [43] (they will be reviewed below), but these black holes are too close to the Planck scale, which again raises the issue of applicability of the classical description. The latter applies if $|P| \gg 1$, but spherical symmetry is then lost.

An important contribution to the analysis of non-spherical black holes was made in Refs. [44, 45] by considering perturbations around a magnetically charged Reissner-Nordström (RN) black hole. For simplicity, the analysis was performed not within the full electroweak theory but within its truncated version, not containing the Z and Higgs bosons. However, the remaining field content of the model was sufficient to show that the perturbative black hole hair has a shape similar to a Platonic solid with only discrete symmetries [45]. We were able to carry out a similar analysis within the full electroweak theory and obtained similar results that will be presented below.

A qualitative description of magnetically charged black holes with electroweak hair was given by Maldacena [46]. These black holes become hairy and lose spherical symmetry because of electroweak condensation — a change in the vacuum structure induced by very strong magnetic fields [47–50]. The EWS theory admits the magnetically charged RN solution supporting a Coulombian magnetic field $B = P/r^2$, with the massive Z, W, and Higgs fields being in the vacuum. If the horizon radius $r_h$ is large, then $|B(r_h)| = |P|/r_h^2$ is small, and the configuration is stable. If $r_h$ is small and $|B(r_h)| > m_W^2/e$, where $m_W$ is the W-boson mass, then the configuration no longer minimizes the energy. The latter can be lowered by developing a condensate of the massive fields around the horizon.

The condensate carries superconducting currents, which form loops along the horizon and produce segments of magnetic vortices orthogonal to the horizon. Due to their mutual repulsion, these vortices are expected to be maximally separated from each other, forming a “corona” around the black hole horizon. This breaks the rotational symmetry.

If $|B(r_h)| > m_H^2/e$, where $m_H$ is the Higgs-boson mass, then, in addition to the massive condensate, the horizon should be surrounded by a region where the Higgs field approaches zero, and the full electroweak symmetry is restored. In the far field, all massive fields approach their vacuum values, and the magnetic field becomes spherically symmetric.

Following this scenario, various aspects of magnetic black holes have been discussed, but always at a phenomenological level [51–56]. Explicit solutions describing non-spherically symmetric hairy black holes have not been studied, as this would require 3D numerical sim-
ulations – a challenging task. However, our analysis of perturbations of the RN background suggests a simpler option.

Although the absolute minimum of the condensate energy corresponds to a configuration without continuous symmetries, the energy has another critical point describing an axially symmetric state. In this state, the corona degenerates into two oppositely directed multi-vortices, which is possible for any value of the magnetic charge. The individual vortices within the multi-vortex configuration repel each other, making the axially symmetric state unstable with respect to spreading into a fully fledged corona. However, this state is a stationary point of the energy and corresponds to an axially symmetric solution of the equations, whose construction requires only 2D simulations.

As we shall see below, the mass of extremal and near-extremal axially symmetric hairy black holes is smaller than their charge. Therefore, they cannot shed their hair and transition into RN black holes, as the mass of the latter exceeds the charge. However, they can further lower their energy when the two multi-vortices split into elementary vortices, forming a fully spread corona. In this state, they appear to approach an absolute energy minimum and become stable.

Thus, axially symmetric hairy black holes are close to the stable hairy configuration and can evolve toward it when perturbed.

Following the above scenario, we will explicitly construct axially symmetric hairy black holes, and the main steps of our analysis are as follows. After introducing the EWS theory in Section II, we describe its simplest solutions, including the RN and RN-de Sitter black holes, in Section III. In Section IV, we study linear perturbations around the RN background and identify the threshold values of the event horizon radius at which electroweak condensation begins. The condensate remains spherically symmetric only when the magnetic charge is $P = \pm 1/e$, and in Section V, we review the corresponding non-perturbative hairy black hole solutions, originally obtained in Ref. [43]. The structure of the non-spherically symmetric condensate for $|P| > 1/e$ is analyzed in Section VI by minimizing the energy of the non-linear perturbations around the RN background.

Restricting to the axially symmetric sector, we present in Section VII the ansatz for the gravitational and electroweak fields, along with the essential details: the integral formula for the mass, the removal of the Dirac string singularity, gauge fixing, boundary conditions, the numerical procedure, and the virial relations used to monitor numerical accuracy.

The profiles of the axially symmetric solutions for $P = \pm 2/e$ are illustrated in Section VIII, while the general properties of non-extremal axially symmetric hairy black holes are discussed in Section IX. Section X is pivotal to our discussion: there, we present the extremal hairy solutions, explain the vortex structure forming the corona, and describe the phase transition when the horizon symmetry changes from spherical to oblate. We also discuss the “weakness of gravity” leading to the relation $M < |Q|$, provide estimates for the mass, radius, and magnetic field of the “maximally hairy” black hole, and include many other important details.

Finally, Section XI contains a brief summary of results, a discussion of possible formation mechanisms for hairy electroweak black holes in the early universe, and potential methods for their detection. Many technical subtleties are addressed in the six Appendices.

We try to be as explicit as possible and present all the essential aspects of our analysis, aiming at a reader who may not necessarily be expert both in General Relativity and electroweak theory. We explain how to introduce dimensionless variables, define the black hole mass, charge, currents, and quadrupole moments, and handle the constraints. We also de-
scribe how to switch to the regular gauge, fix the residual gauge freedom, establish the virial relations, and compute the condensate energy, among other key steps. Although our text is long, we believe it is worth keeping everything in one place rather than splitting it into separate parts, as this would not necessarily improve readability.

A short and highly condensed version of this text can be found in [57]. The results presented below generalize the previous analysis of electroweak monopoles in flat space [58, 59], where the electroweak condensate is generated by the magnetic field of a pointlike hypercharge. Electroweak monopoles closely resemble extremal hairy black holes, but their mass is infinite because the energy of the pointlike magnetic charge diverges. Gravity provides a natural cutoff by creating an event horizon, which regularizes the charge's energy and renders the total mass finite.

## II. EINSTEIN-WEINBERG-SALAM THEORY

Denoting dimensionful quantities by boldfaced letters and dimensionless ones by ordinary letters, the action of the bosonic part of the gravity-coupled electroweak theory is

$$ \boldsymbol{S} = \frac{1}{c} \int \left( \frac{\boldsymbol{c}^4}{16\pi\boldsymbol{G}} \boldsymbol{R} + \boldsymbol{L}_{\text{WS}} \right) \sqrt{-g} d^4\boldsymbol{x}. \qquad (2.1) $$

Here, $c$ is the speed of light, $\boldsymbol{G}$ is Newton's constant, and $\boldsymbol{R}$ is the Ricci scalar. The electroweak Lagrangian is

$$ \boldsymbol{L}_{\text{WS}} = -\frac{1}{4} \boldsymbol{W}_{\mu\nu}^a \boldsymbol{W}^{a\mu\nu} - \frac{1}{4} \boldsymbol{B}_{\mu\nu} \boldsymbol{B}^{\mu\nu} - (\boldsymbol{D}_{\mu} \boldsymbol{\Phi})^{\dagger} \boldsymbol{D}^{\mu} \boldsymbol{\Phi} - \boldsymbol{\lambda} (\boldsymbol{\Phi}^{\dagger} \boldsymbol{\Phi} - \boldsymbol{\Phi}_0^2)^2. \qquad (2.2) $$

Here, the SU(2) and U(1) field strengths are, respectively,

$$ \boldsymbol{W}_{\mu\nu}^a = \partial_\mu \boldsymbol{W}_\nu^a - \partial_\nu \boldsymbol{W}_\mu^a + \boldsymbol{g}\epsilon_{abc}\boldsymbol{W}_\mu^b \boldsymbol{W}_\nu^c, \quad \boldsymbol{B}_{\mu\nu} = \partial_\mu \boldsymbol{B}_\nu - \partial_\nu \boldsymbol{B}_\mu, \qquad (2.3) $$

with $\partial_\mu = \partial/\partial\boldsymbol{x}^\mu$. The complex dublet scalar Higgs field and its covariant derivative are

$$ \boldsymbol{\Phi} = \begin{pmatrix} \phi_1 \\ \phi_2 \end{pmatrix}, \quad \boldsymbol{D}_\mu \boldsymbol{\Phi} = \left( \partial_\mu - \frac{i g'}{2} \boldsymbol{B}_\mu - \frac{i g}{2} \tau_a \boldsymbol{W}_\mu^a \right) \boldsymbol{\Phi}, \qquad (2.4) $$

where $\tau_a$ are the Pauli matrices. Here, $g$ and $g'$ are the gauge couplings, $\boldsymbol{\lambda}$ is the Higgs self-coupling, and $\boldsymbol{\Phi}_0$ is the Higgs field vacuum expectation value. Their values are known from experiments [60].

Denoting $g_0 = \sqrt{g^2 + g'^2}$, one can introduce dimensionless quantities:

$$ W_{\mu}^{a} = \frac{g}{\Phi_{0}} W_{\mu}^{a} \qquad B_{\mu} = \frac{g'}{\Phi_{0}} B_{\mu}, \qquad \Phi = \frac{\boldsymbol{\Phi}}{\Phi_{0}}, \qquad \beta = \frac{8\boldsymbol{\lambda}}{g_{0}^{2}}, \qquad g = \frac{g}{g_{0}}, \qquad g' = \frac{g'}{g_{0}}, \qquad x^{\mu} = \frac{\boldsymbol{x}^{\mu}}{l_{0}}. \qquad (2.5) $$

Here, the length scale and the corresponding mass scale are

$$ l_0 = \frac{1}{g_0 \Phi_0} = 1.53 \times 10^{-16} \text{ cm}, \qquad m_0 c^2 = \frac{\hbar c}{l_0} = \hbar c g_0 \Phi_0 = 128.9 \text{ GeV}. \qquad (2.6) $$

The rescaled gauge couplings are expressed in terms of the Weinberg weak mixing angle:

$$ g' \equiv \sin \theta_w \approx \sqrt{0.223}, \qquad g \equiv \cos \theta_w. \qquad (2.7) $$
The rescaled Higgs self-coupling is given by

$$ \beta \approx 1.88. \qquad (2.8) $$

The electron charge and fine structure constant are

$$ e = \hbar c g_0 \times e, \quad e \equiv gg' = \sin \theta_w \cos \theta_w = 0.416, \quad \alpha = \frac{e^2}{4\pi\hbar c} = \frac{e^2}{4\pi} \hbar c g_0^2 \approx \frac{1}{137}. \qquad (2.9) $$

Therefore,

$$ \frac{1}{cg_0^2} = \frac{e^2}{4\pi\alpha} \hbar = 1.89 \hbar, \quad \frac{1}{g_0} = \frac{\hbar c}{e} = \frac{e}{4\pi\alpha} e. \qquad (2.10) $$

The relation to the Planck mass $M_{Pl}$ and Planck length $L_{Pl}$ is

$$ m_0 = \sqrt{\frac{\alpha\kappa}{2e^2}} M_{Pl} = 1.05 \times 10^{-17} \times M_{Pl} \quad l_0 = \sqrt{\frac{2e^2}{\alpha\kappa}} L_{Pl} = 9.46 \times 10^{16} \times L_{Pl}. \qquad (2.11) $$

One defines the dimensionless gravitational coupling,

$$ \kappa = \frac{8\pi G\Phi_0^2}{c^4} = \frac{4e^2}{\alpha} \left(\frac{m_z}{M_{Pl}}\right)^2 = 5.30 \times 10^{-33}, \qquad (2.12) $$

where $m_z$ is the Z-boson mass.

Substituting everything into (2.1), the action becomes $S = \hbar \times S$ with

$$ S = \frac{e^2}{4\pi\alpha} \int \left( \frac{1}{2\kappa} R + \mathcal{L}_{\text{WS}} \right) \sqrt{-g} \, d^4 x, \qquad (2.13) $$

$$ \mathcal{L}_{\text{WS}} = -\frac{1}{4g^2} W_{\mu\nu}^a W^{a\mu\nu} - \frac{1}{4g'^2} B_{\mu\nu} B^{\mu\nu} - (D_{\mu}\Phi)^{\dagger} D^{\mu}\Phi - \frac{\beta}{8} (\Phi^{\dagger}\Phi - 1)^2, $$

where $B_{\mu\nu} = \partial_{\mu}B_{\nu} - \partial_{\nu}B_{\mu}$ and

$$ W_{\mu\nu}^a = \partial_{\mu}W_{\nu}^a - \partial_{\nu}W_{\mu}^a + \epsilon_{abc}W_{\mu}^bW_{\nu}^c, \quad D_{\mu}\Phi = \left(\partial_{\mu} - \frac{i}{2}B_{\mu} - \frac{i}{2}\tau_a W_{\mu}^a\right)\Phi. \qquad (2.14) $$

This is the theory of the interacting gravitational field $g_{\mu\nu}$, the complex-valued Higgs field $\Phi = (\phi_1, \phi_2)^T$, the U(1) hypercharge field $B = B_\mu dx^\mu$, and the SU(2) field $W = T_a W_\mu^a dx^\mu$, where $T_a = \tau_a/2$ are the SU(2) generators.

The action (2.13) is invariant under spacetime diffeomorphisms and also under SU(2)×U(1) gauge transformations:

$$ \Phi \rightarrow \mathcal{U}\Phi, \quad \mathcal{V} \rightarrow \mathcal{U}\mathcal{V}\mathcal{U}^{-1} + i\mathcal{U}\partial_{\mu}\mathcal{U}^{-1}dx^{\mu}, \qquad (2.15) $$

where

$$ \mathcal{V} = \frac{1}{2} (B_{\mu} + \tau_{a} W_{\mu}^{a}) dx^{\mu}, \quad \mathcal{U} = \exp \left( \frac{i}{2} \nu + \frac{i}{2} \theta^{a} \tau_{a} \right), \qquad (2.16) $$

with $\nu$ and $\theta^a$ being functions of $x^\mu$.
A. Field equations

Varying the action (2.13) with respect to $B_\mu$, $W_\mu^a$, and $\Phi$ yields the EW equations:

$$ \nabla^\mu B_{\mu\nu} = g'^2 \frac{i}{2} (\Phi^\dagger D_\nu \Phi - (D_\nu \Phi)^\dagger \Phi), $$

$$ \mathcal{D}^\mu W_{\mu\nu}^a = g^2 \frac{i}{2} (\Phi^\dagger \tau^a D_\nu \Phi - (D_\nu \Phi)^\dagger \tau^a \Phi), $$

$$ D_{\mu} D^{\mu} \Phi - \frac{\beta}{4} (\Phi^{\dagger} \Phi - 1) \Phi = 0, \qquad (2.17) $$

where

$$ \mathcal{D}_{\mu} W_{\alpha\beta}^{a} = \nabla_{\mu} W_{\alpha\beta}^{a} + \epsilon_{abc} W_{\mu}^{b} W_{\alpha\beta}^{c}, \qquad (2.18) $$

with $\nabla_\mu$ being the geometrical covariant derivative with respect to the spacetime metric $g_{\mu\nu}$. Varying the action with respect to the latter yields the Einstein equations:

$$ G_{\mu\nu} = \kappa T_{\mu\nu}, \qquad (2.19) $$

where $G_{\mu\nu} = R_{\mu\nu} - (R/2)g_{\mu\nu}$ is the Einstein tensor, and the energy-momentum tensor is given by

$$ T_{\mu\nu} = \frac{1}{g^2} W_{\mu\sigma}^a W_{\nu}^{a\sigma} + \frac{1}{g'^2} B_{\mu\sigma} B_{\nu}^{\sigma} + (D_{\mu}\Phi)^{\dagger} D_{\nu}\Phi + (D_{\nu}\Phi)^{\dagger} D_{\mu}\Phi + g_{\mu\nu} \mathcal{L}_{\text{WS}}. \qquad (2.20) $$

The vacuum is defined as the configuration with $T_{\mu\nu} = 0$. Modulo gauge transformations, it can be chosen as

$$ W_{\mu}^{a} = B_{\mu} = 0, \qquad \Phi = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \qquad g_{\mu\nu} = \eta_{\mu\nu}, \qquad (2.21) $$

where $\eta_{\mu\nu} = \text{diag}[1, 1, 1, 1]$ is the Minkowski metric.

Allowing for small fluctuations around the vacuum and linearizing the field equations with respect to these fluctuations yields the perturbative mass spectrum of the theory. The spectrum includes the massless photon and graviton, as well as the massive Z, W and Higgs bosons with masses

$$ m_Z = \frac{1}{\sqrt{2}}, \quad m_W = g m_Z, \quad m_H = \sqrt{\beta} m_Z, \qquad (2.22) $$

in units of the mass scale $m_0$.

B. ADM mass

We shall solve the field equations to obtain asymptotically flat solutions with a finite ADM mass. The metric $g_{\mu\nu}$ then approaches the Minkowski metric $\eta_{\mu\nu}$ at spatial infinity. Introducing the Schwarzschild radial coordinate $r$, the mass is determined from the asymptotic form of the $g_{00}$ metric coefficient as $r \to \infty$,

$$ -g_{00} = 1 - \frac{2M}{r} + \dots = 1 - \frac{2\mathbf{GM}}{c^2 r} + \dots, \qquad (2.23) $$
where $\mathbf{r} = l_0 r$, and the dots denote subleading terms. Therefore,

$$ \frac{\mathbf{M}}{\mathbf{m}_0} = \frac{\mathbf{c}^2 l_0}{\mathbf{G} \mathbf{m}_0} M = \frac{8\pi}{\kappa} \times \frac{M}{\hbar c g_0^2} = \frac{e^2}{4\pi\alpha} \times \mathcal{M} = 1.89 \mathcal{M}, \qquad (2.24) $$

where we have used (2.10), and where the rescaled dimensionless mass is

$$ \mathcal{M} = \frac{8\pi}{\kappa} M. \qquad (2.25) $$

### C. Electromagnetic and Z fields, currents, and magnetic charge

There is no unique way to define electromagnetic and Z fields away from the Higgs vacuum, but two special definitions are often considered in the literature.

#### 1. Nambu definition

According to Nambu, the electromagnetic and Z fields can be defined as [61],

$$ \mathcal{F}_{\mu\nu} = \frac{g}{g'} B_{\mu\nu} - \frac{g'}{g} n^a W_{\mu\nu}^a, \quad \mathcal{Z}_{\mu\nu} = B_{\mu\nu} + n^a W_{\mu\nu}^a, \qquad (2.26) $$

where

$$ n^a = \xi^\dagger \tau^a \xi, \quad \xi = \frac{\Phi}{\sqrt{\Phi^\dagger \Phi}}, \quad \xi^\dagger \xi = 1, \quad n^a n^a = 1. \qquad (2.27) $$

The Higgs vacuum can be defined as a state where, modulo gauge transformation, one has

$$ W_{\mu}^{1} = W_{\mu}^{2} = 0, \quad \Phi = \begin{pmatrix} 0 \\ 1 \end{pmatrix}. \qquad (2.28) $$

This choice is compatible with the equations for $W_\mu^1$, $W_\mu^2$, and $\Phi$ contained in Eqs. (2.17). The field strengths then admit potential representations,

$$ \mathcal{F}_{\mu\nu} = \partial_{\mu} \mathcal{A}_{\nu} - \partial_{\nu} \mathcal{A}_{\mu}, \quad \mathcal{Z}_{\mu\nu} = \partial_{\mu} \mathcal{Z}_{\nu} - \partial_{\nu} \mathcal{Z}_{\mu}, \qquad (2.29) $$

where the potentials are given by

$$ \mathcal{A}_{\mu} = \frac{g}{g'} B_{\mu} + \frac{g'}{g} W_{\mu}^{3}, \quad \mathcal{Z}_{\mu} = B_{\mu} - W_{\mu}^{3}. \qquad (2.30) $$

The remaining non-trivial field equations in (2.17) then reduce to:

$$ \nabla^{\mu} \mathcal{F}_{\nu\mu} = 0, \quad \nabla^{\mu} \mathcal{Z}_{\nu\mu} + m_{\nu}^{2} \mathcal{Z}_{\nu} = 0. \qquad (2.31) $$

Therefore, $\mathcal{F}_{\mu\nu}$ satisfies the source-free Maxwell equations. The 2-forms

$$ \mathcal{F} = \frac{1}{2} \mathcal{F}_{\mu\nu} dx^{\mu} \wedge dx^{\nu}, \quad \mathcal{Z} = \frac{1}{2} \mathcal{Z}_{\mu\nu} dx^{\mu} \wedge dx^{\nu} \qquad (2.32) $$

are closed, satisfying

$$ d\mathcal{F} = 0, \quad d\mathcal{Z} = 0. \qquad (2.33) $$

Away from the Higgs vacuum, however, the form $\mathcal{F}$ is not necessarily closed, i.e., $d\mathcal{F} \neq 0$ (similarly for $\mathcal{Z}$), since there is no fundamental reason for the Maxwell equations to hold in this regime. This is particularly useful for our purposes, as it allows us to introduce a smooth magnetic charge density.
## 2. 't Hooft definition

The electromagnetic and Z fields can be defined as closed forms even away from the Higgs vacuum by using the prescription of 't Hooft [62]. The key observation is the identity

$$n^a W_{\mu\nu}^a + g^2 \psi_{\mu\nu} = \partial_{\mu} W_{\nu} - \partial_{\nu} W_{\mu}, \qquad (2.34)$$

where

$$g^2 \psi_{\mu\nu} = -\epsilon_{abc} n^a \mathcal{D}_\mu n^b \mathcal{D}_\nu n^c, \quad W_\mu = n^a W_\mu^a + i (\xi^\dagger \partial_\mu \xi - \partial_\mu \xi^\dagger \xi), \qquad (2.35)$$

with $\mathcal{D}_\mu n^a \equiv \partial_\mu n^a + \epsilon_{abc} W_\mu^b n^c$. This identity allows one to introduce the closed 2-forms

$$F_{\mu\nu} = \mathcal{F}_{\mu\nu} - e\psi_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}, \quad Z_{\mu\nu} = Z_{\mu\nu} + g^2\psi_{\mu\nu} = \partial_{\mu}Z_{\nu} - \partial_{\nu}Z_{\mu}, \qquad (2.36)$$

where the corresponding potentials are given by

$$A_{\mu} = \frac{g}{g'} B_{\mu} - \frac{g'}{g} W_{\mu}, \quad Z_{\mu} = B_{\mu} + W_{\mu}. \qquad (2.37)$$

The fields $\psi_{\mu\nu}$, $F_{\mu\nu}$, $Z_{\mu\nu}$, and $Z_\mu$ are invariant under gauge transformations (2.16), whereas the electromagnetic potential transforms as $A_\mu \to A_\mu + (1/e) \partial_\mu v$. In the unitary gauge, where $\Phi^T = (0, \phi)$, one finds

$$g^2 \psi_{\mu\nu} = W_{\mu}^{1} W_{\nu}^{2} - W_{\nu}^{1} W_{\mu}^{2}, \quad A_{\mu} = \frac{g}{g'} B_{\mu} + \frac{g'}{g} W_{\mu}^{3}, \quad Z_{\mu} = B_{\mu} - W_{\mu}^{3}. \qquad (2.38)$$

This equation is often used in the literature as the definition of the field potentials. However, the unitary gauge is not always suitable when the Higgs field approaches zero. For this reason, it is preferable to use the definition (2.37), which does not require imposing the unitary gauge.

Since $\psi_{\mu\nu}$ in (2.38) is constructed from the massive W-field, it vanishes exponentially fast in the far field zone. Consequently, both $\mathcal{F}_{\mu\nu}$ and $F_{\mu\nu}$ reduce to the same Maxwell tensor at large distances. This implies that the global quantities expressed as surface integrals at spatial infinity remain identical in both the Nambu and 't Hooft descriptions. However, at short distances, the two frameworks provide different interpretations, which, in a certain sense, are dual to each other.

## 3. Current densities and magnetic fields

The approach of 't Hooft is more commonly used in the literature since $A_\mu$ and $Z_\mu$ can serve as fundamental variables, allowing the entire set of electroweak equations to be formulated in terms of these fields (see Eqs.(6.8) below). The same is not possible with the Nambu fields as they do not admit potential, making 't Hooft's formulation more fundamental in this regard.

However, a key advantage of $\mathcal{F}_{\mu\nu}$ is that it allows for the definition of a smooth magnetic charge density, which can sometimes offer valuable insights. In 't Hooft's approach, the magnetic charge is always pointlike; nevertheless, we shall see that it can effectively "mimic" a smooth charge distribution.
Therefore, in what follows, we will consider both $\mathcal{F}_{\mu\nu}$ and $F_{\mu\nu}$, along with their associated currents,

$$4\pi \mathcal{J}^{\mu} = \frac{1}{\sqrt{-g}} \partial_{\nu} (\sqrt{-g} \mathcal{F}^{\mu\nu}), \quad 4\pi J^{\mu} = \frac{1}{\sqrt{-g}} \partial_{\nu} (\sqrt{-g} F^{\mu\nu}). \qquad (2.39)$$

The rescaled $(-g/g')J^{\mu}$ is sometimes referred to as neutral current because it sources the Z field (see Eq. (6.8b)). Introducing the dual tensors,

$$\tilde{\mathcal{F}}^{\mu\nu} = \frac{1}{2\sqrt{-g}} \epsilon^{\mu\nu\alpha\beta} \mathcal{F}_{\alpha\beta}, \quad \tilde{F}^{\mu\nu} = \frac{1}{2\sqrt{-g}} \epsilon^{\mu\nu\alpha\beta} F_{\alpha\beta}, \qquad (2.40)$$

the Nambu definition allows one to define the magnetic current as:

$$4\pi \tilde{\mathcal{J}}^{\mu} = \frac{1}{\sqrt{-g}} \partial_{\nu} (\sqrt{-g} \tilde{\mathcal{F}}^{\mu\nu}), \qquad (2.41)$$

which can be non-trivial. In contrast, within the 't Hooft definition, this current is locally zero.

We will focus on static and purely magnetic systems. Denoting $\mathcal{N} = \sqrt{|g_{00}|}$ and $h = \sqrt{\det(g_{ik})}$, the magnetic field is defined as:

$$\mathcal{B}^i = \tilde{\mathcal{F}}^{0i} = \frac{1}{\mathcal{N}} \frac{1}{2h} \epsilon^{ijk} \mathcal{F}_{jk} \Rightarrow \mathcal{B}_i = \frac{1}{\mathcal{N}} \frac{h}{2} \epsilon_{ijk} \mathcal{F}^{jk} \Rightarrow \mathcal{F}^{ik} = \frac{\mathcal{N}}{h} \epsilon^{ikm} \mathcal{B}_m. \qquad (2.42)$$

The magnetic field in the 't Hooft definition will be also denoted by $\mathcal{B}^i = \tilde{\mathcal{F}}^{0i}$, with no ambiguity, as it will always be clear from the context which definition is being used. Similarly, we define the massive magnetic fields,

$$\mathcal{B}_{Z}^{i} = \tilde{Z}^{0i}, \quad \mathcal{B}_{\tilde{Z}}^{i} = \tilde{\mathcal{Z}}^{0i}, \qquad (2.43)$$

but we introduce a different symbol for the hypermagnetic field:

$$\mathcal{B}^i = \tilde{\mathcal{B}}^{0i}. \qquad (2.44)$$

### 4. Integrated currents

Injecting (2.42) into (2.39) yields the following relations:

$$4\pi\sqrt{-g} \mathcal{J}^i = \epsilon^{ikm} \partial_k (\mathcal{N}^2 \mathcal{B}_m) \Rightarrow 4\pi\sqrt{-g} \epsilon_{kmi} \mathcal{J}^i = \partial_k (\mathcal{N}^2 \mathcal{B}_m) - \partial_m (\mathcal{N}^2 \mathcal{B}_k). \qquad (2.45)$$

Hence, introducing 1-forms $\mathcal{B} = \mathcal{B}_k dx^k$ and $\mathcal{J} = \mathcal{J}_k dx^k$, we obtain:

$$4\pi\mathcal{N} * \mathcal{J} = d(\mathcal{N}^2\mathcal{B}) = d(\mathcal{N} * \mathcal{F}), \qquad (2.46)$$

where the star * denotes the Hodge duality operator on the spacelike 3-surface. Both sides of these relations are 2-forms. Integrating over a 2-surface $\Sigma$ and applying Green's theorem yields the total current through $\Sigma$,

$$4\pi \mathcal{I} \equiv 4\pi \int_{\Sigma} \mathcal{N} * \mathcal{J} = \oint_{\partial\Sigma} \mathcal{N}^2 \mathcal{B}_m dx^m = \frac{1}{2} \oint_{\partial\Sigma} \sqrt{-g} \mathcal{F}^{ik} \epsilon_{ikm} dx^m, \qquad (2.47)$$

which is the Ampère's law relating the electric current to the the circulation of the magnetic field. Replacing everywhere $\mathcal{F}_{ik} \rightarrow F_{ik}$ yields similarly the 't Hooft current $\mathcal{I}$. The dimensionful currents expressed in amperes are:

$$\mathcal{I} = c\Phi_0 \mathcal{I} = \frac{e}{4\pi\alpha} \frac{e}{t_0} \mathcal{I} = \mathcal{I} \times 1.42 \times 10^8 \text{A}, \quad \mathcal{I} = c\Phi_0 \mathcal{I}, \qquad (2.48)$$

where $t_0 = l_0/c = 5.1 \times 10^{-27}$ sec.
## 5. Magnetic charge

The magnetic charge density in the Nambu definition is given by:

$$4\pi \tilde{\mathcal{J}}^0 = \frac{1}{\sqrt{-g}} \partial_{\nu} \left( \sqrt{-g} \tilde{\mathcal{F}}^{0\nu} \right) = \nabla_k \mathcal{B}^k = \frac{e}{2} \frac{1}{\sqrt{-g}} \epsilon^{ijk} \partial_i \psi_{jk}, \qquad (2.49)$$

whereas in the 't Hooft definition, it is locally zero. However, the total magnetic charge, defined via the flux through a closed two-surface at spatial infinity, is the same in both descriptions:

$$P = \frac{1}{4\pi} \oint_{\partial\mathbb{R}^3} \mathcal{F} = \frac{1}{4\pi} \oint_{\partial\mathbb{R}^3} F, \qquad (2.50)$$

because at infinity one has $\mathcal{F} = F$.

Since $F$ is closed, the second integral in (2.50) is determined solely by the topology of $F$ and remains unchanged as the integration surface shrinks. Consequently, within the 't Hooft definition, the magnetic charge is always pointlike, akin to the Dirac monopole. Its density can only have a distributional support, in which case one can introduce a potential such that, locally, $F = dA$.

At the same time, since we are working within the non-Abelian theory, it is meaningful to consider a scenario where the magnetic charge is smoothly distributed in space, as described by the density in Eq. (2.49). To achieve this, we use the Nambu field $\mathcal{F}$, which lacks potential but allows for a well-defined magnetic charge within a finite volume $V$:

$$P(V) = \frac{1}{4\pi} \oint_{\partial V} \mathcal{F} = \frac{1}{4\pi} \oint_{\partial V} (F + e\psi) = P + \frac{e}{4\pi} \oint_{\partial V} \psi. \qquad (2.51)$$

Since $\psi$ is exponentially suppressed at infinity, its flux vanishes in the asymptotic limit, ensuring that $P(V)$ approaches $P$ as $V \to \infty$.

## III. ABELIAN SOLUTIONS

We begin our analysis with the simplest case: spherically symmetric solutions. The assumption of spherical symmetry is quite restrictive, allowing either solutions with known metrics or, for a specific value of the magnetic charge $|P| = 1/e$, simple hairy solutions. However, it is natural to start with the simplest cases. Moreover, studying perturbations around the Reissner-Nordström background provides valuable insight into more general solutions. For these reasons, we shall focus on the spherically symmetric sector for now.

The static and spherically symmetric metric can be parameterized as:

$$ds^2 = -\sigma^2(r)N(r)dt^2 + \frac{dr^2}{N(r)} + r^2(d\vartheta^2 + \sin^2\vartheta d\varphi^2). \qquad (3.1)$$

The static, spherically symmetric, and purely magnetic electroweak fields can be expressed as (see [58] for details):

$$W = f(r) [T_2 d\vartheta - \nu T_1 \sin\vartheta d\varphi] + T_3 \nu \cos\vartheta d\varphi, \quad B = \nu \cos\vartheta d\varphi, \quad \Phi = \begin{pmatrix} 0 \\ \phi(r) \end{pmatrix}. \qquad (3.2)$$
The gauge fields exhibit a Dirac string singularity at $\vartheta = 0, \pi$ due to their non-vanishing azimuthal components. However, as will be explicitly shown below in the more general axially symmetric case, this singularity can be gauged away if the constant $\nu$ is given by:

$$ \nu = -\frac{n}{2}, \quad n \in \mathbb{Z}. \qquad (3.3) $$

Substituting (3.1), (3.2) into the Einstein-Weinberg-Salam equations (2.17), (2.19), the angular variables separate, yielding the following system of ordinary differential equations:

$$ (N\sigma f')' = \sigma \left( \frac{f^2 - 1}{r^2} + \frac{g^2}{2} \phi^2 \right) f, \qquad (3.4a) $$

$$ (r^2 N \sigma \phi')' = \sigma \left( \frac{\beta r^2}{4} (\phi^2 - 1) + \frac{1}{2} f^2 \right) \phi, \qquad (3.4b) $$

$$ m' = \frac{\kappa}{2} \left( N U_1 + U_0 + \frac{\nu^2}{2g'^2 r^2} \right) = \frac{\kappa}{2} r^2 T_{\hat{0}\hat{0}}, \qquad (3.4c) $$

$$ \frac{\sigma'}{\sigma} = \frac{\kappa}{r} U_1, \qquad (3.4d) $$

$$ (\nu^2 - 1) f' = (\nu^2 - 1) f \phi = 0, \qquad (3.4e) $$

where

$$ U_1 = \frac{1}{g^2} f'^2 + r^2 \phi'^2, \quad U_0 = \nu^2 \frac{(f^2 - 1)^2}{2g^2 r^2} + \frac{1}{2} f^2 \phi^2 + \frac{\beta r^2}{8} (\phi^2 - 1)^2. \qquad (3.5) $$

We have introduced the mass function $m(r)$, with $N(r) = 1 - 2m(r)/r$, and denote $T_{\hat{0}\hat{0}} = -T_0^0$ the projection of $T_{\mu\nu}$ onto the timelike vector of the orthonormal tetrad.

These equations admit simple analytical solutions describing Abelian configurations, for which the non-linear part of the SU(2) field strength vanishes.

### A. Reissner-Nordström solution

The simplest solution to Eqs. (3.4) is given by

$$ f = 0, \quad \phi = 1, \quad \sigma = 1, \quad N(r) = 1 - \frac{2M}{r} + \frac{Q^2}{r^2} \quad \text{with} \quad Q^2 = \frac{\kappa}{2} P^2, \quad P = \frac{n}{2e}, \qquad (3.6) $$

where $M$ is an integration constant. This describes the magnetic RN black hole,

$$ ds^2 = -N(r) dt^2 + \frac{dr^2}{N(r)} + r^2(d\vartheta^2 + \sin^2\vartheta d\varphi^2), \\ B = -(n/2) \cos\vartheta d\varphi, \quad W = T_3 B, \quad \Phi^T = (0, 1). \qquad (3.7) $$

The event horizon is located at $r_h = M + \sqrt{M^2 - Q^2}$. The gauge fields exhibit a Dirac string singularity at $\vartheta = 0, \pi$, which can be removed by gauge transformations (2.15) with

$$ U_{\pm} = \exp\left(\pm i\nu\varphi(1 + \tau_3)/2\right), \qquad (3.8) $$
whose effect is

$$B \rightarrow B_{\pm} = B \pm \nu d\varphi, \quad W \rightarrow W_{\pm} = T_3 B_{\pm}. \qquad (3.9)$$

The fields $W_-$, $B_-$ and $W_+$, $B_+$ should be used, respectively, in the northern and southern hemispheres, where they remain regular at the axis. These two locally regular gauges are related in the equatorial region by the gauge transformation

$$\mathcal{U} = \exp(i\nu\varphi(1 + \tau_3)), \qquad (3.10)$$

which is regular, as $2\nu = -n \in \mathbb{Z}$ [59]. However, for notational simplicity, we shall continue using the singular gauge as in (3.7).

We have

$$dB = \frac{n}{2} \sin\vartheta d\vartheta \wedge d\varphi, \qquad (3.11)$$

while the electromagnetic Nambu tensor $\mathcal{F}$ from (2.26) and the 't Hooft tensor $F$ from (2.36) are given by

$$\mathcal{F} = F = \frac{1}{e} dB = dA \quad \text{with} \quad A = \frac{1}{e} B \qquad (3.12)$$

(while $\mathcal{Z} = Z = 0$). This describes the radial hypermagnetic field $B$ and radial magnetic field $\mathcal{B}$,

$$B \equiv B^r = \frac{n}{2r^2}, \quad \mathcal{B} \equiv \mathcal{B}^r = \frac{B}{e} = \frac{P}{r^2}. \qquad (3.13)$$

The magnetic field is the identical to that of a Dirac monopole with magnetic charge

$$P = \frac{n}{2e}, \quad n \in \mathbb{Z}, \qquad (3.14)$$

in agreement with Dirac's quantization condition [63].

The RN solution (3.7) can thus be interpreted as the gravitational extension of the flat space Dirac monopole, to which it reduces in the $\kappa \rightarrow 0$ limit, where we can set $N(r) = 1$.

In what follows, we shall use the term “magnetic charge” to refer not only to $P$ but also to $Q$ defined in (3.6), as well as to the integer $n$.

The mass of the RN black hole is bounded below,

$$M \ge |Q|, \qquad (3.15)$$

where the bound is saturated in the extremal limit, where $r_h = M = |Q|$. In this limit, the function $N(r) = (1 - r_h/r)^2$ develops a double zero at the horizon, making it degenerate and causing the surface gravity to vanish.

The hairy black holes studied below approach the RN solution (3.7) in the far-field region. Remarkably, the inequality (3.15) is violated for extremal and near-extremal hairy black holes: we shall see later that their mass $M$ is smaller than their charge $|Q|$.
Bertotti-Robinson solution

The metric of the extremal RN solution, where $N = (1 - r_h/r)^2$, can be rewritten in terms of a new radial coordinate $\eta \in (-\infty, \infty)$ defined by

$$d\eta = \frac{dr}{r\sqrt{N}} \Rightarrow r = (1 + e^{\eta}) r_h. \qquad (3.16)$$

Rescaling the time coordinate by a constant factor, the line element takes the form

$$ds^2 = r_h^2 \left\{ -\frac{e^{2\eta}}{(1+e^{\eta})^2} dt^2 + (1+e^{\eta})^2 (d\eta^2 + d\vartheta^2 + \sin^2\vartheta d\varphi^2) \right\}. \qquad (3.17)$$

The horizon is at $\eta \to -\infty$, where the line element reduces to

$$ds^2 = r_h^2 \{-e^{2\eta} dt^2 + d\eta^2 + d\vartheta^2 + \sin^2\vartheta d\varphi^2\}. \qquad (3.18)$$

This is the Bertotti-Robinson metric [64, 65], a solution valid for all $\eta \in (-\infty, \infty)$, describing the direct product geometry $AdS_2 \times S^2$. The B, W, $\Phi$ fields remain as in (3.7), while the radial magnetic and hypermagnetic fields are given by (3.13) with $r = r_h = |Q|$.

B. Reissner-Nordström-de Sitter solution

Another simple solution of (3.4),(3.5) arises when the Higgs field vanishes:

$$f = 1, \quad \phi = 0, \quad \sigma = 1, \quad N = 1 - \frac{2M}{r} + \frac{g^2 Q^2}{r^2} - \frac{\Lambda}{3} r^2, \qquad (3.19)$$

with the same $Q$ as in (3.6) and

$$\Lambda = \frac{\kappa\beta}{8}. \qquad (3.20)$$

Here, the Higgs field remains in the false vacuum, and its potential acts as an effective cosmological constant, making the spacetime geometry RN-de Sitter (RNdS). The B-field remains identical to that of the RN solution, while the W-field is pure gauge. By performing a gauge transformation (Eq.(7.29) below), one can set $W = \Phi = 0$, making the full electroweak symmetry manifest everywhere.

Since $W_{\mu\nu}^a = 0$, Eq.(2.26) gives

$$\mathcal{F}_{\mu\nu} = \frac{g}{g'} B_{\mu\nu} = \frac{g^2}{e} B_{\mu\nu}. \qquad (3.21)$$

This corresponds to the radial hypermagnetic and magnetic fields:

$$B = \frac{n}{2r^2}, \quad \mathcal{B} = \frac{g^2 B}{e} = \frac{g^2 P}{r^2}. \qquad (3.22)$$

Thus, the magnetic charge is now $g^2 P = g^2 n/(2e) = gn/(2g')$, which does not satisfy the standard Dirac quantization condition. This is expected, as the system is not in the Higgs vacuum and standard electrodynamics does not apply.

However, the hypercharge field $B$ can still be made regular by removing the string singularity using the same gauge transformations as in the RN case. Notably, while the magnetic charge transforms as $P \to g^2 P$ compared to the RN case, in (3.19) we have $Q^2 \to g^2 Q^2$, not $Q^2 \to g^4 Q^2$.
Extremal limit

The function $N(r)$ in (3.19) has three positive roots, corresponding to the inner black hole horizon, the outer black hole horizon $r_h$, and the cosmological horizon $r_c$. The value of $r_h$ is bounded below as

$$r_h \ge r_{\text{ex}} = \sqrt{\frac{1 - \sqrt{1 - 4g^2 \Lambda Q^2}}{2\Lambda}} = g|Q| \left( 1 + \frac{g^2 Q^2 \Lambda}{2} + \dots \right) = g|Q| + \mathcal{O}(n^3 \kappa^{5/2}). \quad (3.23)$$

Furthermore, we have

$$1 - 2\Lambda r_{\text{ex}}^2 = \sqrt{1 - \frac{Q^2}{Q_{\text{m}}^2}} \equiv \lambda^2 \quad \text{with} \quad Q_{\text{m}} = \frac{1}{2g\sqrt{\Lambda}} \approx \frac{1.17}{\sqrt{\kappa}}. \qquad (3.24)$$

The lower bound $r_h = r_{\text{ex}}$ corresponds to the extremal limit, where the two black hole horizons merge, forming a degenerate horizon. In this case, the function $N(r)$ factorizes as

$$N(r) = k^2(r) \left(1 - \frac{r_{\text{ex}}}{r}\right)^2, \quad k^2(r) = 1 - \frac{\Lambda}{3} \left[r^2 + 2r_{\text{ex}}r + 3(r_{\text{ex}})^2\right]. \qquad (3.25)$$

After a coordinate transformation similar to (3.16), the horizon limit $r \to r_{\text{ex}}$ leads to an independent solution of the Einstein equations with an $AdS_2 \times S^2$ geometry, similar to (3.18):

$$ds^2 = r_{\text{ex}}^2 \{-e^{2\lambda\eta} dt^2 + d\eta^2 + d\vartheta^2 + \sin^2\vartheta d\varphi^2\}, \quad B = -\frac{n}{2} \cos\vartheta d\varphi, \quad W = \Phi = 0, \quad (3.26)$$

with $\eta \in (-\infty, \infty)$ and $\lambda$ defined in (3.24). A more detailed discussion of this will be provided in Appendix E.

The positive root of $k^2(r)$ determines the position of the cosmological horizon $r_c$, which is typically much larger than the black hole horizon $r_{\text{ex}}$ due to the smallness of $\kappa$:

$$r_c = \sqrt{\frac{3}{\Lambda} - 2r_{\text{ex}}^2} - r_{\text{ex}} = \sqrt{\frac{24}{\beta\kappa}} + \mathcal{O}(n\sqrt{\kappa}) \gg r_{\text{ex}} = \frac{\sqrt{\kappa|n|}}{2\sqrt{2}g'} + \mathcal{O}(n^3\kappa^{5/2}). \qquad (3.27)$$

However, as the magnetic charge increases and $Q = \sqrt{\kappa/2} P$ approaches $Q_m$, the horizons $r_c$ and $r_{\text{ex}}$ approach each other. In the $Q \to Q_m$ limit, $N(r)$ develops a triple zero, and one obtains

$$P = \frac{n}{2e} = \sqrt{\frac{2}{\kappa}} Q_{\text{m}} \Rightarrow |n| = \frac{4g'}{\sqrt{\beta\kappa}} = 2.60 \times 10^{32}, \quad r_c = r_{\text{ex}} = \frac{1}{\sqrt{2\Lambda}} = \frac{2}{\sqrt{\beta\kappa}}, \qquad (3.28)$$

which corresponds to the dimensionful radius $r_c = r_c l_0 = 3$ cm.

Further increasing the magnetic charge would lead to the appearance of a naked singularity. Therefore, the extremal RNdS black holes exist only for $|Q| \le Q_m$.

Although the extremal (and non-extremal) RNdS solution is not asymptotically flat, we will see in the following that its horizon limit is the same as that of the asymptotically flat extremal hairy black holes.
IV. PERTURBATIONS AROUND THE RN BLACK HOLE

The RN solution (3.7) encapsulates an important feature: when its magnetic field becomes sufficiently large at the horizon, this should trigger the electroweak condensation [47, 48], leading to an instability. As the condensate just begins to form, it can be described perturbatively and viewed as the “hair” that appears around the black hole.

Therefore, we now consider small perturbations around the RN solution (3.7):

$$g_{\mu\nu} \rightarrow g_{\mu\nu} + \delta g_{\mu\nu}, \quad W_{\mu}^{a} \rightarrow W_{\mu}^{a} + \delta W_{\mu}^{a}, \quad B_{\mu} \rightarrow B_{\mu} + \delta B_{\mu}, \quad \Phi \rightarrow \Phi + \delta \Phi. \quad (4.1)$$

Assuming the unitary gauge, where $\Phi^T = (0, \phi_2)$ with real-valued $\phi_2$, and linearizing the field equations with respect to the perturbations, the equation for

$$w_{\mu} = \frac{1}{g} (\delta W_{\mu}^{1} + i \delta W_{\mu}^{2}) \quad (4.2)$$

decouples from the rest and reduce to the complex Proca equation:

$$\mathcal{D}^{\mu} w_{\mu\nu} + i e F_{\nu\sigma} w^{\sigma} = m_{\mathrm{W}}^{2} w_{\nu}. \quad (4.3)$$

Here $w_{\mu\nu} = D_{\mu}w_{\nu} - D_{\nu}w_{\mu}$, with $D_{\mu} = \nabla_{\mu} + ieA_{\mu}$, and $A_{\mu} = B_{\mu}/e$ is the background field. One has $F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}$. We recall that $B_{\mu}dx^{\mu} = -(n/2)\cos\vartheta d\varphi$, where $n \in \mathbb{Z}$.

To separate the angular variables, we set

$$w_{\mu} dx^{\mu} = w_{\mathrm{m}} \quad \text{if} \quad n > 0, \quad \text{and} \quad w_{\mu} dx^{\mu} = \bar{w}_{\mathrm{m}} \quad \text{if} \quad n < 0, \quad (4.4)$$

where

$$w_{\mathrm{m}}(t, r, \vartheta, \varphi) = e^{i\omega t - im\varphi} \psi(r) (\sin\vartheta)^j \left(\tan\frac{\vartheta}{2}\right)^{\mathrm{m}} (d\vartheta - i \sin\vartheta d\varphi). \quad (4.5)$$

Here, the orbital angular momentum is

$$j = \left| \frac{n}{2} \right| - 1 \equiv j(n), \quad |n| \ge 2, \quad (4.6)$$

and $m \in [-j, j]$ is the azimuthal number. The values $n = \pm 1$ are excluded because $w_{\mathrm{m}}$ is then unbounded at $\vartheta = 0, \pi$, and the perturbation theory does not apply.

Injecting (4.5) into (4.3), the problem reduces to a Schrödinger-type equation

$$\left(-\frac{d^2}{dr_*^2} + N(r) \left[m_w^2 - \frac{|n|}{2r^2}\right]\right) \psi(r) = \omega^2 \psi(r). \quad (4.7)$$

Here, $N(r)$ is the metric coefficient in the RN solution (3.6), and the tortoise radial coordinate $r_*$ is defined by $dr_* = dr/N(r)$. If this equation admits negative modes, that is, bounded solutions with $\omega^2 < 0$, then perturbations grow in time, and the background is unstable.

In the flat space limit, where $N(r) = 1$ and the RN solution reduces to the Dirac monopole, the equation admits infinitely many bound states with $\omega^2 < 0$. Indeed, if $\omega = 0$ and $r \ll 1$, one can neglect $m_w^2$ in the potential, and the solution is

$$\psi(r) = \sqrt{r} \cos \left( \sqrt{2|n| - 1} \ln(\sqrt{r}) \right).$$
<table><caption>Table I. Values of $r_h^k$ for which Eq. (4.7) with $|n| = 2$ admits a k-node zero mode.</caption><thead><tr><th>k</th><th>r<sub>h</sub><sup>k</sup></th><th>k</th><th>r<sub>h</sub><sup>k</sup></th><th>k</th><th>r<sub>h</sub><sup>k</sup></th></tr></thead><tbody><tr><td>0</td><td>0.8942</td><td>5</td><td>2.2029 × 10<sup>-8</sup></td><td>10</td><td>3.3758 × 10<sup>-16</sup></td></tr><tr><td>1</td><td>4.306 × 10<sup>-2</sup></td><td>6</td><td>5.8553 × 10<sup>-10</sup></td><td>11</td><td>1.2763 × 10<sup>-16</sup></td></tr><tr><td>2</td><td>1.1723 × 10<sup>-3</sup></td><td>7</td><td>1.5563 × 10<sup>-11</sup></td><td>12</td><td>1.2379 × 10<sup>-16</sup></td></tr><tr><td>3</td><td>3.1180 × 10<sup>-5</sup></td><td>8</td><td>4.1367 × 10<sup>-13</sup></td><td>13</td><td>1.2369 × 10<sup>-16</sup></td></tr><tr><td>4</td><td>8.2879 × 10<sup>-7</sup></td><td>9</td><td>1.0997 × 10<sup>-14</sup></td><td>14</td><td>1.23688 × 10<sup>-16</sup></td></tr></tbody></table>

This solution oscillates infinitely many times as $r \to 0$, hence does not correspond to the ground state. Therefore, there should be eigenstates with $\omega^2 < 0$. It follows that the Dirac monopoles are unstable with respect to developing the electroweak condensate [58], except for $n = \pm 1$, when the perturbation theory does not apply.

In curved space, the range of the radial coordinate $r$ in (4.7) is restricted to the outer black hole region, $r \ge r_h$, which bounds the attractive term $-|n|/(2r^2)$ in the potential. As a result, large black holes are stable, while the small ones exhibit instabilities characterized by negative modes – normalizable solutions of (4.7) with $\omega^2 < 0$.

By varying the black hole size, one can identify the threshold radius $r_h = r_h^0(n)$, for which the first instability just starts to appear as a zero mode – normalizable solution of (4.7) with $\omega = 0$. This should occur when the hypermagnetic field at the horizon becomes sufficiently strong to trigger the electroweak condensation, leading to the condition

$$|B(r_h^0)| = \frac{|n|}{2(r_h^0)^2} \approx m_W^2 = \frac{g^2}{2} \quad \Rightarrow \quad r_h^0 \approx \frac{\sqrt{|n|}}{g} = 1.13 \times \sqrt{|n|}. \qquad (4.8)$$

Numerical solutions to Eq. (4.7) for small $|n|$ yield slightly smaller values of $r_h^0(n)$ (see Table II below). This discrepancy is explained by the fact that the estimate $|B| = m_W^2$ is derived for a field orthogonal to a plane, not to a spherical horizon. However, when $|n| \gg 1$, the horizon becomes large and locally approximates a plane. In this regime, the numerical results recover the expected scaling $r_h^0(n) = \sqrt{|n|}/g + \dots$, where the dots denote terms that vanish as $|n| \to \infty$.

For $|n| = 2$, where $j(n) = 0$, the perturbations in (4.5) are spherically symmetric. Starting from a large value of $r_h$ and gradually decreasing it, the first zero mode appears for $r_h^0(2) \approx 0.89$, with the corresponding wavefunction $\psi_0(r)$ being *nodeless*. As $r_h$ decreases further, this solution becomes a negative mode, with $\omega_0^2 < 0$. However, when $r_h = r_h^1(2) \approx 0.04$, the equation develops a second zero mode $\psi_1(r)$, which has *one node*.

As $r_h$ continues to decrease, $\psi_1(r)$ also becomes a negative mode, alongside the *nodeless* $\psi_0(r)$. Further modes with an increasing number of nodes will continue to develop as the black hole size shrinks. In practice, we obtain these bound state solutions by applying the numerical multiple shooting method (see [41] for details).

Reducing $r_h$ all the way down to the minimal value $r_h = |Q|$, Eq. (4.7) develops a total of 15 bound states. The values $r_h = r_h^k$ for which the k-node zero mode appears are shown in Table I. For $r_h < r_h^k$, each zero mode becomes negative, and the number of instabilities increases. For example, at $r_h = 10^{-5}$, there are four negative modes; their profiles $\psi_k(r_\star)$ and the values of $|\omega_k(r_h)| \times r_h$ for $k = 0, 1, 2$ are shown in Fig. 1.

The eigenvalues $\omega_k^2(r_h) \le 0$ start from zero at $r_h = r_h^k$ and then decrease as $r_h$ gets smaller. However, the product $|\omega| r_h$ rapidly approaches a constant value.
<figure><img src="image_3.png" /><figcaption>Figure 1. Left: the rescaled eigenvalue $|\omega|r_h$ against $r_h$ for the bound state solutions of (4.7) with $k = 0, 1, 2$ nodes. Right: the profiles of the bound state eigenfunctions $\psi_k(r_*)$ with $k = 0, 1, 2$ nodes against the radial tortoise coordinate $r_*$ for $r_h = 10^{-5}$.</figcaption></figure>

<table><caption>Table II. Values $r_h^k(n)$ for which Eq.(4.7) admits a k-node zero mode.</caption><thead><tr><th>n</th><th>2</th><th>4</th><th>6</th><th>10</th><th>20</th><th>40</th><th>100</th><th>200</th></tr></thead><tbody><tr><th>r<sub>h</sub><sup>0</sup>(n)</th><td>0.894</td><td>1.465</td><td>1.927</td><td>2.682</td><td>4.109</td><td>6.164</td><td>10.290</td><td>14.966</td></tr><tr><th>r<sub>h</sub><sup>1</sup>(n)</th><td>0.043</td><td>0.245</td><td>0.508</td><td>1.054</td><td>2.274</td><td>4.191</td><td>8.202</td><td>12.823</td></tr><tr><th>r<sub>h</sub><sup>2</sup>(n)</th><td>0.001</td><td>0.024</td><td>0.089</td><td>0.306</td><td>1.066</td><td>2.627</td><td>6.342</td><td>10.829</td></tr></tbody></table>

In summary, small RN black holes with $|n| = 2$ are unstable with respect to the formation of an electroweak condensate. Each new instability manifests as a zero mode $\psi_k(r)$ at $r_h = r_h^k(2)$, which represents a static deformation of the RN black hole by the condensate that begins to form. This can be viewed as a perturbative approximation of a new static hairy solution that merges with the RN family for $r_h = r_h^k$ but deviates from it for $r_h < r_h^k$. The new solutions are spherically symmetric since $j(n) = 0$ for $|n| = 2$.

The same arguments apply for $|n| > 2$, where $j(n) > 0$ and the solutions (4.5) are no longer spherically symmetric. A zero mode $\psi_k(r)$ that appears for $r_h = r_h^k(n)$ then describes a static deformation of the RN black hole by a non-spherically symmetric condensate. This deformation generalizes to a static, non-perturbative hairy solution for $r_h < r_h^k(n)$. The values of $r_h^k(n)$ for $k = 0, 1, 2$ are given in Table II and shown in Fig.2. The profiles of the nodeless zero mode $\psi_0$ in Fig.2 show that the solution becomes increasingly localized at the horizon when $n$ increases.

It is worth mentioning the analytical formula obtained by solving Eq.(4.7) within the WKB approximation [66]:

$$r_h^k(n) = \frac{\sqrt{|n|}}{g} \left( 1 - \frac{2}{\sqrt{|n|}} \left(k + \frac{1}{2}\right) \right), \quad k = 0, 1, 2, \dots, \quad |n| \gg 1. \quad (4.9)$$

For $n = 200$, this gives values of $r_h^k(n)$ that agree with the numerically obtained values in the last column of Table II, with relative discrepancies of 0.4%, 1%, and 4% for $k = 0, 1, 2$, respectively. The agreement is reasonable.
<figure><img src="image_4.png" /><figcaption>Figure 2. The RN horizon radius $r_h^k(n)$, for which there exists a k-node zero mode $\psi_k(r)$ (left), and the profile of the fundamental zero mode $\psi_0(r)$ (right).</figcaption></figure>

To recapitulate, Eqs.(4.7) admits bound state solutions with $\omega = 0$ for discrete values $r_h^k(n)$ of the black hole size. These zero modes can be used for a perturbative approximation of new hairy solutions that exist for $r_h < r_h^k(n)$. The hairy solutions are spherically symmetric only for $|n| = 2$, while for $|n| > 2$, the spherical symmetry is lost.

## V. SPHERICALLY SYMMETRIC HAIRY SOLUTIONS FOR $n = \pm 2$

As shown above, the RN solution (3.7) with $|n| = 2$ becomes unstable with respect to the formation of a spherically symmetric condensate when the horizon radius is small. This instability is expected to lead to the emergence of more general static and spherically symmetric non-Abelian solutions describing hairy black holes. Such solutions must satisfy Eqs.(3.4).

Notably, unless for $|n| = 2$, Eq.(3.4e) requires that $f' = f\phi = 0$, which only permits the simplest Abelian solutions described above. However, for $|n| = 2$, we have $\nu^2 - 1 = 0$, so Eq.(3.4e) is automatically satisfied, allowing for solutions with non-constant $f(r)$ and $\phi(r)$.

Such solutions were previously reported in Ref.[43]; below, we reproduce them and highlight their most important features. Although the two signs of $n = \pm 2$ correspond to two values of the magnetic charge $P = n/(2e) = \pm 1/e$, the radial amplitudes $N(r), \sigma(r), f(r), \phi(r)$ are the same in both cases.

Eq.(3.4d) can be solved as

$$ \sigma(r) = \exp \left( -\kappa \int_r^\infty (\mathcal{U}_1/r) dr \right), \qquad (5.1) $$

using which the remaining Eqs.(3.4a)-(3.4c) reduce to

$$ (N f')' + \frac{\kappa}{r} \mathcal{U}_1 N f' = \left( \frac{f^2 - 1}{r^2} + \frac{g^2}{2} \phi^2 \right) f, \qquad (5.2a) $$

$$ (r^2 N \phi')' + \kappa r \mathcal{U}_1 N \phi' = \left( \frac{\beta r^2}{4} (\phi^2 - 1) + \frac{1}{2} f^2 \right) \phi, \qquad (5.2b) $$

$$ 1 - (rN)' = \kappa \left( N\mathcal{U}_1 + \mathcal{U}_0 + \frac{1}{2g'^2r^2} \right). \qquad (5.2c) $$
The simplest solution is the RN for which $f = 0$, $\phi = 1$, and $N$ is given by (3.6). Perturbing the solution as $f \to \delta f$, $\phi \to 1 + \delta\phi$, and $N \to N + \delta N$, then linearizing with respect to the static perturbations, yields

$$ (N\delta f')' = \left(-\frac{1}{r^2} + m_w^2\right) \delta f, \quad \frac{1}{r^2} (r^2 N \delta\phi')' = m_H^2 \delta\phi, \quad (r\delta N)' = 0. \quad (5.3) $$

The first of these equations is precisely Eq. (4.7) with $|n| = 2$ and $\omega = 0$. As we know, this equation admits a non-trivial solution for $r_h = r_h^0(2) = 0.894$, and we expect this perturbative zero mode to generalize to a non-perturbative solution of Eqs. (5.2) when $r_h < 0.894$.

Assuming that at the horizon $r = r_h$, we have $N(r_h) = 0$ and $N'(r_h) \equiv N'_h > 0$, the equations imply

$$ \begin{aligned} N'_h &= \frac{1}{r_h} \left( 1 - \kappa \mathcal{U}_{0h} - \frac{\kappa}{2g'^2 r_h^2} \right), & f'_h &= \frac{1}{N'_h} \left( \frac{f_h^2 - 1}{r_h^2} + \frac{g^2}{2} \phi_h^2 \right) f_h, \\ \phi'_h &= \frac{1}{r_h^2 N'_h} \left( \frac{\beta r_h^2}{4} (\phi_h^2 - 1) + \frac{1}{2} f_h^2 \right) \phi_h, \end{aligned} \quad (5.4) $$

where the subscript $h$ denotes the value at $r = r_h$. This determines the boundary conditions at the horizon in terms of $f_h$ and $\phi_h$.

For $r \to \infty$, the perturbations should vanish, and the solutions should approach the RN form. Thus, we have

$$ N = 1 - \frac{2M}{r} + \delta N_N, \quad f = C_f e^{-m_W r} + \delta f_N, \quad \phi = 1 + \frac{C_\phi}{r} e^{-m_H r} + \delta \phi_N. \quad (5.5) $$

Here, $M, C_f, C_\phi$ are integration constants corresponding to the massless Newtonian mode and the two massive linear modes described by (5.3), while $\delta N_N, \delta f_N, \delta \phi_N$ contain the subleading terms and nonlinear corrections. These expressions determine the boundary conditions at large $r$. The values of $\delta N_N, \delta f_N, \delta \phi_N$ can be obtained by converting (5.2) into integral equations and iterating (see [41] for details).

Integrating the equations with the above boundary conditions and adjusting the values of the 5 free parameters $-f_h, \phi_h, M, C_f, C_\phi$ within the numerical multiple shooting method [67] yields global solutions in the region $r > r_h$.

There always exists the trivial RN solution with $f(r) = 0, \phi(r) = 1$, but for $r_h < 0.894$, we also find solutions with non-constant $f(r)$ and $\phi(r)$ that increasingly deviate from the RN solution as $r_h$ decreases.

The horizon radius $r_h$ of these hairy black holes ranges within the limits $r_h \in [r_{ex}, r_h^0 = 0.894]$. The maximal value $r_h^0$ corresponds to the bifurcation with the RN solution, while the minimal value $r_{ex}$ determines the extremal solution for which $N'_h = 0$.

Notably, the formulas (5.4) remain valid even as $N'_h \to 0$, provided that $f_h \to 1, \phi_h \to 0, f'_h \to 0, \phi'_h \to 0$, and $r_h \to r_{ex}$. The value $r_{ex}$ is determined as the root of the equation $N'_h = 0$ in (5.4). Remarkably, it coincides with the radius $r_{ex} \approx g|Q|$ of the extremal RNdS solution in (3.23).

For $|n| = 2$ one finds $r_{ex} \approx 1.09 \times 10^{-16}$, which is close to the Planck length.

Examples of hairy solutions for $r_h = 0.8, r_h = 0.24$, and the extremal case with $r_h = r_{ex}$ are shown in Fig.3. The function $\sigma(r)$ (not shown in the figures) is determined by (5.1) and remains very close to unity since $\kappa$ is very small. As $r_h$ shrinks, the horizon value
<figure><img src="image_5.png" /><figcaption>Figure 3. Profiles of the hairy solutions with $r_h = 0.8$ (left), $r_h = 0.24$ (center), and the extremal solution with $r_h = r_{\text{ex}} \approx 1.09 \times 10^{-16}$ (right). In the extremal case, $N(r_h) = N'(r_h) = 0$.</figcaption></figure>

of the Higgs field $\phi_h$ decreases and approaches zero, while $f_h$ increases and approaches unity, as illustrated in Figs.3,4. This is a manifestation of symmetry restoration in the presence of a strong magnetic field [68]. When $r_h$ is small, the U(1) hypermagnetic field, $|B| = n/(2r^2) = \pm 1/r^2$ (as in (3.22)), becomes very large at the horizon. As a result, the value of $\phi_h$ is driven to zero, along with the value of the SU(2) field strength, which is proportional to $(1 - f_h^2)$. As $f_h$ increases, the SU(2) component of the magnetic charge shifts outside the horizon. This phenomenon can be understood as follows.

## A. The charge

Both the hairy black holes and the RN black hole share the same total magnetic charge measured at infinity, given by (2.50):

$$P = \frac{1}{4\pi} \oint_{S^2} \mathcal{F} = \frac{n}{2e} = \pm \frac{1}{e}. \qquad (5.6)$$

In the RN case, this charge is entirely confined inside the horizon. For the hairy black hole, the charge remains inside the horizon within the 't Hooft description. However, in the Nambu description, the charge splits as $P = P_H + P_h$, where $P_H$ is contained inside the horizon, while $P_h$ is distributed outside.

The exterior charge $P_h$ is determined by integrating the magnetic charge density (2.49):

$$ \begin{aligned} P_h &= \int_{r>r_h} \tilde{\mathcal{J}}^0 \sqrt{-g} d^3 x \equiv \int_{r_h}^{\infty} \rho_{\text{SU}(2)}(r) dr = \frac{e}{8\pi} \int_{r>r_h} \epsilon^{ijk} \partial_i \psi_{jk} d^3 x \\ &= -\frac{e}{8\pi} \oint_{r=r_h} \psi_{jk} dx^j \wedge dx^k = -\frac{e}{4\pi} \oint_{r=r_h} \psi_{\vartheta\varphi} d\vartheta d\varphi = f_h^2 g'^2 P. \end{aligned} \qquad (5.7) $$

Here, we used the fact that the fields (3.2) are already in the unitary gauge. Thus, Eq. (2.38) implies that

$$g^2 \psi_{\vartheta\varphi} = W_{\vartheta}^1 W_{\varphi}^2 - W_{\vartheta}^2 W_{\varphi}^1 = \nu f^2(r) \sin\vartheta, \qquad (5.8)$$

which decays exponentially fast at infinity. The same result can be derived using the relation $P_h = P - P(V)$, where $P(V)$ is given by (2.51) and $V$ is the volume for $r < r_h$.
<figure><img src="image_6.png" /><figcaption>Figure 4. Left: the horizon values $\phi_h = \phi(r_h)$ and $f_h = f(r_h)$ for the fundamental hairy solutions. Center: the mass of the hairy and RN solutions, with the former being always less energetic than the latter. Right: the radial densities of the magnetic charge and of the energy stored in the hair.</figcaption></figure>

Since $f_h$ increases when $r_h$ decreases, the hair charge $P_h$ grows from zero to its maximal value of $g'^2P$ in the extremal limit. On the other hand, the horizon charge $P_H$ decreases from $P$ to $g^2P$ in the extremal limit.

To summarize, within the Nambu description, the charged condensate outside the black hole horizon carries 22% of the total charge in the extremal limit, $P_h = g'^2P \approx 0.22P$. The radial profile of the condensate charge density $\rho_{SU(2)}$, as defined in Eq. (5.7), exhibits a peak at $r \approx 1$, as shown in Fig. 4.

## B. The ADM mass

The ADM mass is defined by Eq. (2.23), where $M = m(\infty)$ is the asymptotic value of the mass function appearing in the relation $N(r) = 1 - 2m(r)/r$. The mass function $m(r)$ satisfies Eq. (3.4c), given by

$$m'(r) = \frac{\kappa}{2} r^2 T_{00} = \frac{\kappa}{2} \left( \frac{1}{2g'^2 r^2} + N U_1 + U_0 \right), \qquad (5.9)$$

and integrating this with the boundary condition $N(r_h) = 0$ yields

$$M = m(\infty) = \frac{r_h}{2} + \frac{\kappa}{4g'^2 r_h} + \frac{\kappa}{2} \int_{r_h}^{\infty} (N U_1 + U_0) dr. \qquad (5.10)$$

For the RN solution (3.6), we have $f = 0$, $\phi = 1$, and $Q^2 = \kappa/(2e^2)$ for $n = \pm 2$. Using this in (5.10) yields the mass of the RN solution:

$$M_{\text{RN}} = \frac{1}{2} \left( r_h + \frac{Q^2}{r_h} \right). \qquad (5.11)$$

As shown in the central panel of Fig. 4, the mass $M$ of hairy solutions is always below $M_{\text{RN}}$ for a given $r_h$, so that hairy black holes are less energetic than the RN black holes of same size.
As explained above, hairy black holes contain the charge $P_h = f_h^2 g'^2 P$ in the hair and $P_H = P - P_h$ inside the horizon. This allows one to split the black hole mass into the “horizon mass” $M_H$ and the “hair mass” $M_h$,

$$M = \frac{1}{2} \left( r_h + (g^2 + (1 - f_h^2)^2 g'^2) \frac{Q^2}{r_h} \right) + M_h \equiv M_H + M_h, \qquad (5.12)$$

where $M_H$ is the mass of the RN black hole with charge $P_H$ and size $r_h$. It turns out that $M_H \gg M_h$, so $M \approx M_H$ with very high precision.

The formula (5.12) will be justified below within a more general context. At present, we simply note that in the extremal limit, where $r_h = r_{ex} \approx g|Q|$ and $f_h = 1$, the horizon mass $M_H$ in (5.12) reduces to the sum of the first two terms on the right in (5.10). Hence, $M = M_H + M_h$, with

$$M_H = \frac{1}{2} \left( r_{ex} + \frac{g^2 Q^2}{r_{ex}} \right) \approx g|Q|, \qquad M_h = \frac{\kappa}{2} \int_{r_h}^{\infty} (N\mathcal{U}_1 + \mathcal{U}_0) dr. \qquad (5.13)$$

Here, the hair mass $M_h$ is determined by the non-linear $f, \phi$ fields forming the condensate, so this represents the mass of the condensate. It is instructive to consider their values rescaled according to (2.25):

$$\mathcal{M}_H = \frac{8\pi}{\kappa} M_H = \frac{8\pi g}{\sqrt{2\kappa}e} + \mathcal{O}(\kappa^{7/2}) = 5 \times 10^{17},$$

$$\mathcal{M}_h = \frac{8\pi}{\kappa} M_h = 4\pi \int_{r_{ex}}^{\infty} (N\mathcal{U}_1 + \mathcal{U}_0) dr \equiv \int_{r_{ex}}^{\infty} \varepsilon_h(r) dr = 15.759. \qquad (5.14)$$

The value of $\mathcal{M}_h$ is obtained from our numerical solution. As a result, the hair mass is negligible as compared to the horizon mass, with $\mathcal{M}_h/\mathcal{M}_H \approx 3 \times 10^{-17}$, meaning $M \approx M_H$.

The latter relation has an important consequence: although the mass of the extremal RN black hole is equal to the charge,

$$M_{\text{RN}} = |Q|, \qquad (5.15)$$

the mass of the extremal hairy solution is smaller than the charge:

$$M = M_H + M_h \approx M_H \approx g|Q| \approx 0.88 \times |Q| < |Q|. \qquad (5.16)$$

Therefore, the hairy black hole is energetically favoured compared to the RN black hole. This remarkable feature can be explained by recalling that the Zeeman energy is negative. The W bosons contained in the condensate interact with the magnetic field of the black hole, which shifts their mass as $m_W^2 \rightarrow m_W^2 - |B|$. This renders the condensate almost massless. Although the condensate carries 22% of the total charge, i.e., $Q_h \equiv \sqrt{\kappa/2} P_h = 0.22 \times Q$, its relative contribution to the total mass is only $3 \times 10^{-17}$. Therefore, the mass-to-charge ratio for the condensate is $M_h/|Q_h| \sim 10^{-16}$.

The radial energy density of the condensate, $\varepsilon_h(r)$, defined in (5.14), has a sharp peak outside the horizon at $r \approx 1$, where the charge density $\rho_{SU(2)}(r)$ is also maximal, as shown in Fig.4. The condensate is gravitationally attracted by the black hole but repelled magnetically, as the horizon charge $P_H$ and the condensate charge $P_h$ have the same sign. Since $M_h/|Q_h| \sim 10^{-16}$, the magnetic repulsion is 16 orders of magnitude stronger than gravitational attraction.
This can be viewed as a manifestation of the weak gravity conjecture [69]. However, the condensate cannot be pushed away to infinity because it is composed of massive fields and encounters the Yukawa potential barrier at large $r$.

Since $M \approx M_H$, the dimensionful mass and horizon size are given by

$$ M = \frac{e^2}{4\pi\alpha} M \times m_0 \approx \frac{g}{\sqrt{\alpha}} M_{\text{Pl}} = 10.32 \times M_{\text{Pl}}, \quad r_{\text{ex}} = r_{\text{ex}} l_0 \approx \frac{g}{\sqrt{\alpha}} L_{\text{Pl}}. \quad (5.17) $$

For the extremal RN black hole with the same magnetic charge, the mass and radius are obtained by dividing by $g$: $M = M_{\text{Pl}}/\sqrt{\alpha} = 11.71 \times M_{\text{Pl}}$ and $r_{\text{ex}} = L_{\text{Pl}}/\sqrt{\alpha}$.

### C. Extremal hairy black hole versus Cho-Maison monopole

For the extremal hairy black holes, one has the horizon values $f_h = 1$ and $\phi_h = 0$, so that the fields correspond to the extremal RNdS solution (3.19),(3.25). Therefore, the function $N(r)$ for the extremal hairy black hole has the same form as in (3.25),

$$ N(r) = k^2(r) \left(1 - \frac{r_{\text{ex}}}{r}\right)^2 \approx \left(1 - \frac{r_{\text{ex}}}{r}\right)^2. \quad (5.18) $$

Here, $r_{\text{ex}} = g|Q| + \mathcal{O}(\kappa^{5/2}) \sim 10^{-16}$ is determined by (3.23), and $k(r)$ interpolates between the same horizon limit as in (3.25), $k(r_{\text{ex}}) = \sqrt{1 - 2\Lambda r_{\text{ex}}^2} = 1 + \mathcal{O}(\kappa^2)$, and $k(\infty) = 1$. Hence, $k(r) \approx 1$. The function $N(r)$ quickly approaches unity, and already for $r \simeq 10^2 \times r_{\text{ex}} \sim 10^{-14}$, the geometry is practically flat. As a result, the massive fields described by $f(r)$ and $\phi(r)$ actually live in a flat geometry, and their profiles practically coincide with those for the flat space Cho-Maison monopole.

The Cho-Maison (CM) monopole [70] is the solution in the flat space limit, when $\kappa = 0$ and $N = \sigma = 1$, and Eqs. (5.2) reduce to

$$ f'' = \left( \frac{f^2 - 1}{r^2} + \frac{g^2}{2} \phi^2 \right) f, \quad (r^2 \phi')' = \left( \frac{\beta r^2}{4} (\phi^2 - 1) + \frac{1}{2} f^2 \right) \phi. \quad (5.19) $$

These equations admit a smooth solution with the boundary conditions for $0 \leftarrow r \rightarrow \infty$:

$$ 1 + \mathcal{O}(r^2) \leftarrow f(r) \rightarrow \mathcal{O}(e^{-m_W r}), \quad \mathcal{O}\left(r^{\frac{\sqrt{3}-1}{2}}\right) \leftarrow \phi(r) \rightarrow \mathcal{O}(e^{-m_H r}). \quad (5.20) $$

This describes a magnetic monopole with the same magnetic charge $P = \pm 1/e$ as for the hairy black holes. The SU(2) part of the charge $g'^2 P = 0.22 \times P$ is smoothly distributed in space in the condensate, while the U(1) part $g^2 P = 0.78 \times P$ is pointlike and located at the origin. The latter carries the energy density $1/(2g'^2r^4)$, which renders the total monopole energy $E_{\text{mon}}$ divergent. However, subtracting the divergent part leaves a finite value:

$$ E_{\text{CM}} = E_{\text{mon}} - 4\pi \int_{0}^{\infty} \frac{dr}{2g'^2 r^2} = 15.759, $$

(see [59] for details). This is the same value as that in Eq. (5.14), where it corresponds to the energy of the condensate around the black hole.

The profiles of $f(r)$ and $\phi(r)$ for the extremal black hole in the region $r > r_{\text{ex}}$ are very close to those for the CM monopole and converge to the latter pointwise in the $\kappa \rightarrow 0$ limit.
<figure><img src="image_7.png" /><figcaption>Figure 5. The amplitudes $f(r)$ and $\phi(r)$ for the fundamental hairy solution with $r_h = 10^{-6}$ and for its first two radial excitations.</figcaption></figure>

Therefore, the extremal hairy black hole can be viewed as the CM monopole harbouring in its center a tiny black hole whose radius is microscopically small compared to the size of the monopole itself.

As seen in Fig.3 (right panel), the near-horizon region where the geometry essentially deviates from Minkowski space, $10^{-16} \le r \le 10^{-14}$, is parametrically small compared to the region containing the condensate, where $r \simeq 1$. The amplitudes $f(r)$ and $\phi(r)$ stay close to their horizon values $f_h = 1, \phi_h = 0$ in the interval that exceeds the horizon size by approximately 12 orders of magnitude, and only for $r \ge 10^{-4}$ they start to change. Therefore, the massive fields forming the condensate actually live in flat geometry and almost do not feel the presence of the black hole. However, the black hole is important because it regularizes the divergent part of the monopole energy by replacing it with the finite term $M_H$ in (5.14), so that the total energy becomes finite.

The hair mass $M_h$ of the extremal black hole is given by the integral in (5.14), where the integration starts at $r = r_{ex} \sim \sqrt{\kappa} \sim 10^{-16}$. The regularized energy of the Cho-Maison monopole $E_{CM}$ is given by the same integral, but the integration starts at $r = 0$, hence $M_h = E_{CM} + O(\sqrt{\kappa})$.

To recapitulate, the extremal hairy solution can be viewed as a superposition of the flat space CM monopole of size $\sim 1$ and a black hole of radius $\sim 10^{-16}$. The tiny black hole almost does not affect the monopole configuration but renders its energy $E_{mon}$ finite by providing the cutoff at the horizon. This can be illustrated as follows:

$$E_{mon} = \lim_{\kappa \to 0} \frac{8\pi g}{\sqrt{2\kappa e}} + 15.759 = \infty \Rightarrow M = \frac{8\pi g}{\sqrt{2\kappa e}} + 15.759 = 5 \times 10^{17} + 15.759. \quad (5.21)$$

When the gravity is taken into account, the divergent part of the monopole energy becomes finite, albeit large, and reduces to the black hole horizon mass $M_H = 8\pi g/(\sqrt{2\kappa e})$. The finite part of the monopole energy $E_{CM} = 15.759$ remains the same up to corrections of order $10^{-16}$ and becomes the mass of the black hole hair $M_h = 15.759$.

## D. Stability of the extremal hairy solution

The extremal hairy black hole with $n = \pm 2$ is expected to be stable. This stability is initially suggested by its similarity with the Cho-Maison monopole in the flat space limit,
which has been proven to be stable [58]. Although a full stability analysis incorporating gravity is not yet available, the similarity of the electroweak fields with and without gravity suggests that the results should align.

Next, since the condensate significantly reduces the black hole mass, the hairy black hole with $M \approx 0.88 |Q|$ cannot transition to a RN black hole, for which $M \ge |Q|$. This would require either to increase the energy or to reduce the magnetic charge, but the latter is conserved due to the field topology. It is conceivable that the hairy black hole could decay into a lower-energy state. However, the perturbative analysis in the next Section suggests that the spherical condensate with $n = \pm 2$ already represents the absolute energy minimum.

Therefore, the extremal hairy black hole should be stable.

### E. The excited solutions

For the sake of completeness, it is important to mention also the excited zero mode solutions of Eq. (4.7), $\psi(r) = \psi_k(r)$ with $k > 0$, which exist for specific horizon values $r_h = r_h^k(n)$ as shown in Tables I and II. These excited solutions can also be promoted to static hairy black holes at the non-perturbative level.

As an example, consider the one-node zero mode solution $\psi_1(r)$, which exists for $r_h^1(2) = 4.3 \times 10^{-2}$. When this solution is used as the input for the numerical solver and the horizon radius is set to $r_h = 4 \times 10^{-2} < r_h^1(2)$, we find a static hairy solution for which $f(r)$ is not monotonic, but instead exhibits a node. Additionally, for the same value of $r_h$, there is also the black hole solution with a monotonic $f(r)$.

Thus, for $r_h = 4 \times 10^{-2}$, there are two distinct hairy solutions: one fundamental and one excited. Similarly, for $r_h < r_h^2(2) = 1.1 \times 10^{-3}$, a second radial excitation is found, and so on.

The excited solutions have not been described in the literature. Their profiles for $r_h = 10^{-6}$ are shown in Fig. 5. However, such solutions are presumably unstable, and we will not discuss them anymore.

## VI. PERTURBATIVE DESCRIPTION OF HAIRY BLACK HOLES WITH $|n| > 2$

We considered above spherically symmetric hairy solutions with magnetic charge $n = \pm 2$, beginning with their perturbative description in terms of Eqs. (4.3)–(4.7), and later passing to the fully non-perturbative description. Now, we aim to extend the analysis to cases where $n > |2|$, when spherical symmetry is no longer preserved. As the first step, we will apply the perturbative approach.

Let us get back to Eqs. (4.3)–(4.7), which describe perturbations of the RN black hole. We choose the horizon radius $r_h = r_h^0(n)$ such that Eq. (4.7) admits a zero mode solution $\psi(r) = \psi_0(r)$. Substituting this into (4.5) yields a static solution of the Proca equation (4.3), which describes the electroweak condensate beginning to form on the RN background:

$$w_{\mu}dx^{\mu} = w_{\mathrm{m}}(r, \vartheta, \varphi) = e^{-im\varphi}\psi(r)(\sin\vartheta)^{j}\left(\tan\frac{\vartheta}{2}\right)^{m}(d\vartheta - i\sin\vartheta d\varphi). \quad (6.1)$$

It is assumed here that $n > 0$, while for $n < 0$, one should replace this expression by its complex conjugate. We have $j = |n/2| - 1$. This solution is not unique since the azimuthal
number $m \in [-j, j]$ can take several values. Thus, the general solution is given by:

$$w_{\mu}dx^{\mu} = \sum_{m} c_{m}w_{m}(r, \vartheta, \varphi), \qquad (6.2)$$

where $c_m$ are constant coefficients. This solution is defined up to an overall sign, which can be flipped by the gauge transformation (2.15) generated by $U = -\tau_3$. The effect of this transformation is $T_1 \rightarrow -T_1$ and $T_2 \rightarrow -T_2$, resulting in $w \rightarrow -w$, while the RN background remains unchanged.

Since $w$ must be a single-valued function of $\varphi$ modulo a gauge transformation, it should either be invariant under $\varphi \rightarrow \varphi + 2\pi$ or change sign. Therefore, since $w_m \sim \exp(-im\varphi)$, the azimuthal number $m$ must take either integer or half-integer values, regardless of the value of $j = |n|/2 - 1$, which can itself be either integer or half-integer.

Additionally, one must have $|m| \le j$, since otherwise the mode (6.1) would diverge at $\vartheta = 0$ or $\vartheta = \pi$. Consequently, the allowed values for $m$ are as follows:

* if $j = 0$, then $m = 0$;

* if $j = \frac{1}{2}$, then $m = \pm \frac{1}{2}$ or $m = 0$;

* if $j = 1$, then $m = \pm \frac{1}{2}$ or $m = 0, \pm 1$;

* if $j = \frac{3}{2}$, then $m = \pm \frac{3}{2}, \pm \frac{1}{2}$ or $m = 0, \pm 1$, etc.

Therefore, $m$ assumes either $2j$ or $2j + 1$ distinct values, among which $m = 0$ is always present.

If $|n| = 2$, then $j = m = 0$, and there is only one coefficient $c_0$ in (6.2), which can be fixed by normalization, hence there is no ambiguity. For $|n| > 2$, $j > 0$, spherical symmetry is lost, and the azimuthal number $m$ can take several values. As a result, there are several coefficients $c_m$ in (6.2). Within the framework of linear perturbation theory, these coefficients are arbitrary, but higher-order perturbation theory imposes conditions on them. Specifically, the choice of these coefficients should minimize the energy of the condensate.

In what follows, we adopt the approach of [45], which involves expanding around the RN solution up to the second order terms and then minimizing the energy. We neglect the backreaction of perturbations on the RN geometry, as it is proportional to $\kappa$.

## A. Electroweak equations in the unitary gauge

Switching to the unitary gauge, where $\Phi^T = (0, \phi)$ with $\phi$ being real-valued, one defines the complex-valued W field and its strength:

$$w_{\mu} = \frac{1}{g}(W_{\mu}^{1} + iW_{\mu}^{2}), \quad w_{\mu\nu} = D_{\mu}w_{\nu} - D_{\nu}w_{\mu}, \qquad (6.3)$$

where $D_{\mu} = \nabla_{\mu} + iW_{\mu}^3$. The electroweak Lagrangian in (2.13) then takes the following form:

$$-\mathcal{L}_{\text{EW}} = \frac{1}{4g'^2}(B_{\mu\nu})^2 + \frac{1}{4g^2}(W_{\mu\nu}^3)^2 + \frac{1}{4}|w_{\mu\nu}|^2 \qquad (6.4) \\ +(\partial_{\mu}\phi)^2 + \frac{1}{4}(g^2|w_{\mu}|^2 + (B_{\mu} - W_{\mu}^3)^2)\phi^2 + \frac{\beta}{8}(\phi^2 - 1)^2,$$

with $B_{\mu\nu} = \partial_{\mu}B_{\nu} - \partial_{\nu}B_{\mu}$, where

$$W_{\mu\nu}^3 = \partial_{\mu}W_{\nu}^3 - \partial_{\nu}W_{\mu}^3 + g^2\psi_{\mu\nu}, \quad \psi_{\mu\nu} = \frac{i}{2}(w_{\mu}\bar{w}_{\nu} - w_{\nu}\bar{w}_{\mu}). \qquad (6.5)$$
The Nambu electromagnetic field and Z-field, as given in (2.26), become

$$ 
\begin{aligned}
e\mathcal{F}_{\mu\nu} &= g^2 B_{\mu\nu} + g'^2 W_{\mu\nu}^3 = eF_{\mu\nu} + e^2\psi_{\mu\nu}, \\
\mathcal{Z}_{\mu\nu} &= B_{\mu\nu} - W_{\mu\nu}^3 = Z_{\mu\nu} - g^2\psi_{\mu\nu}.
\end{aligned} 
\qquad (6.6) $$

Here, $F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}$ and $Z_{\mu\nu} = \partial_{\mu}Z_{\nu} - \partial_{\nu}Z_{\mu}$ are the 't Hooft fields, with the potentials

$$ eA_{\mu} = g^2 B_{\mu} + g'^2 W_{\mu}^3, \quad Z_{\mu} = B_{\mu} - W_{\mu}^3. \qquad (6.7) $$

Notice that $\psi_{\mu\nu}$ in (6.5) is the same as in (2.38), while $F_{\mu\nu}$ is the same as in (2.36).

The electroweak equations (2.17) then take the form:

$$ \nabla^{\mu} F_{\nu\mu} = 4\pi J_{\nu}, \qquad (6.8a) $$

$$ \nabla^{\mu} Z_{\nu\mu} + \frac{\phi^2}{2} Z_{\nu} = -4\pi \frac{g}{g'} J_{\nu}, \qquad (6.8b) $$

$$ \mathcal{D}^{\mu} w_{\mu\nu} + i (g^2 \psi_{\nu\sigma} + gg' F_{\nu\sigma} - g^2 Z_{\nu\sigma}) w^{\sigma} = \frac{g^2 \phi^2}{2} w_{\nu}, \qquad (6.8c) $$

$$ \nabla^{\mu} \nabla_{\mu} \phi = \frac{1}{4} (g^2 w_{\mu} \bar{w}^{\mu} + Z_{\mu} Z^{\mu}) \phi + \frac{\beta}{4} (\phi^2 - 1) \phi. \qquad (6.8d) $$

The gauge covariant derivative can be written as

$$ \mathcal{D}_{\mu} = \nabla_{\mu} + i(eA_{\mu} - g^2 Z_{\mu}), \qquad (6.9) $$

and the current is

$$ 4\pi J_{\nu} = e \nabla^{\sigma} \psi_{\sigma\nu} + e \Im(\bar{w}^{\sigma} w_{\sigma\nu}). \qquad (6.10) $$

One has $\mathcal{F}_{\mu\nu} = F_{\mu\nu} + e\psi_{\mu\nu}$, whose divergence yields

$$ 4\pi \mathcal{J}_{\nu} = \nabla^{\mu} \mathcal{F}_{\nu\mu} = \nabla^{\mu} F_{\nu\mu} + e \nabla^{\mu} \psi_{\nu\mu} = 4\pi J_{\nu} - e \nabla^{\mu} \psi_{\mu\nu} = e \Im(\bar{w}^{\sigma} w_{\sigma\nu}), \qquad (6.11) $$

hence the 't Hooft current and Nambu current are related to each other as follows,

$$ J_{\nu} = \mathcal{J}_{\nu} + \frac{e}{4\pi} \nabla^{\sigma} \psi_{\sigma\nu}. \qquad (6.12) $$

## B. Condensate energy

Eqs.(6.8) admit the solution

$$ A_{\mu} dx^{\mu} = -\frac{n}{2e} \cos\vartheta d\vartheta, \quad Z_{\mu} = w_{\mu} = 0, \quad \phi = 1, \qquad (6.13) $$

assuming the background geometry to be the RN given by (3.6),(3.7).

Linearizing the equations around this solution, the only bounded solution for the first order corrections of $A_{\mu}$, $Z_{\mu}$, $\phi$ is the trivial one; hence, these fields remain unchanged in the first perturbative order. However, the first order correction for $w_{\mu}$ can be non-trivial since it fulfils the same Proca equation as in (4.3), whose bounded solution is given by Eq.(6.1) and
(6.2) above. This solution has only two non-vanishing components, $w_ϑ$ and $w_φ = -iw_ϑ \sin \vartheta$ (assuming that $n > 0$), implying the following relations:

$$|w_μ|^2 \equiv w_μ \bar{w}^μ = \frac{2}{r^2} |\psi(r)|^2 \times \Theta(\vartheta, \varphi), \qquad (6.14) \\ \Theta(\vartheta, \varphi) = (\sin \vartheta)^{2j} \sum_{k,m} c_k c_m \left(\tan \frac{\vartheta}{2}\right)^{k+m} e^{i(k-m)\varphi}.$$

The tensor $\psi_{\mu\nu} = -\psi_{\nu\mu}$, defined in (6.5), also has only two components,

$$\psi_{\vartheta\varphi} = -\psi_{\varphi\vartheta} = -\frac{r^2}{2}|w_μ|^2 \sin\vartheta = -|\psi(r)|^2 \Theta(\vartheta, \varphi) \sin\vartheta. \qquad (6.15)$$

The nonzero components of the strength tensor $w_{\mu\nu}$ are

$$w_{r\vartheta} = -w_{\vartheta r} = \partial_r w_{\vartheta}, \quad w_{r\varphi} = -w_{\varphi r} = \partial_r w_{\varphi}, \qquad (6.16)$$

while

$$w_{\vartheta\varphi} = \partial_{\vartheta} w_{\varphi} - (\partial_{\varphi} + ieA_{\varphi}) w_{\vartheta} = w_{\varphi\vartheta} = 0. \qquad (6.17)$$

As a result,

$$(\psi_{\mu\nu})^2 = \frac{1}{2}|w_\mu|^4, \quad |w_{\mu\nu}|^2 = 2N(r) \left|\frac{\psi'(r)}{\psi(r)}\right|^2 |w_\mu|^2 \equiv 2\Omega^2(r)|w_\mu|^2. \qquad (6.18)$$

It follows that $\Im(\bar{w}^\sigma w_{\sigma\nu}) = 0$, hence the Nambu current vanishes, $J_\mu = 0$, while the 't Hooft current (6.10) becomes

$$4\pi J_\mu = e \nabla^\sigma \psi_{\sigma\mu}. \qquad (6.19)$$

At this point, it is important to emphasize the difference between the Nambu tensor $\mathcal{F}_{\mu\nu}$ and the 't Hooft tensor $F_{\mu\nu}$. The electric current $J^\mu$ associated with $\mathcal{F}_{\mu\nu}$ vanishes in the second perturbative order, but the magnetic current $\tilde{J}^\mu$ does not vanish and determines the magnetic charge density. Conversely, for the tensor $F_{\mu\nu}$ the situation is reversed:

$$\mathcal{J}^\mu = \frac{1}{4\pi} \nabla_\sigma \mathcal{F}^{\mu\sigma} = 0, \qquad \tilde{\mathcal{J}}^0 = \frac{1}{4\pi} \nabla_\sigma \tilde{\mathcal{F}}^{0\sigma} = \frac{e}{8\pi\sqrt{-g}} \epsilon^{ijk} \partial_i \psi_{jk}, \\ J^\mu = \frac{1}{4\pi} \nabla_\sigma F^{\mu\sigma} = \frac{e}{4\pi} \nabla_\sigma \psi^{\sigma\mu}, \qquad \tilde{J}^0 = \frac{1}{4\pi} \nabla_\sigma \tilde{F}^{0\sigma} = 0. \qquad (6.20)$$

Therefore, $F_{\mu\nu}$ defines the electric current, whereas $\mathcal{F}_{\mu\nu}$ gives rise to the magnetic current. However, it should be noted that the electric current $J^\mu$ starts deviating from zero at higher perturbative orders, while the magnetic current $\tilde{J}^\mu$ vanishes at all orders.

The non-zero components of $J_\mu$ are

$$J_\vartheta = \frac{e}{8\pi \sin \vartheta} \partial_\varphi |w_\mu|^2, \quad J_\varphi = -\frac{e \sin \vartheta}{8\pi} \partial_\vartheta |w_\mu|^2. \qquad (6.21)$$

This current is tangent to the horizon and orthogonal to the gradient of $|w_\mu|^2$. Restricting to $r = r_h$, the condition $|w_\mu|^2 = \text{const.}$ determines lines on the horizon that are orthogonal to the gradient of $|w_\mu|^2$. Therefore, the current $J^\mu$ is tangent to the level lines of $|w_\mu|^2$.
The current, being quadratic in $w_μ$, acts as a source for the second order correction $f_{μν} = ∂_μ f_ν - ∂_ν f_μ$ to the background electromagnetic field, leading to

$$F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu} + f_{\mu\nu} \equiv F^{(0)}_{\mu\nu} + f_{\mu\nu}, \qquad (6.22)$$

with the same $A_μ$ as in (6.13). The current also induces nonzero $Z_μ$ and a deviation $δφ = φ - 1$. Expanding Eqs. (6.8) to second order gives

$$\nabla^{\mu} f_{\nu\mu} = e \nabla^{\sigma} \psi_{\sigma\nu} \equiv 4\pi J_{\nu}^{(f)}, \qquad (6.23a)$$

$$\nabla^{\mu} Z_{\nu\mu} + \frac{1}{2} Z_{\nu} = -g^2 \nabla^{\sigma} \psi_{\sigma\nu} \equiv 4\pi J_{\nu}^{(Z)}, \qquad (6.23b)$$

$$-\nabla^{\mu} \nabla_{\mu} \delta\phi + \frac{\beta}{2} \delta\phi = -\frac{g^2}{4} |w_{\mu}|^2 \equiv 4\pi J^{(\phi)}, \qquad (6.23c)$$

whose solution is quadratic in $c_m$. Substituting this solution into (6.4) yields the condensate energy density, containing terms up to fourth order in $c_m$:

$$\mathcal{E}_{\text{cond}} \equiv -\mathcal{L}_{\text{EW}} = \frac{1}{4}(gF_{\mu\nu} + g'Z_{\mu\nu})^2 + \frac{1}{4}(g'F_{\mu\nu} - gZ_{\mu\nu} + g\psi_{\mu\nu})^2 \\ + \frac{1}{4}|w_{\mu\nu}|^2 + \frac{g^2}{4}|w_{\mu}|^2(1 + 2\delta\phi) + \frac{1}{4}(Z_{\mu})^2 + (\partial_{\mu}\delta\phi)^2 + \frac{\beta}{2}\delta\phi^2. \qquad (6.24)$$

This determines the condensate energy

$$E_{\text{cond}} = \int \mathcal{E}_{\text{cond}} \sqrt{-g} \, d^3 x. \qquad (6.25)$$

Integrating by parts, using Eqs. (6.23), and omitting the zeroth-order contribution yields

$$E_{\text{cond}} = E_w + E_{fZ\phi}, \qquad (6.26)$$

where

$$\begin{aligned} E_w &= \int \left( \frac{e}{2} F^{(0)}_{\mu\nu} \psi^{\mu\nu} + \frac{g^2}{4} (\psi_{\mu\nu})^2 + \frac{1}{4} |w_{\mu\nu}|^2 + \frac{g^2}{4} |w_{\mu}|^2 \right) \sqrt{-g} \, d^3 x \\ &= \int \left( \frac{g^2}{8} |w_{\mu}|^4 + \frac{\Omega^2(r) + m_w^2(r)}{2} |w_{\mu}|^2 \right) \sqrt{-g} \, d^3 x, \\ E_{fZ\phi} &= 4\pi \int \left( -\frac{1}{2} f^{\nu} J_{\nu}^{(f)} - \frac{1}{2} Z^{\nu} J_{\nu}^{(Z)} - \delta\phi J^{(\phi)} \right) \sqrt{-g} \, d^3 x. \end{aligned} \qquad (6.27)$$

Here $Ω^2(r)$ is defined in (6.18), while

$$m_w^2(r) = m_w^2 - |B(r)| = \frac{g^2}{2} - \frac{|n|}{2r^2}, \qquad (6.28)$$

where the last term arises from the $F_{\mu\nu}\psi^{\mu\nu}$ contribution to (6.27). The modulus $|n|$ appears because flipping the sign of $n$ changes the signs of both $F_{\mu\nu}$ and $\psi_{\mu\nu}$. This explicitly demonstrates that the Zeeman interaction with the background field effectively reduces the W-boson mass.
The $E_w$ term can be computed using the already known first order solution for $w_μ$, while $E_{fZφ}$ requires solving Eqs. (6.23) to determine the second order corrections $f_μ, Z_μ, δφ$. These can be obtained through the following steps.

First, expand the function $Θ(ϑ, φ)$ defined in (6.14) over spherical harmonics:

$$ \Theta(\vartheta, \varphi) = \sum_{J=0}^{2j} \sum_{M=-J}^{J} C_{JM} Y_{JM}(\vartheta, \varphi), \qquad (6.29) $$

with $j = |n/2| - 1$. As explained above after Eq. (6.2), the azimuthal numbers k, m in Eq. (6.14) can either be all integers or all half-integers. Therefore, their difference, k - m, is always an integer, which implies that the azimuthal number M in (6.29) must be an integer. Since $2j$ is also an integer, it follows that $J$ must be an integer as well. Thus, Eq. (6.29) involves sums over integer values of $J$ and $M$.

The next step is to introduce

$$ X^{(I)}(r, \vartheta, \varphi) = \sum_{J=0}^{2j} \sum_{M=-J}^{J} R_{JM}^{(I)}(r) Y_{JM}(\vartheta, \varphi), \qquad (6.30) $$

where the index $I$ runs over $f, Z, φ$, and the radial functions satisfy the equation

$$ \left( \frac{1}{q} \frac{d}{dr} qN \frac{d}{dr} - \frac{J(J+1)}{r^2} - m^2 \right) R_{JM}^{(I)} = a \frac{\psi^2(r)}{r^2} C_{JM}, \qquad (6.31) $$

with the following parameters:

$$ \begin{aligned} I &= f: \quad q = 1, \quad a = -e, \quad m^2 = 0; \\ I &= Z: \quad q = 1, \quad a = g^2, \quad m^2 = m_Z^2 = 1/2; \\ I &= \phi: \quad q = r^2, \quad a = g^2/2, \quad m^2 = m_H^2 = \beta/2. \end{aligned} \qquad (6.32) $$

Since $f_{μν} = ∂_μ f_ν - ∂_ν f_μ$, Eqs. (6.23) are fulfilled by

$$ f_{\vartheta} = \frac{1}{\sin \vartheta} \partial_{\varphi} X^{(f)}, \quad f_{\varphi} = -\sin \vartheta \partial_{\vartheta} X^{(f)}, \qquad (6.33) $$

$$ Z_{\vartheta} = \frac{1}{\sin \vartheta} \partial_{\varphi} X^{(Z)}, \quad Z_{\varphi} = -\sin \vartheta \partial_{\vartheta} X^{(Z)}, \quad \delta\varphi = X^{(\phi)}. $$

Injecting this solution to $E_{fZφ}$ in (6.27) gives, after some algebra,

$$ E_{fZ\phi} = \frac{1}{2} \sum_{J,M} \bar{C}_{JM} \int_{r_h}^{\infty} |\psi(r)|^2 \left\{ \frac{J(J+1)}{r^2} \left( g^2 R_{JM}^{(Z)} - eR_{JM}^{(f)} \right) + g^2 R_{JM}^{(\phi)} \right\} dr, \qquad (6.34) $$

and one should stress that the formula contains the complex conjugated $\bar{C}_{JM}$ and not $C_{JM}$. As a result, the problem reduces to solving the radial equations (6.31). It follows that

$$ f_{\vartheta\varphi} = \sin\vartheta \sum_{J,M} J(J+1) R_{JM}^{(f)}(r) Y_{JM}(\vartheta, \varphi), \qquad (6.35) $$

hence

$$ \oint \frac{1}{2} f_{ik} dx^i \wedge dx^k = \oint f_{\vartheta\varphi} d\vartheta d\varphi = 0. \qquad (6.36) $$
Therefore, the total magnetic flux of $F_{\mu\nu}$, given by (6.22), is the same as for the background:

$$ \oint \frac{1}{2} F_{ik} dx^i \wedge dx^k = \frac{2\pi n}{e}, \qquad (6.37) $$

but the magnetic field is no longer purely radial and homogeneous on the sphere, since components of $f_{ik}$ show a general dependance on $\vartheta, \varphi$, while the components $f_{r\vartheta}, f_{r\varphi}$ do not vanish.

As we shall see below, the fields $f_{\mu\nu}$ and $Z_{\mu\nu}$ create $|n| - 2$ vortices orthogonal to the horizon, forming the black hole corona.

### C. Energy minimization

The condensate energy $E_w + E_{fZ\phi}$ should be minimized with respect to the coefficients $c_m$. Since computing the $E_{fZ\phi}$ term requires additional effort, we leave this for future work and focus only on minimizing $E_w$. Already, this gives the expected result: a lattice of vortices on the horizon, forming the electroweak corona.

This suggests that the minima of $E_w$ coincide with those of $E_w + E_{fZ\phi}$ or close to them. Although we cannot immediately prove this, it seems natural that the $E_w$ term should be dominant because it contains only $w_\mu$, which is the leading correction to the background configuration, while $E_{fZ\phi}$ contains contributions from the subleading second-order corrections for the $f, Z, \phi$ fields. Therefore, we minimize

$$ \langle |w_\mu|^4 \rangle \equiv \int |w_\mu|^4 \sqrt{-g} \, d^3 x, \qquad (6.38) $$

while keeping the norm fixed:

$$ \langle |w_\mu|^2 \rangle \equiv \int |w_\mu|^2 \sqrt{-g} \, d^3 x = \text{const.} \qquad (6.39) $$

This is equivalent to minimizing

$$ E_w = \text{const.} \times \langle |w_\mu|^4 \rangle + \text{const.} \times \langle |w_\mu|^2 \rangle. \qquad (6.40) $$

Using (6.14), one finds that $\langle |w_\mu|^2 \rangle = \text{const.} \times E_2$ with

$$ E_2 = \sum_{m \in [-j, j]} A_{j,m} c_m^2, \qquad (6.41) $$

where

$$ A_{j,m} = \int_0^\pi (\sin \vartheta)^{2j+1} \left( \tan \frac{\vartheta}{2} \right)^{2m} d\vartheta = 2^{2j+1} \frac{\Gamma(j+1+m)\Gamma(j+1-m)}{\Gamma(2j+2)}. \qquad (6.42) $$

Similarly, after some algebra, one finds that $\langle |w_\mu|^4 \rangle = \text{const.} \times E_4$ with

$$ E_4 = \sum_{k,m,l} A_{2j,k+1} c_m c_k c_l c_{k+1-l}. \qquad (6.43) $$
<figure><img src="image_8.png" /><figcaption>Figure 6. Left: the horizon distribution of the W-condensate $\bar{w}^{\mu}w_{\mu}$ minimizing the energy for $n = 10$. The level lines coincide with the current flow, forming loops around vortices that start at the horizon in the radial direction (dark spots). The current flows in the opposite direction around the yellow regions, where the condensate is maximal and where the vortices return to the horizon. Right: the same distribution when all vortices merge into two oppositely directed multi-vortices.</figcaption></figure>

It is assumed here that $c_m = 0$ if $|m| > j$.

Therefore, omitting constant factors, the problem reduces to the minimization of

$$E_c = E_4 + \mu (E_2 - 1), \qquad (6.44)$$

with respect to $c_m$, where $\mu$ is the Lagrange multiplier, and with $E_2$ and $E_4$ expressed by Eqs.(6.41),(6.43).

Assuming, for example, that $n = 10$, hence $j = 4$ and $m = -4, -3, \dots, 4$, there are 9 coefficients $c_m$ to find. Their values corresponding to the absolute minimum of $E_c$ determine the function $|w_{\mu}|^2$ in (6.14), whose horizon profile is shown in Fig.6 (left panel). The figure shows the level lines, which coincide with the flow lines of the current $J^{\mu}$.

As seen in Fig.6, the flow lines form loops. A loop of current generates a magnetic field in the orthogonal direction, so the loops in Fig.6 encircle radially directed vortices located at the points where the condensate $|w_{\mu}|^2$ vanishes (dark spots). There are 4 vortices on the visible side of the surface in Fig.6 (left panel), and the same number on the other side, so altogether there are 8 vortices ($8 = |n| - 2$) homogeneously distributed over the horizon. The condensate exhibits, in this case, the same discrete symmetries as an octahedron.

The difference $|n| - 2$ arises because, for $|n| = 2$, the black hole is spherically symmetric, and the condensate produces no currents, meaning there are no vortices. We checked that minimization of (6.44) gives one vortex for $n = 3$, two vortices for $n = 4$, and so on. These vortices repel each other and form a lattice on the horizon such that their mutual separations are maximal. These solutions correspond to the global energy minimum.

The vortices start at the horizon where the condensate $|w_{\mu}|^2$ vanishes (dark spots in Fig.6), but, as we shall now see, they eventually return to the horizon in the regions where the condensate $|w_{\mu}|^2$ is maximal (yellow spots in Fig.6).
## D. Structure of vortices

The vortices carry non-zero magnetic and Z fluxes, but their net flux through the horizon is zero. Therefore, they start at the horizon and then return, being attached to the horizon at both ends, similarly to the coronal loops in the Sun's atmosphere.

The vortices leave the horizon in the dark regions, which are encircled by an electric current flowing in the anticlockwise direction (for $n > 0$). They return back in the yellow regions, where the current flows in a clockwise direction. The fluxes through the dark and yellow regions mutually compensate, as can be seen as follows.

One notices that the term containing the derivatives in the radial equation (6.31) vanishes at the horizon because the first derivative of $R_{JM}^{(f)}(r)$ vanishes there due to the boundary conditions. This means that $f_{r\vartheta} = f_{r\varphi} = 0$ at the horizon.

Eq.(6.31) then implies that for $J > 0$, one has $J(J+1)R_{JM}^{(f)}(r_h) = e\psi^2(r_h)C_{JM}$. Injecting this to (6.35) (where $J > 0$) and using (6.29), yields the horizon value of the magnetic field:

$$f_{\vartheta\varphi} = e\psi^2(r_h) \sin\vartheta \sum_{J>0,M} C_{JM} Y_{JM}(\vartheta, \varphi) = e\psi^2(r_h) \sin\vartheta \{\Theta(\vartheta, \varphi) - C_{00}Y_{00}\}, \quad (6.45)$$

where $C_{00} = \oint \Theta(\vartheta, \varphi)Y_{00} d^2\Omega$ with $d^2\Omega = \sin\vartheta d\vartheta d\varphi$.

Therefore, the condensate creates a magnetic field that is purely radial at the horizon:

$$\delta\mathcal{B}^r = \frac{\epsilon^{r\vartheta\varphi}}{\sqrt{-g}} f_{\vartheta\varphi} = \frac{n}{|n|} \times \frac{e}{2} \left( |w_\mu|^2 - \frac{1}{4\pi} \oint |w_\mu|^2 d^2\Omega \right), \quad (6.46)$$

where we have used Eq.(6.14) and restored the dependence on the sign of $n$.

The total flux of $\delta\mathcal{B}^r$ through the horizon is zero, but locally this field is non-zero. Specifically, $\delta\mathcal{B}^r > 0$ where the condensate density $|w_\mu|^2$ is large, which results in the enhancement of the background magnetic field $\mathcal{B}^r = n/(2e r_h^2)$. On the other hand, $\delta\mathcal{B}^r < 0$ near the vortex centers where $|w_\mu|^2 = 0$, leading to the weakening of the background field. This phenomenon corresponds to the anti-screening property of the condensate [47, 48]; see Appendix C.

Therefore, the sign of $\delta\mathcal{B}^r$ is determined by the value of the condensate density $|w_\mu|^2$, whose profile is shown in Fig.6. The vortices leave the horizon in regions marked by dark colours in Fig.6, where $\delta\mathcal{B}^r < 0$ (for $n > 0$) and carry away a negative magnetic flux. They return to the horizon in the yellow regions, bringing back the same flux.

The total magnetic flux carried by a vortex can be obtained by integrating $\delta\mathcal{B}^r$ over a dark region where it is sign definite. The flux of the massive Z field can be computed similarly and will have an opposite sign, because of the opposite sign of the source term for the Z field in Eq.(6.31).

The magnetic charge density $\sim \epsilon^{ijk}\partial_i\psi_{jk}$ is also maximal in the yellow regions in Fig.6, together with the condensate density $|w_\mu|^2$. Therefore, these regions are magnetically charged when considered within the Nambu approach. There are 8 vortices and 6 charged yellow regions in Fig.6 (left panel).

Each vortex center (dark region) is surrounded by three yellow regions. Since none of the latter is favored, each vortex, carrying a negative (for $n > 0$) outgoing magnetic flux, must split into three flux jets flowing towards the three neighbouring yellow regions. Each yellow region is surrounded by four vortices and, hence, absorbs incoming jets from four neighbouring vortices. Therefore, there are 8 vortices splitting into 24 jets, which together form the black hole corona.
## E. Axial symmetry

The axially symmetric configuration with $c_m \propto \delta_{m0}$ is also a stationary point of the energy (6.44). It can be viewed as a system in which all elementary vortices merge into two oppositely directed multi-vortices, generated by two oppositely directed azimuthal currents. The system is invariant under $\vartheta \to \pi - \vartheta$, except for the current density, which changes sign.

The condensate $|w_\mu|^2$, the radial magnetic field $F_{\vartheta\varphi}$, and the magnetic charge density $\sim \epsilon^{ijk}\partial_i\psi_{jk}$ attain maximal values at the equator, where the azimuthal current changes sign, as shown in the right panel of Fig.6 (see also (6.21)).

However, the Hessian matrix has negative eigenvalues in this case, indicating that such solutions are unstable and should decay into non-axially symmetric configurations. This is because the elementary vortices forming a multi-vortex repel each other.

Setting $c_m = c_0\delta_{m0}$ and minimizing the energy yields $c_0 = 1/\sqrt{2}$. Eq.(6.14) then implies that

$$ \Theta(\vartheta, \varphi) = \frac{1}{2} (\sin \vartheta)^{2j}. \qquad (6.47) $$

Inserting this Eqs.(6.45),(6.46), we obtain the horizon value of the magnetic field (compare with Eq. (C.15) in Appendix C):

$$ \delta\mathcal{B}^r = e \frac{n}{|n|} \frac{\psi^2(r_h)}{2r_h^2} \times \left( (\sin\vartheta)^{2j} - \frac{\sqrt{\pi}\Gamma(j+1)}{2\Gamma(j+3/2)} \right), \qquad (6.48) $$

where $j = |n|/2 - 1 \ge 0$. This describes two multi-vortices made of $|n| - 2$ elementary vortices: one along the z-axis in the northern hemisphere (the dark spot in Fig.6, right panel), and the other one along the negative direction of the z-axis in the southern hemisphere. Each vortex starts at the horizon in the polar region, extends outside, spreads out like a fountain, and then returns back to the horizon in the equatorial region.

## F. Special cases

Let us consider the lowest values of $n$.

1. $|n| = 2$

In the case where the condensate is spherically symmetric, the currents vanish, and consequently, there are no vortices. This leads to $\delta\mathcal{B}^r = 0$ (see (6.48) with $j = 0$), meaning the magnetic field is purely radial and homogeneous on the sphere.

The energy $E_c$ in (6.44) depends only on one coefficient $c_0$, and the minimization yields a unique solution: $c_0 = 1/\sqrt{2}$, corresponding to the global minimum.

Thus, it is plausible that non-perturbative generalizations of this solution – the hairy black holes with $|n| = 2$, are stable, as discussed in Section V D.

2. $|n| = 4$

If we do not assume the axial symmetry, the condensate corresponding to the global energy minimum would be determined by three coefficients $c_0, c_{\pm 1}$. However, having obtained
these coefficients by minimizing the energy, they can be transformed to $c_m \propto \delta_{m0}$ by a global rotation. This means that the global minimum remains axially symmetric.

In this case, there are $|n| - 2 = 2$ elementary vortices, which position themselves at opposite sides of the sphere due to their mutual repulsion. Since this perturbative solution is stable, it is plausible that its non-perturbative generalizations – the hairy black holes with $|n| = 4$, would also be stable. These solutions will be explored in Section VIII.

3. $|n| = 3$

In this case there are two axially symmetric solutions. The first solution contains a single vortex along the positive z-axis, say. This solution it is not invariant under $\vartheta \to \pi - \vartheta$, since it has a preferred direction for the vortex. Consequently, it must support a nonzero dipole moment. Despite this, it represents the global minimum of the energy.

The second solution is invariant under $\vartheta \to \pi - \vartheta$, and it contains two oppositely directed half-vortices. While this solution is symmetric, it is energetically less favorable than the first one and does not correspond to the global energy minimum.

These two solutions illustrate the balance between symmetry and energy minimization. The first solution is less symmetric but stable, while the second solution is more symmetric but less stable.

### G. Conclusion of the perturbative analysis

Summarizing the above discussion, minimum-energy hairy solutions with higher magnetic charge $|n|$ can have at most only discrete symmetries. As a result, constructing them at the non-perturbative level would require full 3D numerical simulations. However, in the special case of axial symmetry, the magnetic charge can take arbitrary values, but the complexity of the problem is reduced. In this case, only 2D simulations are needed, although they do not necessarily yield the global energy minimum.

Therefore, in what follows, we shall restrict to the axially symmetric sector.

## VII. AXIAL SYMMETRY

Choosing the spacetime coordinates as $x^\mu = (t, r, \vartheta, \varphi)$, we assume the invariance under the action of two symmetry generators:

$$ \iota = \frac{\partial}{\partial t}, \qquad \varsigma = \frac{\partial}{\partial \varphi}, \qquad (7.1) $$

where the timelike Killing vector $\iota$ is assumed to be hypersurface-orthogonal. Additionally, we impose invariance under reflections across the equatorial plane, $\vartheta \to \pi - \vartheta$.

### A. Static and axially symmetric gravitational fields

The line element is parameterized as

$$ ds^2 = -e^{2U}N(r)dt^2 + e^{-2U}dl^2, \qquad dl^2 = e^{2K}\left[\frac{dr^2}{N(r)} + r^2d\vartheta^2\right] + e^{2S}r^2\sin^2\vartheta d\varphi^2, \qquad (7.2) $$
where U, K, S depend on r, ϑ and N(r) is an auxiliary function. It is important to note that the radial coordinate r is not the same as the Schwarzschild coordinate r used in previous sections. The radial coordinate is defined up to diffeomorphisms $r \to \tilde{r} = \tilde{r}(r)$, which lead to the transformations $U \to \tilde{U}$, $N \to \tilde{N}$, $K \to \tilde{K}$, $S \to \tilde{S}$, where

$$ \frac{dr}{r\sqrt{N}} = \frac{d\tilde{r}}{\tilde{r}\sqrt{\tilde{N}}}, \qquad e^U \sqrt{N} = e^{\tilde{U}} \sqrt{\tilde{N}}, \qquad (7.3) $$

$$ r\sqrt{N}e^K = \tilde{r}\sqrt{\tilde{N}}e^{\tilde{K}}, \qquad r\sqrt{N}e^S = \tilde{r}\sqrt{\tilde{N}}e^{\tilde{S}}. \qquad (7.4) $$

This gauge freedom can be fixed by choosing a specific form for the auxiliary function $N(r)$. For example, one can set $N(r) = 1$.

However, even after fixing the radial coordinate r, the line element still admits a large gauge symmetry. Indeed, if $N(r) = 1$, then the r, ϑ part of the metric can be represented as

$$ e^{2K-2U}(dr^2 + r^2 d\vartheta^2) = e^{2K-2U}|d\varpi|^2, \qquad (7.5) $$

where the complex variable is $\varpi = r e^{i\vartheta}$. Rewriting this in terms of a new complex variable $\tilde{\varpi}$ such that $\varpi = \varpi(\tilde{\varpi})$, the line element retains the same structure, up to the replacement

$$ K \to \tilde{K} = K + \ln \left| \frac{d\varpi}{d\tilde{\varpi}} \right|. \qquad (7.6) $$

This reveals that the r, ϑ sector of the geometry enjoys a conformal symmetry generated by holomorphic functions $\varpi(\tilde{\varpi})$. However, it turns out that this symmetry is ultimately fixed by the boundary conditions given in Eq. (7.7) below, which permit only the identity transformation $\varpi = \tilde{\varpi}$.

For non-extremal black hole solutions we choose $N(r) = 1 - r_H/r$, where $r = r_H$ corresponds to the position of the event horizon. To ensure that the geometry is asymptotically flat and regular outside the event horizon, we impose the following boundary conditions:

$$ \text{infinity, } \underline{r = \infty} : U = K = S = 0; $$

$$ \text{horizon, } \underline{r = r_H} : |\partial_r U| < \infty, |\partial_r K| < \infty, |\partial_r S| < \infty; $$

$$ \text{symmetry axis, } \underline{\vartheta = 0, \pi} : \partial_{\vartheta} U = \partial_{\vartheta} K = \partial_{\vartheta} S = 0; $$

$$ \text{equator, } \underline{\vartheta = \pi/2} : \partial_{\vartheta} U = \partial_{\vartheta} K = \partial_{\vartheta} S = 0. \qquad (7.7) $$

It turns out that imposing these conditions automatically ensures that at the symmetry axis one has

$$ \underline{\vartheta = 0, \pi} : K = S, \qquad (7.8) $$

which is needed for the absence of the conical singularity.

The functions U, K, S satisfy second-order differential equations, which are obtained by taking suitable linear combinations of the Einstein equations (2.19):

$$ \left( \frac{1}{2} e^{S-2U} r^2 (Ne^{2U})'_{r} \sin \vartheta \right)'_{r} + (e^S U'_{\vartheta} \sin \vartheta)'_{\vartheta} = -\kappa (T^0_0 - \frac{1}{2} T) \sqrt{-g}, \qquad (7.9a) $$

$$ \left( r\sqrt{N} \left( r\sqrt{N}e^S \right)'_{r} \sin \vartheta \right)'_{r} + (e^S \sin \vartheta)''_{\vartheta\vartheta} = \kappa (T^r_r + T^\vartheta_\vartheta) \sqrt{-g}, \qquad (7.9b) $$

$$ r\sqrt{N} \left( r\sqrt{N} K'_{r} \right)'_{r} + K''_{\vartheta\vartheta} + Nr^2 (U'_{r})^2 + (U'_{\vartheta})^2 + r^2 N' U'_{r} + \frac{1}{2} (r^2 N'_{r})'_{r} = r^2 e^{2K-2U} \kappa T^\varphi_\varphi, \qquad (7.9c) $$
with $T = T^\mu_\mu$. At large distances, where the theory reduces to Einstein-Maxwell, one has $T_r^r + T_\vartheta^\vartheta = 0$, which allows for an exact integration of the S-equation (7.9b). The solution which remains bounded for any $\vartheta$ is given by

$$r\sqrt{N}e^S = A_1e^\xi + A_2e^{-\xi} \quad \text{with} \quad \xi = \int \frac{dr}{r\sqrt{N}}. \qquad (7.10)$$

The integration constants $A_1, A_2$ can be adjusted to ensure that $S \to 0$ as $r \to \infty$. At shorter distances, where $T_r^r + T_\vartheta^\vartheta \neq 0$, the solution can no longer be obtained analytically. Instead, the S-equation must be solved numerically together with the other field equations in (7.9). Among the remaining Einstein equations (2.19), there are also constraints

$$C_1 \equiv r\sqrt{-g}(G^r_r - \kappa T^r_r) = 0, \quad C_2 \equiv \frac{\sqrt{-g}}{r}(G^r_\vartheta - \kappa T^r_\vartheta) = 0, \qquad (7.11)$$

which do not involve second derivatives of U, K, but only those of S. These equations can be solved for $\partial_r K$ and $\partial_\vartheta K$, giving

$$\partial_r K = K_1, \quad \partial_\vartheta K = K_2, \qquad (7.12)$$

where $K_1, K_2$ are functions that do not depend on K. The function K can be obtained by integrating these two equations, as in the vacuum case [71], provided that the compatibility condition,

$$\partial_\vartheta K_1 = \partial_r K_2, \qquad (7.13)$$

is satisfied. This condition follows from equations (7.9a) and (7.9b), while equation (7.9c) is automatically fulfilled as a consequence of (7.12).

In practice, however, it is preferable to treat the two equations (7.11) as constraints rather than solving them explicitly. Specifically, the two non-trivial Bianchi identities $\nabla_\mu(G^\mu_\nu - \kappa T^\mu_\nu) = 0$ for $\nu = r, \vartheta$ can be rewritten as

$$\partial_\vartheta C_2 = -N \partial_r C_1 - \frac{1}{2} N' C_1 + \dots, \quad \partial_\vartheta C_1 = r^2 \partial_r C_2 + r C_2 + \dots, \qquad (7.14)$$

where the dots denote terms that vanish if the second-order equations (7.9) are satisfied. It follows that if the constraints $C_1 = 0, C_2 = 0$ hold at the symmetry axis, then they will be satisfied everywhere. At the symmetry axis, we have

$$\vartheta = 0, \pi : C_1 = -r e^S \partial_\vartheta (K - 2S), \quad C_2 = r N e^S \partial_r (K - S), \qquad (7.15)$$

which vanish owing to the boundary conditions (7.7). Thus, by solving the three second order equations (7.9) with the boundary conditions (7.7), the two constraint equations (7.11) are automatically satisfied.

## B. The integral formula for the ADM mass

The ADM mass $M$ in (2.23) can be expressed as a surface integral over a sphere at spatial infinity, following Komar [72]:

$$M = -\frac{1}{8\pi} \oint_{S_\infty} \star d(\iota_\mu dx^\mu). \qquad (7.16)$$
Here $\iota_\mu$ are the covariant components of the timelike Killing vector, and the star denotes the 4-dimensional Hodge dual. Using Gauss's theorem, the surface integral at infinity can be transformed into a volume integral over the 3-space, plus a surface integral over the event horizon. This gives

$$M = \frac{k_H A_H}{4\pi} + \frac{1}{4\pi} \int_{r>r_H} R_{\hat{0}\hat{0}} \sqrt{-g} d^3 x. \qquad (7.17)$$

The first term represents the contribution from the event horizon, where $k_H$ is the surface gravity and $A_H$ is the horizon area. In the second term one has $R_{\hat{0}\hat{0}} = R_{\mu\nu}\iota^\mu\iota^\nu/|\iota^\alpha\iota_\alpha| = -R_0^0$. The surface gravity is expressed by

$$k_H^2 = -\frac{1}{2} \nabla^\mu \iota_\nu \nabla_\mu \iota^\nu \Big|_{r=r_H} . \qquad (7.18)$$

Using the line element (7.2), one obtains:

$$k_H = \frac{N'}{2} e^{2U-K} \Big|_{r=r_H}, \qquad A_H = 2\pi r_H^2 \int_0^\pi e^{K-2U+S} \sin\vartheta \, d\vartheta \Big|_{r=r_H}. \qquad (7.19)$$

The surface gravity should be constant at the horizon, which follows from the second constraint equation in (7.12). At $r = r_H$, where $N = 0$, this equation simplifies to

$$\partial_\vartheta K|_{r=r_H} = 2\partial_\vartheta U|_{r=r_H} \quad \Rightarrow \quad \partial_\vartheta (2U - K)|_{r=r_H} = 0. \qquad (7.20)$$

Thus, $2U - K$ is independent of $\vartheta$ at the horizon, ensuring that $k_H$ remains constant.

The mass formula (7.17) can be easily verified since $R_{\hat{0}\hat{0}}\sqrt{-g}$ corresponds to the right-hand side of Eq. (7.9a). Therefore, replacing it by the left-hand side of the same equation, integrating, and using the boundary condition (7.7) yields

$$\frac{1}{4\pi} \int_{r>r_H} R_{\hat{0}\hat{0}} \sqrt{-g} d^3 x = \frac{1}{4\pi} \int_0^{2\pi} d\varphi \int_0^\pi \frac{1}{2} e^{S-2U} r^2 (Ne^{2U})'_r \sin\vartheta \, d\vartheta \Big|_{r=r_H}^{r=\infty} \\ = M - \frac{1}{4\pi} \frac{N'(r_H)}{2} r_H^2 \int_0^{2\pi} d\varphi \int_0^\pi e^S \sin\vartheta \, d\vartheta \Big|_{r=r_H} = M - \frac{k_H A_H}{4\pi}. \qquad (7.21)$$

Here we used the fact that at large $r$, one has $S, U \to 0$ and $Ne^{2U} = 1 - 2M/r + O(1/r^2)$, along with the condition that $2U - K$ remains constant at the horizon.

### C. Axially symmetric electroweak fields

Since the two Killing vectors in Eq. (7.1) commute with each other, there is a gauge where the fields do not depend on $t, \varphi$. Additionally, we require the electroweak fields to be purely magnetic and their energy density to be invariant under reflections across the equatorial plane, $\vartheta \to \pi - \vartheta$. We choose the fields as follows:

$$W = T_2(F_1 dr + F_2 d\vartheta) + \nu (T_3 F_3 - T_1 F_4) d\varphi, \quad B = \nu Y d\varphi, \quad \Phi = \begin{pmatrix} \phi_1 \\ \phi_2 \end{pmatrix}, \qquad (7.22)$$
where $F_1, F_2, F_3, F_4, Y, \phi_1, \phi_2$ are 7 real-valued functions of r, $\vartheta$, and the constant is given by

$$ \nu = -n/2, \quad n \in \mathbb{Z}. \qquad (7.23) $$

This is the same condition as in (3.3), and it will be justified shortly.

The SU(2) gauge field in (7.22) corresponds to the static and purely magnetic ansatz of Rebbi and Rossi [73]. The fields (7.22) retain their form under gauge transformations (2.15) generated by $\mathcal{U} = \exp\{i\chi(r, \vartheta)T_2\}$, whose effect is

$$ \begin{aligned} F_1 &\rightarrow F_1 + \partial_r \chi, \quad F_2 \rightarrow F_2 + \partial_\vartheta \chi, \quad Y \rightarrow Y, \\ F_3 &\rightarrow F_3 \cos \chi - F_4 \sin \chi, \quad F_4 \rightarrow F_4 \cos \chi + F_3 \sin \chi, \\ \phi_1 &\rightarrow \phi_1 \cos(\chi/2) + \phi_2 \sin(\chi/2), \quad \phi_2 \rightarrow \phi_2 \cos(\chi/2) - \phi_1 \sin(\chi/2). \end{aligned} \qquad (7.24) $$

Injecting (7.2),(7.22) into (2.20) yields components of the energy-momentum tensor $T_{\mu\nu}$, shown in Appendix A. Injecting into (2.26),(2.39) yields $\mathcal{J}^\varphi$ and $\mathcal{J}^\varphi$ as the only non-zero components of the electric currents. Hence, the currents are purely azimuthal. Injecting into (2.41) yields only one non-trivial component of the magnetic current, corresponding to the magnetic charge density $\mathcal{J}^0$. The energy-momentum tensor and the currents are gauge-invariant.

Modulo gauge transformations (7.24), the zero energy configuration is

$$ F_1 = F_2 = F_4 = \phi_1 = 0, \quad \phi_2 = 1, \quad F_3 = Y = \text{const.} \equiv Y_\infty, \qquad (7.25) $$

which retains its form under gauge transformations generated by $\mathcal{U} = \exp\{iC\nu\varphi(1 + \tau_3)/2\}$ with a constant $C$, whose effect is $Y_\infty \rightarrow Y_\infty + C$.

The electroweak Lagrangian is

$$ L_{\text{WS}} = \int \mathcal{L}_{\text{WS}} \sqrt{-g} \, d^3 x = \int T_0^0 \sqrt{-g} \, d^3 x = -2\pi \int (a_1 + a_2 + a_3 + a_4 + a_5) \, dr \, d\vartheta, \qquad (7.26) $$

where $a_1, \dots, a_5$ are given by (A.2) in Appendix A. Varying the Lagrangian with respect to the seven amplitudes $F_1, F_2, F_3, F_4, Y, \phi_1, \phi_2$ yields the electroweak equations. Together with the Einstein equations (7.9), this gives a closed system of 10 second-order partial differential equations, whose solutions describe static, axially symmetric, and purely magnetic EWS fields. Before solving these equations, some preliminary work remains to be done.

## D. Removing the string singularity

The gauge fields in (7.22) contain the Dirac string singularity at the symmetry axis because their $\varphi$-components do not vanish there. This singularity can be removed by setting

$$ \begin{aligned} F_1 &= -\frac{H_1(r, \vartheta)}{r\sqrt{N}}, \quad F_2 = H_2(r, \vartheta), \quad F_3 = \cos\vartheta + H_3(r, \vartheta)\sin\vartheta, \\ F_4 &= H_4(r, \vartheta)\sin\vartheta, \quad Y = \cos\vartheta + y(r, \vartheta)\sin\vartheta, \end{aligned} \qquad (7.27) $$

and then performing the gauge transformation generated by

$$ \mathcal{U}_{\pm} = e^{-i\nu\varphi T_3} e^{-i\vartheta T_2} e^{\pm i\nu\varphi/2}, \qquad (7.28) $$
which brings the SU(2) field to the form

$$W = T_{\varphi} \left( -\frac{H_1}{r\sqrt{N}} dr + (H_2 - 1) d\vartheta \right) + \nu \left( T_r H_3 + T_{\vartheta} (1 - H_4) \right) \sin\vartheta d\varphi, \quad (7.29)$$

where the angle-dependent SU(2) generators are $T_r = e_r^a T_a$, $T_{\vartheta} = e_{\vartheta}^a T_a$, $T_{\varphi} = e_{\varphi}^a T_a$ with

$$e_r^a = (\sin\vartheta \cos(\nu\varphi), \sin\vartheta \sin(\nu\varphi), \cos\vartheta), \quad e_{\vartheta}^a = \partial_{\vartheta} e_r^a, \quad e_{\varphi}^a = \frac{1}{\nu \sin\vartheta} \partial_{\varphi} e_r^a. \quad (7.30)$$

This field is $\varphi$-dependent but regular at the symmetry axis, owing to the boundary conditions defined in Eq. (7.35) below (see [59] for details). The U(1) and Higgs fields become

$$B_{\pm} = \nu (\cos \vartheta \pm 1 + y \sin \vartheta) d\varphi, \quad (7.31)$$

$$\Phi_{\pm} = e^{\pm i \nu \varphi / 2} \begin{pmatrix} (\phi_1 \cos(\vartheta/2) - \phi_2 \sin(\vartheta/2)) e^{-i \nu \varphi / 2} \\ (\phi_1 \sin(\vartheta/2) + \phi_2 \cos(\vartheta/2)) e^{+i \nu \varphi / 2} \end{pmatrix},$$

where the upper “+” and lower “-” signs correspond to two local gauges used, respectively, in the southern and northern hemispheres. These two gauges are related to each other in the equatorial region by the U(1) transformation with $\mathcal{U} = \exp(i\nu\varphi)$. When expressed in these two local gauges, the U(1) and Higgs fields are everywhere regular.

The above formulas suggest that $\nu = -n/2$ should be integer. However, it can be half-integer as well. In that case, the $\varphi$-dependent part of $T_r$ changes sign under $\varphi \to \varphi + 2\pi$, but the correct sign can be restored by the global gauge transformation generated by $\mathcal{U} = \tau_3$ or $\mathcal{U} = -\tau_3$, whose effect is

$$T_1 \to -T_1, \quad T_2 \to -T_2, \quad (7.32)$$

so that the field (7.29) does not change. Similarly, if $\nu$ is half-integer, then the lower component of $\Phi_+$ and the upper component of $\Phi_-$ change sign under $\varphi \to \varphi + 2\pi$, but the correct sign can be restored by the gauge transformations generated by $\mathcal{U} = \tau_3$ and $\mathcal{U} = -\tau_3$, respectively. The fields (7.31) then remain invariant. As a result, all integer values of $n$ in (7.23) are allowed.

To summarize, by passing to the parameterization in terms of the seven amplitudes $H_1, H_2, H_3, H_4, y, \phi_1$, and $\phi_2$, subject to the boundary conditions to be specified below, the electroweak fields become regular at the symmetry axis.

## E. Fixing the gauge

The equations for $H_1, H_2, H_3, H_4, y, \phi_1$ and $\phi_2$ admit pure gauge solutions due to the residual gauge invariance (7.24). Such zero modes should be removed by fixing the gauge, as otherwise the differential operators in the equations are not invertible. The gauge can be fixed by setting to zero the covariant divergence of the two-vector $F_1 dr + F_2 d\vartheta$ in (7.22), computed with respect to the $r, \vartheta$ part of the metric. This requires that

$$r\sqrt{N} \partial_r H_1 = \partial_\vartheta H_2 \quad \Rightarrow \quad \partial_\xi H_1 = \partial_\vartheta H_2, \quad \text{where} \quad \xi = \int_{r_H}^r \frac{dr}{r\sqrt{N}}. \quad (7.33)$$
This gauge condition insures good numerical convergence, although it gives rise to a spurious long-range mode in the solutions at large r [59]. The condition (7.33) still allows for residual gauge transformations (7.24) with the gauge parameter $\chi$ subject to

$$ (\partial_{\xi\xi}^2 + \partial_{\vartheta\vartheta}^2) \chi = 0. \qquad (7.34) $$

Since $\xi \in [0, \infty)$, $\vartheta \in [0, \pi]$, this equation possesses bounded solutions which generate gauge transformations $H_1 \to H_1 - \partial_\xi \chi$ and $H_2 \to H_2 + \partial_\vartheta \chi$. At the same time, as seen in (7.27), $H_1$ must vanish at the horizon where $\xi = 0$ since otherwise the gauge field amplitude $F_1$ diverges. As a result, the only allowed solution of (7.34) is $\chi = \text{const.}$, but this is also ruled out since otherwise the value of the Higgs field components would change according to (7.24), whereas we impose at infinity $\phi_1 = 0$ and $\phi_2 = 1$.

Therefore, the gauge is completely fixed by the condition (7.33). Using this condition, all equations assume a manifestly elliptic form.

## F. Boundary conditions and numerical procedure

The axially symmetric EWS fields are the line element in (7.2) and the electroweak fields in (7.29),(7.31) expressed in terms of 10 functions U, K, S, $H_1$, $H_2$, $H_3$, $H_4$, y, $\phi_1$, $\phi_2$ which satisfy 10 elliptic equations. Solutions of these equations also satisfy the gravitational constraints (7.11) owing to the boundary conditions described by Eq.(7.35) below.

We use the compact radial variable $x \in [0, 1]$ such that

$$ r = \sqrt{r_H^2 + x^2/(1-x)^2}, $$

and since we assume reflection invariance of the energy density under $\vartheta \to \pi - \vartheta$, we can restrict the range of the angular variable to $\vartheta \in [0, \pi/2]$. This assumption implies that U, K, S, $H_2$, $H_4$, and $\phi_2$ remain invariant under $\vartheta \to \pi - \vartheta$, while $H_1$, $H_3$, y, and $\phi_1$ change sign (see Ref.[59] for details). Regularity at the horizon and along the symmetry axis, together with asymptotic flatness, then imposes the following boundary conditions:

$$ \begin{array}{ll} \underline{\text{axis, } \vartheta = 0:} & \partial_{\vartheta} = 0 \text{ for U, S, } H_2, H_4, \phi_2; \\ & K - S = H_1 = H_3 = y = \phi_1 = 0; \end{array} $$

$$ \begin{array}{ll} \underline{\text{equator, } \vartheta = \pi/2:} & \partial_{\vartheta} = 0 \text{ for U, K, S, } H_2, H_4, \phi_2; \\ & H_1 = H_3 = y = \phi_1 = 0; \end{array} $$

$$ \begin{array}{ll} \underline{\text{horizon, } x = 0:} & H_1 = 0 : \partial_x = 0 \text{ for the rest} \\ \underline{\text{infinity, } x = 1:} & \phi_2 = 1; 0 \text{ for the rest.} \end{array} \qquad (7.35) $$

It turns out that these boundary conditions automatically ensure that at the symmetry axis one has

$$ \underline{\vartheta = 0:} \quad H_2 = H_4, \quad \partial_{\vartheta} K = 0, \qquad (7.36) $$

as needed for the regularity.

We solve the field equations with these boundary conditions to determine the components of the “state vector”

$$ \Psi = [H_1, H_2, H_3, H_4, y, \phi_1, \phi_2, U, K, S], \qquad (7.37) $$
which are functions of r, $\vartheta$. We solve the equations with the FreeFem numerical solver based on the finite element method [74]. This solver uses the weak form of differential equations obtained by transforming them into integral equations, and the functions to be determined are expanded over elementary basis functions obtained by triangulating the integration domain. To handle the non-linearities, we use the Newton-Raphson algorithm. To estimate the numerical accuracy of the solutions, we use the virial relation that should hold on-shell.

## G. Virial relations

Consider the 1-form generating a spatial dilatation,

$$ \zeta_{\mu} dx^{\mu} = r e^{K-U} dr. \qquad (7.38) $$

The tensor $\nabla_{\mu}\zeta_{\nu} + \nabla_{\nu}\zeta_{\mu}$ is diagonal and describes the variation of the spacetime metric induced by the dilatation. The virial relation can be expressed in terms of $\zeta^{\mu} = g^{\mu\nu}\zeta_{\nu}$:

$$ v_m \equiv \int_{r>r_H} T^\mu_\nu \nabla_\mu \zeta^\nu \sqrt{-g} \, d^3x = 0. \qquad (7.39) $$

Since $T^\mu_\nu \nabla_\mu \zeta^\nu = \nabla_\mu(T^\mu_\nu \zeta^\nu) - (\nabla_\mu T^\mu_\nu)\zeta^\nu$ and the energy-momentum tensor is conserved, the volume integral can be transformed to surface integrals over a two-sphere at infinity and over the horizon,

$$ v_m = \oint_{r \to \infty} T^\mu_\nu \zeta^\nu d^2S_\mu - \oint_{r=r_H} T^\mu_\nu \zeta^\nu d^2S_\mu. \qquad (7.40) $$

As $T^\mu_\nu$ is bounded for $r \ge r_H$ and approaches zero at infinity, the surface integral at infinity vanishes, while the integral over the horizon is zero because $\zeta^\nu \sim N$. As a result, the integral in (7.39) vanishes, although its integrand does not.

It is instructive to see what happens in the spherically symmetric case when the metric reduces to (3.1), while the energy-momentum tensor is diagonal, so that we have $T^\mu_\nu = \text{diag}[T_0^0(r), T_r^r(r), T_\vartheta^\vartheta(r), T_\varphi^\varphi(r)]$. This can correspond to any matter type, not only to the electroweak fields. The virial relation (7.39) reduces to:

$$ v_m = 4\pi \int_{r_H}^{\infty} \sigma r^2 \left( N T_k^k + \frac{r N'}{2} T_r^r + \frac{r}{2\sigma^2} (\sigma^2 N)' T_0^0 \right) dr = 0, \qquad (7.41) $$

where $T_k^k = T_r^r + T_\vartheta^\vartheta + T_\varphi^\varphi$ is the sum of the principal pressures. Taking into account the conservation condition $\nabla_\mu T^\mu_\nu = 0$, it is not difficult to see that this is equivalent to:

$$ v_m = 4\pi \int_{r_H}^{\infty} (r^3 \sigma N T_r^r)' dr = 0, \qquad (7.42) $$

which indeed holds, since at the lower limit $N$ vanishes, while for $r \to \infty$, the components $T^\mu_\nu$ decay faster that $1/r^3$ (otherwise the energy would diverge).

In flat space, where $N = \sigma = 1$, (7.41) reduces to:

$$ v_m = \int T_k^k \sqrt{-g} \, d^3x = 0, \qquad (7.43) $$
which should be fulfilled for any finite energy soliton configuration [75]. This is equivalent to the classical scaling relation of Derrick [76]. Therefore, (7.39) provides the generalization of Derrick’s argument for black holes.

This generalisation is not unique, as $\zeta_\mu$ can be multiplied by a scalar function without changing the value of the virial integral. Other forms of virial relation can be found in the literature (see, e.g. [77]).

A similar virial relation in the gravity sector can be obtained by replacing $T^\mu_\nu$ with the Einstein tensor $G^\mu_\nu$:

$$v_g \equiv \int_{r>r_H} G^\mu_\nu \nabla_\mu \zeta^\nu \sqrt{-g} \, d^3x = 0. \qquad (7.44)$$

We adopt (7.39) and (7.44) as virial relations to test the quality of our numerical solutions. These relations are satisfied for all our solutions, with a precision depending on the numbers of discretization points $N_x$ and $N_\vartheta$ along the x, $\vartheta$ axes (these numbers determine the triangulation pattern for the FreeFem solver). Taking $N_x = 130$ and $N_\vartheta = 30$ typically yields: $v_m \sim 10^{-7}$, $v_g \sim 10^{-8}$.

## VIII. HAIRY SOLUTIONS WITH $n = \pm 4$

Before considering axially symmetric solutions with $|n| > 2$, we have reproduced the spherically symmetric solutions with $|n| = 2$ within the axially symmetric setting. This analysis, described in Appendix B, is needed to clarify the choice of the radial gauge. Here are the essential points.

### A. Choice of the radial gauge

The axially symmetric fields in Eqs. (7.2), (7.22), and (7.27) depend on the radial coordinate $r$. This is not the same as the Schwarzschild coordinate $r$, and it is defined up to reparametrizations. This gauge freedom is fixed by specifying the auxiliary function $N(r)$ in the line element (7.2).

For the non-extremal solutions, we choose

$$N(r) = 1 - \frac{r_H}{r}, \qquad r_H \in [0, r_H^0], \qquad (8.1)$$

where $r = r_H$ corresponds to the position of the event horizon, and we refer to $r_H$ as the “black hole size”. However, since $r$ is not the Schwarzschild coordinate, the parameter $r_H$ does not have a direct physical meaning, although its value determines the horizon area. When $r_H \to 0$, the solutions approach the extremal limit with a non-zero horizon area. When $r_H \to r_H^0$, the solutions merge with the RN solution, and we have $r_H^0 = r_+^0 - r_-^0$, where $r_+^0 = r_h^0$ and $r_-^0$ are the two RN horizons (see Eq. (B.20)).

The RN solution (3.6), when transformed to the axially symmetric form (7.2), (8.1), is expressed by Eq. (B.3) in Appendix B.

The horizon values of the metric functions U, K, S grow without bounds in the extremal limit $r_H \to 0$. Therefore, to study the extremal limit, we use a different gauge in which everything is finite at the horizon. This gauge is defined by:

$$N(r) = k^2(r) \times \left(1 - \frac{r_H}{r}\right)^2, \qquad k^2(r) = 1 - \frac{\Lambda}{3} \left[ r^2 + 2 r_H r + 3 r_H^2 \right] \times \frac{1 + r_H^4}{1 + r^4}, \qquad (8.2)$$
<figure><img src="image_9.png" /><figcaption>Figure 7. The SU(2) amplitudes $H_1, H_2, H_3, H_4$ for the hairy black hole solution with $\kappa = 10^{-3}$, $r_H = 0.7$, and $n = 4$ as functions of $\ln(r/r_H)$ and $\vartheta$.</figcaption></figure>

with $r_H = r_{ex}$, which is the horizon radius value for the extremal RNdS solution in Eq.(3.23). The horizon is located at $r = r_H$, and it is degenerate, meaning the surface gravity vanishes. This gauge choice is justified by the fact that, for some extremal hairy solutions (phase I), the horizon geometry coincides with that of the extremal RNdS geometry (3.7),(3.25). The radial coordinate $r$ then matches the Schwarzschild coordinate $r$ at the horizon, where the function $N(r)$ coincides with $N(r)$ from Eq.(3.25).

The gauge choices (8.1) and (8.2) are justified in Appendix B for spherically symmetric hairy solutions with $|n| = 2$, but they also apply to axially symmetric solutions with $|n| > 2$.

## B. Structure of hairy solutions for $n = \pm 4$

The spherically-symmetric solutions with $n = 2$ can be used as an initial guess to construct axially symmetric solutions with $n > 2$. Of course, $n$ should be integer to avoid the line singularities in the gauge fields, but the latter are not seen by the equations, which can
<figure><img src="image_10.png" /><figcaption>Figure 9. The metric amplitude S and the difference K − S for the hairy black hole solution with $\kappa = 10^{-3}$, $r_H = 0.7$, $n = 4$.</figcaption></figure>

the geometry is close to the RN with magnetic charge $P = 2/e$. We therefore choose an intermediate value, $r_H = 0.7$, and show in Figs.7–9 the profiles of the 10 amplitudes $H_k, y, \phi_1, \phi_2, U, K, S$ as functions of $\ln(r/r_H)$ and $\vartheta$. The symmetry axis is at $\vartheta = 0$, and the equator is at $\vartheta = \pi/2$.

The functions $H_2, H_4$, and $\phi_2$, which do not vanish in the spherically symmetric case $n = 2$ (see Eq.(B.2)), remain essentially the same and exhibit only weak $\vartheta$-dependence for $n = 4$. On the other hand, the functions $H_1, H_3, y, \phi_1$, which vanish for $n = 2$, are no longer zero for $n > 2$. They exhibit strong $\vartheta$-dependence, with a pronounced maximum or minimum between the symmetry axis and the equatorial plane, and are antisymmetric under $\vartheta \to \pi - \vartheta$. Their $\vartheta$-dependence remains significant even at the horizon, except for $H_1$, which vanishes there.

The metric function U, K, S almost do not depend on the angle $\vartheta$ and are close to the values defined by Eq.(B.3) in the RN case. In particular, one has $K \approx S$. The difference $K-S$ is maximal at the equator and can be seen as a measure of the deviation from spherical symmetry. For $\vartheta = 0$, one has $K = S$ and also $\partial_\vartheta K = \partial_\vartheta S = 0$. In view of Eq.(7.15), this guarantees that the two gravitational constraints $C_1$ and $C_2$ vanish at the symmetry axis, and Eq.(7.14) then implies that they should vanish everywhere, which is indeed confirmed with a good precision. All three metric functions U, K, S show a weak $\vartheta$-dependence at the horizon, but as required by Eq.(7.20), the combination $2U-K$, which determines the surface gravity,

$$k_H = \frac{N'}{2} e^{2U-K} \Big|_{r=r_H} = \frac{1}{2r_H} e^{2U-K} \Big|_{r=r_H}, \qquad (8.3)$$

should be constant at the horizon. This is indeed confirmed by the numerics (see Fig.13) and shows that the constraints are fulfilled.

For the solution with $r_H = 0.7$ displayed in Figs.7–9, one has $(2U - K)|_{r=r_H} = -0.0377$. Injecting this to (8.3) yields the value $k_H = 0.6878$.

When the horizon size decreases, the surface gravity first increases but then reaches a maximum value and start decreasing, approaching zero in the extremal limit (see Fig.13).
<figure><img src="image_11.png" /><figcaption>Figure 10. The profiles of U, S − U, K − S, and $H_1$ for the extremal hairy black hole with $\kappa = 10^{-3}$ and $n = 4$.</figcaption></figure>

In the extremal limit, $r_H \to 0$, the factor $1/r_H$ in (8.3) diverges, but its growth is damped by the exponent $e^{2U-K}$, which approaches zero because the horizon value of $2U - K$ becomes large and negative. Therefore, as in the spherically symmetric case, the horizon values of U and K grow without bound as the extremal limit is approached. The same conclusion holds for the S amplitude.

To analyze the extremal limit, we choose $N(r)$ according to (8.2). Then, $N'(r_H) = 0$, and the surface gravity vanishes without requiring large values of the metric amplitudes at the horizon. The profiles of the extremal solution with $n = 4$ and $\kappa = 10^{-3}$ are shown in Figs. 10–12. The most notable change compared to the non-extremal solution in Figs. 7–9 is in the profiles of U and S − U. They show a non-monotonic dependence on $r$ and become very similar to those for the extremal spherically symmetric solution with $n = 2$ shown in Fig.24 in Appendix B (for the non-extremal solution, one has $S − U \approx −U$).

Another important change is that all amplitudes are now constant at the horizon: S − U, K − S, $H_1$, $H_3$, y, $\phi_1$, and $\phi_2$ all vanish, $H_2$ and $H_4$ equal to unity, whereas U assumes a
<figure><img src="image_12.png" /><figcaption>Figure 11. The profiles of $H_2$, $H_3$, $H_4$, and $y$ for the extremal hairy black hole with $\kappa = 10^{-3}$ and $n = 4$.</figcaption></figure>

constant value. Therefore, the electroweak fields are in the false vacuum, $W_\mu^a = \Phi = 0$ (see Eqs.(7.29),(7.31)), and the horizon geometry is perfectly spherical, as will be shown in the next section. Thus, the horizon geometry coincides with that of the extremal RNdS solution (3.7),(3.25). This motivates the gauge choice (8.2) with $r_H = r_{ex}$ given by (3.23).

We have described above hairy solutions for $n = 2$ and $n = 4$. The $n = 3$ solution is qualitatively the same as the one for $n = 4$. However, one may also wonder what happens when $n = 1$. The perturbation theory developed in Sections IV and VI is valid only for $|n| > 1$, since for $|n| = 1$, the perturbations are unbounded. This implies that there are no hairy solutions with $|n| = 1$ bifurcating from the RN family. However, there could be hairy solutions that do not intersect with the RN family. To detect them, we start from the hairy solutions with $n = 2$ and then iteratively decrease $n$ down to $n = 1$. The numerical procedure converges, but the properties of the result, such as the virial relation, indicate that this is not a true solution but rather a numerical artifact. Therefore, it seems that there are no hairy solutions with $n = \pm 1$ within the axially symmetric ansatz.
<figure><img src="image_13.png" /><figcaption>Figure 12. The Higgs amplitudes $\phi_1$ and $\phi_2$ for the extremal hairy black hole with $\kappa = 10^{-3}$ and $n = 4$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Left: the angular dependence of 2U–K, U, and S–U at the horizon for the hairy black hole solution with $\kappa = 10^{-3}$, $r_H = 0.7$, and $n = 4$. Right: the surface gravity as a function of $r_H$ for the hairy solutions with $\kappa = 10^{-3}$, $n = 4$, and for the RN solution. The maximal value of the surface gravity, and hence of the Hawking temperature $T = k_H/(2\pi)$, is higher for the hairy solutions than for the RN solution.</figcaption></figure>

## IX. PROPERTIES OF HAIRY BLACK HOLES

Let us now provide a more general description of hairy solutions for higher values of the magnetic charge. Non-extremal hairy black hole form a two-parameter family, characterized by an integer $n$ that determines the magnetic charge and a real parameter $r_H$ that sets the horizon size. At large distances from the horizon, the fields asymptotically approach those of the RN solution with magnetic charge $P = n/(2e)$. To begin with, we assume that the charge is not too large, such that $|Q| = \sqrt{\kappa/2} |P| \ll Q_m$, where $Q_m$ is defined in (3.24).

We emphasize once again that the size parameter $r_H$ does not have a direct physical
<figure><img src="image_15.png" /><figcaption>Figure 14. As the horizon size decreases, the hair first emerges (left), then extends further (center), and eventually, a bubble of symmetric phase forms around the horizon (right).</figcaption></figure>

meaning but merely serves as a label for the solutions. Its value determines the event horizon area $A_H$ according to (7.19), and as $r_H$ decreases, the horizon area shrinks.

When the size parameter reaches its maximal value, $r_H = r_H^0(n)$, the hairy solution reduces to the RN black hole of charge $Q$ and horizon radius $r_h = r_h^0(n)$. The values of $r_h^0(n)$ for small $n$ are shown in Table II (for $\kappa = 5.3 \times 10^{-33}$), while for large $|n|$, we find that $r_h^0(n) \approx \sqrt{|n|}/g$.

The relation between the maximal value $r_H^0(n)$ and the Schwarzschild radial coordinate $r_h^0(n)$ of the RN horizon is given by

$$r_H^0(n) = r_h^0(n) - \frac{Q^2}{r_h^0(n)} = r_h^0(n) + \mathcal{O}(\kappa). \qquad (9.1)$$

For this value of $r_H$, the strength of the hypermagnetic field at the horizon is just enough to trigger the electroweak condensation. In this regime, the condensate can be described perturbatively by linearizing the W-equations around the RN background, which leads to Eq. (4.3). The solution to this equation, given by Eqs. (4.5), (4.7), approaches zero outside the horizon more rapidly as $|n|$ increases (see Fig. 2). Hence, close to the bifurcation point and for large $|n|$, the electroweak condensate is strongly localized near the horizon, as illustrated schematically in Fig. 14 (left panel).

When the size $r_H$ becomes smaller than $r_H^0(n)$, the condensate begins to backreact, changing the geometry so that it is no longer purely RN. In this regime, the electroweak condensate extends further, spanning the interval $r \in [r_H, r_H^0]$, where the hypermagnetic field satisfies the condition $|B| > m_W^2$. The condensate terminates at approximately $r \approx r_H^0(n) \approx r_h^0(n)$, where the hypermagnetic field reduces to the threshold value $|B| \approx m_W^2$ (Fig. 14, central panel).

When $r_H$ decreases below a critical value $r_H^b$, where the hypermagnetic field reaches the second threshold $|B(r_H^b)| \approx m_H^2$, the massive condensate ceases to expand and remains confined within the region $r_H^b < r < r_H^0$, where the field strength satisfies $m_H^2 > |B| > m_W^2$. In the near-horizon region $r_H < r < r_H^b$, the hypermagnetic field becomes strong enough, $|B| > m_H^2$, to suppress the W and Higgs fields. As a result, the full electroweak gauge symmetry is restored in this region (see Fig. 14, right panel). This bubble of symmetric phase expands as $r_H$ decreases further, reaching its maximal extent in the $r_H \to 0$ limit. In this extremal limit, the horizon surface gravity vanishes, yet the horizon area remains finite.
## A. The horizon oblateness

The hairy solutions are not spherically symmetric because the condensate is compressed toward the equatorial plane. Since the condensate backreacts and influences the spacetime geometry, the event horizon itself is not perfectly spherical. As a result, there is no uniquely defined event horizon radius.

However, the horizon area $A_H$ remains a well-defined quantity and can be computed using Eq. (7.19). This allows us to define an effective average event horizon radius as

$$r_h = \sqrt{\frac{A_H}{4\pi}}. \qquad (9.2)$$

As illustrated in Fig. 15, the average radius $r_h$ (and consequently the horizon area $A_H$) is a monotonic function of $r_H$. In the extremal limit $r_H \to 0$, the radius and area approach finite values. Fig. 15 also displays the surface gravity $k_H$, computed from Eq. (7.19), which tends to zero in the extremal limit, as expected.

Since the horizon is not spherically symmetric, one can define the equatorial radius and polar radius as follows:

$$\begin{aligned} r_h^{\text{eq}} &= \sqrt{g_{\varphi\varphi}}|_{r=r_H, \vartheta=\pi/2} = r_H e^{S-U}|_{r=r_H, \vartheta=\pi/2}, \\ r_h^{\text{pl}} &= \frac{1}{\pi} \int_0^\pi \sqrt{g_{\vartheta\vartheta}}|_{r=r_H} d\vartheta = \frac{r_H}{\pi} \int_0^\pi e^{K-U}|_{r=r_H} d\vartheta. \end{aligned} \qquad (9.3)$$

This allows us to introduce the horizon oblateness, which is defined as:

$$\delta = \frac{r_h^{\text{eq}}}{r_h^{\text{pl}}} - 1. \qquad (9.4)$$

As seen in Fig. 15, the oblateness vanishes in the RN limit as expected. Interestingly, it also vanishes in the extremal limit, a feature that will be further analyzed below.

## B. The charge

Within the Nambu description, the total magnetic charge of the black hole splits into two components:

$$P = P_H + P_h = \frac{n}{2e}, \qquad (9.5)$$

where $P_h$ is the charge contained in the condensate outside the horizon. It is computed in the same way as in the spherically symmetric case in Eq. (5.7):

$$\begin{aligned} P_h &= \int_{r>r_H} \tilde{J}^0 \sqrt{-g} d^3 x = \frac{e}{8\pi} \int_{r>r_H} \epsilon^{ijk} \partial_i \psi_{jk} d^3 x = -\frac{e}{8\pi} \oint_{r=r_H} \psi_{jk} dx^j \wedge dx^k \\ &= -\frac{e}{4\pi} \oint_{r=r_H} \psi_{\vartheta\varphi} d\vartheta d\varphi \equiv \gamma g'^2 P \equiv \lambda P. \end{aligned} \qquad (9.6)$$

Here, $\psi_{\mu\nu}$ is defined in (2.35). The remaining charge, $P_H = P - P_h$, is contained inside the black hole.
The ratio $\lambda = P_h/P$ vanishes in the RN limit $r_H \to r_H^0(n)$ and reaches its maximal value of $g'^2 \approx 0.22$ in the extremal limit $r_H \to 0$ (provided that the total charge $P$ is not too large, as discussed in Section X below). In this case we have $P_h \approx 0.22 \times P$ (see Fig.15). We also introduce the parameter $\gamma = P_h/(g'^2P) = \lambda/g'^2$, which varies from zero to one. In the spherically symmetric case, this parameter corresponds to $\gamma = f_h^2$, as seen in Eq.(5.7).

In summary, hairy black holes exist for $r_H \in [0, r_H^0]$. They lose their hair in the RN limit $r_H \to r_H^0$, where all the charge is contained inside the black hole. In contrast, they become maximally hairy in the extremal limit $r_H \to 0$, with approximately 22% of the magnetic charge contained in the hair outside the horizon.

### C. The ADM mass

The ADM mass $M$ can be determined either from the asymptotic expansion of the metric coefficient, $g_{00} = -1 + 2M/r + \mathcal{O}(1/r^2)$, or from the integral expression derived using the Einstein equations in (7.17),

$$M = \frac{k_H A_H}{4\pi} + \frac{\kappa}{8\pi} \int_{r>r_H} (-T_0^0 + T_k^k) \sqrt{-g} \, d^3x, \qquad (9.7)$$

where the surface gravity $k_H$ and the horizon area $A_H$ are given by Eq.(7.19). Both methods yield consistent results with the same level of precision as for the virial relation.

As in the spherically symmetric case, the total mass can be decomposed as $M = M_H + M_h$, where $M_H$ represents the mass of the RN black hole with the same area $A_H$ and charge $P_H$. Some clarifications are necessary at this point.

The horizon radius of the RN solution with the magnetic charge $P$ and mass $M$ is given by $r_h = M + \sqrt{M^2 - Q^2}$, where $Q^2 = (\kappa/2)P^2$. Expressing $M$ in terms of $r_h$, we obtain

$$M = \frac{r_h}{2} + \frac{Q^2}{2r_h} = \frac{r_h}{2} + \frac{\kappa P^2}{4r_h}. \qquad (9.8)$$

On the other hand, the total magnetic charge is carried by both the U(1) gauge field $B_\mu$ and the SU(2) field $W_\mu^a$, leading to a decomposition: $P = P_{U(1)} + P_{SU(2)}$, where $P_{U(1)} = g^2 P$ and $P_{SU(2)} = g'^2 P$. Importantly, the square of the total charge $P^2$ in Eq.(9.8) can be represented in the form

$$P^2 = \frac{1}{g^2} P_{U(1)}^2 + \frac{1}{g'^2} P_{SU(2)}^2. \qquad (9.9)$$

This formula is crucial: although $P$ is the sum of two charges, $P^2$ should not be computed as the square of the sum. Instead, one must use Eq.(9.9) because the two charges are of different nature, and only the square of each charge contributes to the energy.

The magnetic charge of the hairy black hole splits as $P = P_H + P_h$, where $P_H$ is the charge inside the black hole, and $P_h$ is the hair charge defined by Eq.(9.6). The hair charge is always a part of the non-Abelian charge $P_{SU(2)}$, leading to the relations

$$P_h = \gamma P_{SU(2)}, \quad P_H = P_{U(1)} + (1 - \gamma) P_{SU(2)}, \qquad (9.10)$$

where $\gamma \in [0, 1]$ is the same parameter as that introduced in Eq.(9.6). Comparing with Eq.(9.9), the square of the charge contained inside the black hole is given by

$$P_H^2 = \frac{1}{g^2} P_{U(1)}^2 + \frac{(1 - \gamma)^2}{g'^2} P_{SU(2)}^2 = (g^2 + (1 - \gamma)^2 g'^2) P^2. \qquad (9.11)$$
<figure><img src="image_16.png" /><figcaption>Figure 15. Parameters of the non-extremal solutions with $n = 10$ and $\kappa = 10^{-3}$. As $r_H \to 0$, the solutions become extremal, while for $r_H \to 2.66$, they lose their hair and merge with the RN solution. The mass $M$ of extremal and near extremal solutions is less than the charge $Q$.</figcaption></figure>

Using the definition (9.2) of the average horizon radius $r_h$ and comparing with Eq.(9.8), the horizon mass is

$$M_H = \frac{r_h}{2} + \frac{\kappa P_H^2}{4r_h} = \frac{1}{2} \left( r_h + (g^2 + (1-\gamma)^2 g'^2) \frac{Q^2}{r_h} \right). \qquad (9.12)$$

This expression matches the formula in Eq.(5.12) derived in the spherically symmetric case, upon replacing $\gamma \to f_h^2$. In the RN limit, where $\gamma = 0$, the expression reduces to Eq.(9.8). In the extremal limit, where $\gamma = 1$, the average radius is $r_h = r_{ex} \approx g|Q|$, yielding

$$M_H = \frac{1}{2} \left( r_h + g^2 \frac{Q^2}{r_h} \right) \approx g|Q|. \qquad (9.13)$$

As seen in Fig.15, the mass $M$ increases with $r_H$. For $\kappa = 10^{-3}$, one finds that $M \approx M_H$ (when $|Q| \ll Q_m$). However, the relative contribution of the hair mass, $M_H/M_H$, grows with $|n|$. This contribution reaches its maximum in the extremal limit, where the black hole is maximally hairy, and vanishes in the RN limit.

## D. Quadrupole moments

Far away from the horizon, all massive fields approach their vacuum values, and one has

$$B_{\mu} = W_{\mu}^{3} \equiv eA_{\mu}, \quad W_{\mu}^{1} = W_{\mu}^{2} = 0, \quad \Phi = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \qquad (9.14)$$

thus the theory reduces to the electrovacuum:

$$\mathcal{L} = \frac{1}{2\kappa} R - \frac{1}{4} F_{\mu\nu} F^{\mu\nu}, \qquad (9.15)$$

where $F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}$. This allows us to define the gravitational $Q_G$ and magnetic $Q_M$ quadrupole moments, following the methods developed in [78–81].
The hairy solutions approach an electrovacuum configuration at large distances, which, to leading order, coincides with a magnetic RN black hole of mass M and charge $Q = \sqrt{\kappa/2} P$. However, the asymptotic behavior also includes additional subleading terms. These terms encode deviations from the pure RN solution and give rise to the multipole moments.

Choosing the gauge where $N(r) = 1$, the line element in Eq. (7.2) simplifies to:

$$ds^2 = -e^{2U}dt^2 + e^{-2U}dl^2, \qquad (9.16)$$

where the 3-metric is:

$$dl^2 \equiv h_{ik}dx^i dx^k = e^{2K}(dr^2 + r^2d\vartheta^2) + r^2e^{2S}\sin^2\vartheta d\varphi^2. \qquad (9.17)$$

The purely magnetic Maxwell field $F_{\mu\nu}$ can be expressed in terms of a magnetic potential $\Psi$ as:

$$\sqrt{\frac{\kappa}{2}} F_{ik} = e^{-2U} \sqrt{h} \epsilon_{iks} \partial^s \Psi . \qquad (9.18)$$

In the far field region, the hairy solutions reduce to:

$$e^{2U} = e^{2\delta U} \frac{x^2 - \mu^2}{(x + M)^2}, \qquad e^K = e^{\delta K} \left(1 - \frac{\mu^2}{4r^2}\right), \\ e^S = e^{\delta S} \left(1 - \frac{\mu^2}{4r^2}\right), \qquad \Psi = \frac{Q}{x + M} + \delta\Psi, \qquad (9.19)$$

where $\delta U, \delta K, \delta S, \delta \Psi$ are the subleading terms describing deviations from the RN background. The variable $x$ and the radial coordinate $r$ are related by $x = r + \mu^2/(4r)$, where $\mu^2 = M^2 - Q^2$. As seen in Fig. 15 (left panel), $\mu^2$ is not necessarily positive.

In the leading order, the expressions in (9.19) correspond to the RN solution expressed in isotropic coordinates. By solving the equations of the electrovacuum theory (9.15) in power series, we obtain the deviations:

$$\begin{aligned} \delta S &= -\frac{s_2}{2r^2} + \dots, \\ \delta K &= s_2 \frac{2 \sin^2 \vartheta - 1}{2r^2} + \dots, \\ \delta U &= \frac{Q_g (1 - 3 \cos^2 \vartheta) - Ms_2/3}{2r^3} + \dots, \\ \delta \Psi &= \frac{Q_m (1 - 3 \cos^2 \vartheta) + Qs_2/3}{2r^3} + \dots, \end{aligned} \qquad (9.20)$$

where $Q_g, Q_m$ and $s_2$ are integration constants, and the dots indicate higher-order terms. These constants can be determined from our numerical solutions, keeping in mind that the numerical solutions were obtained in a different radial gauge where $N(r) \neq 1$.

Introducing the Ernst potential [82, 83], defined as

$$\mathcal{E} = e^{2U} - \Psi^2, \qquad (9.21)$$

the gravitational $Q_G$ and magnetic $Q_M$ quadrupole moments are determined by the asymptotic behaviour of

$$\xi = \frac{1 - \mathcal{E}}{1 + \mathcal{E}}, \qquad q = \frac{2\Psi}{1 + \mathcal{E}} \qquad (9.22)$$
near spatial infinity [80, 81]. To achieve this, one must switch to the Weyl coordinates $(\rho, z)$, where the 3-metric in (9.17) becomes

$$dl^2 = e^{2K} (dr^2 + r^2 d\vartheta^2) + e^{2S} r^2 \sin^2 \vartheta d\varphi^2 = e^{2K(\rho,z)} (d\rho^2 + dz^2) + \rho^2 d\varphi^2, \quad (9.23)$$

where the function K can be expressed in terms of K, S. The next step is to consider the asymptotic expansions at the symmetry axis $(\rho = 0, z \to \infty)$:

$$\xi = \sum_{k \ge 0} \frac{a_k}{z^{k+1}}, \qquad q = \sum_{k \ge 0} \frac{b_k}{z^{k+1}}. \qquad (9.24)$$

The gravitational and magnetic quadrupole moments are then given by $Q_G = -a_2$ and $Q_M = -b_2$. Therefore, the task is to compute the coefficients $a_2$ and $b_2$ in the series (9.24).

To pass to the Weyl coordinate, we notice that since $T_r^r + T_\vartheta^\vartheta = 0$ in the electrovacuum, Eq.(7.10) with $N(r) = 1$ yields the exact solution for the amplitude S in (9.17),

$$e^S = 1 + \frac{const}{r^2} = 1 - \frac{\mu^2 + 2s_2}{4r^2}, \qquad (9.25)$$

where we have also used (9.19),(9.20). This implies that the relation between the Weyl coordinates $\rho, z$ in (9.23) and the original coordinates $r, \vartheta$ in (9.17) is provided by the complex Joukowski transformation:

$$\rho + iz = \varpi - \frac{\mu^2 + 2s_2}{4\varpi}, \quad \text{with} \quad \varpi = ire^{-i\vartheta}. \qquad (9.26)$$

Using this relation, the 3-metric in (9.17) can be transformed to the Weyl coordinates. After computing the potentials $\xi$ and $q$, expanding them according to (9.24), and using (9.19),(9.20), we obtain the expressions for the gravitational and magnetic quadrupole moments:

$$Q_G = -Q_g - \frac{2}{3} M s_2, \quad Q_M = Q_m - \frac{2}{3} Q s_2. \qquad (9.27)$$

It is important to note that the radial coordinate $r$ used in (9.20) corresponds to the gauge where $N(r) = 1$, while our numerical solutions were obtained in a different gauge where $N(r) = 1 - r_H/r$. Taking all this into account, we can determine the quadrupole moments, as shown in Fig.15.

The quadrupole moments vanish for large $r_H$ when the solution reduces to the RN. They grow as $r_H$ decreases and the solution becomes more and more hairy. They reach maximal values in the extremal limit.

## E. Condensate, vector fields and currents

In the unitary gauge, the electroweak condensate around the black hole is described by a complex vector $w^\mu$ that satisfies Eq.(6.8c). The charged condensate induces an electric current $J^\mu$, defined by Eqs.(6.5),(6.10), which flows along the horizon. The current, in turn, generates the electromagnetic and Z fields governed by (6.8a),(6.8b), forming a corona of vortices orthogonal to the horizon.
However, using the unitary gauge may lead to complications when the Higgs field approaches zero (see Appendix B in [59]), which occurs in the extremal limit. To avoid this issue, we do not impose the unitary gauge and use the fields $B_{\mu\nu}$, $W_{\mu\nu}^a$, $\Phi$ obtained from our numerical solutions. Then we compute the unit vector

$$n^a = \frac{(\Phi^\dagger \tau^a \Phi)}{(\Phi^\dagger \Phi)}, \qquad (9.28)$$

whose components remain finite when $\Phi$ approaches zero. The condensate density can be characterized by the tensor

$$\psi_{\mu\nu} = -\frac{1}{g^2} \epsilon_{abc} n^a D_\mu n^b D_\nu n^c, \qquad (9.29)$$

as defined in Eq. (2.35). Using $n^a$, $\psi_{\mu\nu}$, one can compute the electromagnetic and Z fields. These fields are described by $\mathcal{F}_{\mu\nu}$ and $\mathcal{Z}_{\mu\nu}$, as defined in Eq. (2.26) following Nambu, or by $F_{\mu\nu}$ and $Z_{\mu\nu}$, as defined in Eq. (2.36) following t'Hooft.

The divergences of $\mathcal{F}_{\mu\nu}$ and $F_{\mu\nu}$ yield, according to Eq. (2.39), the currents $J^\mu$ and $J^\mu$, respectively. The definitions of Nambu and 't Hooft lead to two distinct yet dual interpretations of the near-field structure of solutions. However, global parameters such as the total charge and quadrupole moment are determined by the far-field behaviour. Since the difference between $\mathcal{F}_{\mu\nu} - F_{\mu\nu} = e\psi_{\mu\nu}$ decays exponentially at infinity, these global quantities remain the same in both descriptions.

### 1. 't Hooft description

Let us consider first the fields $F_{\mu\nu}$ and $Z_{\mu\nu}$, and the associated current $J^\mu$. According to Eq. (6.8b), the neutral current is given by $J_Z^\mu = -(g/g')J^\mu$.

As discussed in Section VI within the perturbative approach, the current $J^\mu$ flows along the horizon, generating vortices orthogonal to it. This behaviour persists at the non-perturbative level as well. In the axially symmetric case under consideration, the current is purely azimuthal.

The electromagnetic tensor has the structure

$$\frac{1}{2} F_{\mu\nu} dx^\mu \wedge dx^\nu = P \sin\vartheta d\vartheta \wedge d\varphi + \Delta F_{r\varphi} dr \wedge d\varphi + \Delta F_{\vartheta\varphi} d\vartheta \wedge d\varphi, \qquad (9.30)$$

where the first term on the right describes the radial magnetic field surrounding the black hole in the absence of the condensate. The remaining two terms, involving $\Delta F_{r\varphi}$ and $\Delta F_{\vartheta\varphi}$, are generated by the condensate and are obtained from the numerical solution. These terms describe the magnetic field,

$$\Delta B = \frac{1}{\sqrt{-g}} (\Delta F_{\vartheta\varphi} \partial_r + \Delta F_{\varphi r} \partial_\vartheta). \qquad (9.31)$$

The azimuthal current density $J^\varphi$ is given by

$$4\pi\sqrt{-g} J^\varphi = \partial_r(\sqrt{-g} F^{\varphi r}) + \partial_\vartheta(\sqrt{-g} F^{\varphi \vartheta}). \qquad (9.32)$$

The current density changes sign across the equatorial plane and thus can be viewed as describing two distinct currents: one negative (for $n > 0$) in the upper hemisphere and one
<figure><img src="image_17.png" /><figcaption>Figure 16. Profiles of the 't Hooft magnetic field $\Delta\mathcal{B}$ (left), the 't Hooft massive magnetic field $\mathcal{B}_Z$ (center), and the Nambu magnetic field $\Delta\mathcal{B}$ (right) for the non-extremal hairy black hole with $n = 10$, $r_H = 1.6$ and $\kappa = 10^{-3}$. In the first two panels, the blue and red (online) spots indicate, respectively, regions of positive and negative azimuthal currents, $J^\varphi$ and $J_Z^\varphi$. In the third panel, the green spot marks the location of the maximal magnetic charge density $\mathcal{J}^0$. The field strength is represented by the arrow colour.</figcaption></figure>

positive in the lower hemisphere. As a result, two loops of oppositely directed azimuthal currents emerge. Together, these currents support the field $\Delta\mathcal{B}$, whose numerical profile is shown in Fig.16 (left panel). This field is purely solenoidal, winding around regions where the currents attain extremal values. The red spot in the upper half-space marks a negative extremum of the upper current, while the blue spot in the lower half-space indicates a positive extremum of the lower current.

The field $\Delta\mathcal{B}$ ressembles that of a finite segment of a vortex oriented upwards in the upper half-space and a similar vortex oriented downwards in the lower half-space. Each vortex originates at the horizon in the polar region, extends outwards, and eventually returns back to the horizon in the equatorial region.

When superposed with the background magnetic field $P \sin\vartheta \, d\vartheta \wedge d\varphi$, the field $\Delta\mathcal{B}$ exhibits an anti-screening effect: it enhances the background field in the equatorial region, where the condensate density $\sqrt{\psi_{\mu\nu}\psi^{\mu\nu}}$ is maximal, and screens the background field along the axis, where the condensate density vanishes. While the vortices do not change the total magnetic flux $4\pi P$, they redistribute it, making it inhomogeneous.

Integrating the current density (9.32) over a slice $r \in [r_H, \infty)$, $\vartheta \in [0, \pi/2]$, with $\varphi = \text{const}$
(compare with (2.47)) yields the total current in the upper half-space:

$$ 
\begin{aligned}
4\pi I_+ &= 4\pi \int_{r_H}^{\infty} dr \int_{0}^{\pi/2} d\vartheta \sqrt{-g} J^{\varphi} \\
&= \int_{0}^{\pi/2} d\vartheta \sqrt{-g} F^{\varphi r} \Bigg|_{r=r_H}^{r=\infty} + \int_{r_H}^{\infty} dr \sqrt{-g} F^{\varphi \vartheta} \Bigg|_{\vartheta=0}^{\vartheta=\pi/2}.
\end{aligned} 
\quad (9.33) $$

A similar computation in the lower half-space yields the current $I_- = -I_+$.

The dependence of $I_+$ on $r_H$ is shown in Fig.15. It vanishes when the solutions lose hair and bifurcate with the RN configuration, and it reaches its maximal amplitude in the extremal limit.

The direction of the currents depend on the sign of $n$:

$$ I_{\pm} = \mp \frac{n}{|n|} |I_{+}|. \quad (9.34) $$

Thus, for $P = n/(2e) > 0$, the current $I_+$ is always negative (see Fig.15), meaning it flows opposite to the $\varphi$-direction, whereas the current $I_-$ is positive, flowing along the $\varphi$-direction. This result is consistent with the perturbative limit (see Eq.(6.21) and Appendix C).

The central panel of Fig.22 shows the profile of the massive magnetic field,

$$ \mathcal{B}_Z = \frac{1}{\sqrt{-g}} (Z_{\vartheta\varphi} \partial_r + Z_{\varphi r} \partial_{\vartheta}). \quad (9.35) $$

This time, there is no need to subtract the value in the absence of the condensate because it is zero. The profile is similar to that of $\Delta\mathcal{B}$, except that it flows in the opposite direction, since the $Z$ field is sourced by $J_Z^\mu = (-g/g')J^\mu$. The vortex structure is now even more evident because the field $\mathcal{B}_Z$ is directed outwards in the polar regions, resembling two vortices oriented upwards and downwards. The vortices wind around the currents and return back to the horizon in the equatorial region. The field strength decreases rapidly some distance away from the horizon, and the arrows in Fig.16 become black where the field strength decreases by a factor of $\sim 10^{-3}$. Both $\Delta\mathcal{B}$ and $\mathcal{B}_Z$ are purely radial at the horizon, and Fig17 (left panel) shows their horizon value, along with those for the current and condensate densities.

Although $J_Z^\mu = (-g/g')J^\mu$, meaning the extrema of the two currents coincide, the fields $\Delta\mathcal{B}$ and $\mathcal{B}_Z$ in Fig.16 do not wind around exactly the same regions. This discrepancy arises because each field equation contains an additional term:

$$ 
\begin{aligned}
\partial_r(\sqrt{-g} \Delta F^{\varphi r}) + \partial_{\vartheta}(\sqrt{-g} \Delta F^{\varphi \vartheta}) &= 4\pi\sqrt{-g} J^{\varphi} + \partial_{\vartheta}(\sqrt{-g} P \sin \vartheta g^{\vartheta \vartheta} g^{\varphi \varphi}), \\
\partial_r(\sqrt{-g} Z^{\varphi r}) + \partial_{\vartheta}(\sqrt{-g} Z^{\varphi \vartheta}) &= 4\pi\sqrt{-g} J_Z^{\varphi} - \sqrt{-g} \frac{(\Phi^{\dagger}\Phi)}{2} Z^{\varphi}.
\end{aligned} 
\quad (9.36) $$

To summarize, within the 't Hooft description, the electroweak condensate generates azimuthal currents that flow along the horizon. These currents create a magnetic field, which forms two vortices orthogonal to the horizon, anchored in the polar and equatorial regions. These vortices cause the magnetic field to become inhomogeneous, leading to the formations of the black hole "corona".
<figure><img src="image_18.png" /><figcaption>Figure 17. The horizon values of the 't Hooft fields, current density, and condensate for the non-extremal hairy solution with $n = 10$, $r_H = 1.6$ and $\kappa = 10^{-3}$ (left) and for the extremal solution with $n = 10$ (right). The hatted indices denote projections on vectors of the orthonormal tetrad.</figcaption></figure>

## 2. Nambu description

In the Nambu case, the massless and massive magnetic fields are given by:

$$ \mathcal{B}^i = \frac{1}{2\sqrt{-g}} \epsilon^{ijk} \mathcal{F}_{jk}, \qquad \mathcal{B}_Z^i = \frac{1}{2\sqrt{-g}} \epsilon^{ijk} \mathcal{Z}_{jk}. \qquad (9.37) $$

These fields are not purely solenoidal, as their divergences do not vanish (see Eq. (2.49)):

$$ \nabla_k \mathcal{B}^k = \frac{e}{2\sqrt{-g}} \epsilon^{ijk} \partial_i \psi_{jk}, \qquad \nabla_k \mathcal{B}_Z^k = -\frac{g^2}{2\sqrt{-g}} \epsilon^{ijk} \partial_i \psi_{jk}. \qquad (9.38) $$

The Nambu field generated by the condensate, $\Delta\mathcal{B}$, is defined similarly to the previous case, with the substitution $F \to \mathcal{F}$ in Eqs. (9.30), (9.31).

As seen in the right panel of Fig. 16, the field $\Delta\mathcal{B}$ emanates from the region where the magnetic charge density is maximal (green spot in the figure) and then converges to the horizon. Thus, taking into account the rotational invariance, the magnetic charge $P_h$ is distributed over a ring – a toroidal region outside the event horizon. This ring is especially pronounced in the extremal limit, as represented by the green area in Fig. 21 below, where the level surfaces of $\sqrt{-g} \tilde{J}^0$ are shown. The magnetic field created by the ring produces a negative flux through the horizon, ensuring that the latter contains only a part of the total magnetic charge.

The Nambu field contains also a solenoidal part associated with the electric current, defined by $\nabla_k \mathcal{F}^{ik} = 4\pi \mathcal{J}^i$. However, a similar relation for a $\mathcal{Z}$ current, analogous to that in (6.8), cannot be obtained since this would require a potential for $\mathcal{B}_Z$, which does not exist since $\nabla_k \mathcal{B}_Z^k \neq 0$. A potential for $\mathcal{B}$ does not exist either, but the current $\mathcal{J}$ is defined without the use of a potential.

The current density forms two rings of oppositely directed currents, which “sandwich” the magnetically charged ring (see Fig. 21). The effect of the Nambu currents is subleading compared to that of the magnetic charge density. As seen in Fig. 15, the integrated current $\mathcal{I}_+$ is much weaker than the 't Hooft current $I_+$.
One should note that the curve $\mathcal{I}_+(r_H)$ in Fig.15 differs from that shown in Fig.3 of our previous report [57]. The difference arises because $\mathcal{I}_+$ in [57] was evaluated using the surface integral, the same one as in the first line in (9.33) but with $J^\varphi \to \mathcal{J}^\varphi$, which contains second derivatives, leading to large numerical errors. Using instead the contour integrals as in the second line in (9.33) increases the precision, yielding $\mathcal{I}_+(r_H)$ shown in Fig.15. On the contrary, to compute the 't Hooft current $I_+(r_H)$, we always use the contour integral.

Summarizing, within the 't Hooft approach, all the magnetic charge is contained inside the black hole. However, there are strong azimuthal currents in the condensate that generate a corona made of vortices. In contrast, within the Nambu description, a portion of the magnetic charge is contained in a ring outside the horizon, which generates a magnetic field which "imitate" the vortex structure, as can be seen by comparing the left and right panels in Fig.16.

# X. EXTREMAL SOLUTIONS

The extremal hairy black holes are particularly interesting because, as argued by Malda-cena, the presence of the corona should enhance the evaporation rate, causing non-extremal black holes to quickly relax to the extremal state with zero Hawking temperature [46]. For this reason, we discuss the extremal solutions in more detail.

Hairy black hole configurations are squashed toward the equatorial plane. However, as shown in Fig.15, the horizon oblateness $\delta$ vanishes in the extremal limit, indicating that the horizon becomes spherical – provided the black hole charge charge is not too large.

Specifically, the extremal solutions exhibit two distinct phases depending on the value of their charge, $Q = \sqrt{\kappa/2} \times n/(2e) \approx 0.85\sqrt{\kappa} n$, in comparison to the critical value,

$$Q_{\star} \approx 0.62 Q_{\text{m}} = \frac{0.62}{2g\sqrt{\Lambda}} \approx \frac{0.72}{\sqrt{\kappa}}. \qquad (10.1)$$

If $|Q| < Q_\star$ (or equivalently, $|n| < 0.85/\kappa$), the horizon remains spherical, and its oblateness $\delta$ vanishes. In this case, the horizon geometry coincides with that of the extremal RNdS solution (3.23). We refer to this regime as phase I.

For $|Q| > Q_\star$, the oblateness becomes nonzero, $\delta > 0$, and the horizon takes on an oblate shape, so the horizon geometry is no longer that of the RNdS solution. We refer to this regime as phase II.

Thus, the horizon oblateness $\delta$ serves as the order parameter for this phase transition. The parameter $Q_m$ in (10.1) is the same as in (3.24), but it does not represent the maximum charge for extremal hairy solutions. Instead, in (10.1) and throughout the discussion below, $Q_m$ is used as a natural scale for the charge.

The non-extremal solutions were obtained in the gauge (8.1), where $N(r) = 1 - r_H/r$. However, the fields in this gauge become unbounded in the extremal limit. To study the extremal case, we therefore switch to the radial gauge (8.2), where $N(r) \sim (1 - r_H/r)^2$. In this gauge, the horizon becomes degenerate, and the surface gravity vanishes. As explained in Appendix B, in the extremal case the size parameter $r_H$ no longer distinguishes between different solutions, but rather corresponds to different radial gauges of the same solution.

In phase I, where the horizon geometry coincides with that of the extremal RNdS solution, we use the gauge (8.2) with $r_H = r_{\text{ex}}$, where $r_{\text{ex}}$ is the radius of the extremal RNdS solution, determined by (3.23). With this choice, the radial coordinate $r$ coincides with the
Schwarzschild coordinate $r$ at the horizon, and the function $N(r)$ matches the near-horizon behavior of $N(r)$ for the extremal RNdS solution in (3.25). The horizon area is then given by $A_H = 4\pi r_{ex}^2$. Since the value of $r_H$ is fixed, the only remaining essential physical parameter for extremal solutions is the charge $n$.

In phase II, the horizon geometry is no longer spherical. Moreover, the value of $r_{ex}$, as defined by Eq.(3.24), becomes complex if the charge is large, for $|Q| > Q_m$. Therefore, we abandon the condition $r_H = r_{ex}$ and instead use the gauge (8.2) with $r_H = 0.2$.

Below, we provide a detailed description of extremal hairy solutions in phases I and II.

### A. Phase I – spherical horizon, $|Q| < Q_*$

If the total magnetic charge of the extremal black hole is not too large such that $|n| < 0.85/\kappa$, then the hypermagnetic field is strong enough to create around the horizon a bubble of false vacuum where $W_{\mu\nu}^a \approx 0$, $\Phi \approx 0$ and the electroweak symmetry is restored.

These field values are the same as those as for the extremal RNdS solution (3.25). Moreover, as shown in Fig.15, the horizon oblateness vanishes in the extremal limit, indicating that the horizon geometry is spherically symmetric – just as in the extremal RNdS case. Therefore, the field configuration approaches the extremal RNdS solution at the horizon, with radius $r_{ex}$ determined by (3.23), and we accordingly set $r_H = r_{ex}$ in Eq.(8.2).

Since $r_H = r_h = r_{ex} \approx g|Q| \approx 0.75\sqrt{\kappa}|n|$, thus, according to (3.22), the horizon value of the hypermagnetic field is

$$|B(r_h)| = \frac{|n|}{2r_h^2} \approx \frac{0.89}{\kappa|n|}, \qquad (10.2)$$

therefore, $|B(r_h)| > m_H^2$ in phase I.

The bubble of symmetric phase extends up to $r \approx r_b$, where $|B(r_b)| \approx m_H^2$, which occurs at $r_b \approx \sqrt{|n|/\beta} \approx 0.73\sqrt{|n|}$. This corresponds to the yellow area in the left panel of Fig.20 below. The Higgs field is very close to zero in this region and, as explained in Appendix E, it can be approximated by

$$\phi_1 \approx \left( \frac{r - r_{ex}}{r_{ex}} \right)^{\sigma/\lambda} S_1(\vartheta), \quad \phi_2 \approx \left( \frac{r - r_{ex}}{r_{ex}} \right)^{\sigma/\lambda} S_2(\vartheta), \qquad (10.3)$$

where $\lambda, \sigma, S_1, S_2$ are defined by (E.5), (E.13)–(E.15). Let us now analyze the situation from two different perspectives: the approaches of ’t Hooft and of Nambu.

#### 1. Structure of vortices forming the corona

Within the ’t Hooft definition, the electromagnetic and $Z$ fields are defined by Eqs.(2.35)–(2.37). They admit potentials (2.37) whose non-zero components are only azimuthal, hence $F_{j\varphi} = \partial_j A_\varphi$ and $Z_{j\varphi} = \partial_j Z_\varphi$. Subtracting from $F_{ik}$ the background contribution, that is $P \sin\vartheta d\vartheta \wedge d\varphi = -d(P \cos\vartheta d\varphi)$, and using the definition of the magnetic field in Eq.(2.42) yields

$$\sqrt{-g} \Delta \mathcal{B}^i = \epsilon^{ij\varphi} \partial_j (A_\varphi + P \cos \vartheta), \quad \sqrt{-g} \mathcal{B}_Z^i = \epsilon^{ij\varphi} \partial_j Z_\varphi. \qquad (10.4)$$
<figure><img src="image_19.png" /><figcaption>Figure 18. Electroweak corona around the black hole: the level lines of $\Delta A_{\varphi} = A_{\varphi} + P \cos\vartheta$ (left) and of $Z_{\varphi}$ (right) in “cartesian” coordinates $x = \pm r \sin\vartheta$, $z = r \cos\vartheta$ for the extremal hairy solution with $n = 10$ and $\kappa = 10^{-3}$. They coincide, respectively, with the field lines of $\sqrt{-g} \Delta B^k$ and $\sqrt{-g} B_Z^k$, and with those of $\Delta B^k$ and $B_Z^k$. The arrows show the direction of the fields.</figcaption></figure>

It follows that field lines of $\sqrt{-g} \Delta B^k$ and $\sqrt{-g} B_Z^k$ coincide with the level lines of $\Delta A_{\varphi} = A_{\varphi} + P \cos\vartheta$ and $Z_{\varphi}$ shown in Fig.18. They also coincide with lines of $\Delta B^k$ and $B_Z^k$. A disadvantage of this representation is that the density of field lines reduces in the polar regions where the factor $\sqrt{-g}$ approaches zero, whereas the actual values of the fields $\Delta B^k$ and $B_Z^k$ are not small at the poles, as seen in Fig.16.

Nevertheless, Fig.18 clearly shows the structure of the black hole corona, consisting of two vortices: one in the upper hemisphere (coloured in blue) and the other in the lower hemisphere (coloured in red).

Each vortex starts in the polar region, extends outwards, and then returns to the horizon in the equatorial region. Therefore, each vortex forms loops attached to the horizon by both ends. This is reminiscent of the coronal loops in the solar atmosphere formed by the prominences. However, unlike the latter, each vortex is not a localized jet but spreads out in all azimuthal directions, like a fountain, before returning to the horizon.

The vortices carry fluxes, but since they leave the horizon and then return, the total flux through the horizon is zero (excluding the background contribution $4\pi P$).

Each vortex carries a magnetic dipole moment, and the vortices in the upper and lower hemispheres can be viewed as two magnetic dipoles oriented against each other so that the total dipole moment vanishes (but not the quadrupole moment).

The vortices mutually repel, which is why the field lines are squeezed when approaching the equatorial region.
## 2. Vortex fluxes and the similarity with W-strings

The fluxes carried by vortices can be evaluated analytically for the extremal solutions. Specifically, the near-horizon behavior of the Higgs field expressed by Eq. (10.3) implies that the unit vector $n^a = (\Phi^\dagger \tau^a \Phi)/(\Phi^\dagger \Phi)$ depends on $\vartheta$ but is independent of $r$ (within the lowest perturbative order).

Using this information in (2.26)–(2.36), one can evaluate the non-zero components of the condensate density and magnetic field at the horizon (see Appendix E):

$$ -e\psi_{\vartheta\varphi} = g'^2 \frac{n}{2e} \times \frac{|n| \left(\sin \frac{\vartheta}{2} \cos \frac{\vartheta}{2}\right)^{|n|-1}}{\left[\left(\cos \frac{\vartheta}{2}\right)^{|n|} + \left(\sin \frac{\vartheta}{2}\right)^{|n|}\right]^2}, \qquad F_{\vartheta\varphi} = g^2 \frac{n}{2e} \sin\vartheta - e\psi_{\vartheta\varphi}. \qquad (10.5) $$

It is worth noting that, introducing $t = \tan(\vartheta/2)$, one has

$$ g^2 \psi_{\vartheta\varphi} = n \frac{d}{d\vartheta} \left( \frac{1}{1 + t^{|n|}} \right). \qquad (10.6) $$

The derivative of $F_{\vartheta\varphi}$ determines the horizon value of the current density $J_\varphi$ (see Eq. (E.23)). The flux of $F_{\vartheta\varphi}$ through the horizon is $4\pi P = 4\pi \times n/(2e)$.

One can split the electromagnetic field as in (9.30),

$$ F_{\vartheta\varphi} = P \sin\vartheta + \Delta F_{\vartheta\varphi}, \qquad (10.7) $$

where the first term on the right can be associated with the charge contained inside the black hole, while the second term can be interpreted as the magnetic field created by the condensate. One has

$$ \Delta F_{\vartheta\varphi} = -\frac{g'^2}{e} Z_{\vartheta\varphi}, \qquad Z_{\vartheta\varphi} = \frac{n}{2} \sin\vartheta + g^2 \psi_{\vartheta\varphi}. \qquad (10.8) $$

The total fluxes of $\Delta F_{\vartheta\varphi}$ and $Z_{\vartheta\varphi}$ through the horizon vanish, but locally they are not zero because $\Delta F_{\vartheta\varphi}$ and $Z_{\vartheta\varphi}$ are sign definite in the polar regions where $0 < \vartheta < \vartheta_0$ and $\pi - \vartheta_0 < \vartheta < \pi$ (see Fig. 17, right panel). They change sign at $\vartheta = \vartheta_0$, $\vartheta = \pi - \vartheta_0$, where

$$ Z_{\vartheta\varphi}(\vartheta_0) = 0. \qquad (10.9) $$

To analyze this condition, we use (10.6) and (10.8), which yields

$$ Z_{\vartheta\varphi} = n \frac{d}{d\vartheta} \left( \frac{1}{1 + t^{|n|}} - \frac{1}{1 + t^2} \right) = n(1 + t^2) \left( \frac{t}{(1 + t^2)^2} - \frac{|n|}{2} \frac{t^{|n|-1}}{(1 + t^{|n|})^2} \right), \qquad (10.10) $$

hence $Z_{\vartheta\varphi}$ changes sign at the point where

$$ t = \left( \sqrt{\frac{2}{|n|}} \frac{1 + t^{|n|}}{1 + t^2} \right)^{2/(|n| - 2)}. \qquad (10.11) $$

This equation can be solved by iterations, which gives a value $t = t_0(|n|)$. In the large $|n|$ limit, one has $t_0 \to 1$, whereas $(t_0)^{|n|} \to 0$. It follows then from (10.11) that $2|n|t_0^{|n|} \to 1$, therefore

$$ t_0 = 1 - \frac{\ln(2|n|)}{|n|} + \dots \Rightarrow \vartheta_0 = \frac{\pi}{2} - \frac{\ln(2|n|)}{|n|} + \dots, \qquad (10.12) $$
the dots denoting subleading terms. As a result, the point $\vartheta_0$ where the Z field changes sign approaches equator for large $|n|$. The flux through the northern polar region is

$$ \Psi_Z = 2\pi \int_0^{\vartheta_0} Z_{\vartheta\varphi} d\vartheta = 2\pi n c(n), \qquad (10.13) $$

where

$$ c(n) = \frac{t_0^2 - t_0^{|n|}}{(1 + t_0^2)(1 + t_0^{|n|})}. \qquad (10.14) $$

This is a monotonous function of $|n|$, such that $c(2) = 0$, $c(4) \approx 0.15$, and $c(n) \to 1/2$ for $|n| \to \infty$. This gives the fluxes leaving the horizon in the polar region:

$$ \Psi_Z = +4\pi c(n) \times \frac{n}{2}, \qquad \Psi_F = -4\pi g'^2 c(n) \times \frac{n}{2e} = -4\pi g'^2 c(n)P. \qquad (10.15) $$

Integrating over the equatorial region where $\vartheta_0 < \vartheta < \pi/2$ yields the same values, up to the sign, because fluxes leaving the horizon are the opposite of those entering the horizon.

Similarly, one obtains the fluxes in the southern hemisphere, which are the same due to the invariance of the system under the reflection $\vartheta \to \pi - \vartheta$.

The vortex fluxes arise from the electroweak condensation triggered by the background magnetic field. Using the value of the background flux through the northern hemisphere, $\Psi_F^0 \equiv 2\pi P = \pi n/e$, one can define the magnetic susceptibility as the ratio of the flux carried by the vortex to the background flux,

$$ \Upsilon(n) = \Psi_F / \Psi_F^0 = -2g'^2 c(n). \qquad (10.16) $$

The factor $c(n)$ has a geometrical origin because the vortices exist on top of the horizon which is not flat but spherical. The horizon radius increases with $|n|$, and therefore, for large $|n|$, the horizon locally approaches a plane, whereas $c(n) \to 1/2$. Therefore, for large $|n|$, one has $\Upsilon(n) \to -g'^2$ (diamagnetism), whereas the fluxes carried by each vortex are

$$ \Psi_Z = \pi n, \qquad \Psi_F = -\pi n \times \frac{g'^2}{e}. \qquad (10.17) $$

It is worth comparing these fluxes with those for the other known electroweak vortices: Z-strings [84], W-strings [85], and B-strings (see Appendix F). Even without requiring $|n|$ to be large, one obtains from (10.15) the ratio of the fluxes $\Psi_F/\Psi_Z = -g'^2/e$. This coincides with the corresponding ratio for W-strings (see Eq. (F.7)). Therefore, the vortices forming the black hole corona resemble W-strings.

One should, however, note that the similarity is not complete because W-strings interpolate along the orthogonal direction between the symmetric phase and the broken phase, while for the horizon vortices, the orthogonal direction is along the horizon, where everything is in the symmetric phase. We shall now discuss this aspect.

### 3. The origin of vortices within the symmetric phase

The non-extremal hairy black holes carry the hypermagnetic field that creates an electroweak condensate, giving rise to vortices. The condensate is made of fields that are massive because the gauge symmetry is broken. The spherical symmetry is broken as well.
For the extremal solutions, both the electroweak symmetry and spherical symmetry are restored at the horizon. The geometry there is perfectly spherical, the SU(2) field strength and the Higgs field vanish, while the hypercharge field is spherically symmetric, $B = (n/2) \sin\vartheta d\vartheta \wedge d\varphi$. However, the F and Z fields are non-spherical and form vortices, whose horizon values are described by the above formulas. According to our numerics, their amplitudes are maximal at the horizon, where the condensate density also reaches its maximal value (see Fig.19). How can this be, and what plays the role of the condensate responsible for vortex formation?

In fact, a similar issue has already been analyzed in flat space in Refs.[86, 87]. It was found that, assuming the unitary gauge and approaching the symmetric phase, where $|B| > m_H^2$ and $\Phi = W_{\mu\nu}^a = 0$, the SU(2) field potential $W_{\mu}^a$ is not zero but corresponds to a non-trivial “twisted” pure gauge. It is this pure gauge that plays the role of a massless condensate, rendering the F, Z fluxes inhomogeneous and creating vortices even in the symmetric phase.

This is precisely what happens in our case. At the horizon, where $W_{\mu\nu}^a = 0$, Eqs.(2.26), (2.36) yield

$$F_{ik} = \frac{g^2}{e} B_{ik} - e \psi_{ik}, \quad Z_{ik} = B_{ik} + g^2 \psi_{ik}, \qquad (10.18)$$

where the condensate field $\psi_{ik}$ is not spherically symmetric (see Eq.(10.5)), although the horizon geometry is.

The tensor $\psi_{\mu\nu}$ is given by the gauge invariant expression (2.35). Let us first consider the gauge (7.29) – (7.31), in which $W_{\mu}^a$ vanishes at the horizon. Then Eq.(2.35) reduces at the horizon to

$$g^2 \psi_{ik} = -\epsilon_{abc} n^a \partial_i n^b \partial_k n^c, \qquad (10.19)$$

where $n^a$ is the $\Phi \to 0$ limit of $(\Phi^\dagger \tau^a \Phi)/(\Phi^\dagger \Phi)$. Using $\Phi$ given by Eq.(7.31) with $\phi_1, \phi_2$ from (10.3) yields

$$n^a = \frac{(S_1)^2 - (S_2)^2}{(S_1)^2 + (S_2)^2} e_r^a + \frac{2S_1S_2}{(S_1)^2 + (S_2)^2} e_\vartheta^a, \qquad (10.20)$$

with unit vectors $e_r^a, e_\vartheta^a$ defined in (7.30) and $S_1(\vartheta), S_2(\vartheta)$ given by (E.13),(E.14). Inserting this into (10.19) and using the relation $\nu = -n/2$ yields exactly the same components of the condensate density $\psi_{\vartheta\varphi} = -\psi_{\varphi\vartheta}$ as in (10.5).

One notices that $\psi_{ik}$ in (10.19) is proportional to the volume form on the two-sphere spanned by the unit vector $n^a$. When $\vartheta, \varphi$ run over the horizon, the unit vector $n^a$ spans the two-sphere $(-\nu)$ times, where the minus sign arises because $n^a$ runs in the opposite direction, since one has $n^a = \mp \delta_3^a$ for $\vartheta = 0, \pi$, respectively. This corresponds to the “twist” of the condensate.

The horizon flux of $\psi_{ik}$ is therefore proportional to $(-\nu) \times 4\pi$,

$$e \oint \psi_{\vartheta\varphi} d\vartheta d\varphi = e \times \frac{4\pi\nu}{g^2} = -g'^2 \times 4\pi \frac{n}{2e} = -4\pi g'^2 P. \qquad (10.21)$$

This reproduces Eq.(E.20) in Appendix E and ensures that the flux of $F_{\vartheta\varphi}$, defined according to (10.18), is equal to $4\pi P$, while the flux of $Z_{\vartheta\varphi}$ is zero.

However, the gauge in which $W_{\mu}^a$ vanishes at the horizon is not unitary. The Higgs field is given in this gauge by Eq.(7.31), where both the lower and upper components do not
vanish, although they approach zero at the horizon. To pass to the unitary gauge where $\Phi^T = (0, \phi)$, one can first apply the transformation inverse to that in (7.28), which brings the fields back to the form (7.22), where $\Phi^T = (\phi_1, \phi_2)$.

Next, one applies another gauge transformation generated by $\exp(2i\alpha T_2)$, whose effect on the upper component of the Higgs field is $\phi_1 \rightarrow \tilde{\phi}_1 = \phi_1 \cos\alpha + \phi_2 \sin\alpha$. Requiring that $\tilde{\phi}_1 = 0$ yields

$$ \tan(\alpha) = -\frac{\phi_1}{\phi_2} = -\frac{S_1(\vartheta)}{S_2(\vartheta)} = \frac{2}{|\nu|+1} \frac{S'_2(\vartheta)}{S_2(\vartheta)} = t \frac{t^{|\nu|-1}-1}{t^{|\nu|+1}+1}, \quad (10.22) $$

with $t = \tan(\vartheta/2)$ and we used the formula (10.3) for the Higgs field in the horizon vicinity, and also Eqs.(E.13),(E.14) from Appendix E.

As a result, the gauge transformation bringing the fields (7.29) – (7.31) to the unitary gauge is generated by

$$ \mathcal{U} = e^{i(\vartheta+2\alpha)T_2} e^{i\nu\varphi T_3} \quad (10.23) $$

(we have omitted the U(1) part). The horizon value of the SU(2) field in this gauge is

$$ \begin{aligned} T_a W^a_\mu dx^\mu &= i \mathcal{U} \partial_\mu \mathcal{U}^{-1} dx^\mu \\ &= T_2 d(2\alpha + \vartheta) + \nu \{\cos(2\alpha + \vartheta)T_3 - \sin(2\alpha + \vartheta)T_1\} d\varphi. \end{aligned} \quad (10.24) $$

Since in the unitary gauge we have $n^a = -\delta^a_3$, the formula (2.35) reduces to (2.38) and yields

$$ g^2 \psi_{\vartheta\varphi} = W^1_{\vartheta} W^2_{\varphi} - W^1_{\varphi} W^2_{\vartheta} = -\nu (\cos(2\alpha + \vartheta))' = n \frac{d}{d\vartheta} \left( \frac{1}{1 + t^{|n|}} \right), \quad (10.25) $$

which precisely agrees with (10.6). Therefore, the condensate density $\psi_{ik}$ responsible for the vortex formation within the symmetric phase can be viewed as a non-trivial pure gauge, in agreement with [86, 87].

To recapitulate, in the gauge defined by Eqs.(7.29),(7.31), one has $W^a_\mu = 0$ at the horizon, but the condensate is supported by the non-trivial angular dependence of the unit vector $n^a$ in Eq.(10.20). Upon passing to the unitary gauge, the unit vector $n^a$ becomes trivial, while $W^a_\mu$ becomes the pure gauge configuration given by Eq.(10.24), inheriting the non-trivial angular dependence.

### 4. Agreement between the different descriptions

One more issue remains to be addressed. The equatorial region, where the vortices converge back to the horizon, $\vartheta_0 \le \vartheta \le \pi - \vartheta_0$, is notably narrow: its width decreases as $\ln(|n|)/|n|$ when $|n|$ increases (see (10.12)). At the same time, the fluxes (10.17) scale as $|n|$, hence the vortex fields in this region grow as $n^2/\ln(|n|)$. As illustrated in Fig.18, the density of field lines becomes very large in the equatorial region. This has interesting consequences.

Horizon vortices bear a resemblance to W-strings. However, the latter can terminate, and when they do, the trapped magnetic flux $\Psi_F$ spreads outward from the termination point, generating a radial magnetic field that mimics a magnetic monopole with charge proportional to $g'^2/e$ [61]. Based on this, Maldacena proposed that the black hole corona consists of vortices attached to the horizon at one end and terminating on monopoles some
distance away, as depicted in Fig.4b of [46]. This interpretation is supported by Nambu's approach, which suggests that a portion of the magnetic charge resides outside the horizon.

At the same time, we know that the corona rather consists of vortex segments attached to the horizon at both ends, forming loops. This might appear to contradict the interpretation of [46]. However, it turns out that the two interpretations are not mutually exclusive and can, in fact, be reconciled.

Specifically, let us compute the outgoing flux through the equatorial region of the total magnetic field, which includes both the background contribution and the vortex field, i.e., $F_{\vartheta\varphi} = P \sin\vartheta + \Delta F_{\vartheta\varphi}$ (see (10.7)). Assuming that $|n| \gg 1$ and using (10.12), (10.17), we obtain

$$ \begin{aligned} \Psi_{\text{equator}} &= 2\pi \int_{\vartheta_0}^{\pi-\vartheta_0} (P \sin\vartheta + \Delta F_{\vartheta\varphi}) d\vartheta = 4\pi P \cos(\vartheta_0) - 2\Psi_F \\ &= 4\pi P \left( g'^2 + \frac{\ln(2|n|)}{|n|} + \dots \right) \approx 4\pi g'^2 P, \end{aligned} \quad (10.26) $$

where the dots denote subleading terms. Since the vortex field $\Delta F_{\vartheta\varphi}$ is very large, it provides the dominant contribution, while the effect of the background field $P \sin\vartheta$ is subleading.

As a result, the total flux emitted by the narrow equatorial region corresponds to that of a magnetic monopole with charge $g'^2 P \propto g'^2/e$, aligning with the interpretation of [46]. Due to their strong confinement in the equatorial region, the vortex fields mimic an “effective” magnetic charge $g'^2 P$ at the horizon, even though within the ’t Hooft framework, all of the charge remain confined inside the horizon.

This effective charge is not localized at a single point but is instead distributed along a ring formed by the intersection of the horizon and the equatorial plane, as the flux $\Psi_{\text{equator}}$ remains uniform along the azimuthal direction. One can therefore interpret the equatorial region of the event horizon as being effectively charged. Remarkably, this aligns with the Nambu perspective, which posits that the $g'^2 P$ portion of the total magnetic charge $P$ resides in a ring outside the horizon (see below). As a result, what the Nambu approach considers to be the actual charge $g'^2 P$ is, in the ’t Hooft approach, effectively mimicked by the distribution of field lines, producing an equivalent charge of the same strength. The only distinction lies in their locations: in the ’t Hooft approach, the charge is located at the horizon, whereas in the Nambu approach, it resides some distance away – a point that will be discussed below. While the latter perspective aligns more closely with that of [46], it appears to be a feature of the Nambu framework, which is less fundamental than the ’t Hooft description.

It is instructive to revisit the perturbative picture shown in Fig.6. The right panel illustrates the axially symmetric solution, where the yellow equatorial zone appears charged within the Nambu approach. We now understand that this picture extends beyond the perturbative regime, showing that the equatorial zone corresponds to the region of the horizon that is effectively charged within the ’t Hooft framework.

The left panel in Fig.6 depicts a non symmetric minimum energy solution. At present, we have not tried to extend this solution at the non-perturbative level. However, if such a generalization exists, then the six yellow zones in the figure should likewise become effectively charged regions of the horizon, corresponding to monopoles. Consequently, the non-perturbative extension of the non-symmetric solution in Fig.6 should feature eight vortices, interpolating between to the dark zones, where the field lines leave the horizon, and the six monopoles, where the field lines converge back.
<figure><img src="image_20.png" /><figcaption>Figure 19. The norms of the condensate $\sqrt{\psi_{\mu\nu}\psi^{\mu\nu}}$ (left), of the Z field $\sqrt{Z_{\mu\nu}Z^{\mu\nu}}$ (center), and the azimuthal current density $J^\hat{\varphi}$ (right) for the extremal solutions with $n = 10$ and $\kappa = 10^{-3}$.</figcaption></figure>

These observations establish a connection between the 't Hooft perspective, according to which all the magnetic charge is confined within the black hole, and the Nambu perspective, which suggests that part of the charge resides outside. Furthermore, they provide a link between the interpretation of vortices proposed in [46] and our findings.

### 5. Summary of the 't Hooft description

Let us recapitulate everything said above. Within the 't Hooft description, the electroweak condensate around the black hole supports currents flowing along the horizon. The condensate and the current attain their maximal values at the horizon (see Fig. 19). These currents create vortices attached to the horizon at both ends, resembling the coronal loops in the Sun's atmosphere (see Fig. 18). The vortices originate in the polar regions of the horizon, extend outside, spread out like fountains, and eventually return back to the horizon in the equatorial region. These vortices carry a magnetic flux proportional to $4\pi g'^2/e$, similar to W-strings. The narrow equatorial region, where the vortices converge to the horizon, can be viewed as carrying an effective charge $g'^2P$, since the field strengths are significantly high there. The vortices collectively form the black hole corona, while the condensate and corona constitute together the black hole hair.

### 6. The magnetic ring, ADM mass, and the weakness of gravity

Within the Nambu description, the electroweak condensate around the black hole is magnetically charged. For extremal solutions in phase I, the condensate carries the charge $g'^2P$, which represents 22% of the total charge, the remaining 78% being contained inside the black hole (see Section IX B above).

Although the 't Hooft current density is maximal at the horizon, the Nambu current density and the magnetic charge density both vanish there. However, as seen in Fig. 20,
<figure><img src="image_21.png" /><figcaption>Figure 20. The norm of the Higgs field $\sqrt{\Phi^\dagger\Phi}$ (left), the Nambu magnetic charge density $\sqrt{-g}\tilde{\mathcal{J}}^0$ (center), and the Nambu electric current density $\mathcal{J}^\phi$ (right) for the extremal solutions with $n = 40$ and $\kappa = 10^{-4}$.</figcaption></figure>

the magnetic charge density $\sqrt{-g}\tilde{\mathcal{J}}^0$ and the Nambu current density $\mathcal{J}^\phi$ become non-zero and achieve maximal values some distance away from the horizon, where the Higgs field deviates from zero. Therefore, one may say that, in the 't Hooft description, the charge and current are located at the horizon, whereas in the Nambu description, they are situated some distance away from it.

The 3D profiles are shown in Fig.21, where the central ring (green online) corresponds to the level surfaces of the magnetic charge density. This ring is sandwiched between two others (blue and red online) containing superconducting currents $\mathcal{I}_+$ and $\mathcal{I}_- = -\mathcal{I}_+$ represented in Fig.21 by level surfaces of $\sqrt{-g}\mathcal{J}^\phi$. The Nambu currents are, however, relatively weak, the absolute values of $\mathcal{I}_\pm$ being about two orders of magnitude smaller than those for the 't Hooft currents $I_\pm$ (see Fig.22).

The condensate terminates where $|B| \approx m_w^2$, at $r_c \approx \sqrt{|n|}/g \approx 1.13\sqrt{|n|}$. Still farther away, the configuration approaches the RN solution (3.6),(3.7) whose ADM mass turns out to be smaller than the charge,

$$M < |Q|. \qquad (10.27)$$

We have already observed this phenomenon in the spherically symmetric case (Section V B). To repeat, the hair carries 22% of the total charge, hence its charge is $Q_h = 0.22 \times Q$, but its mass $M_h$ is small due to the negative Zeeman energy. This happens because the condensate contains charged W bosons whose interaction with the magnetic field of the black hole shifts their mass as $m_w^2 \rightarrow m_w^2 - |B| [48]$. As a result, one has $M_h < |Q_h|$, which explains (10.27).

If $|n|$ is not very large such that $|Q| \ll Q_*$, then the mass-to-charge ratio for the hair is very small, $M_h/Q_h \sim \sqrt{\kappa} \ll 1$. For example, for $n = 2$, one has $M_h/Q_h \sim 10^{-16}$ (see the estimate below Eq.(5.16)). The horizon contribution to the mass is then dominant, $M = M_H + M_h \approx M_H$, where $M_H$ is determined by (9.13) with $r_h = r_{ex} \approx g|Q|$. One obtains in this case the same estimate as in (5.16),

$$M \approx r_h \approx g|Q| = 0.88 \times |Q|. \qquad (10.28)$$

As $|n|$ increases, the ratio $M_h/Q_h$ also grows, and the estimate (10.28) ceases to be valid, but the inequality $M < |Q|$ continues to hold for all extremal solutions.