# BPS Skyrmions of Generalized Skyrme Model in Higher Dimensions

Emir Syahreza Fadhilla1,2<EMAIL_ADDRESS>Bobby Eka Gunara1
<EMAIL_ADDRESS>(Corresponding author) Ardian Nata Atmaja2
<EMAIL_ADDRESS>1Theoretical High Energy Physics Research Division,
Institut Teknologi Bandung, Jl. Ganesha 10 Bandung 40132, Indonesia.  2
Research Center for Quantum Physics, National Research and Innovation Agency
(BRIN), Kompleks PUSPIPTEK Serpong, Tangerang 15310, Indonesia.

###### Abstract

In this work we consider the higher dimensional Skyrme model, with spatial
dimension $d>3$, focusing on its BPS submodels and their corresponding
features. To accommodate the cases with a higher topological degree, $B\geq
1$, a modified generalized hedgehog ansatz is used where we assign an integer
$n_{i}$ for each rotational plane, resulting in a topological degree that
proportional to product of these integers. It is found via BPS Lagrangian
method that there are only two possible BPS submodels for this spherically
symmetric ansatz which shall be called as BPS Skyrme model and scale-invariant
model. The properties of the higher dimensional version of both submodels are
studied and it is found that the BPS Skyrmions with $B\geq 1$ exist in the
first submodel but there is only $B=1$ BPS Skyrmion in the second submodel. We
also study the higher dimensional version of self-duality conditions in terms
of strain tensor eigenvalues and find that, in general, the scale-invariant
model has a stronger self-duality condition than the BPS Skyrme model.

###### Contents

  1. I Introduction
  2. II The Skyrme Model
    1. II.1 Spacetime Background and the suitable ansatz
    2. II.2 The BPS Submodels
  3. III BPS Bounds and BPS Type Solutions Via Bogomolnyi Equations
    1. III.1 BPS Bounds of The First and Second Submodel
    2. III.2 BPS Solutions
      1. III.2.1 The First Submodel
      2. III.2.2 The Second Submodel
  4. IV Conclusions and Outlooks
  5. A On The Polyspherical Coordinate
    1. A.1 Even $d$ Case
    2. A.2 Odd $d$ Case
  6. B Bogomolnyi equation from Equation of Motions
    1. B.1 Field Equation of The First BPS Submodel
    2. B.2 Field Equation of The Second BPS Submodel
  7. C Components of Energy-Momentum Tensor via Eigenvalues of Strain Tensor

##  I Introduction

The Skyrme model was initially designed to be a model for strongly interacting
particles. The construction of this model is based on the non-linear sigma
model which is slightly modified by the addition of a higher-order term known
as the Skyrme term [1, 2]. This modification is necessary because in the flat
four-dimensional spacetime the two terms scale in opposite ways under
coordinate scaling $\vec{x}\rightarrow\mu\vec{x}$ which means that the
solution of the Skyrme model field equations is stable under the conditions of
Derrick-type stability [3, 4]. The Skyrme model, although eventually failed to
become the general theory for strong interactions, has many uses, one of which
is to be an effectively good model for strong interactions in low energy
regions. This low energy scheme is quite interesting to be studied extensively
because there is no scheme through QCD that can predict hadronic phenomena in
low energy regions [5]. The growing interest in the Skyrme model emerged after
the discovery of the correspondence on large $N_{c}$ limit between such
effective field theory in low energy regime and the general formulation of the
QCD [6, 7].

Although Skyrme formulated this model as a model for strong interactions,
currently Skyrmion is a phenomenon that can be encountered in many other non-
linear systems, especially at low spatial dimensions. Examples of skyrmion
phenomena that can be found in solids varied from the superfluid polariton
[8], Bose-Einstein ferromagnets [9], antiferromagnets [10], magnetic thin
films [11], to chiral nematic liquid crystal [12]. In addition, the phenomenon
of a skyrmion is also quite common in optics such as evanescent waves [13].
This raises a new challenge in the generalization of the Skyrme model in any
spacetime dimension. In fact, we have an exotic matter in higher spacetime
dimensions that is constructed in the same manner, namely skyrmion string that
belongs to the $D$-brane solitons [14, 15]. The Skyrmion string is a Skyrmion-
like string which ends on domain walls since the domain walls are located on
the potential term’s vacua and supported by
$\pi_{d-1}(S^{d-1})=\mathbb{Z}$111It is conventional to have potential terms
with the vacuum expectation value of the skyrmion field minimizes the
potential, i.e. become the vacua of the potential terms.. The original
construction of this skyrmion can only be carried out up to $1+6$ dimensional
spacetime since the Skyrme action contains only up to sextic order term. A
similar object can also be found in the holography theory, known as Skyrme
black brane where a Skyrme model coupled with gravity in asymptotically five-
dimensional AdS spacetime [16]. Another brane-like object in $1+6$ dimensions
that has a skyrmionic solution is the Skyrme brane which has a slightly
different origin from the previously mentioned skyrmion string [17]. A skyrme
brane constructed from a skyrme action up to quartic derivative term (also
known as the Skyrme term) governing the dynamics of the Skyrme field as a map
from three-dimensional submanifold of $1+6$ dimensional spacetime (the
$3$-brane) to a $S^{3}$ manifold represented as $SU(2)$-valued field. It is a
gravitating model with seven-dimensional Einstein-Skyrme action possessing a
spherically symmetric 3-brane along the transverse directions in the bulk. In
contrast with the previous D-brane soliton that lives on a completely flat
higher-dimensional spacetime.

To construct such a higher dimensional object, there is a problem with the
Derrick’s stability conditions because the higher the spacetime dimensions,
the more higher-order terms must be added. The consistent formulation for this
high-order Skyrme model has been carried out in [18] and it is found in this
way that terms consistent with the earlier Skyrme formulation can only be
carried out up to sextic order for four-dimensional space-time. This
generalized model in four-dimensional spacetime with sextic term have been
studied in [19, 20, 21, 22], and the same sextic model have also been studied
on five-dimensional spacetime in [23]. One can also construct a higher-order
Skyrme model through Sutcliffe truncation on quartic Yang-Mills action in an
eight-dimensional spacetime by removing all vector meson terms from the gauge
field expansion, leading to a Skyrme model that obeys Derrick’s stability
condition [24]. It is known that this procedure can produce a higher-order
Skyrme model in spacetime dimensions of multiple of four by applying Sutcliffe
truncation on a $4k$-th order Yang-Mills action. Another approach for the
extension of the Skyrme model in higher dimensions can be done by adding
higher-order terms similar to the one demonstrated in [18] up to $2d$-th order
term which is proportional to the squared norm of topological current and this
model is originally proposed as an extension of $O(d+1)$ sigma model [25]. In
this work, we start our discussion by referring to this $O(d+1)$ sigma model
construction and motivate it with the generalization of Manton’s construction
of static four-dimensional Skyrme model through invariants of strain tensor
for higher-dimensional cases in section II. Here we use the strain tensor
eigenvalue combinations as the representation of higher-order Skyrme terms
which have been proposed in [26] to describe the geometrical properties of
four-dimensional spacetime Skyrmion as a harmonic map of a flat spatial
manifold to a three-dimensional sphere manifold.

As the number of terms grows for a higher dimensional model, it is of interest
to find the term which has special qualitative features. Some examples of this
approach can be found in [23] which shows that omitting quartic term does not
give any significant qualitative feature in the five-dimensional spacetime
model and in [27, 28] which shows that the quartic term is necessary for black
hole solution to exist. Here, we discuss which term gives a BPS submodel for
an arbitrary dimensional Skyrme model. The BPS Skyrme submodel itself is a
model within the generalized Skyrme model which has the maximum symmetry of
the map from spatial space to target space leading to a linear lower bound of
energy in the topological degree of the map. It is known that such a lower
bound exists in four and three-dimensional spacetime Skyrme model [29, 30,
31]. A five-dimensional BPS Skyrme submodel for a specific spherically
symmetric system also has been found in [32]. Our main focus in this work is
to find the conditions needed for a generalized Skyrme model in arbitrary
spatial dimension $d$ to satisfy the BPS limit. To do so, we use the extension
of the BPS Lagrangian method for the Skyrme model which has been demonstrated
for four-dimensional spacetime in [33]. With the conditions for the BPS
Skyrmion at hand, we can deduce the BPS bound on the Skyrmion energy and the
corresponding first-order differential field equation (Bogomolnyi equation)
which could be exploited to derive some features of the BPS Skyrmion. There
are two known higher dimensional BPS Skyrme Models, firstly proposed in [25]
and, in this work, we prove that in spherically symmetric case there are no
other submodels that could saturate the BPS limit, this result is demonstrated
in section II.2. In the lower spacetime dimensional case the saturation of the
BPS limit implies that the total energy is linear to the topological degree
(called baryon number), mostly related to the phenomenology of baryonic
particles. The linear energy as a function of the baryon number is a well-
established fact in nuclear physics [20], but for higher-dimensional objects
such as $D$-brane solitons and strings, such BPS property have a richer
implication even though we have not really understood the physical mechanism
behind it. We also explore some of these physical properties related to the
BPS sector of the Skyrme model as a model for our higher dimensional object.

In this paper we use early alphabets, $a,b,c,\dots$, as scalar multiplets
indices and mid alphabets, $i,j,k,\dots$ as spacetime indices. Both indices
sets take value from $0$ to $d$ where $d+1$ is the spacetime dimensionality.
The Einstein summation convention is always to be understood for any repeating
indices for both tensors and scalar multiplets. i represents the imaginary
number basis, i.e. $\textbf{i}=\sqrt{-1}$ and the combinatorial notation
$C^{p}_{q}$ reads ”$p$ choose $q$” and explicitly given by

$C^{p}_{q}\equiv\frac{p!}{(p-q)!q!}~{},$ (1)

for $q\leq p$ and zero for $q>p$. Because the construction of the theory
relies heavily on combinations of some objects such as eigenvalues and winding
numbers, it is useful for us to simplify the notations by exploiting
elementary symmetric polynomials. The $k$-th elementary symmetric polynomials
from $n$ variables are defined as

$e_{k}(x_{1},x_{2},\dots,x_{n})\equiv\sum_{i_{1}\neq i_{2}\neq\dots\neq
i_{k}}^{n}\frac{1}{k!}x_{i_{1}}x_{i_{2}}\dots x_{i_{k}}~{},$ (2)

and we can calculate the expression of $e_{k}(x_{1},x_{2},\dots,x_{n})$ for
each $k$ recursively via Girard-Newton formula, namely

$ke_{k}=\sum_{i=1}^{k}(-1)^{i-1}e_{k-1}p_{i}$ (3)

for all $n\geq 1$ and $n\geq k\geq 1$ where we have defined $p_{k}\equiv
p_{k}(x_{1},x_{2},\dots,x_{n})$ as

$p_{k}\equiv\sum_{i=1}^{n}x^{k}_{i}~{}.$ (4)

Throughout this paper, we are going to consider the Minkowskian spacetime
manifold which has a compact submanifold, given by

$ds^{2}=-dt^{2}+dr^{2}+r^{2}d\Omega_{d-1}^{2}~{}.$ (5)

The compact submanifold with metric $d\Omega_{d-1}^{2}$ resembles a $S^{d-1}$
which allows $N=\left\lfloor\frac{d}{2}\right\rfloor$ rotational planes.
Explicit choices of the coordinate systems on the manifold of
$d\Omega_{d-1}^{2}$ are described in a more detailed manner in later sections.

##  II The Skyrme Model

Let us first consider the standard Skyrme model action functional in $3+1$
spacetime dimensions given by

$\mathcal{S}_{4}=\int
d^{4}x\sqrt{-g}\left(\gamma_{1}g^{ij}~{}\mathcal{B}_{i}^{a}\mathcal{B}_{j}^{a}+\gamma_{2}g^{ij}g^{kl}~{}\mathcal{C}^{a}_{ik}\mathcal{C}^{a}_{jl}\right)~{},$
(6)

with
$\mathcal{B}^{a}_{i}=\frac{1}{2\textbf{i}}\text{Tr}\left(U^{\dagger}\sigma^{a}\partial_{i}U\right)$
and $\mathcal{C}^{a}_{ij}=\varepsilon^{abc}B^{b}_{i}B^{c}_{j}$ where
$U=\phi^{0}I+\textbf{i}(\phi^{1}\sigma^{1}+\phi^{2}\sigma^{2}+\phi^{3}\sigma^{3})$
is an $SU(2)$ valued chiral field originally proposed by T. Skyrme [1, 2] and
$\varepsilon^{abc}$ is Levi-Civita symbol. The vector
$\sigma^{a}=(\sigma^{1},\sigma^{2},\sigma^{3})$ are Pauli matrices which are
the generators of $SU(2)$, while $(\phi^{a})$ with $a=0,1,2,3$, are real
scalar fields satisfying $O(4)$ model condition, $\phi^{a}\phi^{a}=1$, see
also for example, [23, 3, 25].

