Black holes with electroweak hair

Romain Gervalle, Mikhail S Volkov

► To cite this version:

Romain Gervalle, Mikhail S Volkov. Black holes with electroweak hair. Physical Review Letters, 2024, 133 (17), pp.171402. 10.1103/PhysRevLett.133.171402. hal-04634772

HAL Id: hal-04634772

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

Submitted on 10 Jun 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Black Holes with Electroweak Hair

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

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

We construct static and axially symmetric magnetically charged hairy black holes in the gravity-coupled Weinberg-Salam theory. Large black holes merge with the Reissner-Nordström (RN) family, while the small ones are extremal and support a hair in the form of a ring-shaped electroweak condensate carrying superconducting W currents and up to 22% of the total magnetic charge. The extremal solutions are asymptotically RN with a mass below the total charge, $M < |Q|$, due to the negative Zeeman energy of the condensate interacting with the black hole magnetic field. Therefore, they cannot decay into RN black holes. As their charge increases, they show a phase transition when the horizon symmetry changes from spherical to oblate. At this point they have the mass typical for planetary size black holes of which ≈ 11% is stored in the hair. Being obtained within a well-tested theory, our solutions are expected to be physically relevant.

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, electric and/or magnetic charge [1]. The conjecture was supported by the uniqueness theorems in the electrovacuum sector (see [2] for a review) and by the no-hair theorems proven for a number of special cases (see, e.g., [3, 4]). The first broadly recognized counterexample to the conjecture was discovered in 1989 [5]: static black holes supporting a nontrivial Yang-Mills field without a charge (see [6] for details). This finding was generalized in many ways, leading to a plethora of black holes supporting a Skyrme hair [7], a Higgs hair [8, 9], stringy hair [10], hair without spherical symmetry [11], and many others; see [12] for a review. Some more recent examples include hairy black holes in the Horndeski theory [13], black holes supporting “spinning clouds” of ultralight bosons [14], black holes with a massive graviton hair [15], etc. [16].

Nowadays hairy black holes have become commonplace, but one may wonder if they really exist in nature? Unfortunately, all known solutions were obtained within either simplified models or exotic models relying on the existence of yet undiscovered particles and fields. Therefore, their physical relevance is not immediately obvious. To be relevant, a solution should be obtained within a well-tested physical theory, which is the standard model (SM) of fundamental interactions coupled to the general relativity (GR). The former contains the QCD sector where quantum effects are dominant, but in the electroweak (EW) 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, such solutions do exist but are difficult to construct since magnetically charged hairy EW black holes

are not spherically symmetric, unless in an unphysical limit of the theory [9], or for the low value of the magnetic charge [17]. Nonspherical black holes were considered perturbatively in theories which are similar to the EWS but without the Z boson [18], or without both Z and Higgs bosons [19]. A fully nonperturbative analysis within the EWS theory has been missing up to now.

Maldacena has qualitatively described such black holes [20]. They are not spherically symmetric because their strong magnetic field creates the EW condensate of vortices [21–24] forming a “corona” around the black hole. This happens in the region where the hypermagnetic field B falls within the interval $m_H^2 > B > m_W^2$. If $B > m_H^2$ at the horizon, then the black hole is surrounded in addition by a bubble of false vacuum. In the far field all massive fields assume vacuum values and the magnetic field becomes spherically symmetric.

Following this scenario, various aspects of magnetic black holes have been discussed, but always at the qualitative level [25–29]. Therefore, in what follows we confirm the scenario by explicitly constructing the solutions in the special case of axial symmetry when the EW condensate reduces to rings in the equatorial plane. Analyzing the general case would require 3D simulations. We sketch below only essential points and systematically refer to the Supplemental Material [30] for details; an extended version will be presented separately [31].

**The theory.** We consider the bosonic part of the EWS theory containing gravity described by the metric $g_{\mu\nu}$, a complex doublet Higgs field $\Phi$, a U(1) hypercharge field $B = B_\mu dx^\mu$, and an SU(2) field $W = T_a W_\mu^a dx^\mu$ with $T_a = \tau_a/2$ where $\tau_a$ are Pauli matrices. The action is

$$S = \frac{e^2}{4\pi\alpha} \int \left( \frac{1}{2\kappa} R - \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 \right) \sqrt{-g} d^4 x. \quad (1)$$

Here R is the Ricci scalar, $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^b W_\nu^c$ are the U(1) and SU(2) gauge field strengths, the Higgs covariant derivative is $D_\mu \Phi = (\partial_\mu - \frac{i}{2} B_\mu - \frac{i}{2} \tau_a W_\mu^a) \Phi$, and the Higgs self-coupling is $\beta \approx 1.88$. The two coupling constants

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

<footnote>† mvolkov@univ-tours.fr</footnote>
$g = \cos \theta_w$ and $g' = \sin \theta_w \approx \sqrt{0.22}$ determine the electron charge $e = gg'$ and $\alpha \approx 1/137$ is the fine structure constant. Everything is dimensionless, the Z, W, and Higgs masses are $m_Z = 1/\sqrt{2}$, $m_W = gm_Z$, $m_H = \sqrt{\beta}m_Z$ in units of the mass scale 128.9 GeV, whose Compton wavelength $1.53 \times 10^{-16}$cm sets the length scale. The gravity coupling is $\kappa = (4e^2/\alpha)(m_Z/M_{Pl})^2 = 5.30 \times 10^{-33}$ where $M_{Pl}$ is the Planck mass [30].

The electromagnetic field is $e\mathcal{F}_{\mu\nu} = g^2 B_{\mu\nu} - g'^2 n_a W_{\mu\nu}^a$ with $n_a = (\Phi^\dagger \tau_a \Phi)/(\Phi^\dagger \Phi)$ and the electric current is $4\pi\mathcal{J}^\mu = \nabla_\nu \mathcal{F}^{\mu\nu}$ [32]. Assuming the spacetime to be static and choosing the coordinates $x^\mu = (x^0, x^k)$ such that the timelike Killing vector has components $\delta_0^\mu$, the dual tensor $\tilde{\mathcal{F}}^{\mu\nu}$ defines the magnetic field $\mathcal{B}^k = \tilde{\mathcal{F}}^{0k}$ and the magnetic charge density $4\pi\tilde{\mathcal{J}}^0 = \nabla_k \mathcal{B}^k$. The hypermagnetic field is $\mathcal{B}^k = \tilde{\mathcal{B}}^{0k}$. The closed part of $\mathcal{F}_{\mu\nu}$ defines similarly a current $\mathcal{J}^\mu$ [30].

**Abelian solutions.** Varying the action gives the field equations [30], whose simplest solution corresponds to a vacuum geometry with $B = W = 0$ and $\Phi^{\text{tr}} = (0, 1)$. Choosing the gauge fields within the Cartan subalgebra, gives the Abelian solution with the RN geometry,

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

with $N(r) = 1 - 2M/r + Q^2/r^2$ and $Q^2 = (\kappa/2) P^2$ where $P = n/(2e)$ and $n \in \mathbb{Z}$. The gauge fields show the Dirac string singularity at $\vartheta = 0, \pi$, which can be removed by gauge transformations. The radial hypermagnetic field is $B = n/(2r^2)$ while the magnetic field $\mathcal{B} = P/r^2 = B/e$ corresponds to the Dirac monopole. The solution (2) describes the magnetic RN black hole of charge $P$. The event horizon is at $r_h = M + \sqrt{M^2 - Q^2}$.

Another simple solution is obtained by setting in (2) $W = \Phi = 0$ while keeping the same $B$ and choosing $N(r) = 1 - 2M/r + g^2Q^2/r^2 - (\Lambda/3)r^2$ with $\Lambda = \kappa\beta/8$. This time the magnetic field is $\mathcal{B} = g^2P/r^2 = g^2B/e$; hence, the solution describes the RN-de Sitter (RNdS) black hole of charge $g^2P$. In the extremal limit the two black hole horizons merge, then

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

where $1 - 2\Lambda r_{\text{ex}}^2 = \sqrt{1 - Q^2/Q_m^2}$ with $Q_m = 1/(2g\sqrt{\Lambda})$.

**Perturbative analysis.** Perturbing all fields around the RN solution (2) and linearizing with respect to perturbations, the field $w_\mu = (\delta W_\mu^1 + i \delta W_\mu^2)/g$ decouples from the rest (in the unitary gauge) and fulfills the equations which admit the following solution [30]:

$$w_\mu dx^\mu = e^{i\omega t} \psi(r) (\sin\vartheta)^{j+1} \sum_{m \in [-j, j]} c_m z^{m-1} dz. \quad (4)$$

Here $j = |n|/2 - 1$ is the orbital angular momentum, $|n| > 1$ (if $|n| = 1$ then $w_\mu$ is unbounded at $\vartheta = 0, \pi$),

$z = \tan(\vartheta/2) \exp(-in\varphi/|n|)$, while $c_m$ are arbitrary coefficients. The function $\psi(r)$ fulfills the 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 (5)$$

