Rolling Skyrmions and the Nuclear Spin-Orbit Force

We compute the nuclear spin-orbit coupling from the Skyrme model. Previous attempts to do this were based on the product ansatz, and as such were limited to a system of two well-separated nuclei. Our calculation utilises a new method, and is applicable to the phenomenologically important situation of a single nucleon orbiting a large nucleus. We find that, to second order in perturbation theory, the coefficient of the spin-orbit coupling induced by pion field interactions has the wrong sign, but as the strength of the pion-nucleon interactions increases the correct sign is recovered non-perturbatively.


Introduction
The spin-orbit coupling is an important ingredient in nuclear structure theory. Its presence implies that it is energetically favourable for the spin and orbital angular momentum of a nucleon to be aligned, particularly if this nucleon is moving close to the surface of a larger nucleus. This explains the phenomenon of magic numbers, and it is important in the description of halo nuclei, to name just two examples. Unlike the spin-orbit force encountered in the study of electron shells of an atom, the nuclear spin-orbit force is not merely a relativistic effect but is caused by the strong interaction physics of nuclei.
The Skyrme model is an effective description of QCD, and a candidate model of nuclei with a topologically conserved baryon number. It successfully accounts for phenomena such as the stability of the alpha-particle, the long-range forces between nuclei, and quantum numbers of excited states of very light nuclei. Some of the recent successes of the model include reproducing the excited states of oxygen-16 [1] and carbon-12 [2], nuclear binding energies of the correct magnitude [3], accurately modelling neutron stars [4] and a geometric explanation for certain magic nuclei [5].
However, one of the challenges in analysing the Skyrme model has been accounting for the spin-orbit coupling. There have been several attempts to calculate the spinorbit term in the nucleon-nucleon potential [6,7,8,9,10]. Most of these calculations were only valid for large separations and were also perturbative, and so corresponded to calculations taking into account one-and two-pion exchange. Almost all obtained the nucleon-nucleon spin-orbit coupling with the wrong sign, although [6] obtained the correct sign by introducing additional mesons in the model.
The conventional description of the spin-orbit force is in the framework of relativistic mean field theory [11], which couples nucleons to several mesons (including the pion, σ, ρ and ω). An interesting perspective was put forward by Kaiser and Weise [12]: they argued that the spin-orbit coupling receives several contributions, including a wrong-sign contribution from pion exchange; this is compensated by other effects, including meson exchange and three-body forces. This seems to be related to the sign problem in the Skyrme model.
In this article we investigate in a novel way how a short-range spin-orbit coupling arises in the Skyrme model. Unlike relativistic mean field theory, our calculation is non-relativistic and incorporates pions but no other mesons. Our calculations are for a somewhat simplified model, but we hope this model captures the essence of the effect. Our main discovery is that the sign of the spin-orbit coupling is wrong at weak coupling, where a perturbative approach would be valid. However, the sign is correct when the coupling between the nucleon and the surface of the nucleus with which it interacts is stronger.
A key property of a Skyrmion, distinguishing it from an elementary nucleon, is that it has orientational degrees of freedom. It is a spherical rigid rotor. After quantisation [13], the basic states are nucleons with spin 1 2 , but there are also excited states with spin 3 2 corresponding to Delta resonances, and further states of higher spin and higher energy that play no significant role. The states simultaneously have isospin quantum numbers (isospin 1 2 for the nucleons and 3 2 for the Deltas). In our model, a dynamical Skyrmion interacts quantum mechanically with a background multi-Skyrmion field modelling the nuclear surface. The interaction involves a potential that depends on the Skyrmion orientation and its position, and the potential has a strength parameter that we consider as adjustable. When the parameter is small, a perturbative treatment works. However, the spin-orbit coupling has the wrong sign in this regime. When the parameter is larger (but not too large), the spin-orbit coupling for the Skyrmion has the correct sign.
Indeed, in this latter regime, a better approximation to the Skyrmion wavefunction is to say that the orientation has its probability concentrated near the minimum of the orientational potential, with this minimum varying with the Skyrmion's location on the surface. The quantum state is now close to the classical picture of a Skyrmion rolling over the nuclear surface, maintaining a minimal orientational potential energy. This classical rolling motion gives the correct sign for the spin-orbit coupling. In earlier work, Halcrow and one of the present authors investigated a model of this type [14], but they only treated the case of a disc interacting with another disc in two dimensions. When the potential is strong, the model becomes a quantised version of cog wheels rolling around each other. Here we do better, by treating a realistic three-dimensional Skyrmion interacting quantum mechanically with a nuclear surface. However, we still need to make various approximations. For example, we assume the height of the Skyrmion above the surface is fixed.
Our analysis is based on the following well-known interpretation of the phenomenological spin-orbit coupling. Consider a nucleon near the surface of a spherical nucleus. Suppose that in addition to the usual kinetic terms, the hamiltonian for the nucleon contains a term of the form a S. N × P , where a is a parameter, S is the spin of the nucleon, P is its momentum, and N is an inward-pointing vector normal to the surface, which may be interpreted as the gradient of the density of nuclear matter. Since the position vector r of the nucleon equals −r N /| N |, this term equals −(a| N |/r) S. L, where L = r × P is the orbital angular momentum of the nucleon. This is the usual form of the spin-orbit coupling.
In order to give the correct magic numbers, the spin-orbit coupling must prefer spin and angular momentum to be aligned rather than anti-aligned, so the parameter a needs to be positive. The advantage for us of the formula (1.1) is that it applies when the nucleon is interacting with an essentially flat nuclear surface, as in the model we will discuss below. We will refer to (1.1) as the spin-momentum coupling. Note that N is essential here, and implies that there is no coupling for an isolated nucleon, nor for a nucleon deeply embedded inside a nucleus. There are two practical difficulties with this approach: the first is that the interaction between Skyrmions and multi-Skyrmions is poorly understood at short distances, and the second is that the complicated spatial structure of known multi-Skyrmions with finite baryon number would make the calculations laborious. We solve the first of these problems by working in the lightly bound version of the Skyrme model [15], for which multi-Skyrmions and their interactions are accurately captured by a point particle description, although the particles still have orientational degrees of freedom. We solve the second problem by supposing that the multi-Skyrmion representing the core of the nucleus is large, and approximating its surface by a plane. Since Skyrmions in the lightly bound model naturally arrange themselves to sit at vertices of an FCC lattice, this surface has a high degree of symmetry, making the calculation tractable.
In the next section we review the 2D toy model of [14], but in a modified and simplified form. Here the dynamical, Skyrmion-like object is a coloured disc, and it moves in the background of a straight, periodically coloured rail, rather than around a larger coloured disc as in [14]. The potential depends on the colour difference between the disc and the rail at their closest points. The translational and rotational motion of the disc is quantised, and we compare the result of a perturbative treatment, valid θ x X Figure 1: Coloured disc on a fixed coloured rail. One period of the rail colouring is shown.
when the potential is weak but which leads to a spin-momentum coupling of the wrong sign, with a non-perturbative approach that can deal with stronger coupling but is still algebraically straightforward. The price to pay for working non-perturbatively is that we must assume that the moment of inertia of the disc is small; in our perturbative calculation, no such assumption is necessary. The strong coupling result gives the correct sign for the spin-momentum coupling. In the later sections we perform similar calculations in the more realistic 3D setting. Here, the Skyrmion is visualised as a coloured sphere moving relative to a coloured surface, and the potential again depends on the colour difference at the closest points. The calculations can be done by hand, exploiting the assumed lattice symmetries of the (planar) nuclear surface, but are nevertheless considerably more complicated. The reader may wish to skip the details here.

