# Two Skyrmion Dynamics with 𝜔 Mesons

## I Introduction

In the large $`N_C`$ or classical limit of QCD, nucleons may be identified as classical solitons of the scalar-isovector $`SU(N_f)`$ pion field. The simplest theory that manifests these solitons is the non-linear sigma model. However the solitons of this theory are not stable against collapse. The first attempt to provide a stable theory was by Skyrme (long before QCD)  who introduced a fourth order term. This term does indeed lead to stabilized solitons that are called skyrmions, and there is a vast body of work on their properties and on how to quantize them . Unfortunately the fourth order term introduces numerical instabilities that make complex dynamical calculation nearly impossible . It is possible to stabilize the non-linear sigma model without the fourth order, or Skyrme term as it is called, by coupling the baryon current to the $`\omega `$ meson field . This not only provides stability, but does so in the context of more reasonable physics. It is also possible to stabilize the solitons by introducing a $`\rho `$ meson field, with a gauge coupling, or with both the $`\omega `$ and $`\rho `$, but the $`\rho `$ adds a great deal to the numerical complications. Hence in this first, exploratory work, we stabilize with the $`\omega `$ only. We continue to refer to the solitons as skyrmions.

The Skyrme approach, either in its fourth order or $`\omega `$ stabilized form, has much to recommend it as a model of low energy strong interaction physics. This low energy or long wave length domain is notoriously difficult to describe in the context of standard QCD. The Skyrme approach includes chiral symmetry, baryons (and therefore also antibaryons), pions, the one pion exchange potential, and, with quantization, nucleons and deltas. The idea of having nucleons arise naturally from an effective theory of pions and vector mesons is especially attractive as a path from QCD to such effective Lagrangians is better understood.

It has been shown that the Skyrme model can give a good account, with very few parameters, of the low energy properties of nucleons . In the baryon number two system, the Skyrme approach can describe the principal features of the nucleon-nucleon static potential . Most of these problems have also been successfully studied in more traditional nuclear physics forms.

One problem that has not yielded easily to traditional physics approaches and is out of the reach of standard lattice QCD is low energy annihilation. This is a problem ideally suited to the Skyrme approach. We have already shown that a general picture of post-annihilation dynamics including branching ratios can be obtained from the Skyrme approach . The nucleon-antinucleon potential can be as well . What remains is a full modeling of annihilation from start to finish. The only attempts we know of to do that in the Skyrme model had numerical problems associated with the Skyrme term  . We propose to study skyrmion dynamics using $`\omega `$ stabilization, thus avoiding the usual numerical problems. As a prelude to studying annihilation, we have studied scattering in the baryon number two system. It is that work we report here.

In the Skyrme model the three components of the pion field, $`\pi _1,\pi _2,\pi _3`$, can be aligned with the spatial directions, $`x,y,z`$, providing a correspondence between the two spaces. This simple alignment is called a hedgehog. A rotation of the pion field with respect to space is called a grooming. The energy of a single free skyrmion is independent of grooming, but the interaction between two skyrmions depends critically, as we shall see, on their relative grooming. In this paper we present the results of calculations of skyrmion-skyrmion scattering (with $`\omega `$ meson stabilization) at medium energy for a variety of impact parameters and groomings. We find rich structure. Some of the channels have simple scattering, but some display radiation, scattering out of the scattering plane, orbiting and ultimately capture to a bound $`B=2`$ state. These calculations are numerically complex and require considerable computational resources and computing time, but they are numerically stable. We know of no other calculations of $`B=2`$ skyrmion scattering that show the phenomena we find. Furthermore our success here bodes well for extending the method to annihilation.

In the next section we briefly review the model we use, presenting both the Lagrangian for the $`\omega `$ stabilized non-linear sigma model, and our equations of motion. In Section 3 we discuss our numerical strategies and methods. Section 4 presents our results, mostly in graphical form, and Section 5 deals with conclusions and outlook. The reader interested only in results can go directly to Sections 4 and 5.

## II Formalism

### A Model

Our starting point is the non-linear $`\sigma `$ model Lagrangian,

$$_\sigma =\frac{1}{4}f_\pi ^2\mathrm{t}r\left(_\mu 𝒰^\mu 𝒰^{}\right)+\frac{1}{2}m_\pi ^2f_\pi ^2\mathrm{t}r\left(𝒰1\right),$$
(1)   
where the $`SU(2)`$ field $`𝒰`$ is parameterized by the three real pion fields<sup>*</sup><sup>*</sup>* To avoid confusion, throughout this paper we reserve the plain letter $`\pi `$ for the pion field, and use the symbol $`\underset{¯}{\pi }=3.1415\mathrm{}`$ for the mathematical constant. $`\{\pi _k\}_{k=1,3}=\stackrel{}{\pi }`$,

$$𝒰=\mathrm{exp}\left(i(\stackrel{}{\tau }\stackrel{}{\pi })\right)=\mathrm{cos}\pi +i\mathrm{sin}\pi (\stackrel{}{\tau }\widehat{\pi }).$$
(2)   
Here, $`\{\tau _k\}_{k=1,3}=\stackrel{}{\tau }`$ are the Pauli matrices (in flavor space). We identify the baryon current with

$$B^\mu =\frac{1}{24\underset{¯}{\pi }^2}ϵ^{\mu \nu \alpha \beta }\mathrm{t}r\left[\left(𝒰^{}_\nu 𝒰\right)\left(𝒰^{}_\alpha 𝒰\right)\left(𝒰^{}_\beta 𝒰\right)\right].$$
(3)   
The full Lagrangian also contains the free $`\omega `$ and the interaction part,

$$=_\sigma +\frac{3}{2}g\omega _\mu B^\mu \frac{1}{2}_\mu \omega _\nu \left(^\mu \omega ^\nu ^\nu \omega ^\mu \right)+\frac{1}{2}m^2\omega _\mu \omega ^\mu .$$
(4)   
We take $`f_\pi =62\mathrm{MeV}`$ following  and $`g=\frac{m}{\sqrt{2}f_\pi }`$ from VMD where $`m=770\mathrm{MeV}`$ is the vector meson mass. This gives a skyrmion mass of $`791\mathrm{MeV}`$.

### B Choice of Dynamical Variables