with $dr_* = dr/N(r)$. If $N(r) = 1$, this equation always admits normalizable solutions with $\omega^2 < 0$; hence, Dirac monopoles in flat space are always unstable with respect to the formation of a W condensate [33]. If $N(r)$ corresponds to the RN solution, the instability is absent for a large horizon radius $r_h$, but settles in for $r_h = r_h^0(n)$ when the equation admits a solution with $\omega = 0$ describing a static deformation of the black hole by a condensate. Therefore, $r_h^0$ is the bifurcation point where the first hairy black hole solutions deviate from the RN family.

Solving Eq. (5) numerically yields $r_h^0(n) = 0.89, 2.68, 10.29$ for $n = 2, 10, 100$, respectively. For $|n| \gg 1$ one has $r_h^0(n) \approx 1.13 \sqrt{|n|} = \sqrt{|n|}/g$; hence, $B(r_h^0) \approx m_w^2$. If $|n| = 2$ then $j = 0$ and the hair is spherically symmetric, but if $|n| > 2$ then $j > 0$ and the hair is not spherical.

The coefficients $c_m$ are fixed by going to higher orders of perturbation theory and minimizing the condensate energy [30]. This gives the result shown in Fig. 1 (left): a “platonic” distribution of radial vortex pieces forming the corona. The vortices extend along the radial direction and are encircled by loops of current $4\pi J_\nu = e\Im[\nabla^\sigma(\bar{w}_\sigma w_\nu)]$ tangent to the horizon. According to the EW anti-Lenz law [21–24], the vortex fields vanish at the vortex center where the condensate vanishes and become maximal where the condensate is maximal.

For large $|n|$ the condensate has at most only discrete symmetries, although the background magnetic field creating it is spherically symmetric. It is interesting that a similar phenomenon of symmetry breaking was observed in liquid helium, where the minimal energy defects created by applying an external field do not share the same symmetry with the field applied [34–36].

The choice $c_m = \delta_{m0}$ is also a critical point of the energy, although not the global minimum. It corresponds to the axially symmetric configuration in which all vortices

<figure><img src="image_2.png" /><figcaption>Figure 1. 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 radial vortex pieces (dark spots). The condensate and vortex fields are maximal in the yellow regions. Right: the same when all vortices merge into two oppositely directed multivortices.</figcaption></figure>
merge into two oppositely directed multivortices supported by azimuthal currents; see Fig.1 (right).

Summarizing, the perturbative analysis shows zero modes describing the hair that starts to grow. The hairy solutions are not spherically symmetric for $|n| > 2$, but they can be axially symmetric, and these we are able to construct also at the nonperturbative level.

**Axial symmetry.** We choose the fields as follows:

$$ds^2 = -e^{2U}N(r)dt^2 + e^{-2U}dl^2, \quad (6)$$

$$dl^2 = e^{2K} \left[ \frac{dr^2}{N(r)} + r^2 d\vartheta^2 \right] + e^{2S} r^2 \sin^2 \vartheta d\varphi^2,$$

$$W = T_2(F_1 dr + F_2 d\vartheta) - \frac{n}{2}(T_3 F_3 - T_1 F_4) d\varphi,$$

$$B = -(n/2) Y d\varphi, \quad \Phi^{\text{tr}} = (\phi_1, \phi_2).$$

Here U, K, S, $F_1$, $F_2$, $F_3$, $F_4$, Y, $\phi_1$, $\phi_2$ are 10 real functions of r, $\vartheta$ subject to certain boundary conditions [30]. The gauge fields W, B show the Dirac string singularity that can be gauged away if $n \in \mathbb{Z}$ [30]. They also have a residual U(1) gauge invariance that can be fixed by imposing a gauge condition [30]. The radial coordinate r (not the same as the Schwarzschild coordinate r) is defined up to reparametrizations, its choice is specified by fixing the auxiliary function N(r). For nonextremal solutions we set $N(r) = 1 - r_H/r$ where the black hole “size” $r_H$ does not have a direct physical meaning but labels the solutions. We obtain solutions for the 10 functions in (6) for values of the gravity coupling $10^{-10} < \kappa < 10^{-2}$ and then extrapolate to the physical value $\kappa \sim 10^{-33}$. We use the FreeFem numerical solver [37].

**Hairy black holes.** These black holes have the same charge as in the RN case, $P = n/(2e)$. Choosing a value $|n| > 2$ (such that $Q = \sqrt{\kappa/2} P \ll Q_m$) determines the size $r_H = r_H^0(n)$ for which the RN solution starts to develop a hair [30]. Next, we decrease the size $r_H$ and the horizon shrinks, but the massive hair grows longer spanning the interval $r \in [r_H, r_H^0]$. When $r_H$ descends below a certain value $r_H^b$ such that $B(r_H^b) \approx m_H^2$, the massive hair does not grow anymore and rests confined in the region $r_H^b < r < r_H^0$ where $m_H^2 > B > m_W^2$. In the near-horizon region, $r_H < r < r_H^b$, the hypermagnetic field becomes too strong, $B > m_H^2$, driving to zero W and $\Phi$ and creating a bubble of symmetric phase around the horizon

<figure><img src="image_3.png" /><figcaption>Figure 2. When the horizon size decreases, the massive hair emerges (left), grows longer (center), then a bubble of symmetric phase appears around the horizon (right).</figcaption></figure>

(see Fig.2). The bubble expands as $r_H$ decreases further, and in the $r_H \to 0$ limit the horizon surface gravity approaches zero but the horizon area remains finite.

Using the magnetic current density, one can compute the charge contained in the hair, $P_h = \int_{r>r_H} \tilde{J}^0 \sqrt{-g} d^3x$, while the rest $P_H = P - P_h$ is inside the black hole. Defining $\lambda = P_h/P$, we observe that $\lambda \to 0$ at the bifurcation point while in the extremal limit one has $\lambda \to g'^2$; hence, $P_h = g'^2 P = 0.22 \times P$ (see Fig.3).

Summarizing, hairy black holes exist for $r_H \in [0, r_H^0]$. They lose hair in the RN limit $r_H \to r_H^0$ when all the charge is inside the black hole, and they become maximally hairy in the extremal limit $r_H \to 0$ when 22% of the magnetic charge moves to the hair.

The mass $M$ is determined from the asymptotic expansion $g_{00} = -1 + 2M/r + \mathcal{O}(1/r^2)$ or from the formula

$$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. \quad (7)$$

Here $k_H = (1/2) N' e^{2U-K}|_{r=r_H}$ is the surface gravity and $A_H = 2\pi r_H^2 \int_0^\pi e^{K+S-2U} \sin\vartheta d\vartheta|_{r=r_H}$ is the horizon area, and $T_\nu^\mu$ is the stress-energy tensor. The total mass splits as $M = M_H + M_h$ where the horizon contribution $M_H$ is defined as the mass of the RN black hole with the same area $A_H$ and with the charge $P_H$ [30]. The mass $M$ decreases with decreasing $r_H$ and, unless $|n|$ is very large, one always has $M \approx M_H$. The relative contribution of the hair mass $M_h/M_H$ grows with $n$, it is maximal in the extremal limit and vanishes in the RN limit (see Fig.3).

We define the horizon radius as $r_h = \sqrt{A_H/(4\pi)}$, but since the horizon is nonspherical, one should define also the equatorial radius $r_h^{\text{eq}} = \sqrt{g_{\varphi\varphi}(r_H, \pi/2)}$ and the polar radius $r_h^{\text{pl}} = (1/\pi) \int_0^\pi \sqrt{g_{\vartheta\vartheta}(r_H, \vartheta)} d\vartheta$, which determine the horizon oblateness $\delta = r_h^{\text{eq}}/r_h^{\text{pl}} - 1$.

Far away from the horizon the massive fields decay and the theory reduces to electrovacuum; hence, one can define [30] the gravitational $Q_G$ and magnetic $Q_M$ quadrupole moments following the approach of [38–41]. As seen in Fig.3, the quadrupole moments grow as $r_H$ decreases, reaching maximal values in the extremal limit.

<figure><img src="image_4.png" /><figcaption>Figure 3. Parameters of nonextremal solutions with $n = 10$, $\kappa = 10^{-3}$. For $r_H \to 0$ they become extremal, for $r_H \to 2.66$ they lose hair and become RN.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. The extremal solutions contain a small charged black hole inside a bubble of symmetric phase, surrounded by a ring-shaped EW condensate supporting 22% of the total magnetic charge and two opposite superconducting W currents. These create pieces of two magnetic multivortices along the positive and negative z directions. Farther away the condensate disappears and the magnetic field becomes radial.</figcaption></figure>

The horizon oblateness $\delta$ also grows at first but then decreases and vanishes in the extremal limit, when the horizon geometry becomes perfectly spherical, although the hair is oblate.

**Extremal solutions.** Extremal black holes have zero surface gravity. The value of their charge $Q_* \approx 0.6 Q_m$ with $Q_m = 1/(2g\sqrt{\Lambda})$ separates two phases, with the horizon oblateness $\delta$ being the order parameter.

**Phase I:** $|Q| < Q_*$, $\delta = 0$. The solutions interpolate between the extremal RNdS described by Eq.(3) at the horizon and the RN given by Eq.(2) in the asymptotic region. We obtain them by choosing the auxiliary metric function $N(r) \propto (1 - r_{ex}/r)^2$ [30]. The horizon radius is $r_{ex} = g|Q| + \mathcal{O}(|n|^{3/5/2})$ and the area $A_H = 4\pi r_{ex}^2$.