Disc on a rail
We start with a two-dimensional toy model of spin-momentum coupling, rather similar to what was analysed in [14]. Consider a vertical disc at a fixed height above a fixed, straight rail. The disc can move along the rail and also rotate. Both the edge of the disc and the rail are coloured, and the potential energy is a periodic function of the colour difference at their closest points. When the colours match, the potential energy is lowest. Let us assume that the disc is coloured so that for the potential to remain at its lowest value as the disc moves classically, the disc needs to roll along the rail. See Figure 1. This model is similar to a cog on a rack rail, which can only roll, but not slip. Classically there is spin-momentum coupling, as the (clockwise) spin of a rolling cog is a positive multiple of its linear momentum. Let X be a linear coordinate along the rail. The colour χ along the rail is an angular field variable, and as with an ordinary angle we assume χ takes any real value and identify values that differ by 2π. We suppose that χ = X, so the colour is periodic along the rail, with period 2π. Let the disc have radius 1 and assume that when it is in its standard orientation, the colour is the same as the angle around the disc measured from the bottom in an anticlockwise direction, i.e. the colour is χ at angle χ.
Suppose now that the position and orientation of the disc are (x, θ), where x is the location of the centre of the disc, projected down to the X-axis, and θ is the angle by which the disc is rotated clockwise relative to its standard orientation. The bottom of the disc then has colour θ, and the rail under this point has colour x. We suppose the potential energy of the disc in this configuration is We next introduce some dynamics. Suppose the disc has unit mass, and moment of inertia Λ, so the Lagrangian for its motion is The equations of motion arë Note that as the potential only depends on x−θ, there is a conserved quantityẋ+Λθ. One solution of the equations is x = µt, θ = µt for any constant µ -this is rolling motion.
The conjugate momenta to x and θ are

3) and the Hamiltonian is
with conserved quantity p + s. We now quantise. Stationary wavefunctions are of the form Ψ(x, θ), and the momentum and spin operators are The stationary Schrödinger equation is where the operator on the left hand side is the Hamiltonian (2.4) expressed in terms of the momentum and spin operators. The configuration space of the disc has first homotopy group Z, so wavefunctions can acquire a phase when θ → θ + 2π. Bearing in mind that we are modelling a fermionic nucleon interacting with a large nucleus, we choose this phase to be π. Wavefunctions then have a Fourier expansion a superposition of half-integer spin states. The free motion, in the absence of the potential, has separately conserved momentum p and spin s, and the basic stationary state is where p is arbitrary and s is half-integer. This state has energy We now suppose that Λ is small, so that 1 Λ is large compared to V 0 and to p 2 . The expressions we derive later will only be valid provided p 2 1 Λ . In this regime, the low energy states are those with s = ± 1 2 . This is physically what we are interested in. Spin 3 2 nucleons (i.e. Delta resonances) have energy about 300 MeV greater than spin 1 2 nucleons, and spin-orbit energies are much less than this, of order 1 MeV. So we mostly neglect the small parts of the wavefunction with s = ± 3 2 or larger. Because of the restriction to n = ±1 states, i.e. those with s = ± 1 2 , the wavefunction reduces to Ψ(x, θ) = ψ 1 (x)e i 1 2 θ + ψ −1 (x)e −i 1 2 θ . (2.10) A stationary state like this is not strictly compatible with the Schrödinger equation, because the potential couples it to s = ± 3 2 states. We can deal with this by calculating the matrix form of the Hamiltonian restricted to this subspace of wavefunctions. Recall that there is the conserved quantity p + s. This implies that if ψ 1 (x) = e ipx then ψ −1 (x) = Ae ip x , where p = p + 1, for some amplitude A. Momentum p itself is not a good label for states, but instead we can use r = p + s, where r takes any value in the range (−∞, ∞). The wavefunction (2.10) becomes, for a definite value of r, Alternatively, the crystal momentum k could be defined to be p mod 1 and the (first) Brillouin zone to be − 1 2 ≤ k ≤ 1 2 , but because of the restricted range of spins, we do not need the formalism of Bloch states mixing momentum p with all its integer shifts.
We now work with basis states 1 2π e i(r− 1 2 )x e i 1 2 θ and 1 2π e i(r+ 1 2 )x e −i 1 2 θ . These are normalised in {0 ≤ x ≤ 2π , 0 ≤ θ ≤ 2π}. The matrix elements of the Hamiltonian (2.4), or equivalently the operator on the left of (2.6), are where the diagonal terms are kinetic contributions. The upper off-diagonal term comes from the matrix element of the potential and the lower off-diagonal term is the same, by hermiticity. The potential makes no contribution to the diagonal terms.
It is now convenient to express the energy eigenvalues E of H 2×2 as E = 1 2 ε + 1 8Λ . The matrix with eigenvalues ε is 14) and the eigenvalue equation det( H 2×2 − ε1) = 0 reduces to with solutions ε ± (r) = r 2 + 1 4 ± r 2 + V 2 0 . (2.16) The spectrum has two branches, the lower branch ε − (r) and the upper branch ε + (r), and is symmetric under r → −r. When V 0 = 0 the spectrum simplifies to ε(r) = (r ± 1 2 ) 2 , whose graph consists of two intersecting parabolas, with minima at r = − 1 2 and r = 1 2 , and a crossover at r = 0. We are mainly interested in low energy states on the lower branch, near the minima of ε − (r). There is an important bifurcation at a critical strength of the potential, , ε − has two minima at r = ± 1 4 − V 2 0 and a local maximum at r = 0. For V 0 > 1 2 , there is just one minimum at r = 0; here p = ± 1 2 , so the crystal momentum k is located on the boundary of the Brillouin zone. The upper branch ε + (r) has simpler behaviour, as it just has a minimum at r = 0 for all positive V 0 . Figure 2 shows graphs of the eigenvalue spectrum for two typical values of V 0 .
Recall the form of the (unnormalised) wavefunction (2.11). We evaluate A using the condition that 1 A is the eigenvector of the matrix (2.14). On the lower branch of the spectrum and on the upper branch Note that |A| 2 = 1 at r = 0 on both branches, so spins ± 1 2 are superposed there with equal probability. On the lower branch, the total (unnormalised) wavefunction at r = 0 is 2 cos 1 2 (x − θ), so the highest probability occurs for θ = x, where the disc is oriented so as to minimise the potential energy. This is compatible with a rolling motion.
Except in cases where |A| is very close to 0, or much larger than 1, the quantum states of the disc cannot be thought of as having a definite momentum p or spin s, because the potential strongly superposes states where these have different values. So to consider the correlation between the momentum and spin, we work with their expectation values P and S .
The expectation value of the spin is where A is given by expressions (2.17) and (2.18), respectively, on the lower and upper branches. The expectation value of momentum follows immediately, as p + s = r for both contributing states in (2.11), so Graphs of P and S as functions of r are shown in Figure 3. They are plotted together with ε − for states on the lower branch, for the typical values of V 0 we selected before; for states on the upper branch, they are plotted together with ε + .  Interesting to note is that P vanishes wherever E (or equivalently ε) is stationary with respect to r, as one can see from the graphs. This is because and the right hand side is the matrix form of the momentum operator. Taking expectation values gives the result. (One also needs to use the identity d dr Ψ|Ψ + Ψ| d dr Ψ = 0 for normalised states.) When the potential is relatively weak, such that V 0 < 1 2 , then P passes through 0 at the non-zero minima of ε − on the lower energy branch. On the other hand S does not change sign near here. The signs of P and S are therefore not strongly correlated for these low energy states, and we conclude that for weak potentials there is no significant spin-momentum coupling. Near r = 0, where ε − has a local maximum, P and S have opposite signs, so momentum and spin are anticorrelated. This is the opposite of the classical correlation of momentum and spin for a rolling motion. Similarly, on the upper energy branch, the expectations of momentum and spin have opposite signs for all r, so they are anticorrelated.
When the potential is stronger, such that V 0 > 1 2 , we find the correlation we are seeking. Here, the low energy states on the lower branch are near r = 0, and we see that P and S have the same sign. It is straightforward to estimate these quantities analytically for small r. They are and for V 0 > 1 2 both their slopes with respect to r are positive. In fact, because P is zero only at r = 0 when V 0 > 1 2 , there is a spin-momentum correlation of the desired sign for all r, on the lower branch. On the other hand, the momentum and spin are anticorrelated for all r on the upper branch.
The conclusion is that the potential has to be quite strong to achieve the spinmomentum coupling for quantum states that mimics the classical phenomenon of rolling motion for a cog. As in the usual model of spin-orbit coupling for a spin 1 2 particle, there are two states, a lower energy state with a positive correlation, and a higher energy state with an anticorrelation.