The traditional approach to numerical simulations involving skyrmions uses the cartesian decomposition of the unitary matrix $`𝒰`$,

$`𝒰=\mathrm{\Phi }𝐈_2+\stackrel{}{\mathrm{\Psi }}\stackrel{}{\tau };\mathrm{\Phi }=\mathrm{cos}\pi ,\stackrel{}{\mathrm{\Psi }}=\widehat{\pi }\mathrm{sin}\pi `$ (5)   
The quadruplet of real numbers $`(\mathrm{\Phi },\stackrel{}{\mathrm{\Psi }})`$ is constrained by the unitarity condition $`𝒰𝒰^{}=𝐈_2\mathrm{\Phi }^2+\stackrel{}{\mathrm{\Psi }}\stackrel{}{\mathrm{\Psi }}=1`$, also known as the chiral condition. The Lagrangian is usually written in terms of $`(\mathrm{\Phi },\stackrel{}{\mathrm{\Psi }})`$ and the chiral condition is imposed using a Lagrange multiplier field. The four coordinates of $`𝒰`$ are similar to the cartesian coordinates of a point in $`𝐑^4`$ confined to the surface of a hypersphere of unit radius.

An attractive idea is to use an approach which ensures naturally the unitarity of $`𝒰`$ at all times. Such a method is used successfully in lattice QCD in the context of the hybrid Monte-Carlo algorithm . There, the dynamical variables are $`𝒰`$ itself and the Hermitian matrix $`\dot{𝒰}𝒰^{}`$. In the present work, we employ the parameterization of $`𝒰`$ that follows from (2). We will use $`\pi `$ and $`\widehat{\pi }`$ and their time derivatives $`\dot{\pi },\dot{\widehat{\pi }}`$ as dynamical variables. This implements only in part the principle mentioned above, since the unit vector $`\widehat{\pi }`$ is subject to conditions similar to the four-dimensional vector $`(\mathrm{\Phi },\stackrel{}{\mathrm{\Psi }})`$. However, one has better geometrical intuition for the former. The connection to $`(\mathrm{\Phi },\stackrel{}{\mathrm{\Psi }})`$ is straightforward via (5). The Lagrangian in terms of $`\pi ,\widehat{\pi }`$ is

$``$ $`=`$ $`_\sigma +_{int}+_\omega `$ (6)   
$`_\sigma `$ $`=`$ $`{\displaystyle \frac{1}{2}}f_\pi ^2\left(_\mu \pi ^\mu \pi +\mathrm{sin}^2\pi _\mu \widehat{\pi }^\mu \widehat{\pi }\right)+m_\pi ^2f_\pi ^2(\mathrm{cos}\pi 1)`$ (7)   
$`_{int}`$ $`=`$ $`{\displaystyle \frac{3g}{8\underset{¯}{\pi }^2}}ϵ^{\mu \nu \alpha \beta }\mathrm{sin}^2\pi \omega _\mu _\nu \pi \left[\widehat{\pi }\left(_\alpha \widehat{\pi }\times _\beta \widehat{\pi }\right)\right]`$ (8)   
$`_\omega `$ $`=`$ $`{\displaystyle \frac{1}{2}}_\mu \omega _\nu \left(^\mu \omega ^\nu ^\nu \omega ^\mu \right)+{\displaystyle \frac{1}{2}}m_{vec}^2\omega _\mu \omega ^\mu .`$ (9)   
### C Equations of Motion

We wish to obtain the Euler-Lagrange equations for $`\pi ,\widehat{\pi }`$. The problem is that $`\widehat{\pi }`$ is a unit vector, $`\widehat{\pi }\widehat{\pi }=1`$, so its components cannot be treated as true coordinates. One may still write down the Euler-Lagrange equations, by considering the action

$`S={\displaystyle d^3x𝑑t(\pi (x,t),\widehat{\pi }(x,t),\omega _\alpha (x,t);_\mu \pi (x,t),_\mu \widehat{\pi }(x,t),_\mu \omega _\alpha (x,t))}.`$ (10)   
The Lagrange equation for $`\pi `$ is obtained by requiring that $`S`$ be stationary with respect to a small variation $`\delta \pi (x,t)`$ and the corresponding variations $`\delta (_\mu \pi )=_\mu (\delta \pi )`$. The variation of $`S`$ which has to vanish for any $`\delta \pi `$, is

$`\delta S={\displaystyle d^3x𝑑t\left\{\frac{}{\pi }\delta \pi +\frac{}{(_\mu \pi )}\delta (_\mu \pi )\right\}}={\displaystyle d^3x𝑑t\left\{\frac{}{\pi }_\mu \left[\frac{}{(_\mu \pi )}\right]\right\}\delta \pi }.`$ (11)   
Therefore the quantity in curly brackets in the last equality of (11) has to be identically zero, since $`\delta \pi `$ is completely arbitrary. Of course this leads exactly to the usual form of the Euler-Lagrange equations. We may repeat formally the same steps for $`\widehat{\pi }`$, leading to

$`\delta S={\displaystyle d^3x𝑑t\left\{\frac{}{\widehat{\pi }}_\mu \left[\frac{}{(_\mu \widehat{\pi })}\right]\right\}\delta \widehat{\pi }}.`$ (12)   
Here, partial differentiation with respect to the vector $`\widehat{\pi }`$ or its derivatives gives the vector obtained by differentiating with respect to the corresponding components of $`\widehat{\pi }`$ or its derivatives. If we want $`\delta S`$ to vanish for any $`\delta \widehat{\pi }`$, we do not need to require that the quantity in curly brackets vanish identically. This is because $`\delta \widehat{\pi }`$ is not completely arbitrary. Both $`\widehat{\pi }`$ and $`\widehat{\pi }^{}=\widehat{\pi }+\delta \widehat{\pi }`$ have to be unit vectors. The variation of $`\widehat{\pi }\widehat{\pi }`$ must vanish,

$`0=\delta \left(\widehat{\pi }\widehat{\pi }\right)=2\delta \widehat{\pi }\widehat{\pi }.`$ (13)   
In other words, $`\delta \widehat{\pi }(x,t)\widehat{\pi }(x,t)`$ for any $`(x,t)`$. Therefore, the necessary condition for the stationarity of $`S`$ is simply

