# Non-Abelian Wormholes Threaded by a Yang-Mills-Higgs Field beyond the BPS
Limit

Xiao Yan Chew<EMAIL_ADDRESS>Department of Physics Education, Pusan
National University, Busan 46241, Republic of Korea  Research Center for
Dielectric and Advanced Matter Physics, Pusan National University, Busan
46241, Republic of Korea  Kok-Geng Lim<EMAIL_ADDRESS>University of
Southampton Malaysia, 79200 Iskandar Puteri, Johor, Malaysia

###### Abstract

We construct numerically the symmetric non-Abelian wormholes which are
supported by a phantom field in the Einstein-Yang-Mills-Higgs theory beyond
Bogomol’nyi-Prasad-Sommerfield (BPS) limit where the Higgs self-interaction
constant $\lambda$ is non-vanishing. Analogous to the BPS limit, the probe
limit is the Yang-Mills-Higgs field in the background of the Ellis wormhole
when the gravity is switched off. In the presence of gravity, the wormhole
solutions possess the Yang-Mills-Higgs hair where families of hairy wormholes
solutions emerge from the Ellis wormhole when the gravitational coupling
constant increases. In contrast to the BPS limit, the properties of wormholes
change drastically when the gravitational strength approaches a critical value
for a fixed $\lambda$. The hairy wormholes possess two types of double throat
configurations. The first type is wormholes that develop the double throat
when the gravitational strength almost approaches the critical value for lower
$\lambda$, whereas the second type is wormholes that exhibit the double throat
for a certain range of gravitational strength for higher $\lambda$. These two
types of double throat configurations can coexist for a certain range of
$\lambda$ where it is a transition process for first type double throat
disappears gradually and the second type double throat becomes dominant.

##  I Introduction

The solitonic magnetic monopole with finite energy is an inevitable outcome of
the non-Abelian SU(2) Yang-Mills-Higgs (YMH) theory upon spontaneously broken
by the Higgs field to a residual symmetry U(1) Rossi (1982); Shnir (2005). The
non-Abelian magnetic monopole discovered by ’t Hooft and Polyakov is a three-
dimensional topological soliton t Hooft (1974); Polyakov (1974), whereby its
topological charge corresponds to the magnetic charge. Exact monopole Prasad
and Sommerfield (1975) and axially symmetric multimonopole solutions Rebbi and
Rossi (1980); Forgács et al. (1981) are found in the Bogomol’nyi-Prasad-
Sommerfield (BPS) limit Bogomol’Nyi (1976) where the Higgs potential vanishes,
and the Higgs field becomes massless. These BPS solutions satisfy the first-
order Bogomol’nyi equation and their mass saturated at the lower Bogomol’nyi
bound Bogomol’Nyi (1976). Beyond the BPS limit, when the Higgs potential is
non-vanishing and the Higgs field becomes massive, the non-BPS solutions of
monopole Bogomol’nyi and Marinov (1976); Bais and Primack (1976); Teh et al.
(2010) and multimonopole Kleihaus et al. (1998) no longer fulfil the
Bogomol’nyi equation but only can be numerically obtained by solving the
second-order Euler-Lagrange equations of YMH theory. Likewise, the monopole-
antimonopole pair Kleihaus and Kunz (1999), vortex ring and monopole-
antimonopole chain Kleihaus et al. (2003); Kleihaus et al. (2004a) also do not
satisfy the Bogomol’nyi equation and they are the non-BPS, saddle-points
solutions in flat space which possess finite energy.

The coupling of the YMH model with gravity or known as Einstein-Yang-Mills-
Higgs (EYMH) theory gives rise to a branch of gravitating monopole solutions
Lee et al. (1992); Breitenlohner et al. (1992, 1995); Brihaye et al. (1998);
Lue and Weinberg (1999); Brihaye et al. (1999); Brihaye et al. (2000a, b);
Brihaye and Hartmann (2002). The mass of the gravitating monopole strictly
decreases as the gravitational strength increases. The solutions of
gravitating monopole, which includes the static axially symmetric gravitating
monopole solutions Hartmann et al. (2001a, b) end up as an extremal Reissner-
Nordstrom black hole when the gravitational strength reaches a critical value.
Note that the black holes in the EYMH system are dubbed as the “black hole
within monopole” Lee et al. (1992); Brihaye et al. (1998) because they possess
the non-abelian gauge field outside the event horizon, which is a
counterexample to the “no hair” conjecture for black holes. Similar to the YMH
theory, the EYMH theory also reported various gravitating soliton solutions
such as gravitating monopole and antimonopole pairs Kleihaus and Kunz (2000a,
b); Hartmann et al. (2001b); Van der Bij and Radu (2002); Paturyan and
Tchrakian (2004); Paturyan et al. (2005); Kleihaus et al. (2004b) and vortex
rings Kleihaus et al. (2005a, b).

Besides, some wormhole-like structures in EYMH theory have been reported in
Hajicek (1983a, b); Degen (1987). Recently the authors have constructed
numerically the solutions of symmetric wormholes in the EYMH system for the
BPS limit Chew and Lim (2020). The throat of the wormholes is supported by a
phantom field which can violate the null energy condition in order to prevent
the collapse of the throat for the construction of traversable wormholes in
general relativity. A phantom field is a real-valued scalar field that has an
opposite sign for the kinetic term in the Lagrangian. It has been used to
construct the classic example of traversable wormhole, which is known as the
Ellis wormhole Ellis (1973, 1979); Bronnikov (1973). However, the Ellis
wormhole possesses the unstable radial modes Gonzalez et al. (2009a, b); Torii
and Shinkai (2013). The Ellis wormhole has been generalized to the higher-
dimensional case Torii and Shinkai (2013), the slowly rotating case with
perturbative method Kashargin and Sushkov (2008a, b), the rapidly rotating
case in four dimensions Kleihaus and Kunz (2014); Chew et al. (2016) and five
dimensions with equal angular momenta Dzhunushaliev et al. (2013a), in the
modified gravity, e.g., the scalar-tensor theory Chew et al. (2018) and f(R)
gravity Karakasis et al. (2021). Recently, the Ellis wormhole has been
considered in the bouncing universe Huang et al. (2020) and generalized in
asymptotic anti-de Sitter Blázquez-Salcedo et al. (2020).

In the BPS limit, the EYMH wormholes admit the Ellis wormhole as the trivial
solution, which is analogous to Schwarzschild black hole when the gauge fields
vanish Chew and Lim (2020). Analogous to particle-like solutions of EYMH, the
corresponding wormholes possess a probe limit where only the gauge field is
present in the background solution of Ellis wormhole when the gravitational
field strength vanishes. When the gravity present, a branch of hairy wormholes
emerge from the Ellis wormhole where the wormholes gain the mass and its
phantom scalar charge increases. The wormholes possess a double-throat
configuration when the gravitational field strength exceeds a value.

In addition, another non-abelian wormhole in Einstein-Yang-Mills (EYM) and
phantom field system has been obtained numerically Hauser et al. (2014).
Similarly, these hairy wormholes solutions also possess a probe limit without
the presence of gravity. These hairy wormholes solutions possess a sequence of
solutions, which are labelled by the node number $k$ of the gauge field
function. They are analogous to the Bartnik-McKinnon solution, which is the
regular and spherically symmetric solutions of the EYM system Bartnik and
McKinnon (1988).

On the other hand, since wormholes can be a candidate for the black hole
mimicker, then several astrophysical signatures of wormholes have been
proposed to search for their existence in near future, for example, the shadow
Nedkova et al. (2013); Gyulchev et al. (2018); Amir et al. (2019); Narzilloev
et al. (2021); Benavides-Gallego et al. (2021); Jusufi et al. (2021); Bambi
and Stojkovic (2021), the gravitational lensing Abe (2010); Toki et al.
(2011); Takahashi and Asada (2013); Cramer et al. (1995); Perlick (2004);
Tsukamoto et al. (2012); Bambi (2013), the accretion disk around the wormhole
Zhou et al. (2016); Deligianni et al. (2021a, b), and the ringdown phase in
the emission of gravitational waves Blázquez-Salcedo et al. (2018).

Since EYMH wormholes for the BPS limit in our previous work Chew and Lim
(2020) give rise to new and interesting phenomena due to the presence of the
non-Abelian field. In this paper, our motivation is to continue our
investigation by numerically obtaining the symmetric wormhole solutions in
EYMH with the presence of Higgs potential and analyze their properties. Thus
our paper is organized as follows. In section II, we briefly introduce the
EYMH theory and present the equation of motions. Subsequently, we introduce
the geometrical properties of the wormhole. We then obtain the ordinary
differential equations (ODEs) from the equation of motions and discuss the
global charges of the wormholes. In section III, we exhibit and discuss our
numerical results. In section IV, we conclude our research work and discuss
the possible outlook from this present work.

##  II Theoretical Framework

###  II.1 Einstein-Yang-Mills-Higgs Theory

The Einstein gravity couples with a phantom field $\psi$ and a gauge field
$A_{\mu}$ in SU(2) YMH theory in the Einstein-Hilbert action,

$S=\int d^{4}x\sqrt{-g}\left[\frac{R}{16\pi
G}+\mathscr{L}_{\text{ph}}+\mathscr{L}_{\text{YMH}}\right]\,,$ (1)

where the Lagrangian of phantom field and YMH Hartmann et al. (2001b) are
respectively, given by

$\mathscr{L}_{\text{ph}}=\frac{1}{2}\partial_{\mu}\psi\partial^{\mu}\psi\,,\quad\mathscr{L}_{\text{YMH}}=-\frac{1}{2}\text{Tr}\left(F_{\mu\nu}F^{\mu\nu}\right)-\frac{1}{4}\text{Tr}\left(D_{\mu}\Phi
D^{\mu}\Phi\right)-\frac{\lambda}{8}\text{Tr}(\Phi^{2}-\upsilon^{2})^{2}\,,$
(2)

where $\lambda$ is the Higgs field self-interaction constant and $\upsilon$ is
the vacuum expectation value of the Higgs field. The covariant derivative of
the Higgs field and the gauge field strength tensor are given respectively by

$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+i\left[A_{\mu},A_{\nu}\right]\,,\quad
D_{\mu}\Phi=\partial_{\mu}\Phi+i\left[A_{\mu},\Phi\right]\,,$ (3)

