<figure><img src="image_22.png" /><figcaption>Figure 21. The extremal solutions in phase I contain a small charged black hole inside a bubble of the symmetric phase (grey online), surrounded by a ring-shaped electroweak condensate supporting 22% of the total magnetic charge (green) and two opposite superconducting currents (blue and red). Farther away the condensate disappears and the magnetic field becomes coulombian.</figcaption></figure>

The smallness of the ratio $M_h/|Q_h|$ can be interpreted as a manifestation of the weak gravity conjecture [69]: gravity is the weakest force because the magnetic repulsion between the condensate and the black hole is stronger than the gravitational attraction. However, the condensate cannot be pushed away to infinity due to the Yukawa potential barrier at large distances, as the condensate consists of massive fields.

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

We use in our numerics values of $\kappa$ ranging between $10^{-10}$ and $10^{-2}$ since otherwise, it is challenging to achieve numerical convergence. Extrapolating our results to the physical value $\kappa \sim 10^{-33}$ then yields the following picture of extremal solutions with $|Q| \ll Q_*$.

Their horizon size $r_{\text{ex}} \approx 0.75 \sqrt{\kappa} |n|$ is parametrically small compared to the size of the region where the massive fields change, $r_b \approx 0.73 \sqrt{|n|} < r < r_c \approx 1.13 \sqrt{|n|}$. Therefore, the tiny black hole does not gravitationally affect the electroweak configuration but only supports the strong hypermagnetic field B, which is necessary to create the electroweak condensate. It also regularizes the energy of this field to a finite value, $M_H \approx g|Q| \sim \sqrt{\kappa}, |n|$. The massive condensate resides far from the horizon, where the geometry is nearly flat. The condensate carries mass $M_h \sim \kappa, |n|^{3/2}$ (as explained below after Eq. (10.34)), meaning the ratio $M_h/\kappa$ remains finite as $\kappa \to 0$. This leads to the following consequence.

The dimensionful ADM mass in units of the electroweak mass scale is expressed by (2.24), (2.25),

$$ M = \frac{e^2}{4\pi\alpha} \mathcal{M} \times \mathbf{m}_0, \qquad (10.29) $$

where

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

Therefore,

$$ M = M_H + M_h, \qquad (10.31) $$

where the horizon contribution is

$$ M_H = \frac{e^2}{4\pi\alpha} \mathcal{M}_H \times \mathbf{m}_0 = \frac{eg|n|}{\sqrt{2\kappa}\alpha} \mathbf{m}_0 = 5.1 |n| \mathbf{M}_{\text{Pl}}. \qquad (10.32) $$
Notice that this diverges when $\kappa \to 0$. At the same time, the hair contribution,

$$ M_h = \frac{e^2}{4\pi\alpha} M_h \times m_0, \qquad (10.33) $$

remains finite when $\kappa \to 0$ because $M_h$ reduces in this limit precisely to the regularized mass of the flat space electroweak monopole: the Cho-Maison monopole for $|n| = 2$ [70] and its generalizations for $|n| > 2$ [59]. We have already discussed in Section V C the correspondence between the Cho-Maison monopole and the extremal hairy black hole with $|n| = 2$. It turns out that a similar correspondence applies to higher charge solutions as well.

As a result, the extremal hairy solution with $|Q| \ll Q_*$ can be viewed as the flat space monopole with the same charge $Q$ harboring in the center a tiny charged black hole that almost does not affect the electroweak fields. The total mass of the flat space monopole is infinite because of the Coulombian singularity of the hypermagnetic field B at the origin, but gravity renders the mass finite due to the cutoff at the horizon. Therefore, the extremal hairy black holes provide “a classical UV completion” for the electroweak monopoles.

## 8. Increasing the charge

The size of the black hole is small compared to the size of the surrounding condensate only if the charge is not too large: $|Q| \approx 0.85\sqrt{\kappa}|n| \ll Q_* \approx 0.72/\sqrt{\kappa}$. Increasing $|n|$, the size of the vacuum bubble $r_b \approx 0.73\sqrt{|n|}$ increases and the condensate size $r_c \approx 1.13\sqrt{|n|}$ grows as well, but the horizon radius $r_{ex} \approx 0.75\sqrt{\kappa}|n|$ grows faster, and when $n$ becomes of the order of $1/\kappa$, the black hole starts absorbing the false vacuum region where the Higgs field is very small.

The horizon mass and hair mass are

$$ M_H \approx r_{ex} \approx 0.75\sqrt{\kappa}|n|, \qquad M_h \approx 0.10\kappa|n|^{3/2}, \qquad (10.34) $$

where the latter follows from the estimate obtained in Minkowski space, because the rescaled hair mass $M_h = (8\pi/\kappa)M_h$ essentially coincides with the regularized energy $E_{reg}$ of the flat space monopole given by Eq. (5.24) in [59]. Taking gravity into account, our numerics confirm that $M_h \propto |n|^{3/2}$ if $|Q| < Q_*$.

As a result, although for small $|n|$ the horizon mass $M_H$ is totally dominant as compared to the hair mass $M_h$, the latter increases faster with growing $|n|$ and eventually overtakes the horizon mass. This means that the ratio $M_h/|Q_h|$ increases, but according to our numerics it always remains less than unity. The estimate (10.28) is no longer valid when $Q$ approaches $Q_*$, but the inequality $M < |Q|$ remains valid for all extremal solutions.

On the other hand, the horizon value of the hypermagnetic field, $|B(r_h)| \sim 1/|n|$ (see (10.2)), decreases with growing $|n|$, and when $|Q| \to Q_*$ then $|B(r_h)| \to m_H^2$ and the Higgs field starts deviating from zero at the horizon. The solution then enters the next phase in which the electroweak symmetry is broken everywhere. Analyzing the condition $|B(r_h)| = |n|/(2r_h^2) = m_H^2$ together with Eq. (3.24) determining the value of $r_h = r_{ex}$ yields

$$ Q_* = \frac{2\sqrt{2}g'}{g(g'^2 + \beta)} \frac{1}{\sqrt{\kappa}} = \frac{0.7203}{\sqrt{\kappa}} = 0.6157 Q_m, \qquad (10.35) $$

which reproduces (10.1) and agrees with our numerics.
## B. Phase II – oblate horizon, $|Q| > Q_*$

If $|Q| < Q_*$ then the horizon is spherical and one has $\delta = 0$, but for $|Q| > Q_*$ one has $\delta > 0$, the horizon geometry is no longer spherical and is not described by the extremal RNdS solution (3.19),(3.25). Therefore, the horizon symmetry changes through the transition point $|Q| = Q_*$.

The horizon remains degenerate, and we continue to use the radial gauge defined by Eq.(8.2). However, since the horizon is no longer spherical, there is no longer any reason to enforce similarity with the RNdS case by imposing $r_H = r_{ex}$ in Eq.(8.2). Moreover, the expression for $r_{ex}$ in Eq.(3.24) becomes complex-valued when $|Q| > Q_m$. For these reasons, we fix $r_H = 0.2$ for all solutions in phase II.

Close to the transition point for $|Q| > Q_*$ one has

$$ \delta \propto (|Q| - Q_*)^s. \qquad (10.36) $$

This behavior resembles a second-order phase transition or possibly a crossover, given the rather large value of the critical exponent: according to our numerics, one finds $s \approx 10.8$ for $\kappa = 10^{-2}$.

Additionally, the fraction of the magnetic charge contained in the hair, which was constant in the previous phase, $\lambda = P_h/P \approx 0.22$, starts to decrease for $|Q| > Q_*$ (see Fig.22, left panel). As the charge increases, the black hole gradually absorbs the condensate and becomes less hairy; the ratio $M/|Q|$ increases, the geometry approaches that of the extremal RN solution, and ultimately merges with it when the horizon size overtakes the condensate size. This occurs at $r_h = |Q| = r_c \equiv r_h^0(n) = \sqrt{|n|}/g$, corresponding to the maximal charge

$$ Q_{\max} = 2g'\sqrt{\beta} Q_m = 1.29 Q_m, \qquad (10.37) $$

which in turn corresponds to the maximal value $|n| = 8g'^2/\kappa = 1.78/\kappa$. This estimate of the maximal charge is supported by our numerics. No hairy solutions exist beyond this charge; black holes with $|Q| > Q_{\max}$ are thus RN.

The estimate (10.28) should be reformulated as an inequality valid for all extremal hairy solutions,

$$ 0.88 < M/|Q| \le 1. \qquad (10.38) $$

Here the lower bound corresponds to the minimal charge when $|n| = 2$, whereas the upper bound is reached when the charge is maximal, $|Q| = Q_{\max}$.

## C. Maximally hairy black holes

The black hole is maximally hairy for $|Q| \approx Q_* \approx 0.62 Q_m$ when the ratio $M_h/M$ is maximal (see Fig.22, the insertion in the right panel). This corresponds to the following values of the black hole charge, size and mass:

$$ |n| \approx \frac{0.85}{\kappa} \approx 1.6 \times 10^{32}, \quad r_h \approx 1.41 \text{ cm}, \quad M \approx 2 \times 10^{25} \text{ kg}, \quad \frac{M_h}{M} \approx 0.11. \qquad (10.39) $$

These size and mass are typical for the planetary mass black holes. The hair account for $\approx 11\%$ of the total mass, which means the black hole is surrounded by $\approx 2 \times 10^{24}$ kg of
electroweak condensate (the Earth mass is $M_⊕ ≈ 6 × 10^{24}$ kg). These estimates for the mass and radius are obtained for the physical value of $κ ∼ 10^{-33}$ by noting that for maximally hairy black holes one has $M ≈ |Q|$ while $r_h = r_{ex}$ is given by (3.23). The value ≈ 11% is obtained numerically for $10^{-3} < κ < 10^{-2}$, but it seems to be stable when varying $κ$, therefore we expect it to apply for the physical value of $κ$ as well.

It is also interesting to estimate the condensate size. The condensate terminates where $|B| ≈ m_w^2$ at $r_c ≈ 1.13\sqrt{|n|}$, which corresponds to the dimensionful value

$$r_c ≈ 2.20 \text{ cm.} \qquad (10.40)$$

One may wonder how massive fields can extend over macroscopic distances; the key point is that they form a Bose-Einstein condensate. A familiar example is the condensate in a superconductor, which also can be macroscopically large despite involving massive fields. When the electroweak condensate terminates, the corona does not vanish abruptly but instead decays as $\sim Q_M/r^4$, due to the presence of a magnetic quadrupole moment $Q_M$ carried by the magnetic field.

One might think that the magnetic charge with $|n| = 1.6 × 10^{32}$ is extremely large. However, as explained in Appendix D, a charge of this magnitude placed at the center of the Sun would create a magnetic field of $≈ 10^3$ Gauss at the Sun surface. This is 1000 times larger than the average proper magnetic field of the Sun, but 2-4 times smaller than the field in sunspots.

## D. Summary

Summarizing, the extremal hairy solutions show two phases depending on value of their charge $Q ≈ 0.85\sqrt{κ}n$. If the charge is small then the hypermagnetic field B is strong enough to restore the full electroweak symmetry at the horizon where $W_{\muν}^a$ and $Φ$ vanish. The horizon geometry is then perfectly spherical, but the magnetic fields $B$ and $B_Z$ keep a non-trivial dependence on the polar angle $ϑ$ and form a corona of vortex loops attached to the horizon. These vortices do not change the magnetic flux $4πP$ through the horizon but render it inhomogeneous, highly increasing the field strength in the narrow equatorial region. Although all the charge $P$ is contained inside the black hole in the 't Hooft description, the total flux can be interpreted as a superposition of two contributions: the flux produced by an effective charge $g'^2P$ distributed along a ring at the horizon in the equatorial plane, and the remaining flux $g^2P$ distributed more or less homogeneously over the entire horizon.

When the charge increases, the horizon value of B decreases. If the charge exceeds the critical value $Q_* ≈ 0.62 Q_M$ where $Q_M ≈ 1.17/\sqrt{κ}$, then B becomes insufficient to create a region of symmetric phase around the horizon. As a result, the Higgs field starts deviating from zero at the horizon. The horizon geometry is then no longer spherical. However, the B-field is still strong enough to create a condensate around the horizon. The condensate disappears and the black hole loses its hair when the charge approaches the maximal value $Q_{max} ≈ 1.29 Q_M$.