$`\left\{{\displaystyle \frac{}{\widehat{\pi }}}_\mu \left[{\displaystyle \frac{}{(_\mu \widehat{\pi })}}\right]\right\}\widehat{\pi }.`$ (14)   
The above statement is easily turned into an equation by subtracting the projection of the left hand side onto $`\widehat{\pi }`$,

$`\left\{{\displaystyle \frac{}{\widehat{\pi }}}_\mu \left[{\displaystyle \frac{}{(_\mu \widehat{\pi })}}\right]\right\}\widehat{\pi }\left\{{\displaystyle \frac{}{\widehat{\pi }}}_\mu \left[{\displaystyle \frac{}{(_\mu \widehat{\pi })}}\right]\right\}=0.`$ (15)   
The Euler-Lagrange equations we obtain finally are

$`_\nu ^\nu \pi `$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}2\pi _\mu \widehat{\pi }^\mu \widehat{\pi }m_\pi ^2\mathrm{sin}\pi +{\displaystyle \frac{3g}{8\underset{¯}{\pi }^2f_\pi ^2}}ϵ^{\mu \nu \alpha \beta }\mathrm{sin}^2\pi _\mu \omega _\nu \left[\widehat{\pi }\left(_\alpha \widehat{\pi }\times _\beta \widehat{\pi }\right)\right]`$ (16)   
$`_\mu ^\mu \widehat{\pi }`$ $`=`$ $`\widehat{\pi }\left(\widehat{\pi }_\mu ^\mu \widehat{\pi }\right){\displaystyle \frac{2\mathrm{cos}\pi }{\mathrm{sin}\pi }}_\mu \pi ^\mu \widehat{\pi }{\displaystyle \frac{3g}{4\underset{¯}{\pi }^2f_\pi ^2}}ϵ^{\mu \nu \alpha \beta }_\mu \omega _\nu _\alpha \pi \left(_\beta \widehat{\pi }\times \widehat{\pi }\right)`$ (17)   
$`_\nu ^\nu \omega ^\mu `$ $`=`$ $`m_{vec}^2\omega ^\mu {\displaystyle \frac{3g}{8\underset{¯}{\pi }^2}}ϵ^{\mu \nu \alpha \beta }\mathrm{sin}^2\pi _\nu \pi \left[\widehat{\pi }\left(_\alpha \widehat{\pi }\times _\beta \widehat{\pi }\right)\right].`$ (18)   
## III Computation

Our calculation is based on the equations of motion (16), using $`(\pi ,\widehat{\pi },\omega _\mu ;\dot{\pi },\dot{\widehat{\pi }},\dot{\omega }_\mu )`$ as variables. This choice leads to two problems which have to be tackled by our discretization scheme. First, we have a coordinate singularity at $`\pi =0`$ and $`\pi =\underset{¯}{\pi }`$. Here, $`\widehat{\pi }`$ is not defined. The situation is similar to that of angles in polar coordinates when the radius vanishes. While the equations of motion are correct for any non-zero $`\pi `$, in the vicinity of the coordinate singularity small changes of $`\stackrel{}{\pi }`$ are translated into very large variations of $`\widehat{\pi }`$. A scheme based on the equations of motion (16) breaks down around the ‘poles’ of the hypersphere due to the large discretization errors involved.

One way out is to rotate the coordinate system in $`SU(2)`$ so as to avoid the problem regions. We introduce a new field $`𝒱`$ obtained by rotation with a fixed $`U_0SU(2)`$,

$`𝒰=U_0𝒱;𝒱=\mathrm{exp}\left(i\stackrel{}{\tau }\stackrel{}{\sigma }\right);U_0=\mathrm{exp}\left(i\stackrel{}{\tau }\stackrel{}{\mathrm{\Theta }}\right).`$ (19)   
Substituting (19) into our Lagrangian (1)(4), we find that only the pion mass term is modified since everywhere else $`𝒰`$ is combined with $`𝒰^{}`$ or its derivatives so $`U_0`$ drops out. The equations of motion can then be derived in an identical fashion to (16). We cite them here for the sake of completeness,

$`_\nu ^\nu \sigma `$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}2\pi _\mu \widehat{\sigma }^\mu \widehat{\sigma }m_\pi ^2\left(\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\sigma +(\widehat{\mathrm{\Theta }}\widehat{\sigma })\mathrm{sin}\mathrm{\Theta }\mathrm{cos}\sigma \right)+{\displaystyle \frac{3g}{8\underset{¯}{\pi }^2f_\pi ^2}}ϵ^{\mu \nu \alpha \beta }\mathrm{sin}^2\sigma _\mu \omega _\nu \left[\widehat{\sigma }\left(_\alpha \widehat{\sigma }\times _\beta \widehat{\sigma }\right)\right]`$ (20)   
$`_\mu ^\mu \widehat{\sigma }`$ $`=`$ $`\widehat{\sigma }\left(\widehat{\sigma }_\mu ^\mu \widehat{\sigma }\right){\displaystyle \frac{2\mathrm{cos}\sigma }{\mathrm{sin}\sigma }}_\mu \sigma ^\mu \widehat{\sigma }{\displaystyle \frac{3g}{4\underset{¯}{\pi }^2f_\pi ^2}}ϵ^{\mu \nu \alpha \beta }_\mu \omega _\nu _\alpha \sigma \left(_\beta \widehat{\sigma }\times \widehat{\sigma }\right)m_\pi ^2\left(\widehat{\mathrm{\Theta }}\widehat{\sigma }(\widehat{\mathrm{\Theta }}\widehat{\sigma })\right)\mathrm{sin}\mathrm{\Theta }\mathrm{sin}\sigma `$ (21)   
$`_\nu ^\nu \omega ^\mu `$ $`=`$ $`m_{vec}^2\omega ^\mu {\displaystyle \frac{3g}{8\underset{¯}{\pi }^2}}ϵ^{\mu \nu \alpha \beta }\mathrm{sin}^2\sigma _\nu \sigma \left[\widehat{\sigma }\left(_\alpha \widehat{\sigma }\times _\beta \widehat{\sigma }\right)\right].`$ (22)   
In our calculation, we first rotate the field to be updated, together with the surrounding fields, so that the corresponding $`(\sigma ,\widehat{\sigma })`$ is comfortably away from the coordinate singularities. Then we apply the discretized equations of motion (see below) derived from (20) and finally we rotate back the updated fields. We stress that switching to rotated fields and (20) amounts to a mere changing of the coordinate system. The content of the equations is identical irrespective of the coordinate system. The difference is in the discretization which in the vicinity of the ‘poles’ leads to large errors which are avoided in the rotated frame. The choice of $`\stackrel{}{\mathrm{\Theta }}`$ is largely arbitrary. For simplicity we choose it so that $`\stackrel{}{\sigma }`$ is always on the ‘equator’ of the $`SU(2)`$ hypersphere.