The configurations harbor at the center a small black hole containing 78% of the total magnetic charge and supporting a strong B field which creates a bubble of symmetric phase where $W \approx \Phi \approx 0$. The bubble spans the region of size $r_b \approx \sqrt{|n|/\beta} \approx 0.73\sqrt{|n|}$ where $B > m_H^2$. Farther away, where $m_H^2 > B > m_W^2$, the massive fields deviate from zero and form a ring-shaped condensate, as shown in Fig.4. The central ring (green online) contains the remaining 22% of the charge whose distribution is described by the magnetic charge density $\tilde{J}^0$. This ring is sandwiched between two others (blue and red online) containing equal but oppositely directed superconducting currents $I_+$ and $I_-$ (in units of $1.4 \times 10^8$ A) which are fluxes of $\mathcal{J}^\varphi$ in the upper and lower half-spaces. Fluxes of $J^\varphi$, $I_\pm$, look similar but have a much larger amplitude (a more detailed description will be given in [31]).

The condensate terminates where $B \approx m_W^2$ at $r_c \approx \sqrt{|n|/g} \approx 1.13\sqrt{|n|}$, which determines the corona size. Still farther away the configuration approaches the RN solution (2) whose mass is less than the charge,

$$M < |Q|. \qquad (8)$$

This is because, although the hair carries the charge $Q_h = 0.22 \times Q$, its mass $M_h$ is small due to the negative Zeeman energy of the condensate interacting with the magnetic field of the black hole, which shifts the W mass as $m_W^2 \to m_W^2 - |B|$ [22]. As a result, the mass-to-charge ratio for the hair is small, $M_h/|Q_h| \sim \sqrt{\kappa} \ll 1$, which can be viewed as a manifestation of the weak grav-

ity conjecture [42]. The condensate is magnetically repelled by the black hole stronger than attracted gravitationally, but it cannot fly away because it has to obey the Yukawa law. Since the hair mass is small, one has $M = M_H + M_h \approx M_H = (r_{ex}+g^2Q^2/r_{ex})/2 \approx g|Q| < |Q|$.

Extrapolating our results to the physical value $\kappa \sim 10^{-33}$ yields the following picture. The horizon size $r_{ex} \propto Q \propto \sqrt{\kappa}$ is parametrically small as compared to the size of the hairy region $\sim \sqrt{|n|}$. Therefore, the tiny black hole does not affect the electroweak hair but only creates the B field which produces the condensate and accounts for the horizon part of the mass, $M_H = g|Q| \propto \sqrt{\kappa}$. The condensate carries the mass $M_h \sim \kappa$ and lives far away from the horizon, in the region where the geometry is almost flat. The total mass in units of the electroweak mass scale is $e^2/(4\pi\alpha) \times \mathcal{M}$, where [30]

$$\mathcal{M} \equiv \frac{8\pi}{\kappa} \times M = \frac{4\pi g |n|}{\sqrt{2\kappa} e} + \frac{8\pi}{\kappa} M_h \equiv \mathcal{M}_H + \mathcal{M}_h. \quad (9)$$

The horizon contribution $\mathcal{M}_H$ corresponds to the mass $5.1|n| M_{Pl}$, which diverges when $\kappa \to 0$, but the hair contribution $\mathcal{M}_h$ remains then finite and reduces precisely to the regularized mass of the flat space EW monopoles: the Cho-Maison monopole [43] for $|n| = 2$ and its generalizations for $|n| > 2$ [44]. As a result, extremal solutions can be viewed as flat space monopoles harboring in the center a tiny RNdS black hole. The total mass of the flat space EW monopoles is infinite; hence, it costs infinite energy to create them, but gravity renders the mass finite due to the cutoff at the horizon. Therefore, the extremal black holes provide the “classical UV completion” for the EW monopoles.

When increasing the charge $Q \propto |n|$, the bubble size and hair length scale as $\sqrt{|Q|}$, but the horizon size $r_{ex} \propto |Q|$ grows faster and the black hole absorbs the bubble. The hair mass $M_h$ grows faster than the horizon mass $M_H$ increasing the ratio $M/|Q|$. The horizon value of the B field is $\propto Q/r_{ex}^2 \propto 1/Q$; hence, it decreases, and when $B < m_H^2$, the Higgs field deviates from zero at the horizon, the bubble is totally absorbed and the solution enters next phase.

<figure><img src="image_6.png" /><figcaption>Figure 5. The parameters of extremal solutions (left) and the ratio $M/|Q|$ for the hairy solutions (right) for $\kappa = 10^{-2}$. The insertion in the right panel shows the fraction of the mass stored in the hair for the extremal solutions.</figcaption></figure>
Phase II: $|Q| > Q_\star$, $\delta > 0$. The horizon geometry changes from spherical to oblate and is no longer described by the extremal RNdS solution (3) [30]. One has $\delta \propto (|Q| - Q_\star)^s$ at the transition point, which resembles a second order phase transition (the critical exponent is rather large, $s \approx 10.8$ if $\kappa = 10^{-2}$). In addition, the fraction of the hair charge which was constant in the previous phase, $\lambda = P_h/P = 0.22$, starts to decrease (see Fig.5). The black hole absorbs the corona and becomes less hairy, the ratio $M/|Q|$ increases, the geometry approaches the extremal RN and finally merges with it when the horizon size overtakes the corona size, then $|Q| = r_h^0(n) = \sqrt{|n|/g}$. This corresponds to the maximal charge $Q_{max} = 2g'\sqrt{\beta}Q_m = 1.29Q_m$. No hairy solutions exist beyond this value of charge.

The black hole is maximally hairy for $Q \approx Q_\star$ ($|n| \approx 1.5 \times 10^{32}$) when the ratio $M_h/M$ is maximal (see Fig.5). The horizon size is then $\approx 1.37$ cm and the total mass $\approx 2 \times 10^{25}$ kg is typical for the planetary mass black holes, of which $\approx 11\%$ (if $\kappa = 10^{-2}$) is stored in the hair.

**Summary.** We obtain the hairy solutions within the region corresponding to the shaded area in Fig.5 (right panel): below the upper curve corresponding to the bifurcation with the RN branch and above the lower curve corresponding to the extremal solutions. Solutions for $|n| > 2$ are new. The corona greatly enhances the evaporation rate [20]; hence, nonextremal holes should quickly relax to the extremal state where $M \le |Q|$, therefore they cannot decay into RN black holes.

However, they can still lower their energy by developing a “hedgehog” of radial vortices forming the corona. In the axially symmetric case the corona degenerates into two oppositely directed multivortices which are unstable due to the repulsion between the elementary vortices. The axially symmetric corona can lower its energy via splitting into elementary vortices with maximal mutual separations, similarly to what is shown in Fig.1. The perturbative analysis indicates that such a “spread” configuration of vortices achieves the absolute energy minimum [30]. Therefore, the corresponding “platonic” black holes are expected to be stable, although nonaxially symmetric solutions describing such black holes are not yet known and their systematic stability analysis is still to be carried out.

Since they are described by well-tested theories, the hairy EW black holes are expected to be physically relevant. They could probably originate from primordial black holes. Suppose that the fluctuating hypermagnetic field in the ambient EW plasma becomes at some moment mostly orthogonal to the black hole horizon, or that a piece of a magnetic vortex gets attached to the horizon. This will create a flux through the horizon, hence a charge, to be compensated by the opposite flux created on the other black hole(s). The oppositely charged black holes will not necessarily annihilate, being pushed apart by the cosmic expansion, or maybe they form bound systems stabilized by a scalar repulsion [45]. Such a mechanism for the charge generation should be verified, but if

true, it would imply that no new physics beyond the SM is needed for monopoles. Therefore, the extremal hairy EW black holes are the plausible candidates for magnetic monopoles which may perhaps exist in nature.

Such black holes can show various properties, for example they should catalyze the proton decay [20]. They can be detected when captured by a neutron star, causing a sudden change of the star’s rotation period [28]. Their phenomenology can be quite rich and has been discussed in Refs. [25–29].

**Acknowledgements.** We thank Eugen Radu, Maxim Chernodub, and Julien Garaud for discussions.

# SUPPLEMENTAL MATERIAL

We give below many technical details.

## S.1. DIMENSIONFUL AND DIMENSIONLESS QUANTITIES

We denote dimensionful quantities by boldfaced letters and the dimensionless ones by ordinary letters. The action of the bosonic part of the gravity-coupled electroweak theory is

$$ \mathcal{S} = \frac{1}{c} \int \left( \frac{c^4}{16\pi G} R + L_{\text{WS}} \right) \sqrt{-g} \, d^4 x, \quad (1.1) $$

where the electroweak Lagrangian is usually chosen in the form

$$ L_{\text{WS}} = -\frac{1}{4} \mathbf{W}_{\mu\nu}^a \mathbf{W}^{a\mu\nu} - \frac{1}{4} B_{\mu\nu} B^{\mu\nu} - (D_{\mu} \mathbf{\Phi})^{\dagger} D^{\mu} \mathbf{\Phi} - \lambda (\mathbf{\Phi}^{\dagger} \mathbf{\Phi} - \mathbf{\Phi}_0^2)^2. \quad (1.2) $$