The mass of all extremal hairy solutions is smaller than the charge, $M/|Q| ≤ 1$, which can be related to the Zeeman interaction between the W bosons in the condensate and the black hole magnetic field.

Various parameters of the extremal solutions against the charge are shown in Fig.22. One can see that the quadrupole moments and 't Hooft currents attain their maximal values in
<figure><img src="image_23.png" /><figcaption>Figure 22. The parameters of extremal solutions (left) and the existence diagram for hairy solutions (right) for $\kappa = 10^{-2}$. The upper violet curve on the right corresponds to solutions bifurcating with the RN while the lower green curve corresponds to extremal solutions. The blue shaded area between these two curves is the domain of existence of non-extremal hairy solutions. The scale for the charge is set by $Q_m \approx 1.17/\sqrt{\kappa}$ defined in (3.24).</figcaption></figure>

phase I. They vanish for the minimal value of the charge $|n| = 2$ when the configuration is spherically symmetric, and for the maximal charge when the hair disappears. The horizon oblateness and the Nambu currents attain maximal values in phase II, while the fraction of the hair charge $\lambda = P_h/P$ is constant in phase I and decreases down to zero in phase II. As seen in Fig.22, one has $r_h/|Q| \approx M/|Q| \le 1$.

The domain of existence for these hairy solutions is shown in the right panel of Fig.22, which presents the ratio $M/|Q|$ against the magnetic charge $|Q| = \sqrt{\kappa}/2|n|/(2e)$. The hairy solutions occupy the shaded region between two curves: below the upper curve corresponding to the bifurcation with the RN branch, and above the lower curve corresponding to the extremal hairy solutions. Spherically symmetric solutions with $|n| = 2$ correspond to the left edge of the region, they were previously known [43], whereas all solutions with $|n| > 2$ are new. Descending from the upper curve to the lower one at fixed $Q$ increases the amount of charge stored in the condensate and decreases the ratio $M/|Q|$. As shown in the figure, the inequality $M/|Q| \le 1$ holds not only for extremal solutions but also for near-extremal ones. For the extremal solutions in Fig.22, one finds $0.91 \le M/|Q| \le 1$, which slightly differs from the estimate in Eq.(10.38) because the figure is obtained not for the physical value of the gravitational coupling $\kappa \sim 10^{-33}$, but for $\kappa = 10^{-2}$.

XI. CONCLUDING REMARKS

We have constructed solutions describing static and axially symmetric black holes that support a magnetic field strong enough to change the structure of the electroweak vacuum. This results in the formation of an electroweak condensate outside the black hole horizon. The condensate carries currents that flow along the horizon and generate a corona of fountain-like vortex loops around the horizon.

According to Maldacena [46], the presence of the corona should significantly enhance the evaporation rate. As a result, non-extremal black holes are expected to rapidly relax
toward the extremal state, which correspond to the lower edge of the existence region where $M \le |Q|$. In this extremal state, the black holes cannot shed their hair to evolve into RN black holes, for which $M \ge |Q|$. In principle, they could decay by emitting pairs of hairy black holes with $n = \pm 2$, which minimize the ratio $M/|Q|$, but such a process would require quantum tunnelling and is therefore expected to be strongly suppressed.

However, axially symmetric black holes can reduce their energy by developing a fully
spread corona of vortices, which may retain at most only discrete symmetries. In the axially
symmetric case, the corona consists of two oppositely directed multi-vortices. Each multi-
vortex is composed of many elementary vortices that repel each other, making the configur-
ation inherently unstable (except for $|n| = 4$). This instability is evident at the perturbative
level, as shown in Section VI, where the axially symmetric configuration corresponds to a
saddle point of the energy. In contrast, the absolute energy minimum corresponds to a
configuration of elementary vortices that maximize their mutual separation, as illustrated
in Fig.6 (left panel). Non-perturbative generalizations of this configuration are expected to
describe stable hairy black holes surrounded by an electroweak corona. However, the cor-
responding solutions are currently unknown, except for the spherically symmetric case with
$n = \pm 2$, which contains no vortices, and the axially symmetric case with $n = \pm 4$, where the
vortices are already elementary.

The transition from the axially symmetric state to a fully spread corona should not only
reduce the energy, but also increase the entropy. One may expect that, for very large $|n|$,
the condensate and the associated vortices in the corona form a gas-like or liquid-like phase
[50]. The ground state of such a system is likely to be highly degenerate, with a high density
of low-lying excitations above it. As a result, there should be a nontrivial entropy associated
with the black hole hair.

As we have seen above, the black hole mass and charge split into two parts, associated
respectively with the black hole itself and with the surrounding hair. It is therefore natural
to assume that the entropy also splits as $S = S_H + S_h$, where the first term is determined
by the horizon area. As shown in Fig. 22, for extremal hairy black holes one finds $r_h < |Q|$,
whereas for extremal RN black holes $r_h = |Q|$, so the horizon entropy satisfies $S_H < S_{RN}$.

On the other hand, as shown in Fig.22, the ratio $r_h/|Q|$ is close to unity. It is therefore
possible that including the hair entropy may reverse the inequality to $S_H + S_h > S_{RN}$.
This would imply that extremal hairy black holes with a fully spread corona are not only
energetically preferred over RN black holes, but are also favoured entropically, which would
enhance their formation rate. However, for this to be true, the hair entropy $S_h$ should
be sufficiently large. It is expected to be maximal for maximally hairy black holes, where
$r_h \sim Q_* \sim 1/\sqrt{\kappa} \sim \sqrt{|n|}$, implying that $S_H \sim S_{RN} \sim r_h^2/\kappa \sim |n|/\kappa \sim n^2$. Therefore, $S_h$
should also scale as $n^2$. Using, for instance, the ideal gas approximation for the hair entropy
gives $S_h \sim |n|\ln|n|$, which is not large enough. On the other hand, the vortices in the
corona are not point particles but extended objects (fountains) that can become mutually
entangled in a complex way, as we observed in the perturbative description of Section VI.
This should increase their entropy. Since we cannot currently give a reliable estimate for
$S_h$, we leave this issue for future analysis.

Our solutions provide the first example of hairy black holes described by well-tested
physical theories – the Standard Model and General Relativity – and may therefore be
physically relevant (the $n = 2$ solutions of [43] are too close to the Planck scale to be
considered realistic). Such black holes could plausibly originate from primordial black holes
acquiring charge from the early electroweak plasma.
Suppose that the fluctuating hypermagnetic field in the ambient plasma becomes, at some moment, predominantly orthogonal to the black hole horizon. This would induce a net flux through the horizon, and hence generate a charge. Since the total magnetic charge is conserved, the flux generated on one black hole must be compensated by an opposite flux on another black hole (or several others). The resulting oppositely charged black holes would attract one another, but they might not necessarily annihilate, as cosmic expansion can drive them apart. It is also possible that oppositely charged black holes form bound systems stabilized by scalar repulsion, as suggested in [88].

Such a mechanism for magnetic charge generation remains to be verified, but if confirmed, it would have significant implications. Magnetic monopoles are typically associated with Grand Unified Theories (GUT), which involve physics beyond the Standard Model – the existence of which remains an open question. However, if primordial black holes can indeed acquire magnetic charge through a process similar to the one described above, this would suggest that no new physics beyond the Standard Model is required to account for monopoles. In this scenario, the magnetic charge does not require a conventional carrier particle but is instead associated with the magnetic flux threading the black hole event horizon. From this perspective, extremal hairy electroweak black holes naturally emerge as compelling candidates for magnetic monopoles – objects that may, in principle, exist in Nature without invoking physics beyond the Standard Model.

These black holes can exhibit a variety of physical effects. For instance, they are expected to catalyze proton decay [46]. They could potentially be detected if captured by a neutron star, leading to a sudden change in the star’s rotation period [54]. Their phenomenology is quite rich and has been explored in Refs. [51–55].

Finally, it is worth noting that there are numerous directions in which our results can be extended and generalized. Most obviously, one should explicitly construct non-axially symmetric solutions with $|n| > 4$ and demonstrate that they indeed minimize the energy. It would also be natural to explore their stationary generalizations. Crucially, one should clarify the physical mechanism by which black holes can acquire magnetic charge. In addition, it is important to account for the role of fermions, verify the enhancement of Hawking evaporation, examine the Rubakov-Callan effect, further develop the phenomenology of these objects, etc.

**Acknowledgements.** We thank Eugen Radu, Maxim Chernodub, Julien Garaud, and Christos Charmousis for valuable discussions. M.S.V. is grateful to the Theory Division of CERN for its hospitality during the final stage of this work.

## Appendix A. ENERGY-MOMENTUM TENSOR

In this Appendix, we explicitly present the components of the electroweak energy-momentum tensor in the axially symmetric case. Substituting (7.2), (7.22) into (2.20) yields the following diagonal components of the energy-momentum tensor:

$$ 
\begin{aligned}
-T_0^0 \sqrt{-g} &= a_1 + a_2 + a_3 + a_4 + a_5, \\
-T_r^r \sqrt{-g} &= -a_1 + a_2 - a_3 + a_4 + a_5, \\
-T_\vartheta^\vartheta \sqrt{-g} &= a_1 - a_2 - a_3 + a_4 + a_5, \\
-T_\phi^\phi \sqrt{-g} &= -a_1 - a_2 + a_3 + a_4 - a_5,
\end{aligned} 
\quad (A.1) $$
where $\sqrt{-g} = e^{2K+S-2U}r^2 \sin \vartheta$, and