The second important numerical issue follows from the unit vector nature We will continue to use $`\stackrel{}{\pi }`$ when referring to the pion field, even though we perform our updates using the rotated pion fields $`\stackrel{}{\sigma }`$. of $`\widehat{\pi }`$, $`\stackrel{}{\pi }\stackrel{}{\pi }=1`$. This and the corresponding constraints on $`\dot{\widehat{\pi }}`$ and $`\ddot{\widehat{\pi }}`$, namely $`\widehat{\pi }\dot{\widehat{\pi }}=0`$, $`\ddot{\widehat{\pi }}\widehat{\pi }+\dot{\widehat{\pi }}\dot{\widehat{\pi }}=0`$ are consistent with the equations of motion but are violated by the discretized equations. Therefore they need to be imposed by projecting out the spurious components of $`\dot{\widehat{\pi }}`$ and $`\ddot{\widehat{\pi }}`$ at every step in the update. There are similar conditions for the spatial derivatives, which also have to be taken into account. This is of course a rather technical point, but it is a reminder of the fact that our approach does not completely conform to the idea of using a minimal set of coordinates. That would be achieved for instance by trading $`\widehat{\pi }`$, which describes a direction in three dimensions, for the two angles that define that direction. That would have introduced another set of coordinate singularities to be avoided. The opposite strategy would be to apply the same considerations we followed for $`\widehat{\pi }`$ in deriving the Euler-Lagrange equations and imposing the unit length constraint, to the four-dimensional unit vector $`(\mathrm{\Phi },\stackrel{}{\mathrm{\Psi }})`$. That also remains an open possibility, along with using unitary and Hermitian matrices as dynamical variables and of course, the traditional path employing Lagrange multipliers. We have not found it necessary to employ any of these notions at this time. Besides the fact that it leads to a reasonably stable calculation, our choice of variables has the advantage of fairly simple equations of motion. The decomposition of $`\stackrel{}{\pi }`$ into its length and unit vector follows quite naturally, and leads to great simplification in the much more complicated case of including a vector-isovector $`\rho `$ meson field.

The main idea of our numerical scheme is the following. The discretized time evolution for the fields themselves is quite straightforward, given the knowledge of their time-derivatives or ‘velocities’. The evolution of the velocities is the core of our procedure. The velocities are assigned to half-integer timesteps. Thus, one timeslice contains the fields at a given time and the velocities half a timestep earlier. The time evolution of the velocities follows from solving the equations of motion for the second time derivatives. The latter are written in terms of the retarded velocity (at time $`(t\mathrm{\Delta }t/2)`$) and the one at $`(t+\mathrm{\Delta }t/2)`$. The rest of the equations of motion contain the local fields and their spatial derivatives, (all defined at time $`t`$), but they also contain the velocities as it is evident from (20). We approximate the velocities at $`t`$ by the average of their values at $`(t\pm \mathrm{\Delta }t/2)`$. This leads to an implicit equation for the updated velocities. We solve the implicit equation iteratively to cubic order in $`\mathrm{\Delta }t`$ which exceeds the order of accuracy in which it was derived. Finally, the new velocities are used to update the fields.

Use of a flux-conservative form for the equations of motion might have lead to improved accuracy. However, this way we can use the known conservation laws, in particular, energy and baryon number conservation, as a check for the accuracy of our simulations. The long-term stability of the calculation is largely improved by the feedback mechanism built into our implicit scheme. Calculations performed using an explicit scheme involving a second timeslice shifted with one half timestep give virtually identical results as far as the time evolution of the fields and energy and baryon number conservation are concerned. Part of the calculations presented below – those not involving complicated physical situations – have been in fact performed with this earlier version of our code. When pushed beyond $`25\mathrm{fm}/\mathrm{c}`$ in time (a few thousand timesteps), these calculations typically develop instabilities which build up fast to destructive levels. By contrast, the feedback calculations ran without exception up to 5000 timesteps or more without becoming unstable. This is achieved without adding an explicit numerical viscosity term to our equations.

The calculation based on the algorithm described above was implemented on a three-dimensional grid of points. The physical size of the box used for the calculations shown in this paper was $`18\times 10\times 10`$ fermi. We take advantage of the four-fold spatial symmetry of the problem , therefore our mesh covers only one quadrant of the physical box. We used a fixed lattice spacing of $`0.1\mathrm{fm}`$ therefore our lattice has $`91\times 101\times 51`$ regular points. In addition to these, we have a layer of unequally spaced points on each outside wall of the box simulating an absorbing boundary. Thus our full mesh has $`101\times 121\times 61`$ points. An indication of the better intrinsic stability of the model Lagrangian employed here is the fact that our calculations are fairly stable out to $`50\mathrm{fm}/\mathrm{c}`$ in time with timesteps of $`0.1`$ to $`0.4\mathrm{fm}/\mathrm{c}`$, corresponding to a CFL ratio of $`0.1`$ to $`0.4`$. This is larger than in the early works of Verbaarschot et al  ($`0.05`$) and of the Caltech group  ($`0.075`$ to $`0.013`$), and comparable to the more recent calculations of , which employ fourth-order spatial differences, in contrast with our second-order spatial difference scheme. In the absence of radiation, the total energy is conserved to within $`3\%`$ for typically $`20\mathrm{fm}`$. Due to the emission of radiation which eventually leaves the box it is harder to assess the degree of energy conservation for those processes which involve (quasi)bound states and have to be followed for a longer time. We can estimate the amount of energy being radiated by calculating the energy contained in a sphere of given radius ($`2\mathrm{fm}`$) around the origin, large enough to contain most of the field. The radiated energy first leaves this sphere and then the box. Excluding loss through radiation identified in this manner, the total energy is conserved to within $`4\%`$ in all the calculations presented below, including the long runs ($`40\mathrm{}60\mathrm{fm}/\mathrm{c}`$) involving bound states. A check for consistency is to reverse the arrow of time at some point in the calculation, which should lead back close to the initial state. We performed this check successfully on one case with nontrivial dynamics but little radiation.