Here $\mathbf{\Phi}_0$ is the Higgs field vacuum expectation value and

$$ \begin{aligned} \mathbf{W}_{\mu\nu}^a &= \partial_\mu \mathbf{W}_\nu^a - \partial_\nu \mathbf{W}_\mu^a + g\epsilon_{abc}\mathbf{W}_\mu^b \mathbf{W}_\nu^c, \\ B_{\mu\nu} &= \partial_\mu B_\nu - \partial_\nu B_\mu, \\ D_\mu \mathbf{\Phi} &= \left( \partial_\mu - \frac{ig'}{2} B_\mu - \frac{ig}{2} \tau_a \mathbf{W}_\mu^a \right) \mathbf{\Phi}, \end{aligned} \quad (1.3) $$

with $\partial_\mu = \partial/\partial x^\mu$. Denoting $g_0 = \sqrt{g^2 + g'^2}$ one can pass to dimensionless quantities by setting

$$ \begin{aligned} g &= g_0 g, & g' &= g_0 g', & \mathbf{W}_\mu^a &= \frac{\mathbf{\Phi}_0}{g} \mathbf{W}_\mu^a, \\ B_\mu &= \frac{\mathbf{\Phi}_0}{g'} B_\mu, & \mathbf{\Phi} &= \mathbf{\Phi}_0 \mathbf{\Phi}, & \lambda &= \frac{\beta}{8} g_0^2, \end{aligned} \quad (1.4) $$

while the spacetime coordinates $x^\mu = l_0 x^\mu$ with the length scale $l_0 = 1/(g_0 \mathbf{\Phi}_0)$. The corresponding mass scale is $m_0 = \hbar/(cl_0) = (\hbar/c) g_0 \mathbf{\Phi}_0$.

One has $g' = \sin\theta_w \approx \sqrt{0.223}$ and $g = \cos\theta_w$. The dimensionful electron charge is $e = \hbar c g_0 \times e$ with
$e \equiv gg' = \sin\theta_w \cos\theta_w \approx 0.416$ and the fine structure constant

$$
\begin{align*}
\alpha &= \frac{e^2}{4\pi\hbar c} = \frac{e^2}{4\pi} \frac{\hbar c g_0^2}{\hbar c g_0^2} \approx \frac{1}{137} \Rightarrow \\
\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. \tag{1.5}
\end{align*}
$$

Injecting everything into (1.1) and introducing the di-
mensionless gravitational coupling,

$$
\kappa = \frac{8\pi\mathbf{G}\mathbf{\Phi}_0^2}{\mathbf{c}^4}, \qquad (1.6)
$$

the action becomes $\mathcal{S} = \hbar \times \mathcal{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, \quad (1.7)
$$

and

$$
\begin{equation}
\begin{split}
\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.
\end{split}
\tag{1.8}
\end{equation}
$$

This is the action in Eq.(1) in the main text.

The Z-boson mass is $m_Z = 1/\sqrt{2}$ whose dimensionful version

$$
m_Z = \frac{m_0}{\sqrt{2}} = \frac{\hbar g_0 \Phi_0}{\sqrt{2} c} = 91.18 \text{ GeV}/c^2 \quad (1.9)
$$

determines the mass scale $m_0 = 128.9 \text{ GeV}/c^2$, and the
length scale $l_0 = 1.53 \times 10^{-16} \text{ cm}$, hence

$$
\kappa = \frac{4 e^2}{\alpha} \left( \frac{m_z}{\mathbf{M}_{\mathrm{Pl}}} \right)^2 = 5.30 \times 10^{-33}. \quad (1.10)
$$

The dimensionful version of the ADM mass is determ-
ined from the asymptotic form of the g<sub>00</sub> metric coeffi-
cient,

$$
-g_{00} = 1 - \frac{2M}{r} + \dots = 1 - \frac{2GM}{c^2 r} + \dots, \quad (1.11)
$$

where $\mathbf{r} = l_0 \mathbf{r}$ hence

$$
\begin{align}
\frac{\mathbf{M}}{m_0} &= \frac{c^2 l_0}{G m_0} M = \frac{8\pi}{\kappa} \times \frac{M}{\hbar c g_0^2} \nonumber \\
&= \frac{e^2}{4\pi\alpha} \times \frac{8\pi}{\kappa} M \equiv \frac{e^2}{4\pi\alpha} \times \mathcal{M}. \tag{1.12}
\end{align}
$$

S.2. FIELD EQUATIONS

Varying the action (1.7) with respect to the EW fields
gives the equations,

$$
\begin{align}
\nabla^{\mu} B_{\mu\nu} &= g'^2 \frac{i}{2} (\Phi^{\dagger} D_{\nu} \Phi - (D_{\nu} \Phi)^{\dagger} \Phi), \nonumber \\
\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), \nonumber \\
D_{\mu} D^{\mu} \Phi - \frac{\beta}{4} (\Phi^{\dagger} \Phi - 1) \Phi &= 0, \tag{2.1}
\end{align}
$$

with $\mathcal{D}_{\mu} W_{\alpha\beta}^{a} = \nabla_{\mu} W_{\alpha\beta}^{a} + \epsilon_{abc} W_{\mu}^{b} W_{\alpha\beta}^{c}$ where $\nabla_{\mu}$ is 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.2)
$$

with the energy-momentum tensor

$$
\begin{equation}
\begin{split}
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}}.
\end{split}
\tag{2.3}
\end{equation}
$$

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, \quad \Phi = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad g_{\mu\nu} = \eta_{\mu\nu}, \quad (2.4)
$$

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 the fluctu-
ations, gives the perturbative mass spectrum containing
the massless photon, massless graviton, and the massive
Z, W and Higgs bosons with dimensionless masses

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

S.3. GAUGE TRANSFORMATIONS

The action (1.7) is invariant under SU(2)×U(1) gauge
transformations

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

with

$$
W = \frac{1}{2} (B_{\mu} + \tau^{a} W_{\mu}^{a}) dx^{\mu}, \quad \mathcal{U} = \exp \left( \frac{i}{2} \Sigma + \frac{i}{2} \tau^{a} \theta^{a} \right),
$$

where $\Sigma$ and $\theta^a$ are functions of $x^\mu$.

The ansatz for the electroweak fields in Eq.(6) in the
main text is

$$
\begin{align}
W &= T_2 (F_1 dr + F_2 d\vartheta) + \nu (T_3 F_3 - T_1 F_4) d\varphi, \nonumber \\
B &= \nu Y d\varphi, \tag{3.2} \\
\Phi &= \begin{pmatrix} \phi_1 \\ \phi_2 \end{pmatrix}, \nonumber
\end{align}
$$

with $\nu = -n/2$. It keeps its form under gauge transform-
ations generated by $\mathcal{U} = \exp\{i\chi(r, \vartheta)T_2\}$, whose effect is

$$
\begin{align}
F_1 &\rightarrow F_1 + \partial_r \chi, \quad F_2 \rightarrow F_2 + \partial_\vartheta \chi, \quad Y \rightarrow Y, \nonumber \\
F_3 &\rightarrow F_3 \cos \chi - F_4 \sin \chi, \nonumber \\
F_4 &\rightarrow F_4 \cos \chi + F_3 \sin \chi, \nonumber \\
\phi_1 &\rightarrow \phi_1 \cos(\chi/2) + \phi_2 \sin(\chi/2), \nonumber \\
\phi_2 &\rightarrow \phi_2 \cos(\chi/2) - \phi_1 \sin(\chi/2). \tag{3.3}
\end{align}
$$
The fields in (3.2) are singular at the symmetry axis. To remove the singularity, one sets

$$
\begin{align}
F_1 &= -\frac{H_1(\mathbf{r}, \vartheta)}{r\sqrt{\mathbf{N}}}, & F_2 &= H_2(\mathbf{r}, \vartheta), \nonumber \\
F_3 &= \cos \vartheta + H_3(\mathbf{r}, \vartheta) \sin \vartheta, & F_4 &= H_4(\mathbf{r}, \vartheta) \sin \vartheta, \nonumber \\
Y &= \cos \vartheta + y(\mathbf{r}, \vartheta) \sin \vartheta. & & \tag{3.4}
\end{align}
$$

The gauge transformation generated by

$$
\bar{U}_{\pm} = e^{-i\nu\varphi T_3} e^{-i\vartheta T_2} e^{\pm i\nu\varphi/2} \quad (3.5)
$$

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 (T_r H_3 + T_{\vartheta} (1 - H_4)) \sin \vartheta \, d\varphi, \quad (3.6)
$$

where T<sub>r</sub> = T<sub>1</sub> sin ϑ cos(νϕ) + T<sub>2</sub> sin ϑ sin(νϕ) + T<sub>3</sub> cos ϑ
and T<sub>ϑ</sub> = ∂<sub>ϑ</sub>T<sub>r</sub>, T<sub>ϕ</sub> = ∂<sub>ϕ</sub>T<sub>r</sub>/(ν sin ϑ). This field is ϕ-
dependent but regular at the symmetry axis, owing to
the boundary conditions defined in Eq.(6.1) below (see
[44] for details). The U(1) and Higgs fields become

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

$$
\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 hemi-spheres. These two gauges are re-
lated to each other in the equatorial region by the U(1)
transformation with $\tilde{U} = \exp(i\nu\varphi)$ (see [44] for details).
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 $\tilde{U} = \tau_3$ or
$\tilde{U} = -\tau_3$ whose effect is

$$
T_1 \rightarrow -T_1, \quad T_2 \rightarrow -T_2, \qquad (3.8)
$$

so that the field (3.6) 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 global gauge transformation generated by $\tilde{U} = \tau_3$ and $\tilde{U} = -\tau_3$, respectively. The fields (3.7) then remain invariant. As a result, all integer values of $n$ are allowed.

S.4. DIFFEOMORPHISMS

The line element in Eq.(6) of the main text

$$
\begin{align}
ds^2 &= -e^{2U} \mathrm{N}(\mathbf{r}) dt^2 + e^{-2U} dl^2, \nonumber \\
dl^2 &= e^{2K} \left[ \frac{dr^2}{\mathrm{N}(\mathbf{r})} + r^2 d\vartheta^2 \right] + e^{2S} r^2 \sin^2 \vartheta d\varphi^2, \tag{4.1}
\end{align}
$$

keeps its form under diffeomorphisms $r \to \tilde{r} = \tilde{r}(r)$ whose
effect is $U \to \tilde{U}, N \to \tilde{N}, K \to \tilde{K}, S \to \tilde{S}$, where

$$
\begin{align}
\frac{dr}{r\sqrt{\tilde{N}}} &= \frac{d\tilde{r}}{\tilde{r}\sqrt{\tilde{N}}}, & e^U \sqrt{\tilde{N}} &= e^{\tilde{U}} \sqrt{\tilde{N}}, \nonumber \\
r\sqrt{\tilde{N}}e^K &= \tilde{r}\sqrt{\tilde{N}}e^{\tilde{K}}, & r\sqrt{\tilde{N}}e^S &= \tilde{r}\sqrt{\tilde{N}}e^{\tilde{S}}. \tag{4.2}
\end{align}
$$

This symmetry can be fixed by making a specific choice for the auxiliary function N(r), for example one can set N(r) = 1. We choose N(r) = 1−r<sub>H</sub>/r for the non-extremal solutions, then r = r<sub>H</sub> corresponds to the position of the event horizon. However, since r is not the Schwarzschild coordinate, the parameter r<sub>H</sub> does not have a direct physical meaning, although its value determine the horizon area. One has r<sub>H</sub> ∈ [0, r<sub>H</sub><sup>0</sup>] and when r<sub>H</sub> → 0, the solutions approach the extremal limit with a non-zero horizon area, while for r<sub>H</sub> → r<sub>H</sub><sup>0</sup> they merge with the RN solution and one has r<sub>H</sub><sup>0</sup> = r<sub>+</sub> − r<sub>−</sub> where r<sub>±</sub> = M ± √(M<sup>2</sup> − Q<sup>2</sup>) are the two RN horizons.

For the extremal solutions we make the choice N(r) = k(r) × (1 − r<sub>H</sub>/r)<sup>2</sup> with

$$
k(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} . \quad (4.3)
$$

The horizon is again at r = r<sub>H</sub> and it is degenerate so that the surface gravity vanishes. The parameter r<sub>H</sub> again does not have a direct physical meaning and is not related anymore to the horizon area, so that, for example, one can set r<sub>H</sub> = 1. However, for the extremal solutions in phase I which have a spherical horizon it is convenient to choose r<sub>H</sub> = r<sub>ex</sub> where r<sub>ex</sub> is the radius of the extremal RNdS solution described by Eq.(3) in the main text. In this case the radial coordinate r coincides at the horizon with the Schwarzschild coordinate r and the function N(r) has the same limit as N(r) for the extremal RNdS solution. For extremal solutions in phase II one can set r<sub>H</sub> = const.

S.5. CHARGE AND CURRENTS

The electromagnetic and Z fields are defined according
to Nambu [32],

$$
\begin{align}
\mathcal{F}_{\mu\nu} &= \frac{g}{g'} B_{\mu\nu} - \frac{g'}{g} n^a W_{\mu\nu}^a, \notag \\
\mathcal{Z}_{\mu\nu} &= B_{\mu\nu} + n^a W_{\mu\nu}^a, \tag{5.1}
\end{align}
$$

where $n^a = \Phi^\dagger \tau^a \Phi / (\Phi^\dagger \Phi)$. The 2-form

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

is closed in the Higgs vacuum where $\mathcal{F}_{\mu\nu}$ reduces to the
usual electromagnetic field tensor. Away from the Higgs
vacuum the form $\mathcal{F}$ is not necessarily closed since in this
case there is no reason for Maxwell equations to apply. Defining the dual field tensor,

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

the conserved electric current and magnetic densities can be defined as

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

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

The flux of $\mathcal{F}$ through a two-sphere $S^2$ at spatial infinity defines the conserved magnetic charge,

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

Since $\mathcal{F}_{\mu\nu}$ consists of two parts determined by the contribution of $B_{\mu\nu}$ and of $W_{\mu\nu}^a$, the magnetic charge splits,

$$ P = P_{\mathrm{U}(1)} + P_{\mathrm{SU}(2)}, \qquad (5.7) $$

with

$$ \begin{aligned} P_{\mathrm{U}(1)} &= \frac{1}{4\pi} \frac{g}{g'} \oint_{S^2} \frac{1}{2} B_{\mu\nu} dx^{\mu} \wedge dx^{\nu} && (5.8) \\ &= \frac{1}{4\pi} \frac{g}{g'} \oint_{S^2} dB = \frac{g}{g'} \frac{n}{2} = g^2 \frac{n}{2e} = g^2 P, \end{aligned} $$

where we used the expression for the B-field in (3.2),(3.4) and the boundary conditions specified below in (6.1). The integral here does not depend on the radius of the sphere hence this part of the charge is always contained inside the black hole.

The remaining part of the charge,

$$ P_{\mathrm{SU}(2)} = P - P_{\mathrm{U}(1)} = g'^2 P, \qquad (5.9) $$

is not always contained inside the black hole. Its part distributed outside can be obtained by subtracting the flux at the horizon from flux at infinity,

$$ \begin{aligned} 4\pi P_h &= \oint_{r \to \infty} \mathcal{F} - \oint_{r=r_H} \mathcal{F} = \int_{r>r_H} d\mathcal{F} && (5.10) \\ &= \int_{r>r_H} \frac{1}{2} \epsilon^{ijk} \partial_i \mathcal{F}_{jk} d^3 x = \int_{r>r_H} \tilde{\mathcal{J}}^0 \sqrt{-g} d^3 x. \end{aligned} $$

Since the form $\mathcal{F}_{\mu\nu}$ is not necessarily closed, the latter integral can be non-zero. One has

$$ \mathcal{F}_{\mu\nu} = F_{\mu\nu} + e\psi_{\mu\nu}, \qquad (5.11) $$

where $F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}$ is a closed form while the rest is not closed,

$$ g^2 \psi_{\mu\nu} = -\epsilon_{abc} n^a \mathcal{D}_\mu n^b \mathcal{D}_\nu n^c, \qquad (5.12) $$

with $\mathcal{D}_\mu n^a = \partial_\mu n^a + \epsilon_{abc} W_\mu^b n^c$. This tensor is made of massive fields and approaches zero exponentially fast at infinity, therefore,

$$ \begin{aligned} 4\pi P_h &= e \int_{r>r_H} \frac{1}{2} \epsilon^{ijk} \partial_i \psi_{jk} d^3 x \\ &= -e \oint_{S^2, r=r_H} \frac{1}{2} \psi_{ik} dx^i \wedge dx^k. \end{aligned} \qquad (5.13) $$

The tensor $F_{\mu\nu}$ leads to another definition of the electric current,

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

This definition is better adapted for the perturbative analysis carried out in Section S.9 below. Integrating the current densities $\mathcal{J}^\mu$ and $J^\mu$ over a 2-surface gives currents.

In the static case, denoting $\mathcal{N} = \sqrt{-g_{00}}$, the relation (5.4) (similarly for (5.14)) can be represented in the form language as

$$ d(\mathcal{N} * \mathcal{F}) = 4\pi \mathcal{N} * \mathcal{J} \qquad (5.15) $$

where $*$ denotes the Hodge duality operator on the space-like 3-surface. Both sides of this relation are 2-forms. Integrating over a 2-surface $\Sigma$ and applying Green's theorem yields

$$ \frac{1}{4\pi} \int_{\partial\Sigma} \mathcal{N} * \mathcal{F} = \int_{\Sigma} \mathcal{N} * \mathcal{J} \equiv \mathcal{I}. \qquad (5.16) $$

Here on the left is the circulation of the magnetic field along the boundary $\partial\Sigma$, hence on the right is the total current through $\Sigma$.

In the situation considered in the main text the currents are along the azimuthal direction, hence they flow through planes of constant $\varphi$. The above definition then implies that the total azimuthal current in the upper half-space, $\vartheta < \pi/2$, is

$$ \mathcal{I}_+ = \int_{r_H}^{\infty} dr \int_{0}^{\pi/2} d\vartheta \sqrt{-g} \mathcal{J}^{\varphi}; \qquad (5.17) $$

similarly for the current $I_+$ associated with $\mathcal{J}^\varphi$. A similar integral with $\vartheta \in [\pi/2, \pi]$ yields the current $\mathcal{I}_-$ in the lower half-space. Since $\mathcal{J}^\varphi$ is antisymmetric with respect to the equatorial plane, the total current is zero, $\mathcal{I}_+ + \mathcal{I}_- = 0$. The dimensionful current expressed in amperes is

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

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

## S.6. BOUNDARY CONDITIONS AND NUMERICAL PROCEDURE

The axially symmetric EWS fields are the line element in (4.1) and the electroweak fields in (3.2),(3.4) expressed
in terms of 10 functions U, K, S, $H_1$, $H_2$, $H_3$, $H_4$, $y$, $\phi_1$, $\phi_2$ depending on r, $\vartheta$. The fields in (3.2), (3.4) have a residual U(1) gauge invariance (3.3) that we fix by imposing the gauge condition $r\sqrt{N}\partial_r H_1 = \partial_\vartheta H_2$. This gives 10 elliptic equations whose solutions fulfil also the first order gravitational constraints owing to the boundary conditions described by Eq. (6.1) below.

We use the compact radial variable $x \in [0, 1]$ such that $r = \sqrt{r_H^2 + x^2/(1-x)^2}$ and assume the reflection invariance of the energy density under $\vartheta \to \pi - \vartheta$ to restrict the range to $\vartheta \in [0, \pi/2]$. The following boundary conditions at the borders of the integration domain are imposed:

$$ \text{axis, } \vartheta = 0: \quad \partial_{\vartheta} = 0 \text{ for U, S, } H_2, H_4, \phi_2; \\ \text{K-S} = H_1 = H_3 = y = \phi_1 = 0; $$

$$ \text{equator, } \vartheta = \pi/2: \quad \partial_{\vartheta} = 0 \text{ for U, K, S, } H_2, H_4, \phi_2; \\ H_1 = H_3 = y = \phi_1 = 0; $$

$$ \text{horizon, } x = 0: \quad H_1 = 0; \quad \partial_x = 0 \text{ for the rest} \\ \text{infinity, } x = 1: \quad \phi_2 = 1; \quad 0 \text{ for the rest.} \qquad (6.1) $$

It turns out that these boundary conditions automatically imply that at the symmetry axis one has $H_2 = H_4$ and $\partial_\vartheta K = 0$, 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 (6.2) $$

which are functions of r, $\vartheta$. We solve the equations with the FreeFem numerical solver based on the finite element method [37]. This solver uses the weak form of differential equations obtained by transforming them into integral equations, expanding with respect to basis functions obtained by triangulating the integration domain, and handling the non-linearities with the Newton-Raphson procedure. The quality of the solutions is checked by monitoring the virial relation that should hold on-shell,

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

where $\zeta_\mu dx^\mu = r e^{K-U} dr$ and $T^\mu_\nu$ is the EW energy-momentum tensor (2.3). For a domain triangulation by 100 × 60 elements we typically obtain $v \sim 10^{-7}$.

## S.7. HORIZON MASS

The horizon radius of the RN solution with the magnetic charge $P$ and mass $M$ is $r_h = M + \sqrt{M^2 - Q^2}$ where $Q^2 = (\kappa/2)P^2$, hence

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

On the other hand, the total magnetic charge consists of the U(1) part and SU(2) part, $P = P_{U(1)} + P_{SU(2)}$ with $P_{U(1)} = g^2 P$ and $P_{SU(2)} = g'^2 P$ and one has

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

This formula is important: although $P$ is the sum of two charges, $P^2$ is not the square of the sum because the two charges are of different nature and only the square of each charge contributes to the energy.

The magnetic charge of a hairy black hole splits as $P = P_H + P_h$ where $P_H$ is the charge contained inside the black hole and $P_h$ is the charge contained outside in the hair defined by (5.10). The hair charge is always a part of the non-Abelian charge $P_{SU(2)}$, hence

$$ P_h = \tilde{\lambda} P_{SU(2)}, \quad P_H = P_{U(1)} + (1 - \tilde{\lambda}) P_{SU(2)}, \qquad (7.3) $$

where the parameter $\tilde{\lambda}$ varies from zero for black holes which bifurcate with RN till unity for the extremal black holes in phase I. Comparing with (7.2), the square of the charge contained inside the black hole is computed according the rule

$$ \begin{aligned} P_H^2 &= \frac{1}{g^2} P_{U(1)}^2 + \frac{(1 - \tilde{\lambda})^2}{g'^2} P_{SU(2)}^2 \\ &= (g^2 + (1 - \tilde{\lambda})^2 g'^2) P^2. \end{aligned} \qquad (7.4) $$

Therefore, defining the size of a hairy black hole with a horizon area $A_H$ as $r_h = \sqrt{A_H/(4\pi)}$ and comparing with (7.1), the horizon mass is

$$ M_H = \frac{r_h}{2} + \frac{\kappa D_H^2}{4r_h} = \frac{r_h}{2} + \frac{\kappa P^2}{4r_h} (g^2 + (1 - \tilde{\lambda})^2 g'^2). \qquad (7.5) $$

For the hairy solutions bifurcating with the RN branch one has $\tilde{\lambda} = 0$ and the formula reduces to (7.1). For the extremal hairy solutions in phase I one has $\tilde{\lambda} = 1$ and

$$ M_H = \frac{r_h}{2} + \frac{\kappa g^2 P^2}{4r_h} = \frac{r_h}{2} + g^2 \frac{Q^2}{2r_h}. \qquad (7.6) $$

Notice finally that $P_h = \tilde{\lambda} P_{SU(2)} = \tilde{\lambda} g'^2 P \equiv \lambda P$ hence $\tilde{\lambda} = \lambda/g'^2$ where $\lambda = P_h/P$ is used in the main text.

## S.8. QUADRUPOLE MOMENTS

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

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

hence the Einstein-Weinberg-Salam theory reduces to the electrovacuum,

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

where $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$. The solutions approach at large distances the magnetic RN configuration of mass $M$ and charge $Q = \sqrt{\kappa/2} P$. Choosing the gauge where $N(r) = 1$, the line element in (4.1) becomes

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

$$ dl^2 \equiv h_{ik} dx^i dx^k = e^{2K}(dr^2 + r^2 d\vartheta^2) + r^2 e^{2S} \sin^2 \vartheta d\varphi^2, $$
while 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 (8.4) $$

Our solutions in the far field region reduce to

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

with $x = r + \mu^2/(4r)$ and $\mu^2 = M^2 - Q^2$. In the leading order this is the RN solution in the isotropic coordinates, while the subleading terms describe deviations from the RN background,

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

where $Q_g, Q_m, s_2$ are integration constants and the dots denote higher order terms. The values of these constants can be determined from our numerical solutions, taking into account that the latter are obtained in a different radial gauge where $N(r) \neq 1$.

Introducing the Ernst potential,

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

the gravitational $Q_G$ and magnetic $Q_M$ quadrupole moments defined within the Geroch-Hansen approach [38, 39] are determined by the behaviour of

$$ \xi = \frac{1 - \mathcal{E}}{1 + \mathcal{E}}, \qquad q = \frac{2\Psi}{1 + \mathcal{E}} \qquad (8.8) $$

near spatial infinity. This requires [40, 41] passing to the Weyl coordinates where

$$ dl^2 = e^{2K(\rho,z)} (d\rho^2 + dz^2) + \rho^2 d\varphi^2, \qquad (8.9) $$

(the function $K$ can be expressed in terms of $K, S$) and considering the asymptotic expansions at the symmetry axis, $\rho = 0$, for $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 (8.10) $$

from where $Q_G = -a_2$, $Q_M = -b_2$. Transforming the 3-metric in (8.3) to the Weyl coordinates, computing the potentials $\xi$ and $q$ and expanding them according to (8.10) yields

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

## S.9. BEYOND AXIAL SYMMETRY

Although this is not the main subject of this work, we consider also non-axially symmetric solutions by adopting the perturbative approach of [19]. The idea is to expand around the RN solution up to the second order terms and then minimize the energy.

Passing to the unitary gauge where the Higgs field $\Phi^{\text{tr}} = (0, \phi)$ with real-valued $\phi$, defining the complex-valued vector field $w_\mu = (W_\mu^1 + iW_\mu^2)/g$ with the strength $w_{\mu\nu} = D_\mu w_\nu - D_\nu w_\mu$ where $D_\mu = \nabla_\mu + iW_\mu^3$, the EW Lagrangian (1.8) assumes the form

$$ \begin{aligned} -\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 & (9.1) \\ & + (\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. \end{aligned} $$

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

$$ \begin{aligned} W_{\mu\nu}^3 &= \partial_\mu W_\nu^3 - \partial_\nu W_\mu^3 + g^2 \psi_{\mu\nu}, \\ \psi_{\mu\nu} &= \frac{i}{2} (w_\mu \bar{w}_\nu - w_\nu \bar{w}_\mu). \end{aligned} \qquad (9.2) $$

The electromagnetic and Z-fields (5.1) then become

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

where $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$ and $Z_{\mu\nu} = \partial_\mu Z_\nu - \partial_\nu Z_\mu$ with

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

The WS equations (2.1) assume the form

$$ \begin{aligned} \nabla^\mu F_{\nu\mu} &= 4\pi J_\nu, & \nabla^\mu Z_{\nu\mu} + \frac{\phi^2}{2} Z_\nu &= -4\pi \frac{g}{g'} J_\nu, \\ 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, \\ \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. \end{aligned} \qquad (9.5) $$

Here the gauge covariant derivative can be written in the form $D_\mu = \nabla_\mu + i(eA_\mu - g^2 Z_\mu)$ and the current is

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

<figure><img src="image_7.png" /><figcaption>Figure 6. The RN horizon size $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>
It should be emphasized that $F_{\mu\nu}$ has the potential, $A_\mu$, while $F_{\mu\nu}$ given by (5.1) and used in the main text to define $J_\nu$, $\tilde{J}_\nu$ is not always closed,

$$F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu} + e\psi_{\mu\nu}. \quad (9.7)$$

This is exactly the same relation as in (5.11) because the tensor $\psi_{\mu\nu}$ defined in (5.12) reduces in the unitary gauge precisely to $\psi_{\mu\nu}$ defined in (9.2).

Taking the divergence of (9.7) and comparing with (9.6) yields

$$\begin{aligned} 4\pi 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}). \end{aligned} \quad (9.8)$$

hence

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

Eqs.(9.5) admit the solution

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

assuming the background geometry to be RN. Expanding Eqs.(9.5) around this solution yields in the first order of perturbation theory

$$D^{\mu}w_{\mu\nu} + ieF_{\nu\sigma}w^{\sigma} = \frac{g^2}{2}w_{\nu}, \quad (9.11)$$

with $D_\mu = \nabla_\mu + ieA_\mu$ where $A_\mu$ is the same as in (9.10) (the fields $Z_\mu$, $\phi$ also do not change in the first order). This equation admits a solution of the form

$$\begin{aligned} w_{\mu}dx^{\mu} &= e^{i\omega t} \sum_{m \in [-j, j]} c_m w_m(r, \vartheta, \varphi), \\ w_m &= \psi(r) (\sin \vartheta)^j \left(\tan \frac{\vartheta}{2}\right)^m e^{-im\varphi} (d\vartheta - i \sin \vartheta d\varphi), \end{aligned} \quad (9.12)$$

if $n > 1$, while for $n < -1$ one has to replace $w_m \to \bar{w}_m$. The value $|n| = 1$ is excluded since $w_m$ is then unbounded at the symmetry axis. One has $j = |n|/2 - 1$. The coefficients $c_m$ are arbitrary and can be chosen to be real. The function $\psi(r)$ fulfills the 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 (9.13)$$

with $dr_* = dr/N(r)$ where $N(r) = 1 - 2M/r + Q^2/r^2$ and $Q^2 = (\kappa n^2)/(8e^2)$. This equation admits normalizable solutions with $\omega^2 < 0$ (negative modes) if the horizon radius $r_h = M + \sqrt{M^2 - Q^2}$ is small. Therefore, small RN black holes are unstable with respect to the formation of the charged W-condensate, whereas large black holes are stable.

Varying $r_h$, one can detect critical values $r_h^k(n)$ for which instabilities just settle in as zero modes, that is normalizable solutions with $\omega = 0$. These solutions, $\psi_k(r)$,

are characterized by $k = 0, 1, 2, \dots$ nodes and they exist only for discrete values of the event horizon radius $r_h = r_h^k(n)$ shown in Fig.6. The functions $\psi_k(r)$ are finite at the horizon and decay exponentially fast at large $r$; see Fig.6. Since $\omega = 0$, the solution (9.12) with $\psi(r) = \psi_k(r)$ describes static deformations of the RN black hole, the value $k = 0$ corresponding to the fundamental mode and $k > 0$ being the radial excitations.

Setting in (9.12) $\omega = 0$, the solution can be transformed to the form expressed by Eq.(4) of the main text, while (9.13) reduces to Eq.(5) of the main text. This solution is defined up to an overall sign which can be flipped by the gauge transformation generated by $\mathcal{U} = -\tau_3$, whose effect is $T_1 \to -T_1$ and $T_2 \to -T_2$, hence $w_\mu \to -w_\mu$ while the RN background does not change. It follows that there are two options: either $w_\mu$ does not change under $\varphi \to \varphi + 2\pi$ or it flips sign. Therefore, since $w_m \sim \exp(-im\varphi)$, the azimuthal number m can assume either only integer or only half-integer values. For example, if $j = 1$ then the two options are either $m = \pm 1/2$ or $m = 0, \pm 1$.

There is always an option containing the value $m = 0$ and it is possible to choose $c_m = \delta_{m0}$ since coefficients $c_m$ are arbitrary within the linear perturbation theory. However, they are no longer arbitrary when the higher order corrections are taken into account, although the possibility to choose $c_m = \text{const.} \times \delta_{m0}$ always exists.

One has $\Im(\bar{w}^\sigma w_{\sigma\nu}) = 0$ for the solution (9.12) hence $4\pi J_\mu = e \nabla^\sigma \psi_{\sigma\mu}$ is tangential to the sphere, with non-zero components

$$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, \quad (9.14)$$

where

$$|w_{\mu}|^2 \equiv w_{\mu} \bar{w}^{\mu} = \frac{2}{r^2} |\psi(r)|^2 \times \Theta(\vartheta, \varphi), \quad (9.15)$$

$$\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}.$$

