Date: ]

# Constrained instanton approximation of Skyrmions with massive pions

Alberto García Martín-Caro<EMAIL_ADDRESS>Departamento de Física
de Partículas, Universidad de Santiago de Compostela and Instituto Galego de
Física de Altas Enerxias (IGFAE) E-15782 Santiago de Compostela, Spain

([)

###### Abstract

We present the idea of using the holonomy along a line of a constrained
instanton solution (an approximate solution of the Euclidean equations of
motion subject to a constraint) of an $\rm SU(2)$ Yang-Mills-Higgs theory to
approximate the Skyrmion solution in a chiral model with massive pions. The
fact that the gauge field acquires a nonzero mass due to the Higgs mechanism
implies that a constrained instanton decays exponentially far from its center,
and so does the Skyrmion configuration that it generates via the Atiyah-Manton
construction. This is precisely the desired behavior at large distances for
the true Skyrmion solutions when the pion mass is included.

##  I Introduction

The vacuum structure of an $\rm SU(2)$ Yang-Mills theory in $3+1$ dimensions
is nontrivial, in the sense that it presents a discrete set of degenerate
classical minima. These minima are associated to static (in the temporal
gauge), pure gauge configurations of the form
$A_{\mu}=\Omega^{-1}(\mathbf{x})\partial_{\mu}\Omega(\mathbf{x})$, with ${\rm
SU(2)}\ni\Omega(\mathbf{x})\rightarrow 1$ as $|\mathbf{x}|\rightarrow\infty$.
Mathematically, these configurations can be seen as continuous mappings
$\mathbb{R}^{3}\cup\\{\infty\\}\sim S^{3}\rightarrow{\rm SU(2)}\sim S^{3}$,
which are classified by their third homotopy group,
$\pi_{3}(S^{3})=\mathbb{Z}$. In other words, vacuum gauge field configurations
are separated in different homotopy classes, labeled by an integer
$n[\Omega]$, which can be obtained as an integral

$n[\Omega]=-\frac{1}{24\pi^{2}}\int_{S^{3}}d^{3}x\epsilon^{ijk}\operatorname{Tr}\left(\Omega^{-1}\partial_{i}\Omega\Omega^{-1}\partial_{j}\Omega\Omega^{-1}\partial_{k}\Omega\right).$
(1)

Vacuum configurations with different topological degree cannot be continuously
deformed into each other without generating non-vacuum gauge fields, so these
vacua are separated by a potential barrier in the quantum theory. The
tunneling transition between classically degenerate vacua of a quantum theory
is generally described by instantons, which are solutions of the Euclidean
field equations, i.e absolute minimizers of the Euclidean action

$S_{\text{E}}=\frac{1}{4g^{2}}\int\Tr(F^{a}_{\mu\nu}F^{a}_{\mu\nu})d^{4}x,$
(2)

with
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+\epsilon^{abc}A^{b}_{\mu}A^{c}_{\nu}$
is the gauge field strength tensor, and
$A_{\mu}=-ig\tfrac{\tau^{a}}{2}A^{a}_{\mu}$ is the Lie algebra-valued Yang-
Mills field.

Instanton solutions are localised in spacetime, and can also be classified by
their homotopy class when seen as mappings on $\mathbb{R}^{4}$ that asymptote
to pure gauge configurations on its boundary, $S^{3}_{\infty}$ (so that their
action remains finite). Therefore, we can also think of Yang-Mills instantons
as topological solitons in the Euclidean version of the theory, with
topological charge (instanton number) given by

$n=-\frac{1}{16\pi^{2}}\int d^{4}x\Tr\,(F_{\mu\nu}\tilde{F}_{\mu\nu}).$ (3)

These tunnelling solutions are usually necessary in order to understand
various non-perturbative processes in the Standard Model and other
(super)Yang-Mills theories. In particular, they play an important role in
Quantum Chromodynamics (QCD), where they have been used to solve the axial
$\rm U(1)$ problem ’t Hooft (1986) or to explain chiral symmetry breaking
Diakonov (1996).

Further, $\rm SU(2)$ Yang-Mills instantons are also related to solitons in the
Skyrme model, a relativistic, $3+1$ dimensional chiral effective field theory
first proposed by Tony Skyrme Skyrme (1961). The Skyrme model describes the
dynamics of self-interacting pion fields, which are collected into an $\rm
SU(2)$-valued field $U(t,\mathbf{x})$. Baryons are described in this model as
topological solitons, called Skyrmions, whose topological charge is identified
with baryon number. In a seminal paper Atiyah and Manton (1989), M. Atiyah and
N. Manton proposed that a (static) Skyrme field configuration $U(\vec{x})$
with baryon number $k$ may be obtained as the holonomy of an $\rm SU(2)$ Yang-
Mills $k$-instanton configuration along the Euclidean time direction, formally

$U(\mathbf{x})=\mathcal{P}\exp(-\int^{\infty}_{-\infty}A_{4}(\mathbf{x},x_{4})dx_{4}),$
(4)

with $\mathcal{P}$ denoting path ordering of the exponential. Such a
construction is almost gauge invariant, since the effect of a gauge
transformation of $A_{\mu}$ by $g(x)$ is just a conjugation of $U$ by the
asymptotic value $g(\infty)$. This can be understood as an isospin rotation of
the corresponding Skyrme field.

Moreover, the moduli space of $\rm SU(2)$ instanton configurations of charge
$k$, including the choice of gauge at infinity, is a connected manifold of
dimension $8k$ Atiyah, Hitchin, and Singer (1977); Jackiw and Rebbi (1977).
Hence, the set of Skyrme fields with $B=k$ obtained via the Atiyah-Manton
construction will be an $8k-1$ dimensional manifold, as a global time
translation of the instanton leaves the Skyrme field unchanged. This manifold
is sometimes called the _moduli space of $B=k$ instanton-generated Skyrmions_,
since even though the Skyrme fields so-constructed are not exact solutions to
the Skyrme field equations, some are good approximations to minimal energy
Skyrmions and other field configurations describing their low energy dynamics
in that topological sector. The Atiyah-Manton approximation was further
developed in Leese and Manton (1994); Hosaka, Oka, and Amado (1991); Sutcliffe
(1994), and it was also extended in order to include instantons at finite
temperature, Cork (2018); Cork and Halcrow (2022) and skyrmions from
gravitational instantons Dunajski (2013). Its application to describe low
energy interactions between nuclei was studied shortly after it was first
proposed Atiyah and Manton (1993); Hosaka, Oka, and Amado (1991), and has been
put forward recently in Halcrow and Winyard (2021); Halcrow and Harland
(2022).

The useful idea that Skyrmions can be seen as instantons in a higher
dimensional space has been thoroughly studied in the holographic context. In
Sutcliffe Sutcliffe (2010), it is explicitly shown that the Skyrme model
appears, together with a Kaluza-Klein tower of vector mesons, after a
dimensional reduction of a $4+1$ dimensional Yang-Mills theory. This
phenomenon is also realized in several holographic QCD models. In particular,
the Sakai-Sugimoto model, a holographic dual of QCD with $N_{f}$ massless
quarks proposed in Sakai and Sugimoto (2005). The low energy effective action
of such a model consists of a five-dimensional Yang-Mills and Chern-Simons
theory on a curved background. In this model, the massless pion and an
infinite tower of massive (axial-)vector mesons are interpreted as Kaluza-
Klein states associated with the holographic direction, and baryons are
identified as D4-branes wrapped on a non-trivial four-cycle in the D4
background. Such a D4-brane is realized as an instanton configuration in the
worldvolume gauge theory. Also, the pion effective action obtained from the KK
reduction of the $4+1d$ Yang-Mills theory is precisely that of the Skyrme
model, in which baryons appear as solitons. It can be shown that the baryon
number of a Skyrmion, i.e the winding number carried by the pion field, is
equivalent to the instanton number in the Yang-Mills theory. In this way, the
Sakai-Sugimoto model presents two dual descriptions of baryons, as wrapped
D-branes (instantons) and as Skyrmions in the KK-reduced theory.

Apart from QCD-like theories, instanton-mediated processes are also important
in more complicated gauge theories presenting the Higgs mechanism, such as the
electroweak sector of the Standard Model, where they are associated to the
violation of baryon number Espinosa (1990). Such transitions become relevant
only at very high energies and have been proposed as a plausible source of
baryon number violation leading to baryogenesis in the early universe Rubakov
and Shaposhnikov (1996). However, the Higgs mechanism poses a slight
complication concerning instantons in these gauge theories. Indeed, there are
no solutions to the Euclidean field equations, i.e. no exact minima of the
Euclidean action in sectors with nonzero topological degree. The reason, as we
will discuss below, is that the action for instanton-like configurations, with
a nontrivial Higgs field, depends on the instanton size and decreases as it
tends to zero. To evaluate the functional integral in that case one introduces
a constraint that fixes the size of the configuration ’t Hooft (1976), then
minimizes the action under this constraint and finally integrates over all
possible values of the instanton size. The constrained solutions that one uses
to calculate instantonic contribution to the path integral in the Higgs phase
are called _constrained instantons_.

##  II The Skyrme model and the Atiyah-Manton construction

The Lagrangian of the Skyrme model is given by

$\mathscr{L}=-\frac{f_{\pi}}{16}\Tr(L_{\mu}L^{\mu})+\frac{1}{32e^{2}}\Tr([L_{\mu},L_{\nu}]^{2})+\frac{1}{8}m_{\pi}^{2}f_{\pi}^{2}\Tr(U-\bm{1}),$
(5)

where $L_{\mu}=U^{\dagger}\partial_{\mu}U$ is the associated left-invariant
Maurer-Cartan form. The lagrangian depends on the three parameters
$f_{\pi},m_{\pi}$ and $e$, namely, the pion decay constant, the pion mass and
the (dimensionless) Skyrme parameter. The first two constants are common to
every chiral effective theory involving pions, and can be fixed by low energy
constraints on their properties. The latter is associated to the fourth-order
Skyrme term, which allows for stable soliton solutions. Such solutions are
called Skyrmions, and correspond to field configurations that minimize the
static energy functional

$E=-\int\Big{[}\frac{1}{2}\Tr(L_{i}L_{i})+\frac{1}{16}\Tr([L_{i},L_{j}]^{2})+m^{2}\Tr(U-\bm{1})\Big{]}d^{3}x,$
(6)

where we have changed our energy and lenght units to $f_{\pi}/4e$ and
$2/ef_{\pi}$, so that the dimensionless pion mass parameter is related to the
pion mass $m_{\pi}$ by $m=2m_{\pi}/ef_{\pi}$.

Furthermore, Skyrmions are topological solitons, since any finite energy
configuration must satisfy the boundary condition
${U(\mathbf{x})\rightarrow\bm{1}}$ as $|\mathbf{x}|\rightarrow\infty$. Thus,
static skyrmions are maps $U:\mathbb{R}^{3}_{\infty}\sim S^{3}\mapsto{\rm
SU(2)}\sim S^{3}$, classified by their third homotopy group,
$\pi_{3}(S^{3})=\mathbb{Z}$, into homotopy classes labelled by their
corresponding topological degree, which can be written in terms of $L_{\mu}$
as

$B=\int B_{0}d^{3}x\in\mathbb{Z},\qquad
B_{\mu}=\frac{1}{24\pi^{2}}\varepsilon_{\mu\nu\rho\sigma}\Tr(L^{\nu}L^{\rho}L^{\sigma}).$
(7)

In this work, we will be primarily interested in the instanton-generated $B=1$
Skyrmion. Indeed, the $B=1$ Skyrmion solution is well known, and can be
written using the so-called hedgehog configuration

$U(\mathbf{x})=e^{if(\rho)\mathbf{n}^{a}\sigma_{a}},\qquad\text{where}\qquad\rho=|\mathbf{x}|,\quad\mathbf{n}=\frac{\mathbf{x}}{\rho},$
(8)

and $f(\rho)$ is the radial profile of the Skyrmion, which, for the $B=1$
case, must satisfy $f(0)=\pi$ and $f(\infty)=0$. Introducing the hedgehog
ansatz into the static energy functional, one obtains

$E=4\pi\int\left(\rho^{2}f^{\prime 2}+2\left(f^{\prime
2}+1\right)\sin^{2}f+\frac{\sin^{4}f}{\rho^{2}}+2m^{2}\rho^{2}(1-\cos
f)\right)d\rho,$ (9)

On the other hand, the (Euclidean, $4d$) pure Yang-Mills action presents
instanton solutions which, in the ’t Hooft gauge, read

$A_{\mu}(x)=\frac{1}{2}\sigma_{\mu\nu}\partial_{\nu}\log\alpha(x),$ (10)

with $\alpha$ a solution of the Laplace equation in $\mathbb{R}^{4}$,
$\partial^{2}\alpha=0$, and $\sigma_{\mu\nu}$ the $\mathfrak{so(4)}$ Lie
algebra generators, defined as 111we are using the conventions in Espinosa
(1990)

$\sigma_{\mu\nu}=\frac{1}{4i}(\sigma_{\mu}\sigma_{\nu}^{\dagger}-\sigma_{\mu}^{\dagger}\sigma_{\nu}),\qquad\text{with}\quad\sigma_{\mu}=(\bm{\tau},i\bm{1}_{2}).$
(11)

In particular, the $k=1$-instanton is given by

$\alpha(x)=1+\frac{\lambda^{2}}{|x-X|^{2}}=1+t,\quad
t\doteq\frac{\lambda^{2}}{|x-X|^{2}},$ (12)

with $\lambda\in\mathbb{R}$, $X\in\mathbb{R}^{4}$ the moduli parametrizing the
size and position of the instanton in Euclidean space, respectively. These
$1$-instantons generate, after calculating the corresponding holonomy as in
(4), a seven-dimensional manifold of B = 1 Skyrme hedgehogs, the coordinates
being the position, orientation and scale size, and whose radial profile is
given by

$f(\rho)=\pi\left[1-\left(1+\frac{\lambda^{2}}{\rho^{2}}\right)^{-1/2}\right].$
(13)

Although instantons are scale invariant, the energy of the associated Skyrmion
configurations depends on the particular value of the (in principle,
arbitrary) parameter $\lambda$. Therefore, its value must be fixed so that the
energy of the resulting Skyrme field is minimized, which reduces the moduli
space by one dimension. It turns out that, for the massless case, the value of
this scale is $\lambda^{2}=2.11$, with an energy of $E^{\rm
inst}_{m=0}=1.243\times 12\pi^{2}$, which is only $1\%$ above that of the true
Skyrmion solution.

##  III Constrained instantons in Yang-Mills-Higgs theories

We now consider an $\rm SU(2)$ Yang-Mills-Higgs theory on four dimensional
Euclidean space $\mathbb{R}^{4}$, given by the Lagrangian

$\mathscr{L}=\frac{1}{g^{2}}\quantity{\frac{1}{4}F_{\mu\nu}^{a}F_{\mu\nu}^{a}+\kappa\quantity[(D_{\mu}\phi)^{\dagger}D_{\mu}\phi+\frac{1}{4}(\phi^{\dagger}\phi-\mu^{2})]},$
(14)

where
$D_{\mu}=\partial_{\mu}-i\frac{\tau^{a}}{2}A^{a}_{\mu}\equiv\partial_{\mu}+A_{\mu}$
is the associated covariant derivative, $\phi$ an $\rm SU(2)$ Higgs doublet
$g$ and $\kappa$ the gauge and Higgs couplings, respectively. Finally,
$\mu^{2}$ is the vacuum expectation value of the Higgs field in the
spontaneously broken phase, in which the gauge field acquires a mass
$m=\sqrt{\tfrac{\kappa\mu^{2}}{2}}$ due to the Higgs mechanism.

For $\mu=0$, the Yang-Mills-Higgs field equations present instanton solutions
which are the same as in the pure Yang-Mills case, with $\phi=0$. On the other
hand, for $\mu\neq 0$, the nontrivial vacuum expectation value of the Higgs
field breaks conformal symmetry, and there cannot exist an exact instanton
solution of the coupled field equations. This can be easily seen by applying
Derrick’s scaling argument Affleck (1981). Indeed, rescaling any instanton
field configuration with $\phi\neq 0$ as

$A_{\mu}(x)\rightarrow aA_{\mu}(ax),\quad\phi(x)\rightarrow\phi(ax),$ (15)

the action becomes $(y=ax)$:

$S=\int
d^{4}y\frac{1}{g^{2}}\quantity{\frac{1}{4}F_{\mu\nu}^{a}(y)F_{\mu\nu}^{a}(y)+\frac{\kappa}{a^{2}}\quantity[(D_{\mu}\phi(y))^{\dagger}D_{\mu}\phi(y)+\frac{1}{4a^{2}}(\phi^{\dagger}\phi-\mu^{2})]},$
(16)

and since all terms in the integrand are positive, the action can always be
made smaller by taking the limit $a\rightarrow\infty$. Since the expression
for $aA_{\mu}(ax)$ is exactly equal to $A_{\mu}(x)$ with $\lambda$ replaced by
$\lambda/a$, we get that the action-minimizing instanton has zero size.
However, the nontrivial topology of the space of field configurations is not
affected by the Higgs mechanism, so that instanton-like field configurations
with nonzero topological degree still exist as approximate solutions of the
field equations. These are the so-called _constrained instantons_ , and have
been used to approximately calculate nonperturbative contributions to the path
integral in the standard model and (super) Yang-Mills theories at the Higgs
phase ’t Hooft (1976); Nielsen (2005); Rubakov and Shaposhnikov (1996).

The constrained instanton idea is based on the existence of a finite number of
destabilizing directions in the gauge field space, along which the action
always decreases, so that one introduces constraints that prevent deformations
in these directions. Then, one just minimizes the action for fields subject to
such constraints. For the spontaneously broken phase of $\rm SU(2)$ Yang-
Mills-Higgs theory, the destabilizing direction is that parametrized by the
rescaling parameter $a$ (note that this would correspond to a zero mode in the
pure Yang-Mills phase), so one needs to impose a constraint that fixes the
instanton size $\lambda$ at some value.222Also, the constraint must allow
finiteness of the action, which is a nontrivial demand for choosing a valid
constraint, see Nielsen and Nielsen (2000). The contribution to the path
integral coming from these configurations is then obtained by integration over
all (fixed) values of $\lambda$. Remarkably, the Higgs mechanism resolves the
infrared divergence problem due to large instantons dominating the path
integral measure on an otherwise scale-symmetric field theory such as pure
Yang-Mills, since these become suppressed on the Higgs phase ’t Hooft (1976);
Espinosa (1990); Shifman (2022).

The Yang-Mills-Higgs equations coming from the Lagrangian (14) are

$\displaystyle
D_{\mu}F_{\mu\nu}^{a}+\frac{\kappa}{2}\quantity(-i\phi^{\dagger}\tau^{a}\stackon[-4.2pt]{\partial}{\,\smash{\stackon[-2.1pt]{\mkern
0.0mu\mathchar 382\relax}{\rotatebox{180.0}{$\mkern 0.0mu\mathchar
382\relax$}}}}_{\mu}\phi-A^{a}_{\mu}\phi^{\dagger}\phi)$ $\displaystyle=0,$
(17) $\displaystyle D^{2}\phi-\phi(\phi^{\dagger}\phi-\mu^{2})$
$\displaystyle=0.$ (18)

One could naively try to find instantonic solutions to the Yang-Mills-Higgs
equations. The finite action condition in the broken phase implies that gauge
fields must become pure gauge configurations asymptotically, and the Higgs
field must tend to a point in the corresponding vacuum manifold,

$A_{\mu}(x)\rightarrow\Omega\partial_{\mu}\Omega^{-1},\quad\phi(x)\rightarrow\Omega\phi_{\rm
vac},\quad x\rightarrow\infty,\quad\Omega(x)\in{\rm SU(2)}$ (19)

with $\phi_{\rm vac}$ a general base point of the Higgs potential vacuum
manifold. Without loss of generality, we may choose $\phi_{\rm
vac}=(0,\mu)^{\intercal}$. Therefore, a general ansatz for the (spherically
symmetric) $k=1$ instanton would be of the form

$A_{\mu}=\tilde{\xi}(x)\Omega\partial_{\mu}\Omega^{-1},\qquad\phi(x)=(1-\chi(x))\Omega\phi_{\rm
vac},$ (20)

where $\Omega=i\frac{\sigma_{\mu}x^{\mu}}{\sqrt{x_{\mu}x_{\mu}}}$ is an
element of the first nontrivial homotopy class in $\pi_{3}(\rm SU(2))$, and
$\tilde{\xi},\chi$ are functions of the radial euclidean coordinate $r=|x|$
which satisfy the conditions ${\tilde{\xi},\chi\rightarrow 1,0}$ as
$r\rightarrow\infty$. This ansatz corresponds to the generalization of the
BPST instanton solution Belavin  _et al._ (1975) in the regular gauge to the
Higgs phase, and was used in Klinkhamer (1993) to find a constrained instanton
solution using a specific constraint. In particular, for $\mu=0$, i.e. the
pure Yang-Mills case, the instanton solution is given in this gauge by
$\tilde{\xi}=r^{2}/(r^{2}+\lambda^{2})$.

However, we will find it useful to gauge-transform (20) into a new form, in
which $\phi\rightarrow\phi_{\rm vac}$. This can be achieved with a gauge
transformation by $\Omega^{-1}$, which brings the ansatz into the so-called
singular gauge:

$A_{\mu}(x)=\frac{x^{\nu}\sigma_{\mu\nu}}{x^{2}}\xi(r)=-\sigma_{\mu\nu}\partial_{\nu}a(r),\quad\phi(x)=(1-\chi(x))\phi_{\rm
vac},$ (21)

where $\xi$ is a function related to the original $\tilde{\xi}$, and
${\xi=a^{\prime}/r}$. With this gauge choice, the boundary conditions
${A_{\mu}(x)\xrightarrow{x\rightarrow\infty}0,\,\phi(x)\xrightarrow{x\rightarrow\infty}(0,\mu)^{\intercal},}$
are satisfied, and eqs. 17 and 18 linearize for $A_{\mu}(x)$ and $\chi(x)$ in
the $x\rightarrow\infty$ limit,

$\displaystyle\quantity[-\delta_{\mu\nu}\partial^{2}+\partial_{\mu}\partial_{\nu}+m^{2}]A_{\mu}$
$\displaystyle=0,\qquad m=\sqrt{\frac{\kappa\mu^{2}}{2}}$ (22)
$\displaystyle(-\partial^{2}+\mu^{2})\chi$ $\displaystyle=0.$ (23)

Any solution of the linearized equations that satisfies the correct boundary
conditions is given by

$\chi\propto G(x;\mu),\qquad
A_{\mu}\propto\sigma_{\mu\nu}\partial_{\nu}G(x;m),$ (24)

with

$G(x;\alpha)=\frac{\alpha}{r}K_{1}(\alpha r)$ (25)

the Green function associated to the four dimensional, Euclidean Klein-Gordon
operator, $-\partial^{2}+\alpha^{2}$, and $K_{n}(z)$ is the modified Bessel
function of degree $n$.

###  III.1 Perturbative construction

Unfortunately, as previously shown using Derrick’s argument, the boundary
conditions we have imposed are not compatible with a regular solution of the
full system of Yang-Mills-Higgs equations. Nonetheless, it was argued in
Affleck (1981) that a perturbative solution can be constructed if we add extra
terms to the right-hand sides of such equations, with coefficients that can be
adjusted order by order in perturbation theory to obtain the desired boundary
conditions. The presence of these extra terms is equivalent to adding a
constraint at the level of the action, so the obtained solutions are
_constrained instantons_. This procedure was further extended in Nielsen and
Nielsen (2000) with the additional requirement that the perturbative solution
yields a finite value for the action functional, and a constrained instanton
solution in the Higgs phase of $\rm SU(2)$ YMH theory was calculated
perturbatively in the parameter $\mu$ up to $\order{\mu^{4}}$. The solution
looks like the standard self-dual instanton near its core, but decays
exponentially (instead of polynomially) as $\exp(-\mu r)$ far from the center.
We will now reproduce the construction of such solution, following the same
perturbative method.

We start with the singular gauge ansatz (21), and define two real functions
$\alpha=\exp(a)$ and $f=\mu(1-\chi)$ depending on $x$ in terms of the
parameter $t=\tfrac{\lambda^{2}}{r^{2}}$. We also expand these functions as

$\alpha=\sum\limits_{n=0}^{\infty}\alpha_{2n},\qquad
f=\sum\limits_{n=0}^{\infty}f_{2n+1},$ (26)

where the subscript on each function indicates the lowest power of $\mu$
appearing on it. Introducing this ansatz into the Yang-Mills-Higgs Lagrangian,
the equations of motion for $\alpha$ and $f$ are (${}^{\prime}=d/dt$) Nielsen
and Nielsen (2000):

$(\alpha^{-3}t^{3}\alpha^{\prime\prime})^{\prime}=\frac{\kappa\lambda^{2}}{8}f^{2}\alpha^{-3}\alpha^{\prime},\quad\alpha^{2}f^{\prime\prime}-\frac{3}{4}(\alpha^{\prime})^{2}f=\frac{\alpha^{2}\lambda^{2}f}{8t^{3}}(f^{2}-\mu^{2}),$
(27)

Figure 1: Perturbative solution of the constrained instanton profile at
different orders in the Higgs vev parameter, for $\mu^{2}=m^{2}=0.1$. The
large $x$ limit towards it needs to converge is also plotted in green.

which can be solved by quadrature at each order. The required boundary
conditions plus finiteness of the action require that these equations are
modified order by order by adding certain terms at the rhs. These terms (and
the corresponding constraint) are computed at each order in $\mu$, and are
unique up to $\order{\mu^{4}}$ Nielsen and Nielsen (2000). This can be done
iteratively, starting with the zeroth-order instanton solution,
$\alpha_{0}(t)=1+t$, and substituting in the equation for $f$, which yields
$f_{1}=\mu/\sqrt{1+t}$. This solution already satisfies the boundary condition
for $\phi$, so no modification of the equation for $f$ is needed at this
order. In fact, the equation for $f$ does not need to be modified at any order
Nielsen and Nielsen (2000). On the other hand, at second order, the
modification of the equation for $\alpha$ turns out to be unique, (i.e. the
corresponding constraint is unique up to $\order{\mu^{4}}$), and therefore,
one can solve for $\alpha_{2}$ unambiguously,

$\alpha_{2}(t)=\frac{\lambda^{2}m^{2}}{4}\quantity[\log{\frac{\lambda^{2}m^{2}}{4t}}+2\gamma-1],$
(28)

where $\gamma$ is Euler’s constant, appearing in the small $m$ expansion of
the modified Bessel function in (25). In order to find the next order term in
$\alpha$, one needs to solve first for the Higgs field at order three. As we
are not interested in its explicit form, let us skip it and just give the
result for the next order term in $\alpha$, i.e. the $\order{\mu^{4}}$ term.
It’s found in Nielsen and Nielsen (2000) to be

$\displaystyle\alpha_{4}(t)=\frac{\lambda^{4}m^{4}}{64t}\Bigg{\\{}\left(-2t^{2}+4\gamma-5\right)+2\Bigg{[}t\log(t)\Big{(}\frac{\mu^{2}}{m^{2}}(4\gamma+3)+$
$\displaystyle
2\left(\frac{\mu^{2}}{m^{2}}-1\right)\log\left(\frac{\lambda^{2}\mu^{2}}{4}\right)-2\log\left(\frac{m^{2}}{\mu^{2}}\right)-4\gamma+1\Big{)}+$
$\displaystyle+$ $\displaystyle
t(t+1)\log\left(\frac{1}{t}+1\right)\left(4\log\left(\frac{\lambda^{2}m^{2}}{4t}\right)+(t+8\gamma-9)-\frac{\mu^{2}}{m^{2}}(t-5)\right)+$
$\displaystyle+$
$\displaystyle\log\left(\frac{\lambda^{2}m^{2}}{4t}\right)\Bigg{]}-4t\Phi\left(\frac{1}{t}\right)\left(\frac{\mu^{2}}{m^{2}}(1-2t)+4t+1+2\frac{\mu^{2}}{m^{2}}t^{2}\right)\Bigg{\\}},$
(29)

with

$\Phi(x)=\int_{0}^{x}\frac{\log(1+u)}{u}du.$ (30)

For simplicity, in the following we will make the Higgs field vev equal to the
mass (i.e. we fix the Higgs coupling constant $\kappa=2$ ). Doing so, the
expression above simplifies to

$\displaystyle\alpha_{4}(t)$
$\displaystyle=\frac{\lambda^{2}m^{2}}{32t}\Bigg{\\{}\log\left(\frac{\lambda^{2}m^{2}}{4t}\right)-4t(t+1)\Phi\left(\frac{1}{t}\right)+2\gamma-\frac{5}{2}+$
$\displaystyle+$ $\displaystyle
4t\left[(t+1)\log\left(\frac{1}{t}+1\right)\left(\log\left(\frac{\lambda^{2}m^{2}}{4t}\right)+2\gamma-1\right)+\log(t)\right]\Bigg{\\}}.$
(31)

Furthermore, the authors in Nielsen and Nielsen (2000) found a recurrence
relation for obtaining the leading term on each higher order component
$\alpha_{2n}$

$\frac{d}{dt}t^{3}\frac{d^{2}\alpha_{n}}{dt^{2}}\simeq\frac{\lambda^{2}m^{2}}{4}\frac{d\alpha_{n-2}}{dt},$
(32)

whose solutions coincide, at each $n$, with the $n$-th order term in the small
$m$ expansion of $G(x;m)$. The convergence towards $G(x,m)$ of $\alpha(t)$ at
large distances for different orders in $\mu^{2}$ is shown in (1). The
$\alpha_{6}$ and $\alpha_{8}$ are taken as the corresponding solutions to
(32). The perturbative solution is unsatisfactory for our purposes, because it
will always diverge for some value of $r$ at any given truncation order (as we
can see in fig. 1). After computing the corresponding holonomy, this
translates into a Skyrme field that does not satisfy the correct boundary
conditions at infinity. A simple way around this problem is just to consider
the asymptotic solution, namely

$\alpha(x)=1+\frac{\lambda^{2}m}{r}K_{1}(mr),$ (33)

which does of course yield well-behaved Skyrmions at infinity, as we will show
in the next section.

###  III.2 Non-perturbative approach

Apart from the perturbative method, an alternative approach for constructing
constrained instanton configurations in Yang-Mills-Higgs theory was used in
Wang (1994), in which instead of fixing the constraint and trying to solve the
equation for the corresponding constrained instanton, the reverse method is
followed. One chooses a particular functional form for the constrained
instanton, and the corresponding constraint may be systematically calculated
afterwards due to the freedom in choosing the constraint in the first place.
The only requirement in choosing the constrained instanton shape a priori is
that it satisfies some constraint-independent boundary conditions at the
origin and far from its center. Indeed, if one considers again the BPST-like
ansatz, the (dimensionless) functions $\xi,\chi$ can depend on $r=|x|$ only
through two dimensionless combinations, $\lambda/r$ and $mr$. We may expand
these functions both for small $(r\ll\lambda)$ and large $(r\gg m^{-1})$
distances in terms of these two combinations, as Espinosa (1990); Wang (1994):

$\displaystyle\xi(r)$
$\displaystyle=\xi_{0}(\tfrac{\lambda}{r})+(mr)^{2}\xi_{1}(\tfrac{\lambda}{r})+\cdots=$
$\displaystyle=(\tfrac{\lambda}{r})^{2}\xi^{0}(mr)+(\tfrac{\lambda}{r})^{4}\xi^{1}(mr)+\cdots,$
(34) $\displaystyle\chi(r)$
$\displaystyle=\chi_{0}(\tfrac{\lambda}{r})+(mr)^{2}\ln(mr)\chi_{1}(\tfrac{\lambda}{r})+\cdots=$
$\displaystyle=(\tfrac{\lambda}{r})^{2}\chi^{0}(mr)+(\tfrac{\lambda}{r})^{4}\chi^{1}(mr)+\cdots.$
(35)

Thus for recovering the pure Yang-Mills solution in the $m\rightarrow 0$ limit
we require

$\xi_{0}(\tfrac{\lambda}{r})=\frac{2\lambda^{2}}{r^{2}+\lambda^{2}},\qquad\chi_{0}(\tfrac{\lambda}{r})=1-\sqrt{\frac{r^{2}}{r^{2}+\lambda^{2}}},$
(36)

where $\chi_{0}$ above corresponds to the scalar field profile in the
instanton background without backreaction ’t Hooft (1976).

On the other hand, the linearized equations at large distances impose

$\xi^{0}(mr)=m^{2}r^{2}K_{2}(mr),\qquad\chi^{0}(mr)=\tfrac{1}{2}mrK_{1}(mr).$
(37)

A pair of functions proposed by Wang in Wang (1994) satisfying eqs. 34, 35, 36
and 37 are

$\xi(r)=\frac{\lambda^{2}m^{2}K_{2}\left(mr\right)}{[1+\lambda^{2}m^{2}K_{2}\left(mr)/2\right]},\quad\chi(r)=\quantity[1-\left(\frac{r^{2}}{r^{2}+\lambda^{2}mrK_{1}\left(mr\right)}\right)^{\frac{1}{2}}].$
(38)

Obviously, the expressions in eq. 38 are associated with a particular
constraint at the level of the action, which could be calculated, at least
formally. Of course, they are not unique, and a different pair of functions
satisfying the correct boundary conditions could have been chosen.

##  IV Skyrmions from constrained instanton configurations

Despite its success in describing the moduli space of massless pion Skyrmions,
the instanton approximation does not work in the massive pion case. This can
be understood by noting that the asymptotic decay of the instanton generated
profile function is $f\sim\lambda^{2}/(2\rho^{2})$, which has the correct form
for massless pions, but decays too slowly for the massive pion case,
ultimately generating configurations with infinite energy. Thus this form of
the instanton holonomy method is not applicable to massive pions.

However, as we have seen in the previous section, we are still able to find
instanton-like configurations in Yang-Mills-Higgs theory that preserve their
topological properties and nevertheless present an exponential radial decay
due to the Higgs mechanism. Hence, constrained instantons are natural
candidates for generating exponentially decaying Skyrme hedgehogs via the
Atiyah-Manton construction. Indeed, from the ansatz in the previous section we
get

$A_{4}(\mathbf{x},x_{4})=-in^{a}\tau_{a}\quantity(\rho\frac{\lambda^{2}}{(x_{4}^{2}+\rho^{2})^{2}}\frac{\alpha^{\prime}}{\alpha}),$
(39)

with $\alpha(t)$ given by (33). The profile function generated by the
constrained instanton approximation can therefore be written

$f(\rho)=\int_{0}^{\infty}\frac{\rho\lambda^{2}}{(x_{4}^{2}+\rho^{2})^{2}}\frac{\alpha^{\prime}(t)}{\alpha(t)}dx_{4}.$
(40)

On the other hand, if we choose to use the nonperturbative configuration
generated by the functions (38) for constructing the constrained instanton,
the profile function of the corresponding Skyrmion can be directly obtained as

$f(\rho)=\int_{0}^{\infty}\frac{\rho\lambda^{2}m^{2}K_{2}(m\sqrt{x_{4}^{2}+\rho^{2}})}{2(x_{4}^{2}+\rho^{2})[1+\lambda^{2}m^{2}K_{2}(m\sqrt{x_{4}^{2}+\rho^{2}})/2]}dx_{4}.$
(41)

As pointed out in Espinosa (1990), the Lagrangian eq. 14 is invariant under
independent global rotations of $A_{\mu}$ and $\phi$, even after substituting
the ansatz (21). Thus, in the constrained instanton the gauge and Higgs field
orientations are not correlated. Furthermore, the Higgs field orientation is
completely fixed once we choose a specific point in the vacuum manifold,
$\phi$ as the vacuum state in the spontaneously broken phase. On the other
hand, the orientation of $A_{\mu}$ remains a symmetry of the Yang-Mills-Higgs
action and therefore has to be treated as a collective coordinate along with
the position of the constrained instanton. Of course, these collective
coordinates will be inherited by the Skyrmion field configurations generated
by constrained instanton configurations, so that the structure of the
instanton-generated Skyrmion moduli space is preserved even after the addition
of the Higgs field.

Obviously, the constrained instanton solution has lost its scale invariance
due to a finite Higgs vev, so the scale parameter $\lambda$ doesn’t represent
a zero mode anymore, but a quasi-zero mode, since it parametrizes a family of
local minima of the action when the constraint is imposed. Therefore, we still
have the freedom of choosing the value that better fits the true Skyrmion
energy. Also, even for a given, specific constraint, its value will not be
unique, but will depend on the value of the other free parameter, namely, the
gauge field mass $m$, which coincides with the (dimensionless) pion mass
parameter of the Skyrme field. Therefore, there’s a one parameter family of
Skyrmion solutions, each with different values of $m$, and the best value of
$\lambda$ for each one can be understood as a function $\lambda(m)$.

At the end of the day, constrained instantons also generate a seven
dimensional Skyrmion moduli space, as fixing the Higgs vev at infinity does
not break translation invariance, nor global gauge symmetry of the gauge field
solution.

Figure 2: Solution of the Skyrme field profile for $m\in\\{0,0.5,3\\}$, and
their corresponding best approximation using instantons (massless) and
constrained instantons (massive)

As we can see from fig. 2, the similarity between the true solution in the
massive pion case and the constrained instanton approximation is astounding,
and even better than the regular instanton approximation in the massless case.
Indeed, the difference in energies between both cases is smaller than $1\%$,
for instance, $E_{m=0.5}=1.433\times 12\pi^{2}$ and $E^{\rm
c.inst}_{m=0.5}=1.438\times 12\pi^{2}$. In this case, we have
$\lambda^{2}(m=0.5)=1.47$. In table 1 we show the energies for each case in
units of $12\pi$.

$m$ | 0 | 0.5 | 3
---|---|---|---
$E_{\rm real}$ | $1.232$ | $1.433$ | $2.170$
$E^{\rm pert}_{\rm inst}$ | $1.243$ | $1.438$ | $2.250$
$E^{\rm n.pert}_{\rm inst}$ | $1.243$ | $1.435$ | $2.224$
Table 1: Energies of $B=1$ hedgehog obtained by numerically solving the Skyrme
field equations ($E_{\rm real}$) and by the constrained instanton
approximation, both using the perturbative configuration ($E^{\rm pert}_{\rm
inst}$) and Wang’s proposal ($E^{\rm n.pert}_{\rm inst}$).

Now an interesting question is whether we can say anything about higher baryon
number skyrmions. In the pure Yang-Mills case, the self duality of instantons
allows for the construction of the most general $k$ instanton configuration
via the powerful method of ADHM local data Atiyah  _et al._ (1978). In Cork
and Halcrow (2022), this method was used to describe multi-Skyrmion moduli
spaces in terms of the ADHM data of the corresponding instanton approximation.

In our case, however, self duality is broken, and the problem of finding
expressions for $k>1$ constrained instantons, considerably harder. We could
try to construct a $k$ constrained instanton as a superposition of $k$
individual 1-instantons, with arbitrary widths and positions, using a
generalization of the ’t Hooft ansatz to our case, i.e. taking

$\alpha(x)=1+\sum\limits_{i=1}^{k}\frac{\lambda_{i}^{2}m}{|x-X_{i}|}K_{1}(m|x-X_{i}|),$
(42)

where $\lambda_{i}$, $X_{i}$ correspond to the sizes and euclidean positions
of each constrained instanton. Unfortunately, this ansatz does not reproduce
the full picture, as it does not allow for different orientations of the
individual instantons. A modification of the ’t Hooft ansatz that includes the
relative orientation between instantons and does not modify the topological
degree was proposed in Park  _et al._ (2002). It’s based on the introduction
of $3k$ additional parameters, namely, the Euler angles of a $SU(2)$ rotation
matrix $R_{n}$. Then, the fourth component of the gauge field is given by

$A_{4}(\mathbf{x},x_{4})=-i\frac{1}{\alpha}\sum\limits_{i=1}^{k}R_{i}\tau_{a}\partial_{a}\quantity[\frac{\lambda_{i}^{2}m}{|x-X_{i}|}K_{1}(m|x-X_{i}|)]R^{\dagger}_{i},$
(43)

which yields the correct number of collective coordinates on each topological
sector. To what extent the gauge field configuration given in eq. 43 does
actually reproduce the moduli space of multi-Skyrmions in the massive pion
case after calculating its holonomy along $x_{4}$ is out of the scope of this
paper. Nevertheless, it is an interesting calculation that may allow us to
study the effect of a nonzero pion mass in the Skyrmion-Skyrmion interactions.
Such configuration, in the $k\rightarrow\infty$ limit could be used as well to
approximate the (half-)Skyrmion crystals, since it is the constrained-
instanton version of that proposed in Park  _et al._ (2002). Indeed, a nonzero
pion mass provides a natural cutoff for the individual instantons, which decay
exponentially, and the holonomy can be calculated directly without needing to
impose an artificial cutoff on instanton tails.

##  V Conclusions and outlook

We have presented an approximation of Skyrmions in a model with massive pions
using the holonomy of a constrained instanton of a Yang-Mills-Higgs theory in
4d Euclidean space. Constrained instantons are not even true minima of the
Yang-Mills-Higgs action (due to the presence of a destabilizing scaling mode),
and hence they are not unique. Nevertheless, they preserve their topological
structure, so may still be useful to generate nontrivial Skyrmions via the
Atiyah-Manton construction. Indeed, we proposed two of these configurations,
namely, the linearized solutions of the equations of motion in their
asymptotic form and the constrained instanton found in Wang (1994), which are
both easy to write and still allow to describe the correct asymptotic
exponential behavior of the corresponding Skyrmion configuration.

Therefore, we have shown that the relation between Euclidean Yang-Mills theory
and the Skyrme model can be extended to account for massive pions by adding a
scalar field with a Higgs potential to the Yang-Mills action. An interesting
open question is whether this fact is a mere mathematical curiosity, or, on
the contrary, it can be deduced from a holography-based argument, in the
spirit of Sutcliffe’s simple model. If this were the case, our finding may
lead to new ideas on how to construct more realistic, bottom-up holographic
duals for QCD that are able to describe explicit chiral symmetry breaking in
the boundary theory.

Obviously, having a short, semi-analytical approximation to the $k=1$ Skyrmion
with massive pions can also be useful to study some dynamical aspects of the
model without the need to rely in numerical methods to solve the corresponding
Euler Lagrange equation. We would like to remark that a different method of
generating Skyrmions with massive pions via computing the holonomy of
instantons along circles in $\mathbb{R}^{4}$ was proposed in Atiyah and
Sutcliffe (2005). The field configurations obtained this way correspond to
Skyrmions with massless pions in hyperbolic space, which are shown to be a
good approximation to Skyrmions in flat space with a nonzero pion mass, when
the curvature parameter is related to the pion mass in a certain way. This
idea results also in a very good approximation to the real solutions, but a
drawback of this method is that the moduli space structure is not preserved,
since the position of the skyrmion is related to a fixed point in the base
manifold and it is not a free parameter (as there is not translation symmetry
in hyperbolic space). However, the instanton approximation is most useful for
studying low energy dynamical processes involving two or more Skyrmions, since
it allows to (approximately) describe the corresponding moduli space. In
principle, the extension of our approach to describe a multi-skyrmion moduli
space is possible, although not straightforward. A first step would be to find
constrained multi-instanton solutions, and study the corresponding moduli
space. In this paper, we have proposed a simple configuration that is based on
’t Hooft multi-instanton configuration, which is just a superposition of
oriented constrained instantons. We believe this is a promising path towards a
more realistic treatment of nuclear dynamics using the Skyrme model in the
instanton approximation, since the effects of a nonzero pion mass should play
a nontrivial role in Skyrmion-Skyrmion interactions.

## Acknowledgements

I would like to thank the organizers of the conferences SIG X in Krakow and
Geometric Models of Nuclear Matter in Kent, for providing a friendly
environment that allows the exchange of fresh ideas in the soliton community.
I’d like to thank also J. Evslin for introducing to me the concept of
constrained instantons, and to A.Wereszczynsky and M. Huidobro for reading the
manuscript and providing useful comments. This work has been done thanks to
the funding of my predoctoral research activity from the Spanish Ministry of
Science, Innovation and Universities, and the European Social Fund. This work
has also received financial support from Xunta de Galicia (Centro singular de
investigación de Galicia accreditation 2019-2022), by European Union ERDF, and
by the “María de Maeztu” Units of Excellence program MDM-2016-0692 and the
Spanish Research State Agency.

## References

  * ’t Hooft (1986) G. ’t Hooft, “How Instantons Solve the U(1) Problem,” Phys. Rept. 142, 357–387 (1986).
  * Diakonov (1996) D. Diakonov, “Chiral symmetry breaking by instantons,” Proc. Int. Sch. Phys. Fermi 130, 397–432 (1996), arXiv:hep-ph/9602375 .
  * Skyrme (1961) T. H. R. Skyrme, “A Nonlinear field theory,” Proc. Roy. Soc. Lond. A 260, 127–138 (1961).
  * Atiyah and Manton (1989) M. Atiyah and N. Manton, “Skyrmions from instantons,” Physics Letters B 222, 438–442 (1989).
  * Atiyah, Hitchin, and Singer (1977) M. F. Atiyah, N. J. Hitchin, and I. M. Singer, “Deformations of Instantons,” Proc. Nat. Acad. Sci. 74, 2662–2663 (1977).
  * Jackiw and Rebbi (1977) R. Jackiw and C. Rebbi, “Degrees of freedom in pseudoparticle systems,” Physics Letters B 67, 189–192 (1977).
  * Leese and Manton (1994) R. A. Leese and N. S. Manton, “Stable instanton generated Skyrme fields with baryon numbers three and four,” Nucl. Phys. A 572, 575–599 (1994).
  * Hosaka, Oka, and Amado (1991) A. Hosaka, M. Oka, and R. D. Amado, “Skyrmions and their interactions using the Atiyah-Manton construction,” Nucl. Phys. A 530, 507–531 (1991).
  * Sutcliffe (1994) P. M. Sutcliffe, “Instanton moduli and topological soliton dynamics,” Nucl. Phys. B 431, 97–118 (1994), arXiv:hep-th/9408168 .
  * Cork (2018) J. Cork, “Skyrmions from calorons,” JHEP 11, 137 (2018), arXiv:1810.04143 [hep-th] .
  * Cork and Halcrow (2022) J. Cork and C. Halcrow, “ADHM skyrmions,” Nonlinearity 35, 3944–3990 (2022), arXiv:2110.15190 [hep-th] .
  * Dunajski (2013) M. Dunajski, “Skyrmions from gravitational instantons,” Proc. Roy. Soc. Lond. A 469, 20120576 (2013), arXiv:1206.0016 [hep-th] .
  * Atiyah and Manton (1993) M. F. Atiyah and N. S. Manton, “Geometry and kinematics of two skyrmions,” Commun. Math. Phys. 153, 391–422 (1993).
  * Halcrow and Winyard (2021) C. Halcrow and T. Winyard, “A consistent two-skyrmion configuration space from instantons,” JHEP 12, 039 (2021), arXiv:2103.15669 [hep-th] .
  * Halcrow and Harland (2022) C. Halcrow and D. Harland, “Nucleon-nucleon potential from instanton holonomies,” (2022), arXiv:2208.04863 [hep-th] .
  * Sutcliffe (2010) P. Sutcliffe, “Skyrmions, instantons and holography,” JHEP 08, 019 (2010), arXiv:1003.0023 [hep-th] .
  * Sakai and Sugimoto (2005) T. Sakai and S. Sugimoto, “Low energy hadron physics in holographic QCD,” Prog. Theor. Phys. 113, 843–882 (2005), arXiv:hep-th/0412141 .
  * Espinosa (1990) O. Espinosa, “High-Energy Behavior of Baryon and Lepton Number Violating Scattering Amplitudes and Breakdown of Unitarity in the Standard Model,” Nucl. Phys. B 343, 310–340 (1990).
  * Rubakov and Shaposhnikov (1996) V. A. Rubakov and M. E. Shaposhnikov, “Electroweak baryon number nonconservation in the early universe and in high-energy collisions,” Usp. Fiz. Nauk 166, 493–537 (1996), arXiv:hep-ph/9603208 .
  * ’t Hooft (1976) G. ’t Hooft, “Computation of the quantum effects due to a four-dimensional pseudoparticle,” Phys. Rev. D 14, 3432–3450 (1976).
  * Note (1) We are using the conventions in Espinosa (1990).
  * Affleck (1981) I. Affleck, “On constrained instantons,” Nuclear Physics B 191, 429–444 (1981).
  * Nielsen (2005) N. K. Nielsen, “Instanton constraints in supersymmetric gauge theories. II. N = 2 Yang-Mills theory,” (2005), arXiv:hep-th/0503120 .
  * Note (2) Also, the constraint must allow finiteness of the action, which is a nontrivial demand for choosing a valid constraint, see Nielsen and Nielsen (2000).
  * Shifman (2022) M. Shifman,  _Advanced Topics in Quantum Field Theory: A Lecture Course_, 2nd ed. (Cambridge University Press, 2022).
  * Belavin  _et al._ (1975) A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Y. S. Tyupkin, “Pseudoparticle Solutions of the Yang-Mills Equations,” Phys. Lett. B 59, 85–87 (1975).
  * Klinkhamer (1993) F. R. Klinkhamer, “Existence of a new instanton in constrained Yang-Mills Higgs theory,” Nucl. Phys. B 407, 88–114 (1993), arXiv:hep-ph/9306208 .
  * Nielsen and Nielsen (2000) M. Nielsen and N. K. Nielsen, “Explicit construction of constrained instantons,” Phys. Rev. D 61, 105020 (2000), arXiv:hep-th/9912006 .
  * Wang (1994) M.-y. Wang, “On constrained instanton valleys,” (1994), arXiv:hep-ph/9502333 .
  * Atiyah  _et al._ (1978) M. F. Atiyah, N. J. Hitchin, V. G. Drinfeld, and Y. I. Manin, “Construction of Instantons,” Phys. Lett. A 65, 185–187 (1978).
  * Park  _et al._ (2002) B.-Y. Park, D.-P. Min, M. Rho, and V. Vento, “Atiyah-Manton approach to skyrmion matter,” Nucl. Phys. A 707, 381–398 (2002), arXiv:nucl-th/0201014 .
  * Atiyah and Sutcliffe (2005) M. Atiyah and P. Sutcliffe, “Skyrmions, instantons, mass and curvature,” Physics Letters B 605, 106–114 (2005).


Safari Escape

# Multi-Agent Path Planning with Asymmetric Interactions In Tight Spaces

V. Modi1 0000-0002-9350-494X Y. Chen 10000-0001-7547-9587, A. Madan
10000-0001-7547-9587, S. Sueda 20000-0003-4656-498X, and D. I. W. Levin
10000-0001-7079-1934
1University of Toronto, Toronto, Canada 2Texas A&M University, College
Station, TX

###### Abstract

By starting with the assumption that motion is fundamentally a decision making
problem, we use the world-line concept from Special Relativity as the
inspiration for a novel multi-agent path planning method. We have identified a
particular set of problems that have so far been overlooked by previous works.
We present our solution for the global path planning problem for each agent
and ensure smooth local collision avoidance for each pair of agents in the
scene. We accomplish this by modeling the trajectories of the agents through
2D space and time as curves in 3D. Global path planning is solved using a
modified Djikstra’s algorithm to ensure that initial trajectories for agents
do not intersect. We then solve for smooth local trajectories using a quasi-
Newton interior point solver, providing the trajectory curves with a radius to
turn them into rods. Subsequently, resolving collision of the rods ensures
that no two agents are in the same spatial position at the same time. This
space-time formulation allows us to simulate previously ignored phenomena such
as highly asymmetric interactions in very constrained environments. It also
provides a solution for scenes with unnaturally symmetric agent alignments
without the need for jittering agent positions or velocities.

<ccs2012> <concept> <concept_id>10010405</concept_id> <concept_desc>Applied
computing</concept_desc> <concept_significance>500</concept_significance>
</concept> </ccs2012> [500]Computing Methodologies Applied computing

##  1 Introduction

On a hot dry day in Kruger National Park, an empty truck idles on the side of
a road. Sam, the driver of the truck, has wandered a hundred meters off the
road in an attempt to take a picture of a tree with ten baboons. The baboons
suddenly jump out of the tree and charge towards Sam. Sam panics and starts
running back to the truck; however, coincidentally an elephant wandering in
the area is on a path perpendicular to Sam’s, between him and the truck. How
will Sam get back to the truck safely while outrunning the baboons and
avoiding collision with the elephant? What are the paths of the twelve agents
in the scene: Sam, the elephant and the ten baboons?

The problem above poses many challenges to a multi-agent path planning
algorithm. First Sam must anticipate collisions ahead of time, in order to
move quickly and efficiently to their truck. Second, the three groups of
agents, Sam, the baboons and the elephant have dramatically different masses
and behaviors. For instance while Sam seeks to avoid all animals, the massive
elephant is untroubled, and will stubbornly continue on its path. Finally,
geographic features such as additional trees and ponds can lead to a highly
constrained environment. Existing state-of-the-art crowd simulation methods
struggle to compute anticipatory agent paths in constrained environments when
asymmetric interactions are involved (Table 1) making them ill-suited for
application in planning problems such as the example given above.

We propose a new multi-agent path planning algorithm well suited for these
problems. Our model directly optimizes the space-time trajectories of all
agents which allows for per-agent physical and psychological characteristics
and smooth anticipatory trajectories. A novel, differentiable space-time
repulsive energy ensures collision free trajectories. Using our approach, Sam
arrives safely at the truck, escaping the baboons and avoiding the elephant.

Figure 1:  Sam needs to escape ,from the baboons while avoiding the elephant,
whereas the elephant is unconcerned with the other agents in the scene.

##  2 Related Work

Successful multi-agent path planning requires an algorithm to both correctly
model the behavior of independent agents as well as their interactions. A
common approach is to apply a dynamics model based on Newton’s second law of
motion which is integrated over time to produce plausible agent trajectories.
Psychological and social characteristics of the group can be incorporated into
the dynamics equations as social forces [POSB05, HM98, WLJT17]. Intra-agent
forces, such as those that handle collision avoidance are added through a
variety of means, and we can partition the space of successful approaches
based on locality of their models in both space and time. Approaches such as
[GKLM11] incorporate psychological factors into an underlying dynamics model,
but it is impossible to tell apart the psychological traits of the agents by
simply observing the simulation. Our method makes the impact of agent
characteristics on the trajectory very obvious, while providing a standalone
local and global dynamics model for the scene.

In local methods, agent decision making requires information only local in
space and in time. Local collision avoidance methods [VGLM11, GNCL14, WLP16]
compute collision response using local information in space and time (the
planning horizon can be as small as a single time step). These methods
struggle to generate smooth anticipatory collision responses. Implicit Crowds
[KSNG17] attempts to overcome this difficulty by introducing a time-to-
collision potential to the crowd dynamics. This gives agents richer space-time
information on which to act; however this energy is effectively local,
computed from the current state of the system (position and velocity). NH-TTC
[DKG20] improves upon prior work by utilizing longer planning horizon,
geometrically represented by curves in space. Intersection checks between
these curves allow agents to react to collisions likely to occur in the near
future. Similarly, vision based methods, such as [OPOD10, DMC*17] provide an
anticipatory collision avoidance model, but with extremely limited path
planning capabilities and no guarantee of smooth agent motion. Data-driven
aproaches such as [CC14] use an underlying state-action graph created via
external data to generate trajectories. However, these trajectories are highly
data-dependent, do not factor in environmental constraints, seem to operate in
sparse crowds and must be used in conjunction with some other higher level
global path planner.

An alternative approach to more local methods is to extend the collision
response globally in space using a fluid like pressure solve [Hug02, TCP06,
NGCL09]. However because these methods only consider the configuration of the
system (position and velocity of each agent) at a single time, their ability
to produce smooth anticipatory collision response, especially in a sparser
setting, is reduced. Continuum Crowds stands out as one of the only methods
that incorporates both a global planner through Djikstra’s search and local
collision avoidance through a fluid-like pressure solve. Another option is to
handle local collision avoidance using RVO or Social Forces, as done by the
data-driven method [TYK*09]. Optionally, one might use [HKHL13], another data-
driven method which uses a purely stochastic collision avoidance method. None
of these approaches solve the problems of handling tight environmental
constrains or asymmetric agent interactions. In fact, an adjacent field of
research, that involves measuring the "correctness" of various crowd models (
[GVL*12, WJO*14]) also fails to test for asymmetric agents and behavior in
complex environments.

Figure 2: The repulsive curves energy forces agents to be maximally far away
from one another which is not an intuitive behavior.

Finally, while not strictly designed for multi-agent path planning, Repulsive
Curves [YSC21] can potentially be used for agent planning. While it yields
impressive results in 3D curve untangling, when applied to multi-agent path
planning it has several drawbacks. First, in unconstrained environments,
Repulsive Curves will maximally repell agents away as shown in Figure 2.
Secondly, like other previous methods, there is no straightforward way to
model asymmetric agents. Lastly, Repulsive Curves uses random "jitter" to
ensure that no trajectories initially overlap. From our own experiments, we
have observed that this is not sufficient to guarantee non-overlapping initial
trajectories. Meanwhile, our space-time Djikstra’s approach is guaranteed to
create initially non-intersecting trajectories, but is not sufficient for
overcoming the other limitations of Repulsive Curves.

In this paper we present a multi-agent path planner that computes its response
globally in space and time. Rather than formulating our multi-agent planner as
an initial value problem, we instead take a space-time optimization approach.

Space-time trajectory optimization [WK88] was initially applied to keyframe
interpolation, with subsequent methods such as [PSE*00] allowing manual
editing of object trajectories and [SKF08] allowing interactive motion
correction and synthesis using graph search methods. The concept of using
space-time graph search methods for collision avoidance is furthered in
[LLKP11] for a single agent and [SKH*11] for small groups. For larger
situations, standard space-time approaches would be computationally
prohibitive. Our method builds on these prior methods by expanding the look-
ahead globally and ensuring smooth paths by collision resolution on the entire
trajectories rather than using a limited space-time graph search approach.

Our method treats each agent as an individual space time curve with only three
requirements: a starting position and time and an ending position (at
indeterminate time). We use a globally supported, differentiable LogSumExp
smooth distance in space-time to guarantee collision free trajectories and
solve the resulting problem using an interior-point technique. For small to
medium scale crowds, our method outperforms current state-of-the-art methods
in simple scenarios (3(a) and 3(b)) and more complicated, constrained
environments. Our agents exhibit anticipatory behaviors such as slowing and
waiting to let others pass and makes no assumptions about agents having
identical mass or other physical properties. This enables planning of
intricate scenes with environmental constraints such as the one described in
our introduction.

(a) 3 agents

(b) 8 agents

Figure 3: On the left, circles denote the start position of the agents and starts denote the desired end position. Only our method lets agents reach their desired ends without any locking. On the right, each agent wants to traverse to the directly opposite end of the circle. Our method leads to the most efficient trajectories. Features | Ours | (1) RVO | (2) IC | (3) NH-TTC | (4) RC | (5) CC | (6) STPL | (7) AABS
---|---|---|---|---|---|---|---|---
Asymmetric Interactions | Y | N | P | P | N | N | N | N
Extremely Constrained Environments | Y | N | N | N | Y | N | P | Y
Smooth Local Interactions | Y | N | Y | N | Y | N | Y | N
Intuitive Control Parameters | Y | Y | N | P | N | N | N | Y
Multiple Agents | Y | Y | Y | Y | Y | Y | N | Y
Code Available Online | Y | Y | Y | Y | Y | N | N | N
Table 1: Compare (1) RVO ( [VGLM11]), (2) Implicit Crowds ([KSNG17]), the
brand new time-to-collision method (3), NH-TTC ([DKG20]), (4) Repulsive Curves
([YSC21]), (5) Continuum Crowds ([TCP06]), (6) Space Time Planning With
Parameterized Locomotion Control ([LLKP11]), and (7) Modular Framework for
Adaptive Agent Base Steering ( [SKH*11]). Y - feature is available, N -
feature is not possible, P - feature might be possible, but not demonstrated.

##  3 Method

At a high level we are influenced by the notion of correlated equilibrium. In
a correlated equilibrium solution to a non-cooperative game, an "oracle"
chooses a strategy for each player, and no player has any reason to deviate
from the chosen strategy assuming others do not deviate either. The result is
an equilibrium solution which maximizes collective utility.

Our method acts as the oracle of the scene and plans agent trajectories in a
way that collectively maximizes the utility of the entire scene. Additionally,
unlike previous approaches which ‘pre-set’ trajectories for agents (often done
manually), we allow our agents to find their own utility-maximizing
trajectories. Lastly, real life agents have different sizes, masses, and
personalities, which affect the agents’ utilities, and therefore its path as
well. Our local-global path planning approach allows for these nuances.

Figure 4:  Top: A path through space can be represented as a curve embedded in
$R^{2}$. Every path has a scalar cost (utility) value. For a cost function
that minimizes distance travelled, $c_{1}$ and $c_{2}$ are inoptimal paths
from the start to the goal (red X), but $c_{3}$ is optimal. Bottom: A path
through space and time can be represented as a curve in $R^{3}$. A space-time
rod (as shown here) is simply a 3D curve with a collision radius.

###  3.1 Paths Map To Utility

As shown in Figure 4, for an agent, each path from the start location to the
end location through space maps to some scalar utility value. For example, if
the agent’s utility minimized distance traveled, a straight line from start to
end would prove to be the optimal path. In Figure 4, the 2D x-y plane
describes the domain of spatial motion for agent paths. We make the additional
observation that agents do not simply move through space, but also move
through time which is denoted by the vertical z-axis. So, in 3D space-time,
the agent’s motion through space ($x,y$) and time ($t$) is described by a 3D
curve which corresponds to a utility. For any given moment in time, the
projection of the curve onto the x-y plane gives us the agent’s location.

Let us consider a simple case where there are no other agents in the scene and
the terrain is flat and free of obstacles. The agent, indexed henceforth by
subscript $a$, wishes to maximize its utility (by minimizing cost,
$\bar{\Psi}$); the tilde symbol indicates the variable or function is
continuous, not discrete. Lower case variables indicate a single agent. Upper
case indicates the variable aggregates all agents. The motion of the agent
through space and time denoted by the agent’s 3D space-time curve
$\tilde{\mathbf{s}}_{a}$ embedded in $R^{x,y,t}$ requires constraints.
Fortunately all our constraints are linear, so we lump them together into
$\tilde{B}$ for now. Put together, the full optimization for a single agent
$a$ is

$\displaystyle f_{a}^{*}$ $\displaystyle=\min\;\int_{\tilde{\mathbf{s}}_{a}}$
$\displaystyle\tilde{\psi}(\tilde{\mathbf{s}}_{a})\text{d}\tilde{\mathbf{s}}_{a}$
(1a) $\displaystyle s.t.$
$\displaystyle\tilde{B}\tilde{\mathbf{s}}_{a}\leq\mathbf{0}.$ (1b)

The cost, Equation 1a, for agent $a$ is integrated over the trajectory of the
agent $(\tilde{\mathbf{s}}_{a})$. The specific nature of the cost function
determines the agents’ behavior based on its characteristics. An agent might
have a preferred walking speed, or a stubbornness factor, or a radius of
comfort, all of which (and more) can be encoded into components of
$\tilde{\psi}$. In order to solve this optimization, we must first discretize
the agent’s trajectory into the discrete curve $\mathbf{s}_{a}$ shown in
Figure 5. We also discretize the cost $\tilde{\psi}$ into separate intra-agent
costs, which solely affect the path of one individual, and interaction costs,
which can affect multiple agents.

###  3.2 Discretizing

Our agent’s cost function takes in the agent’s discrete 3D space-time curve
$\mathbf{s}_{a}$ as input and outputs a scalar cost for that path, $f_{a}$. We
descretize the agent’s continuous $\tilde{\mathbf{s}}$ into a piecewise linear
curve described by $n+1$ nodes
$\mathbf{s}_{a}=[(x_{a}^{0},y_{a}^{0},t_{a}^{0}),\allowbreak...,\allowbreak(x_{a}^{i},y_{a}^{i},t_{a}^{i}),\allowbreak....,\allowbreak(x_{a}^{n},y_{a}^{n},t_{a}^{n})]$
for nodes $i=0..n$ connected sequentially by edges. We must also discretize
our constraints. First, even though time is a variable in our formulation,
Equation 2b ensures the agent cannot move backwards in time. Second, the agent
has a start location $(x_{a}^{0},y_{a}^{0})$ and start time $t_{a}^{0}$
denoted in Equation 2c. Third, Equation 2d sets a goal or an end location
$(x_{a}^{n},y_{a}^{n})$. Lastly, the agent cannot take an infinite amount of
time, so Equation 2e bounds the agent by a max time $T_{a}^{max}$. Putting all
of this together, we can re-write the optimization for agent $a$ with the
discrete generalized cost function $\psi$ as

$\displaystyle f_{a}^{*}$ $\displaystyle=\min\;\sum_{i=0}^{n}$
$\displaystyle\psi(\mathbf{s}_{a})$ (2a) $\displaystyle s.t.$ $\displaystyle
t_{a}^{i}\leq t_{a}^{i+1}$ (2b)
$\displaystyle(x_{a}^{0},y_{a}^{0},t_{a}^{0})=\mathbf{(x_{a},y_{a},t_{a})^{start}}$
(2c) $\displaystyle(x_{a}^{n},y_{a}^{n})=\mathbf{(x_{a},y_{a})^{end}}$ (2d)
$\displaystyle t_{a}^{n}\leq T_{a}^{max}.$ (2e)

Now with a template for our optimization problem, we can replace the
generalized cost function with specific discrete costs. Our specific cost
functions are derived from observation of real behavior. These behavioral
observations are divided into two categories, intra-agent costs, and
interactions costs. Intra-agent cost functions only look at one agent’s path
at a time. Interaction cost functions include avoiding collisions with other
agents, collisions with static obstacles such as walls or furniture, grouping
behavior within friends, asymmetric behavior based on the mass, size of the
agents, or other characteristics.

Figure 5: This space-time rod describes a discrete space-time trajectory curve
made up of many nodes. Our agent moves along the projection of the rod from
its start to end, taking into account the boundary constraints of the problem.

###  3.3 Intra-Agent Costs

Since these costs apply to a single agent, we calculate intra-agent costs for
arbitrary agent $a$ and later we show how to sum the costs over the entire
scene. The path of our agent, $\mathbf{s}_{a}$, is comprised of of $n+1$ nodes
indexed by $i=1..n$. Each node
$(x^{i}(\mathbf{s}_{a}),y^{i}(\mathbf{s}_{a}),t^{i}(\mathbf{s}_{a}))$ is
comprised of the agent’s spatial $(x,y)$ coordinates and time $(t)$
coordinates. For all cost functions, agent $a$ also has constant weighting
terms, $K_{a}$, as well as constant characteristics such as mass, $m_{a}$, and
preferred end time, $T_{a}^{p}$.

####  3.3.1 Intra-Agent Kinetic Cost

Our agent $a$ will try to minimize energy expenditure. Given an agent’s space-
time curve $\mathbf{s}_{a}$, the kinetic energy cost of the agent over the
curve can be written as

$\displaystyle
C_{a}^{K}(\mathbf{s}_{a},K_{a}^{K})=K_{a}^{K}\sum_{i=0}^{n}\frac{1}{2}m_{a}\frac{(\Delta\mathbf{x}_{a}^{i})^{T}(\Delta\mathbf{x}_{a}^{i})}{(\Delta
t_{a}^{i})^{2}}(\Delta t_{a}^{i})$ (3)
$\displaystyle=K_{a}^{K}\sum_{i=0}^{n}\frac{1}{2}m_{a}\frac{(x^{i+1}(\mathbf{s}_{a})-x^{i}(\mathbf{s}_{a}))^{2}+(y^{i+1}(\mathbf{s}_{a})-y^{i}(\mathbf{s}_{a}))^{2}}{t^{i+1}(\mathbf{s}_{a})-t^{i}(\mathbf{s}_{a})}$
(4)

where Equation 3 is the kinetic energy $\frac{1}{2}mv^{2}$ for curve segment
$e_{a}^{i}=[x_{a}^{i},y_{a}^{i},t_{a}^{i},x_{a}^{i+1},y_{a}^{i+1},t_{a}^{i+1}]$
integrated over total time travelled $\Delta t_{a}^{i}=t_{a}^{i+1}-t^{i}$.
This simplifies into Equation 4 where $K_{a}^{K}$ is the agent-wise weighting
coefficient, $m_{a}$ is the mass of the agent and $\mathbf{s}_{a}$ is the
discretized path curve (our input variable) made up of $n+1$ nodes.

####  3.3.2 Intra-Agent Acceleration Cost

In order to penalize acceleration in agent $a$’s trajectory, the acceleration
cost

$\displaystyle
C_{a}^{A}(\mathbf{s}_{a},K_{a}^{A})\;=K^{A}_{a}\sum_{i=1}^{n-1}\frac{1}{2}(\theta^{i}_{a})^{2}$
(5)

measures the curvature of $\mathbf{s}_{a}$ through the discretized
acceleration cost where angle $\theta$ is the angle between two piece-wise
linear segments of the curve computed using the stable arctan function
$arctan2$. The angle is $\theta_{a}^{i}\;=\allowbreak
arctan2(\frac{\|(\mathbf{s}_{a}^{i+1}-\mathbf{s}_{a}^{i})\times(\mathbf{s}_{a}^{i}-\mathbf{s}_{a}^{i-1})\|}{(\mathbf{s}_{a}^{i+1}-\mathbf{s}_{a}^{i})^{T}(\mathbf{s}_{a}^{i}-\mathbf{s}_{a}^{i-1})})$
where $\mathbf{s}_{a}^{i}=[x^{i},y^{i},t^{i}]$ for each node in the curve.

####  3.3.3 Intra-Agent Preferred End Time Cost

Any agent $a$ has a preferred end time, $T_{a}^{p}$ at which they expect to
reach the end position. Sometimes this is the same as the max end time
$T_{a}^{max}$ by when the agent is required to be at the end positions, but
sometimes the preferred end time $T^{p}_{a}$ can be sooner. Deviation from the
preferred end time is modeled as a quadratic cost

$\displaystyle
C_{a}^{T}(\mathbf{s}_{a},T_{a}^{p},K_{a}^{T})\;=K_{a}^{T}\frac{1}{2}(t_{a}^{n}-T_{a}^{p})^{2}$
(6)

incentivizing the agents to arrive at their preferred end time $T_{a}^{p}$ by
keeping the actual end time for each agent, $t_{a}^{n}$, close to $T^{p}_{a}$.

####  3.3.4 Regularizing Cost

We find that adding a regularizing term to penalize extremely short time
segments improves the overall quality of the paths by reducing near
instantaneous motions in time. To penalize very fast agent motion, we use the
regularizing term

$\displaystyle
C_{a}^{R}(\mathbf{s}_{a},K_{a}^{R})\;=K_{a}^{R}\sum_{i=0}^{n}\left(\frac{\frac{t^{n}_{a}}{n}}{t^{i+1}_{a}-t^{i}_{a}}\right).$
(7)

This is important because it allows us to feed the solver an initial space-
time curve such as the ones in Figure 10 with many near instantaneous agent
motions, and the solver is able to optimize the final space-time curve to a
much more reasonable trajectory.

###  3.4 Interaction Cost

Interactions include agent-environment interactions and agent-agent
interactions for which we must introduce a new arbitrary agent indexed by $b$
where $a\neq b$. All our interactions depend on each agent’s ‘radius of
comfort’ (or collision radius), denoted by $r_{a},r_{b}$ for agents $a$ and
$b$. The collision radius might be based on size, or a combination of size and
personality: for example, people stay away from angry people. We encode these
agent characteristics into our cost functions by extruding circle with radius
$r_{a}$ along the agent’s path $\mathbf{s}_{a}$ thus forming a rod. In a scene
with multiple agents, as long as no two rods are intersecting, the scene is
collision free as shown in Figure 6. The collision radius ensures that no two
agents are in the same place at the same time. So even though the forces
between the two space-time rods might be symmetric, the agents’ response to
these forces (change in path) leads to asymmetric interactions. For static
object interactions, the environment boundary is encoded into boundary
vertices $b_{v}$ and boundary elements $b_{e}$. Let us examine how we detect
which rods are intersecting and how we deal with our three interaction types:
agent-agent collisions, agent-agent groupings (friendships), and agent-
environment collisions.

Figure 6:  Resolving intersections between agents space-time rods resolves
collisions in the scene since no agent shares the same $x,y,t$ location at any
point in the scene.

####  3.4.1 Collision Interaction Cost

Figure 7: Tunneling artifact.

Given a pair of agents $a,b$ with 3D space-time curves
$\mathbf{s}_{a},\mathbf{s}_{b}$ and collision radii $r_{a},r_{b}$, we densely
sample each space-time curve uniformly. Next, we calculate the minimum smooth
distances between the upsampled centerlines using the method described in
subsection 3.5. As long as the smooth minimum distance d between the upsampled
rods
$u_{a}=\text{upsample}(\mathbf{s}_{a}),u_{b}=\text{upsample}(\mathbf{s}_{b})$
is further apart than $r_{a}+r_{b}$, collisions will not occur. We implement
this non-linear constraint on the agent paths as a log barrier energy:

$\displaystyle
C_{a,b}^{C}(\text{d}(\alpha,\text{upsample}(\mathbf{s}_{a}),\text{upsample}(\mathbf{s}_{b})),r_{a},r_{b},K_{a,b}^{C})$
$\displaystyle=-K_{a,b}^{C}\log(-(r_{a}+r_{b})+\text{d})).$ (8)