Now, we introduce the strain tensor which is defined as $D\equiv JJ^{T}$ where
$J$ is the Jacobian of generalized harmonic map $\phi=(\phi^{0},\vec{\phi})$.
In the recipe given in [26], we can rewrite the two terms in the Lagrangian of
the static case as the invariants of the harmonic map, namely
$\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}$ and
$\lambda_{1}^{2}\lambda_{2}^{2}+\lambda_{1}^{2}\lambda_{3}^{2}+\lambda_{2}^{2}\lambda_{3}^{2}$
with $\lambda_{1}^{2},\lambda_{2}^{2},\lambda_{3}^{2}$ are the eigenvalues of
$D$. We can observe that the use of the strain tensor eigenvalues reduces the
problem of finding the BPS bound into an algebraic problem. Thus, we can
exploit this method to generalize the static Skyrme model

Let us discuss the details of this strain tensor for higher dimensional static
Skyrme model as follows. For a map from a $d$ dimensional Euclidean space,
$\mathbb{R}^{d}$, to a $d$ dimensional sphere, $S^{d}$, we have a $d\times d$
Jacobian matrix $J$, leading to a $d\times d$ strain tensor. Suppose we have
the eigenvalue of $D$ from a $d$ dimensional Skyrme Model as
$\lambda_{1}^{2},\lambda_{2}^{2},\dots,\lambda_{d}^{2}$. The corresponding
Lagrangian can be written in the most general way as

$\mathcal{L}=-\gamma_{0}V-\sum_{n=1}^{d}\gamma_{n}\mathcal{L}_{2n}~{},$ (7)

where $\gamma_{i}\in\mathbb{R}$ are coupling constants, $V\equiv V(\phi)$ is
an arbitrary potential and $\mathcal{L}_{2n}$ are combinations of eigenvalues
which can be written in elementary symmetric polynomials, namely

$\displaystyle\mathcal{L}_{2}$ $\displaystyle\propto$ $\displaystyle
e_{1}(\lambda_{1}^{2},\lambda_{2}^{2},\dots,\lambda_{d}^{2})~{},$
$\displaystyle\mathcal{L}_{4}$ $\displaystyle\propto$ $\displaystyle
e_{2}(\lambda_{1}^{2},\lambda_{2}^{2},\dots,\lambda_{d}^{2})~{},$
$\displaystyle\vdots$ $\displaystyle\mathcal{L}_{2d}$ $\displaystyle\propto$
$\displaystyle
e_{d}(\lambda_{1}^{2},\lambda_{2}^{2},\dots,\lambda_{d}^{2})~{}.$ (8)

In static four dimensional spacetime case ($d=3$), this gives the generalized
Skyrme model which is just the traditional model from Skyrme added with BPS-
Skyrme term
($\mathcal{L}_{6}\propto\lambda_{1}^{2}\lambda_{2}^{2}\lambda_{3}^{2}$)
proposed in [20].
By following the prescription given in four dimensional spacetime model we can
construct $\mathcal{L}$ using

$\displaystyle\mathcal{L}_{2n}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{(n!)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{n}}_{i_{n}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{n}}_{j_{n}]}g^{i_{1}j_{1}}\dots
g^{i_{n}j_{n}}~{},$ (9)

so that the Lagrangian of the model mentioned above can be written as

$\mathcal{L}=-\gamma_{0}V-\sum_{n=1}^{d}\frac{\gamma_{n}}{\left(n!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{n}}_{i_{n}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{n}}_{j_{n}]}g^{i_{1}j_{1}}\dots
g^{i_{n}j_{n}}~{},$ (10)

where we have defined
$\phi^{a}_{i}=\partial_{i}\phi^{a}=\frac{\partial\phi^{a}}{\partial x^{i}}$
for simplicity. This is just an $O(d+1)$ model as the generalization proposed
in [23] if we take the potential to be $V=1-\phi^{0}$. As the summary, the
action functional for $d+1$ dimensional Skyrme model takes the form

$\mathcal{S}_{d+1}=-\int
d^{d+1}x\sqrt{-g}\left[\gamma_{0}V+\sum_{n=1}^{d}\frac{\gamma_{n}}{\left(n!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{n}}_{i_{n}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{n}}_{j_{n}]}g^{i_{1}j_{1}}\dots
g^{i_{n}j_{n}}\right]~{}.$ (11)

In fact, the action functional (11) does accommodate time dependencies on the
fields, but we restrict our model to be spherically symmetric and static.
Because the model is static, then $\phi^{a}_{0}=0$ is implied for every
possible value of $a$ which means that the strain tensor effectively lives on
the $d$ dimensional spatial manifold. This fact reduces the representation of
strain tensor into a $d\times d$ matrix, hence, consistent with the
construction for $\phi:\mathbb{R}^{d}\rightarrow S^{d}$ we have shown above.
The topological charge density of this higher dimensional spacetime model is
given by (up to a normalization constant. See [25, 23])

$\rho\propto\varepsilon^{i_{1}\dots i_{d}}\varepsilon^{a_{1}\dots
a_{d+1}}\phi^{a_{1}}_{i_{1}}\dots\phi^{a_{d}}_{i_{d}}\phi^{d+1}~{}.$ (12)

In the next subsection, we cast the above model into the combination of the
strain tensor $D$’s eigenvalue scheme using so called the generalized hedgehog
ansatz. The field equations for Skyrmion can be found by substituting the
corresponding ansatz to the dynamical equations of the Skyrme model came from
the variation of $\phi^{a}$ on the Skyrme action functional, namely

$\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\left[\nabla_{k}\left(\sum_{n=1}^{d}\frac{2\gamma_{n}}{\left(n!\right)^{2}}g^{i_{1}j_{1}}\dots
g^{i_{n}j_{n}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{n}}_{i_{n}]}\frac{\partial}{\partial\phi^{b}_{k}}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{n}}_{j_{n}]}\right)-\frac{\partial}{\partial\phi^{b}}\gamma_{0}V\right]=0~{}.$
(13)

The dynamical field equations above can be simplified by exploiting the Levi-
Civita symbol and strain tensor, which lead to the following field equations

$\displaystyle\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\left[\sum_{n=1}^{d}\frac{2n\gamma_{n}}{n!(d-n)!}\nabla_{j}\left(\phi^{b,i}H_{i}^{j}(n,d)\right)-\frac{\partial}{\partial\phi^{b}}\gamma_{0}V\right]=0~{},$
(14)

Here, we define a new tensor

$H_{i}^{j}(n,d)\equiv\varepsilon_{i~{}i_{2}\dots i_{n}k_{n+1}\dots
k_{d}}\varepsilon^{j~{}j_{2}\dots j_{n}k_{n+1}\dots
k_{d}}\prod_{m=2}^{n}D_{j_{m}}^{i_{m}}$ (15)

to simplify the computational process, with $D_{j}^{i}$s are components of
strain tensor $D$. $H_{i}^{j}(n,d)$ has special features, namely, it is a
diagonal tensor if the strain tensor is diagonal. This is the desired property
which are going to be exploited more in the calculation of field equations for
each submodels.

###  II.1 Spacetime Background and the suitable ansatz

It is necessary to choose in which coordinate the calculations are done
because different coordinate requires different ansatz to reduce the
complexity of the resulting field equation (without any loss of generality).
Here, we consider polyspherical coordinate system222This coordinate system is
an extended version of spherically symmetric Euclidean space that allows more
than one independent rotational planes. The details of polyspherical
coordinate system’s properties and the way we can construct such
parametrization for $\mathbb{R}^{d}$ are given in a more detailed manner in
Appendix A. with $N=\left\lfloor\frac{d}{2}\right\rfloor$ azimuthal coordinate
$\varphi_{i}$’s, namely

$ds^{2}=-dt^{2}+dr^{2}+\sum\limits_{i=1}^{N}r^{2}\left(d\mu_{i}^{2}+\mu_{i}^{2}d\varphi_{i}^{2}\right)$
(16)

with $\sum_{i}\mu_{i}^{2}=1$ for even $d$, and

$ds^{2}=-dt^{2}+dr^{2}+r^{2}d\mu_{0}^{2}+\sum\limits_{i=1}^{N}r^{2}\left(d\mu_{i}^{2}+\mu_{i}^{2}d\varphi_{i}^{2}\right)$
(17)

with $\mu_{0}^{2}+\sum_{i}\mu_{i}^{2}=1$ for odd $d$. The $\mu_{i}$ are called
cosine coordinates and we choose to reparameterize them as follows

$\displaystyle\mu_{0}$ $\displaystyle=$ $\displaystyle\cos\theta_{0}~{},$
$\displaystyle\mu_{1}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\cos\theta_{1}~{},$ $\displaystyle\mu_{2}$
$\displaystyle=$ $\displaystyle\sin\theta_{0}\sin\theta_{1}\cos\theta_{2}~{},$
$\displaystyle\vdots$ $\displaystyle\mu_{N-1}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\sin\theta_{1}\dots\cos\theta_{N-1}~{},$
$\displaystyle\mu_{N}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\sin\theta_{1}\dots\sin\theta_{N-1}~{}.$ (18)

the even $d$ case is simply recovered by choosing $\theta_{0}=\frac{\pi}{2}$
and it is straightforward to prove that the constraint
$\mu_{0}^{2}+\sum_{i}\mu_{i}^{2}=1$ is satisfied.

The spherically symmetric Skyrme field ansatz is a coordinate system on the
target space, $S^{d}$, given by

$\phi=(\cos\xi,\vec{n}\sin\xi)$ (19)

with $\vec{n}\cdot\vec{n}=n^{a}n^{a}=1$ and $\xi\equiv\xi(r)$ interpolates
between $\pi$ and zero. In this form of ansatz we have the following
expression of strain tensor (in matrix form)

$D=\frac{\sin^{2}\xi}{r^{2}}M+\left(\xi^{\prime}{}^{2}-\frac{\sin^{2}\xi}{r^{2}}\right)\frac{\tilde{\textbf{x}}\tilde{\textbf{x}}^{T}}{r^{2}}$
(20)

where components of matrix $M$ is given by

$M_{ij}=\frac{\partial\tilde{x}^{a}}{\partial
x^{i}}\frac{\partial\tilde{x}^{a}}{\partial x^{j}}$ (21)

with $x^{i}$s are coordinates in the standard Euclidean space, defined in (85)
for even $d$ and (90) for odd $d$, and we have defined a new vector
$\tilde{\textbf{x}}\equiv r\vec{n}$, hence $\tilde{x}^{a}$ represents
component of $\tilde{\textbf{x}}$. here, $i$ runs from $1$ to $d$,
corresponding to spatial coordinate only. In order to diagonalize the
resulting strain tensor of this map we need to employ the same symmetry
features of space-time spherical submanifold’s coordinate to the target space
coordinate by choosing $\vec{n}$ field to be parameterized by the cosine and
azimuthal coordinates. By exploiting the normalization constraint of
$\mu_{i}$’s we arrive at the following ansatz of $\vec{n}$

$\vec{n}_{\text{even}}=\left(\mu_{1}\cos(g_{1}),~{}\mu_{1}\sin(g_{1}),~{}\dots,~{}\mu_{N}\cos(g_{N}),~{}\mu_{N}\sin(g_{N})\right)~{}.$
(22)

for even $d$ and

$\vec{n}_{\text{odd}}=\left(\mu_{0},~{}\mu_{1}\cos(g_{1}),~{}\mu_{1}\sin(g_{1}),~{}\dots,~{}\mu_{N}\cos(g_{N}),~{}\mu_{N}\sin(g_{N})\right)~{}.$
(23)

for odd $d$ where we have defined the functions $g_{i}\equiv n_{i}\varphi_{i}$
with $n_{i}\in\mathbb{Z}$ are the winding numbers of each rotational plane.
The constants $n_{i}$ is crucial to construct Skyrmion with non-unit
topological degree because substitution of (22) and (23) to (12) gives the
following topological charge

$B=\frac{1}{C_{d}}\int\sqrt{-g}\rho~{}d^{d}x=n_{1}n_{2}\dots n_{N}$ (24)

where $C_{d}$ is a normalization constant, given that $B$ must be equal to one
when we consider a unit-topological case with all winding numbers equal to
one. We can see that the well-known property of Skyrmion topological charge is
satisfied where $B\in\mathbb{Z}$. In fact, this is inevitable because we
consider a model which belongs to the homotopy class
$\pi_{d}(S^{d})=\mathbb{Z}$ [3]. It is important to note that $n_{i}$s
themselves are not topological invariants and only their product in the form
of topological charge is a topological invariant. Here, $n_{i}$s only describe
how we parametrize topological charge in poly spherical coordinate.

We can view the static Skyrme field ansatz as a map $\mathbb{R}^{d}\rightarrow
S^{d}$ given by

$(r,\theta_{0},\dots,\theta_{N-1},\varphi_{1},\dots,\varphi_{N})\rightarrow(\xi,\theta_{0},\dots,\theta_{N-1},n_{1}\varphi_{1},\dots,n_{N}\varphi_{N})~{}.$
(25)

Hence, we can calculate the Jacobian matrix and its corresponding strain
tensor explicitly. The resulting eigenvalues of strain tensor within this
choice of ansatz are

$\lambda_{2N+1}^{2}=\xi^{\prime}{}^{2},~{}\lambda_{N+1}^{2}=\dots=\lambda_{2N}^{2}=\frac{\sin^{2}\xi}{r^{2}},\text{
and}~{}\lambda_{k}^{2}=n_{k}^{2}\frac{\sin^{2}\xi}{r^{2}}~{}\text{for
}k=1,2,\dots,N~{},$ (26)