where $A_{\mu}=\frac{1}{2}\tau^{a}A_{\mu}^{a}$ and $\Phi=\phi^{a}\tau^{a}$
with $\tau^{a}$ is the Pauli matrices. Since we only construct spherically
symmetric wormholes, then we also employ the spherically symmetric Ansatz in a
purely magnetic gauge field $(A_{t}=0)$ for the gauge and Higgs field Brihaye
et al. (1998),

$A_{\mu}dx^{\mu}=\frac{1-K(\eta)}{2}\left(\tau_{\varphi}d\theta-\tau_{\theta}\sin\theta
d\varphi\right)\,,\quad\Phi=H(\eta)\tau_{\eta}\,.$ (4)

The action is varied with respect to the metric $g_{\mu\nu}$, which yields the
Einstein equation,

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\beta\left(T_{\mu\nu}^{\text{ph}}+T_{\mu\nu}^{\text{YMH}}\right)\,,$
(5)

where $\beta=8\pi G$, the stress-energy tensor for phantom field
$T_{\mu\nu}^{\text{ph}}$ and YMH $T_{\mu\nu}^{\text{YMH}}$ are respectively,
given by

$\displaystyle T_{\mu\nu}^{\text{ph}}$
$\displaystyle=\frac{1}{2}g_{\mu\nu}\partial_{\alpha}\psi\partial^{\alpha}\psi-\partial_{\mu}\psi\partial_{\nu}\psi\,,$
(6) $\displaystyle T_{\mu\nu}^{\text{YMH}}$
$\displaystyle=\text{Tr}\left(\frac{1}{2}D_{\mu}\Phi
D_{\nu}\Phi-\frac{1}{4}g_{\mu\nu}D_{\alpha}\Phi
D^{\alpha}\Phi\right)+2\text{Tr}\left(g^{\alpha\beta}F_{\mu\alpha}F_{\nu\beta}-\frac{1}{4}g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\right)-\frac{\lambda}{8}g_{\mu\nu}\text{Tr}(\Phi^{2}-\upsilon^{2})^{2}\,.$
(7)

The equations of motion for the matter fields are

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}\partial^{\mu}\psi\right)=0\,,\quad
D_{\mu}F^{\mu\nu}=\frac{i}{4}\left[\Phi,D^{\nu}\Phi\right]\,,\quad
D_{\mu}D^{\mu}\Phi=\lambda(\Phi^{2}-\eta^{2})\Phi\,.$ (8)

###  II.2 The Geometrical Structure of Wormhole

A globally regular wormhole spacetime can be constructed by employing the
quasi-isotropic line element,

$ds^{2}=-F_{0}(\eta)dt^{2}+F_{1}(\eta)\left[d\eta^{2}+h(\eta)\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\right]\,,$ (9)

where $h(\eta)=\eta^{2}+\eta_{0}^{2}$ with $\eta_{0}$ as the throat parameter.
The wormhole spacetime possesses two asymptotically flat regions in the limit
$\eta\rightarrow\pm\infty$. Here we define the function $R(\eta)^{2}$ as the
circumferential radius of the wormhole,

$R(\eta)^{2}=F_{1}h\,.$ (10)

Note that $R(\eta)$ should not contain zero for a globally regular wormhole
solution. When $R$ contains a local minimum, which is known as the throat of
wormhole $\eta_{\text{th}}$, then the wormhole possesses a minimal surface
area at the throat, $A_{\text{th}}=4\pi R(\eta_{\text{th}})^{2}$. However, if
$R$ contains a local maximum, then it is an equator of the wormhole
$\eta_{\text{eq}}$, which corresponds to the maximal surface area of the
wormhole, $A_{\text{eq}}=4\pi R(\eta_{\text{eq}})^{2}$. The equator of the
wormhole is normally sandwiched between two throats.

The wormholes that possess either a throat or an equator can be determined by
the following conditions, respectively

$\displaystyle 1.\quad R^{\prime}(\eta_{\text{th}})=0\,,\ $ and
$\displaystyle\ R^{\prime\prime}(\eta_{\text{th}})>0\,,$ (11) $\displaystyle
2.\quad R^{\prime}(\eta_{\text{eq}})=0\,,\ $ and $\displaystyle\
R^{\prime\prime}(\eta_{\text{eq}})<0\,.$ (12)

When $R^{\prime}(\eta_{\text{crit}})=R^{\prime\prime}(\eta_{\text{crit}})=0$,
the geometry of wormhole is in a transition state because the circumferential
radius forms a turning point at some value of the radial coordinate
$\eta_{\text{crit}}$, such that the double throat and the equator can
simultaneously exist, this also implies that there is a transition can occur
from the single throat configuration to the double throat configuration
Dzhunushaliev et al. (2014a); Hoffmann et al. (2017).

In this paper, we only consider the metric functions symmetric w.r.t. the
coordinate $\eta=0$, thus the circumferential radius of the wormhole at
$\eta=0$ is assumed to be either a throat or an equator, which implies $R$
should have an extremum at $\eta=0$ by demanding

$R^{\prime}(0)=0\quad\Rightarrow\quad\frac{\left(hF^{\prime}_{1}+2\eta
F_{1}\right)}{2R}\Bigg{|}_{\eta=0}=0\,,$ (13)

where we have to set $F^{\prime}_{1}(0)=0$. In particular, if a wormhole only
contains a single throat, then the throat must be located at $\eta=0$ with the
minimal surface area $A_{\text{th}}=4\pi R(0)^{2}=4\pi F_{1}(0)\eta_{0}^{2}$.

###  II.3 Ordinary Differential Equations (ODEs)

A set of second-order and nonlinear ODEs is obtained for the metric functions
and gauge fields by substituting Eqs. (9) and (4) into the Einstein equation
Eq. (5) and equations of motion for the gauge fields Eq. (8),

$\displaystyle
F^{\prime\prime}_{1}+\frac{2\eta}{h}F^{\prime}_{1}-\frac{3F^{\prime
2}_{1}}{4F_{1}}+\frac{\eta_{0}^{2}F_{1}}{h^{2}}$
$\displaystyle\qquad\qquad\qquad=\beta\frac{F_{1}}{2}\psi^{\prime
2}-\beta\frac{(K^{2}-1)^{2}+2hK^{\prime 2}+h^{2}F_{1}H^{\prime
2}+2hF_{1}H^{2}K^{2}}{2h^{2}}-\beta\frac{\lambda}{4}F_{1}^{2}\left(H^{2}-\upsilon^{2}\right)^{2}\,,$
(14)
$\displaystyle\left(\frac{F^{\prime}_{1}}{2F_{1}}+\frac{\eta}{h}\right)\frac{F^{\prime}_{0}}{F_{0}}+\frac{F^{\prime
2}_{1}}{4F^{2}_{1}}+\frac{\eta}{hF_{1}}F^{\prime}_{1}-\frac{\eta_{0}^{2}}{h^{2}}$
$\displaystyle\qquad\qquad\qquad=-\frac{\beta}{2}\psi^{\prime
2}+\beta\frac{-(K^{2}-1)^{2}+2hK^{\prime 2}+h^{2}F_{1}H^{\prime
2}-2hF_{1}H^{2}K^{2}}{2h^{2}F_{1}}-\beta\frac{\lambda}{4}F_{1}(H^{2}-\upsilon^{2})^{2}\,,$
(15) $\displaystyle
F^{\prime\prime}_{0}+\left(-\frac{F^{\prime}_{0}}{2F_{0}}+\frac{\eta}{h}\right)F^{\prime}_{0}+\frac{F_{0}}{F_{1}}F^{\prime\prime}_{1}+\left(-\frac{F^{\prime}_{1}}{F_{1}}+\frac{\eta}{h}\right)\frac{F_{0}F^{\prime}_{1}}{F_{1}}+\frac{2F_{0}\eta_{0}^{2}}{h^{2}}$
$\displaystyle\qquad\qquad\qquad=\beta F_{0}\psi^{\prime
2}-F_{0}\left[\frac{-(K^{2}-1)^{2}+h^{2}F_{1}H^{\prime
2}}{h^{2}F_{1}}\right]-\beta\frac{\lambda}{2}F_{0}F_{1}(H^{2}-\upsilon^{2})^{2}\,.\,$
(16) $\displaystyle K^{\prime\prime}$
$\displaystyle+\frac{1}{2}\left(\frac{F^{\prime}_{0}}{F_{0}}-\frac{F^{\prime}_{1}}{F_{1}}\right)K^{\prime}=\frac{K(K^{2}-1+hF_{1}H^{2})}{h}\,,$
(17) $\displaystyle H^{\prime\prime}$
$\displaystyle+\frac{1}{2}\left(\frac{F^{\prime}_{0}}{F_{0}}+\frac{F^{\prime}_{1}}{F_{1}}+\frac{4\eta}{h}\right)H^{\prime}-\frac{2K^{2}}{h}H=\lambda
F_{1}(H^{2}-\upsilon^{2})H\,.$ (18)

where the prime denotes the derivative of the functions w.r.t. the radial
coordinate $\eta$. The equation of motion Eq. (8) for the phantom field yields
a first-order integral, $\psi^{\prime}=D/(h\sqrt{F_{0}F_{1}})$ where $D$ is
the scalar charge of the phantom field. Thus, the term $\psi^{\prime 2}$ in
Eqs. (40)-(42) can be replaced by $\psi^{\prime}=D/(h\sqrt{F_{0}F_{1}})$. We
subtract Eqs. (40) and (42) with Eq. (41) to obtain the final set of ODEs,

$\displaystyle F^{\prime\prime}_{0}$
$\displaystyle=\frac{F^{\prime}_{0}}{2}\left(\frac{F^{\prime}_{0}}{F_{0}}-\frac{F^{\prime}_{1}}{F_{1}}-\frac{4\eta}{h}\right)+\beta
F_{0}\frac{(K^{2}-1)^{2}+2hK^{\prime 2}}{h^{2}F_{1}}-\frac{1}{2}\beta
F_{0}F_{1}\lambda(H^{2}-\upsilon^{2})^{2}\,,$ (19) $\displaystyle
F^{\prime\prime}_{1}$ $\displaystyle=\frac{F^{\prime
2}_{1}}{2F_{1}}-\frac{3\eta
F^{\prime}_{1}}{h}-\left(\frac{F^{\prime}_{1}}{2}+\frac{F_{1}}{h}\right)\frac{F^{\prime}_{0}}{F_{0}}-\frac{\beta\left[(K^{2}-1)^{2}+2hF_{1}H^{2}K^{2}\right]}{h^{2}}-\frac{\beta\lambda}{2}F_{1}^{2}(H^{2}-\upsilon^{2})^{2}\,,$
(20) $\displaystyle K^{\prime\prime}$
$\displaystyle=\left(\frac{F^{\prime}_{1}}{F_{1}}-\frac{F^{\prime}_{0}}{F_{0}}\right)\frac{K^{\prime}}{2}+\frac{K\left(K^{2}-1+hF_{1}H^{2}\right)}{h}\,,$
(21) $\displaystyle H^{\prime\prime}$
$\displaystyle=-\frac{1}{2}\left(\frac{F^{\prime}_{0}}{F_{0}}+\frac{F^{\prime}_{1}}{F_{1}}+\frac{4\eta}{h}\right)H^{\prime}+\frac{2HK^{2}}{h}+\lambda
F_{1}(H^{2}-\upsilon^{2})H\,,$ (22)