where $K_{a,b}^{C}$ is the pair-wise weighting coefficient on the energy. The
intuition behind a log barrier energy is explained in Figure 8. Upsampling the
trajectories prevents tunneling artifacts (Figure 7) between edges during
collision detection. The collision cost increases exponentially as the minimum
distance between the two agent rods approaches the sum of the collision radii,
so collisions are exponentially penalized. By summing this cost with the
kinetic energy term, which incentivizes agents to move at a constant velocity,
we get smooth local collision resolution as shown in Figure 6.

We employ two methods for sparsifying interactions. First, we use a 3D tree
structure to store agent trajectory nodes (analogous to a Bounded Volume
Hierarchy BVH) for broad-phase collision detection between trajectories. The
tree-based broad phase reduces collision detection costs from $O(n^{2})$ to
$O(nlogn)$. Second, if two agents are known to be far apart, we entirely avoid
the collision detection step between those two agents thus reducing costs to
$O(n)$ in the best case. Using a combination of tree based broad phase along
with manual denotation of interacting agents, we find that our method scales
nearly linearly $O(n)$as the number of agents in the scene increases.

Figure 8:  The log barrier energy quickly tends towards infinity as the
minimum distance approaches the collision radius.

####  3.4.2 Grouping Interaction Cost

While the collision resolution term pushes agent paths apart, we introduce a
grouping term which pulls paths together