for odd $d$ cases and

$\lambda_{2N}^{2}=\xi^{\prime}{}^{2},~{}\lambda_{N+1}^{2}=\dots=\lambda_{2N-1}^{2}=\frac{\sin^{2}\xi}{r^{2}},\text{
and}~{}\lambda_{k}^{2}=n_{k}^{2}\frac{\sin^{2}\xi}{r^{2}}~{}\text{for
}k=1,2,\dots,N~{},$ (27)

for even $d$ cases. This is the expected form of eigenvalues for higher
dimensional case which reduce to the one found in $d=4$ case [32]. From here,
calculating the explicit form of effective Lagrangian by combinations of
strain tensor eigenvalues is straightforward via equation (8). It is found
that in five dimensional Skyrme model, there exist an interesting property,
known as topological degeneracy, where some multi-soliton states of BPS
skyrmion could have different energy states with the same value of topological
charge as a result of several possible configurations of winding numbers. This
property bring up a question on how winding numbers $n_{i}$ fundamentally
affect skyrmion properties and the eigenvalues we have for this arbitrary
dimensional Skyrme model that depend directly to the winding numbers is a good
starting point to do so.

###  II.2 The BPS Submodels

From equations (26), (27) and (8) we can construct the effective Skyrme
Lagrangian of the form

$\mathcal{L}_{eff}=-\left[(\xi^{\prime})^{2}\sum_{n=1}^{d}\gamma_{n}\frac{\sin^{2(n-1)}\xi}{r^{2(n-1)}}{}_{d}K_{n-1}+\gamma_{0}V+\sum_{n=1}^{d}\frac{\gamma_{n}\sin^{2n}\xi}{r^{2n}}{}_{d}K_{n}\right]~{},$
(28)

where we have defined functions of winding number
${}_{d}K_{n}\equiv{}_{d}K_{n}(n_{1},\dots n_{N})~{},$ which also depend on
spatial dimensions $d$, given by

  * •

For odd $d$ cases

${}_{d}K_{n}=\sum_{i=0}^{n}C^{(d-1)/2}_{n-i}e_{i}(n_{1}^{2},n_{2}^{2},\dots,n_{(d-1)/2}^{2})=\sum_{i=0}^{n}C^{N}_{n-i}e_{i}(n_{1}^{2},n_{2}^{2},\dots,n_{N}^{2})~{}.$
(29)

  * •

For even $d$ cases

${}_{d}K_{n}=\sum_{i=0}^{n}C^{(d/2)-1}_{n-i}e_{i}(n_{1}^{2},n_{2}^{2},\dots,n_{d/2}^{2})=\sum_{i=0}^{n}C^{N-1}_{n-i}e_{i}(n_{1}^{2},n_{2}^{2},\dots,n_{N}^{2})~{}.$
(30)

The construction of ${}_{d}K_{n}$ from eigenvalues of strain tensor is
straightforward but could be very cumbersome for higher values of $d$, hence
we carry the calculations only for some special cases where their explicit
form is necessary. These functions satisfies

${}_{d}K_{n}\geq C^{d}_{n}-C^{d-1}_{n-1}~{}.$ (31)

with equality is achieved only if all winding numbers are equal to one,
$n_{1}=\dots=n_{N}=1$. To find the BPS type solution of this model we can use
BPS Lagrangian method for Skyrme model which can be found on [33, 34]. Since
our metric satisfies $\sqrt{-g}\propto r^{d-1}$, the method suggests a BPS
Lagrangian $\mathcal{L}_{BPS}=-\frac{Q}{r^{d-1}}\xi^{\prime}$ with $Q\equiv
Q(\xi)$ does not depend explicitly on the coordinates. At the BPS limit
$\mathcal{L}_{eff}-\mathcal{L}_{BPS}=0$ we have

$\left[\sum_{n=1}^{d}\gamma_{n}\frac{\sin^{2(n-1)}\xi}{r^{2(n-1)}}{}_{d}K_{n-1}\right](\xi^{\prime})^{2}-\left[\frac{Q}{r^{d-1}}\right]\xi^{\prime}+\left[\gamma_{0}V+\sum_{n=1}^{d}\frac{\gamma_{n}\sin^{2n}\xi}{r^{2n}}{}_{d}K_{n}\right]=0~{},$
(32)

which is the quadratic equation of $\xi^{\prime}$. Solving (32) for
$\xi^{\prime}$ and then imposing condition for which we have only single
solution, it gives

$\frac{Q^{2}}{4r^{2(d-1)}}-\left[\sum_{n=1}^{d}\gamma_{n}\frac{\sin^{2(n-1)}\xi}{r^{2(n-1)}}{}_{d}K_{n-1}\right]\left[\gamma_{0}V+\sum_{n=1}^{d}\frac{\gamma_{n}\sin^{2n}\xi}{r^{2n}}{}_{d}K_{n}\right]=0~{}.$
(33)

Equation (33) is nothing but a polynomials in $r$, thus, by setting every
corresponding coefficient to be zero, we have a set of equations to determine
each $\gamma_{n}$ that satisfies the BPS limit. These equations are

$\displaystyle
0={}_{d}K_{k-1}\gamma_{0}\gamma_{k}V+\sum_{n+m=k}\gamma_{n}\gamma_{m}~{}{}_{d}K_{n-1}~{}{}_{d}K_{m}~{}~{}~{}$
$\displaystyle\text{For }k=2,3,\dots,d~{},$ (34) $\displaystyle
0=\sum_{n+m=k}\gamma_{n}\gamma_{m}~{}{}_{d}K_{n-1}~{}{}_{d}K_{m}~{}~{}~{}$
$\displaystyle\text{For }k=d+1,d+2,\dots,2d~{},$ (35)
$\displaystyle\sin^{2(d-1)}\xi\left[\gamma_{0}\gamma_{d}V+\sum_{n+m=d}\gamma_{n}\gamma_{m}~{}{}_{d}K_{n-1}~{}{}_{d}K_{m}\right]$
$\displaystyle=$ $\displaystyle\frac{Q^{2}}{4}~{}.$ (36)

Since ${}_{d}K_{n}$ are positive definite quantities, we can deduce that for
all $d\in\mathbb{N}$ there exists only one BPS submodel with
$\gamma_{1}=\dots=\gamma_{d-1}=0$. This submodel is a higher dimensional
extension of BPS-Skyrme model in four spacetime dimensions [20]. On the other
side, there is another BPS submodel if we restrict our case for only even $d$
with coupling constants
$\gamma_{0}=\dots=\gamma_{\frac{d}{2}-1}=\gamma_{\frac{d}{2}+1}=\dots=\gamma_{d}=0$,
i.e. only $\gamma_{\frac{d}{2}}$ is non-zero. This submodel is also a higher
dimensional extension of the model in five spacetime dimensions introduced in
[23]. Both submodels in $d+1$ dimensional spacetime are constructed from the
scale invariant hierarchy and the scalar norm of the topological current (or
topological charge density). Solutions and properties of both BPS submodels
are discussed in the next sections.

##  III BPS Bounds and BPS Type Solutions Via Bogomolnyi Equations

###  III.1 BPS Bounds of The First and Second Submodel

Firstly, let us consider the first submodel that exist for any value of
spatial dimension $d\in\mathbb{N}$. In this case we only have one BPS
submodels as mentioned above. The Lagrangian of this submodel is given by

$\mathcal{L}_{I}=-\gamma_{0}V-\frac{\gamma_{d}}{\left(d!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d}}_{i_{d}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d}}_{j_{d}]}g^{i_{1}j_{1}}\dots
g^{i_{d}j_{d}}~{},$ (37)

which gives the following equation of motions for the fields

$\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\left[\nabla_{k}\left(\frac{2\gamma_{d}}{\left(d!\right)^{2}}g^{i_{1}j_{1}}\dots
g^{i_{d}j_{d}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d}}_{i_{d}]}\frac{\partial}{\partial\phi^{b}_{k}}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d}}_{j_{d}]}\right)-\frac{\partial}{\partial\phi^{b}}\gamma_{0}V\right]=0~{}.$
(38)

The BPS inequality can be found by introducing the energy functional of the
static case which satisfies $E=\int\sqrt{-g}dx^{d}\mathcal{L}$. Recall that
the second term in (37) is just the determinant of the strain tensor, thus our
energy functional takes the following form

$E=\int\sqrt{-g}d^{d}x\left[\gamma_{d}\text{det}D+\gamma_{0}V\right]=\int\sqrt{-g}d^{d}x\left[\left(\sqrt{\gamma_{d}\text{det}D}\mp\sqrt{\gamma_{0}V}\right)^{2}\pm
2\sqrt{\gamma_{0}\gamma_{d}V\text{det}D}\right]~{}.$ (39)

In the second term we can transform the spacetime invariant measure
$d^{d}x\sqrt{-g}$ to the target space invariant measure by using
$\sqrt{\text{det}D}\sqrt{-g}~{}d^{d}x=\sqrt{-g}\left|\det{J}\right|d^{d}x=\sqrt{h}~{}d^{d}\eta$
where $\eta^{a}$’s are the coordinates and $h$ is the metric tensor
determinant of the target manifold. The integration on the entire target
manifold is actually $B$ fold of integration on the entire $S^{d}$, where $B$
is the topological degree of $\phi$, hence we state the resulting BPS
inequality

$E\geq
2B\sqrt{\gamma_{0}\gamma_{d}}\text{Vol}S^{d}\left<\sqrt{V}\right>_{S^{d}}~{},$
(40)

with $\text{Vol}S^{d}$ is the total volume of compact manifold $S^{d}$ and
$\left<a\right>_{S^{d}}$ denotes the average of function $a$ on $S^{d}$. The
straightforward transformation of base space invariant measure (Volume form)
to the target space measure implies that the derivative term of this submodel
is invariant under all volume-preserving diffeomorphisms on $S^{d}$, which has
been shown in [35]. Because this model belong to the same class of BPS Skyrme
model in [36] and generally (for higher dimensions) in [35], then we can
conclude that the integrability is guaranteed. The nonzero potential term
itself provides a symmetry braking into a smaller group which leaves the $\xi$
unchanged, represented by the constraint $\vec{n}\cdot\vec{n}=1$. This
restricts the symmetries further for only volume-preserving diffeomorphisms on
$S^{d-1}$, which again, implies that this model is a straightforward extension
of the model in [19].

The BPS bound which restricts the lower bound of energy to be linear in
topological degree is the main feature of this model. This BPS bound given in
(40) strongly constrains the energy spectrum of this type of soliton which
belongs to the class of $D$-brane soliton. The bound is satisfied by
generalized Skyrme as well since adding the other terms only makes the energy
higher. Furthermore, the result (40) implies that for every $d\in\mathbb{N}$
there is at least one class of Skyrmion (or a submodel) that is non-
interacting, provided that the BPS bound is saturated. This fact not only
applies to the compacton solution but also the regular ones, as we show in a
later subsection that regular solutions exist and are well-behaved.

Now, let us consider the special case of even $d$. For this case there exists
two type of BPS submodels. The first one is from the general $d\in\mathbb{N}$
case which has been previously described above and the second one is obtained
by taking only $\gamma_{\frac{d}{2}}$ to be non-zero. The Lagrangian of the
second submodel is given by

$\mathcal{L}_{II}=-\frac{\gamma_{d/2}}{\left((d/2)!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d/2}}_{i_{d/2}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d/2}}_{j_{d/2}]}g^{i_{1}j_{1}}\dots
g^{i_{d/2}j_{d/2}}$ (41)

that implies the following field equations (Euler-Lagrange equations)

$\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\nabla_{k}\left(g^{i_{1}j_{1}}\dots
g^{i_{d/2}j_{d/2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d/2}}_{i_{d/2}]}\frac{\partial}{\partial\phi^{b}_{k}}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d/2}}_{j_{d/2}]}\right)=0$
(42)

Here, we take a slightly different approach to find the BPS inequality of this
submodel. We know that from the construction of Lagrangian terms via strain
tensor eigenvalues we have

$\mathcal{L}_{d}\propto
e_{d/2}(\lambda_{1}^{2},\dots,\lambda_{d}^{2})\propto\lambda_{1}^{2}\dots\lambda_{\frac{d}{2}}^{2}+\dots+\lambda_{\frac{d}{2}+1}^{2}\dots\lambda_{d}^{2}$
(43)

There are $C^{d}_{d/2}$ terms in (43), hence by employing arithmetic-geometric
means inequality, which reads $\sum_{i}w_{i}f_{i}\geq\prod_{i}f_{i}^{w_{i}}$
with $\sum_{i}w_{i}=1$, to (43) with identical $w_{i}=\frac{1}{C^{d}_{d/2}}$
for every term we can proceed as follow