It follows from (9.14) that $J^\mu$ is orthogonal to the gradient of $|w_\mu|^2$, which is in turn orthogonal to the level lines of $|w_\mu|^2$. Therefore, the latter are parallel to $J^\mu$.

The current is quadratic in $w_\mu$, $c_m$ and sources second order correction for the electromagnetic field, $f_{\mu\nu}$, second order $Z_\mu$ and second order $\delta\phi = \phi - 1$. Expanding (9.5) up to second order terms yields

$$\nabla^{\mu} f_{\nu\mu} = e \nabla^{\sigma} \psi_{\sigma\nu}, \quad (9.16a)$$

$$\nabla^{\mu} Z_{\nu\mu} + \frac{1}{2} Z_{\nu} = -g^2 \nabla^{\sigma} \psi_{\sigma\nu}, \quad (9.16b)$$

$$\nabla^{\mu} \nabla_{\mu} \delta\phi - \frac{\beta}{2} \delta\phi = \frac{g^2}{4} |w_{\mu}|^2. \quad (9.16c)$$

Solving these equations and injecting to (9.1), yields the
result containing terms up to fourth order,

$$ -\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 (9.17) $$

Here $F_{\mu\nu} = F^{0}_{\mu\nu} + f_{\mu\nu}$ with $F^{0}_{\mu\nu}$ being the zeroth-order field corresponding to (9.10). Since the fields are static and purely magnetic, the energy density is $\mathcal{E} = -\mathcal{L}_{\text{EW}}$.