$\displaystyle
C_{a,b}^{G}(\text{d}(\alpha,\text{upsample}(\mathbf{s}_{a}),\text{upsample}(\mathbf{s}_{b})),r_{a},r_{b},K_{a,b}^{G})$
$\displaystyle=K_{a,b}^{G}(\text{d}-(r_{a}+r_{b}))^{2}.$ (9)

This quadratic grouping term allows us to model scenarios in which friends
going in the same direction tend to stick together as shown in Figure 16, or
scenarios in which an agent needs to deliver a message to another agent who is
out of the way from his final destination such as Figure 18. The weighting
term $K_{a,b}^{G}$ controls the desire of the agents $a,b$ to group together.

####  3.4.3 Agent-Environment Interaction

We surpass previous methods in three ways in terms of agent-environment
interactions. Firstly, our method works on agents traversing highly
constrained environments such as 13(a), 13(b), or 17(a). Secondly, our method
allows agents to explore routes when multiple paths are available and choose
the most optimal path using a novel modification to Djikstra’s algorithm.
Thirdly, our path planning step is topology aware, so this method will work on
more fascinating terrain.

Environment maps are stored as 2D mesh files with vertices and faces. We pre-
compute the vertices and edges around the boundaries of the map and any static
obstacles into $b_{v}$ and $b_{e}$ as a one-time preprocessing step. Next, we
compute the minimum smooth distance between the map boundary edges and the
agent rod. Luckily since static obstacles are fixed in space for all time, we
can ignore the third dimension and only compute the 2D minimum distance,
$\text{d}(\alpha,\mathbf{s}_{a},b_{v})$ from the rod’s projection onto the map
and the map’s boundary edges. Next we pass the smooth distance d into the log
barrier function