$\displaystyle\lambda_{1}^{2}\dots\lambda_{\frac{d}{2}}^{2}+\dots+\lambda_{\frac{d}{2}+1}^{2}\dots\lambda_{d}^{2}$
(44) $\displaystyle=$ $\displaystyle
C^{d}_{d/2}\left(\frac{1}{C^{d}_{d/2}}\lambda_{1}^{2}\dots\lambda_{\frac{d}{2}}^{2}+\dots+\frac{1}{C^{d}_{d/2}}\lambda_{\frac{d}{2}+1}^{2}\dots\lambda_{d}^{2}\right)$
$\displaystyle\geq$ $\displaystyle
C^{d}_{d/2}(\lambda_{1}^{2})^{\frac{C^{d-1}_{(d/2)-1}}{C^{d}_{d/2}}}\dots(\lambda_{d}^{2})^{\frac{C^{d-1}_{(d/2)-1}}{C^{d}_{d/2}}}$
$\displaystyle=$ $\displaystyle
C^{d}_{d/2}\lambda_{1}\dots\lambda_{d}=C^{d}_{d/2}\sqrt{\text{det}D}$

From here, it is straightforward to conclude that this BPS submodel satisfies
a BPS inequality of the form

$E\geq B~{}\gamma_{d/2}C^{d}_{d/2}\text{vol}S^{d}$ (45)

with BPS limit saturation happens if and only if
$\lambda_{1}^{2}=\dots=\lambda_{d}^{2}$, i.e. all the eigenvalues are
identical. This is the self-duality condition in terms of strain tensor
eigenvalues. In fact, this property can be derived directly from the usual
self-dual equation for $O(d+1)$ Skyrme-sigma (see, [25]) model given by

$\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d/2}}_{i_{d/2}]}=\frac{1}{((d/2)!)^{2}}\varepsilon_{i_{1}\dots
i_{d}}\varepsilon^{a_{1}\dots
a_{d+1}}\phi^{a_{(d/2)+1}}_{[j_{(d/2)+1}}\dots\phi^{a_{d}}_{j_{d}]}\phi^{a_{d+1}}g^{i_{(d/2)+1}j_{(d/2)+1}}\dots
g^{i_{d}j_{d}}.$ (46)

that reduces into $\frac{C^{d}_{(d/2)}}{2}$ equations which can be expressed
as follows

$\lambda_{k_{1}}\dots\lambda_{k_{d/2}}\pm\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{d}}=0,$
(47)

with all the indices $k_{i}$ must take different values,
$k_{i}\in{1,2,\dots,d}$. We can show that equation (47) is equivalent to

$\lambda_{1}^{2}=\dots=\lambda_{d}^{2},$ (48)

as follows. Let us choose two arbitrary $\lambda_{k_{i}}$ from the left-hand-
side and $\lambda_{k_{j}}$ from the right-hand-side of
$\lambda_{k_{1}}\dots\lambda_{k_{i}}\dots\lambda_{k_{d/2}}=\mp\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{j}}\dots\lambda_{k_{d}}$.
By switching $\lambda_{k_{i}}$ to the right-hand-side and $\lambda_{k_{j}}$ to
the left-hand-side we have
$\lambda_{k_{1}}\dots\lambda_{k_{j}}\dots\lambda_{k_{d/2}}=\mp\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{i}}\dots\lambda_{k_{d}}$,
and since all $k_{i}$s are different, this equation is also allowed by (47).
Dividing the first equation with the switched equation gives

$\frac{\lambda_{k_{i}}}{\lambda_{k_{j}}}=\frac{\lambda_{k_{j}}}{\lambda_{k_{i}}},$
(49)

that implies $\lambda_{k_{i}}^{2}=\lambda_{k_{j}}^{2}$. Since $k_{i}$ and
$k_{j}$ are arbitrary, the equation applies to all $k_{i}\in{1,2,\dots,d}$.
This shows the equivalence between self-dual equation and BPS limit saturation
conditions in terms of strain tensor eigenvalues.

Lastly, it is beneficial to show that (45) is indeed the true lower bound of
the energy. To do this, consider the following squared quantity from the self-
duality equation

$\left(\lambda_{k_{1}}\dots\lambda_{k_{d/2}}\pm\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{d}}\right)^{2}=\left(\lambda_{k_{1}}\dots\lambda_{k_{d/2}}\right)^{2}+\left(\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{d}}\right)^{2}\pm
2\lambda_{1}\lambda_{2}\dots\lambda_{d},$ (50)

where we have used the fact that all $k_{i}$s must take different values that
leads to the last term in the right-hand side become two times the product of
all $\lambda_{i}$s. Both the $\pm$ sign in the left-hand and right-hand sides
take the same value simultaneously with positive sign represents anti-soliton
(with negative topological degree), and vice versa. Thus, we proceed with
negative sign for the above equation since we are only dealing with soliton
(with positive topological degree). Next, summing all of the above
$\frac{C^{d}_{d/2}}{2}$ possible terms gives

$\sum_{k_{1}\neq k_{2}\neq\dots
k_{d}}\left(\lambda_{k_{1}}\dots\lambda_{k_{d/2}}-\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{d}}\right)^{2}=e_{d/2}(\lambda_{1}^{2},\dots,\lambda_{d}^{2})-C^{2}_{d/2}\lambda_{1}\lambda_{2}\dots\lambda_{d}.$
(51)

Notice that the first term in the right-hand side is proportional to
$\mathcal{L}_{d}$, thus we can write the total energy as

$E=\gamma_{d/2}\int\sqrt{-g}\left[\sum_{k_{1}\neq k_{2}\neq\dots
k_{d}}\left(\lambda_{k_{1}}\dots\lambda_{k_{d/2}}-\lambda_{k_{(d/2)+1}}\dots\lambda_{k_{d}}\right)^{2}+C^{2}_{d/2}\lambda_{1}\lambda_{2}\dots\lambda_{d}\right]~{}d^{d}x.$
(52)

We can see that if we impose the self-dual equation to the above equation then
the summation vanishes, which leave us only the lower-bound term that have
been previously found via arithmetic-geometric inequality.

###  III.2 BPS Solutions

####  III.2.1 The First Submodel

The first BPS submodel within the ansatz of this work for the static case
($\omega=0$) is given by

$E=\int\sqrt{-g}\left[\gamma_{d}B^{2}\frac{\sin^{2(d-1)}\xi}{r^{2(d-1)}}\xi^{\prime
2}+\gamma_{0}V\right]d^{d}x~{}$ (53)

where we have substituted the value of topological charge, $B=n_{1}n_{2}\dots
n_{N}$. The BPS limit is satisfied when
$\sqrt{\gamma_{d}\text{det}D}\mp\sqrt{\gamma_{0}V}=0$ which leads to the
Bogomolny equation of this submodel, given by

$B\frac{\sin^{d-1}\xi}{r^{d-1}}\xi^{\prime}\mp\sqrt{\frac{\gamma_{0}}{\gamma_{d}}}\sqrt{V}=0~{}.$
(54)

which implies that the cut of radius now scales as $r_{cut-off}\propto
B^{\frac{1}{d}}$.333Similar to the property of BPS Skyrme model in $d=3$
spatial dimensions This property of cut-off radius is expected from the
scaling property of the model itself. This can be physically interpreted as
the volume of the first type of BPS Skyrmions is proportional to its
topological degree.

The Bogomolnyi equation for this model gives an identical solution with the
$B=1$ solution in (56) and (57) if we rescale the radial coordinate as
$r\rightarrow r/B^{\frac{1}{d}}$, hence, there is no new qualitative feature
of this non-unit topological charge case. Thus, from this point on, we
restrict ourselves to the case of $B=1$ without any loss of generality which
satisfies the following Bogomolnyi equation

$\frac{\sin^{d-1}\xi}{r^{d-1}}\xi^{\prime}\mp\sqrt{\frac{\gamma_{0}}{\gamma_{d}}}\sqrt{V}=0~{}.$
(55)

We can immediately see that vacuum solution is the trivial solution of this
submodel. The non-trivial solution needs a specific form of potential, in this
work we use the pion-mass type potential to the power of $s$ which takes the
form $V=(1-\phi_{0})^{s}=(1-\cos\xi)^{s}$ with $s\in\mathbb{R}^{+}$.
Integrating equation (55) with mentioned pion-mass type potential gives

$\frac{\sin^{d-s}\left(\frac{\xi}{2}\right)}{d-s}~{}{}_{2}\text{F}_{1}\left(1-\frac{d}{2},~{}\frac{d-s}{2};~{}\frac{2+d-s}{2};~{}\sin^{2}\left(\frac{\xi}{2}\right)\right)\mp\sqrt{\frac{2^{s}\gamma_{0}}{d^{2}\gamma_{d}}}\left(\frac{r}{2}\right)^{d}+c_{d}=0~{},$
(56)

for $d\neq s$ or

$\frac{\cos^{d}\left(\frac{\xi}{2}\right)}{d}~{}{}_{2}\text{F}_{1}\left(1,~{}\frac{d}{2};~{}\frac{2+d}{2};~{}\cos^{2}\left(\frac{\xi}{2}\right)\right)\mp\sqrt{\frac{2^{d}\gamma_{0}}{d^{2}\gamma_{d}}}\left(\frac{r}{2}\right)^{d}+c^{\prime}_{d}=0~{},$
(57)

for special case where $d=s$ with ${}_{2}\text{F}_{1}$ is the hypergeometric
function while $c_{d}$ and $c^{\prime}_{d}$ are integration constants.
Equation (56) holds for all $d\in\mathbb{N}$ because the hypergeometric series
from ${}_{2}\text{F}_{1}$ always terminate at finite term. From equation (56)
we can find the cut -off value of radius $r$ in which $\xi$ takes vacuum
value. This cut-off radius is

$r_{cut-
off}=\begin{cases}2\left[\frac{d}{d-s}\sqrt{\frac{\gamma_{d}}{2^{s}\gamma_{0}}}\frac{\Gamma\left(\frac{d}{2}\right)\Gamma\left(\frac{2+d-s}{2}\right)}{\Gamma\left(d-\frac{s}{2}\right)}\right]^{\frac{1}{d}}&\text{for}~{}0<s<d\\\
\infty&\text{for}~{}s\geq d\end{cases}~{},$ (58)

which implies that Skyrmion distributed in a smaller radius for higher
dimensions and the solution becomes a regular solution for $s\geq d$. It is
also can be observed that the cut-off radius is asymptotic to 1 as $d$ goes to
infinity for compacton solutions. This property ensures the existence of
positive finite energy of compacton solutions in BPS Skyrme model for any
dimensions. Since equations (56) and (57) are transcendental in $\xi$ we use
numerical method to find the solutions for some special cases. The first case
is the ordinary pion-mass potential where $s=1$. This case is a higher
dimensional extension of the four-dimensional case in [20]. We can see in
Figure 1 that all solutions are of compacton type for $d\geq 2$, in agreement
with the special four-dimensional spacetime case.

Figure 1: Numerical solution for equation (55) with pion-mass type potential

The next special case is $s=d$ where the solution starts to become a regular
solution that extends to infinity. We demonstrate this feature in Figure 2
where we plot the solutions in $d=5$ spatial manifold for $s=1,2,\dots,8$. In
Figure 3 we show the solution profiles for some value of $d$. We can observe
that the solutions become closer as $d$ goes to infinity which is the same
feature we have found in the compacton solutions above. Furthermore, there are
no intersections between solutions with different value of $d$ and all
solutions are monotonically decreasing, i.e. $\xi_{d}\geq\xi_{d+1}\geq 0$ for
$r\in[0,\infty)$ with $\xi_{d}$ is the solution of $\xi(r)$ in $d$ spatial
dimension. From these mentioned features and further support by the fact that
energy for $d=1$ case is finite, we can conclude that the energy, in general,
is positive and finite for all values of $d$.

Figure 2: Numerical solutions for equation (55) with pion-mass type potential
to the power $s=1,2,\dots,8$ for $d=5$ spatial manifold. The red curves
represent $s=1,2,3,4$ cases, and the blue curves represent $s=5,6,7$ cases. We
can see that for higher $s$ the profile of $\xi$ decays slower as
$r\rightarrow\infty$ Figure 3: Numerical solutions for equation (55) with
pion-mass type potential to the power of $s=d$.

Now, consider the energy-momentum tensor of this submodel which is given by

$\displaystyle T_{ij}$ $\displaystyle=$ $\displaystyle
2\frac{d\gamma_{d}}{\left(d!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d-1}}_{i_{d-1}}\phi^{a_{d}}_{i]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d-1}}_{j_{d-1}}\phi^{a_{d}}_{j]}g^{i_{1}j_{1}}\dots
g^{i_{d-1}j_{d-1}}$ (59) $\displaystyle-
g_{ij}\gamma_{0}V-g_{ij}\frac{\gamma_{d}}{\left(d!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d}}_{i_{d}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d}}_{j_{d}]}g^{i_{1}j_{1}}\dots
g^{i_{d}j_{d}}$

Substituting the Bogomolnyi equation to (59) gives

$T_{00}=\varepsilon,~{}T_{0j}=T_{ij}=0~{}\text{for all values
of}~{}i,j~{}\text{except 0,}$ (60)

where $\varepsilon$ is the energy density. Here, we can observe that the
energy density does not have a simple expression in $r$, but the asymptotic
expression is useful in energy estimate. Suppose that we choose a large
positive real $R$ such that for $r>R$ the energy density can be approximated
by its asymptotic form, namely