with Eq. (41) is expressed as

$\displaystyle D^{2}$
$\displaystyle=\frac{2h^{2}F_{0}F_{1}}{\beta}\left[-\left(\frac{F^{\prime}_{1}}{2F_{1}}+\frac{\eta}{h}\right)\frac{F^{\prime}_{0}}{F_{0}}-\frac{F^{\prime
2}_{1}}{4F^{2}_{1}}-\frac{\eta}{hF_{1}}F^{\prime}_{1}+\frac{\eta_{0}^{2}}{h^{2}}\right.$
$\displaystyle\qquad\qquad\left.+\beta\frac{-(K^{2}-1)^{2}+2hK^{\prime
2}+h^{2}F_{1}H^{\prime
2}-2hF_{1}H^{2}K^{2}}{2h^{2}F_{1}}-\beta\frac{\lambda}{4}F_{1}(H^{2}-\upsilon^{2})^{2}\right]\,.$
(23)

to monitor the quality of the numerical solutions with the condition
$D^{2}=const$.

Since we only consider the wormhole solutions with the metric functions
symmetric w.r.t. $\eta=0$, thus we solve Eqs. (19)-(22) numerically from
$\eta=0$ to the infinity by Colsys and Matlab package bvp4c. Both packages
solve the boundary value problems for systems of nonlinear coupled ODEs,
Colsys is based on the Newton-Raphson method with adaptive mesh refinement and
error estimation for the solutions Ascher et al. (1979) whereas bvp4c
implements the three-stage Lobatto IIIa collocation formula Kierzenka and
Shampine (2001). In the numerics, we impose the eight boundary conditions at
$\eta=0$ and $\eta=\infty$. First, we let the the metric functions possess the
extremum at $\eta=0$ by requiring the first order derivative of the metric
functions vanish at the throat, $F^{\prime}_{0}(0)=F^{\prime}_{1}(0)=0$. The
metric functions approach Minkowski spacetime at the infinity,
$F_{0}(\infty)=F_{1}(\infty)=1$. We impose the following boundary conditions
$K(0)=1$, $H(0)=0$, $K(\infty)=0$, $H(\infty)=1$ for the gauge fields.
Furthermore, we compactify the radial coordinate $\eta$ by
$\eta=\eta_{0}\tan\left(\pi x/2\right)$ in the numerics. We also scale some
parameters by the throat parameter $\eta_{0}$,

$\eta\rightarrow\eta_{0}\eta\,,\quad
h\rightarrow\eta_{0}^{2}h\,,\quad\beta\rightarrow\eta^{2}_{0}\beta\,,\quad
H\rightarrow\frac{H}{\eta_{0}}\,,\quad\upsilon\rightarrow\frac{\upsilon}{\eta_{0}}\,.$
(24)

Besides, the second-order derivative of $R$ at $\eta=0$ is given by

$R^{\prime\prime}(0)=\frac{2F_{1}+hF^{\prime\prime}_{1}}{2R}\Bigg{|}_{\eta=0}=\frac{2F_{1}-\frac{1}{2}F^{2}_{1}\beta\lambda\eta^{2}_{0}\upsilon^{4}}{2R}\,,$
(25)

where we have used $F^{\prime\prime}_{1}(0)$ from the ODEs. In the BPS limit,
$R(0)$ always remains as the throat since $R^{\prime\prime}(0)>0$.

###  II.4 Global Charges

(a)  (b)  (c)

Figure 1: (a) The solutions of Hooft-Polyakov monopole: $K(x)$ (solid line)
and $H(x)$ (dash-dotted line) for several values of $\lambda$ in the
compactified coordinate, $x=R/(1+R)$. (b) The mass of Hooft-Polyakov monopole
as a function of $\lambda$. (c) The gauge fields $H(x)$ (dash-dotted line) and
$K(x)$ (solid line) in the probe limit of hairy wormholes for several values
of $\lambda$ in the compactified coordinate $x$.

The mass of wormhole can be read off directly from the asymptotic expansion of
the metric at $\eta\rightarrow\infty$,

$F_{0}\rightarrow 1-\frac{2GM}{\eta}=1-\frac{2\mu}{\bar{\eta}}\,,$ (26)

where $\mu$ is the mass parameter, which is given by

$\mu=\frac{\beta M}{8\pi\eta_{0}}\,.$ (27)

One also can obtain the expression for the mass via the Komar integral Hauser
et al. (2014),

$M=M_{\text{th}}+\frac{1}{4\pi
G}\int_{\Sigma}R_{\mu\nu}\xi^{\mu}n^{\nu}dV=M_{\text{th}}+\frac{1}{4\pi G}\int
R^{0}_{0}\sqrt{-g}d^{3}x\,,$ (28)

where $\Sigma$ is the spacelike hypersurface ($0\leq\eta<\infty$), $n^{\nu}$
is a normal vector on $\Sigma$, $\xi^{\mu}=(1,0,0,0)$ is a timelike Killing
vector, $dV$ is the volume element on $\Sigma$. The term $M_{\text{th}}$ is
the contribution of the throat to the mass,

$M_{\text{th}}=\frac{\kappa A_{\text{th}}}{4\pi G}\,,$ (29)

where $\kappa$ is the surface gravity at the throat, which is given by Chew
and Lim (2020)

$\kappa=\frac{F^{\prime}_{0}}{2\sqrt{F_{0}F_{1}}}\,.$ (30)

Eq. (30) shows that $\kappa$ and $M_{\text{th}}$ vanish for symmetric
wormholes with only a single throat but remains finite for wormholes with
double throat configuration.

The charge of the phantom field is given by $D^{2}$ and the magnetic charge
for the non-abelian gauge fields is given by Corichi et al. (2000); Ashtekar
et al. (2001)

$\mathscr{P}^{\text{YMH}}=\frac{1}{4\pi}\oint\sqrt{\sum_{i}\left(F^{i}_{\theta\varphi}\right)^{2}}d\theta
d\varphi=|P|\,,$ (31)

where the integral is evaluated at the spatial infinity, yielding $P=0$ for
the hairy wormholes Hauser et al. (2014).

##  III Results and Discussions

###  III.1 Probe Limit

(a)  (b)

Figure 2: The global charges of wormhole solutions for several values of
$\lambda$: (a) The scaled mass $\mu/\alpha^{2}$ versus the scaled
gravitational coupling constant $\alpha$. (b) The logarithmic of scaled scalar
charge $\log_{10}(\alpha^{2}D^{2})$ versus the scaled gravitational coupling
constant $\alpha$. Note that $\beta=2\alpha^{2}$.

Here we fix the vacuum expectation value $\upsilon=1$ in all numerics
computation. In this part, we discuss the probe limit of hairy wormholes. But
first, let us recall that YMH theory possesses a probe limit which is known as
Hooft-Polyakov monopole. In the BPS limit, it possesses an exact solution,

$K(R)=\frac{R}{\sinh(R)}\,,\quad H(R)=\text{coth}(R)-\frac{1}{R}\,.$ (32)

However, the solutions of Hooft-Polyakov monopole beyond the BPS limit can
only be obtained numerically which are shown in Fig. 1(a). The mass of Hooft-
Polyakov monopole as shown in Fig. 1(b) increases monotonically from the unity
which corresponds to the BPS limit.

When the gravity is switched off $(\beta=0)$, the metric Eq. (9) in the probe
limit of hairy wormholes is the massless Ellis wormhole
$(F_{0}(\eta)=F_{1}(\eta)=1)$ in which the spacetime is symmetric and has two
identical asymptotically flat regions. The phantom field $\psi$ is given
exactly by

$\psi=\frac{D}{\eta_{0}}\left[\arctan\left(\frac{\eta}{\eta_{0}}\right)-\frac{\pi}{2}\right]\,.$
(33)

Hence, we obtain the pure YMH equations in the background of the Ellis
wormhole because the YMH field doesn’t contribute to the Einstein equation,

$\displaystyle K^{\prime\prime}$
$\displaystyle=\frac{K(K^{2}-1+hH^{2})}{h}\,,$ (34) $\displaystyle
H^{\prime\prime}$
$\displaystyle=-\frac{2\eta}{h}H^{\prime}+\frac{2K^{2}}{h}H+\lambda
F_{1}(H^{2}-1)H\,.$ (35)

The above ODEs are solved numerically and shown in Fig. 1(c).

(a)  (b)  (c)  (d)

Figure 3: The gauge field $K(x)$ (solid line) and Higgs field $H(x)$ (dash-
dotted line) in the compactified coordinate $x$ for the wormhole solutions by
varying $\beta$ with fixed $\lambda$: (a) $\lambda=0.2$; (b) $\lambda=1$; (c)
$\lambda=30$ and (d) $\lambda=70$.

(a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)

Figure 4: The geometrical properties of the wormhole solutions. (a) and (b)
are the location of throats $x_{\text{th}}$ and equator $x_{\text{eq}}$ in the
compactified coordinate $x$ versus $\beta$, respectively for $\lambda=0.2$
(Type I), $\lambda=30$ (Type II); (c) and (d) are the second derivative of the
circumferential radius of throats and equator versus $\beta$, respectively for
$\lambda=0.2,30$; (e) and (f) are the size of throats and equator versus
$\beta$, respectively for $\lambda=0.2,30$; (g) and (h) are the
circumferential radius $R$ in the radial coordinate $\eta$ with several values
of $\beta$, respectively for $\lambda=0.2,30$;

(a)  (b)

Figure 5: The isometric embedding for the wormholes solutions: (a) Type I
double throat $\lambda=0.2$ and $\beta=6.15$; (b) Type II double throat
$\lambda=30$ and $\beta=0.4$.

(a)  (b)