$\displaystyle
C_{a}^{M}(\text{d}(\alpha,\text{upsample}(\mathbf{s}_{a}),b_{v}),r_{a})=-K^{M}_{a}log(-r_{a}+d)$
(10)

where $K^{M}_{a}$ is the energy coefficient. Like in the agent-agent collision
function (Equation 8), we can use the collision radius $r_{a}$ for each agent
since the initial agent paths are sufficiently far enough away from static
obstacles such that $d>r_{a}$. The coarseness of the map mesh does not impact
the smoothness of agent path; however, Djikstra’s algorithm requires two
vertices on the map which correspond to the agent’s start and end positions to
generate the initial path. Therefore, we set the nearest vertices to the given
start and end points as the boundary conditions for the solve.

###  3.5 Smooth Min Distance Implementation

The absolute minimum distance between agent paths is not a differentiable
function. Therefore, we use a differentiable LogSumExp smooth minimum distance
function to smoothly approximate the distance between two upsampled space-time
curve vertices $u_{a}$ and $u_{b}$. See the supplementary material for
derivatives of the function. The smooth minimum distance function is

$\displaystyle
d(\alpha,u_{a},u_{b})=\frac{-1}{\alpha}\log\left(\sum_{i=0}^{n}\sum_{j=0}^{n}e^{-\alpha\|\mathbf{s}_{a}^{i}-\mathbf{s}_{b}^{j}\|}\right)$
(11)