$\varepsilon\approx\frac{2\gamma_{0}}{2^{s}}\xi^{2s}\approx\frac{2\gamma_{0}}{2^{s}}\left(2^{s/2}\sqrt{\frac{\gamma_{d}}{\gamma_{0}}}\frac{d}{s-d}\right)^{\frac{2s}{s-d}}\frac{1}{r^{\frac{2sd}{s-d}}}~{}.$
(61)

Note that we have used the asymptotic expression of (55) to find the
asymptotic solution of $\xi$ and we only consider the case with $s>d$ where
asymptotic vacuum is slower to be reached for higher $s$.444For the cases with
$s\in[1,d]$ the profiles decay faster as $r\rightarrow\infty$, thus they have
finite total energy if the cases with higher $s$ have finite total energy. If
$R$ is arbitrary, then the total energy is always smaller than the following
integral

$\displaystyle\int_{0}^{R}\sqrt{-g}\varepsilon~{}d^{d}x+\frac{2\gamma_{0}}{2^{s}}\left(2^{s/2}\sqrt{\frac{\gamma_{d}}{\gamma_{0}}}\frac{d}{s-d}\right)^{\frac{2s}{s-d}}~{}\text{Vol}S^{d-1}\int_{R}^{\infty}r^{d-1}\left(\frac{1}{r^{\frac{2sd}{s-d}}}\right)dr~{}.$
(62)

To have finite energy the second term above must be finite which implies that
$d-1-\frac{2sd}{s-d}<-1$ or, in other words, the total energy is finite for
every $s\in(d,\infty)$. In a conclusion, the BPS submodel with unit
topological degree given in (37) is equipped with the pion-mass potential to
the power $s$ and static Skyrme field ansatz has finite total energy which is
bounded from below by (40).

####  III.2.2 The Second Submodel

The second submodel has the following energy functional

$E=\gamma_{d/2}\int\sqrt{-g}\left[{}_{d}K_{(d/2)-1}\xi^{\prime
2}+{}_{d}K_{d/2}\frac{\sin^{2}\xi}{r^{2}}\right]\frac{\sin^{d-2}\xi}{r^{d-2}}~{}d^{d}x~{}.$
(63)

In general, the Bogomolnyi equation of this submodel is non-self-dual which
takes the form

$\xi^{\prime}\pm\sqrt{\frac{{}_{d}K_{d/2}}{{}_{d}K_{(d/2)-1}}}\frac{\sin\xi}{r}=0.$
(64)

We can see that this equation is, indeed, non self-dual since the self-duality
equation (46) is not satisfied for $d>2$ unless we take $B=1$, but the
solution of (64) is the skyrmion solution of second submodel which satisfies
the field equations as demonstrated in appendix B.2. Thus, the resulting
solutions are higher topological degree skyrmions that admit spherical
symmetry.

The Bogomolnyi equation of this second submodel satisfies (42) and can be
derived directly from effective lagrangian via BPS Lagrangian method. The
solution of (64) is given by

$\xi=2\arctan\left(\left(\frac{r_{0}}{r}\right)^{\sqrt{\frac{{}_{d}K_{d/2}}{{}_{d}K_{(d/2)-1}}}}\right)~{}.$
(65)

This solution admits topological degeneracy where different solutions and
energies exist for a single value of the topological degree.555We take the
positive signed solution of (64) in order to place the potential vacuum at
spatial infinity, resulting a Skyrmion with positive topological degree. In
fact, the negative signed solution also exist which is known as anti-Skyrmion
with quite similar properties with ordinary Skyrmion This is because we can
choose different configurations of $n_{1}$ to $n_{N}$ without changing the
value of the topological degree. In contrast with the topological degree which
is just a product of all winding numbers, $n_{1}$ to $n_{N}$, the solution and
energy highly depend on the value of ${}_{d}K_{n}$ which gives different
values for different configurations of winding numbers [32].

The energy momentum tensor of the second submodel is

$\displaystyle T_{ij}$ $\displaystyle=$
$\displaystyle\frac{d\gamma_{d/2}}{\left((d/2)!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d/2}}_{i]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d/2}}_{j]}g^{i_{1}j_{1}}\dots
g^{i_{(d/2)-1}j_{(d/2)-1}}$ (66) $\displaystyle-
g_{ij}\frac{\gamma_{d/2}}{\left((d/2)!\right)^{2}}\phi^{a_{1}}_{[i_{1}}\dots\phi^{a_{d/2}}_{i_{d/2}]}\phi^{a_{1}}_{[j_{1}}\dots\phi^{a_{d/2}}_{j_{d/2}]}g^{i_{1}j_{1}}\dots
g^{i_{d/2}j_{d/2}}~{}.$

If we substitute the first order field equation (64) the trace of energy-
momentum tensor of this submodel becomes

$T_{ij}g^{ij}=T_{0}^{0}=-\varepsilon=-2^{d+1}\gamma_{d/2}~{}{}_{d}K_{d/2}~{}\left(\frac{r^{\alpha-1}r_{0}^{\alpha}}{r_{0}^{2\alpha}+r^{2\alpha}}\right)^{d}~{},$
(67)

because the radial component of energy-momentum tensor vanishes, $T_{1}^{1}=0$
and the angular component, $T_{2}^{2},\dots,T_{d}^{d}$, are non-zero but the
summation of all angular component vanishes, i.e. $\sum_{i=2}^{d}T^{i}_{i}=0$.
In (67) we have defined
$\alpha\equiv\sqrt{\frac{{}_{d}K_{d/2}}{{}_{d}K_{(d/2)-1}}}$ to simplify the
expression. Now, consider the expression of total energy of this second
submodel which is given by

$E=2\gamma_{d/2}\sqrt{{}_{d}K_{(d/2)-1}~{}{}_{d}K_{d/2}}\text{Vol}S^{d}.$ (68)

We can see that this energy is higher than the BPS bound due to the non-self-
duality of the solution. As a result, there are possible repulsive
interactions between $B=1$ Skyrmions in a multi-soliton configuration where
the total topological degree is more than one, which is a similar property
found in $d=4$ case [32].

Interestingly, for $d=2$ the Bogomolnyi equation (64) becomes self-dual
because it satisfies

$\lambda_{1}^{2}=\xi^{\prime
2}=\lambda_{2}^{2}=n_{1}^{2}\frac{\sin^{2}\xi}{r^{2}}=\frac{{}_{2}K_{1}(n_{1})}{{}_{2}K_{0}(n_{1})}\frac{\sin^{2}\xi}{r^{2}}~{}.$
(69)

This lower dimensional situation cannot be found in a higher-dimensional case,
at least, within the same construction as has been shown here. Hence, $d=2$ is
the only case where we can find stable non-interacting BPS Skyrmions of the
second BPS submodel for the case of $B>1$. The stability of this higher
topological degree case is further supported by the fact that $d=2$ is the
only case where all components of the energy-momentum tensor vanishes except
the $0,0$ component, $T^{0}_{0}$. This came from the fact that $d=2$
dimensional space has only one azimuthal coordinate and no polar coordinate,
which leads to only one winding number. The second submodel of $d=2$ we have
discussed here is well-known as $O(3)$ sigma model in $2+1$ dimensional
spacetime and it has been proven in the previous section that this sigma model
is a BPS submodel of Baby-Skyrme model, proposed in [30, 31].

Now, let us consider the case where Bogomolnyi equation coincides with self-
duality equation.666The self-dual conditions in terms of strain tensor
eigenvalues is given by
$\lambda_{1}^{2}=\lambda_{2}^{2}=\dots=\lambda_{d}^{2}$. Substituting the
ansatz from (II.1) to the self-dual equation for $d>2$ gives
$\xi^{\prime}{}^{2}=\frac{\sin^{2}\xi}{r^{2}}=n_{1}^{2}\frac{\sin^{2}\xi}{r^{2}}=\dots=n_{N}^{2}\frac{\sin^{2}\xi}{r^{2}}$
that can only be satisfied for $B=1$. In this submodel, this can only be done
by taking

$\lambda_{1}^{2}=\lambda_{2}^{2}=\dots=\lambda_{d}^{2}~{},$ (70)

which lead to a $B=1$ BPS Skyrmion unless $d=2$. The Corresponding Bogomolny
equation becomes

$\xi^{\prime}\pm\frac{\sin\xi}{r}=0~{}.$ (71)

Features and solutions of this submodel are well-known and has been discussed
in detail in [23, 32], namely

$\xi=2\tan^{-1}\left(\frac{r_{0}}{r}\right)$ (72)

where $r_{0}$ is an arbitrary scaling coming from integration constant. This
BPS skyrmion represent another class of $D$-brane soliton, firstly proposed in
[37] which is constructed from a generalized Nambu-Goto action and is expected
to be scale invariant. This fact can be easily seen from the energy functional
(or static part of Lagrangian) which is invariant under spatial coordinate
transformation given by $\vec{x}\rightarrow\mu\vec{x}$.

The resulting BPS bound in this model given in (45), again, proves that in any
arbitrary $d$ we always have a $D$-brane soliton (or specifically, skyrmion)
which saturates the linear energy, but in a quite different manner from the
previous submodel above. Firstly, we could easily notice that the Lagrangian
does not obey the symmetry under all volume-preserving diffeomorphisms on
$S^{d}$ unless the self-duality condition is satisfied. For the case of
hedgehog ansatz, this can only be satisfied for $B=1$ skyrmion, implying a
spherical symmetry. Several studies have been done in $B>1$ hedgehog setups
and it is found that $B>1$ cases cannot saturate the bound (45), resulting in
a repulsive configuration that cannot be stable [32]. Secondly, in this
submodel there is no potential term to break the symmetry, thus the target
space symmetry of $d$ dimensional sphere is conserved, provided that self-
duality conditions are satisfied. It is also interesting that in hedgehog
ansatz, the resulting skyrmion profile does not depend on spatial dimensions
$d$ which implies perfectly identical solution and features of skyrmion.

By substituting the Bogomolnyi equation (71) to (66) we arrive at the
following results

$T_{00}=\varepsilon,~{}T_{0j}=T_{ij}=0~{}\text{for all values
of}~{}i,j~{}\text{except 0,}$ (73)

where $\varepsilon$ is the energy density of this second submodel, given by

$\varepsilon=\gamma_{d/2}C^{d}_{d/2}2^{d}\left(\frac{r_{0}}{r_{0}^{2}+r^{2}}\right)^{d}~{}.$
(74)

The finiteness of total energy is guaranteed because in this submodel we do
not have a potential term. We can observe that the energy density takes a non-
zero value at $r=0$, namely

$\varepsilon(r=0)=2^{d}C^{d}_{d/2}\frac{\gamma_{d/2}}{r_{0}^{d}}$ (75)

which depends only on the scaling $r_{0}$ and coupling constant
$\gamma_{d/2}$. This property implies that for $B=1$ skyrmions are distributed
close to the coordinate origin. Because dependency of $\gamma_{d/2}$ can be
easily omitted by scaling of the coordinate on the Lagrangian level, hence we
are left with only $r_{0}$ as the only arbitrary constant777Usually, the
integration constants of Bogomolnyi equation is not arbitrary due to the
specific boundary conditions, but in this special submodel the boundary
conditions are satisfied for all size modulus $r_{0}\in\mathbb{R}$ which is
the consequence of scale invariance..

##  IV Conclusions and Outlooks

We have considered the higher dimensional Skyrme model and deduced that we can
make a slight modification on generalized hedgehog ansatz by introducing
winding numbers $n_{1}$ to $n_{N}$ for every azimuthal coordinate to
accommodate higher values of topological degree. The winding numbers
themselves are non-topological but the product of all winding numbers is
proportional to the topological degree of the model. This form of ansatz is an
extension of similar forms found in $d=2$ and $d=4$ skyrmions [30, 23]. In
this work, we consider only the static case and we are going to address the
more restrictive time-dependent cases in future works.

We have shown via BPS Lagrangian method that, in general, there exist at most
two independent BPS Skyrme submodels within the static Skyrme field ansatz.
The first BPS submodel is the model with only $\mathcal{L}_{2d}$ term and
potential term which can be found in any $d\in\mathbb{N}$ dimensional spatial
manifold. The second BPS submodel is the model with only $\mathcal{L}_{d}$
term which can only be found in an even-dimensional spatial manifold. Both
have been proposed by different hierarchies in [25].

The first BPS submodel has both compacton and regular solutions for the pion-
mass potential to the power of $s$ and lowers the bound of energy which is
proportional to the topological degree as expected. This submodel is the
higher dimensional extension of the known BPS Skyrme model in four and five
spacetime dimensions [20, 32]. The profile of the solution of this submodel
highly depends on the form of potential and the dimensions of the target
space. We found that the total energy is always finite for the pion-mass
potential to the power of $s$ with arbitrary power $s$, and both compacton and
regular solutions are non-interacting if the BPS bound is saturated.

The second submodel only has a regular solution which is similar to the one
found in the five-dimensional spacetime case [23, 32]. The resulting
Bogomolnyi equation is similar for every value of $d$ up to a constant which
depends on both $d$ and the configuration of winding numbers. In contrast with
the Bogomolnyi equation, the energy density has a distinctive feature that
depends only on the topological degree. The energy density of the second
submodel is concentrated at the origin for $B=1$. Other important features of
this type of BPS Skyrmion are that the solution and energy functional of this
submodel are scale-invariant and BPS limit saturation happens if and only if
all strain tensor eigenvalues are identical, i.e.