Figure 6: (a) The location of throats and equator of hybrid type I + II double
throat configuration for $\lambda=2$. (b) The surface gravity $\kappa$ at the
throat of wormholes for several values of $\lambda$.

(a)  (b)  (c)  (d)

Figure 7: The metric function $F_{1}(x)$ in the compactified coordinate $x$
for the wormhole solutions by varying $\beta$ with fixed $\lambda$: (a)
$\lambda=0.2$; (b) $\lambda=1$; (c) $\lambda=30$ and (d) $\lambda=70$.

(a)  (b)  (c)  (d)  (e)  (f)

Figure 8: The violation of null energy condition in the compactified
coordinate $x$ for the wormhole solutions by varying $\beta$ with fixed
$\lambda$: (a) $\lambda=0.2$; (b) $\lambda=1$; (c) $\lambda=30$; (d)
$\lambda=70$; (e) $\lambda=2$ with the inset of zoom in for small postive
values of NEC at $x=0$; (f) $\lambda=2$ with the inset of zoom in for the peak
near $x=1$.

###  III.2 With Backreaction

Here we exhibit our numerical results for hairy wormholes beyond the BPS limit
($\lambda\neq 0$) by fixing the Higgs self-interaction parameter $\lambda$ for
the range $[0,100]$ while varying the gravitational coupling constant $\beta$.
The hairy wormholes solutions should exist for all positive real values of
$\lambda$. The hairy wormholes in the BPS limit could take any real positive
values of $\beta$ Chew and Lim (2020), this is in contrast to the YM wormhole
where the limiting configuration is the extremal Reissner-Nordstrom black hole
for higher nodes Hauser et al. (2014). Recall that when the gauge fields don’t
present in the probe limit, the Ellis wormhole has the analogue of
Schwarzschild solution for the black holes Hauser et al. (2014); Chew and Lim
(2020). Recall also that the Ellis wormhole is massless, the circumferential
radius of throat and scalar charge are unity. With fixed $\lambda$, when we
increase $\beta$ from zero, the families of solutions for hairy wormholes
emerge from the Ellis wormhole, thus the properties of these hairy wormholes
behave differently with the Ellis wormholes when $\beta$ increases.

Fig. 2(a) shows that the hairy wormholes gain the mass when scaled
gravitational strength $\alpha$ ($\beta=2\alpha^{2}$) increases from zero for
several fixed values of $\lambda$. In the BPS limit, the mass of hairy
wormholes increase monotonically as $\alpha$ increases and the mass is always
positive. For a small value of the Higgs self-interaction $\lambda$, the mass
of wormholes increases to a maximum value when $\alpha$ increases from zero.
However, their mass subsequently drops very sharply to the negative value as
$\alpha$ approaches a critical value. For a large value of $\lambda$, the mass
of wormholes decreases to the negative value monotonically from zero as
$\beta$ increases from zero and then decreases very sharply when $\beta$
approaches the critical value. Thus, we see that the mass of wormholes for the
non-BPS limit and BPS limit are quite different.

Fig. 2(b) shows that the scaled scalar charge $\alpha^{2}D^{2}$ of the phantom
field increases monotonically from unity as scaled gravitational strength
$\alpha$ increases from zero for several fixed values of $\lambda$. Then
$\alpha^{2}D^{2}$ increases very sharply as $\alpha$ approaches a critical
value. However, in the BPS limit $\alpha^{2}D^{2}$ increases monotonically as
$\alpha$ increases. Analogous to the mass, the scaled phantom charge of
wormholes for beyond BPS limit and BPS limit are also quite different. Note
that in Figs. 2(a) and (b) the domain of existence for scaled gravitational
strength $\alpha$ shrinks for the wormholes beyond BPS limit when $\lambda$
increases, hence the global charges can change dramatically even for a smaller
value of $\alpha$ for very large values of $\lambda$. Furthermore, the
precision of wormhole solutions drop significantly when the gravitational
strength approaches the critical value, thus we could not generate the
solutions and study their properties beyond that critical value.

Figs. 3(a)-(d) exhibit the profiles of gauge field $K(x)$ and Higgs field
$H(x)$ in the compactified coordinate $x$ for $\lambda=0.2,1,30,70$,
respectively with several values of $\beta$. The two fields behave
monotonically. The Higgs field $H$ increases monotonically from zero to reach
its asymptotic value $H=1$ at the infinity while the gauge field $K$ decreases
monotonically from unity to reach its asymptotic value $K=0$ at the infinity.
As $\beta$ approaches the critical value, $H$ increases very sharply and
becomes very steep (with a very large gradient) near the asymptotically flat
region, this implies that $H$ might diverge at the asymptotically flat region.
However, $K$ approaches unity when $\beta$ reaches the critical value.
Nevertheless in the BPS limit, the gauge field $K$ decays faster while the
profile of gauge field $H$ doesn’t vary too much when $\beta$ increases Chew
and Lim (2020).

Let’s turn our discussion to the geometry of the hairy wormholes with non-
vanishing $\lambda$. In the BPS limit ($\lambda=0$) when $\beta$ exceeds a
value the hairy wormholes possess a double throat configuration where an
equator is sandwiched by a throat at the radial coordinate $\eta=0$ and
another throat which is near the asymptotic flat region Chew and Lim (2020).
When $\lambda>0$, the hairy wormholes possess two different types of double
throat configurations, which are type I for the range $0<\lambda<3.0283$ and
type II for the range $1.8814\leq\lambda\leq 100$ as shown in Fig. 4 (a) and
(b), respectively. Interestingly, type I and II double throat configurations
can coexist in the range of $1.8814\leq\lambda<3.0283$ where we will discuss
this in detail in the later paragraphs.

We exhibit the wormhole solutions with $\lambda=0.2$ for type I double throat
configuration as depicted in Figs. 4(a), (c), (e) and (g). In Fig. 4(a), the
wormholes only have a single throat at $x=0$ represented by the yellow square
on the horizontal axis when $\beta$ increases from zero. When $\beta$
approaches a critical value, the wormholes possess a double throat in which
the throat at $x=0$ remains as the throat, but another throat $x_{\text{th}}$
(purple curve) simultaneously exists near the asymptotically flat region as
shown in the inset of Fig. 4(a). We observe that there is an equator
$x_{\text{eq}}$ (green curve) in the inset of Fig. 4(a) lies in between that
two throats as shown inside the inset in Fig. 4(a). Note that at the beginning
of the formation of the double throat, $x_{\text{th}}$ and $x_{\text{eq}}$ are
very close to each other, they move away from each other and move toward the
spatial infinity $x=1$ as $\beta$ increases. When $\beta$ almost reaches the
critical value, both $x_{\text{th}}$ and $x_{\text{eq}}$ approach each other
again. In Fig. 4(c), as $\beta$ increases from zero before approaches to the
critical value, the second derivative of circumferential radius
$R^{\prime\prime}(0)$ in the compactified coordinate $x$ for the throat at
$x=0$ decreases from its maximum value to a very small value which is close to
zero. When type I double throat appears, $R^{\prime\prime}(x_{\text{eq}})$
(green curve) is always negative and decreases very sharply for the equator
$x_{\text{eq}}$ and $R^{\prime\prime}(x_{\text{th}})$ (purple curve) is always
positive and increases very sharply for another throat $x_{\text{th}}$.

In Fig. 4(e), the size of the throat $R(0)$ (yellow curve) at radial
coordinate $x=0$ increases from unity to a maximum value and then decreases to
nearly zero as $\beta$ increases from zero to the critical value. In the inset
of Fig. 4(e), the size of another throat $R(x_{\text{th}})$ (purple curve) and
the size of equator $R(x_{\text{eq}})$ (green curve) are very close to each
other at the beginning of the formation of the double throat because they are
very close to each other. Then $R(x_{\text{th}})$ and $R(x_{\text{eq}})$
decrease and their difference become larger when $\beta$ approaches the
critical value.

In Fig. 4(g), the wormhole solutions with $\beta=0.5,1.0,2.0,3.0$ clearly only
possess a single throat at $x=0$. Although the solutions with
$\beta=5.0,6.0,6.1$ still possess a single throat, their circumferential radii
are being deformed for $0\leq\log_{10}(\eta)\leq 2$ which is a process in
developing the double throat configuration. Hence, the solutions with
$\beta=6.15,6.194$ (dashed line) which approach to critical value possess the
double throat configuration because the circumferential radius $R$ contains
two local minima and one local maximum. The isometric embedding of the type I
double throat configuration for $\lambda=0.2$ can be visualized in Fig. 5(a).

For type II double throat configuration, we exhibit the wormhole solutions
with $\lambda=30$ as depicted in Figs. 4(b), (d), (f) and (h). As shown in
Fig. 4(b), the wormholes only possess a single throat at $x=0$ for small
values of $\beta$, since $R^{\prime\prime}(0)>0$ (yellow curve) in Fig. 4(d).
As $\beta$ increases to a value in which $R^{\prime\prime}(0)$ decreases to
zero in Fig. 4(d), this implies that the wormholes are in a transition from a
single throat configuration to a double throat configuration. Hence, $x=0$
becomes an equator (green dot in Fig. 4(b)) with $R^{\prime\prime}(0)<0$
(green curve) in Fig. 4(d), and simultaneously another throat $x_{\text{th}}$
(purple curve) in Fig. 4(b) with $R^{\prime\prime}(x_{\text{th}})>0$ (purple
curve) in Fig. 4(d) is developed nearby $x=0$. As $\beta$ continues to
increase, the throat $x_{\text{th}}$ moves away from the equator $x=0$ to a
maximum distance and then moves back toward the equator $x=0$ at a certain
distance, while in Fig. 4(d) $R^{\prime\prime}(0)$ (green curve) decreases
from zero to a minimum value then increases again to zero and
$R^{\prime\prime}(x_{\text{th}})$ increases from zero to a maximum value then
decreases again to zero. Hence, both
$R^{\prime\prime}(0)=R^{\prime\prime}(x_{\text{th}})=0$, then this again
implies the wormholes are in another transition from a double throat
configuration to a single throat configuration where in Fig. 4(b) the throat
$(x_{\text{th}})$ disappears and the equator $x=0$ changes back as a throat
until $\beta$ reaches to the critical value. Moreover, we find that when
$\lambda$ increases, type II double throat appears earlier for a small value
of $\beta$.