Perturbation theory
We have just seen that the spin-momentum correlation has the desired form only when the potential is quite strong. Nevertheless, it is of some interest to calculate what happens in perturbation theory. When the potential is weak, we can calculate the energy spectrum to second order in perturbation theory, treating V 0 as small. The perturbative result overlaps what we have already calculated, and we can allow for the possibility that the moment of inertia Λ is not small. This is a useful check on our calculations, both for the two-dimensional disc, and later, when we consider three-dimensional Skyrmion dynamics.
When V 0 = 0, the eigenstates of the Hamiltonian are Ψ 0 (x, θ) = e ipx e isθ , with definite momentum and spin, and energy E = 1 2 p 2 + 1 2Λ s 2 . Low energy states are those with s = ± 1 2 and p 0. These are near the centre of the Brillouin zone. Let us focus on the states with s = 1 2 (the results are similar for s = − 1 2 ), whose energy is Recall that when the potential is included, there is still the good quantum number r = p + s, so the states that we are focussing on have r = p + 1 2 1 2 . The effect of the cosine potential −V 0 cos(x − θ), at leading order, is to mix the unperturbed state Ψ 0 = e ipx e i 1 2 θ with states where p is shifted by ±1, i.e. the states e i(p+1)x e −i 1 2 θ and e i(p−1)x e i 3 2 θ , whose unperturbed energies are 1 2 (p + 1) 2 + 1 8Λ and 1 2 (p − 1) 2 + 9 8Λ , respectively. The potential has no diagonal matrix element, so the energy is unchanged to first order in V 0 .
The eigenfunction of the Hamiltonian to first order in V 0 , for the fixed value p + 1 2 of r, is where the denominators of the coefficients are proportional to differences between the energies of the unperturbed states. The energy of the state Ψ, to second order in V 0 (found either by acting with the Hamiltonian, or by using the standard formula) is .