where the $\alpha$ term controls the numerical sensitivity of the distance
which changes depending on the size of the environment mesh and the distances
between the agents’ path curves. A higher $\alpha$ leads to a more accurate
minimum distance, but reduces numerical stability. Since this smooth minimum
distance function is an underestimation of the true minimum distance, it is
possible to get negative distances for very close objects.

We use 1 to update $\alpha$ during the solve to guarantee a usable (real, non-
negative) smooth minimum distance. Our initial $\alpha_{0}$ for the scene is
as low as possible, but large enough to guarantee that both our log barrier
interaction costs are real and defined. If the initial $\alpha_{0}$ is too
small, the distance underestimation is too great, and the log barrier costs
are either complex or undefined, then the interior point solver will throw an
error. We must choose a large enough $\alpha_{0}$ that the heuristic will
provide an $\alpha$ that guarantees a real smooth distance values where
$D>r_{a}+r_{b}$ to ensure real log barrier costs. For subsequent iterations,
any real, non-negative smooth minimum distance is fine.

Algorithm 1 A heuristic to find a usable $\alpha$ to be used within the agent-
collision and map interaction cost, gradient and hessian functions.

function ALPHAHEURISTIC($\alpha_{0},\mathbf{s}_{a},\mathbf{s}_{b}$)