We construct the initial state by first numerically solving the field equations in the spherically symmetric ansatz. This gives us the radial functions for the spherically symmetric static skyrmion (hedgehog). We then place a boosted hedgehog configuration in the simulation box. The skyrmion is boosted towards its mirror image implemented via the boundary conditions at the symmetry walls of the box, which corresponds only to one quadrant of the physical region being simulated. One can dial all the relative groomings discussed below as well as the corresponding skyrmion-antiskyrmion configurations by simply choosing the appropriate signs in the boundary conditions for the various field components . The problem of the overlap of the tails of the two skyrmions in the initial state is not easy to solve. Instead, we chose the initial configuration with a large separation ($`9\mathrm{fm}`$) between the two centers, so that the overlap becomes truly negligible.

Our calculations have been performed on clusters of $`612`$ IBM SP-2 computers. Depending on the number of processors and the timestep, one $`20\mathrm{fm}/\mathrm{c}`$ calculation takes half a day to two days to complete.

## IV Results

For skyrmion-skyrmion scattering with finite impact parameter, there are four relative groomings that are distinct. The first is no grooming at all. This is the hedgehog-hedgehog channel. The second grooming we consider is a relative rotation of $`\underset{¯}{\pi }`$ about an axis perpendicular to the scattering plane (the plane formed by the incident direction and the impact parameter). The third is a relative grooming by $`\underset{¯}{\pi }`$ about the direction of motion. For zero impact parameter this channel is known to be repulsive.  The fourth grooming consists of a relative rotation by $`\underset{¯}{\pi }`$ about the direction of the impact parameter. In the limit of zero impact parameter the second and fourth groomings become equivalent, they then correspond to a rotation of $`\underset{¯}{\pi }`$ about an axis normal to the incident direction. In this zero impact parameter case, this channel is known to be attractive.  All of our scattering studies are done at a relative velocity of $`\beta =v/c=0.5`$, which corresponds to a center of mass kinetic energy of 230 MeV. In order to keep our study relatively small, we have not studied the effect of varying the incident energy. As the main means of presentation we chose energy density contour plots. The baryon density plots are very similar to the energy density plots.

### A Head-on collisions

Let us begin, for simplicity, with the case of zero impact parameter. For the hedgehog-hedgehog channel (HH) and the repulsive channel (rotation by $`\underset{¯}{\pi }`$ about the incident direction) symmetry dictates that the scattering can only be exactly forward (the skyrmions passing through on another) or exactly backward (the skyrmions bouncing back off each other). We find that in fact the scattering is backward in both the HH and the repulsive channels.

We illustrate this type of scattering in Figure 1 with contour plots of the energy density in the $`xy`$ plane Throughout the paper we adopt the following convention. The direction of the initial motion is $`x`$, the impact parameter – if nonzero – points in the $`y`$ direction and $`z`$ is the direction perpendicular to the \[initial\] plane of motion, $`xy`$. For zero impact parameter the choice of $`y`$ and $`z`$ is of course arbitrary. for the repulsive channel. Unless otherwise specified we keep the same choice for the energy contour levels for all similar plots that will follow, namely, the first contour at $`5\mathrm{M}\mathrm{e}\mathrm{V}/\mathrm{fm}^3`$, and the others equally spaced at $`100\mathrm{M}\mathrm{e}\mathrm{V}/\mathrm{fm}^3`$. The head-on scattering in the repulsive channel has nothing surprising. The process is reminiscent of two tennis balls bouncing off each other. The skyrmions start compressing as soon as they touch. They slow down, compress and stop, and then expand and move off in the direction opposite to the one they came in along. The collision is practically elastic. We were unable to detect any energy loss through radiation and the velocities of the topological centers before and after collision are practically the same.

For the hedgehog-hedgehog channel we again find backward scattering with an evolution very similar to that shown in Figure 1, hence we show no figure.

In the attractive channel (rotation by $`\underset{¯}{\pi }`$ normal to the incident direction) the skyrmions scatter at $`90^0`$ along an axis perpendicular to the plane formed by the incident direction and the grooming axis. This right angle scattering is well known.  It proceeds through the axially symmetric $`B=2`$ configuration . The skyrmions lose their individual identity in this process.

In Figure 2 we show the energy contours for head-on collision in the attractive channel. At the midpoint of this process, shown in in the fourth frame of Figure 2, one can clearly see the torus-shaped configuration. It is situated in the plane defined by the incoming and outgoing directions and perpendicular to the grooming axis. The bulk of the energy density avoids the center of the doughnut as it shifts from the incoming direction to the one perpendicular to it. To illustrate this point better we plot the middle frames of Figure 2 in perspective in Figure 3.<sup>§</sup><sup>§</sup>§ We remind the reader that we plot level contours of the total energy density in the median plane of our three-dimensional system, as opposed to three-dimensional surfaces of constant energy.

In Figure 4 we plot the coordinates of the center of one skyrmion (defined as the point where the pion field amplitude is exactly $`\underset{¯}{\pi }`$) versus. time. The only nonzero coordinate is $`x`$, initially. Then, after the right-angle scattering, $`y`$ is the only non-zero coordinate. Straight lines indicate uniform motion. That is the case both before and after the collision. The slopes of the $`y`$ line is noticeably smaller than that of the $`x`$ line before collision, in other words, the outgoing velocity is slightly smaller. This is due to a genuine physical process, radiation, rather than to a numerical artifact, because there is no decrease in the total energy of the system. In some of the processes we discuss below, there are stronger examples of slowing down, accompanied by detectable radiation.

As we have seen above, the scattering direction is determined by the incident direction and the grooming direction. When the grooming direction is normal to the incident direction, a torus is formed in the plane normal to the grooming direction. In the presence of a non-zero impact parameter, it matters whether the grooming direction is parallel or normal to that impact parameter. In both cases there is a tendency to the formation of a torus normal to the grooming direction, but the ultimate evolution of the scattering is different in the two cases. We will refer to the case where the grooming is normal to the impact parameter as attractive (1) and attractive (2) when they are parallel.