$\lambda_{1}^{2}=\lambda_{2}^{2}=\dots=\lambda_{d}^{2}~{}.$ (76)

This is the self-duality condition in terms of strain tensor eigenvalues.

In a more general case where a wider set of symmetry is considered with a
higher topological degree, we expect that BPS features of the first submodel
are maintained because the key term, which is the $\mathcal{L}_{2d}$ term, is
just a product of all eigenvalues. This is further supported by the expression
of the energy which can only be altered if we choose another kind of compact
manifold other than $S^{d}$ as the target space. In contrast with the first
one, the second submodel needs to satisfy a strong condition where all
eigenvalues are the same. This might lead to smaller sets of possible target
spaces where BPS solutions exist for this submodel. It is interesting to point
out some, or all, target space symmetries that are compatible with this second
submodel, this topic will be addressed in future works.

There exists non-self-dual solutions of the second submodel for higher
topological degree but these solutions are found to be non-BPS since they do
not saturate the BPS bound given in (45). Some aspects of these $B>1$
Skyrmions will be discussed elsewhere.

##  Appendix A On The Polyspherical Coordinate

As we know that a standard unit $d-1$ dimensional sphere manifold is defined
as

$S^{d-1}=\\{\vec{x}\in\mathbb{R}^{d},\|\vec{x}\|=1\\},$ (77)

where $\mathbb{R}^{d}$ is a $d$ dimensional Euclidean manifold and
$\|\vec{x}\|$ is the Euclidean norm of $\mathbb{R}^{d}$, defined as

$\|\vec{x}\|\equiv\sqrt{\sum_{i=1}^{d}\left(x^{i}\right)^{2}}.$ (78)

We define a rotational plane to be a $\mathbb{R}^{2}=\mathbb{R}^{+}\times
S^{1}$ manifold equipped with coordinate system $\\{\rho,\varphi\\}$ and
metric $ds_{2}^{2}=d\rho^{2}+\rho^{2}d\varphi^{2}$. Thus, a $S^{d-1}$ has
$N=\left\lfloor\frac{d}{2}\right\rfloor$ independent rotational planes due to
the fact that we can decompose $\mathbb{R}^{d}$ into a product of some
$\mathbb{R}^{2}$s. For even $d$ the decomposition is given by

$\mathbb{R}^{d}=\underbrace{\mathbb{R}^{2}\times\mathbb{R}^{2}\times\dots\times\mathbb{R}^{2}}_{N=\frac{d}{2}},$
(79)

and for odd $d$ we have the following decomposition

$\mathbb{R}^{d}=\underbrace{\mathbb{R}^{2}\times\mathbb{R}^{2}\times\dots\times\mathbb{R}^{2}}_{N=\frac{d-1}{2}}\times\mathbb{R}.$
(80)

In a conclusion, we have been able to parametrize a $d$ dimensional Euclidean
manifold into patches of rotational plane coordinates
$ds_{2}^{2}=d\rho^{2}+\rho^{2}d\varphi^{2}$. Let us explain the details of
this parametrization for both odd and even $d$ cases in the subsections below.

###  A.1 Even $d$ Case

From the decomposition given in (79) we have the following metric for even $d$
case,

$\displaystyle
ds^{2}=\sum_{i=1}^{d}\left(dx^{i}\right)^{2}=\sum_{i=1}^{N}\left(d\rho_{i}^{2}+\rho_{i}^{2}d\varphi_{i}^{2}\right).$
(81)

with $N=\frac{d}{2}$. Suppose that we have $\mu_{i}$ such that,
$\rho_{i}\equiv r\mu_{i}$, relating the radial coordinate of one rotational
plane to the radial coordinate of full space, then we have the following
relation

$\displaystyle
d\rho_{i}^{2}+\rho_{i}^{2}d\varphi_{i}^{2}=\mu_{i}^{2}dr^{2}+r^{2}d\mu_{i}^{2}+r^{2}\mu_{i}^{2}d\varphi_{i}^{2}+2r\mu_{i}drd\mu_{i},$
(82)

with no summation for repeated indices. Recall, that
$r=\|\vec{x}\|=\sqrt{\sum_{i=0}^{N}\rho_{i}^{2}}=\sqrt{r^{2}\sum_{i=0}^{N}\mu_{i}^{2}}$,
hence, we need to impose $\sum_{i=1}^{N}\mu_{i}^{2}=1$ which leads to
$\sum_{i=0}^{N}\mu_{i}d\mu_{i}=0$. Substituting these properties of $\mu_{i}$
to the metric we have

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left(d\rho_{i}^{2}+\rho_{i}^{2}d\varphi_{i}^{2}\right)$
(83) $\displaystyle=$ $\displaystyle
dr^{2}\sum_{i=1}^{N}\mu_{i}^{2}+r^{2}\sum_{i=1}^{N}d\mu_{i}^{2}+r^{2}\sum_{i=1}^{N}\mu_{i}^{2}d\varphi_{i}^{2}+2rdr\sum_{i=1}^{N}\mu_{i}d\mu_{i}$
$\displaystyle=$ $\displaystyle
dr^{2}+r^{2}\left(\sum_{i=1}^{N}d\mu_{i}^{2}+\sum_{i=1}^{N}\mu_{i}^{2}d\varphi_{i}^{2}\right).$

This is the desired form of metric where we have $N$ azimuthal coordinates,
$\varphi_{i}$, each of them living in a rotational plane.

It is conventional to take $r\geq 0$ and $0\leq\varphi_{i}<2\pi$ for every
$i$. Thus, we only need to specify the form of $\mu_{i}$ in terms of polar
coordinates, $\theta_{i}$, that represent the relative embedding position of
rotational planes in the full space, $\mathbb{R}^{d}$. One of possible
parametrization for $\mu_{i}$ is given by

$\displaystyle\mu_{1}$ $\displaystyle=$ $\displaystyle\cos\theta_{1}~{},$
$\displaystyle\mu_{2}$ $\displaystyle=$
$\displaystyle\sin\theta_{1}\cos\theta_{2}~{},$ $\displaystyle\vdots$
$\displaystyle\mu_{N-1}$ $\displaystyle=$
$\displaystyle\sin\theta_{1}\dots\cos\theta_{N-1}~{},$ $\displaystyle\mu_{N}$
$\displaystyle=$ $\displaystyle\sin\theta_{1}\dots\sin\theta_{N-1}~{}.$ (84)

From the definition, we know that $\mu_{i}=\rho_{i}/r$ and $0\leq\rho_{i}\leq
r$. Thus, $0\leq\mu_{i}\leq 1$ which means that every polar coordinate must
satisfy $0\leq\theta_{i}\leq\frac{\pi}{2}$. This completes our parametrization
of $S^{d-1}$ in polyspherical coordinate for even $d$ case.

Obtaining the inverse transformation to the standard Euclidean coordinate
$x^{i}$ is straightforward. Firstly, we need to divide the $x^{i}$s into two
sets with one set contains coordinate of every horizontal axes of submanifold
$\mathbb{R}^{2}$ and the other set contains the vertical axes coordinate.
Thus, we have the following relations

$\displaystyle x^{i}$ $\displaystyle=$ $\displaystyle
r\mu_{i}\cos\varphi_{i},$ (85) $\displaystyle x^{2i}$ $\displaystyle=$
$\displaystyle r\mu_{i}\sin\varphi_{i}.$ (86)

where $i$ runs from $1$ to $N$.

###  A.2 Odd $d$ Case

We can apply the same method from previous subsection for decomposition (80)
with extra coordinate $\rho_{0}$ for submanifold $\mathbb{R}$ such that the
metric for odd $d$ case is given by

$\displaystyle
ds^{2}=\sum_{i=1}^{d}\left(dx^{i}\right)^{2}=d\rho_{0}^{2}+\sum_{i=1}^{N}\left(d\rho_{i}^{2}+\rho_{i}^{2}d\varphi_{i}^{2}\right).$
(87)

with $N=\frac{d-1}{2}$. The constraint for $\mu_{i}$ now becomes
$\mu_{0}^{2}+\sum_{i=1}^{N}\mu_{i}^{2}=1$, which leads to
$\mu_{0}d\mu_{0}+\sum_{i=1}^{N}\mu_{i}d\mu_{i}=0$. Again, by substituting the
constraint for $\mu_{i}$s and its implications to the metric, we arrive at the
following expression

$\displaystyle ds^{2}$ $\displaystyle=$ $\displaystyle
r^{2}d\mu_{0}^{2}+2r\mu_{0}drd\mu_{0}+\mu_{0}^{2}dr^{2}+\sum_{i=1}^{N}\left(d\rho_{i}^{2}+\rho_{i}^{2}d\varphi_{i}^{2}\right)$
(88) $\displaystyle=$ $\displaystyle
dr^{2}\left(\mu_{0}^{2}+\sum_{i=1}^{N}\mu_{i}^{2}\right)+r^{2}\left(d\mu_{0}^{2}+\sum_{i=1}^{N}d\mu_{i}^{2}\right)+r^{2}\sum_{i=1}^{N}\mu_{i}^{2}d\varphi_{i}^{2}+2rdr\left(\mu_{0}d\mu_{0}+\sum_{i=1}^{N}\mu_{i}d\mu_{i}\right)$
$\displaystyle=$ $\displaystyle
dr^{2}+r^{2}d\mu_{0}^{2}+r^{2}\left(\sum_{i=1}^{N}d\mu_{i}^{2}+\sum_{i=1}^{N}\mu_{i}^{2}d\varphi_{i}^{2}\right).$

We can adopt a similar form of $\mu_{i}$ from even $d$ case with additional
$\mu_{0}$ term that gives us the following result

$\displaystyle\mu_{0}$ $\displaystyle=$ $\displaystyle\cos\theta_{0}~{},$
$\displaystyle\mu_{1}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\cos\theta_{1}~{},$ $\displaystyle\mu_{2}$
$\displaystyle=$ $\displaystyle\sin\theta_{0}\sin\theta_{1}\cos\theta_{2}~{},$
$\displaystyle\vdots$ $\displaystyle\mu_{N-1}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\sin\theta_{1}\dots\cos\theta_{N-1}~{},$
$\displaystyle\mu_{N}$ $\displaystyle=$
$\displaystyle\sin\theta_{0}\sin\theta_{1}\dots\sin\theta_{N-1}~{}.$ (89)

Consider, $\mu_{0}=\frac{\rho_{0}}{r}$ and $-r\leq\rho_{0}\leq r$. This leads
to $-1\leq\mu_{0}\leq 1$ which means that $0\leq\theta_{0}\leq\pi$. The domain
of the rest of the $\theta_{i}$s have been specified in previous subsection
and that completes our parametrization of $S^{d-1}$ in polyspherical
coordinate for odd $d$ case.

The inverse transformation to standard Euclidean coordinate for this case is
similar with the one we found for even $d$ case, namely

$\displaystyle x^{i}$ $\displaystyle=$ $\displaystyle
r\mu_{i}\cos\varphi_{i},$ (90) $\displaystyle x^{2i}$ $\displaystyle=$
$\displaystyle r\mu_{i}\sin\varphi_{i},$ (91) $\displaystyle x^{d}$
$\displaystyle=$ $\displaystyle r\mu_{0}.$ (92)

where $i$ runs from $1$ to $N$. The only modification we need to accommodate
$\mu_{0}$ is the last equation (92) and we can choose to place at the $d$-th
coordinate, $x^{d}$.

##  Appendix B Bogomolnyi equation from Equation of Motions

In this part of the paper we are going to show explicitly how we can reproduce
the same differential equation of field $\xi(r)$ from both BPS Lagrangian
method and dynamical field equation (14)

###  B.1 Field Equation of The First BPS Submodel

Firstly, let us focus on the field equation of the first submodel which is
recast in this following form

$\displaystyle\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\left[\frac{2d\gamma_{d}}{d!}\nabla_{j}\left(\phi^{b,i}H_{i}^{j}(d,d)\right)-\frac{\partial}{\partial\phi^{b}}\gamma_{0}V\right]=0~{}.$
(93)