$\alpha\leftarrow\alpha_{0}$

$D\leftarrow d(\alpha,\mathbf{s}_{a},\mathbf{s}_{b})$

while $D\leq 0$ do

$\alpha\leftarrow\alpha+0.1*\alpha_{0}$

$D\leftarrow d(\alpha,\mathbf{s}_{a},\mathbf{s}_{b})$

end while

return $\alpha$

end function

###  3.6 Initial Paths

Using an interior point solver requires that we find feasible initial
trajectories in which the agents do not collide with other agents or any
obstacles. The supplementary material shows our failed experiments in
generating initial paths based on shortest distances. Some global path planner
is necessary to generate feasible initial paths. Initially intersecting space-
time rods causes NaNs in the log barrier interaction costs which makes the
interior point method intractable. In order to generate initial plausible
paths to feed into the optimization, we sequentially run a modified Djikstra’s
algorithm for all agents traversing the environment mesh inspired by [TCP06].

Djikstra finds the distance-weighted least costly path along the edges of the
map mesh from the start location to the end location for a given agent. We
avoid obstacles (environment boundary) by giving these edges, $b_{e}$, a
prohibitively high traversal cost. We only consider map edges that have a
distance from the map boundary at least greater than the agent’s collision
radius.

Figure 9:  Running Djikstra’s algorithm to set initial agent paths on the map.
Red map edges indicate a higher traversal cost. Since we want to set a
feasible initial path for each agent, we keep them away from the map edges as
well as away from other agent’s initial paths.