In Fig. 4(f), the size of the throat $R(0)$ (yellow curve) at $x=0$ increases
from unity as $\beta$ increases from zero. At the beginning of the formation
of type II double throat configuration where the throat at $x=0$ becomes the
equator with equator size $R(0)$ (green curve) and another throat
$x_{\text{th}}$ with throat size $R(x_{\text{th}})$ (purple curve)
simultaneously appears nearby $x=0$, hence $R(0)=R(x_{\text{th}})$. As $\beta$
continues to increase, the equator size $R(0)$ increases to a maximum value
then decreases while $R(x_{\text{th}})$ also increases to a maximum value
where $x_{\text{th}}$ is farthest from $x=0$ and then decreases. Note that the
difference between $R(x_{\text{th}})$ between $R(0)$ is very small. When the
wormholes transform back from double throat configuration to the single throat
configuration, $R(0)$ and $R(x_{\text{th}})$ coincide again and then $R(0)$
decreases toward zero when $\beta$ reaches the critical value.

In Fig. 4(h), the wormhole solutions with $\beta=0.05,0.1$ only possess a
single throat at $x=0$ but $\beta=0.15,0.2,0.3,0.4,0.5$ (dashed line) possess
type II double throat configuration with one local minimum and one local
maximum. Although the wormholes with $\beta=0.8,1.0,1.45$ only possess a
single throat, we see that their circumferential radius has been slightly
deformed and eventually the circumferential radius satisfies the relation
$|R|\rightarrow|\eta|$ at infinity. The isometric embedding for type II double
throat configuration with $\lambda=30$ can be visualized in Fig. 5(b).

Previously we have mentioned that type I and II double throat configurations
can coexist for $1.8814\leq\lambda<3.0283$, here we denote the hybrid of type
I and type II as type I + II configuration and illustrate them by choosing
$\lambda=2$ in Fig. 6(a). In fact, this type I + II configuration serves as a
transition state where the type I double throat configuration can transform
into the type II double throat configuration gradually. Initially, the hairy
wormholes only possess type I double throat configuration for $\lambda<1.8814$
when $\beta$ approaches a critical value. When $\lambda=1.8814$, the type II
double throat configuration starts to appear in addition to type I for a
certain range of gravitational strength which is far lesser than that critical
value. The example $\lambda=2$ has shown that the type I + II double throat
configurations coexist when type II appears for $0.8806\leq\beta\leq 1.0074$
in the left inset of Fig. 6(a). As $\lambda$ increases from $1.8814$ to
$3.0283$, the type I double throat configuration is diminishing gradually
whereas the type II configuration is enlarged increasingly. We observe that in
this process, the range of $\beta$ for the appearance of type I double throat
configuration decreases and the locations of the new throat and equator are
approaching infinity. All these factors combined contribute to the
disappearance of type I double throat configuration. Eventually, when
$\lambda=3.0283$, the type I double throat configuration disappears completely
and the double throat configuration is now dominated by type II. Hence, the
varying of Higgs self-interaction value in the Higgs potential gives rise to
the transition of type I to type II via the intermediate state of type I + II
configuration.

The surface gravity $\kappa$ of wormholes is depicted in Fig. 6(b). Recall
that $\kappa$ vanishes when the wormholes only possess a single throat at
$x=0$. For type I double throat configuration, $\kappa$ only assumes a finite
value when $\beta$ approaches the critical value where the wormholes develop
another throat $x_{\text{th}}$ near the asymptotically flat region. As
$\lambda$ increases, $\kappa$ decreases very sharply. For type II double
throat configuration, $\kappa$ also only assumes finite value when double
throat configuration exists. Thus, starting from the beginning of the
formation of double throat configuration, $\kappa$ decreases from zero to a
minimum value and then increases to zero again where the wormholes transform
back to their original single throat configuration. Likewise, in the BPS
limit, the surface gravity of the single throat wormhole vanishes but the
double throat wormhole assumes a positive finite value.

Figs. 7(a)-(d) exhibit the metric function $F_{1}(x)$ in the compactified
coordinate $x$ for $\lambda=0.2,1,30,70$, respectively. When $\beta=0$,
$F_{1}(x)=1$ for Ellis wormholes. As $\beta$ increases, $F_{1}(0)$ assumes the
maximum value and $F_{1}(x)$ decreases monotonically to its asymptotic value
$F_{1}(1)=1$. After $\beta$ exceeds a value, $F_{1}$ develops another local
minimum which is near the asymptotically flat region while $F_{1}(0)$ still
retain as the local maximum. As $\beta$ increases and then approaches the
critical value, both local maximum $F_{1}(0)$ and local minimum decrease near
to zero and the function $F(x)$ increases very steeply from that local minimum
toward the asymptotic value $F(1)=1$. In the BPS limit Chew and Lim (2020),
the $F_{1}$ function is strictly declining from the maximum value at $\eta=0$
to the asymptotic value.

Lastly we address the violation of scaled null energy condition (NEC) in the
compactified coordinate $x$ demonstrated in Figs. 8(a)-(e) for
$\lambda=0.2,1,30,70,2$, respectively. The closed form of the expression for
the NEC can be found in our previous work Chew and Lim (2020). For all cases
of $\lambda$, the violation of scaled NEC at the throat of wormholes $x=0$
initially is maximum when $\beta=0$ where the corresponding wormhole is Ellis
wormhole (black curve). As $\beta$ increases, the violation of scaled NEC
decreases but then increases again and the scaled NEC becomes most violated
when $\beta$ approaches the critical value at the end. For $\lambda=0.2,1$,
the NEC at $x=0$ can be minimally violated because it reaches a very small
number which is close to zero. Besides, we observe that the appearance of
double throat configuration can affect the violation of NEC. For
$\lambda=0.2,1$, when $\beta$ approaches the critical value, both cases
develop a small peak indicating the NEC is satisfied at the location adjacent
to the asymptotically flat region (refer to the insets in Figs. 8(a)-(b)),
since the type I double throat appears in that limit of $\beta$. Nevertheless,
the peak disappears for large values of $\lambda$ because type I double throat
doesn’t exist anymore. Similarly, for the case of $\lambda=30$ and $70$, the
violation of scaled NEC at the throat $x=0$ decreases when $\beta$ increases
from zero, we see that the scaled NEC is satisfied at $x=0$ when the type II
double throat configuration appears with the throat at $x=0$ becomes the
equator. As $\beta$ continues to increase and type II double throat
configuration starts to diminish, the scaled NEC decreases and becomes
negative, hence it is violated again when the wormholes become single throat
again until $\beta$ approaches the critical value. Figs. 8(e)-(f) show the NEC
of $\lambda$=2 for type I + II double throat configuration which is related to
Fig. 6(a). Since the wormholes possess two types of double throat, they
certainly exhibit the phenomena of NEC due to these double throats. When the
wormholes possess a type II double throat in the range of $0.8806\leq\beta\leq
1.0074$ (left inset of Fig. 6(a)), the NEC at $x=0$ assumes a small positive
value which is slightly above the zero (inset of Fig. 8(e)). Similarly, a
small peak of positive NEC is noticed (inset of Fig. 8(f)) when the wormholes
exhibit a type I double throat near the asymptotically flat region (right
inset of Fig. 6(a)). However, in the BPS limit Chew and Lim (2020), the
violation of NEC decreases with the increase of $\beta$.

###  III.3 Junction Condition

We follow the approach of Kanti et. al. Kanti et al. (2012) to evaluate the
Einstein-matter field equations at $\eta=0$ for studying the junction
condition of wormholes,

$\left<G^{\mu}\,_{\nu}-\beta
T^{\mu}\,_{\nu}\right>=s^{\mu}\,_{\nu}\,,\quad\left<D_{\mu}F^{\mu\nu}-\frac{i}{4}\left[\Phi,D^{\nu}\Phi\right]\right>=\tilde{s}^{\nu}\,,\quad\left<D_{\mu}D^{\mu}\Phi-\lambda\left(\Phi^{2}-\upsilon^{2}\right)\Phi\right>=s_{2}\,.$
(36)

Here we have denoted by $s^{\mu}\,_{\nu}$ the stress-energy tensor for the
matter at the throat, by $\tilde{s}^{\nu}$ the source term for the SU(2)
vector fields and by $s_{2}$ the source term for the Higgs field. To obtain
the left-hand side of the junction conditions, one has to integrate the system
of EYMH equations across the boundary $l=0$, i.e., to evaluate the expressions

$\displaystyle\left<G^{\mu}\,_{\nu}-\beta T^{\mu}\,_{\nu}\right>$
$\displaystyle=\frac{1}{2}\lim_{L\rightarrow
0}\int_{-L}^{L}\left(G^{\mu}\,_{\nu}-\beta T^{\mu}\,_{\nu}\right)dl\,,$ (37)
$\displaystyle\left<D_{\mu}F^{\mu\nu}-\frac{i}{4}\left[\Phi,D^{\nu}\Phi\right]\right>$
$\displaystyle=\frac{1}{2}\lim_{L\rightarrow
0}\int_{-L}^{L}\left(D_{\mu}F^{\mu\nu}-\frac{i}{4}\left[\Phi,D^{\nu}\Phi\right]\right)dl\,,$
(38)
$\displaystyle\left<D_{\mu}D^{\mu}\Phi-\lambda\left(\Phi^{2}-\upsilon^{2}\right)\Phi\right>$
$\displaystyle=\frac{1}{2}\lim_{L\rightarrow
0}\int_{-L}^{L}\left(D_{\mu}D^{\mu}\Phi-\lambda\left(\Phi^{2}-\upsilon^{2}\right)\right)dl\,,$
(39)

where $dl=\sqrt{F_{1}}d\eta$.

Hence, the EYMH equations read