From the $O(d+1$ model constraint, we have $\phi^{a}\phi^{a}=1$ which implies
that

$\displaystyle\phi^{a}\phi^{a}_{i}=0~{}.$ (94)

From this equation, we can prove another useful identity from the derivative
of $\phi^{a}\phi^{a}_{i}$, given by
$\nabla_{j}(\phi^{a}\phi^{a}_{i})=\phi^{a}_{i}\phi^{a}_{j}+\phi^{a}\nabla_{j}\phi^{a}_{i}$,
which leads to

$\phi^{a}\nabla_{j}\phi^{a}_{i}=-\phi^{a}_{i}\phi^{a}_{j}=-D_{ij}~{}.$ (95)

Applying these two identities to (93) gives

$\displaystyle\frac{2\gamma_{d}}{(d-1)!}\left[\nabla_{j}\left(\phi^{c,i}H_{i}^{j}(d,d)\right)+\phi^{c}H_{i}^{j}(d,d)D^{i}_{j}\right]-\gamma_{0}\left[\frac{\partial}{\partial\phi^{c}}V-\phi^{c}\phi^{b}\frac{\partial}{\partial\phi^{b}}V\right]=0~{}.$
(96)

Recall that we choose to work specifically on potential which takes the form
$V\equiv V(\xi)$ and $\phi^{a}\phi^{a}=1$ where the only field which does not
contain angular factor is $\phi^{0}=\cos\xi$, then the summation
$\phi^{b}\frac{\partial}{\partial\phi^{b}}V$ is equal to
$\phi^{0}\frac{\partial}{\partial\phi^{0}}V(\xi)=-V^{\prime}(\xi)\cot\xi$
because the other term for $b\neq 0$ vanishes. Upon substitution of the ansatz
to equations (93) we arrive at a single equation for $\xi$, namely

$\frac{2\gamma_{d}}{(d-1)!}\left[-\frac{1}{r^{d-1}}\frac{d}{dr}\left(r^{d-1}\sin\xi\xi^{\prime}\left(B^{2}\frac{\sin^{2(d-1)}\xi}{r^{2(d-1)}}\right)\right)+d!\cos\xi
B^{2}\frac{\sin^{2(d-1)}\xi}{r^{2(d-1)}}\xi^{\prime}{}^{2}\right]+\gamma_{0}V^{\prime}(\xi)\sin\xi=0$
(97)

and all other equations only differ by some angular factors. Equation (97) is
equivalent to Euler-Lagrange equations derived from effective Lagrangian. With
some manipulations on field equation (97) we can simplify it to become

$B\frac{\sin^{d-1}\xi}{r^{d-1}}\frac{d}{dr}\left(B\frac{\sin^{d-1}\xi}{r^{d-1}}\xi^{\prime}\right)=\frac{\gamma_{0}}{2\gamma_{d}}V^{\prime}(\xi)~{}.$
(98)

Multiplying both sides with $\xi^{\prime}$ gives

$\frac{1}{2}\frac{d}{dr}\left(B\frac{\sin^{d-1}\xi}{r^{d-1}}\xi^{\prime}\right)^{2}=\frac{\gamma_{0}}{2\gamma_{d}}\frac{d}{dr}V(\xi)~{},$
(99)

which implies a Bogomolnyi equation of the form

$B\frac{\sin^{d-1}\xi}{r^{d-1}}\xi^{\prime}\mp\sqrt{\frac{\gamma_{0}}{\gamma_{d}}}\sqrt{V(\xi)}=0~{}.$
(100)

As a conclusion, it is proven that we can rederive the Bogomolnyi equation of
the first submodel from its full field equation, implying a consistent
formulation of BPS Lagrangian method.

###  B.2 Field Equation of The Second BPS Submodel

Again, let us consider the field equation for this submodel, given by

$\displaystyle\left(\delta^{cb}-\phi^{c}\phi^{b}\right)\frac{d\gamma_{d/2}}{\left((d/2)!\right)^{2}}\nabla_{j}\left(\phi^{b,i}H_{i}^{j}(d/2,d)\right)=0~{}.$
(101)

By exploiting the identity given in (94) we can simplify the equation above to
become

$\nabla_{j}\left(\phi^{c,i}H_{i}^{j}(d/2,d)\right)-\phi^{c}H_{i}^{j}(d/2,d)\nabla_{j}\phi^{b,i}=0~{}.$
(102)

Now consider only the second term of the equation above. From equation (95) we
can recast it into the following form

$-H_{i}^{j}(d/2,d)\nabla_{j}\phi^{b,i}=H_{i}^{j}(d/2,d)D^{i}_{j}=\left((d/2)!\right)^{2}\left({}_{d}K_{(d/2)-1}\xi^{\prime
2}+{}_{d}K_{d/2}\frac{\sin^{2}\xi}{r^{2}}\right)\frac{\sin^{d-2}\xi}{r^{d-2}}~{}.$
(103)

If we substitute the ansatz into (102), then we arrive at $d+1$ equations
which are equivalent and differs only by some angular factor. Thus, the
equations become, effectively, single equation of $\xi$, given by

$\displaystyle-\frac{1}{r^{d-1}}\frac{d}{dr}\left(r^{d-1}\sin\xi\xi^{\prime}\left((d/2)!{}_{d}K_{(d/2)-1}\frac{\sin^{d-2}\xi}{r^{d-2}}\right)\right)$
$\displaystyle+\cos\xi\left((d/2)!\right)^{2}\left({}_{d}K_{(d/2)-1}\xi^{\prime
2}+{}_{d}K_{d/2}\frac{\sin^{2}\xi}{r^{2}}\right)\frac{\sin^{d-2}\xi}{r^{d-2}}=0~{}.$
(104)

Equation (B.2) is equivalent to the Euler-Lagrange equation derived from the
effective Lagrangian. Now, Suppose that the corresponding Bogomolnyi equation
of this second submodel is given by

$\xi^{\prime}=-\alpha\frac{\sin\xi}{r}~{},$ (105)

then, by substituting this expression into (B.2) we arrive at

$\frac{2^{d}d}{r}\left(-1+\frac{r_{0}^{2\alpha}}{r^{2\alpha}}\right)\left(\frac{1}{\frac{r^{a}}{r_{0}^{a}}+\frac{r_{0}^{a}}{r^{a}}}\right)\left(\alpha^{2}{}_{d}K_{(d/2)-1}-{}_{d}K_{d/2}\right)=0~{}.$
(106)

Since the equation above must be satisfied for all $r$ then we need to take
$\left(\alpha^{2}{}_{d}K_{(d/2)-1}-{}_{d}K_{d/2}\right)=0$ which lead to
$\alpha=\sqrt{\frac{{}_{d}K_{d/2}}{{}_{d}K_{(d/2)-1}}}$. We can observe that
we arrive at the same expression for $\alpha$ from BPS Lagrangian method, thus
we can conclude that the ansatz is consistent for both Bogmolnyi equation and
full field equation from Euler-Lagrange equation.

##  Appendix C Components of Energy-Momentum Tensor via Eigenvalues of Strain
Tensor

The strain tensor given in (20) possesses a set of eigenvalues of given by
(27) for even $d$ and (26) for odd $d$. These simplified version of
eigenvalues actually came from

$\lambda_{1}^{2}=g^{11}\xi^{\prime
2},~{}~{}~{}\text{and}~{}\lambda_{k}^{2}=g^{kk}\sin^{2}\xi\partial_{k}\vec{n}\cdot\partial_{k}\vec{n}~{}\text{for}~{}k=2,\dots,d~{}.$
(107)

with no summation on repeating indices. In equation (107) we can see the
reason why choosing a suitable ansatz for a particular coordinate is important
to reduce the complexity of the resulting effective model. The main rule for
the choice of ansatz which we have used above is that we need the ansatz to
satisfy proportionality $\partial_{k}\vec{n}\cdot\partial_{k}\vec{n}\propto
g_{kk}$, such that the quantity
$g^{kk}\partial_{k}\vec{n}\cdot\partial_{k}\vec{n}$ does not depend on any
angular coordinates.

We can exploit the eigenvalues given in (107) to calculate the energy-momentum
tensor more simply. Firstly, recall that the energy-momentum tensor is given
by

$T_{ij}=-2\frac{\partial\mathcal{L}}{\partial g^{ij}}+g_{ij}\mathcal{L},$
(108)

and Skyrme model Lagrangian expressed in terms of eigenvalues is given by

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle-\left[\gamma_{0}V+\sum_{i=1}^{d}\gamma_{i}e_{i}(\lambda_{1}^{2},\dots,\lambda_{d}^{2})\right]$
(109) $\displaystyle=$
$\displaystyle-\left[\gamma_{0}V+\sum_{i=1}^{d}\gamma_{i}e_{i}(g^{11}\xi^{\prime
2},g^{22}\sin^{2}\xi,\dots,n_{N}^{2}g^{dd}\sin^{2}\xi\sin^{2}\theta_{0}\dots\sin^{2}\theta_{N-1})\right].$

Because all of the arguments of symmetric polynomials $e_{i}$ contains
$g^{kk}$, then the term $\frac{\partial\mathcal{L}}{\partial g^{ij}}$
implicitly removes all the term which are not proportional to $g^{ij}$. It is
easy to prove that all the $T_{ij}$ for $i\neq j$ vanishes, hence there only
$T_{kk}$ components are left. This $T_{kk}$ components can be easily
calculated as the full effective Lagrangian subtracted by two times the terms
which are proportional to $g^{kk}$ which means that the eigenvalues which are
proportional to $g^{kk}$, $\lambda_{k}^{2}$, changes their sign to $(-)$. In
other words, we have the following form of non-vanishing components of energy-
momentum tensor, given by

$T_{kk}=-g_{kk}\left[\gamma_{0}V+\sum_{i=1}^{d}\gamma_{i}e_{i}(\lambda_{1}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})\right]~{}.$
(110)

It is important to note that the expression of the energy-momentum tensor in
(110) might not be the most general one, since we restrict our case to the
ansatz satisfying the properties mentioned above to suit the spacetime metric.

Let us focus on the case of BPS submodels. The spatial diagonal component of
the energy-momentum tensor for the first submodel is given by

$T^{(1)}_{kk}=-g_{kk}\left[\gamma_{0}V+\gamma_{d}e_{d}(\lambda_{1}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})\right]=-g_{kk}\left[\gamma_{0}V-\gamma_{d}\text{det}D\right]~{}.$
(111)

Substitution of the Bogomolnyi equation,
$\sqrt{\gamma_{0}V}\pm\sqrt{\gamma_{d}\text{det}D}=0$ to $T^{(1)}_{kk}$ above
gives

$T^{(1)}_{kk}=0~{}\text{for}~{}k=1,\dots,d~{}.$ (112)

Thus, for every $B$ the spatial components of energy-momentum tensor of the
first submodel are equal to zero. For the second submodel, the spatial
diagonal components of energy-momentum tensor is given by

$T^{(2)}_{kk}=-g_{kk}\gamma_{d/2}~{}e_{d/2}(\lambda_{1}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})~{}.$
(113)

If we impose the self-duality condition for the second submodel where all
eigenvalues are identical, then the equation becomes
$T^{(2)}_{kk}=-g_{kk}\gamma_{d/2}\lambda_{k}^{2}\left(C^{d-1}_{d/2}-C^{d-1}_{(d/2)-1}\right)=0$,
which implies that, if BPS limit is saturated then we always have vanishing
spatial components of the energy-momentum tensor. This is not the case for
non-self-dual Skyrmion with $B>1$, implying an unstable Skyrmion solution. We
demonstrate this property as follows: Suppose that we take the first
eigenvalues to be $\xi^{\prime}{}^{2}$ and factorize it from $e_{d/2}$ or
$e_{N}$,888we take $\frac{d}{2}=N$ for this special case from now on to
simplify notations. then we arrive at the following relation

$\displaystyle T^{(2)}_{kk}$ $\displaystyle=$ $\displaystyle-
g_{kk}\gamma_{N}e_{N}(\lambda_{1}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})$
(114) $\displaystyle=$ $\displaystyle-
g_{kk}\gamma_{N}\left[\lambda_{1}^{2}e_{N-1}(\lambda_{2}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})+e_{N}(\lambda_{2}^{2},\dots,-\lambda_{k}^{2},\dots,\lambda_{d}^{2})\right]~{}.$

We know that the rest of the eigenvalues are proportional to
$\frac{\sin^{2}\xi}{r^{2}}$, hence it is useful to factorize it and then
divide the case of energy-momentum tensor into polar components corresponding
to $\theta_{k}$ and azimuthal components corresponding to $\varphi_{k}$. Since
elementary symmetric polynomials satisfies the following identity

$e_{n}(x_{1},\dots,x_{m})=\sum_{i=0}^{n}e_{n-i}(x_{1},\dots,x_{k})e_{i}(x_{k+1},\dots,x_{m})~{},$
(115)

then we can divide the eigenvalues related to $N-1$ polar coordinates which
are proportional to 1 and the eigenvalues related to $N$ azimuthal coordinates
which are proportional to their corresponding winding numbers $n_{k}$ and
substitute it into (115) after factorizing the $\frac{\sin^{2}\xi}{r^{2}}$
factor. This process lead us to the following form of polar components

$\displaystyle T^{(2)}_{22}$ $\displaystyle=$
$\displaystyle-r^{2}\xi^{\prime}{}^{2}\frac{\sin^{d-2}\xi}{r^{d-2}}\sum_{i=0}^{N-1}e_{N-1-i}(\underbrace{1,\dots,-1,\dots,1}_{N-1})e_{i}(n_{1}^{2},\dots,n_{N}^{2})$
(116)
$\displaystyle-r^{2}\frac{\sin^{d}\xi}{r^{d}}\sum_{i=0}^{N}e_{N-i}(\underbrace{1,\dots,-1,\dots,1}_{N-1})e_{i}(n_{1}^{2},\dots,n_{N}^{2})$
$\displaystyle=$
$\displaystyle-r^{2}\xi^{\prime}{}^{2}\frac{\sin^{d-2}\xi}{r^{d-2}}\sum_{i=0}^{N-1}\left(C^{N-2}_{N-1-i}-C^{N-2}_{N-2-i}\right)e_{i}(n_{1}^{2},\dots,n_{N}^{2})$
$\displaystyle-r^{2}\frac{\sin^{d}\xi}{r^{d}}\sum_{i=0}^{N}\left(C^{N-2}_{N-i}-C^{N-2}_{N-1-i}\right)e_{i}(n_{1}^{2},\dots,n_{N}^{2})~{},$
$\displaystyle g^{22}T^{(2)}_{22}$ $\displaystyle=$ $\displaystyle
g^{33}T^{(2)}_{33}=g^{44}T^{(2)}_{44}=\dots=g^{NN}T^{(2)}_{NN}~{},$ (117)