Furthermore, we ensure that every agent has a collision-free initial path by
making it impossible for an agent’s initial path to pass too close to the
initial path of another agent. For each agent, upon finding an initial path,
we wipe out all edges connected to the path within the agent’s collision
radius. For very constrained scenes, it becomes impossible to ensure collision
free initial trajectories since agents often run out of traversible edges from
$[x^{0},y^{0}]$ to $[x^{n},y^{n}]$. Since agents move through both space and
time, we extrude our environment map into the $t$ dimension and create a pre-
set number of layers from $t^{0}$ to $T^{max}$ as shown in Figure 9. Solving
this 3D space-time Djikstra’s algorithm makes it a lot more feasible to find
collision-free inital paths for all agents. Sometimes, the coarseness of the
environment mesh, or the low number of layers in the 3D space-time map makes
it impossible for Djikstra’s algorithm to find a feasible initial path. In
this case, either the environment map would need to be subdivided, more layers
would need to be added, or the scene itself might be infeasible given the
$T^{max}$ time constraints and the number and size of agents involved.

Figure 10:  Top is the initial space-time curve output from Djikstra’s path
finding algorithm. Bottom is the final result of the solver with optimal
paths. Agents are allowed to find their own optimal end-times while avoiding
collisions in tightly constrained bottlenecks such as this.

Additionally, although the space-time Djikstra’s algorithm will provide
feasible initial paths, they will have many kinks and sharp turns (Figure 10)
indicating unnatural motion and thus cannot be used in the simulation
directly. These get smoothed out during the optimization process resulting in
much more realistic motion. It is possible that this proposed method to
generate initial paths will affect the state of the final trajectories. For
example, if an initial path passes through the left of an obstacle, it might
be stuck in a local minimum even though an optimal trajectory might be through
the right side. This drawback is somewhat mitigated by using Djikstra’s search
to pre-compute the least costly paths. Either way, whether or not the solver
produces a globally minimum outcome, the trajectories are guaranteed to be
smooth, viable and collision free.

####  3.6.1 Optimization

Now aggregating over all the agents in the scene, $\|A\|$, we put intra-agent
cost functions Equations 4-7 together with interaction costs Equation 8-10 and
construct an optimization problem for our scene. We aggregate boundary
conditions for each agent and denote them with capital letters to show that
they apply to the entire scene. The optimization problem

$\displaystyle f^{*}=\min\;$
$\displaystyle\sum_{a=1}^{\|A\|}C_{a}^{K}(\cdot)+C_{a}^{A}(\cdot)+C_{a}^{T}(\cdot)+C_{a}^{R}(\cdot)+C_{a}^{M}(\cdot)+$
(12a) $\displaystyle\sum_{\begin{subarray}{c}a,b\epsilon\\\ \\{1..\|A\|\\}\\\
|a\neq b\end{subarray}}C_{a,b}^{G}(\cdot)+C_{a,b}^{C}(\cdot)$ (12b)
$\displaystyle s.t.$ $\displaystyle T^{i}\leq T^{i+1}$ (12c)
$\displaystyle(X^{0},Y^{0},T^{0})=(X,Y,T)^{start}$ (12d)
$\displaystyle(X^{n},Y^{n})=(X,Y)^{end}$ (12e) $\displaystyle T^{n}\leq
T^{max}$ (12f)

is solved using a quasi-newton interior point method. An outer loop over this
optimization is required for a proper log-barrier method as explained by
[NW06]. The full pseudocode overview of our method’s outer loop, 2,
illustrates this. The supplementary material for our method include the
subroutines, the calculations for the gradients and hessian approximations of
our costs, intuition for controlling agent behavior through function weights.

Algorithm 2 A full pseudocode overview of our method with subroutines provided
in the supplementary material.

global input variables

$OuterIts\geq 1$, outer loop iterations

$\mu,c=0.75$, log-barrier coeff and its decrement factor

$\alpha_{0},cutoff=0.2$, initial alpha, Hess sparsifying cutoff

$T^{max}$, max time constraint for agents

$b_{e},b_{v}$, environment edges and vertices.

$K$, cost function weights for agents

$R$, collision radius for agents

$T$, preferred end times

$M$, agent masses

end global input variables

interior point solver parameters

$MaxIts$, max iterations

$B$, trajectory boundary conditions

Stop if $StepSize>10^{-2}(meters)$, norm of solver step

Stop if $FirstOrderOpt>10^{-2}(meters)$, first order optimality criteria

end interior point solver parameters

@COSTS, @GRADS, @HESS (supplementary material)

function FINDOPTIMALPATHS

$[\mathbf{S},B]\leftarrow DJIKSTRASPREPROCESS(R,T,b_{v})$

$\mathbf{S}\leftarrow$ increment $t$ by $\epsilon$ to ensure $t_{i+1}>t_{i}$

while OuterIts>0 do

$\mathbf{S}\leftarrow IPSOLVER(\mathbf{S},B,@COSTS,@GRADS,@HESS)$

$\mu\leftarrow c\mu$

$OuterIts\leftarrow OuterIts-1$

end while

return $\mathbf{S}$

end function

Figure 11: (left) we varying number of agents with 300 DOFs each. (middle) we
vary the number of DOFs/agent from 90 to 540 on a fixed scene with 8 agents.
(right) we vary DOFs in the environment boundary from 2520 to 20160 on the
corn maze mesh with 3 agents of 600 DOFs each in the scene. Scaling tests show
near linear performance in the worst-case and linear performance in the best
case.

##  4 Results