Let us now look at scattering in each grooming as a function of impact parameter. In each of the four cases we study impact parameters of $`0.4\mathrm{fm}`$, $`0.8\mathrm{fm}`$, $`1.6\mathrm{fm}`$, and $`2.8\mathrm{fm}`$. In all cases we find that the scattering at $`2.8\mathrm{fm}`$ is “routine” and hence we do not go to larger impact parameter.

### B Simple scattering: the hedgehog-hedgehog and attractive (1) channels

We begin with two groomings that lead to simple dynamics. Consider first the hedgehog-hedgehog channel. In Figure 5 we plot the trajectories of the topological center of one of the skyrmions for zero impact parameter and for each of the four non-zero impact parameters. Here and in the other similar plots, we define the center as the point where the norm of the pion field reaches the value of $`|\stackrel{}{\pi }|=\underset{¯}{\pi }`$. This is a good indicator of the global movement of a skyrmion, especially when the two colliding objects are somewhat separated.The topological center does not necessarily coincide with the center of mass. While the latter is insensitive to internal oscillations of the skyrmion, the topological center oscillates. We cut off the trajectories in Figure 5 after $`t=17\mathrm{fm}/\mathrm{c}`$. We see normal scattering trajectories corresponding to a repulsive interaction. This channel is indeed known to be mildly repulsive. For the smallest impact parameter, $`b=0.4\mathrm{fm}`$, the scattering is in the backward direction (recall that for $`b=0`$ the scattering angle is $`180^o`$) and gradually turns forward for increasing impact parameter. This is exactly what one might expect, since this channel is the most similar to point-particle scattering. The interaction between the skyrmions is central in the HH channel, because it is independent of the direction of the relative position vectors of the topological centers.

The large-angle scattering for small impact parameters probes the interaction of the soft core of the skyrmions. We illustrate this in Figure 6 with energy contour plots from the $`b=0.4\mathrm{fm}`$ case. Just like in the head-on case, the skyrmions compress as they touch. They slow down and then accelerate and proceed in the outgoing direction. Internal oscillations of the skyrmions  can be observed after the collision, therefore this process is not entirely elastic.

The other channel that exhibits only simple scattering is the attractive (1) grooming – where the grooming direction is perpendicular to the plane of motion.

We have studied the same five impact parameters here. In Figure 7 we show the trajectory of the skyrmion centers for each of those impact parameters. For $`b=0`$ and this grooming, the scattering angle is $`90^o`$. For large impact parameters the scattering is nearly forward. Hence in between, we should see something in between, which is precisely what the figure shows. Note that in this channel, in contrast to the HH channel, the trajectories begin their scattering by curving toward the scattering center, as one would expect for an attractive channel. The attraction is seen to act practically unperturbed in the scattering with $`b=2.8\mathrm{fm}`$.

Again we illustrate one scattering process in more detail. In Figure 8 we plot energy contours for the $`b=0.8\mathrm{fm}`$ scattering. The most important feature is the existence in the third and fourth frames of a configuration reminiscent of the doughnut-shaped intermediate state identified in the $`b=0`$ case for this grooming (Figures 2 and 3). We are further away from axial symmetry than in Figure 2, but the relatively large void in the middle is clearly identifiable. It is important to point out that even at this mutual distance of $`0.8\mathrm{fm}`$, the attractive interaction is strong enough to start forming the $`B=2`$ configuration. In the present case it is ripped apart by momentum. As we shall see that is not always the case in all groomings. Another feature to mention is the presence of internal oscillations. They can be clearly seen in the last few frames of Figure 8 deforming the outgoing skyrmions. In the case of $`b=1.6\mathrm{fm}`$ trajectory, where the attraction is too weak to lead to a $`B=2`$ configuration, the topological centers are attracted toward each other but then they bounce back and oscillate transversally.

### C Orbiting, capture and radiation: the repulsive channel

For the remaining two groomings, where the grooming direction is in the scattering plane, the dynamics is more complex. Let us first consider grooming about the incident direction, referred to as the repulsive channel. For the smaller impact parameters $`b=0.4\mathrm{fm}`$ and $`b=0.8\mathrm{fm}`$, we find a remarkable type of trajectory. The scattering begins as repulsive as we see in Figure 9, but as the skyrmions pass one another, they find themselves groomed by a rotation of $`\underset{¯}{\pi }`$ about an axis normal to the line joining them. This is the most attractive configuration. Hence as they pass one another they begin to attract. The two skyrmions now, feeling this attraction, couple and begin to orbit. In order to make the motion clearer in Figure 9, we have decreased the spatial scale of our plot and increased the length of time shown on the orbiting trajectories, which are cut off at $`28.0\mathrm{fm}/\mathrm{c}`$.

In Figure 10 we show energy density contour plots for $`b=0.8\mathrm{fm}`$ to illustrate our point. Notice that the plots cover a longer time period than in the previous cases. As the skyrmions orbit, they radiate. This radiation carries off energy and the skyrmions couple more strongly. The radiation and coupling last for a long time. The first energy level we plot is $`5\mathrm{MeV}/\mathrm{fm}^3`$, comparable to the amplitude of the radiation. The first burst of radiation can be seen in the second and third frames of Figure 10. It appears to consist of the spinning off of a region of high local energy density. A later burst is seen in the last frame. These bursts accompany the orbiting movement of the now bound skyrmions.

The radiation takes away some of the angular momentum and a significant part of the energy. The total energy radiated is greater than the incident kinetic energy, hence the two are now in a bound state. Thus for this choice of parameters, ($`b=0.8\mathrm{fm}`$, $`v=0.5c`$, repulsive grooming) we have the phenomena of orbiting, radiation and capture. We assume that this $`B=2`$ system will eventually find its way to axial-symmetric torus shaped bound state at rest. To reach that state they will have to radiate more energy and the remaining angular momentum. We have followed the orbiting for a time of more than $`60`$ fm/c. We observed continued orbiting with slowly decreasing amplitude, as shown in Figure 11, and continuous energy loss through bursts radiation. We are not able, numerically, to follow the state to the very end.