The procedure now is to solve Eqs.(9.16), compute the energy $E = \int \mathcal{E} \sqrt{-g} d^3x$, and then minimize it with respect to the coefficients $c_m$ by keeping fixed the norm $\int w_\mu \bar{w}^\mu \sqrt{-g} d^3x$. Eqs.(9.16) can be solved with Green's functions, as was done in [19] where a simplified theory without $Z$, $\Phi$ was considered. However, we postpone this for future work since our main goal at present is to study axially-symmetric solutions. To get some idea of more general solutions, we use an approximate procedure, which nevertheless seems to capture the essential points. The key observation is that the energy density always contains the fourth order term $(\psi_{\mu\nu})^2 \sim |w_\mu|^4$. Therefore, we minimize

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

under the condition that the norm should be fixed,

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

One can also argue as follows. Eq. (9.16a) is fulfilled by

$$ f_{\mu\nu} = -e \psi_{\mu\nu}. \qquad (9.20) $$

This solution is not totally satisfactory since the form $\psi_{\mu\nu}$ is not closed, but we ignore this. Then the non-zero components are

$$ f_{\vartheta\varphi} = -f_{\varphi\vartheta} = e \frac{n}{|n|} |\psi(r)|^2 \Theta(\vartheta, \varphi) \sin \vartheta, \qquad (9.21) $$