$\displaystyle G^{t}\,_{t}-\beta T^{t}\,_{t}$
$\displaystyle=\frac{F^{\prime\prime}_{1}}{F^{2}_{1}}+\frac{2\eta}{hF^{2}_{1}}F^{\prime}_{1}-\frac{3F^{\prime
2}_{1}}{4F^{3}_{1}}+\frac{\eta_{0}^{2}F_{1}}{h^{2}F_{1}}-\beta\frac{1}{2F_{1}}\psi^{\prime
2}+\beta\frac{(K^{2}-1)^{2}+2hK^{\prime 2}+h^{2}F_{1}H^{\prime
2}+2hF_{1}H^{2}K^{2}}{2F^{2}_{1}h^{2}}$
$\displaystyle\quad+\beta\frac{\lambda}{4}\left(H^{2}-\upsilon^{2}\right)^{2}\,,$
(40) $\displaystyle G^{\eta}\,_{\eta}-\beta T^{\eta}\,_{\eta}$
$\displaystyle=\left(\frac{F^{\prime}_{1}}{2F_{1}}+\frac{\eta}{h}\right)\frac{F^{\prime}_{0}}{F_{0}F_{1}}+\frac{F^{\prime
2}_{1}}{4F^{3}_{1}}+\frac{\eta}{hF^{2}_{1}}F^{\prime}_{1}-\frac{\eta_{0}^{2}}{h^{2}F_{1}}+\frac{\beta}{2F_{1}}\psi^{\prime
2}$ $\displaystyle\quad-\beta\frac{-(K^{2}-1)^{2}+2hK^{\prime
2}+h^{2}F_{1}H^{\prime
2}-2hF_{1}H^{2}K^{2}}{2h^{2}F^{2}_{1}}+\beta\frac{\lambda}{4}F_{1}(H^{2}-\upsilon^{2})^{2}\,,$
(41) $\displaystyle G^{\theta}\,_{\theta}-\beta T^{\theta}\,_{\theta}$
$\displaystyle=\frac{F^{\prime\prime}_{0}}{2F_{0}F_{1}}+\left(-\frac{F^{\prime}_{0}}{4F^{2}_{0}F_{1}}+\frac{\eta}{2hF_{0}F_{1}}\right)F^{\prime}_{0}+\frac{1}{2F^{2}_{1}}F^{\prime\prime}_{1}+\left(-\frac{F^{\prime}_{1}}{F_{1}}+\frac{\eta}{h}\right)\frac{F^{\prime}_{1}}{2F^{2}_{1}}+\frac{\eta_{0}^{2}}{h^{2}F_{1}}-\beta\frac{\psi^{\prime
2}}{2F_{1}}$ $\displaystyle\quad+\frac{-(K^{2}-1)^{2}+h^{2}F_{1}H^{\prime
2}}{2h^{2}F^{2}_{1}}+\beta\frac{\lambda}{4}(H^{2}-\upsilon^{2})^{2}\,,\,$ (42)
$\displaystyle G^{\phi}\,_{\phi}-\beta T^{\phi}\,_{\phi}$
$\displaystyle=\sin^{2}\theta\left(G^{\theta}\,_{\theta}-\beta
T^{\theta}\,_{\theta}\right)\,,$ (43) $\displaystyle
D_{\mu}F^{\mu\theta}-\frac{i}{4}\left[\Phi,D^{\theta}\Phi\right]$
$\displaystyle=\left[-\frac{K^{\prime\prime}}{2hF^{2}_{1}}-\frac{1}{4hF^{2}_{1}}\left(\frac{F^{\prime}_{0}}{F_{0}}-\frac{F^{\prime}_{1}}{F_{1}}\right)K^{\prime}+\frac{K(K^{2}-1+hF_{1}H^{2})}{2h^{2}F^{2}_{1}}\right]\tau_{\phi}\,,$
(44) $\displaystyle
D_{\mu}F^{\mu\phi}-\frac{i}{4}\left[\Phi,D^{\phi}\Phi\right]$
$\displaystyle=-\left[-\frac{K^{\prime\prime}}{2hF^{2}_{1}}-\frac{1}{4hF^{2}_{1}}\left(\frac{F^{\prime}_{0}}{F_{0}}-\frac{F^{\prime}_{1}}{F_{1}}\right)K^{\prime}+\frac{K(K^{2}-1+hF_{1}H^{2})}{2h^{2}F^{2}_{1}}\right]\frac{\tau_{\theta}}{\sin\theta}\,,$
(45) $\displaystyle
D_{\mu}D^{\mu}\Phi-\lambda\left(\Phi^{2}-\upsilon^{2}\right)\Phi$
$\displaystyle=\left[\frac{H^{\prime\prime}}{F_{1}}+\frac{1}{2F_{1}}\left(\frac{F^{\prime}_{0}}{F_{0}}+\frac{F^{\prime}_{1}}{F_{1}}+\frac{4\eta}{h}\right)H^{\prime}-\frac{2K^{2}}{hF_{1}}H-\lambda(H^{2}-\upsilon^{2})H\right]\tau_{\eta}\,.$
(46)

Next, we evaluate the left-hand sides of the above contributions

$\displaystyle\left<G^{t}\,_{t}-\beta T^{t}\,_{t}\right>$
$\displaystyle=\frac{F^{\prime}_{1}}{F_{1}}\,,$ (47)
$\displaystyle\left<G^{\eta}\,_{\eta}-\beta T^{\eta}\,_{\eta}\right>$
$\displaystyle=0\,,$ (48) $\displaystyle\left<G^{\theta}\,_{\theta}-\beta
T^{\theta}\,_{\theta}\right>$
$\displaystyle=\frac{F^{\prime}_{0}}{2F_{0}}+\frac{F^{\prime}_{1}}{2F_{1}}\,,$
(49) $\displaystyle\left<G^{\phi}\,_{\phi}-\beta T^{\phi}\,_{\phi}\right>$
$\displaystyle=\sin^{2}\theta\left<\left(G^{\theta}\,_{\theta}-\beta
T^{\theta}\,_{\theta}\right)\right>\,,$ (50)
$\displaystyle\left<D_{\mu}F^{\mu\theta}-\frac{i}{4}\left[\Phi,D^{\theta}\Phi\right]\right>$
$\displaystyle=-\frac{K^{\prime}}{2hF_{1}}\tau_{\phi}\,,$ (51)
$\displaystyle\left<D_{\mu}F^{\mu\phi}-\frac{i}{4}\left[\Phi,D^{\phi}\Phi\right]\right>$
$\displaystyle=\frac{K^{\prime}}{2hF_{1}}\frac{\tau_{\theta}}{\sin\theta}\,,$
(52)
$\displaystyle\left<D_{\mu}D^{\mu}\Phi-\lambda\left(\Phi^{2}-\upsilon^{2}\right)\Phi\right>$
$\displaystyle=H^{\prime}\tau_{\eta}\,,$ (53)

where the right-hand sides should be evaluated at the throat $\eta=0$.
Therefore we need to take into account the expansion at $\eta=0$ by assuming
the following forms,

$\displaystyle F_{0}(\eta)$
$\displaystyle=F_{00}+F_{02}\eta^{2}+F_{03}\eta^{3}+F_{04}\eta^{4}+......\,,$
(54) $\displaystyle F_{1}(\eta)$
$\displaystyle=F_{10}+F_{12}\eta^{2}+F_{13}\eta^{3}+F_{14}\eta^{4}+......\,,$
(55) $\displaystyle K(\eta)$
$\displaystyle=1+K_{1}\eta+K_{2}\eta^{2}+K_{3}\eta^{3}+K_{4}\eta^{4}+......\,,$
(56) $\displaystyle H(\eta)$
$\displaystyle=H_{1}\eta+H_{2}\eta^{2}+H_{3}\eta^{3}+H_{4}\eta^{4}+......\,.$
(57)

We substitute them into Eqs. (19)-(22), some of the leading terms are given by

$\displaystyle F_{02}$ $\displaystyle=\frac{\beta
F_{00}\left(4K^{2}_{1}-\lambda\eta^{2}_{0}F^{2}_{10}\upsilon^{4}\right)}{4\eta^{2}_{0}F_{10}}\,,$
(58) $\displaystyle F_{12}$ $\displaystyle=-\frac{1}{4}\beta\lambda
F^{2}_{10}\upsilon^{4}\,,$ (59) $\displaystyle K_{2}$ $\displaystyle=0\,,$
(60) $\displaystyle H_{2}$ $\displaystyle=0\,,$ (61) $\displaystyle F_{03}$
$\displaystyle=0\,,$ (62) $\displaystyle F_{13}$ $\displaystyle=0\,,$ (63)
$\displaystyle K_{3}$ $\displaystyle=-\frac{K_{1}\left(\beta
K^{2}_{1}-2F_{10}\right)}{6\eta^{2}_{0}F_{10}}\,,$ (64) $\displaystyle H_{3}$
$\displaystyle=-\frac{H_{1}\left(2\beta
K^{2}_{1}-\beta\lambda\eta^{2}_{0}\upsilon^{4}F^{2}_{10}+2\eta^{2}_{0}\lambda\upsilon^{2}F^{2}_{10}\right)}{12\eta^{2}_{0}F_{10}}\,,$
(65) $\displaystyle F_{04}$ $\displaystyle=\frac{\beta
K^{2}_{1}F_{00}\left(F_{10}+\beta
K^{2}_{1}\right)}{6\eta^{4}_{0}F^{2}_{10}}-\frac{\beta\lambda\upsilon^{2}F_{00}\left(\beta\upsilon^{2}K^{2}_{1}-F_{10}\upsilon^{2}-\eta^{2}_{0}H^{2}_{1}F_{10}\right)}{12\eta^{2}_{0}}+\frac{\beta^{2}\lambda^{2}\upsilon^{8}F_{00}F^{2}_{10}}{48}\,,$
(66) $\displaystyle F_{14}$
$\displaystyle=-\frac{\beta\left(3K^{2}_{1}+\eta^{2}_{0}H^{2}_{1}F_{10}\right)}{6\eta^{4}_{0}}+\frac{\beta\lambda\upsilon^{2}F_{10}\left(4\upsilon^{2}F_{10}+2\eta^{2}_{0}H^{2}_{1}F_{10}+\beta\upsilon^{2}K^{2}_{1}\right)}{24\eta^{2}_{0}}+\frac{\beta^{2}\lambda^{2}\upsilon^{8}F^{3}_{10}}{48}\,,$
(67) $\displaystyle K_{4}$
$\displaystyle=\frac{\eta^{2}_{0}H^{2}_{1}F_{10}+3K^{2}_{1}}{12\eta^{2}_{0}}\,,$
(68) $\displaystyle H_{4}$ $\displaystyle=\frac{H_{1}K_{1}}{3\eta^{2}_{0}}\,.$
(69)

The expansion shows that $F^{\prime}_{0}(0)=F^{\prime}_{1}(0)=0$, whereas
$K^{\prime}_{0}(0)$ and $H^{\prime}(0)$ are finite. Therefore we only need to
include the source terms for the vector field and for the Higgs field. Thus we
assume

$\tilde{s}^{\nu}=g^{\nu\mu}\tilde{s}_{\mu}\,,\quad\tilde{s}_{\mu}dx^{\mu}=\frac{\upsilon}{2}\left(\tau_{\phi}d\theta-\tau_{\theta}\sin\theta
d\phi\right)\,,\quad s_{2}=s\tau_{\eta}\,.$ (70)