$$
\begin{align*}
a_1 &= \mathrm{Ne}^{2\mathrm{U}-\mathrm{S}} \left[ \frac{1}{2g'^2} (\partial_r Y)^2 + \frac{1}{2g^2} ((\mathcal{D}_r F_3)^2 + (\mathcal{D}_r F_4)^2) \right] \frac{\nu^2}{\sin \vartheta}, \\
a_2 &= \mathrm{Ne}^{2\mathrm{U}-\mathrm{S}} \left[ \frac{1}{2g'^2} (\partial_\vartheta Y)^2 + \frac{1}{2g^2} ((\mathcal{D}_\vartheta F_3)^2 + (\mathcal{D}_\vartheta F_4)^2) \right] \frac{\nu^2}{r^2 \sin \vartheta}, \\
a_3 &= \mathrm{Ne}^S \left[ \frac{1}{2g^2} e^{2\mathrm{U}-2\mathrm{K}} (\partial_\vartheta F_1 - \partial_r F_2)^2 + r^2 (\mathcal{D}_r \phi_1)^2 + r^2 (\mathcal{D}_r \phi_2)^2 \right] \sin \vartheta, \\
a_4 &= \mathrm{Ne}^S \left[ (\mathcal{D}_\vartheta \phi_1)^2 + (\mathcal{D}_\vartheta \phi_2)^2 + \frac{\beta r^2}{8} e^{2\mathrm{K}-2\mathrm{U}} (\phi_1^2 + \phi_2^2 - 1)^2 \right] \sin \vartheta, \\
a_5 &= \mathrm{Ne}^{2\mathrm{K}-\mathrm{S}} \left[ ((F_3 + Y)\phi_1 - F_4\phi_2)^2 + ((F_3 - Y)\phi_2 + F_4\phi_1)^2 \right] \frac{\nu^2}{4 \sin \vartheta}. \tag{A.2}
\end{align*}
$$

Here, the following abbreviations have been introduced:

$$
\begin{align*}
\mathcal{D}_r F_3 &= \partial_r F_3 + F_1 F_4, & \mathcal{D}_\vartheta F_3 &= \partial_\vartheta F_3 + F_2 F_4, & \mathcal{D}_r F_4 &= \partial_r F_4 - F_1 F_3, & \mathcal{D}_\vartheta F_4 &= \partial_\vartheta F_4 - F_2 F_3, \\
\mathcal{D}_r \phi_1 &= \partial_r \phi_1 - \frac{F_1}{2} \phi_2, & \mathcal{D}_\vartheta \phi_1 &= \partial_\vartheta \phi_1 - \frac{F_2}{2} \phi_2, & \mathcal{D}_r \phi_2 &= \partial_r \phi_2 + \frac{F_1}{2} \phi_1, & \mathcal{D}_\vartheta \phi_2 &= \partial_\vartheta \phi_2 + \frac{F_2}{2} \phi_1.
\end{align*}
$$

The only off-diagonal components are $T_{r\vartheta} = T_{\vartheta r}$ with

$$
\begin{equation}
\begin{split}
T_{\vartheta}^{r} \sqrt{-g} = \mathrm{Ne}^{2\mathrm{U}-\mathrm{S}} & \left[ \frac{1}{g'^2} \partial_r Y \partial_\vartheta Y + \frac{1}{g^2} (\mathcal{D}_r F_3 \mathcal{D}_\vartheta F_3 + \mathcal{D}_r F_4 \mathcal{D}_\vartheta F_4) \right] \frac{\nu^2}{\sin \vartheta} \\
& + 2\mathrm{r}^2 \mathrm{Ne}^S (\mathcal{D}_r \phi_1 \mathcal{D}_\vartheta \phi_1 + \mathcal{D}_r \phi_2 \mathcal{D}_\vartheta \phi_2) \sin \vartheta.
\end{split}
\tag{A.3}
\end{equation}
$$

The energy-momentum tensor is gauge-invariant. Its diagonal components are invariant under the reflection $\vartheta \to \pi - \vartheta$, while the off-diagonal components then change sign.

## Appendix B. AXIALLY-SYMMETRIC FORM OF THE $n = \pm 2$ SOLUTIONS

Before considering axially symmetric solutions with $|n| > 2$, it is instructive to first reproduce the spherically symmetric solutions with $|n| = 2$ within the axially symmetric formalism. A spherically symmetric solution is described by Eqs.(3.1),(3.2), with functions $N(r)$, $\sigma(r)$, $f(r)$, $\phi(r)$ depending on the Schwarzschild radial coordinate $r$.

The same solution, expressed within the axially-symmetric formalism, is described by the line element (7.2) and electroweak fields (7.22),(7.27) parameterized by functions N, U, K, S, $H_k$, y, $\phi_1$, $\phi_2$. These functions depend on the radial coordinate $r$, distinct from $r$.

Notice that $N(r)$ is obtained by solving the equations, whereas $N(r)$ is the auxiliary function whose choice determines the radial gauge. For example, one can choose $N(r) = 1$.

The relation between the coordinates $r$ and $r$ is determined by

$$ \frac{dr}{r\sqrt{N(r)}} = \frac{dr}{r\sqrt{N(r)}}, \qquad (B.1) $$

integrating which yields $r = r(r)$. One has then

$$
\begin{align}
e^{2\mathrm{U}} &= \frac{N(r)}{\mathrm{N}(r)} \sigma^2(r), & e^{\mathrm{K}-\mathrm{U}} &= e^{\mathrm{S}-\mathrm{U}} = \frac{r}{r}, \nonumber \\
H_2 = H_4 &= f(r), & H_1 = H_3 &= y = \phi_1 = 0, & \phi_2 = \phi(r). \tag{B.2}
\end{align}
$$
A. Non-extremal solutions

For non-extremal solutions the metric coefficient g$_{00}$ shows a simple zero at the horizon, which can be ensured by choosing the auxiliary function N(r) to be

$$N = 1 - \frac{r_H}{r}.$$

For the non-extremal RN solution (3.6), Eqs.(B.1),(B.2) then yield

$$\begin{gathered} r = r - r_{-}, \quad r_H = r_+ - r_{-}, \quad e^U = \frac{r}{r + r_{-}}, \\ K = S = H_1 = H_2 = H_3 = H_4 = y = \phi_1 = 0, \quad \phi_2 = 1, \end{gathered} \qquad (B.3)$$

where $r_{\pm} = M \pm \sqrt{M^2 - Q^2}$ are the two RN horizons. One has

$$r_+ r_- = Q^2 \equiv \frac{\kappa}{2e^2} = r_-^2 + r_H r_- \quad \Rightarrow \quad r_- = \frac{1}{2} \left( \sqrt{r_H^2 + 4Q^2} - r_H \right). \qquad (B.4)$$

In the extremal limit, when $M \to Q$, one has $r_H \to 0$, while the horizon value $U_H = \ln(r_H/r_+)$ diverges.

<figure><img src="image_24.png" /><figcaption>Figure 23. The metric amplitude U and the difference S – U for the non-extremal spherically symmetric hairy black hole with $n = 2$, $\kappa = 10^{-3}$, and $r_H = 0.02$ as functions of $\ln(r/r_H)$ and $\vartheta$.</figcaption></figure>

For the non-extremal hairy solutions one can represent $N(r)$ in (3.1) as

$$N(r) = \frac{1}{\alpha^2(r)} \left(1 - \frac{r_h}{r}\right), \qquad (B.5)$$

where $\alpha(r_h) \equiv \alpha_h$ and $\alpha(\infty) = 1$. Integrating then Eq.(B.1) yields the relation between r and r, expressed by

$$X(r, r_H) = (X(r, r_h))^{\alpha(r)} \times \mathcal{F}(r) \quad \text{with} \quad \mathcal{F}(r) = \exp \left( \int_r^\infty \alpha'(r) \ln (X(r, r_h)) dr \right), \quad (B.6)$$
where the following definition has been used:

$$X(x, x_1) = \frac{\sqrt{x} + \sqrt{x - x_1}}{2} \qquad (B.7)$$

This implies that at infinity, where $r, r \to \infty$, one has $r/r \to 1$ and U, S, K $\to 0$. At the horizon, where $r = r_H$ and $r = r_h$, one has

$$\sqrt{r_H} = 2 \left( \frac{\sqrt{r_h}}{2} \right)^{\alpha_h} \times \mathcal{F}(r_h), \qquad (B.8)$$

hence $e^{S_H - U_H} = r_h/r_H \neq 1$. With some algebra, one finds

$$\lim_{r \to r_h} \frac{dN(r)}{dr} = \frac{1}{\alpha_h^2 r_h}, \qquad \lim_{r \to r_h} \frac{dN(r(r))}{dr} = \frac{\alpha_h^2}{r_h}, \qquad \Rightarrow \qquad e^{U_H} = \frac{\sigma_h}{\alpha_h^2}. \qquad (B.9)$$

Applying the surface gravity formula (7.19) separately within the spherically symmetric and axially symmetric formalisms yields

$$k_H = \frac{\sigma_h}{2\alpha_h^2 r_h} = \frac{1}{2r_h} e^{U_H} = \frac{1}{2r_H} e^{2U_H - K_H}. \qquad (B.10)$$

In the extremal limit, $r_h \to r_{ex} > 0$ when the horizon becomes degenerate and $k_H \to 0$, the value $\sigma_h$ approaches a non-zero limit, so $\alpha_h \to \infty$. It follows that

$$r_h \to r_{ex}, \qquad r_H \to 0, \qquad e^{U_H} \to 0, \qquad e^{K_H - U_H} = e^{S_H - U_H} = \frac{r_h}{r_H} \to \infty. \qquad (B.11)$$

Having computed U, S, ... from (B.2), the same values are obtained by solving the axially symmetric equations numerically. Fig.23 shows the profiles of U and S - U for $\kappa = 10^{-3}$ and $r_H = 0.02$. One finds that $r_h = r_H e^{S_H - U_H} = 0.058 > r_{ex} = 0.047$, while $r_H < r_{ex}$. As $r_H$ decreases further and approaches zero, $U_H$ shifts more and more to negative values, while $S_H - U_H$ grows, as suggested by (B.11).

## B. Extremal solutions

Since $U_H$ becomes unbounded when approaching the extremal limit, this case should be considered separately. In this limit, the horizon becomes degenerate, and the function $N(r)$ takes structure shown in (5.18). Equations (B.1) and (B.2) then imply that the amplitudes U, K, S will be finite at the horizon and approach zero at infinity, provided that $N(r)$ has the same structure as $N(r)$:

$$N(r) = k^2(r) \left(1 - \frac{r_H}{r}\right)^2, \qquad N(r) = k^2(r) \left(1 - \frac{r_{ex}}{r}\right)^2, \qquad (B.12)$$

where $r_{ex}$ is defined in (3.23) with $k(r_H) = k(r_{ex}) \equiv k$, while $k(\infty) = k(\infty) = 1$.

It is convenient to set

$$r_H = r_{ex}, \qquad (B.13)$$
so that the radial coordinate r coincides with the Schwarzschild coordinate r at the horizon. Integrating Eq. (B.1) then yields

$$ \frac{\ln(r - r_H)}{k(r)} - \frac{\ln(r - r_H)}{k(r)} = \int_{r}^{\infty} \ln(r - r_H) \frac{k'(r)}{k^2(r)} dr - \int_{r}^{\infty} \ln(r - r_H) \frac{k'(r)}{k^2(r)} dr. \quad (B.14) $$

The right-hand side of this relation approaches zero when $r, r \to \infty$, so that $r/r \to 1$. According to (B.2), this implies that $U, K, S \to 0$. When $r, r \to r_H$, the right-hand side of (B.14) approaches a finite value, which we denote by $C$. Hence, close to the horizon, we have

$$ \frac{r - r_H}{r - r_H} = \frac{\sqrt{N}}{\sqrt{N}} = e^{kC}. \quad (B.15) $$

Using Eq. (B.2), this implies

$$ e^{U_H} = e^{S_H} = e^{K_H} = \sigma_h e^{kC}, \quad (B.16) $$

so that the metric amplitudes are now finite at the horizon.

Applying the above formulas to the extremal RN solution yields a very simple result:

$$ N = \left(1 - \frac{r_H}{r}\right)^2, \quad U = K = S = H_1 = H_2 = H_3 = H_4 = y = \phi_1 = 0, \quad \phi_2 = 1. \quad (B.17) $$

For the extremal hairy black hole, the function $k(r)$ takes the form (3.25) at the horizon. Hence, by introducing an additional cutoff factor, we set

$$ 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_H^4}, \quad (B.18) $$

so that $k^2(r_H) = k^2(r_{ex})$ and $k^2(\infty) = 1$. Integrating the axially symmetric equations

<figure><img src="image_25.png" /><figcaption>Figure 24. The metric amplitude U and the difference S − U for the extremal spherically symmetric hairy black hole with $n = 2$ and $\kappa = 10^{-3}$, in which case $r_H = r_{ex} = 0.047$.</figcaption></figure>

numerically then yields profiles of U and S − U, shown in Fig.24. These profiles are quite different from those for the non-extremal solution in Fig.23, but they agree with the result in (B.16). In particular, $S − U = \ln(r/r)$ now vanishes both at the horizon and at infinity.
### C. Fixing the radial gauge – summary

The upshot is as follows. The spherically symmetric hairy solutions can be labelled by the Schwarzschild radius of their event horizon, $r_h \in [r_{ex}, r_h^0]$. The same solutions within the axially symmetric description are labeled by $r_H$, which we call event horizon size, even though it does not have a direct physical meaning. The parameter $r_H$ and the auxiliary function $N(r)$ in (7.2) take the following values:

$$
\begin{align}
\text{Non-extremal case:} \qquad & 0 < r_H \le r_H^0 = r_+^0 - r_-^0, \qquad N(r) = 1 - \frac{r_H}{r}; \nonumber \\
\text{Extremal case:} \qquad & r_H = r_{ex}, \qquad N(r) = k^2(r) \left(1 - \frac{r_H}{r}\right)^2; \tag{B.19}
\end{align}
$$

with $k(r)$ given by (B.18) and $r_{ex}$ defined in (3.23). Here, $r_+^0 \equiv r_h^0$ and $r_-^0$ are the two horizons of the RN solution, to which the hairy black holes reduce at the bifurcation point $r_h = r_h^0$ (note that for the RN solution, $r_H = r_+ - r_-$). Since $r_+ r_- = Q^2$ with $Q$ from (3.6), one has

$$r_H^0 = r_+^0 - r_-^0 = r_h^0 - Q^2/r_h^0. \tag{B.20}$$

We use the same gauge conditions (B.19), (B.20) also for $|n| > 2$.

### D. Choice of $r_H$ and diffeomorphisms

The value of $r_H$ can be changed by a diffeomorphism in the extremal case, while for the non-extremal case this is not possible. Suppose that $r \to r(\tilde{r})$ and consider Eq. (7.3),

$$\frac{dr}{r\sqrt{N(r)}} = \frac{d\tilde{r}}{\tilde{r}\sqrt{\tilde{N}(\tilde{r})}}. \tag{B.21}$$

Integrating this equation yields a relation between $r$ and $\tilde{r}$. In the non-extremal case, one chooses $N = 1 - r_H/r$ and $\tilde{N} = 1 - \tilde{r}_H/\tilde{r}$. Performing the integration then shows that the only possibility is $r = \tilde{r}$ and $r_H = \tilde{r}_H$. Therefore, the value of $r_H$ cannot be changed by a diffeomorphism.

In the extremal case, $N$ is chosen as in Eq. (B.12), and similarly for $\tilde{N}$, with the replacement $r \to \tilde{r}$, $r_H \to \tilde{r}_H$. In this case the integration of (B.21) yields a relation similar to that in (B.14):

$$\frac{\ln(r - r_H)}{k(r)} + \int_{r_H}^{r} \ln(r - r_H) \frac{k'(r)}{k^2(r)} dr = \frac{\ln(\tilde{r} - \tilde{r}_H)}{\tilde{k}(\tilde{r})} + \int_{\tilde{r}_H}^{\tilde{r}} \ln(\tilde{r} - \tilde{r}_H) \frac{\tilde{k}'(\tilde{r})}{\tilde{k}^2(\tilde{r})} d\tilde{r}. \tag{B.22}$$

This establishes a non-trivial correspondance between $r$ and $\tilde{r}$ for any given $r_H$ and $\tilde{r}_H$. In the vicinity of the horizon, it simplifies to

$$\tilde{r} - \tilde{r}_H = (r - r_H)^{\tilde{k}/k}(1 + \dots), \tag{B.23}$$

where $k \equiv k(r_H)$, $\tilde{k} \equiv \tilde{k}(\tilde{r}_H)$, and the dots denote subleading term.

Assuming that $r_H = r_{ex}$, which is the radius of the RNdS solution defined in Eq. (3.24), the metric amplitudes U, K, S remain everywhere regular. However, since $k \neq \tilde{k}$, it follows from
the relation $e^U\sqrt{N} = e^{\tilde{U}}\sqrt{\tilde{N}}$ (see Eq. (7.3)) that the amplitude $\tilde{U}$ becomes unbounded at the horizon. This implies that, while the value of $r_H$ can indeed be changed by a diffeomorphism in the extremal case, the resulting gauge choice leads to a singularity at the horizon.

At the same time, the singularity is very mild and located in a narrow region near the horizon. We have verified that choosing in our numerical code $r_H = r_{ex}$, in which case everything is regular, or selecting a different value such as $r_H = 1$, in which case there is a singularity, returns the same values of the physical parameters such as the total mass, horizon area, etc.

Therefore, we adopt the following strategy. As described in Section X, extremal hairy solutions exhibit two distinct phases depending on their charge. For extremal solutions in phase I, where the horizon remains spherical, as for the RNdS solution, we always choose $r_H = r_{ex}$ to ensure regularity at the horizon.

In contrast, for extremal solutions in phase II, the horizon is still degenerate but no longer spherical. In this regime, there is no compelling reason to preserve similarity with the RNdS case. Furthermore, hairy solutions in phase II exist for $|Q| > Q_m$, in which case the value of $r_{ex}$ defined by Eq. (3.24) becomes complex. Therefore, we abandon the condition $r_H = r_{ex}$ and instead set $r_H = 0.2$ for all solutions in phase II. Numerical solutions obtained in this gauge exhibit no irregularities, and the virial relation is satisfied. We therefore conclude that this gauge choice is acceptable. Moreover, we observe a smooth transition between phase I and phase II.

## Appendix C. ANTI-SCREENING PROPERTY OF THE CONDENSATE

It is instructive to consider simple models illustrating the anti-screening property of the electroweak condensate. This property is the opposite of the well-known screening within the conventional Abelian model of superconductivity defined by Eq. (F.1) in Appendix F. This model admits vortices consisting of a magnetic field and a complex scalar [89, 90]. Their magnetic field is maximal at the vortex center, where scalar describing the superconducting condensate vanishes. However, the magnetic field decays to zero away from the center, where the condensate approaches a constant value. Therefore, the field is screened by the condensate: it is zero where the condensate is maximal.

As we shall now see, the structure of vortices appearing via the electroweak condensation is just the opposite: their magnetic field is maximal where the condensate is maximal, so that the field is enhanced by the condensate [47, 48].

In the linear approximation, the electroweak condensate is described by a complex vector $w_\mu$ subject to Eq. (4.3),

$$ \mathcal{D}^\mu w_{\mu\nu} + ieF_{\nu\sigma}w^\sigma = m_W^2 w_\nu, \qquad (C.1) $$

with $w_{\mu\nu} = \mathcal{D}_\mu w_\nu - \mathcal{D}_\nu w_\mu$ and $\mathcal{D}_\mu = \nabla_\mu + ieA_\mu$, where $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$ is the background electromagnetic field. In some special cases, this equation can be solved by

$$ w_{\mu\nu} = 0, \qquad ieF_{\nu\sigma}w^\sigma = m_W^2 w_\nu. \qquad (C.2) $$

The condensate generates the current $J_\nu$ and produces a second order correction $f_{\mu\nu}$ to the background electromagnetic field $F_{\mu\nu}$. One has (see (6.5),(6.10),(6.23))

$$ 4\pi J_\nu = \frac{ie}{2} \nabla^\sigma (w_\sigma \bar{w}_\nu - w_\nu \bar{w}_\sigma), \qquad \nabla^\mu f_{\nu\mu} = 4\pi J_\nu. \qquad (C.3) $$
1. Magnetic field along z-axis.

Let us consider a constant magnetic field inside a cylinder of radius $R$ in flat space,

$$ds^2 = -dt^2 + d\rho^2 + \rho^2 d\varphi^2 + dz^2, \quad A_\varphi = \frac{B}{2} \rho^2, \quad \mathcal{B}^z = \frac{1}{\rho} \epsilon^{z\rho\varphi} F_{\rho\varphi} \equiv B > 0, \quad (C.4)$$

where $\rho \in [0, R]$. The equations (C.2) are solved by

$$B = \frac{m_w^2}{e}, \quad w_0 = w_z = 0, \quad w_\rho = \frac{i}{\rho} w_\varphi = \exp(-i\varphi) \exp\left(-\frac{m_w^2}{4}\rho^2\right), \quad (C.5)$$

hence, the magnetic flux in the cylinder is equal to $\pi m_w^2 R^2 / e$ and condensate density is

$$w(\rho) \equiv w^\mu \bar{w}_\mu = 2|w_\rho|^2 = 2 \exp\left(-\frac{m_w^2}{2}\rho^2\right). \quad (C.6)$$

This gives rise to the current with the only non-vanishing component in the azimuthal direction,

$$4\pi J_{\varphi} = -e\rho (|w_{\rho}|^2)'_{\rho} = em_w^2 \rho^2 |w_{\rho}|^2, \quad (C.7)$$

whose sign is opposite to that of the Larmor precession. Therefore, the field $\delta\mathcal{B}$ created by the current rather enhances the background field instead of screening it.

The Maxwell equation (C.3) reduces to

$$\frac{1}{\rho} (\rho f^{\varphi\rho})'_{\rho} = 4\pi J^{\varphi} = -\frac{e}{\rho} (|w_{\rho}|^2)'_{\rho} \quad \rightarrow \quad \delta\mathcal{B}^z = \frac{1}{\rho} f_{\rho\varphi} = e|w_{\rho}|^2 + \text{const.}, \quad (C.8)$$

where the constant is fixed by the requirement that the field $\delta\mathcal{B}_z$ should not change the magnetic flux, hence

$$\int_0^R \rho \delta\mathcal{B}^z d\rho = 0 \quad \Rightarrow \quad \delta\mathcal{B}^z = e \left( \frac{1}{2} w(\rho) - \frac{2 - w(R)}{m_w^2 R^2} \right). \quad (C.9)$$

This field is maximal at $\rho = 0$ where the condensate density $w(\rho)$ is maximal and where $\delta\mathcal{B}^z > 0$, therefore the background field is enhanced by the condensate (anti-screening). At the cylinder surface, for $\rho = R$, the condensate is minimal and one has $\delta\mathcal{B}^z < 0$.

2. Magnetic field on a sphere.

Let us consider a radial magnetic field in a spacetime with the direct product geometry,

$$ds^2 = -dt^2 + dr^2 + R^2(d\vartheta^2 + \sin^2\vartheta d\varphi^2), \quad A_\varphi = -P \cos\vartheta, \quad \mathcal{B}^r = \frac{F_{\vartheta\varphi}}{R^2 \sin\vartheta} = \frac{P}{R^2}, \quad (C.10)$$

where $r \in (-\infty, \infty)$ while $R = \text{const}$. The magnetic flux through the 2-sphere is $4\pi P$, and we assume for notational simplicity that $P > 0$.

The equations (C.2) are solved by

$$R^2 = \frac{eP}{m_w^2}, \quad w_0 = w_r = 0, \quad w_\vartheta = \frac{i}{\sin\vartheta} w_\varphi = (\sin\vartheta)^j, \quad j = eP - 1. \quad (C.11)$$
The condensate density is $w(\vartheta) = w^\mu \bar{w}_\mu = (2/R^2)|w_\vartheta|^2$.

The current has only one component in the azimuthal direction,

$$4\pi J_{\varphi} = -\frac{e}{R^2} \sin\vartheta (|w_{\vartheta}|^2)'_{\vartheta} = -\frac{2ej}{R^2} \cos\vartheta |w_{\vartheta}|^2. \qquad (C.12)$$

The sign in (C.12) agrees with that in (6.21), but in the upper hemisphere, for $0 \le \vartheta \le \pi/2$, it is opposite to that in (C.7) because $|w_\vartheta|^2$ increases with $\vartheta$, whereas $|w_\rho|^2$ in (C.7) decreases with $\rho$. This produces the second order correction for the magnetic field,

$$\delta B^r = \frac{f_{\vartheta\varphi}}{R^2 \sin\vartheta} = \frac{e}{R^2} |w_{\vartheta}|^2 + \text{const.}, \qquad (C.13)$$

where the integration constant is fixed by the requirement that the total magnetic flux should not change, $\oint \delta B_r d^2\Omega = 0$ (with $d^2\Omega = \sin\vartheta d\vartheta d\varphi$), hence

$$\frac{R^2}{e} \delta B^r = (\sin \vartheta)^{2j} - \frac{\sqrt{\pi} \Gamma(j+1)}{2 \Gamma(j+3/2)}. \qquad (C.14)$$

Therefore, the total field is

$$\frac{P}{R^2} + \delta B^r = \frac{m_w^2}{e} + \frac{e}{R^2} \left( (\sin \vartheta)^{2j} - \frac{\sqrt{\pi} \Gamma(j+1)}{2 \Gamma(j+3/2)} \right). \qquad (C.15)$$

For $j = 0$, this reduces to the background field $P/R^2 = m_w^2/e$. For $j > 0$, the total field is maximal and positive at the equator $\vartheta = \pi/2$, where the condensate is maximal, and it is minimal and negative at the poles $\vartheta = 0, \pi$, where the condensate vanishes (see Fig.25). Therefore, the background magnetic field is again enhanced by the condensate and the total field is maximal where the condensate is maximal.

<figure><img src="image_26.png" /><figcaption>Figure 25. The magnetic field $\delta B^r$ in (C.14), created by the condensate (C.11) with $j = 4$, the corresponding azimuthal current $4\pi J_\varphi$ in (C.12), and the condensate density ($R^2/2$) w are shown against $\vartheta$. Also shown the background magnetic $P/R^2 = m_w^2/e$, which corresponds to the threshold of the electroweak condensation.</figcaption></figure>

The current flow lines are parallel to $\vartheta = \text{const.}$ on the sphere, encircling two vortices centered at the poles $\vartheta = 0, \pi$. The condensate vanishes at the vortex centres, where the
total magnetic field $P/R^2 + \delta B^r$ is minimal, while the condensate is maximal at the equator, where the total magnetic field is maximal. Therefore, these vortices resemble “hollow tubes”. This agrees with the analysis in Sec. VI.

The solution (C.11)–(C.13) is still valid if we set $R = Q = \sqrt{\kappa/2} P$ and replace the geometry in (C.10) with the $AdS_2 \times S^2$ geometry (3.18) of Bertotti-Robinson. The latter describes the near-horizon limit of the extremal RN geometry. Therefore, Eqs. (C.11)–(C.13) describe the condensate on the extremal RN black hole, hence the bifurcation of the hairy black holes with the extremal RN black holes. This bifurcation corresponds to the rightmost point on the existence diagram in Fig. 22 (right panel), which relates to the maximal value of the charge $Q_{\text{max}} = 1.29 Q_{\text{m}} = 2g'\sqrt{\beta} Q_{\text{m}} = \sqrt{8/\kappa} g'/g$. The same value is obtained by setting $R = Q$ in (C.11), indeed,

$$Q^2 = \frac{\kappa}{2} P^2 = \frac{eP}{m_w^2} \Rightarrow P = \frac{2e}{\kappa m_w^2} \Rightarrow Q = \sqrt{\frac{\kappa}{2}} P = \sqrt{\frac{8}{\kappa}} \frac{g'}{g}. \quad (C.16)$$

## Appendix D. ESTIMATES OF THE MAGNETIC FIELD STRENGTH

It is instructive to evaluate the dimensionful magnetic field expressed in Teslas (T),

$$\left[ \frac{\mathcal{B}}{c} \right] = \frac{1 \text{ kg}}{1 \text{ Coulomb} \times 1 \text{ sec}} = 1 \text{ T} = 10^4 \text{ Gauss}; \quad (\text{D.1})$$

($[e\mathcal{B}]$ has the dimension of force). Consider the magnetic field corresponding to the onset of the W-condensation, $\mathcal{B} = m_w^2/e$. Its dimensionful version is

$$\frac{\mathcal{B}}{c} = \frac{1}{c} \frac{(m_w c^2)^2}{\hbar c e} = 1.097 \times 10^{20} \frac{1 \text{ kg}}{1 \text{ Coulomb} \times 1 \text{ sec}} \approx 1.1 \times 10^{20} \text{ T}, \quad (\text{D.2})$$

where we used the numerical values $e = 1.6 \times 10^{-19}$ Coulombs, $m_w c^2 = 80.4$ GeV, $1 \text{ GeV} = 1.78 \times 10^{-27} \text{ kg} \times c^2$, $c = 3 \times 10^8 \text{ m/sec}$ and $\hbar = 1.05 \times 10^{-34} \text{ kg m}^2/\text{sec}$.

Consider now the magnetic field of the Dirac monopole, $\mathcal{B} = n/(2er^2)$. Its dimensionful version is $\mathcal{B} = (\Phi_0/l_0)\mathcal{B}$ whereas $r = r/l_0$ while $\Phi_0 l_0 = 1/g_0 = e\mathbf{e}/(4\pi\alpha) = e\hbar\mathbf{c}/e$ (see Eqs. (2.6), (2.10)), hence

$$\frac{\mathcal{B}}{c} = n \times \frac{\hbar}{2e} \times \frac{1}{r^2} = \frac{n}{r^2} \times 3.28 \times 10^{-16} \text{ T} \times \text{m}^2. \quad (\text{D.3})$$

Setting here $r = 1.13 \times \sqrt{|n|} l_0$ according to (4.8) with $l_0 = 1.53 \times 10^{-18} \text{ m}$ yields again the same value as in (D.2). Setting $r$ to be the Sun's radius, $r = R_{\odot} \approx 6.96 \times 10^8 \text{ m}$, and choosing $n = 1.6 \times 10^{32}$, which corresponds to the charge of the maximally hairy black hole, yields

$$\frac{\mathcal{B}}{c} = n \times \frac{\hbar}{2e} \times \frac{1}{R_{\odot}^2} = n \times 6.77 \times 10^{-34} \text{ T} \approx 0.1 \text{ T} \approx 10^3 \text{ Gauss}. \quad (\text{D.4})$$

Therefore, a magnetic charge with $n = 1.6 \times 10^{32}$ placed at the center of the Sun would create at the surface a magnetic field of $\approx 10^3$ Gauss. For comparison, the average magnetic field at the Sun's surface is $\sim 1$ Gauss, but that in the sunspots can reach $2-4 \times 10^3$ Gauss.
Appendix E. HORIZON LIMIT OF EXTREMAL BLACK HOLES

The extremal hairy solutions in phase I reduce, in the vicinity of the horizon, to the extremal RNdS solution (3.1),(3.2),(3.19),(3.23)–(3.25), up to small corrections that vanish at the horizon. To determine these corrections, it is convenient to transform the solution to coordinates similar to those used in (3.17) for the extremal RN case.

The geometry of the extremal hairy solutions is described by the line element (7.2), with N(r) given by (8.2). Let us perform the coordinate transformation (7.3), assuming that in the new gauge, we have $\tilde{N} = 1/\tilde{r}^2$, and denote the new radial coordinate as $\tilde{r} = \eta$. This amounts to setting

$$d\eta = \frac{dr}{r\sqrt{N}} = \frac{dr}{k(r)(r-r_H)} \Rightarrow \eta - \eta_0 = \frac{1}{k(r)} \ln \frac{r-r_H}{r_H} + \int_{r_H}^{r} \frac{k'}{k^2(r)} \ln \frac{r-r_H}{r_H} dr. \quad (E.1)$$

As $r \to r_H$, we have $\eta \to -\infty$.

Re-expressing the metric amplitudes U, K, S in (7.2) in terms of three other functions, $\delta U, \delta K, \delta S$, as

$$e^{2U} = r_H^2 \eta^2 e^{2\lambda\eta+2\delta U}, \quad \eta^2 e^{2K-2U} = r_H^2 e^{2\delta K}, \quad \eta^2 e^{2S-2U} = r_H^2 e^{2\delta S}, \quad (E.2)$$

with a constant $\lambda$, the line element takes the form

$$ds^2 = r_H^2 \{-e^{2\lambda\eta}e^{2\delta U}dt^2 + e^{2\delta K}(d\eta^2 + d\vartheta^2) + e^{2\delta S}\sin^2\vartheta d\varphi^2\}. \quad (E.3)$$

We also re-express the electroweak amplitudes in (7.29),(7.31) as

$$\begin{aligned} H_1 &= \delta H_1, \quad H_2 = 1 + \delta H_2, \quad H_3 = \delta H_3, \quad H_4 = 1 + \delta H_4, \\ y &= \delta y, \quad \phi_1 = \delta \phi_1, \quad \phi_1 = \delta \phi_2. \end{aligned} \quad (E.4)$$

In these formulas, the deviations $\delta U, \dots, \delta \phi_2$ are not necessarily small, but they approach zero at the horizon, where $\eta \to -\infty$. Setting them to zero, $\delta U = \delta S = \dots = \delta \phi_2 = 0$, all field equations are satisfied provided that the constants $\lambda$ and $r_H$ obey the relations

$$r_H = r_{\text{ex}}, \quad \lambda^2 = 1 - \frac{\kappa\beta}{4} r_{\text{ex}}^2 = \sqrt{1 - \frac{\beta\kappa^2 n^2}{16g'^2}}. \quad (E.5)$$

These are precisely the same relations as for the extremal RNdS solution in (3.24). The argument of the square root is positive for $|Q| < Q_m$, which applies to all extremal hairy solutions in phase I, where $|Q| \le 0.6 Q_m$. This yields

$$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 (E.6)$$

which is an exact solution for any $\eta \in (-\infty, \infty)$. This is the same as the horizon limit of the extremal RNdS solution in (3.26).

In the vicinity of the horizon, the deviations $\delta U, \dots, \delta \phi_2$ are expected to be small, rather than zero. Therefore, one can expand the field equations around the limiting solution (E.6), assuming $\eta$ to be large and negative. This yields, in the first perturbative order, linear equations for the deviations. These equations split into three independent groups.
A. The first group of perturbations

This group contains the three metric amplitudes as well as $\delta y$, whose dependence on $\eta$ and $\vartheta$ can be separated as follows:

$$ \delta U = S_U(\vartheta)e^{\lambda\eta}, \quad \delta K = S_K(\vartheta)e^{\lambda\eta}, \quad \delta S = S_S(\vartheta)e^{\lambda\eta}, \quad \delta y = S_y(\vartheta)e^{\lambda\eta}. \qquad (E.7) $$

The functions $S_U(\vartheta)$, $S_K(\vartheta)$, $S_S(\vartheta)$, $S_y(\vartheta)$ satisfy equations, whose simplest solution is $\vartheta$-independent,

$$ S_K = S_S = \text{const.}, \quad S_U = -\frac{\lambda^2 + 2}{3\lambda^2} S_S, \quad S_y = 0. \qquad (E.8) $$

The condition $S_K = S_S$ implies that the geometry remains spherically symmetric in the first perturbative order. One can set $S_S = 1$ by adjusting the integration constant $\eta_0$ in (E.1) and then absorbing the constant $e^{\lambda\eta_0}$ by redefining the time coordinate. As a result, the Schwarzschild radius of the 2-sphere, $r = \sqrt{g_{\vartheta\vartheta}}$, is no longer constant:

$$ r = r_{\text{ex}} e^{\delta S}, \quad \text{where} \quad \delta S = e^{\lambda\eta} = \ln \frac{r}{r_{\text{ex}}} \approx \frac{r}{r_{\text{ex}}} - 1. \qquad (E.9) $$

The same results can be obtained by directly applying (E.1) to the line element of the extremal RNdS solution (3.25).

B. The second group of perturbations

This group includes the two Higgs amplitudes, which can be expressed as:

$$ \delta\phi_1 = S_1(\vartheta)e^{\sigma\eta}, \quad \delta\phi_2 = S_2(\vartheta)e^{\sigma\eta}, \qquad (E.10) $$

where $\sigma$ is a constant. The equations governing $S_1$ and $S_2$ are given by:

$$ \left( l(l+1) + \frac{d^2}{d\vartheta^2} + \cot\vartheta \frac{d}{d\vartheta} - \frac{\nu^2}{\sin^2\vartheta} + \frac{3\nu^2-1}{4} \right) S_1 - \left( \frac{d}{d\vartheta} + \frac{1-\nu^2}{2} \cot\vartheta \right) S_2 = 0, \\ \left( l(l+1) + \frac{d^2}{d\vartheta^2} + \cot\vartheta \frac{d}{d\vartheta} - \frac{\nu^2+1}{4} \right) S_2 + \left( \frac{d}{d\vartheta} + \frac{1+\nu^2}{2} \cot\vartheta \right) S_1 = 0, \qquad (E.11) $$

where

$$ \nu = -\frac{n}{2}, \quad l(l+1) = \sigma^2 + \lambda\sigma + \frac{\beta r_{\text{ex}}^2}{4}. \qquad (E.12) $$

This defines the eigenvalue problem to determine $l$, whose solution is

$$ l = \frac{\sqrt{1 + 2|\nu| - 1}}{2}, \quad S_1(\vartheta) = -\frac{2}{|\nu| + 1} \frac{d}{d\vartheta} S_2(\vartheta), \qquad (E.13) $$

where

$$ \left( \frac{d^2}{d\vartheta^2} - |\nu| \cot \vartheta \frac{d}{d\vartheta} + \frac{1-\nu^2}{4} \right) S_2 = 0 \quad \Rightarrow \quad S_2 = \left( \sin \frac{\vartheta}{2} \right)^{|\nu|+1} + \left( \cos \frac{\vartheta}{2} \right)^{|\nu|+1}. \qquad (E.14) $$
As a result, the exponent $\sigma$ in (E.10) is given by:

$$ \sigma = \frac{-\lambda + \sqrt{\lambda^2 + |n| - \beta r_{\text{ex}}^2}}{2}, \qquad (\text{E.15}) $$

and

$$ \phi_1 \propto \phi_2 \propto e^{\sigma\eta} = \left(\ln \frac{r}{r_{\text{ex}}}\right)^{\sigma/\lambda} \approx \left(\frac{r - r_{\text{ex}}}{r_{\text{ex}}}\right)^{\sigma/\lambda}. \qquad (\text{E.16}) $$

When the magnetic charge $|n|$ increases, $\sigma$ grows as $\sqrt{|n|}$, while $\lambda$ decreases. Therefore, the Higgs field approaches zero faster and faster in the vicinity of the horizon, and the size of the false vacuum region around the horizon increases.

If $|n| \ll 1/\sqrt{\kappa}$, then the horizon size $r_{\text{ex}} \ll 1$, and the spacetime geometry at the electroweak scale, where $r \sim 1$, is essentially flat. In that case, the size of the vacuum region can be estimated to be $\sim \sqrt{|n|}$, as expressed by Eq. (5.24) in [59].

The relation (E.16) implies that the components of the unit vector $n^a = (\Phi^\dagger\tau^a\Phi)/(\Phi^\dagger\Phi)$ used in (2.26) to define the electromagnetic field do not depend on the radial coordinate:

$$ n^1 = \frac{2S_1S_2}{(S_1)^2 + (S_2)^2}, \qquad n^2 = 0, \qquad n^3 = \frac{(S_1)^2 - (S_2)^2}{(S_1)^2 + (S_2)^2}, \qquad (\text{E.17}) $$

where we used $\Phi$ in the gauge (7.22). These relations are exact at the horizon, while in the vicinity of the horizon, they may receive corrections only in higher perturbative orders.

This allows us to determine the horizon value of the electromagnetic field. Since the perturbations $\delta H_1, \delta H_2, \delta H_3, \delta H_4$ approach zero at the horizon, it follows from (7.22),(7.27) that the SU(2) gauge field strength vanishes at the horizon, $W_{\mu\nu}^a = 0$. The Nambu field $\mathcal{F}_{\mu\nu}$, defined in (2.26), then reduces to the U(1) hypercharge field, up to a constant factor.

Therefore, the only non-zero components of the Nambu field at the horizon are $\mathcal{F}_{\vartheta\varphi} = -\mathcal{F}_{\varphi\vartheta}$, where

$$ \mathcal{F}_{\vartheta\varphi} = \frac{g}{g'} B_{\vartheta\varphi} = \frac{g}{g'} \frac{n}{2} \sin\vartheta = g^2 \frac{n}{2e} \sin\vartheta = g^2 P \sin\vartheta. \qquad (\text{E.18}) $$

This corresponds to the purely radial magnetic field of total flux $4\pi g^2 P$, created by the magnetic charge $g^2 P$ inside the black hole. The remaining part of the charge $g'^2 P$ is contained in the black hole hair.

The 't Hooft field, $F_{\mu\nu} = \mathcal{F}_{\mu\nu} - e\psi_{\mu\nu}$, defined in (2.36), includes the condensate tensor $\psi_{\mu\nu}$, which is constructed from the unit vector $n^a$ according to (2.35). Using $n^a$ in (E.17) and the SU(2) field $W_\mu^a$ in the gauge (7.22), we find that the only components of the condensate which do not vanish at the horizon are $\psi_{\vartheta\varphi} = -\psi_{\varphi\vartheta}$, where

$$ -e\psi_{\vartheta\varphi} = g'^2 \frac{n}{2e} \times \frac{|n| (\sin \frac{\vartheta}{2} \cos \frac{\vartheta}{2})^{|n|-1}}{[(\cos \frac{\vartheta}{2})^{|n|} + (\sin \frac{\vartheta}{2})^{|n|}]^2}. \qquad (\text{E.19}) $$

It is instructive to compare this with the formula $\psi_{\vartheta\varphi} \sim (\sin\vartheta)^{|n|-1}$, which follows from (6.14),(6.15) in the axially symmetric case (when $c_k = \delta_{0k}$ in (6.14)). This formula was obtained using perturbation theory at the bifurcation with the RN solution. In contrast,
the relation in (E.19) is non-perturbative and corresponds to the fully non-linear solution restricted to the event horizon. According to (6.14),(6.15), the value of $|\psi_{\vartheta\varphi}|$ (or more precisely, $|\psi_{\vartheta\varphi}|/\sin\vartheta$) can be interpreted as a measure of the condensate strength. Remarkably, both the perturbative and non-perturbative expressions for the condensate exhibit the same qualitative behavior: they reach their maximum at the equator and vanish at the poles.

The flux of $\psi_{\mu\nu}$ is given by

$$ -e \oint \psi_{\vartheta\varphi} d\vartheta d\varphi = 4\pi g'^2 P. \qquad (E.20) $$

As a result, the horizon value of the 't Hooft field is the sum of (E.18) and (E.19),

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

whose flux is $4\pi P$. Since the 2-form $F_{\mu\nu}$ is closed, its flux is independent of the integration surface. As a result, the flux though a 2-sphere at spatial infinity is the same as that though the horizon, implying that the entire magnetic charge is contained within the black hole. In contrast, for the Nambu field (E.18), the flux through the horizon is $4\pi g^2 P$, while the flux at infinity remains $4\pi P$. The difference is $4\pi g'^2 P$, corresponds to the magnetic charge $g'^2 P$ distributed outside the black hole.

The Z field, defined by (2.36), has non-zero components

$$ Z_{\vartheta\varphi} = -Z_{\varphi\vartheta} = B_{\vartheta\varphi} + g^2 \psi_{\vartheta\varphi}. \qquad (E.22) $$

Its horizon flux vanishes, as expected, since the 2-form $Z_{\mu\nu}$ is also closed, and its flux at infinity is zero because this field is massive.

Let us now consider the current densities. Taking the divergence of the Nambu field (E.18) gives zero, meaning that the Nambu current vanishes at the horizon. The magnetic charge density in (2.49) is proportional to the radial derivative of $\psi_{\vartheta\varphi}$, which is also zero at the horizon. Therefore, both the Nambu current density $J^\mu$ and the magnetic charge density $\tilde{J}^0$ vanish at the horizon, although they deviate from zero at some distance away, as indicated by our numerics.

On the other hand, taking the divergence of the 't Hooft field (E.21) gives a non-zero value of the azimuthal current density at the horizon:

$$ 4\pi r_{\text{ex}}^2 J_{\varphi} = - \sin\vartheta \left( \frac{F_{\vartheta\varphi}}{\sin\vartheta} \right)'_{\vartheta} = - \frac{g' n |n|}{g 2^{|n|}} \sin\vartheta \left( \frac{(\sin\vartheta)^{|n|-2}}{\left[ (\cos\frac{\vartheta}{2})^{|n|} + (\sin\frac{\vartheta}{2})^{|n|} \right]^2} \right)'_{\vartheta}. \qquad (E.23) $$

According to our numerics (see Fig.19, right panel) the current density is maximal at the horizon and starts to decrease as we move away from it. Notice that the formulas (E.18)–(E.23) do not contain $\kappa$.

### C. The third group of perturbations

This group contains four SU(2) amplitudes that we re-express as:

$$ \delta H_1 = e^{q\eta} h_1(\vartheta), \quad \delta H_2 = e^{q\eta} h_2(\vartheta), \quad \delta H_3 = e^{q\eta} h_3(\vartheta), \quad \delta H_4 = e^{q\eta} h_4(\vartheta), \qquad (E.24) $$
where $q$ is a constant. The corresponding equations for the perturbations are:

$$ 
\begin{aligned}
& \left( \frac{d^2}{d\vartheta^2} + \cot \vartheta \frac{d}{d\vartheta} - \frac{\nu^2}{\sin^2 \vartheta} + q^2 \right) h_1 + q \cot \vartheta h_2 + q\nu^2 h_3 - q\nu^2 \cot \vartheta h_4 = 0, \\
& \left( \frac{d^2}{d\vartheta^2} - \frac{\nu^2}{\sin^2 \vartheta} + q(q + \lambda) \right) h_2 + \lambda \frac{d}{d\vartheta} h_1 - \nu^2 \left( \frac{d}{d\vartheta} + 2 \cot \vartheta \right) h_3 \\
& \qquad + \nu^2 (\cot \vartheta \frac{d}{d\vartheta} + \cot^2 \vartheta - 1) h_4 = 0, \quad (\text{E.25}) \\
& \left( \frac{d^2}{d\vartheta^2} + \cot \vartheta \frac{d}{d\vartheta} - \frac{1}{\sin^2 \vartheta} + q(q + \lambda) - 1 \right) h_3 - \lambda h_1 - \cot \vartheta h_2 + (2 \frac{d}{d\vartheta} + \cot \vartheta) h_4 = 0, \\
& \left( \frac{d^2}{d\vartheta^2} + \cot \vartheta \frac{d}{d\vartheta} - \frac{1}{\sin^2 \vartheta} + q(q + \lambda) - 1 \right) h_4 + \lambda \cot \vartheta h_1 + \cot^2 \vartheta h_2 - (2 \frac{d}{d\vartheta} + \cot \vartheta) h_3 = 0.
\end{aligned}
$$

By using the transformations:

$$ h_4(\vartheta) = \cos \vartheta \tilde{h}_3(\vartheta) + \sin \vartheta \tilde{h}_4(\vartheta), \quad h_3(\vartheta) = -\sin \vartheta \tilde{h}_3(\vartheta) + \cos \vartheta \tilde{h}_4(\vartheta), \quad (\text{E.26}) $$

the system of equations (E.25) splits into three coupled equations involving $h_1, h_2, \tilde{h}_3$, plus one independent equation for $\tilde{h}_4$:

$$ \left( \frac{d^2}{d\vartheta^2} + \cot \vartheta \frac{d}{d\vartheta} - \frac{1}{\sin^2 \vartheta} + q(q + \lambda) \right) \tilde{h}_4 = 0. \quad (\text{E.27}) $$

The solution to this equation is given by:

$$ \tilde{h}_4 = C_4 P_l^1(\cos \vartheta), \quad q(q + \lambda) = l(l + 1), \quad l = 1, 2, \dots, \quad (\text{E.28}) $$

where $P_l^1(\cos \vartheta)$ are the Legendre polynomials and $C_4$ is an integration constant.

By setting $h_1 = h_2 = \tilde{h}_3 = 0$, and using the equation (E.26), we obtain

$$ h_4 = C_4 \sin \vartheta P_l^1(\cos \vartheta), \quad h_3 = C_4 \cos \vartheta P_l^1(\cos \vartheta), \quad (\text{E.29}) $$

and the value of $q$ is

$$ q = \frac{-\lambda + \sqrt{\lambda^2 + 4l(l+1)}}{2}, \quad (\text{E.30}) $$

where $C_4$ is an integration constant.

The analysis of the remaining equations for $h_1, h_2, \tilde{h}_3$ involves applying the gauge fixing condition (7.33), which now simplifies to:

$$ qh_1(\vartheta) = h'_2(\vartheta). \quad (\text{E.31}) $$

This allows us to eliminate $h_1$ since its equation becomes a differential consequence of the two other equations. The remaining equations reduce to

$$ 
\begin{aligned}
& \left( \frac{d^2}{d\vartheta^2} + \cot \vartheta \frac{d}{d\vartheta} - \frac{1}{\sin^2 \vartheta} + q(q + \lambda) \right) \tilde{h}_3 + \frac{1}{\sin \vartheta} \left( \frac{\lambda}{q} \frac{d}{d\vartheta} + \cot \vartheta \right) h_2 = 0, \\
& \left( (q + \lambda) \frac{d^2}{d\vartheta^2} + q^2(q + \lambda) - \frac{q \nu^2}{\sin^2 \vartheta} \right) h_2 + \frac{q \nu^2}{\sin \vartheta} \left( \frac{d}{d\vartheta} + \cot \vartheta \right) \tilde{h}_3 = 0.
\end{aligned}
\quad (\text{E.32}) $$

Unfortunately, we were unable to find an exact solution of these equations, although it is clear that the solution should be expressible in terms of Legendre polynomials. We leave this problem for future work.
Appendix F. ELECTROWEAK VORTICES

The vortices were originally discovered [89, 90] within the Abelian Higgs model, which includes an Abelian vector $C_μ$ and a complex scalar $\phi$ with the Lagrangian

$$ \mathcal{L} = -\frac{1}{4} C_{\mu\nu} C^{\mu\nu} - \overline{D_{\mu} \phi} D^{\mu} \phi - \frac{\beta}{8} (|\phi|^2 - 1)^2. \qquad (F.1) $$

Here $C_{\mu\nu} = \partial_{\mu}C_{\nu} - \partial_{\nu}C_{\mu}$ and $D_{\mu}\phi = (\nabla_{\mu} - iqC_{\mu})\phi$ with $q$ being the electric charge. The equations of this theory,

$$ \nabla^{\mu} C_{\mu\nu} = iq(\overline{\phi}D_{\nu}\phi - \phi\overline{D_{\nu}\phi}), \qquad D_{\mu}D^{\mu}\phi = \frac{\beta}{4}(|\phi|^2 - 1)\phi, \qquad (F.2) $$

admit solutions describing a localized magnetic flux along the $z$ direction. In the orthogonal $\mathbb{R}^2$ plane, the fields approach the Higgs vacuum at infinity, $\phi \to e^{i\alpha}$, $C_\mu \to \partial_\mu\alpha/q$. The phase $\alpha$ should wind an integer number of times for $\phi$ to be single-valued, hence $\alpha = k\varphi$, where $\varphi$ is the azimuthal angle and $k \in \mathbb{Z}$. This determines the magnetic flux carried by the vortex,

$$ \Psi_C = \oint_{\partial\mathbb{R}^2} C_k dx^k = \frac{2\pi k}{q}. \qquad (F.3) $$

These vortex solutions can be embedded into the electroweak theory. There are several possibilities for this, always assuming that $W_\mu^1 = W_\mu^2 = 0$.

1. Z-strings

If $C_\mu$, $\phi$ fulfill (F.2), then the electroweak equations (2.17) will be fulfilled by setting [84]

$$ B_{\mu} = g'^{2}C_{\mu}, \quad W_{\mu}^{3} = -g^{2}C_{\mu}, \quad \phi_{1} = 0, \quad \phi_{2} = \phi, \quad q = \frac{1}{2}. \qquad (F.4) $$

Injecting into (2.38) yields $A_\mu = 0$, $Z_\mu = C_\mu$, hence the vortex contains in this case only the massive Z-field, while the massless F-field vanishes. Using (F.3) determines the fluxes:

$$ \Psi_F = 0, \qquad \Psi_Z = 4\pi k. \qquad (F.5) $$

2. W-strings

Another possibility to fulfill the electroweak equations (2.17) is to set [85]

$$ B_{\mu} = 0, \quad W_{\mu}^{3} = g C_{\mu}, \quad \phi_{1} = \overline{\phi_{2}} = \frac{\phi}{\sqrt{2}}, \quad q = \frac{g}{2}. \qquad (F.6) $$

This time, both the massless and massive fields are non-zero, with fluxes:

$$ \Psi_F = g'^2 \frac{4\pi k}{e}, \qquad \Psi_Z = -4\pi k. \qquad (F.7) $$
### 3. B-strings

Yet one more option is to choose

$$B_{\mu} = g' C_{\mu}, \quad W_{\mu}^{3} = 0, \quad \phi_{1} = \phi_{2} = \frac{\phi}{\sqrt{2}}, \quad q = \frac{g'}{2}, \qquad (F.8)$$

which yields the fluxes:

$$\Psi_F = g^2 \frac{4\pi k}{e}, \quad \Psi_Z = 4\pi k. \qquad (F.9)$$

These possibilities correspond to three physically different types of vortices. They are infinitely long, but within the electroweak theory, there are more general solutions describing semi-infinite or finite vortex pieces. In other words, electroweak vortices can terminate. If a vortex terminates, the magnetic flux $\Psi_F$ trapped inside spreads out from the termination point, producing a radial magnetic field that imitates a magnetic monopole. For example, the extremity of a W-string corresponds to the Nambu monopole, whose charge is proportional to $g'^2/e$ [61].

Comparing the above formulas with those for the horizon vortices, we see that the fluxes in (10.17) reduce to those for W-strings in (F.7) upon setting $n = -4k$. Therefore, the horizon vortices ressemble W-strings; in particular, they have the same value for the ratio $\Psi_F/\Psi_Z$. At the same time, the analogy is not complete since $n/4$ is not always integer for the horizon vortices. Additionally, as explained in Section XA 3, the horizon vortices are made of massless fields, whereas in W-strings, the fields become massive far away from the center.

[1] R. Ruffini and J. Wheeler, *Introducing the black hole*, Phys. Today **24** (1971) 30–41, [1206.3482], [doi:10.1088/1475-7516/2013/01/035].

[2] W. Israel, *Event horizons in static vacuum space-times*, Phys. Rev. **164** (1967) 1776–1779, [doi:10.1103/PhysRev.164.1776].

[3] D. C. Robinson, *Uniqueness of the Kerr black hole*, Phys. Rev. Lett. **34** (1975) 905–906, [doi:10.1103/PhysRevLett.34.905].

[4] P. O. Mazur, *Proof of uniqueness of the Kerr-Newman black hole solution*, J. Phys. A **15** (1982) 3173–3180, [doi:10.1088/0305-4470/15/10/021].

[5] M. Heusler, *Black Hole Uniqueness Theorems*. Cambridge University Press, 1996.

[6] J. E. Chase, *Event horizons in static scalar-vacuum space-times*, Commun. Math. Phys. **19** (1970), no. 4 276–288, [doi:10.1007/BF01646635].

[7] J. D. Bekenstein, *Nonexistence of baryon number for static black holes*, Phys. Rev. D **5** (1972) 1239–1246, [doi:10.1103/PhysRevD.5.1239].

[8] J. D. Bekenstein, *Transcendence of the law of baryon-number conservation in black hole physics*, Phys. Rev. Lett. **28** (1972) 452–455, [doi:10.1103/PhysRevLett.28.452].

[9] J. D. Bekenstein, *Novel “no-scalar-hair” theorem for black holes*, Phys. Rev. D **51** (1995), no. 12 R6608, [doi:10.1103/PhysRevD.51.R6608].

[10] A. E. Mayo and J. D. Bekenstein, *No hair for spherical black holes: Charged and nonminimally coupled scalar field with selfinteraction*, Phys. Rev. D **54** (1996) 5059–5069, [arXiv:gr-qc/9602057], [doi:10.1103/PhysRevD.54.5059].
[11] J. D. Bekenstein, *Black hole hair: 25 - years after*, in *2nd International Sakharov Conference on Physics*, pp. 216-219, 5, 1996. arXiv:gr-qc/9605059.

[12] A. Saa, *New no scalar hair theorem for black holes*, J. Math. Phys. **37** (1996) 2346-2351,
[arXiv:gr-qc/9601021], [doi:10.1063/1.531513].

[13] C. A. R. Herdeiro and E. Radu, *Asymptotically flat black holes with scalar hair: a review*, Int. J. Mod. Phys. D **24** (2015), no. 09 1542014, [arXiv:1504.08209],
[doi:10.1142/S0218271815420146].

[14] M. S. Volkov and D. V. Galtsov, *NonAbelian Einstein Yang-Mills black holes*, [Pis'ma Zh. Eksp. Teor. Fiz. **50** (1989) 312.], [JETP Lett. 50, 346 (1989)].

[15] M. S. Volkov, *Hairy black holes in the XX-th and XXI-st centuries*, arXiv:1601.08230.

[16] H. P. Kuenzle and A. K. M. Masood-ul Alam, *Spherically symmetric static SU(2) Einstein Yang-Mills fields*, J. Math. Phys. **31** (1990) 928-935, [doi:10.1063/1.528773].

[17] P. Bizon, *Colored black holes*, Phys. Rev. Lett. **64** (1990) 2844-2847,
[doi:10.1103/PhysRevLett.64.2844].

[18] S. Droz, M. Heusler, and N. Straumann, *New black hole solutions with hair*, Phys. Lett. B **268** (1991) 371-376, [doi:https://doi.org/10.1016/0370-2693(91)91592-J].

[19] P. Bizon and T. Chmaj, *Gravitating skyrmions*, Phys. Lett. B **297** (1992) 55-62,
[doi:10.1016/0370-2693(92)91069-L].

[20] P. Breitenlohner, P. Forgács, and D. Maison, *Gravitating monopole solutions*, Nucl. Phys. B **383** (1992) 357-376, [doi:https://doi.org/10.1016/0550-3213(92)90682-2].

[21] B. R. Greene, S. D. Mathur, and C. M. O'Neill, *Eluding the no hair conjecture: Black holes in spontaneously broken gauge theories*, Phys. Rev. D **47** (1993) 2242-2259,
[arXiv:hep-th/9211007], [doi:10.1103/PhysRevD.47.2242].

[22] G. Lavrelashvili and D. Maison, *Regular and black hole solutions of Einstein-Yang-Mills dilaton theory*, Nucl. Phys. B **410** (1993) 407-422,
[doi:https://doi.org/10.1016/0550-3213(93)90441-Q].

[23] S. Mignemi and N. R. Stewart, *Charged black holes in effective string theory*, Phys. Rev. D **47** (1993) 5259-5269, [arXiv:hep-th/9212146], [doi:10.1103/PhysRevD.47.5259].

[24] P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, *Dilatonic black holes in higher curvature string gravity*, Phys. Rev. D **54** (1996) 5049-5058, [arXiv:hep-th/9511071],
[doi:10.1103/PhysRevD.54.5049].

[25] T. Torii, H. Yajima, and K.-i. Maeda, *Dilatonic black holes with Gauss-Bonnet term*, Phys. Rev. D **55** (1997) 739-753, [arXiv:gr-qc/9606034], [doi:10.1103/PhysRevD.55.739].

[26] B. Kleihaus and J. Kunz, *Static black-hole solutions with axial symmetry*, Phys. Rev. Lett. **79** (Sep, 1997) 1595-1598, [doi:10.1103/PhysRevLett.79.1595].

[27] B. Kleihaus and J. Kunz, *Rotating hairy black holes*, Phys. Rev. Lett. **86** (2001) 3704-3707,
[arXiv:gr-qc/0012081], [doi:10.1103/PhysRevLett.86.3704].

[28] M. S. Volkov and D. V. Gal'tsov, *Gravitating nonAbelian solitons and black holes with Yang-Mills fields*, Phys. Rept. **319** (1999) 1-83, [arXiv:hep-th/9810070],
[doi:10.1016/S0370-1573(99)00010-1].

[29] U. Nucamendi and M. Salgado, *Scalar hairy black holes and solitons in asymptotically flat space-times*, Phys. Rev. D **68** (2003) 044026, [arXiv:gr-qc/0301062],
[doi:10.1103/PhysRevD.68.044026].

[30] K. G. Zloshchastiev, *On co-existence of black holes and scalar field*, Phys. Rev. Lett. **94** (2005) 121101, [arXiv:hep-th/0408163], [doi:10.1103/PhysRevLett.94.121101].

[31] A. Anabalon, D. Astefanesei, and R. Mann, *Exact asymptotically flat charged hairy black
holes with a dilaton potential, JHEP 10 (2013) 184, [arXiv:1308.1693],
[doi:10.1007/JHEP10(2013)184].

[32] M. Rinaldi, *Black holes with non-minimal derivative coupling*, Phys. Rev. D **86** (2012) 084048, [arXiv:1208.0103], [doi:10.1103/PhysRevD.86.084048].

[33] T. P. Sotiriou and S.-Y. Zhou, *Black hole hair in generalized scalar-tensor gravity: An explicit example*, Phys. Rev. D **90** (Dec, 2014) 124063, [doi:10.1103/PhysRevD.90.124063].

[34] E. Babichev and C. Charmousis, *Dressing a black hole with a time-dependent Galileon*, JHEP **08** (2014) 106, [arXiv:1312.3204], [doi:10.1007/JHEP08(2014)106].

[35] T. Kobayashi and N. Tanahashi, *Exact black hole solutions in shift symmetric scalar-tensor theories*, PTEP **2014** (2014) 073E02, [arXiv:1403.4364], [doi:10.1093/ptep/ptu096].

[36] C. A. R. Herdeiro and E. Radu, *Kerr black holes with scalar hair*, Phys. Rev. Lett. **112** (2014) 221101, [arXiv:1403.2757], [doi:10.1103/PhysRevLett.112.221101].

[37] C. Herdeiro, E. Radu, and H. Rúnarsson, *Kerr black holes with Proca hair*, Class. Quant. Grav. **33** (2016), no. 15 154001, [arXiv:1603.02687], [doi:10.1088/0264-9381/33/15/154001].

[38] W. E. East and F. Pretorius, *Superradiant instability and backreaction of massive vector fields around Kerr black holes*, Phys. Rev. Lett. **119** (2017), no. 4 041101, [arXiv:1704.04791], [doi:10.1103/PhysRevLett.119.041101].

[39] D. D. Doneva, F. M. Ramazanoğlu, H. O. Silva, T. P. Sotiriou, and S. S. Yazadjiev, *Spontaneous scalarization*, Rev. Mod. Phys. **96** (2024), no. 1 015004, [arXiv:2211.01766], [doi:10.1103/RevModPhys.96.015004].

[40] R. Brito, V. Cardoso, and P. Pani, *Black holes with massive graviton hair*, Phys. Rev. D **88** (2013) 064006, [arXiv:1309.0818], [doi:10.1103/PhysRevD.88.064006].

[41] R. Gervalle and M. S. Volkov, *Asymptotically flat hairy black holes in massive bigravity*, Phys. Rev. D **102** (2020), no. 12 124040, [arXiv:2008.13573], [doi:10.1103/PhysRevD.102.124040].

[42] E. Alonso-Monsalve and D. I. Kaiser, *Primordial black holes with QCD color charge*, Phys. Rev. Lett. **132** (2024) 231402, [arXiv:2310.16877], [doi:10.1103/PhysRevLett.132.231402].

[43] Y. Bai and M. Korwar, *Hairy magnetic and dyonic black holes in the Standard Model*, JHEP **04** (2021) 119, [arXiv:2012.15430], [doi:10.1007/JHEP04(2021)119].

[44] K.-M. Lee and E. J. Weinberg, *Nontopological magnetic monopoles and new magnetically charged black holes*, Phys. Rev. Lett. **73** (1994) 1203-1206, [arXiv:hep-th/9406021], [doi:10.1103/PhysRevLett.73.1203].

[45] S. A. Ridgway and E. J. Weinberg, *Static black hole solutions without rotational symmetry*, Phys. Rev. D **52** (1995) 3440-3456, [arXiv:gr-qc/9503035], [doi:10.1103/PhysRevD.52.3440].

[46] J. Maldacena, *Comments on magnetic black holes*, JHEP **04** (2021) 079, [arXiv:2004.06084], [doi:10.1007/JHEP04(2021)079].

[47] J. Ambjorn and P. Olesen, *On electroweak magnetism*, Nucl. Phys. B **315** (1989) 606-614, [doi:10.1016/0550-3213(89)90004-7].

[48] J. Ambjorn and P. Olesen, *Electroweak magnetism: theory and application*, Int. J. Mod. Phys. A **5** (1990) 4525-4558, [doi:10.1142/S0217751X90001914].

[49] M. N. Chernodub, J. Van Doorsselaere, and H. Verschelde, *Magnetic-field-induced superconductivity and superfluidity of W and Z bosons: in tandem transport and kaleidoscopic vortex states*, Phys. Rev. D **88** (2013) 065006, [arXiv:1203.5963], [doi:10.1103/PhysRevD.88.065006].

[50] M. N. Chernodub, V. A. Goy, and A. V. Molochkov, *Phase structure of electroweak vacuum in a strong magnetic field: the lattice results*, Phys. Rev. Lett. **130** (2023), no. 11 111802, [arXiv:2206.14008], [doi:10.1103/PhysRevLett.130.111802].
[51] Y. Bai, J. Berger, M. Korwar, and N. Orlofsky, *Phenomenology of magnetic black holes with electroweak-symmetric coronas*, JHEP **10** (2020) 210, [arXiv:2007.03703], [doi:10.1007/JHEP10(2020)210].

[52] Y. Bai, S. Lu, and N. Orlofsky, *Searching for magnetic monopoles with the Earth's magnetic field*, Phys. Rev. Lett. **127** (2021), no. 10 101801, [arXiv:2103.06286], [doi:10.1103/PhysRevLett.127.101801].

[53] D. Ghosh, A. Thalapillil, and F. Ullah, *Astrophysical hints for magnetic black holes*, Phys. Rev. D **103** (2021), no. 2 023006, [arXiv:2009.03363], [doi:10.1103/PhysRevD.103.023006].

[54] J. Estes, M. Kavic, S. L. Liebling, M. Lippert, and J. H. Simonetti, *Stability and observability of magnetic primordial black hole-neutron star collisions*, JCAP **06** (2023) 017, [arXiv:2209.06060], [doi:10.1088/1475-7516/2023/06/017].

[55] M. D. Diamond and D. E. Kaplan, *Constraints on relic magnetic black holes*, JHEP **03** (2022) 157, [arXiv:2103.01850], [doi:10.1007/JHEP03(2022)157].

[56] Y. M. Cho, S.-W. Kim, and S. H. Oh, *Electroweak primordial magnetic blackhole: cosmic production and physical implication*, arXiv:2408.05531.

[57] R. Gervalle and M. S. Volkov, *Black holes with electroweak hair*, Phys. Rev. Lett. **133** (2024), no. 17 171402, [arXiv:2406.14357], [doi:10.1103/PhysRevLett.133.171402].

[58] R. Gervalle and M. S. Volkov, *Electroweak monopoles and their stability*, Nucl. Phys. B **984** (2022) 115937, [arXiv:2203.16590], [doi:10.1016/j.nuclphysb.2022.115937].

[59] R. Gervalle and M. S. Volkov, *Electroweak multi-monopoles*, Nucl. Phys. B **987** (2023) 116112, [arXiv:2211.04875], [doi:10.1016/j.nuclphysb.2023.116112].

[60] E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, *CODATA recommended values of the fundamental physical constants: 2018**, Rev. Mod. Phys. **93** (2021), no. 2 025010, [doi:10.1103/RevModPhys.93.025010].

[61] Y. Nambu, *String-like configurations in the Weinberg-Salam theory*, Nucl. Phys. B **130** (1977) 505, [doi:10.1016/0550-3213(77)90252-8].

[62] G. 't Hooft, *Magnetic monopoles in unified gauge theories*, Nucl. Phys. B **79** (1974) 276-284, [doi:10.1016/0550-3213(74)90486-6].

[63] P. A. M. Dirac, *Quantised singularities in the electromagnetic field*, Proc. Roy. Soc. Lond. A **133** (1931), no. 821 60-72, [doi:10.1098/rspa.1931.0130].

[64] B. Bertotti, *Uniform electromagnetic field in the theory of general relativity*, Phys. Rev. **116** (1959) 1331, [doi:10.1103/PhysRev.116.1331].

[65] I. Robinson, *A Solution of the Maxwell-Einstein equations*, Bull. Acad. Pol. Sci. Ser. Sci. Math. Astron. Phys. **7** (1959) 351-352.

[66] S. Hod, *Electroweak clouds supported by magnetically charged black holes: Analytic treatment along the existence line*, Phys. Rev. D **110** (2024), no. 10 104045, [arXiv:2501.18702], [doi:10.1103/PhysRevD.110.104045].

[67] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, *Numerical Recipes 3rd Edition: The Art of Scientific Computing*. Cambridge University Press, New York, NY, USA, 3rd ed., 2007.

[68] J. Ambjorn and P. Olesen, *A Condensate solution of the electroweak theory which interpolates between the broken and the symmetric phase*, Nucl. Phys. B **330** (1990) 193-204, [doi:10.1016/0550-3213(90)90307-Y].

[69] N. Arkani-Hamed, L. Motl, A. Nicolis, and C. Vafa, *The String landscape, black holes and gravity as the weakest force*, JHEP **06** (2007) 060, [arXiv:hep-th/0601001], [doi:10.1088/1126-6708/2007/06/060].

[70] Y. M. Cho and D. Maison, *Monopoles in Weinberg-Salam model*, Phys. Lett. B **391** (1997)
360–365, [arXiv:hep-th/9601028], [doi:10.1016/S0370-2693(96)01492-X].

[71] J. L. Synge, ed., *Relativity: The General theory*. North-Holland Publishing Company, Amsterdam, 1960.

[72] A. Komar, *Covariant conservation laws in general relativity*, Phys. Rev. **113** (1959) 934–936, [doi:10.1103/PhysRev.113.934].

[73] C. Rebbi and P. Rossi, *Multi-monopole solutions in the Prasad-Sommerfield limit*, Phys. Rev. D **22** (1980) 2010, [doi:10.1103/PhysRevD.22.2010].

[74] F. Hecht, *New development in freefem++*, J. Numer. Math. **20** (2012), no. 3-4 251–265.

[75] S. Deser, *Absence of static solutions in source-free Yang-Mills theory*, Phys. Lett. B **64** (1976) 463–464, [doi:10.1016/0370-2693(76)90122-2].

[76] G. H. Derrick, *Comments on nonlinear wave equations as models for elementary particles*, J. Math. Phys. **5** (1964) 1252–1254, [doi:10.1063/1.1704233].

[77] C. A. R. Herdeiro, J. a. M. S. Oliveira, A. M. Pombo, and E. Radu, *Deconstructing scaling virial identities in general relativity: Spherical symmetry and beyond*, Phys. Rev. D **106** (2022), no. 2 024054, [arXiv:2206.02813], [doi:10.1103/PhysRevD.106.024054].

[78] R. P. Geroch, *Multipole moments. II. Curved space*, J. Math. Phys. **11** (1970) 2580–2588, [doi:10.1063/1.1665427].

[79] R. Hansen, *Multipole moments of stationary spacetimes*, J. Math. Phys. **15** (1974) 46–52, [doi:10.1063/1.1666501].

[80] G. Fodor, C. Hoenselaers, and Z. Perjes, *Multipole moments of axisymmetric systems in relativity*, J. Math. Phys. **30** (1989) 2252–2257, [doi:10.1063/1.528551].

[81] G. Fodor, E. d. S. C. Filho, and B. Hartmann, *Calculation of multipole moments of axistationary electrovacuum spacetimes*, Phys. Rev. D **104** (2021), no. 6 064012, [arXiv:2012.05548], [doi:10.1103/PhysRevD.104.064012].

[82] F. J. Ernst, *New formulation of the axially symmetric gravitational potential*, Phys. Rev. **167** (1968) 1175–1167, [doi:10.1103/PhysRev.167.1175].

[83] F. J. Ernst, *New Formulation of the axially symmetric gravitational field problem. II*, Phys. Rev. **168** (1968) 1415–1417, [doi:10.1103/PhysRev.168.1415].

[84] T. Vachaspati, *Vortex solutions in the Weinberg-Salam model*, Phys. Rev. Lett. **68** (1992) 1977–1980, [doi:10.1103/PhysRevLett.68.1977]. [Erratum: Phys. Rev. Lett. 69, 216 (1992)].

[85] T. Vachaspati and M. Barriola, *A New class of defects*, Phys. Rev. Lett. **69** (1992) 1867–1870, [doi:10.1103/PhysRevLett.69.1867].

[86] P. Olesen, *Vacuum structure of the electroweak theory in high magnetic fields*, Phys. Lett. B **268** (1991) 389–393, [doi:10.1016/0370-2693(91)91595-M].

[87] J. Van Doorsselaere, *Properties of symmetry restoration in the electroweak vacuum in very high magnetic fields*, Phys. Rev. D **88** (2013), no. 2 025013, [arXiv:1206.6205], [doi:10.1103/PhysRevD.88.025013].

[88] C. A. R. Herdeiro and E. Radu, *Two Schwarzschild-like black holes balanced by their scalar hair*, Phys. Rev. D **107** (2023), no. 6 064044, [arXiv:2302.00016], [doi:10.1103/PhysRevD.107.064044].

[89] A. A. Abrikosov, *On the magnetic properties of superconductors of the second group*, Sov. Phys. JETP **5** (1957) 1174–1182.

[90] H. B. Nielsen and P. Olesen, *Vortex line models for dual strings*, Nucl. Phys. B **61** (1973) 45–61, [doi:10.1016/0550-3213(73)90350-7].