We now have a bound $`B=2`$ configuration, and therefore should expect the appearance of a torus. The grooming in this case is about the $`x`$ axis. The torus should now appear in the $`yz`$ plane, and should rotate about the $`z`$ axis to carry the initial angular momentum. Energy contours in the $`yz`$ plane corresponding to three frames from Figure 10 are shown in Figure 12. We recognize the familiar pattern of the doughnut-shaped $`B=2`$ bound state. These snapshots correspond to situations when the doughnut is aligned with the $`yz`$ plane. The doughnut is strongly deformed, with the two skyrmions clearly distinguishable. This deformation is only slightly alleviated during the $`20\mathrm{fm}/\mathrm{c}`$ (and a full $`360^\mathrm{o}`$ rotation) between the first and the last frame in Figure 12.

The appearance of the $`B=2`$ bound state is even clearer in the $`b=0.4\mathrm{fm}`$ case. The $`xy`$ plane energy contours are very similar to those in the $`b=0.8\mathrm{fm}`$ process. There is again quite some radiation early on, as illustrated in Figure 13. As can be seen from the trajectories in Figure 9, the attraction is strong in this case (once the skyrmions are past the initial repulsion), and the topological centers come very close.

The $`yz`$ plane contour plot in Figure 14, showing the doughnut just after its formation, exhibits more axial symmetry than the corresponding ones in the $`b=0.8\mathrm{fm}`$ case. We conclude that we have an example of two skyrmions merging into a $`B=2`$ axial symmetric configuration. The angular momentum remaining after the initial radiation burst is carried by the rotation of the doughnut around the $`z`$ axis. The individual skyrmions lose their identity early on and the movement of the topological centers is always very close to the symmetry center. Periodically, they move out of the $`xy`$ plane in the $`z`$ direction. We interpret this as oscillations of the torus. We assume that eventually the kinetic energy and the angular momentum will be radiated away.

In the $`b=1.6\mathrm{fm}`$ case the skyrmions first feel some repulsion, which turns to attraction as they pass. This results in small persistent transverse oscillations. For $`b=2.8\mathrm{fm}`$ there is only a very small attractive interaction.

### D Scattering out of the plane of motion: the attractive (2) channel

The scattering in the last remaining channel, with grooming around the direction of the impact parameter, shows the most remarkable behavior of all. The corresponding trajectories are shown in Figure 15. Recall that for zero impact parameter, this grooming leads to right angle scattering in a direction normal to the plane formed by the incident direction and the grooming axis. That direction is now normal to the scattering plane, which is the one that contains both the impact parameter and the incident momenta. Usually one says that there cannot be scattering, for finite impact parameter, out of the scattering plane by angular momentum conservation. However, we have meson radiation (mostly pion but also some $`\omega `$) that can carry off angular momentum, albeit inefficiently.

For the impact parameter of $`0.4\mathrm{fm}`$ this is just what happens. This remarkable trajectory can be better understood by studying the cartesian components of the topological center’s position as a function of time. We take the $`x`$ direction as the incident one, the impact parameter in the $`y`$ direction and the normal to the scattering plane in the $`z`$ direction. The behavior of $`x`$, $`y`$, and $`z`$ as functions of time for the topological center of each of the skyrmions in the case of $`b=0.4\mathrm{fm}`$ is shown in Figure 16. The skyrmions meet, interact, and then go off normal to the scattering plane in the $`z`$ direction. Therefore they disappear from the plot of trajectories in the $`xy`$ plane, Figure 15. They disappear at the point $`x=y=0`$, with their final state trajectory along the $`z`$ axis having no impact parameter, as it must not, by angular momentum conservation. Note that the slope of the $`z`$ trajectory in Figure 16 is less than that of the $`x`$ trajectory, reflecting the energy carried off by the radiation.

The energy contour plots shown in Figure 17 show the skyrmions coming together and attempting to form a doughnut about a skewed axis in the $`xy`$ plane and then flying apart in the $`z`$ direction. Figure 18 shows the last two configurations at higher resolution in energy. This reveals the radiation that carries off the angular momentum in the $`xy`$ plane.

For the smallest non-zero impact parameter we investigated, $`b=0.4\mathrm{fm}`$, there is meson field radiation left behind in the scattering plane that carries off the initial angular momentum. When we go to an impact parameter of $`0.8\mathrm{fm}`$, things change. The cartesian coordinates for the $`b=0.8\mathrm{fm}`$ case are shown in Figure 19. As the skyrmions meet there is some curvature and then an attempt at uniform motion along $`z`$ at $`x=y=0`$. Now there is too much angular momentum for the field to carry away. The skyrmions try to go off normal to the scattering plane, but they have only a brief excursion in that direction while the meson field is radiating. The skyrmions then return to the scattering plane, but by now the field has taken off so much energy that they are bound and they begin to orbit in the $`xy`$ plane, alternating with excursions in the $`z`$ direction, of slowly decreasing amplitude.

In Figure 15 we have truncated the trajectory of the $`b=0.8`$ fm case at the point where the skyrmions first leave the $`xy`$ plane. Presumably their final state would once again be the static torus. Between $`0.4\mathrm{fm}`$ and $`0.8\mathrm{fm}`$ there must be a critical impact parameter dividing the cases of skyrmions that escape normal to the scattering plane and those that are trapped in orbit in that plane. We are investigating this critical impact parameter and the nature of the solutions in its vicinity.

Energy contour plots for the $`b=0.8\mathrm{fm}`$ impact parameter are shown in Figure 20. The first few frames are very similar to those in Figure 17, showing the skyrmions approaching, forming a doughnut practically in the $`yz`$ plane, and attempting to escape in the $`z`$ direction. As we see in Figure 21, there is considerable radiation at this time as the field tries to carry off the angular momentum. The field cannot, and the skyrmions come back first along $`z`$ and then through the doughnut into the $`xy`$ plane, avoiding the center. Note that both the $`z`$ turn around and the subsequent $`x`$ turn around comes at no more than $`3.0`$ fm separation, as it must. At larger separation, the skyrmions would be out of effective force range and not able to come back. When the skyrmions return to the $`xy`$ plane, they have lost so much energy that they are bound. They begin to orbit, and then make another unsuccessful attempt to escape along $`z`$. As can be seen from the time evolution plot, this sequence of orbiting of the topological centers in the $`xy`$ plane alternating with excursions out of the plane in the $`z`$ direction continues for a long time. We interpret this as follows. The two skyrmions have essentially merged into a $`B=2`$ configuration. The residual angular momentum forces the torus to rotate around the $`z`$ axis. The remainder of the momentum with which the two skyrmions came into this configuration (more precisely, its component pointing to the center), while not enough to push the skyrmions out in the perpendicular direction, results in oscillations which deform the doughnut as illustrated in Figure 22. This central component of the momentum allows the skyrmions to escape in the $`b=0.4\mathrm{fm}`$ and the central case.