and for azimuthal components we have

$\displaystyle T^{(2)}_{(N+k)(N+k)}$ $\displaystyle=$ $\displaystyle-
g_{(N+k)(N+k)}\xi^{\prime}{}^{2}\frac{\sin^{d-2}\xi}{r^{d-2}}\sum_{i=0}^{N-1}e_{N-1-i}(\underbrace{1,\dots,1}_{N-1})e_{i}(n_{1}^{2},\dots,-n_{k}^{2},\dots,n_{N}^{2})$
(118) $\displaystyle-
g_{(N+k)(N+k)}\frac{\sin^{d}\xi}{r^{d}}\sum_{i=0}^{N}e_{N-i}(\underbrace{1,\dots,1}_{N-1})e_{i}(n_{1}^{2},\dots,-n_{k}^{2},\dots,n_{N}^{2})$
$\displaystyle=$ $\displaystyle-
g_{(N+k)(N+k)}\xi^{\prime}{}^{2}\frac{\sin^{d-2}\xi}{r^{d-2}}\sum_{i=0}^{N-1}C^{N-1}_{N-1-i}e_{i}(n_{1}^{2},\dots,-n_{k}^{2},\dots,n_{N}^{2})$
$\displaystyle-
g_{(N+k)(N+k)}\frac{\sin^{d}\xi}{r^{d}}\sum_{i=0}^{N}C^{N-1}_{N-i}e_{i}(n_{1}^{2},\dots,-n_{k}^{2},\dots,n_{N}^{2})~{}.$

We can observe that all these spatial components of $B>1$ cases are non-zero,
but unfortunately equation (116) stops to work for $d=2$ where $N=1$. This
peculiarity leads to vanishing spatial components as we show below.

Let us consider the case of $d=2$ which has two eigenvalues,

$\lambda_{1}^{2}=\xi^{\prime
2},~{}~{}~{}\lambda_{2}^{2}=n_{1}^{2}\frac{\sin^{2}\xi}{r^{2}}~{}.$ (119)

The $T_{11}$ and $T_{22}$ components are given by

$\displaystyle T_{11}$ $\displaystyle=$
$\displaystyle-\left[\gamma_{0}V+\gamma_{1}\left(-\xi^{\prime
2}+n_{1}^{2}\frac{\sin^{2}\xi}{r^{2}}\right)-\gamma_{2}\left(n_{1}^{2}\xi^{\prime
2}\frac{\sin^{2}\xi}{r^{2}}\right)\right]~{},$ (120) $\displaystyle T_{22}$
$\displaystyle=$
$\displaystyle-r^{2}\left[\gamma_{0}V+\gamma_{1}\left(\xi^{\prime
2}-n_{1}^{2}\frac{\sin^{2}\xi}{r^{2}}\right)-\gamma_{2}\left(n_{1}^{2}\xi^{\prime
2}\frac{\sin^{2}\xi}{r^{2}}\right)\right]~{}.$ (121)

We can directly see if we consider the first BPS submodel ($\gamma_{1}=0$) and
impose the BPS limit then both $T_{11}$ and $T_{22}$ vanish. This property can
be found in all first BPS submodels for arbitrary $d$. In contrast with the
first submodel, we can observe that the second BPS submodel
($\gamma_{0}=\gamma_{2}=0$) does have vanishing angular components of energy-
momentum tensor for the case of $d=2$ but we have shown in (116) and (118)
that they are non-zero for $d>2$. This property leads to a stable self-dual
BPS Skyrmion within this ansatz for the second BPS submodel which generally
cannot be found for $d$ other than 2.

###### Acknowledgements.

The work in this paper is supported by GTA Research Group ITB, Riset ITB 2022,
and PDUPT Kemenristek-ITB 2022. A. N. A. would like to acknowledge the support
from the ICTP through the Associates Programme (2018-2023). B. E. G. would
also acknowledge the support from the ICTP through the Associates Programme
(2017-2022). E. S. F. also would like to acknowledge the support from BRIN
through the Research Assistant Programme 2022.

## References

  * Skyrme [1961] T. Skyrme, A Nonlinear field theory, Proc. Roy. Soc. Lond. A A260, 127 (1961).
  * Skyrme [1962] T. H. R. Skyrme, A Unified Field Theory of Mesons and Baryons, Nucl, Phys. 31, 556 (1962).
  * Manton and Sutcliffe [2007] N. S. Manton and P. Sutcliffe,  _Topological solitons_ (Cambridge University Press, 2007).
  * Derrick [1964] G. Derrick, Comments on nonlinear wave equations as models for elementary particles, J. Math. Phys. 5, 1252 (1964).
  * Dothan and Biedenharn [1987] Y. Dothan and L. Biedenharn, Old Models Never Die: The Revival of the Skyrme Model, Comments Nucl. Part. Phys. 17, 63 (1987).
  * ’t Hooft [1974] G. ’t Hooft, A Planar Diagram Theory for Strong Interactions, Nucl. Phys. B 72, 461 (1974).
  * Witten [1983] E. Witten, Current Algebra, Baryons, and Quark Confinement, Nucl. Phys. B 223, 433 (1983).
  * Donati  _et al._ [2016] S. Donati, L. Dominici, G. Dagvadorj, D. Ballarini, M. De Giorgi, A. Bramati, G. Gigli, Y. G. Rubo, M. H. Szymańska, and D. Sanvitto, Twist of generalized skyrmions and spin vortices in a polariton superfluid, Proceedings of the National Academy of Sciences 113, 14926 (2016).
  * Khawaja and Stoof [2001] U. Khawaja and H. Stoof, Skyrmion physics in Bose-Einstein ferromagnets, Phys. Rev. A 64, 043612 (2001), arXiv:cond-mat/0107302 .
  * Baskaran [2011] G. Baskaran, Possibility of Skyrmion Superconductivity in Doped Antiferromagnet K2Fe4Se5, (2011), arXiv:1108.3562 [cond-mat.supr-con] .
  * Kiselev  _et al._ [2011] N. S. Kiselev, A. N. Bogdanov, R. Schäfer, and U. K. Rößler, Chiral skyrmions in thin magnetic films: new objects for magnetic storage technologies?, Journal of Physics D: Applied Physics 44, 392001 (2011).
  * Fukuda and Žumer [2011] J. Fukuda and S. Žumer, Quasi-two-dimensional Skyrmion lattices in a chiral nematic liquid crystal, Nature Communications 2 (2011).
  * Tsesses  _et al._ [2018] S. Tsesses, E. Ostrovsky, K. Cohen, B. Gjonaj, N. H. Lindner, and G. Bartal, Optical skyrmion lattice in evanescent electromagnetic fields, Science 361, 993 (2018), https://science.sciencemag.org/content/361/6406/993.full.pdf .
  * Gudnason and Nitta [2015] S. B. Gudnason and M. Nitta, D-brane solitons in various dimensions, Phys. Rev. D 91, 045018 (2015), arXiv:1412.6995 [hep-th] .
  * Prasetyo and Ramadhan [2017] I. Prasetyo and H. S. Ramadhan, Classical defects in higher-dimensional Einstein gravity coupled to nonlinear $\sigma$-models, Gen. Rel. Grav. 49, 115 (2017), arXiv:1707.06415 [gr-qc] .
  * Cartwright  _et al._ [2021] C. Cartwright, B. Harms, and M. Kaminski, Topological or rotational non-Abelian gauge fields from Einstein-Skyrme holography, JHEP 03, 229, arXiv:2010.03578 [hep-th] .
  * Blanco-Pillado  _et al._ [2009] J. J. Blanco-Pillado, H. S. Ramadhan, and N. Shiiki, Skyrme Branes, Phys. Rev. D 79, 085004 (2009), arXiv:0809.0930 [hep-th] .
  * Gudnason and Nitta [2017] S. B. Gudnason and M. Nitta, A higher-order Skyrme model, JHEP 2017 (9).
  * Adam  _et al._ [2010a] C. Adam, J. Sanchez-Guillen, and A. Wereszczynski, A Skyrme-type proposal for baryonic matter, Phys. Lett. B 691, 105 (2010a), arXiv:1001.4544 [hep-th] .
  * Adam  _et al._ [2010b] C. Adam, J. Sanchez-Guillen, and A. Wereszczynski, A BPS Skyrme model: Mathematical properties and physical applications, Acta Phys. Polon. B 41, 2717 (2010b).
  * Adam  _et al._ [2020a] C. Adam, K. Oles, and A. Wereszczynski, The Dielectric Skyrme model, Phys. Lett. B 807, 135560 (2020a), arXiv:2005.00018 [hep-th] .
  * Adam  _et al._ [2020b] C. Adam, A. García Martín-Caro, M. Huidobro, R. Vázquez, and A. Wereszczynski, A new consistent neutron star equation of state from a generalized Skyrme model, Phys. Lett. B 811, 135928 (2020b), arXiv:2006.07983 [hep-th] .
  * Brihaye  _et al._ [2017] Y. Brihaye, C. Herdeiro, E. Radu, and D. H. Tchrakian, Skyrmions, Skyrme stars and black holes with Skyrme hair in five spacetime dimension, JHEP 11, 037, arXiv:1710.03833 [gr-qc] .
  * Nakamula  _et al._ [2017] A. Nakamula, S. Sasaki, and K. Takesue, Atiyah-Manton Construction of Skyrmions in Eight Dimensions, JHEP 03, 076, arXiv:1612.06957 [hep-th] .
  * Arthur  _et al._ [1996] K. Arthur, G. Roche, D. H. Tchrakian, and Y. S. Yang, Skyrme models with selfdual limits: d = 2,3, J. Math. Phys. 37, 2569 (1996).
  * Manton [1987] N. S. Manton, Geometry of skyrmions, Comm. Math. Phys. 111, 469 (1987).
  * Adam  _et al._ [2016] C. Adam, O. Kichakova, Y. Shnir, and A. Wereszczynski, Hairy black holes in the general Skyrme model, Phys. Rev. D 94, 024060 (2016), arXiv:1605.07625 [hep-th] .
  * Gudnason  _et al._ [2016] S. B. Gudnason, M. Nitta, and N. Sawado, Black hole Skyrmion in a generalized Skyrme model, JHEP 09, 055, arXiv:1605.07954 [hep-th] .
  * Faddeev [1976] L. Faddeev, Some Comments on the Many Dimensional Solitons, Lett. Math. Phys. 1, 289 (1976).
  * Piette  _et al._ [1995a] B. M. A. G. Piette, B. J. Schroers, and W. J. Zakrzewski, Dynamics of baby skyrmions, Nucl. Phys. B 439, 205 (1995a), arXiv:hep-ph/9410256 .
  * Piette  _et al._ [1995b] B. M. A. G. Piette, B. J. Schroers, and W. J. Zakrzewski, Multi - solitons in a two-dimensional Skyrme model, Z. Phys. C 65, 165 (1995b), arXiv:hep-th/9406160 .
  * Fadhilla  _et al._ [2021] E. S. Fadhilla, B. E. Gunara, and A. N. Atmaja, BPS Skyrme submodels of the five-dimensional Skyrme model, JHEP 07, 090, arXiv:2010.03226 [hep-th] .
  * Atmaja  _et al._ [2020] A. N. Atmaja, B. E. Gunara, and I. Prasetyo, BPS submodels of the generalized Skyrme model and how to find them, Nucl. Phys. B 955, 115062 (2020), arXiv:1902.02988 [hep-th] .
  * Stȩpień [2019] L. T. Stȩpień, On Bogomolny equations in the Skyrme model, Acta Phys. Polon. B 50, 65 (2019), arXiv:1804.09020 [physics.gen-ph] .
  * Adam  _et al._ [2009] C. Adam, P. Klimas, J. Sanchez-Guillen, and A. Wereszczynski, Pullback of the Volume Form, Integrable Models in Higher Dimensions and Exotic Textures, J. Math. Phys. 50, 022301 (2009), arXiv:0810.1943 [hep-th] .
  * Adam  _et al._ [2013] C. Adam, C. D. Fosco, J. M. Queiruga, J. Sanchez-Guillen, and A. Wereszczynski, Symmetries and exact solutions of the BPS Skyrme model, J. Phys. A 46, 135401 (2013), arXiv:1210.7839 [hep-th] .
  * Dolan and Tchrakian [1988] B. P. Dolan and D. H. Tchrakian, New Lagrangians for Bosonic M-branes With Vanishing Cosmological Constant, Phys. Lett. B 202, 211 (1988).