The two vector equations then yield

$-K^{\prime}=\upsilon\,,$ (71)

while the Higgs equation yields

$H^{\prime}=s\,,$ (72)

with the left-hand sides evaluated at the throat.

##  IV Conclusion and Outlook

We have constructed the symmetric wormholes in which the throat is supported
by the phantom field in the Einstein-Yang-Mills-Higgs (EYMH) system beyond the
Bogomol’nyi-Prasad-Sommerfield (BPS) limit by including the Higgs potential
with the Higgs self-interaction value $\lambda$. The wormhole spacetime is
symmetric with respect to the radial coordinate $\eta=0$. Analogous to the BPS
limit Chew and Lim (2020), the wormholes possess the probe limit which is the
Yang-Mills-Higgs (YMH) field in the background of the Ellis wormhole when the
gravity is absent. In the presence of gravity, the wormholes possess the non-
trivial non-abelian hair where the families of hairy wormholes solutions
emerge from the Ellis wormhole. In the BPS limit, the mass of wormholes and
the scalar charge of the phantom field increase monotonically when the
gravitational strength increases. However, beyond the BPS limit by increasing
$\lambda$, the masses of wormholes initially increase from null to a maximum
positive value but then decrease dramatically to a negative value when the
gravitational strength approaches a critical value for small values of
$\lambda$ but strictly decreases to negative for high values of $\lambda$. The
scaled scalar charge for the phantom field also increases very steeply at the
critical value of gravitational strength compare to the BPS case. Note that
when $\lambda$ increases, the domain of existence of wormholes decreases, thus
the drastic changing of properties of wormholes occurs earlier since the
gravitational strength approaches the critical value earlier. This is in
contrast to the BPS limit where it can take very large values of gravitational
strength Chew and Lim (2020).

The gauge field and Higgs field behave monotonically where the gauge field
assumes a maximum value at $\eta=0$ and decreases to zero at the asymptotic
region while the Higgs field increases from zero at $\eta=0$ to its asymptotic
value. Beyond the BPS limit when the gravitational strength approaches a
critical value, the Higgs field increases very steeply toward its asymptotic
value while the gauge field approaches unity. Similarly, the behaviour of
metric functions also changes very sharply toward their asymptotic values at
infinity. Furthermore, the precision of wormhole solutions drop significantly
when the gravitational strength approaches the critical value, thus we are not
able to study the properties of wormholes beyond that critical value.

The hairy wormholes possess a double throat configuration in the BPS limit
when the gravitational strength exceeds a value. By varying $\lambda$, the
hairy wormholes exhibit richer geometrical structures than their counterpart
in the BPS limit for the double throat configurations, namely, type I, type II
and the hybrid of type I and II. Type I and II double throat configurations
exist for $0<\lambda<3.0283$ and $1.8814\leq\lambda\leq 100$, respectively.
The hybrid of type I and II appears for $1.8814\leq\lambda<3.0283$ where type
I and II configurations can coexist together.

In type I the wormholes only exhibit double throat configuration when the
gravitational strength approaches a critical value where the throat at
$\eta=0$ still remains but simultaneous another throat is formed near the
asymptotically flat region and the equator is sandwiched between them.
However, the wormholes merely possess a single throat when the gravitational
strength before approaching the critical value.

On the contrary, in type II double throat configuration the wormholes start to
develop the double throat configuration for small values of gravitational
strength which is far below the critical value. The throat at $\eta=0$ becomes
the equator and another throat concurrently develops in the vicinity of
$\eta=0$. As gravitational strength further increases, the distance between
that throat and the equator increases to a maximum and then decreases until
that throat disappears. Hence, the wormholes transform back from double throat
to single throat until the gravitational strength approaches the critical
value. Besides, we find that type II double throat appears earlier for small
gravitational strength when $\lambda$ increases. Moreover, the hybrid of type
I and II double throat configuration acts as a transition state from type I to
II, where type I gradually disappears with the increase of $\lambda$.
Therefore, the non-zero of Higgs self-interaction value can give rise to the
interesting geometrical structure of double throat for the wormholes.

The null energy condition (NEC) is most violated at the throat $\eta=0$ when
the gravitational strength vanishes since the Ellis wormhole is the trivial
solution. For small values of $\lambda$, the increase of gravitational
strength initially weakens the violation of NEC at $\eta=0$ to a minimum value
which is close to zero but then the corresponding violation increases again
and the NEC becomes most violated when the gravitational strength approaches
the critical value. We find that the appearance of the double throat can
affect the NEC. When the type I double throat exists, a small peak of NEC
develops near the asymptotically flat region when the gravitational strength
approaches the critical value, hence the NEC is satisfied there. However, the
peak disappears for large value of $\lambda$ because type I double throat
doesn’t exist anymore. Whereas for those wormholes with a higher value of
$\lambda$, as the gravitational strength increases, the violation of NEC at
$\eta=0$ is also weakening and then NEC can become positive, thus NEC can be
satisfied also with the appearance of type II double throat configuration.
However, when type II double throat gradually disappears, the NEC at $\eta=0$
decreases and then becomes negative, hence it is most violated when the
gravitational strength approaches the critical value. When the wormholes
possess the hybrid of type I and II double throats, they indeed exhibit the
phenomena of violation of NEC for small and large values of $\lambda$. In the
BPS limit, the violation of NEC decreases when the gravitational strength
increases.

On the other hand, since the first-order derivatives of YMH fields are finite
at $\eta=0$, thus there will be a discontinuity of YMH equations which are
evaluated at $\eta=0$, therefore we need to introduce the extra source terms
to overcome the discontinuity.

Let us comment on the stability issue of the non-Abelian wormholes. As was
pointed out in our previous paper Chew and Lim (2020), the wormholes supported
by a phantom field are generally unstable against linear perturbations
Gonzalez et al. (2009a, b); Torii and Shinkai (2013); Dzhunushaliev et al.
(2013b, 2014b); Aringazin et al. (2015). Likewise, the particle-like and hairy
black hole solutions of the EYMH system are unstable as well Greene et al.
(1993); Winstanley and Mavromatos (1995); Mavromatos and Winstanley (1996).
Therefore, we conjecture that the EYMH hairy wormholes will inherit these
instabilities and should behave qualitatively unstable. Hence, in the future,
we plan to carry out an exhaustive numerical analysis of full linear stability
for hairy wormholes consistently by perturbing all the functions. However, the
calculation of (un)stable modes would be tedious and very challenging since
the presence of the YMH field makes the calculation becomes non-trivial.
Besides, the unstable modes disappear for sufficiently rapidly rotating
5-dimensional Ellis wormholes with equal angular momenta Dzhunushaliev et al.
(2013a). Since the counterpart EYMH black holes can rotate, then it is
interesting to construct the rotating EYMH wormholes which might be stable
against the perturbation.

Moreover, the static and regular EYMH solutions can also possess only axial
symmetry and need not be spherically symmetric, their counterpart static black
holes also can possess only axially symmetric horizon Hartmann et al. (2001a,
b) which are the counterexamples to Israel’s theorem. Therefore, as a first
step to constructing the rotating wormholes in EYMH, we could consider
constructing the static hairy wormhole solutions with a throat that is also
axially symmetric.

## Acknowledgement

XYC is supported by the National Research Foundation of Korea (Grant No.:
2021R1C1C1008622, 2021R1A4A5031460). We really appreciate having a useful
discussion with Jutta Kunz and Burkhard Kleihaus about the junction condition
of wormhole.