which corresponds to a radial magnetic field of the same direction as the background field. This confirms the EW anti-Lenz law: the magnetic field induced by the condensate enhances the external field and is maximal where the condensate is maximal [21–24].

Injecting (9.20) to (9.17), integrating by parts, using Eqs.(9.16), and omitting the zeroth order order term yields

$$ \mathcal{E} = \frac{g^4}{4}(\psi_{\mu\nu})^2 + \frac{1}{4}|w_{\mu\nu}|^2 + \frac{g^2}{4}|w_{\mu}|^2 \\ - \frac{g^2}{4}\psi^{\mu\nu}Z_{\mu\nu} + \frac{g^2}{4}|w_{\mu}|^2\delta\phi, \qquad (9.22) $$

where

$$ (\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. \qquad (9.23) $$

We do not know what $Z_{\mu\nu}$ and $\delta\phi$ are, but if we simply ignore them and set to zero the second line in (9.22), then the integration gives

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

The same result is obtained if we ignore only $Z_{\mu\nu}$ but use the approximate solution for $\delta\phi$ obtained by neglecting the derivatives in (9.16c), $\delta\phi = -g^2/(2\beta)|w_\mu|^2$. Of course, the energy obtained by genuinely solving the equations (9.16) should have a more complicated structure. The justification for the heuristic prescription (9.24) is that it is simple and gives what is expected – the “platonic” distribution of vortices on the horizon. Therefore, it seems that the prescription (9.24) gives correct results, even though it relies on approximations.

Using (9.15), the second term in (9.24) reduces to $\text{const.} \times E_2$ with

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

where

$$ \begin{aligned} 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)}. \end{aligned} \qquad (9.26) $$