(2.24)
This formula is valid, provided the unperturbed energy differences are not small compared to V 0 . So V 0 must be much less than 1 and p must not approach − 1 2 . The perturbative approach therefore definitely fails for the states near r = 0 that we were considering earlier for fairly strong V 0 . However, it is successful for small p, even if Λ is not small and the last term of the formula (2.24) makes a significant contribution. Therefore, perturbation theory allows us to consider easily the spin 3 2 contribution to low energy states, in contrast to our matrix method, which required this contribution to be negligible.
Let us compare our previous calculation of ε − , as a matrix eigenvalue, with this perturbative estimate. From the expression (2.16), and converting it back to give the energy E as a function of momentum p, we find, to second order in V 0 , that and this agrees with (2.24) provided Λ is small. So the matrix method and perturbation theory agree where they should. The conclusion is that perturbation theory is a good way to find states of the disc in a certain regime, but that regime does not extend to where spin-momentum coupling has the correlation we are seeking. In the following sections we shall investigate the quantised three-dimensional dynamics of a Skyrmion in a background potential. We should expect the matrix method to be more effective than perturbation theory for finding the desired form of spin-momentum coupling. We shall need a model where the potential is fairly strong, and where states of the Skyrmion with spin 3 2 and higher are suppressed, relative to the spin 1 2 states.
A multi-Skyrmion modelling a nucleus with mass number N is described by N Skyrmion-like point particles, each with three positional degrees of freedom and three orientational degrees of freedom. The rotational degrees of freedom could be expressed using an SO(3) matrix, but for quantum mechanical calculations it is more convenient to use an SU(2) matrix q. Throughout this section we will identify the group SU(2) with the group of unit quaternions, making the identifications i = −iσ 1 , j = −iσ 2 , k = −iσ 3 between imaginary quaternions and Pauli matrices. The Lagrangian for the model consists of standard kinetic terms for the positions and orientations, and interaction potentials between pairs of particles (see [15] for the precise form).
The interaction potential is such that the particles tend to arrange themselves into crystals with an FCC lattice structure, with a preferred orientation at each lattice site. In suitable length units the FCC lattice is the set of vectors (x, y, z) ∈ Z 3 such that x + y + z is even. The preferred orientation at lattice site (x, y, z) is i x j y k z .
We want to study the problem of a charge-1 Skyrmion rolling along the surface of a half-filled FCC lattice. See Figure 4. We assume that the lattice sites with x + y + z ≤ −2 are filled with particles in their preferred orientations, and consider a Skyrmion moving freely in the plane The degrees of freedom for this Skyrmion are its position coordinates (x, y, z) and its orientation q ∈ SU(2). Its dynamics can be described by a Lagrangian consisting of a standard kinetic term and a potential function V : Π × SU(2) → R. The kinetic terms are invariant under the group SU(2) I × SU(2) S of isorotations and rotations, with action where we identify vectors x ∈ R 3 with imaginary quaternions xi + yj + zk = −i(xσ 1 + yσ 2 + zσ 3 ). These terms are also invariant under translations ( x, q) → ( x + c, q) and parity transformations ( x, q) → (− x, q).
The potential function V must be invariant under the group of symmetries of the half-filled lattice. This group is generated by the following transformations: The invariance of V under −τ 2 , −ρτ 2 ρ −1 and −ρ 2 τ 2 ρ −2 implies in particular that V is invariant under the translation action of the two-dimensional lattice The transformations listed above acting on (Π/Γ) × SU(2) generate a finite group which is isomorphic to the binary cubic group (the double cover of the cubic group). Note that since τ 2 = −1 when acting on (Π/Γ) × SU (2) we are free to use ρ, σ, τ as a set of generators.
We employ an ansatz for the potential of the form with R(q) the rotation matrix induced by q (i.e. qσ j q −1 = σ i R(q) ij ), and Y ( x) a 3 × 3 matrix-valued function. This ansatz is motivated by the dipole description of Skyrmion interactions; to a good approximation a single Skyrmion interacts with a background field of pions like a triple of orthogonal scalar dipoles, and this dipole interaction has similar q-dependence to our ansatz. Alternatively, one may regard our ansatz as the first two terms in an expansion of V in harmonics on SU (2). Note that this potential satisfies V ( x, −q) = V ( x, q), as required by symmetry. We simplify the ansatz further using Fourier series. Both U and Y are required to be invariant under the lattice Γ, so have Fourier series with summands corresponding to dual lattice vectors. We assume that these Fourier series only contain terms corresponding to the shortest dual lattice vectors; the associated functions are 1 and e ±i a j . x , where The symmetries ρ, σ, τ generate an action of the binary cubic group on the vector spaces occupied by U, Y . Since the ansatz (3.8) is invariant under q → −q, this action descends to an action of the cubic group. Representation theory can be used to find all functions U and Y which are invariant under this action. This calculation involves the irreducible representations of the cubic group: we recall these briefly. Besides the trivial representation A 1 , there is another one-dimensional representation A 2 in which ρ and τ map to 1 and σ maps to −1. There is a unique two-dimensional representation E and two three-dimensional representations T 1 and T 2 ; the first of these is the standard rotational action as the symmetry group of the cube and the second is It can be shown that the functions e ±i a j . x transform in the representation 2T 2 of the cubic group; since this contains no trivial subrepresentations the only allowed form for U is a constant function. Since this constant does not alter differences between energy eigenvalues we set it to zero.
The elements of the group act on matrix-valued functions Y by simultaneously multiplying with matrices from the left and right, and by permuting the Fourier modes. The matrix acting from the left corresponds to the representation A 2 ⊕ E, and that acting from the right corresponds to the representation T 1 . The action on Fourier modes is A 1 ⊕ 2T 2 . Therefore the representation acting on the vector space space contains four copies of A 1 , so the space of allowed potential functions has real dimension four.
This space of allowed potential functions can be parametrised by (U 0 , U 1 ) = (W 0 e iθ 0 , W 1 e iθ 1 ) ∈ C 2 as follows: The values of the constants can be estimated in the lightly bound Skyrme model using its point particle approximation. One calculates a function Y true by adding up the interaction energies between fixed Skyrmions in the planar lattice x + y + z = −2 and the Skyrmion moving freely in the plane x + y + z = 0, and then calculates its Fourier coefficients. The values obtained are With these parameters the truncated Fourier series Y approx given in eq. (3.10) is a good approximation to Y true , in the sense that the ratio of the squares of the L 2 norms of Y true − Y approx and Y true is 0.095. Our final potential V ( x, q) = Tr(R(q)Y approx ( x)) Figure 5: The path of a rolling Skyrmion.
is not exact, even in the point particle description of Skyrmion interactions, but it is analogous to the potential −V 0 cos(x − θ) that we chose for the disc in section 2. We claim that for the parameter values (3.11), the potential given by equations (3.8) and (3.10) induces classical motion similar to a ball rolling on a surface. Consider the situation where a particle moves from (x, y, z, q) = (0, 0, 0, 1) to (x, y, z, q) = (1, −1, 0, ±k). Both of these points are critical points of the potential, and for our parameter set they are minima. We will treat this situation adiabatically, assuming that the mass M of the Skyrmion is much greater than its moment of inertia Λ. If the spatial kinetic energy 1 2 M v 2 is much larger than the energy scale W = W 2 0 + W 2 1 of the potential then the path in space will to a good approximation be a straight line: x(t) = (vt/ √ 2)(1, −1, 0). If the velocity is not too large then, at each time t, q(t) will to a good approximation be the orientation that minimises V ( x(t), q) with respect to variations in q. In this situation the rotational kinetic energy is roughly 1 2 Λv 2 , and the approximation is reliable as long as this is much less than W . Thus our approximation assumes that W/M v 2 W/Λ. We wish to compare this motion with that of a rolling ball. If a ball of radius r rolls with velocity v along a surface with inward-pointing unit normal n its angular velocity will be ω = − n × v/r. For n = (−1, −1, −1)/ √ 3 and v = (v/ √ 2)(1, −1, 0) as above this makes the angular velocity a positive multiple of (1, 1, −2). The angle θ(t) between the angular velocity vector ω(t) = −2q −1q for the path q(t) and the vector (1, 1, −2) measures deviation from rolling motion: acute angles indicate motion similar to rolling, and obtuse angles indicate motion that is opposed to rolling. We have computed q(t) using the adiabatic approximation described above and have hence determined θ(t). The maximum angle along the path is 0.89 ≈ 2π/7, indicating that the motion induced by the potential is similar to that of a rolling ball.
This adiabatic motion of a Skyrmion is illustrated in Figure 5 (see also Figure  4). The orientation of the rolling Skyrmion is illustrated at the start, mid-point, and end of the path. The start and end points are neighbouring lattice sites, and their orientations differ by a rotation of 180 degrees about the red-green axis. The most natural guess for the orientation at the mid-point is a rotation through 90 degrees about the same axis, and there are two possibilities here (depending on whether one rotates clockwise or anticlockwise). Figure 5 shows the orientation for one sense of rotation, but the alternative would have made the Skyrmion's red, white and yellow faces visible at the mid-point. Now observe that just below the Skyrmion at the mid-point there is a nearby Skyrmion in the lattice (white and yellow faces visible). It is straightforward to find the pion dipole fields of this pair of Skyrmions along the line joining them and verify that for the illustrated sense of rotation, the fields are identical at the closest points, implying that the potential energy is minimal. (The associated colouring is predominantly green, but with a small tilt towards white and yellow.) If the sense of rotation had been opposite, the field match would have been less good and the energy greater.
We conclude that the rolling motion illustrated in Figure 5 is along a particularly deep valley in the potential energy landscape, and favoured as a low energy classical motion. Anti-rolling is disfavoured. Figure 5 suggests that to a good approximation the spin vector S for the rolling Skyrmion points in the direction of the red-green axis, from green to red. This spin vector S, the vector N pointing into the halffilled lattice, and the momentum vector P do not form an orthonormal triad ( P is orthogonal to both N and S , and N . S = −1/ √ 3), but their triple scalar product S. N × P is negative. This is what is expected classically if the parameter a in equation (1.1) is positive.