## References

  * Rossi (1982) P. Rossi, Physics Reports 86, 317 (1982).
  * Shnir (2005) Y. M. Shnir, _Magnetic Monopoles_ , Text and Monographs in Physics (Springer-Verlag Berlin Heidelberg, Berlin/Heidelberg, 2005), ISBN 978-3-540-25277-1, 978-3-540-29082-7.
  * t Hooft (1974) G. t Hooft, Nucl. Phys. B 79, 276 (1974).
  * Polyakov (1974) A. M. Polyakov, JETP Letters 20, 194 (1974).
  * Prasad and Sommerfield (1975) M. Prasad and C. M. Sommerfield, Phys. Rev. Lett. 35, 760 (1975).
  * Rebbi and Rossi (1980) C. Rebbi and P. Rossi, Phys. Rev. D 22, 2010 (1980).
  * Forgács et al. (1981) P. Forgács, Z. Horváth, and L. Palla, Phys. Lett. B 99, 232 (1981).
  * Bogomol’Nyi (1976) E. Bogomol’Nyi, Sov. J. Nucl. Phys. 24 (1976).
  * Bogomol’nyi and Marinov (1976) E. Bogomol’nyi and M. S. Marinov, Sov. J. Nucl. Phys. 23, 357 (1976).
  * Bais and Primack (1976) F. Bais and J. Primack, Physical Review D 13, 819 (1976).
  * Teh et al. (2010) R. Teh, K.-M. Wong, and K.-G. Lim, International Journal of Modern Physics A 25, 5731 (2010).
  * Kleihaus et al. (1998) B. Kleihaus, J. Kunz, and D. Tchrakian, Modern Physics Letters A 13, 2523 (1998).
  * Kleihaus and Kunz (1999) B. Kleihaus and J. Kunz, Phys. Rev. D 61, 025003 (1999).
  * Kleihaus et al. (2003) B. Kleihaus, J. Kunz, and Y. Shnir, Phys. Lett. B 570, 237 (2003).
  * Kleihaus et al. (2004a) B. Kleihaus, J. Kunz, and Y. Shnir, Physical Review D 70, 065010 (2004a).
  * Lee et al. (1992) K. Lee, V. Nair, and E. J. Weinberg, Phys. Rev. D 45, 2751 (1992).
  * Breitenlohner et al. (1992) P. Breitenlohner, P. Forgacs, and D. Maison, Nucl. Phys. B 383, 357 (1992).
  * Breitenlohner et al. (1995) P. Breitenlohner, P. Forgacs, and D. Maison, Nucl. Phys. B 442, 126 (1995), eprint gr-qc/9412039.
  * Brihaye et al. (1998) Y. Brihaye, B. Hartmann, and J. Kunz, Phys. Lett. B 441, 77 (1998).
  * Lue and Weinberg (1999) A. Lue and E. J. Weinberg, Phys. Rev. D 60, 084025 (1999), eprint hep-th/9905223.
  * Brihaye et al. (1999) Y. Brihaye, B. Hartmann, J. Kunz, and N. Tell, Phys. Rev. D 60, 104016 (1999), eprint hep-th/9904065.
  * Brihaye et al. (2000a) Y. Brihaye, B. Hartmann, and J. Kunz, Phys. Rev. D 62, 044008 (2000a), eprint hep-th/9911148.
  * Brihaye et al. (2000b) Y. Brihaye, F. Grard, and S. Hoorelbeke, Phys. Rev. D 62, 044013 (2000b), eprint hep-th/9912023.
  * Brihaye and Hartmann (2002) Y. Brihaye and B. Hartmann, Phys. Rev. D 66, 064018 (2002), eprint hep-th/0206004.
  * Hartmann et al. (2001a) B. Hartmann, B. Kleihaus, and J. Kunz, Phys. Rev. Lett. 86, 1422 (2001a), eprint hep-th/0009195.
  * Hartmann et al. (2001b) B. Hartmann, B. Kleihaus, and J. Kunz, Phys. Rev. D 65, 024027 (2001b).
  * Kleihaus and Kunz (2000a) B. Kleihaus and J. Kunz, Phys. Rev. Lett. 85, 2430 (2000a), eprint hep-th/0006148.
  * Kleihaus and Kunz (2000b) B. Kleihaus and J. Kunz, Phys. Lett. B 494, 130 (2000b), eprint hep-th/0008034.
  * Van der Bij and Radu (2002) J. Van der Bij and E. Radu, Int. J. Mod. Phys. A 17, 1477 (2002), eprint gr-qc/0111046.
  * Paturyan and Tchrakian (2004) V. Paturyan and D. Tchrakian, J. Math. Phys. 45, 302 (2004), eprint hep-th/0306160.
  * Paturyan et al. (2005) V. Paturyan, E. Radu, and D. Tchrakian, Phys. Lett. B 609, 360 (2005), eprint hep-th/0412011.
  * Kleihaus et al. (2004b) B. Kleihaus, J. Kunz, and F. Navarro-Lerida, Phys. Lett. B 599, 294 (2004b), eprint gr-qc/0406094.
  * Kleihaus et al. (2005a) B. Kleihaus, J. Kunz, and Y. Shnir, Phys. Rev. D 71, 024013 (2005a), eprint gr-qc/0411106.
  * Kleihaus et al. (2005b) B. Kleihaus, J. Kunz, and U. Neemann, Phys. Lett. B 623, 171 (2005b), eprint gr-qc/0507047.
  * Hajicek (1983a) P. Hajicek, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 386, 223 (1983a), ISSN 00804630\.
  * Hajicek (1983b) P. Hajicek, Journal of Physics A: Mathematical and General 16, 1191 (1983b).
  * Degen (1987) F. Degen, General Relativity and Gravitation 19, 739 (1987).
  * Chew and Lim (2020) X. Y. Chew and K.-G. Lim, Phys. Rev. D 102, 124068 (2020), eprint 2009.13334.
  * Ellis (1973) H. G. Ellis, J. Math. Phys. 14, 104 (1973).
  * Ellis (1979) H. G. Ellis, Gen. Relativ. Gravit. 10, 105 (1979).
  * Bronnikov (1973) K. A. Bronnikov, Acta. Phys. Pol p. B4 (1973).
  * Gonzalez et al. (2009a) J. Gonzalez, F. Guzman, and O. Sarbach, Class. Quant. Grav. 26, 015010 (2009a), eprint 0806.0608.
  * Gonzalez et al. (2009b) J. Gonzalez, F. Guzman, and O. Sarbach, Class. Quant. Grav. 26, 015011 (2009b), eprint 0806.1370.
  * Torii and Shinkai (2013) T. Torii and H.-a. Shinkai, Phys. Rev. D 88, 064027 (2013).
  * Kashargin and Sushkov (2008a) P. Kashargin and S. Sushkov, Gravit. Cosmol. 14, 80 (2008a).
  * Kashargin and Sushkov (2008b) P. Kashargin and S. Sushkov, Phys. Rev. D 78, 064071 (2008b).
  * Kleihaus and Kunz (2014) B. Kleihaus and J. Kunz, Phys. Rev. D 90, 121503 (2014).
  * Chew et al. (2016) X. Y. Chew, B. Kleihaus, and J. Kunz, Phys. Rev. D 94, 104031 (2016).
  * Dzhunushaliev et al. (2013a) V. Dzhunushaliev, V. Folomeev, B. Kleihaus, J. Kunz, and E. Radu, Phys. Rev. D 88, 124028 (2013a).
  * Chew et al. (2018) X. Y. Chew, B. Kleihaus, and J. Kunz, Phys. Rev. D 97, 064026 (2018).
  * Karakasis et al. (2021) T. Karakasis, E. Papantonopoulos, and C. Vlachos (2021), eprint 2107.09713.
  * Huang et al. (2020) H. Huang, H. Lü, and J. Yang (2020), eprint 2010.00197.
  * Blázquez-Salcedo et al. (2020) J. L. Blázquez-Salcedo, X. Y. Chew, J. Kunz, and D.-H. Yeom (2020), eprint 2012.06213.
  * Hauser et al. (2014) O. Hauser, R. Ibadov, B. Kleihaus, and J. Kunz, Phys. Rev. D 89, 064010 (2014).
  * Bartnik and McKinnon (1988) R. Bartnik and J. McKinnon, Phys. Rev. Lett. 61, 141 (1988).
  * Nedkova et al. (2013) P. G. Nedkova, V. K. Tinchev, and S. S. Yazadjiev, Phys. Rev. D 88, 124019 (2013).
  * Gyulchev et al. (2018) G. Gyulchev, P. Nedkova, V. Tinchev, and S. Yazadjiev, Eur. Phys. J. C 78, 544 (2018).
  * Amir et al. (2019) M. Amir, A. Banerjee, and S. D. Maharaj, Ann. Phys. 400, 198 (2019).
  * Narzilloev et al. (2021) B. Narzilloev, D. Malafarina, A. Abdujabbarov, B. Ahmedov, and C. Bambi (2021), eprint 2105.09174.
  * Benavides-Gallego et al. (2021) C. A. Benavides-Gallego, W.-B. Han, D. Malafarina, B. Ahmedov, and A. Abdujabbarov (2021), eprint 2107.07998.
  * Jusufi et al. (2021) K. Jusufi, S. K., M. Azreg-Aïnou, M. Jamil, Q. Wu, and C. Bambi (2021), eprint 2106.08070.
  * Bambi and Stojkovic (2021) C. Bambi and D. Stojkovic, Universe 7, 136 (2021), eprint 2105.00881.
  * Abe (2010) F. Abe, Astrophys. J. 725, 787 (2010).
  * Toki et al. (2011) Y. Toki, T. Kitamura, H. Asada, and F. Abe, Astrophys. J. 740, 121 (2011).
  * Takahashi and Asada (2013) R. Takahashi and H. Asada, Astrophys J. Lett. 768, L16 (2013).
  * Cramer et al. (1995) J. G. Cramer, R. L. Forward, M. S. Morris, M. Visser, G. Benford, and G. A. Landis, Phys. Rev. D 51, 3117 (1995).
  * Perlick (2004) V. Perlick, Phys. Rev. D 69, 064017 (2004).
  * Tsukamoto et al. (2012) N. Tsukamoto, T. Harada, and K. Yajima, Phys. Rev. D 86, 104062 (2012).
  * Bambi (2013) C. Bambi, Phys. Rev. D 87, 107501 (2013).
  * Zhou et al. (2016) M. Zhou, A. Cardenas-Avendano, C. Bambi, B. Kleihaus, and J. Kunz, Phys. Rev. D 94, 024036 (2016).
  * Deligianni et al. (2021a) E. Deligianni, J. Kunz, P. Nedkova, S. Yazadjiev, and R. Zheleva (2021a), eprint 2103.13504.
  * Deligianni et al. (2021b) E. Deligianni, B. Kleihaus, J. Kunz, P. Nedkova, and S. Yazadjiev (2021b), eprint 2107.01421.
  * Blázquez-Salcedo et al. (2018) J. L. Blázquez-Salcedo, X. Y. Chew, and J. Kunz, Phys. Rev. D 98, 044035 (2018).
  * Dzhunushaliev et al. (2014a) V. Dzhunushaliev, V. Folomeev, C. Hoffmann, B. Kleihaus, and J. Kunz, Phys. Rev. D 90, 124038 (2014a).
  * Hoffmann et al. (2017) C. Hoffmann, T. Ioannidou, S. Kahlen, B. Kleihaus, and J. Kunz, Phys. Rev. D 95, 084010 (2017).
  * Ascher et al. (1979) U. Ascher, J. Christiansen, and R. D. Russell, Math. Comput. 33, 659 (1979).
  * Kierzenka and Shampine (2001) J. Kierzenka and L. F. Shampine, ACM Transactions on Mathematical Software (TOMS) 27, 299 (2001).
  * Corichi et al. (2000) A. Corichi, U. Nucamendi, and D. Sudarsky, Phys. Rev. D 62, 044046 (2000), eprint gr-qc/0002078.
  * Ashtekar et al. (2001) A. Ashtekar, A. Corichi, and D. Sudarsky, Class. Quant. Grav. 18, 919 (2001), eprint gr-qc/0011081.
  * Kanti et al. (2012) P. Kanti, B. Kleihaus, and J. Kunz, Phys. Rev. D 85, 044007 (2012).
  * Dzhunushaliev et al. (2013b) V. Dzhunushaliev, V. Folomeev, B. Kleihaus, and J. Kunz, Phys. Rev. D 87, 104036 (2013b).
  * Dzhunushaliev et al. (2014b) V. Dzhunushaliev, V. Folomeev, B. Kleihaus, and J. Kunz, Phys. Rev. D 89, 084018 (2014b).
  * Aringazin et al. (2015) A. Aringazin, V. Dzhunushaliev, V. Folomeev, B. Kleihaus, and J. Kunz, J. Cosmol. Astropart. Phys. 2015, 005 (2015).
  * Greene et al. (1993) B. R. Greene, S. D. Mathur, and C. M. O’neill, Phys. Rev. D 47, 2242 (1993).
  * Winstanley and Mavromatos (1995) E. Winstanley and N. E. Mavromatos, Phys. Lett. B 352, 242 (1995).
  * Mavromatos and Winstanley (1996) N. E. Mavromatos and E. Winstanley, Phys. Rev. D 53, 3190 (1996).