Similarly, after some algebra, the first term in (9.24) reduces to $\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-m}. \qquad (9.27) $$

It is assumed here that $c_m = 0$ if $|m| > j$. Therefore, omitting constant factors, the problem reduces to minimization with respect to $c_m$ of

$$ E = E_4 + \mu (E_2 - 1), \qquad (9.28) $$

where $\mu$ is the Lagrange multiplier.

Assuming for example that $n = 10$ hence $j = 4$ and $m = -4, -3, \dots, 4$, so that there are 9 coefficients $c_m$ to find. Their values corresponding to the absolute minimum of $E$ determine the function $|w_\mu|^2$ in (9.15), whose value at the horizon is shown in Fig.(1) (left panel) of the main text. The figure shows the contour lines which, as explained above, coincide with the flow lines of the current $J^\mu$.

As seen in Fig.1, flow lines form loops (blue online). A loop of current generates a magnetic field in the orthogonal direction, hence the loops in Fig.1 encircle (pieces of) radially directed vortices approximately described by (9.21). The special feature of these EW vortices is that their magnetic field is enhanced by the condensate: it vanishes at the vortex center where the condensate is zero and maximal where the condensate is maximal [23, 24]. In other words, these vortices resemble hollow tubes:
everything is concentrated at the tube surface and nothing is inside.

There are 4 vortices on the visible side of the surface in Fig.1 and the same number on the other side, so that altogether there are $8 = n - 2$ vortices homogeneously distributed over the horizon. The total magnetic flux is $2\pi \times n$, and the difference $|n| - 2$ arises because for $n = 2$ the black hole is spherically symmetric [17] hence there are no vortices at all. We checked that minimization of (9.28) gives one vortex for $n = 3$, two vortices for $n = 4$, etc. These solutions correspond to the global energy minimum.

The axially symmetric configuration with $c_m \propto \delta_{m0}$ is also a stationary point of the energy (9.28), which can be

viewed as a system in which all elementary vortices merge into two oppositely directed collective multi-vortices generated by two oppositely-directed azimuthal currents.

The Hessian matrix has negative eigenvalues in this case, therefore such solutions are unstable and decay into non-axially symmetric ones. The exceptional case is $|n| = 2$ when only $m = 0$ is possible; the corresponding solution is stable in the flat space limit [33], while its gravitating version should very likely to be stable too. Another exceptional case is $|n| = 4$ when the values of $c_0, c_{\pm 1}$ corresponding to the global energy minimum can be transformed to $c_m \propto \delta_{m0}$ by a global rotation. Therefore, the axially symmetric configuration is already in the global energy minimum.

[1] R. Ruffini and J. Wheeler, Phys. Today **24**, 30 (1971), 1206.3482.

[2] M. Heusler, *Black Hole Uniqueness Theorems* (Cambridge University Press, 1996).

[3] J. D. Bekenstein, Phys. Rev. D **5**, 1239 (1972).

[4] A. Saa, J. Math. Phys. **37**, 2346 (1996), arXiv:gr-qc/9601021.

[5] M. S. Volkov and D. V. Galtsov, Pis'ma Zh. Eksp. Teor. Fiz. **50**, 312 (1989) [JETP Lett. 50, 346 (1989)].

[6] M. S. Volkov, (2016), arXiv:1601.08230 [gr-qc].

[7] S. Droz, M. Heusler, and N. Straumann, Phys. Lett. B **268**, 371 (1991).

[8] P. Breitenlohner, P. Forgács, and D. Maison, Nucl. Phys. B **383**, 357 (1992).

[9] B. R. Greene, S. D. Mathur, and C. M. O'Neill, Phys. Rev. D **47**, 2242 (1993), arXiv:hep-th/9211007.

[10] P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Phys. Rev. D **54**, 5049 (1996), arXiv:hep-th/9511071.

[11] B. Kleihaus and J. Kunz, Phys. Rev. Lett. **79**, 1595 (1997).

[12] M. S. Volkov and D. V. Gal'tsov, Phys. Rept. **319**, 1 (1999), arXiv:hep-th/9810070 [hep-th].

[13] T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. D **90**, 124063 (2014).

[14] C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. **112**, 221101 (2014), arXiv:1403.2757 [gr-qc].

[15] R. Brito, V. Cardoso, and P. Pani, Phys. Rev. D **88**, 064006 (2013), arXiv:1309.0818 [gr-qc].

[16] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D **24**, 1542014 (2015), arXiv:1504.08209 [gr-qc].

[17] Y. Bai and M. Korwar, JHEP **04**, 119 (2021), arXiv:2012.15430 [hep-ph].

[18] K.-M. Lee and E. J. Weinberg, Phys. Rev. Lett. **73**, 1203 (1994), arXiv:hep-th/9406021.

[19] S. A. Ridgway and E. J. Weinberg, Phys. Rev. D **52**, 3440 (1995), arXiv:gr-qc/9503035.

[20] J. Maldacena, JHEP **04**, 079 (2021), arXiv:2004.06084 [hep-th].

[21] J. Ambjorn and P. Olesen, Nucl. Phys. B **315**, 606 (1989).

[22] J. Ambjorn and P. Olesen, Int. J. Mod. Phys. A **5**, 4525 (1990).

[23] M. N. Chernodub, J. Van Doorsselaere, and H. Verschelde, Phys. Rev. D **88**, 065006 (2013), arXiv:1203.5963 [hep-ph].

[24] M. N. Chernodub, V. A. Goy, and A. V. Molochkov, Phys. Rev. Lett. **130**, 111802 (2023), arXiv:2206.14008 [hep-lat].

[25] Y. Bai, J. Berger, M. Korwar, and N. Orlofsky, JHEP **10**, 210 (2020), arXiv:2007.03703 [hep-ph].

[26] Y. Bai, S. Lu, and N. Orlofsky, Phys. Rev. Lett. **127**, 101801 (2021), arXiv:2103.06286 [hep-ph].

[27] D. Ghosh, A. Thalapillil, and F. Ullah, Phys. Rev. D **103**, 023006 (2021), arXiv:2009.03363 [hep-ph].

[28] J. Estes, M. Kavic, S. L. Liebling, M. Lippert, and J. H. Simonetti, JCAP **06**, 017 (2023), arXiv:2209.06060 [astro-ph.HE].

[29] M. D. Diamond and D. E. Kaplan, JHEP **03**, 157 (2022), arXiv:2103.01850 [hep-ph].

[30] Supplemental Material attached after the main text.

[31] R. Gervalle and M. S. Volkov, Phys. Rev. D (to be published).

[32] Y. Nambu, Nucl. Phys. B **130**, 505 (1977).

[33] R. Gervalle and M. S. Volkov, Nucl. Phys. B **984**, 115937 (2022), arXiv:2203.16590 [hep-th].

[34] E. Thuneberg, Phys. Rev. Lett. **56**, 359 (1986).

[35] M. Salomaa and G. Volovik, Phys. Rev. Lett. **56**, 363 (1986).

[36] M. M. Salomaa and G. E. Volovik, Rev. Mod. Phys. **59**, 533 (1987), [Erratum: Rev. Mod. Phys. 60, 573–573 (1988)].

[37] F. Hecht, J. Numer. Math. **20**, 251 (2012).

[38] R. P. Geroch, J. Math. Phys. **11**, 2580 (1970).

[39] R. Hansen, J. Math. Phys. **15**, 46 (1974).

[40] G. Fodor, C. Hoenselaers, and Z. Perjes, J. Math. Phys. **30**, 2252 (1989).

[41] G. Fodor, E. d. S. C. Filho, and B. Hartmann, Phys. Rev. D **104**, 064012 (2021), arXiv:2012.05548 [gr-qc].

[42] N. Arkani-Hamed, L. Motl, A. Nicolis, and C. Vafa, JHEP **06**, 060 (2007), arXiv:hep-th/0601001.

[43] Y. M. Cho and D. Maison, Phys. Lett. B **391**, 360 (1997), arXiv:hep-th/9601028.

[44] R. Gervalle and M. S. Volkov, Nucl. Phys. B **987**, 116112 (2023), arXiv:2211.04875 [hep-th].

[45] C. A. R. Herdeiro and E. Radu, Phys. Rev. D **107**, 064044 (2023), arXiv:2302.00016 [gr-qc].