Weak coupling to the potential
In this section and the next we will study the quantum mechanical problem of a Skyrmion interacting with the surface of a half-filled lattice. Since the potential experienced by the Skyrmion is periodic it is natural to analyse this problem using the theory of Bloch waves. Let k ∈ R 3 be a crystal wavevector satisfying k 1 + k 2 + k 3 = 0 and let H k be the Hilbert space of wavefunctions Ψ : The first condition ensures that Ψ is effectively defined in the plane Π rather than all of R 3 . The second condition has the implication that two crystal wavevectors whose difference lies in the reciprocal lattice Γ * generated by a j define the same Hilbert space, so k should be regarded as an element of Π * /Γ * . The natural operators on H k are isospin, spin, and momentum. Spin and isospin are just the infinitesimal versions of the actions described in (3.2): Although the space in which the Skyrmion moves is two-dimensional, it will be convenient to write momentum as a three-vector (due to the three-dimensional origin of the problem). Thus we set noting that P 1 + P 2 + P 3 = 0. Then for a plane wave of the form Ψ( It will be useful in what follows to decompose H k into eigenspaces of | S| 2 . Fix a non-negative integer or half-integer and let η : SU(2) → SU(2 + 1) be the spin irreducible representation of SU(2). If ψ : Π → Mat(2 + 1, C) is a matrix-valued function of x then Ψ( x, q) := Tr(ψ( x)η (q)) (4.6) satisfies | S| 2 Ψ = | I| 2 Ψ = ( + 1)Ψ. Thus this wavefunction describes a particle of total spin and total isospin . The space of all such wavefunctions in H k will be denoted H k . The Peter-Weyl theorem implies that any wavefunction in H k can be decomposed as an infinite sum of wavefunctions of this type: The wavefunction Ψ describing the Skyrmion is required to satisfy the Finkelstein-Rubinstein constraints [13]. These simply state that Ψ is an odd function of q: Ψ( x, −q) = −Ψ( x, q). Functions in H k are odd if is a half-integer and even if is an integer. Thus the Finkelstein-Rubinstein constraints require Ψ to be in the subspace H odd k of H k , where the summation over is restricted to half-integers. This ensures the quantised Skyrmion has half-integer spin.

Outline of perturbation theory
The hamiltonian that we will study is Here M, Λ > 0 are parameters representing the mass and moment of inertia of the Skyrmion, and V is the potential introduced in the previous section. We will construct an effective hamiltonian for the lowest-energy eigenstates using perturbation theory, with the parameters W 0 and W 1 of V treated as small. If V = 0 and k is in the first Brillouin zone (i.e. | k + v| > | k| for all v ∈ Γ * ) then the lowest energy eigenstates in the space H k are clearly of the form Ψ 0 ( x, q) = Tr(ψq)e i k. x , (4.9) with ψ ∈ Mat(2, C). The space of all such wavefunctions has dimension four and will be denoted by K k . The energy of these states is This is minimised by k = 0. In the following calculations we will assume that k is close to 0, discarding terms of O( k 2 ). When the potential V is non-zero, the four degenerate energy levels with energy E 0 will separate. We will study this effect using perturbation theory. Let us review the overall methodology, which generalises the formulae (2.23) and (2.24). We seek an operator I : K k → H k which depends continuously on the parameters U 0 , U 1 in the potential, such that the image under I of the H 0 -invariant subspace K k is Hinvariant, and such that the composition Π K I of I with the orthogonal projection Π K : H k → K k is the identity map. The effective hamiltonian is then defined to be H eff = Π K HI. The operators I and H eff will be constructed as power series in the parameters that appear in the potential.
To zeroth order, I is just the inclusion: I|Ψ 0 = |Ψ 0 + O(V ) for all Ψ 0 ∈ H k . The first order correction to H eff is given by (4.11) The term linear in V vanishes. The reason for this is simple: the only non-zero terms in the Fourier series of V Ψ 0 correspond to plane waves of the form e i( k± a j ). x , as one sees from eqs. (4.9) and (3.10), and these are all L 2 -orthogonal to e i k. x . As a consequence, H eff = H 0 + O(V 2 ). The first order correction to I is given by This satisfies HI|Ψ 0 = E 0 I|Ψ 0 + O(V 2 ), so its image is H-invariant up to terms quadratic in V . The second order correction to H eff is given by In the next subsection we will calculate the action of Π K V (H 0 − E 0 ) −1 V on wavefunctions Ψ 0 of the form (4.9), and thereby evaluate H eff to second order. A reader uninterested in the details of this calculation may skip to the final result, eq. (4.34).

The effective hamiltonian H eff
We begin by analysing V |Ψ 0 , with the potential V given by eqs. (3.8) and (3.10).
From eq. (3.8) we see that V ∈ H 1 0 , and from eq. (4.9) we see that Ψ 0 ∈ H 1 2 k . It follows from the Clebsch-Gordan rules that the excited wavefunction V ( x, q)Ψ 0 ( x, q) will be a sum of terms with spin 1 2 and spin 3 2 . Thus with Π denoting projection onto H k . Applying (H 0 − E 0 ) −1 to the spin 1 2 term gives As the Fourier modes that appear in the excited wavefunction V ( x, q)Ψ 0 ( x, q) are e i( k± a j ). x , the operator | P − k| 2 takes the constant value | a j | 2 = 2π 2 /3 on Π 1 2 V |Ψ 0 , which simplifies this expression. The spin 3 2 term can be analysed in the same way, yielding Thus to compute Π K V (H 0 −E 0 ) −1 V |Ψ 0 we need to compute the following four terms: This can be evaluated with the help of the following identity, which is proved in the appendix: We introduce a vector with the index i − j understood modulo 3. Then applying the identity (4.17) yields To apply the operator Π K V to this expression we multiply the function with V , discard all terms in the Fourier series except e i k. x , and apply Π 1 2 with the help of the identity (4.17). The result is The next term, Π K V k.( P − k)Π 1 2 V |Ψ 0 , can be evaluated using a similar method. The calculation will make use of the identity in which e j are the standard basis vectors for R 3 and is an inward-pointing normal vector of unit length representing the normalised gradient of the nuclear charge density. Since k.( P − k)e i( k± a j ). x = ± k. a j e i( k± a j ). x , we obtain Now j k. a j ( n × e 3 ) l−j simplifies algebraically to π( n × k) l and Tr(σ l ψq)e i k. x = 2S l Ψ 0 , so The remaining two terms will be evaluated indirectly, using the identities In other words, we calculate the contributions from the sum of the spin 1 2 and spin 3 2 excited states and subtract the spin 1 2 contribution. We begin with Π K V 2 |Ψ 0 . The term in the Fourier series of (4.28) The other terms in the Fourier series will be annihilated by Π K , so need not be computed.
By the Clebsch-Gordan rules, . We only need to calculate the piece in H 0 0 ⊕ H 1 0 , because multiplying a spin 1 2 wavefunction with a spin 2 function yields wavefunctions with spin 3 2 and 5 2 , both of which will be annihilated by Π K . We show in the appendix that, for any vectors v, w ∈ R 3 , Therefore the relevant part of | i R ij (q)u i−j | 2 is | u| 2 + | u| 2 = W 2 0 + 2W 2 1 . It follows that and, using our earlier result (4.21), The term Π K V k.( P − k)V |Ψ 0 can be evaluated using similar techniques. The coefficient of e i k. x in the Fourier series of (V k. (4.32) As before, the other terms in the Fourier series are irrelevant. Also as before, we may replace | i R ij (q)u i−j | 2 with W 2 0 + 2W 2 1 . The resulting sum over j is zero, because j a j = 0. Therefore Π K V k.( P − k)V |Ψ 0 = 0 and, by our previous result (4.25), We are now in a position to evaluate the effective hamiltonian. Collecting together the results (4.13), (4.16), (4.21), (4.25), (4.31) and (4.33) gives This hamiltonian, which is analogous to equation (2.24) in the 2D model, contains the sought-after coupling between momentum and spin (1.1). Besides scalars, this is the only term in the hamiltonian, and it is at first sight surprising that no other terms occur. The explanation lies in the symmetries of the lattice: S. n × k is the only term linear in k which is invariant under the action of the binary cubic group. For the parameter set (3.11) the coefficient of the term (1.1) in H eff is negative, which is opposite to what would be expected based on the classical rolling motion of Skyrmions. This is not such a surprise, given what we learnt from the toy model. In the toy model, spin-momentum effects consistent with the classical rolling motion of Skyrmions only occurred for a relatively strong potential, and were inaccessible to perturbation theory. In the next section we investigate stronger potentials.

Strong coupling to the potential
In the previous section we discussed the situation where the potential is small; in this section we discuss the case where the potential is slightly larger. Recall that in the 2D toy model, if the potential was strong the lowest energy Bloch wave had a non-zero crystal wave vector (at r = 0 so k = ± 1 2 ). We expect a similar effect in the 3D model. We begin this section by looking for candidate crystal wave vectors for the ground state, using symmetry as a guide.
Recall that the hamiltonian is invariant under an action of the binary cubic group. The action of this group on wavefunctions induces an action on the space of crystal wavevectors k. The generator τ acts trivially on k, while the generators ρ and σ act on k as multiplication by the matrices The vectors are special because they represent fixed points of the action of the subgroup generated by ρ and τ , namely the binary tetrahedral group (bear in mind that k is only defined up to addition of the reciprocal lattice vectors a j ). These two crystal wavevectors are plausible candidates for the wavefunction of the ground state at strong coupling. Note that they are at the vertices of the first Brillouin zone, as shown in Figure 6. In order to analyse the Hilbert spaces corresponding to these crystal wavevectors it is convenient to apply a rotation to the lattice and the moving Skyrmion: After rotation, the Skyrmion moves in the plane z = 0, and the half-filled lattice of Skyrmions is the region z < 0. The generators of the binary cubic group now act as follows: After rotation the reciprocal lattice vectors are

Perturbation theory in k
We will be interested in eigenfunctions of the hamiltonian whose crystal wavevector is close to k ± . It is enough to analyse just wavevectors close to k + , as the transformation τ swaps k + and k − . First we will identify an orthonormal basis |Ψ 0a ∈ H k + for the eigenspace of the hamiltonian with (degenerate) lowest energy eigenvalue E 0 . Then we will consider nearby wavevectors k = k + + δ k. Perturbing k in this way is mathematically equivalent to perturbing the momentum operator: where P 0 = −i∇ x is the usual momentum operator acting on H k + . Thus nearby wavevectors can be analysed using perturbation theory. The perturbed hamiltonian is We will show below that Ψ 0a | P 0 |Ψ 0b = 0 for reasons of symmetry, so the effective hamiltonian acting on this eigenspace is unchanged to linear order in k. Therefore the perturbed wavefunctions satisfy H|Ψ a = E 0 |Ψ a + O(δ k 2 ). We then compute the matrix elements of the hamiltonian H to second order in δ k: For large enough M the second term on the right dominates the third term, meaning that the lowest energy eigenvalue has a stable local minimum at δ k = 0. Below we will quantify how large M needs to be for this to happen. The expectation value P of P = (P 1 , P 2 ) in the state |Ψ 0 is, as we have already noted, zero. Similarly, group theoretical arguments will show that the expectation value S of S = (S 1 , S 2 , S 3 ) has vanishing planar components (although the component perpendicular to the plane will be non-vanishing). For δ k = 0 we expect these expectation values to be non-zero and correlated. More precisely, we expect P 0 to point in the same direction as n × S , where n = (0, 0, −1) is now the normalised gradient of the nuclear matter density. Equivalently, where S ± := S 1 ± iS 2 and P ± 0 := P 1 0 ± iP 2 0 . It is straightforward to derive expressions for these expectation values within the framework of perturbation theory in δ k. The expectation value of P in a normalised We will show below that for sufficiently large M the second term is negligible and we have that P ≈ δ k. For S we compute This equation and (5.16) are analogues of eqs. (2.22) in the 2D model. In terms of κ = δk 1 + iδk 2 , we have that P + ≈ κ and to leading order. Thus to verify (5.15) it is sufficient to show that This concludes the outline of what we intend to show. In the remainder of this section we verify equations (5.19) and (5.20) by explicit calculation. In the next section we provide an alternative verification based mainly on symmetry.

Truncation of Hilbert space
In order to calculate the eigenstates |Ψ 0a we make a number of simplifying assumptions. First, we assume that the only terms that occur in the spatial Fourier series of Ψ 0a are those with the shortest possible wavevectors, namely Note that these all have the same crystal wavevector; for example, in the case of e 1 and e 2 this is because Second, we assume that the only terms that occur in the expansions of Ψ 0a in harmonics on SU(2) are those corresponding to spin 1 2 . In other words, for 2 × 2 matrices ψ 1a , ψ 2a , ψ 3a . Since these three matrices have altogether 12 degrees of freedom, the eigenstates |Ψ 0a belong to a 12-dimensional subspace of the Hilbert space H k + . These assumptions are justified as long as energies of states in the 12-dimensional subspace are appreciably lower than those in its complement. If the moment of inertia Λ is small then states with spin greater than 1 2 will have much greater energy than the spin 1 2 states considered here, so truncation to spin 1 2 can always be justified by choosing Λ small. To justify the truncation in momentum space, we need to consider the next-shortest wavevectors associated with k + . These are 2 3 ( a 3 − a 2 ), 2 3 ( a 1 − a 3 ) and 2 3 ( a 2 − a 1 ), and their associated kinetic energies are Later we will compare these with the energies of states in the 12-dimensional subspace. The generators r = ρ, τ of the binary tetrahedral group act naturally on wavefunctions H k + via r · Ψ( x, q) = Ψ(r −1 ( x, q)), and these actions fix the 12-dimensional subspace. However, they only define a projective representation and not a true representation, because The binary tetrahedral group is known to be Schur-trivial, meaning that every projective representation can be turned into a true representation by twisting the actions of the group elements. In this case, a true representation is obtained by choosing ρ · Ψ( x, q) = Ψ(ρ −1 ( x, q)) , τ · Ψ( x, q) = ωΨ(τ −1 ( x, q)) .

(5.26)
Here we have introduced ω = e 2πi/3 = − 1 2 + i √ 3 2 , the cube root of unity. We wish to break up the 12-dimensional subspace of the Hilbert space into irreducible subrepresentations of the binary tetrahedral group. To this end, we review these irreducible representations. Besides the trivial representation, there are two further 1-dimensional representations A a with a = 1, 2, given by The binary tetrahedral group can be identified with the subgroup of the group of unit quaternions generated by −1, ρ = − 1 2 (1+i+j+k), τ = i. The standard identification of unit quaternions with SU(2) matrices gives a two-dimensional representation E 3 . There are two further inequivalent representations Finally, there is a three-dimensional representation F given by R : SU(2) → SO (3).
It is straightforward to check that the action of the binary tetrahedral group on the span of e 1 , e 2 , e 3 ∈ H k + is isomorphic to the representation F . The action on the four-dimensional subspace of H 0 consisting of functions of the form Ψ( x, q) = Tr(ψq) is isomorphic to E 1 ⊕ E 2 . This can be seen as follows: the induced action on the 2 × 2 matrix ψ is The matrices acting on the left correspond to the representation A 1 ⊕ A 2 , and those acting on the right correspond to the representation E 3 , so the representation is The action on our 12-dimensional subspace is therefore F ⊗ (E 1 ⊕ E 2 ). This, it turns out, is isomorphic to 2E 3 ⊕ 2E 1 ⊕ 2E 2 . To fully describe the decomposition, we introduce basis vectors f ia , with i = 1, . . . 6 and a = 1, 2: It can be checked that f 1a span an irreducible subrepresentation isomorphic to E 3 , f 2a span a second copy of E 3 , f 3a and f 4a span two copies of E 1 , and f 5a and f 6a span two copies of E 2 .

Hamiltonian matrix and the ground state
Next we need the matrix elements (5.14) for the hamiltonian acting on our truncated Hilbert space. The non-trivial part is the potential. After rotation, the potential given by equations (3.8) and (3.10) becomes where U α = W α e iθα (α = 0, 1) as before and This acts on wavefunctions from our 12-dimensional space by multiplication, and, in order to have a well-defined action, the resulting functions need to be projected back onto the 12-dimensional space. Consider first the action of the functions e i a j . x with j = 1, 2, 3. In the case j = 1 we find that The first and third of these are orthogonal to e 1 , e 2 , e 3 so only the second of these survives projection onto the span of e 1 , e 2 , e 3 . By performing similar computations we find that the actions of e i a j . x are e i a j . x e k = δ j+1,k e k+1 (no sum over k) . (5.44) In this expression, indices i, j, k are to be understood modulo 3.
The effect of multiplying a wavefunction with R ij (q) and projecting back to the 12-dimensional space is described by the identity (4.17). Therefore the action of the functions Tr(A α (x)R(q)) that appear in the potential on the 12-dimensional subspace of the Hilbert space can be computed using equations (5.44) and (4.17), and turns out to be Tr(A α (x)R(q)) · ψ ib = B α;ji ψ jb , (5.45) where B α are 6 × 6 block diagonal matrices of the form The action of the potential function is therefore described by the 6 × 6 block diagonal matrix The latter was computed in the previous section using perturbation theory to be The value of W 2 0 + 2W 2 1 is approximately 1.06. Since we are only interested in energy differences we ignore the term 3/8Λ which occurs in both expressions. Since we have been assuming that Λ is small, the other Λ-dependent term in brackets can be ignored. Thus the state with crystal wavevector k + will have lower energy if This inequality holds for M in the range 4.04 < M < 10.11. Thus for M close to zero (equivalent to small potentials) the state with k = 0 is preferred, but as M increases past the value 4.04 the state with k = k + is preferred. Now we assess the reliability of the approximation that we made by truncating in momentum space. The largest eigenvalue of the 6 × 6 block diagonal matrix that describes the potential is 0.37. Thus the largest energy involved in our calculation is In our truncation of the Hilbert space we neglected states whose energy is bounded below by (5.24). We are justified in neglecting these provided that This means that our approximation is valid for the values of M around 4.04 where the transition between the k = 0 and k = k + states occurs.

Expectation values for spin and momentum
Now we turn our attention to the expectation value of spin and momentum in the ground state. We need to compute matrices describing the action of P 0 and S on the 12-dimensional subspace of the Hilbert space. The action of S is given by It follows that S 3 f ia = (−1) i+1 1 2 f ia for i = 1, . . . , 6 and a = 1, 2. In particular, for the ground states ψ 0a given by (5.48) we have S 3 ψ 0a = 1 2 (−µf 3a − νf 4a ), and the expectation value of S 3 is 1 2 (µ 2 − ν 2 ) < 0. The action of S + = S 1 + iS 2 is described by the 6 × 6 matrix The action of S − is given by the conjugate transpose of this matrix. As the blocks on the diagonal are zero, the expectation values of S 1 and S 2 in the states ψ 0a are zero, so the expected spin points vertically down into the half-filled lattice of Skyrmions, as was previously claimed. The action of P 0 = −i∇ on the functions e 1 , e 2 , e 3 is simply It follows that the action of P + 0 = P 1 0 + iP 2 0 is described by the 6 × 6 block diagonal matrix The action of P − 0 is given by the hermitian conjugate of this matrix. It follows that the expectation value of P in the ground state is zero as claimed.
Using these formulae it is straightforward to verify equations (5.19) and (5.20). We have that The block diagonal structure of the matrix representing H 0 means that the inner products on the left hand side of (5.20) vanish as required. Using these identities and our particular values for U 0 , U 1 , we find that Note that by construction T ij ab = T ji ab . It is straightforward to show using the matrix given earlier for P + 0 that Ψ 0a |P + 0 (H 0 − E 0 ) −1 P + 0 |Ψ 0b = 0 and Ψ 0a |P − 0 (H 0 − E 0 ) −1 P − 0 |Ψ 0b = 0. These two identities imply that T 11 ab = T 22 ab and T 12 ab = −T 21 ab . Altogether, this means that T ij ab is proportional to δ ij . The coefficient can be determined by evaluating So for M < 4.03, the expectation of momentum points in the opposite direction to δ k and for M > 4.03, the region of most interest, they point in the same direction. Notice that the transition occurs at almost exactly the same value of M as where the energy for k = k + drops below that for k = 0. Finally we consider the subleading corrections to the eigenvalue E of H implied by eq. (5.14). This equation can be rewritten in terms of T ij ab as follows: This concludes our verification of spin-momentum coupling based on the crystal wavevector k + . If M > 4.04 then the crystal wavevector k + is preferred over k = 0, and the expectation values of spin and momentum are correlated in the manner predicted by the spin-momentum coupling. Our calculation is reliable as long as M 8.9.

Symmetry arguments
To conclude, we would like to point out that our results in the previous section are robust and insensitive to the details of the choice of potential function. Many of them can be derived using symmetry alone, as we now explain. We begin by analysing the symmetry properties of the operators P and S. Their commutation relations with ρ are as follows: ρS 3 = S 3 ρ , ρS + = ω 2 S + ρ , ρP + = ω 2 P + ρ . (6.1) Since the hamiltonian commutes with the action of the binary tetrahedral group, the eigenspace corresponding to the lowest eigenvalue E 0 forms a representation K of this group. Generically this representation will be irreducible, as was the case in the above calculation. Since the group element −1 acts non-trivially on the Hilbert space, K must be isomorphic to one of the three representations E 3 , E 1 and E 2 introduced above, because −1 acts trivially in all other irreducible representations of the binary tetrahedral group. The commutation relations above show that the images of K under S + and P + are isomorphic to K ⊗ A 2 . Since tensoring with A 2 cyclically permutes the representations E 3 , E 1 and E 2 , these image representations are not isomorphic to K. It follows that they are orthogonal to K. This means that Ψ 0a |S + |Ψ 0b = 0 and Ψ 0a |P + 0 |Ψ 0b = 0 , (6.2) and in particular that S + and P + 0 have zero expectation value in the ground state. The identity (5.20) can be proved similarly. The operators S + (H 0 − E 0 ) −1 P + and P + (H 0 − E 0 ) −1 S + map K onto a representation isomorphic to K ⊗ A 1 , which is again not isomorphic to K, so the inner products in (5.20) have to vanish.
To analyse the identity (5.19) we need the symmetry σ. As has already been noted, σ maps the Hilbert space H k + onto H k − . There is another transformation which swaps k + and k − , namely time reversal T . This acts as 3) The composition σT maps H k + onto H k + . Its commutation relations with S and P are σT P 1 0 = −P 1 0 σT , σT P 2 0 = P 2 0 σT , σT S 1 = S 1 σT , σT S 2 = −S 2 σT . (6.4) Since multiplication with i anticommutes with σT , the transformation σT anticommutes with P ± and commutes with S ± .
The operator that appears in (5.19 When composed with projection onto the eigenspace K it defines a linear map K → K. This map commutes with the action of ρ and τ , so by Schur's lemma it acts as multiplication by a scalar. Since it anticommutes with the action of σT , this scalar must be pure imaginary. Thus symmetry arguments show that an identity similar to (5.19) must hold, with λ ∈ R. However, symmetry arguments alone cannot determine the sign of λ. This is because replacing the potential V with its negative −V changes the sign of λ without altering the symmetry properties. Nevertheless, the sign of λ does seem to be fixed by a few coarse features of the above calculation. Consider again the basis vectors ψ 0a for the lowest-energy eigenspace. Each of these can be written as a sum of three terms: Each summand is an eigenvector of P 0 , so has a definite momentum vector. Each summand also determines a unique spin vector v, such that it is an eigenstate of v. σ acting from the left with eigenvalue 1 2 . The momentum vectors and spin vectors for the summands involving e 1 , e 2 and e 3 are listed below: Note that for each summand, the momentum vector points in the opposite direction to the cross product of n with the spin vector.
The expectation values for momentum and spin are weighted averages of these vectors. In the case δ k = 0 the three summands contribute equally to the wavefunction, and weighted averages are ordinary averages. Since the momentum vectors sum to zero and the unweighted average of the spin vectors is 1 2 (ν 2 − µ 2 ) n, we recover the results derived earlier. When δ k = 0 the momentum eigenvalues get shifted by δ k and the dominant contribution to the wavefunction is from the summand with the shortest wavevector. For example, when δ k points in the direction − and n × S points in the direction of δ k. There are two effects contributing to the expectation value for P : the shift in momentum vectors and the change of weights. For strong potentials the former dominates, and the expectation value for P points in the same direction as the naive momentum δ k (see the discussion around eq. (5.16)). Thus n × S and P point in the same direction, consistent with the spin-momentum coupling.
Note that all of this follows from the correlation between the spin and momentum vectors of the three summands making up ψ 0a , and any vector similar to ψ 0a with µν > 0 would produce the same effect. Thus we expect a similar correlation between spin and momentum for all values of U 0 , U 1 close to those used in our calculation.

Conclusions and further work
In a classical picture, the experimentally observed nuclear spin-orbit coupling arises from a rolling motion of a nucleon over the surface of a larger nucleus. However, understanding why such a rolling motion is energetically preferred remains something of a mystery. We have shown here that for a Skyrmion close to the planar surface of a half-filled lattice of Skyrmions, a rolling motion is energetically favoured by the orientational part of the potential energy. To describe this planar rolling motion, it is convenient to introduce the notion of spin-momentum coupling.
We have next investigated the quantum mechanics of the Skyrmion, first by analysing the hamiltonian describing the Skyrmion interacting with the half-filled lattice of Skyrmions using perturbation theory. A spin-momentum coupling term appears at second order in perturbation theory, but has the wrong sign, at least for the parameter set obtained from the lightly bound Skyrme model. We then calculated spin-momentum coupling at the level of expectation values, and found that the correct sign is recovered non-perturbatively at stronger potential strengths. The change of sign is correlated with a jump in the crystal momentum of the lowest energy state.
Our results were based on a half-filled FCC lattice that has been sliced in the plane x + y + z = 0. There is another natural way to slice the FCC lattice, in a plane parallel to one of the coordinate planes (or x = 0, y = 0 or z = 0). It would be interesting to investigate the spin-momentum coupling in that situation.