We are not able to follow this cycle to completion, but believe that eventually the system will radiate enough energy and angular momentum to settle into the static, axial-symmetric $`B=2`$ torus, similarly to the small impact parameters for repulsive grooming discussed in the previous section. For impact parameters of $`1.6\mathrm{fm}`$ and $`2.8\mathrm{fm}`$ in this channel, the scattering is unremarkable by comparison, as is seen in Figure 15. For the largest impact parameter we notice the medium-range repulsive interaction, just the opposite of the corresponding case in the repulsive grooming discussed above. The two groomings are interchanged at the point where the centers cross the $`x=0`$ plane.

## V Conclusion and Outlook

The phenomenology of two-skyrmion scattering for various groomings as it emerges from our qualitative study, may be summarized as follows. In the absence of grooming , the scattering process is almost elastic at the energy we considered. The kinematics is similar to that of point-particles interacting via a central repulsive interaction. This is to be expected, since the interaction between two hedgehog configurations is central. We know  that at higher velocities one expects to excite vibrational modes of the individual skyrmions. We do see modest indications of that.

The head-on collision in the repulsive channel is quasi-elastic, similar to the hedgehog-hedgehog processes. The remaining collisions involving grooming of $`180^o`$ can be divided into two categories, depending of whether or not the impact parameter, $`b`$, is small enough for the formation of the $`B=2`$ (torus-shaped) bound state. If the impact parameter is large, the collisions are quasi-elastic with a weak attractive or repulsive character depending on the grooming. If the impact parameter is small enough, the collision proceeds through a (sometimes deformed) doughnut configuration in the plane normal to the grooming axis, even if the final state is not bound.

The cleanest example is the attractive (1) channel, with the grooming axis normal to the plane of motion. For zero impact parameter it is well known  that a torus is formed in the $`xy`$ plane, the skyrmions lose their identity, and fly out at right angles with respect to the incident direction and the grooming axis. For a nonzero impact parameter, a deformed torus still appears in the $`xy`$ plane, and the scattering angle decreases continuously from $`90^o`$ as $`b`$ increases.

In the repulsive case (grooming around the direction of motion, $`x`$), the torus is formed close to the $`yz`$ plane. As the two incident skyrmions approach the $`y`$ axis (that of the impact parameter), they find themselves in an attractive configuration, since they are now groomed around an axis ($`x`$) perpendicular to the one ($`y`$) connecting them. They form the doughnut initially in the plane perpendicular to the grooming axis. The configuration carries some of the initial angular momentum by rotating around the $`z`$ axis. The skyrmions do not tend to exit in the perpendicular direction because they came in with very little momentum in the plane of the torus.

In the attractive (1) and (2) cases, the initial momenta have a large component pointing to the center of the torus. This momentum is channeled into the perpendicular direction, always leading to scattering in the attractive (1) case. In the attractive (2) case, with grooming around the impact parameter direction $`y`$, the doughnut is formed close to the $`xz`$ plane. The $`90^o`$ scattering would therefore happen in the $`z`$ direction, but this is strongly limited by the necessity of shedding the angular momentum around $`z`$. Radiation provides a mechanism for this, and for a small enough impact parameter the skyrmions can escape in the $`z`$ direction. Otherwise they go into a rotating doughnut configuration. In this case however, there is radiation and significant momentum in the plane of the torus, which torus also exhibits quadrupole oscillations.

While our investigation is by no means complete, we did identify a significant number of distinct patterns for the outcome of skyrmion-skyrmion collisions. One may define a number of critical configurations which separate the different outcomes, even for fixed velocity. The picture might get even richer by sweeping a range of incident velocities.

The calculations of skyrmion-skyrmion scattering reported in the Section above demonstrate two things. First they show that the $`\omega `$ meson stabilizes the Skyrme model and makes it numerically tractable out to long times and through complex dynamical situations. Second the calculation shows a rich mix of phenomena. These include capture and orbiting, radiation, and excursions out of the scattering plane. Although each of these arises naturally in the model, they have not been seen or demonstrated before, nor have most of them been suggested before. They therefore add to the rich and often surprising mix of results in the Skyrme model. Since this is a model of low energy QCD at large $`N_C`$, the new results give insight into aspects of QCD in the non-perturbative long wave length or low energy domain. This is a region in which our best hope for insight comes from effective theories like the Skyrme model.

The success of these calculations also gives us confidence that the method can be carried over, with only simple changes, to the annihilation problem. In particular the stability of the calculations, thanks both to $`\omega `$ meson stabilization and a number of numerical strategies, suggests that the annihilation calculation will also be stable. For annihilation, meson field (pion and omega) radiation in the final state is all there is and the fact that we clearly see radiation in the skyrmion-skyrmion case is reassuring.

Our work here suggests a number of further avenues. We have already discussed annihilation. It would also be interesting to explore the landscape of skyrmion-skyrmion scattering as a function of energy as well as grooming and impact parameter and also to examine more closely the behavior of the model in the neighborhood of critical parameters where the scattering behavior changes abruptly. Also interesting would be to try to extract information about nucleon-nucleon scattering from the results for skyrmion-skyrmion scattering. We are investigating these questions.

## VI Acknowledgements

We are grateful to Jac Verbaarschot and Yang Lu for many helpful discussions on both physics and numerical issues. We are indebted to Folkert Tangerman and Monica Holbooke for a number of suggestions regarding the calculation. Jac Verbaarschot is also acknowledged for a critical reading of the manuscript.

This work was supported in part by the National Science Foundation. We are very grateful to Prof. R. Hollebeek for making available to us the considerable computing resources of the National Scalable Cluster Project Center at the University of Pennsylvania, which center is also supported by the National Science Foundation.