THe performance of our method relies on (1) the number of agents in the scenes
and (2) the complexity of the agent’s paths through space and time. We push
along each of these axes separately in this section. We show extremely complex
maze-like environments with bottlenecks, as well as large scenes with tens (to
hundreds) of agents. In addition to the comparisons shown in the related
works, we include the obligatory circles of agents Figure 12. Notice that our
agents take smooth, natural trajectories rather than the strange spinning
motions demonstrated by other methods. Next, in Figure 16 we highlight the
several different types of asymmetric interactions supported by our method. We
show the naive case where agents collide. We show standard symmetric collision
between agents. We show a stubborn snake which forces the scared humans to
move out of its path. We show a large elk who is still scared of the humans so
it changes its trajectory just as the human agents change theirs. We show a
large and stubborn elephant that goes straight to its destination while the
humans have to significantly alter trajectories to avoid it. Finally, we show
a scenario where two friends stick closer to each other on their way to their
final destination.

Figure 12: Circles of agents.

In Figure 13 and Figure 14, we show our method works in extremely tight
environments. RVO and Implicit Crowds show that they can navigate
environmental constraints with a lot of manual pre-processing, but nowhere
near as tight as our examples. Meanwhile, NHTTC does not implement
environmental constraints altogether. We show that our method works in
extremely constrained scenarios such as a subway tunnel or a tight warehouse
of robots where each lane is big enough to accommodate only a single agent and
agents have to take turns to pass through to avoid locking.

(a) Maze of roombas.

(b) Subway platform

Figure 13: In both (a) and (b) only one agent can move through the
intersection at a time.

Although our method is built to handle asymmetric interactions and agents in
highly constrained spaces, we also show several other features of our method.
First, our agents have flexible arrival times. In Figure 10 we show a
bottleneck where all the agents simply cannot arrive at their destinations by
their preferred end times. This feature serves as a counterpoint to the idea
of using fixed start and end boundary conditions to model trajectory curves.
In Figure 14 we show a highly constrained bottleneck of agents trying to get
to their seats on an airplane.

Figure 14: Chaos in first class. All the passengers are in incorrect seats.
Figure 15: A large scene with 102 agents moving to opposite sides of a meadow
while a couple of deer avoid them.

Additionally, we show the navigability of agents in larger constrained
environments in 17(a) and 17(b). Agents are able to navigate the mazes while
avoiding collisions with each other at the bottlenecks.

Figure 16: (0) No interactions. (1) Mass (analogous to stubbornness) weighted
asymmetric interactions. (2) Size and mass weighted interactions. (3)
Symmetric interactions.(4) Grouping friends together. (5) Size based
interactions.

(a) Corn maze

(b) Circle maze

Figure 17: (a) large cornmaze with eight agents all trying to get to different
locations through the maze. (b) a circular shaped maze where agents need to
navigate while avoiding collisions.

In Figure 18 we show the flexibility and controllability of our method. We
create a scenario where one agent must meet up with another agent to deliver a
message and then arrive at his end goal faster than the other agents. In
another scenario, we direct the agent to meet up with the second agent, yet
arrive at his end goal at the same time as the other agents. Lastly, we also
provide the simple scene where agents are ignorant of each other. The flash
mob example in Figure 19 shows how our method can be used to position agents
intricately. The corresponding submission video shows how we are able to
control the order of the placement of the agents in the flash mob as well
through the preferred end time parameter. Our method allows careful control
with no manual effort on the part of the user, thus making it useful for games
or other industrial applications.

Figure 18: Left: An agent delivers a message to another agent and rushes to
his goal position. Middle: Two agents meet up briefly before walking to their
separate destinations. Right: No contact between any agents.

(a) Pre-smile flash mob

(b) Post-smile

Figure 19: (a) group of people in a flash mob instructed to form a smily face.
Agents follow unintuitive paths through space and time to maintain the
fluidity of their motion. (b) a smily face created by controlling the motion
of the mob.

And finally, for the sake of completeness we show Figure 15, a large scene
with 104 agents. We include this example to demonstrate that even though our
method is designed for small to medium sized scenes, it can work on simple
scenes with larger numbers of agents with asymmetric interactions and still
scales linearly in time.

###  4.1 Timings

Figure 11 provides three plots, each with different scaling information for
different usage scenarios. The first plot measure the average time per solver
function iteration for an increasing number of agents. Each agent trajectory
has 100 nodes (300 DOFs) and the performance is linear with broad-phase
collision detection enabled. Our second plot measure the average time per
iteration for an increasing number of DOFS for an 8 agent circle. The number
of nodes starts off at 30 and goes all the way to 180 nodes per agent.
Performance is very nearly linear with broad phase collision detection
enabled. Without a broad phase, performance is quadratic. Something to note is
that increasing the number of DOFs per agent provides no additional benefit to
the quality of the simulation since the actual physics of the space-time rods
are not the end result of our method. As long as agents trajectory curves have
enough DOFs to traverse the environment smoothly, the end results will be
good. Our third plot measures the time per iteration for an increasing
complexity in the environment mesh on the corn maze ( 17(a)) example with
three agents and 200 trajectory nodes per agent. Again, performances is nearly
linear.

As mentioned before there is no gold-standard crowd simulation algorithm since
each scenario is so unique and intricate. Rather than focus on timing
performance for large crowds, our method focuses on solving previously
overlooked path planning scenarios with asymmetric interactions for small to
medium sized groups in highly constrained environments. Our asymptotics show
linear to near-linear performance, but there is plenty of room for improvement
in wall-clock-times through optimization and parallelization.

##  5 Conclusion and Future Work

In this paper we show that modeling the motion of agents through space and
time using a 3D curve resolves a number of difficulties with multi-agent path
planning. Assigning physical characteristics such as a radius and mass to the
agent’s trajectory curve lets us intuitively simulate asymmetric interactions
between agents. Our method outputs agent paths that are smooth, can navigate
through highly constrained environments, and are parameterized by intuitive
controls.

In the future we hope to improve our agent model to include limited
environmental perception (rather than the omniscience our agents enjoy
currently) as well as additional dynamics. As mentioned before, performance
depends on two factors: number of agents and their temporal support. While we
can currently push along these axes independently, a future work is to be able
to push along both axes together, i.e. handle complex environments with large
numbers of agents with intricate motion. We also hope to extend the method in
order to handle scenes with conflicting and changing goals, e.g. a predator-
prey scenario. We would also like to extend our space-time approach to fully
3D environments which requires performing our optimization in four dimensions.
Finally, while our focus is on robustness in path planning, there is still
further room for improvement in performance. We are excited to explore fast
space-time multi-agent path planning to scale our approach to the large dense
crowds that continuum approaches excel at, while maintaining our unique
advantages.

## References

  * [CC14] Panayiotis Charalambous and Yiorgos Chrysanthou  “The pag crowd: A graph based approach for efficient data-driven crowd simulation”  In _Computer Graphics Forum_ 33.8, 2014, pp. 95–108  Wiley Online Library
  * [DKG20] Bobby Davis, Ioannis Karamouzas and Stephen Guy  “NH-TTC: A gradient-based framework for generalized anticipatory collision avoidance”  In _Robotics Science and Systems 2020_ , 2020  DOI: 10.15607/RSS.2020.XVI.078
  * [DMC*17] Teófilo Bezerra Dutra et al.  “Gradient-based steering for vision-based crowd simulation algorithms”  In _Computer graphics forum_ 36.2, 2017, pp. 337–348  Wiley Online Library
  * [GKLM11] Stephen J Guy, Sujeong Kim, Ming C Lin and Dinesh Manocha  “Simulating heterogeneous crowd behaviors using personality trait theory”  In _Proceedings of the 2011 ACM SIGGRAPH/Eurographics symposium on computer animation_ , 2011, pp. 43–52
  * [GNCL14] A. Golas, R. Narain, S. Curtis and M.. Lin  “Hybrid Long-Range Collision Avoidance for Crowd Simulation”  In _IEEE Transactions on Visualization and Computer Graphics_ 20.7, 2014, pp. 1022–1034  DOI: 10.1109/TVCG.2013.235
  * [GVL*12] Stephen J Guy et al.  “A statistical similarity measure for aggregate crowd dynamics”  In _ACM Transactions on Graphics (TOG)_ 31.6 ACM New York, NY, USA, 2012, pp. 1–11
  * [HKHL13] Kyunglyul Hyun, Manmyung Kim, Youngseok Hwang and Jehee Lee  “Tiling Motion Patches”  In _IEEE Transactions on Visualization and Computer Graphics_ 19.11, 2013, pp. 1923–1934  DOI: 10.1109/TVCG.2013.80
  * [HM98] Dirk Helbing and Peter Molnar  “Social Force Model for Pedestrian Dynamics”  In _Physical Review E_ 51, 1998  DOI: 10.1103/PhysRevE.51.4282
  * [Hug02] Roger Hughes  “A continuum theory for the flow of pedestrians”  In _Transportation Research Part B: Methodological_ 36, 2002, pp. 507–535  DOI: 10.1016/S0191-2615(01)00015-7
  * [KSNG17] Ioannis Karamouzas, Nick Sohre, Rahul Narain and Stephen J. Guy  “Implicit Crowds: Optimization Integrator for Robust Crowd Simulation”  In _ACM Transactions on Graphics_ 36.4, 2017  DOI: 10.1145/3072959.3073705
  * [LLKP11] Sergey Levine, Yongjoon Lee, Vladlen Koltun and Zoran Popovic  “Space-Time Planning with Parameterized Locomotion Controllers”  In _ACM Transactions on Graphics_ 30, 2011  DOI: 10.1145/1966394.1966402
  * [NGCL09] Rahul Narain, Abhinav Golas, Sean Curtis and Ming C Lin  “Aggregate dynamics for dense crowd simulation”  In _ACM SIGGRAPH Asia 2009 papers_ , 2009, pp. 1–8
  * [NW06] Jorge Nocedal and Stephen J. Wright  “Numerical Optimization”  New York, NY, USA: Springer, 2006
  * [OPOD10] Jan Ondřej, Julien Pettré, Anne-Hélène Olivier and Stéphane Donikian  “A synthetic-vision based steering approach for crowd simulation”  In _ACM Transactions on Graphics (TOG)_ 29.4 ACM New York, NY, USA, 2010, pp. 1–9
  * [POSB05] Núria Pelechano Gómez, Kevin O’Brien, Barry G Silverman and Norman Badler  “Crowd simulation incorporating agent psychological models, roles and communication”  In _First International Workshop on Crowd Simulation_ , 2005
  * [PSE*00] Jovan Popović et al.  “Interactive manipulation of rigid body simulations”  In _Proceedings of the 27th annual conference on Computer graphics and interactive techniques_ , 2000, pp. 209–217
  * [SKF08] Ari Shapiro, Marcelo Kallmann and Petros Faloutsos  “Interactive motion correction and object manipulation”  In _ACM SIGGRAPH 2008 classes_ , 2008, pp. 1–8
  * [SKH*11] Shawn Singh et al.  “A modular framework for adaptive agent-based steering”  In _Symposium on interactive 3D graphics and games_ , 2011, pp. 141–150
  * [TCP06] Adrien Treuille, Seth Cooper and Zoran Popović  “Continuum crowds”  In _ACM Transactions on Graphics (TOG)_ 25.3 ACM New York, NY, USA, 2006, pp. 1160–1168
  * [TYK*09] Shigeo Takahashi et al.  “Spectral-based group formation control”  In _Computer Graphics Forum_ 28.2, 2009, pp. 639–648  Wiley Online Library
  * [VGLM11] Jur Van Den Berg, Stephen J Guy, Ming Lin and Dinesh Manocha  “Reciprocal n-body collision avoidance”  In _Robotics research_ Springer, 2011, pp. 3–19
  * [WJO*14] David Wolinski et al.  “Parameter estimation and comparative evaluation of crowd simulations”  In _Computer Graphics Forum_ 33.2, 2014, pp. 303–312  Wiley Online Library
  * [WK88] Andrew Witkin and Michael Kass  “Spacetime constraints”  In _ACM Siggraph Computer Graphics_ 22.4 ACM New York, NY, USA, 1988, pp. 159–168
  * [WLJT17] Tomer Weiss, Alan Litteneker, Chenfanfu Jiang and Demetri Terzopoulos  “Position-based multi-agent dynamics for real-time crowd simulation”  In _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ , 2017, pp. 1–2
  * [WLP16] David Wolinski, Ming C Lin and Julien Pettré  “WarpDriver: context-aware probabilistic motion prediction for crowd simulation”  In _ACM Transactions on Graphics (TOG)_ 35.6 ACM New York, NY, USA, 2016, pp. 1–11
  * [YSC21] Chris Yu, Henrik Schumacher and Keenan Crane  “Repulsive Curves”  In _ACM Trans. Graph._ 40.2 New York, NY, USA: ACM, 2021