Realistic classical binding energies in the $\omega$-Skyrme model

An omega-meson extension of the Skyrme model - without the Skyrme term but including the pion mass - first considered by Adkins and Nappi is studied in detail for baryon numbers 1 to 8. The static problem is reformulated as a constrained energy minimisation problem within a natural geometric framework and studied analytically on compact domains, and numerically on Euclidean space. Using a constrained second-order Newton flow algorithm, classical energy minimisers are constructed for various values of the omega-pion coupling. At high coupling, these Skyrmion solutions are qualitatively similar to the Skyrmions of the standard Skyrme model with massless pions. At sufficiently low coupling, they show similarities with those in the lightly bound Skyrme model: the Skyrmions of low baryon number dissociate into lightly bound clusters of distinct 1-Skyrmions, and the classical binding energies for baryon numbers 2 through 8 have realistic values.


Introduction
Skyrmions were used to model nuclei even before the birth of Quantum ChromoDynamics (QCD) [1]. The symmetries of hadronic physics at low energies were understood before QCD was an accepted theory of the strong interactions. In fact, QCD contains an extra U (1) symmetry compared to the low-energy chiral Lagrangian, and this caused scepticism until the so-called U (1)-problem was solved by 't Hooft [2]. As a consequence of Derrick's theorem [3], the topological solitons of the Skyrme model [1] need something more than the kinetic term to be stabilised. There is, however, little -if any -phenomenological support for adding the Skyrme term 1 . Starting from just the symmetries of the low-energy hadronic physics, it is possible to include just one more particle into the theory to stabilise the topological solitons, namely the omega vector meson. This was understood already in the seminal paper by Adkins and Nappi [4]. In a full-fledged hadronic physics model, several vector mesons would have to be incorporated. However, if the scope is simply the low-energy effective nuclear spectrum, perhaps a few -or just one -vector meson could be sufficient. The alternative option of including the rho meson instead of the omega meson was considered in a series of papers [5][6][7] and recently also by Sutcliffe and Naya [8][9][10][11]. 2 In the past 36 years, the omega vector meson extension of the chiral Lagrangian as a model for nuclei has not received much attention. Sutcliffe considered the model [17], but only constructed solutions of degree 1 to 4 within the rational map approximation, which approximates the field equations by ODEs. Recently, Speight considered the model with the addition of an (explicit) isospin symmetry breaking term in the form of a derivative coupling of the omega meson field and the pion field [18], but considered only the degree 1 sector where, again, only ODEs need be solved.
There is a good reason for this relative paucity of results: the static field equations in this model are not the Euler-Lagrange equations for the theory's static energy functional, so standard energy minimisation algorithms (based on gradient descent or simulated annealing) do not solve the static problem. The underlying cause for this difficulty is that the vector field representing the omega meson enters the Lagrangian with the "wrong sign." We overcome this obstacle by observing that static solutions solve a constrained energy minimisation problem in which ω 0 (the temporal component of the omega field) is uniquely determined by the Skyrme field. We solve this constrained energy minimisation problem by arrested Newton flow for the Skyrme field, updating ω 0 after each time step by solving the constraint equation. This equation is a linear inhomogeneous PDE which can be efficiently solved via a standard conjugate gradient method. The resulting algorithm, being based on a second order flow, is much faster than comparable heat-flow methods [19], allowing, for the first time, extensive simulation of a wide selection of topological sectors for a range of coupling values.
We find that this omega extended Skyrme model, although very simple -with only two parameters to dial -has regions in parameter space with extremely low binding energies. This addresses one of the usual problems with Skyrme-type models -that they are too strongly bound. In this model, we have a line of vanishing classical binding energy and beyond that even "negatively" bound solutions (that is, they are metastable 3 ). A where g is (without loss of generality) a positive coupling constant, vol M denotes the volume form on (M, η), and ·, · L 2 denotes the L 2 pseudo-inner-product on M defined by its Lorentzian metric (and the metric on N for the first term). The case of direct interest has M = R 1,3 (Minkowski space), N = S 3 , the unit sphere in R 4 and Ω the normalised volume form on S 3 (normalised so that S 3 Ω = 1). Note, in particular, that the baryon current in this formulation is the vector field on M metrically dual to the 1-form B = ϕ * Ω, where denotes the Hodge isomorphism on (M, η), and that this vector field is divergenceless by closure of Ω. We may identify ϕ = (ϕ 0 , ϕ 1 , ϕ 2 , ϕ 3 ) whose components are traditionally named σ = ϕ 0 and π i = ϕ i , i = 1, 2, 3 (the pions). A standard choice of potential is

2)
tunnel over to the global minimum in a finite time, which is exponentially prolonged by the barrier between the two minima.
which gives the pions mass m (in units of the omega mass). With these choices, the action (2.1) coincides with that introduced by Adkins and Nappi, in the normalisation used by Sutcliffe [17].
Returning to the general case, the field equations are obtained by demanding that (ϕ, ω) is a formal critical point of S: for all smooth variations (ϕ s , ω s ) of (ϕ, ω) = (ϕ s , ω s )| s=0 of compact support in M, To proceed further, it is convenient to choose an isometric embedding N ⊂ R k (such an embedding certainly exists; for N = S 3 we may choose the canonical embedding in R 4 ) and to associate to any smooth map ϕ : M → N the (d − 1)-form Ξ ϕ on M valued in ϕ −1 T N defined so that h(Y, Ξ ϕ (X 1 , X 2 , . . . , X d−1 )) = Ω(Y, dϕ(X 1 ), dϕ(X 2 ), . . . , dϕ(X d−1 )), (2.4) for all X 1 , . . . , X d−1 ∈ T p M and Y ∈ T ϕ(p) N . Recall that ϕ −1 T N is the vector bundle over M whose fibre over p ∈ M is the vector space T ϕ(p) N . This bundle will be of some significance in the following. A comprehensive description of it, and the geometric structures it canonically possesses, may be found in ref. [20].
Given a smooth variation (ϕ s , ω s ), we define ε = ∂ s ϕ s | s=0 and α = ∂ s ω s | s=0 . Note that ε is a section of ϕ −1 T N while α is a 1-form on M and that both, by assumption, have support in some compact set K ⊂ M. It follows immediately from eq. (2.1) and the Homotopy Lemma (see, for example ref. [21]) that d ds S(ϕ s , ω s ) where we have used Stokes's Theorem and the facts that, on a Lorentzian (d + 1)-manifold, the coderivative Ω p (M) → Ω p−1 (M) adjoint to d is (−1) p(d+1) d , and = (−1) d(p+1) [22]. This should vanish for all ε ∈ Γ(ϕ −1 T N ) and all α ∈ Ω 1 (M). Hence where P ϕ : R k → T ϕ N denotes 4 the orthogonal projection defined by the isometric embedding N ⊂ R k . These are the field equations for the action S. Note that each term on the left hand side of eq. (2.6), and hence the left-hand side itself, is a section of ϕ −1 T N .
So far, M was an arbitrary Lorentzian manifold. Henceforth, we assume that M = R × X with a product metric η = dt 2 − ζ, where (X, ζ) is a Riemannian d-manifold. We shall denote the Hodge isomorphism on X by * , to distinguish it from the isomorphism on M. Now ω = ω 0 dt + ω X where ω 0 and ω X are curves (parametrised by t) in Ω 0 (X) and Ω 1 (X) respectively. We shall denote by dω 0 and dω X the curves in Ω 1 (X) and Ω 2 (X) obtained by applying dΩ p (X) → Ω p+1 (X) at each fixed t, andω 0 = ∂ t ω 0 ∈ Ω 0 (X), ω X = ∂ t ω X ∈ Ω 1 (X). Similar conventions apply to dϕ andφ, having interpreted ϕ as a curve in C ∞ (X, N ). In this case, the theory enjoys time translation symmetry and hence, by Noether's Theorem, has a conserved energy functional where B 0 = * ϕ * Ω ∈ Ω 0 (X). Note that the quantity is a homotopy invariant of the map ϕ(t, .) : X → N since Ω is closed, and hence is independent of t. For suitable X and N it may be interpreted as the baryon number of the field ϕ.
Our aim is to find static solutions of the field equations, so let us assume that all fields are independent of t. Then ϕ = φ•π, where φ : X → N is a fixed map and π : R×X → X is projection. Furthermore, ϕ * Ω = B 0 dt = ( * φ * Ω)dt. Hence, eqs. (2.6), (2.7) are satisfied by ω = f dt and ϕ = φ • π, where f : X → R, provided ( + 1)f = −g * φ * Ω, (2.11) where is the usual 5 Laplacian on (X, ζ). This is the coupled pair of PDEs we seek to solve. Note that any solution of them has, by virtue of eq. (2.11) (and, if X is noncompact, a suitable decaying boundary condition on and hence energy The term Pϕ( d dϕ) is, up to sign, the tension field of the map ϕ. It can be defined without reference to an embedding N ⊂ R k using the natural connexion on the bundle ϕ −1 T N [20]. The extrinsic formulation is more convenient for our purposes. 5 We use the geometer's sign convention, so We claim that eq. (2.10) is precisely the Euler-Lagrange equation for the energy functional E(φ, f ) subject to the constraint (2.11). To verify this, let (φ s , f s ) be a smooth variation of (φ, f ) satisfying (2.11) for all s. Once again let ε = ∂ s φ s | s=0 ∈ Γ(φ −1 T N ) and α = ∂ s f s | s=0 ∈ C ∞ (X). Then, differentiating eq. (2.11) with respect to the variation parameter yields ( + 1)α = −g * d(φ * ι ε Ω), (2.14) and hence d ds where we have used eq. (2.14) in the last line. Now where, once again, decaying boundary conditions were imposed if X is noncompact. Hence Equation (2.17) has a useful reinterpretation. Given any smooth map φ : X → N , the constraint equation (2.11) uniquely determines the smooth function f : X → R, so we may think of E as a function C ∞ (X, N ) → R, that is, as a functional of φ only. Formally, C ∞ (X, N ) is an infinite dimensional manifold whose tangent space at a map φ is Γ(φ −1 T N ), the vector space of smooth sections of the bundle φ −1 T N . This space carries a natural inner product called the L 2 metric, so that, formally, C ∞ (X, N ) is a Riemannian manifold. In this picture, eq. (2.17) states that the gradient of the function E : C ∞ (X, N ) → R with respect to the L 2 metric is Note that this is, at each fixed φ, a section of φ −1 T N , and hence defines a vector field on C ∞ (X, N ).

Stability and the second variation formula
Since eq. (2.11) uniquely determines f for each given φ, we may interpret E(φ, f ) as a functional of φ only which, in a slight abuse of notation, we will denote E(φ). The static ω-Skyrme model thus defines a natural geometric variational problem for maps φ : (X, ζ) → (N, h) between Riemannian manifolds -to minimise E(φ) in a given homotopy class of maps -analogous to the classical harmonic map problem, where the energy to be extremised is simply the Dirichlet energy, Equation (2.10) is the condition for φ to be a critical point of E(φ), but its solutions are not necessarily local minima: they could be saddle points instead. To distinguish between minima and saddle points of E(φ) we must consider its second variation. The goal of this section is to compute and apply this second variation, exploiting the close analogy with the well-established setting of harmonic maps. To avoid technical issues with boundary conditions, we will assume throughout this section that (X, ζ) is closed.
We begin by briefly recalling the first and second variation formulae for E D (φ). Associated to any smooth map φ : (X, ζ) → (N, h) is a smooth section of φ −1 T N called the tension field, where {e i } is a local orthonormal frame on (X, ζ), and ∇ φ denotes the pullback of the Levi-Civita connexion ∇ N on T N to φ −1 T N . This is the natural connexion on φ −1 T N (recall, this vector bundle over X whose fibre above x ∈ X is T φ(x) N ), constructed from ∇ N . A thorough treatment of its definition and properties is presented in ref. [20, ch. 4]. In the extrinsic formulation used in section 2, τ (φ) = −P φ φ. Given a smooth one-parameter variation φ t of φ = φ 0 : X → N , with infinitesimal generator ε : so φ is a critical point of E D (a harmonic map) if and only if τ (φ) = 0. Consider now an arbitrary two-parameter variation φ s,t of a harmonic map φ = φ 0,0 , with infinitesimal generators ε := ∂ s φ s,t | s=t=0 and ε := ∂ t φ s,t | s=t=0 . Then is a certain second-order linear self-adjoint elliptic differential operator, constructed from ∇ φ and the curvature tensor R of (N, h), called the Jacobi operator. Explicitly, The second variation thus defines a symmetric bilinear form on Γ(φ −1 T N ) called the Hessian. We say that the harmonic map φ is stable if Hess D φ (ε, ε) ≥ 0 for all ε, and unstable otherwise. Determining the stability of a harmonic map thus reduces to a question about the eigenvalues of its Jacobi operator.
How does this generalise to our variational problem? We have already computed the first variation, (2.17), in the notation just introduced. To state the second variation formula requires two more preliminary definitions. First, given a smooth map φ : X → N , we define the linear first-order differential operatorΞ φ : for all x ∈ X, Y ∈ T φ(x) N , X 1 , . . . , X d−1 ∈ T x X. Second, given a smooth map φ : X → N , we define the linear integral operator α φ : Γ(φ −1 T N ) → C ∞ (X) which maps a section ε of φ −1 T N to the solution α of the linear PDE which exists and is smooth and unique by standard elliptic PDE theory. The linear operator α φ maps infinitesimal variations of φ to the corresponding infinitesimal variations of f . That is, given a variation φ t of φ, generated by ε = ∂ t φ t | t=0 , the corresponding variation f t of the solutions of eq. (2.11) has ∂ t f t | t=0 = gα φ (ε). We may now state the second variation formula (the proof, which is rather involved, is presented in Appendix A): Proposition 1 Let φ : X → N and f : X → R satisfy (2.10), (2.11). Let φ s,t be any smooth two-parameter variation of φ = φ 0,0 , f s,t be the corresponding variation of f = f 0,0 , preserving (2.11), ε = ∂ s φ s,t | s=t=0 and ε = ∂ t φ s,t | s=t=0 . Then In direct analogy with harmonic map theory, a critical point is stable if Hess φ (ε, ε) ≥ 0 for all ε, and unstable otherwise. Since α φ is not invertible, the stability question does not easily reduce to a spectral problem. Nonetheless, in an interesting family of special cases we can make significant progress.
Consider the case that (X, ζ) = (N, h), Ω is the volume form on (N, h), V = 0 and φ = Id, the identity map, that is, φ(x) = x. If N = S 3 , this is a simple model of dense nuclear matter with uniform baryon density, whose stability in the conventional Skyrme model was studied by Manton [23]. We will, for the time being, leave X = N general, however. It is well known that Id : X → X is harmonic, so τ (Id) = 0 [24]. Furthermore, * Id * Ω = * Ω = 1, since Ω was chosen to be the volume form. Hence the function f determined by eq. (2.11) is simply the constant function f = −g, so df = 0, and it follows immediately that φ = Id satisfies the Euler-Lagrange equation (2.10): Id is a critical point of E(φ) for all g. As we will see, the stability of Id depends, in general, on the coupling g, however.
The formula for the Hessian given by Proposition 1 simplifies radically in this case. First, since df = 0, the difficult terms involvingΞ Id and ∇ N Ω vanish (actually ∇ N Ω ≡ 0 since the volume form is parallel, so the latter term vanishes even for critical points with nonconstant f ). So, noting that V = 0, and it follows that if Id is stable as a harmonic map, it is also a stable critical point of E(φ). Hence, Id is stable for all g in dimensions d = 1, 2, or if (X, ζ) is Kähler, or if (X, ζ) is Ricci negative, for example [24]. If Id is unstable as a harmonic map (for example, if X = S d , d ≥ 3), things are more interesting: it is an unstable critical point of E(φ) for g ≥ 0 small, but may exhibit a stability transition, as g increases.
To proceed further, we note that the variation section ε is now a section of Id −1 T N ≡ T N ≡ T X, that is a vector field on (X, ζ), which greatly simplifies the Jacobi operator. In fact [24] J Id ε = ε − 2ρε, (3.11) where is the usual Hodge Laplacian on one-forms, is the metric isomorphism T X → T * X defined by ζ (i.e. ( ε)(u) := ζ(ε, u) for all u ∈ T x X), is its inverse, and ρ is the Ricci endomorphism of (X, ζ) (the linear map ρ : where Ric is the usual Ricci curvature tensor). Hence, where ε is an arbitrary smooth vector field on X. Every term in this, except the curvature term, − ε, ρε /2, is manifestly non-negative, so the question of stability of Id is nontrivial only if the Ricci curvature of (X, ζ) is positive somewhere. We shall prove that Hess Id is non-negative when evaluated on the subspace of divergenceless vector fields, and is, for large enough g, also non-negative on the subspace of pure gradients. From this, we can deduce that Id is stable, for g sufficiently large, if (X, ζ) is Einstein.
Proof: For any vector field ε on X, * Id * ι ε Ω = * ι ε Ω = divε, so α Id (ε) satisfies the PDE Hence, for all divergenceless vector fields ε 0 , α Id (ε 0 ) = 0. Further, by a formula of Bochner and Yano [24], where L denotes the Lie derivative, so for all divergenceless vector fields ε 0 , Lemma 3 There exists g 0 ≥ 0 such that, for all g ≥ g 0 and all smooth functions : X → R, Proof: Since X is compact, there exists a constant c > 0 such that, for all u ∈ T x X, Ric(u, u) ≤ cζ(u, u), and hence, for all vector fields ε, ε, ρε ≤ c ε 2 L 2 . Let 0 = λ 0 < λ 1 ≤ λ 2 ≤ λ 3 ≤ · · · be the eigenvalues of the Laplacian (on functions) on (X, ζ) and {f n } be a corresponding L 2 orthonormal basis of eigenfunctions, so f n = λ n f n . Since the sequence (λ n ) diverges to infinity, there exists q ∈ N such that, for all n > q, λ n ≥ 2c. Any function ∈ C ∞ (X) has a unique expansion = ∞ n=0 a n f n in the harmonics {f n }. Now λ n a n f n , Hence since λ 2 n ≥ 2cλ n for all n > q. If g is chosen so that that 4g 2 exceeds all the terms in this finite sum are non-negative, and the claim immediately follows.
Every smooth vector field ε on X uniquely decomposes into gradient and divergenceless components (just apply the Hodge decomposition to the one-form ε), and we have just shown that, for g sufficiently large, Hess Id is non-negative on both the gradient and divergenceless subspaces of Γ(T X). If Hess Id is diagonal with respect to the Hodge decomposition, it follows immediately that Id is stable for g sufficiently large. In particular: Proposition 4 Let (X, ζ) be a closed Einstein manifold. Then there exists g 0 ≥ 0 such that, for all g > g 0 , Id : (X, ζ) → (X, ζ) is a stable critical point of E(φ).
Proof: By Lemmas 2 and 3, there exists g 0 such that, for all g ≥ g 0 and all ε, (3.20) where ε = ε 0 + ∇ is the Hodge decomposition of ε into divergenceless and gradient parts (obtained by decomposing the one-form ε into coclosed and exact parts). Since (X, ζ) is Einstein, ρ = cId where c is a constant. Hence The claim immediately follows.
Proposition 4 covers, in particular, the case of most interest, X = S 3 . A careful recapitulation of the proof of Lemma 3 using the spectrum of the Laplacian for the unit d-sphere reveals that the critical coupling for X = S d , above which Id is stable, is We conclude by establishing a topological lower bound for E(φ). We now revert to the case of general (N, h), Ω and V while maintaining the assumption that X is compact and without boundary.
Proof: By the Cauchy-Schwartz inequality and eq. (2.11), and hence Note that this bound is quadratic in the topological invariant X φ * Ω. So, if N = S 3 and Ω is the (normalised) volume form on N , we see that the ω-Skyrme energy grows at least quadratically with the baryon number, E ≥ const × B 2 . This contrasts with the conventional Skyrme model, where the analogous bound on compact domains is E ≥ const × |B| 4/3 [25].

The numerical method
We seek to find, within a given topological sector, the minimum of E(φ) as defined in eq. (2.13), f being determined by φ using eq. (2.11). To do this, we choose an initial configuration φ(0) and solve Newton's equation for the motion of φ(t) in C ∞ (X, N ) subject to the potential function E : C ∞ (X, N ) → R, that is starting at rest,φ(0) = 0. In practice, we discretise space on a cubic grid and approximate grad E φ using finite differences, then use a 4th order Runge-Kutta scheme to perform the time stepping. This flow will start to roll "downhill", that is, reduce E, but will not, as it stands, converge to a minimum of E. To achieve this, we compare, after each time step, the energies of the new and old configurations. If E(t + δt) > E(t), we arrest the flow, restarting it withφ = 0. This strategy 6 is quite widely used in the study of topological solitons, but does not appear to have received a commonly accepted name. We propose to call it "arrested Newton flow".
In the present case, at each time step, to evaluate grad E φ (and E) we must construct the function f satisfying the constraint equation (2.11). This is a linear inhomogeneous PDE, or rather, having discretised space, a high-dimensional linear system of algebraic equations, so the obvious strategy is to use an off-the-peg linear algebra solver to compute f . This turns out to be inefficient, since such solvers are not iterative, in the sense that they start from scratch, making no use of an initial guess for the solution. For our application, after each time step, φ, and hence the right hand side of eq. (2.11), has changed very little, so we have access to an excellent approximation to f (t + δt), namely f (t). To exploit this feature, we reinterpret eq. (2.11) as the Euler-Lagrange equation for the quadratic functional which we solve by minimising G using an off-the-peg conjugate gradient method starting at f (t) (a particularly efficient choice for quadratic functions). The first application of this method (at t = 0), where we have only a rough guess for f (we use f = −gB 0 ) is quite computationally costly, but after each subsequent time step very few cycles of the conjugate gradient method (typically 0 to 3) are required to correct f to match the new Skyrme field φ to within the tolerance we require.
To illustrate our numerical scheme, we present classical energy minimisers of charges B = 1, 2, . . . , 8 for the coupling and pion mass proposed by Sutcliffe [17]: The calibration chosen by Sutcliffe fixes g by using the experimental value for the pion decay constant and the omega mass (hence fixing the length and energy units) and adjusting g to match the mass of the 4-Skyrmion to that of Helium-4. Fig. 1 shows coloured surfaces of constant baryon density for these solutions. The colouring represents the value of the normalised pion field π/|π| using a standard colouring of the unit sphere, which can be deduced from the picture for B = 1. The 1-Skyrmion is spherically symmetric, while the 2-Skyrmion is stable and has the shape of a torus as usual in the Skyrme-like models -this confirms the stability of the 2-Skyrmion which was an open question in the rational map approach with the same value of the coupling g [17]. The B = 3 topological sector contains the first metastable solution (local, but not global, energy minimiser), which is a baguette-shaped solution of three 1-Skyrmions stacked together horizontally (with the middle one flipped with respect to the outer two), see 3b in fig. 1 7 . It has, nevertheless, higher energy compared with the tetrahedrally symmetric "standard" 3-Skyrmion. The B = 4 Skyrmion is octahedrally symmetric and the B = 5 is dihedrally symmetric, as usual. The B = 6 sector contains a global minimiser of the energy functional with dihedral symmetry (which is the "standard" 6-Skyrmion solution) as well as a metastable solution; it can be interpreted as three 2-Skyrmions (tori) that are stacked on top of each other (with the middle one flipped); this is similar to how a cube is made of two tori (with one of them flipped), but just with an extra torus added in, see fig. 1(6b). In the B = 7 sector the energy functional is minimised by the icosahedrally symmetric Skyrmion as usual. Finally, the B = 8 topological sector contains three solutions. The stable solution is the dihedrally symmetric "standard" 8-Skyrmion with D 6d symmetry, unlike in the standard Skyrme model with a pion mass term (where the stable solution is composed of two B = 4 cubes). Additionally, here, there are two metastable solutions: the first and closest in energy to the minimiser of the energy functional in the B = 8 sector has a slightly smaller amount of symmetry, which we think is D 6 . The last and highest-energy solution in this sector is composed of two cubes, but unlike in the standard Skyrme model, they do not "melt" together, but merely attach to each other and hence look more like a regular crystal than the "standard" solution of the standard Skyrme model with a pion mass term does.
The first two B = 8 Skyrmions depicted in figs. 1(8) and 1(8b) are both approximately described by the rational map [27]: with z = e iϕ tan θ 2 being the coordinate on the Riemann sphere and a ∈ C. If a is real, there is an enhanced symmetry, i.e. D 6d , otherwise it is simply D 6 . The Skyrme field φ obtained by suspending this rational map is [27] where F is some (so far, unspecified) profile function. The standard Skyrme energy of this field depends on a only via which is minimised independently from F (r). The ω-Skyrme energy in the rational map approximation analogously depends only on a via I [17]. The minimum of I(a) is at a = 0.135 [27], but there is a saddle point at a = 0.101i. We think that in the ω-Skyrme theory, this saddle point has become a local minimum (and possibly moved a bit in the aplane). Thus we want to identify the stable and normal B = 8 D 6d symmetric Skyrmion of fig. 1(8) with a = 0.135 and the metastable (local minimum) B = 8 Skyrmion of fig. 1(8b) with a = 0.101i, which has D 6 symmetry.
We have searched extensively for a solution that looks like two cubes attached to each other with and without a twist along the axis that joins them (i.e. the global minimisers in the standard Skyrme model with a pion mass term), but have found -to our surprisethat they only exist as saddle points in the theory and decay into the dihedrally symmetric global minimiser (see the ancillary files for videos of this decay).
To summarise, all global energy minimisers for B = 1 to 8 turn out to have the same symmetries as the global minimisers in the standard Skyrme energy without a pion mass term. This model has a pion mass term and thus differs from the standard Skyrme model with massive pions in having a dihedrally D 6d symmetric fullerene-like B = 8 solution as the global minimiser of the energy functional.
Although the static solutions for the pion fields φ uniquely determine the corresponding omega meson functions f via the constraint (2.11), it will prove instructive to look at the difference between the baryon charge density B 0 and the function f . It is intuitively clear that the two quantities have some similarities and in particular, for large enough level set, they display surfaces of the same topology. Of course the difference between B 0 and f is due to the presence of the Laplace operator in the constraint equation which smooths out f in comparison with B 0 . In particular, this means that the "holes" -well known to reside in Skyrmion solutions -are filled up by said smoothing of the Laplace operator. This in turn has consequences for the energy density, which receives contributions from the omega meson field f and hence also is less "hollow" than the Skyrmion solutions in the standard Skyrme model. We conclude this section by comparing our solutions, obtained by solving the full PDE system, with the approximate solutions obtained by Sutcliffe [17]. These latter were obtained by using the rational map approximation for the pion field for B = 1 through B = 4, where the rational maps have spherical, axial, tetrahedral and cubic (octahedral) symmetries, respectively. The ω 0 = f field was obtained in ref. [17], by expanding it in symmetric harmonics, which are a linear combination of the usual spherical harmonics. The expansion was further truncated to angular quantum numbers l ≤ 10. Similarly, the baryon density was expanded in the same basis as the ω meson. This procedure led to at most 10 ODEs for the ω field and a single ODE for the pion fields. with the energies found in ref. [17] using the rational map approximation. For convenience, we also display the ratio of the energies with respect to that of the 1-Skyrmion.
In the usual Skyrme model without a pion mass term, the precision of the solutions obtained within the rational map approximation is surprisingly good, and the energies for B ≤ 22 are only about 1% higher than the energies of the true solutions (to the full PDEs), see ref. [27]. As can be seen in tab. 1, the accuracy of the rational map approximation is slightly worse in the ω-Skyrme model. Nevertheless, for B = 1, 2, 3, 4 the correct symmetries were predicted using the rational map approximation and their energies were at most 4.2% too large compared with the true solutions. Our results should therefore be regarded as a vindication of Sutcliffe's ingenious approximation.

Collective coordinate quantisation
The question remains, what value of the coupling g best reproduces the physical properties of atomic nuclei for low B? To answer this, we must calibrate the model (choose its length and energy units), and compare its data with experiment. For B = 1, particularly, quantum mechanical effects are an important component of these data, so we must devise a tractable quantisation scheme for our Skyrmions. The traditional approach is "rigid body quantisation", in which the action of the field theory is restricted to the spin-isospin orbit of a degree B classical energy minimiser. Recent studies of the standard Skyrme model suggest that this is, for B > 1, often too restrictive: the Skyrme field should instead be restricted (for each fixed t) to lie in some finite dimensional manifold M of configurations which includes the spin-isospin orbits of the global energy minimiser and all nearby local minima, and field configurations interpolating between these [28][29][30][31][32][33]. In general, determining M is a difficult task, more art than science at present. Note that by choosing M to be the spin-isospin orbit of the B-Skyrmion, we recover rigid body quantisation from the more general picture.
Let us assume that a finite dimensional manifold M of static degree B field configurations has been chosen, and that φ(t) moves slowly in M . As already observed, static fields produce no source for ω X = ω i dx i , so each point φ in M determines a function f = ω 0 , but induces no ω X . Once we allow φ(t) to move slowly in M , it produces a source for ω X of order |φ| so that, even in the approximation of low velocity, the terms in S involving ω X contribute significant terms to the Lagrangian determining slow dynamics in M . This subtlety was already apparent to Adkins and Nappi [4], although they do not give a detailed justification of their proposed resolution of it.
We propose the following procedure: for each φ ∈ M andφ ∈ T φ M , we take ω 0 and ω X to be the fields determined by eq. (2.7). We then substitute φ and ω into the Lagrangian defined by S (eq. (2.1)), keeping only terms up to quadratic order in time derivatives. This gives a Lagrangian L| governing the dynamics of a point moving in M (i.e. a slow curve of Skyrme fields) which can be quantised by standard methods. The Lagrangian defined by S of eq. (2.1) is where ·, · denotes L 2 inner product on X and · the associated norm, B X denotes the spatial part of the baryon current, B X = * φ * ιφΩ and which coincides with eq. (2.13) in the case where (φ, ω 0 ) is a static solution of the model. Assume now that ω X satisfies eq. (2.7). It follows immediately that the form ω = ω 0 dt+ω X is coclosed on M, and hence thatω where δ = (−1) p * d * denotes the coderivative of p-forms on X. Furthermore, the spatial component of eq. (2.7) reads Substituting eq. (5.5) into eq. (5.1) yields where, in the last line, we have used eq. (5.3) and discarded the irrelevant total time derivative.
In principle, formula (5.6) determines L|, the Lagrangian for motion in M . Given a curve φ(t) ∈ M , ω 0 (t) is determined at each time t by eq. (2.11), soω 0 is determined. We work to quadratic order in time derivatives and note that both B X and ω X are of linear order, so only the leading term in ω X is required. Henceω X may be discarded from eq. (5.4) which, given eq. (5.3) reduces to ( + 1)ω X = −gB X = −gφ * ιφΩ. (5.7) Then φ,φ uniquely determine ω X (by solving eq. (5.7)), so every term in L is determined by φ(t).

Quantizing the 1-Skyrmion
Let us apply this formalism to the motion of a B = 1 Skyrmion, where M is its spin-isospin orbit. Then E static is constant, and may be discarded from L|. Since the unit Skyrmion is a hedgehog field, rotation is equivalent to isorotation, and isorotation always leaves ω 0 fixed. Hence, for any curve in M ,ω 0 = 0, and so where ω X is determined by eq. (5.7). To proceed further, we must solve eq. (5.7) approximately. For this purpose we formally invert the operator 1 + yielding If we keep only the leading term, ω X ≈ −gB X , we obtain The curve φ(t) takes the form where E 1 , E 2 , E 3 is any orthonormal frame for T rn R 3 . Choosing E 1 = ∂ r , E 2 = Y /r, E 3 = n × Y /r where Y is a unit vector in T n S 2 one finds, after some routine algebra, . This covering map is an isometry, so the lifted metric is γ = Λ γ 0 where γ 0 is the round metric with radius 2 on SU (2) ≡ S 3 . The quantum Hamiltonian for geodesic flow is

Electric charge radius
The final phenomenological observable that we need is the electric charge radius. Computing this will require us to consider the Noether current associated with isospin symmetry, so it is convenient to revert to the Lorentz covariant setting in which the Skyrme field is regarded as a map on spacetime ϕ : M → SU (2) (rather than a curve φ(t) of maps X → SU (2)). Using the Gell-mann-Nishijima relation, the electric charge Q is given by where I 3 is the isospin and Y is the hypercharge which is given by where B is the baryon number and S is the strangeness quantum number. Since S = 0 for Skyrmions in SU (2) models (meaning 2 light flavors of quarks), we can write the electric charge density as We can construct the isospin density from the vectorial (Noether) current that is defined from the vectorial (isospin) transformation (as opposed to the axial transformation), whose infinitesimal form can be written as ϕ + α · ∆ϕ, (5.23) which in component form can be written as with i, j, k = 1, 2, 3 and α k being infinitesimal parameters and ∆ k the k-th isospin generator. The Noether current corresponding to the above infinitesimal transformation is given by the 1-form As usual with Noether currents, the time component contains the Noether charge, once integrated. The isospin charge density is thus proportional to which corresponds to the third isospin direction. We still have to find a proper normalisation of the current in order to use it for the electric charge density. Since we know that the nucleon with isospin ± 1 2 has electric charge 1 and 0, corresponding to the proton and the neutron, respectively, we can normalise the vectorial Noether current such that it integrates to ± 1 2 : The normalisation constant can thus be obtained simply as Using the baryon charge density B 0 = * ϕ * Ω and inserting the hedgehog Ansatz (5.12), we can finally write the electric charge density as Q ± = − sin 2 F (r)F (r) 4π 2 r 2 ± sin 2 F (r) + g 2 π −4 r −2 sin 4 F (r)F (r) 2 8π ∞ 0 r 2 sin 2 F (r) + g 2 π −4 sin 4 F (r)F (r) 2 dr , (5.30) which can readily be checked to integrate to 1 (0) for the upper (lower) sign, corresponding to the electric charge of the proton (neutron). We can now define the electric charge radius as the weighted integral

Calibration
An appealing point about the ω-Skyrme theory that we study in this paper is that it only contains 2 physical parameters: m ∈ (0, ∞) and g ∈ (0, ∞). m is physically the ratio of the pion mass to the omega meson mass m = mπ mω and g is a coupling constant β multiplied by the ratio of the omega meson mass and the pion decay constant g = βmω Fπ . β is related to the decay ω → 3π and is limited from above by experimental data [4]. The reason the data give only an upper bound on β is that the calculation of the ω decay to 3 pions in the model does not include the resonance ω → ρ + π (since the rho meson is absent from this theory), which enhances the decay rate. The upper bound calculated in Ref. [4] is β ≤ 25.4, whereas the same calculation with updated experimental data reads β ≤ 23.9, where we have used the decay width Γ(ω → 3π) 8.49 MeV, m ω 782.65 MeV F π 184.13 MeV [34]. In result, using the new data we get an upper bound for g ≤ 101.4, if we use the experimental values for m ω and F π . The energy units of the model are F 2 π mω and the length units are 1/m ω .
Physically, the pions are pseudo-Nambu-Goldstone bosons of chiral symmetry breaking in QCD and would be massless if the quarks were all massless. Nevertheless, this physical explanation for their small masses, puts an upper bound on m < 1. Furthermore, the Skyrmions tend to destabilise for m > 1. However, the limit m → 1 is theoretically interesting as it tends to unbind the Skyrmions and hence lower their mutual binding energies, which we shall see shortly. Using the experimental values for the meson masses, m = 0.176.
In the literature, two values of g have been used: g = 98.7 [4] and g = 34.7 [17]. The former value is found by letting F π and g be free parameters and fit the rotational excitation energy of the Skyrmion to the nucleon and Delta masses [4]. Fitting parameters to the Delta in Skyrme-type models, however, is filled with subtleties [35,36]. The latter value of g, on the other hand, is found by setting F π to its experimental value (186 MeV) and fitting the B = 4 Skyrmion mass to that of 4 He [17].

Fitting the nucleon and helium-4 masses
In this paper, we will consider the following calibration based on the idea that in a minimalistic model like the ω-Skyrme theory, we cannot accurately describe all phenomena of hadronic and nuclear physics with only 2 parameters over a large energy range. Hence, if we allow to fit the parameters of the model in order to fit baryonic quantities, disregarding the mesonic observables, then an appropriate list of quantities to fit the model with contains the masses of the nucleon and helium-4 as well as the size of the nucleon. The justification for doing so could either be that the model is incomplete or somewhat equivalently, that the parameters in the effective low-energy field theory have been renormalised.
The two equations for our calibration thus read fitting the mass of helium-4 to that of the 4-Skyrmion and fitting the mass of the nucleon to that of the 1-Skyrmion with the spin quantum correction (5.19), where m B is the static energy of the B-Skyrmion. Eq. (6.1) does not have a quantum correction from the spin, because the ground state of helium-4 is a spin 0, isospin 0 state. The factors of F 2 π /m ω and 1/m ω have been reinstated to convert to physical units (MeV). In principle, these two equations fix (F π , m ω ) in terms of m 1 (g, m), m 4 (g, m) and Λ(g, m). However, there is not always a solution, which we can see by taking the ratio of the two equations where we have used the definition of g. Substituting back into eqs. (6.1)-(6.2), we get There is always a solution if m 1 m 4 < m N m 4 He , however, we would additionally like the size of the nucleon to fit experimental data as well where c 197.3 fm/MeV and the radius of the nucleon as perceived by an electron in scattering experiments, is the electric charge radius given in eq. (5.31).    4 shows the omega mass m ω , the pion mass m π , the pion decay constant F π , the nucleon radius r N and the coupling constant β between the omega meson and the baryon current as functions of the dimensionless coupling constant g for various values of the mass ratio m. First we can see that this fitting procedure yields omega masses in the range ∼ (10, 90) MeV, which is between 1 and 2 orders of magnitude too small. The largest values of the omega mass tend to prefer small values of g. The pion masses are in the range ∼ (3, 80) MeV, which is also too small compared with data. The pion decay constant is in the range of ∼ (14, 95) MeV, which is not much worse than in many other Skyrme-like models, but still at least a factor of 2 too small compared with data. The nucleon radii are in the range ∼ (1.5, 41) fm, which is at least 71% too large compared with data; this is the Achilles heel of this fitting procedure. The coupling constant β is in the range ∼ (0.48, 91); the experimental upper bound is at about 23.7 and there are many solutions that obey this bound for g 33.
The biggest issue here is that the nucleon radius (electric charge radius) is at least 71% too large compared with experimental data.  . Solutions that fit to the nucleon mass and the helium-4 mass. The panels show the omega mass m ω , the pion mass m π , the pion decay constant F π , the nucleon radius r N and the coupling constant β. The figures for r N and β have been cropped so as to better see the viable content.

Fitting the nucleon radius and the helium-4 mass
In this subsection, we will fit the size of the nucleon and the mass of helium-4 to experimental data. The mismatch that naturally will happen now is that the nucleon mass will be larger than its experimental value. Fig. 5 shows the omega mass m ω , the pion mass m π , the pion decay constant F π , the nucleon mass including the spin quantum correction m N and finally the coupling constant β as functions of the dimensionless coupling constant g for various mass ratios m. The omega mass is generally too small in this fitting scheme,  Figure 5. Solutions that fit to the nucleon radius and the helium-4 mass. The panels show the omega mass m ω , the pion mass m π , the pion decay constant F π , the nucleon mass m N and the coupling constant β. The figure for m N has been cropped so as to better see the viable content. but for m 2.5 and large g, its experimental value can be reproduced, but at the price of the nucleon mass being more than 5 times heavier than it should be. The pion mass can be reproduced in this fitting procedure for g 100 for various mass ratios m < 0.9. The pion decay constant is generally larger in this fitting procedure than in the latter and is in the range ∼ (70, 145) MeV and hence always smaller than its experimental value. The nucleon mass is too large and in the range ∼ (1035, 5045) MeV. An issue is that the best values for the nucleon mass is just before the B = 4 Skyrmion becomes unstable; this is problematic because it is one of the most tightly bound Skyrmions. Finally, the coupling constant β is in the range ∼ (7,19) and hence is everywhere smaller than the upper bound from pion scattering.
Ideally we would choose a point in the model parameter space where the nucleon mass -including the spin quantum correction -fits experimental data. Since such a point is absent from the set of solutions, we could consider a less ambitious calibration scheme: we could continue to fit the 4-Skyrmion mass to that of helium-4 and the size of the 1-Skyrmion to that of the nucleon. If we set the classical mass ratio m 1 /m 4 ∼ 1/4, then we are in the right ballpark for a model with small binding energies -provided that the quantum corrections to each of the Skyrmions are roughly proportional to the topological degree. This choice corresponds to the green-dashed and the orange lines in fig. 3. Then the nucleon mass with the spin quantum correction is off and generally (always) too large compared with data. The justification of this lowering of ambition is that we do not really expect the spin quantisation to be the only quantum correction to the Skyrmion energies -especially in a regime where the binding energy is small [37]. The latter is due to the expectation of small binding energies yielding small vibrational frequencies [38]. Since there is no perfect data point (because the nucleon mass with the quantum spin correction is always too large), we will select a point in the parameter space as follows. We notice that although the minimum of the nucleon mass is around m ∼ 0. 25 We will present numerical solutions for g = 14.37, m = 0.176 in the next section. As we will see, they exhibit some striking differences from those obtained previously for Sutcliffe's coupling g = 34.7. (The situation for the Adkins-Nappi coupling g = 98.7 is rather similar to g = 34.7). For the multi-Skyrmion solutions, we begin the numerical calculations with initial configurations which are all made up of 1-Skyrmions placed in various random spatial patterns -generally rotated so as to attract each other. The existence of an attractive channel for m < 1 follows from a point source analysis whose details we postpone until next section. The numerical method described in sec. 4 then evolves the initial configuration using the arrested Newton flow until a local minimum of the energy functional has been obtained.  Fig. 7 shows the numerically obtained multi-Skyrmion solutions for B = 1 through B = 8. Obviously the B = 1 Skyrmion is a spherically symmetric solution. The first surprise is that the B = 2 and B = 3 solutions are delocalised bound states for the chosen calibration. Some insight into this phenomenon will be gained from a study of the inter-Skyrmion interaction energy. The obtained solutions are similar to those found in the pointparticle model [39,40] 8 . The remaining Skyrmion solutions with B = 4 through B = 8 are very similar to those found in sec. 4 for g = 34.7 (the Sutcliffe coupling), showing some universal features of the solutions. Briefly, the B = 4 Skyrmion has octahedral symmetry, the B = 5 Skyrmion has dihedral symmetry, the B = 6 Skyrmion has dihedral symmetry, the B = 6b Skyrmion is metastable and composed of three tori, the B = 7 Skyrmion has icosahedral symmetry, the B = 8 Skyrmion is D 6d symmetric whereas the B = 8b is only D 6 symmetric. Finally the B = 8c Skyrmion is similar to that of sec. 4, i.e. composed by two cubes sitting next to each other. However, for this value of the coupling, g = 14.34, the two cubes have repelled themselves to become a bound state of separated B = 4 cubes.
To summarise, the solutions for B = 2, 3 are like in the point-particle models, whereas the remaining solutions are qualitatively similar to solutions of the standard Skyrme model without pion mass.

Inter-Skyrmion forces
In this section we will compute the forces between widely separated 1-Skyrmions using a point-source formalism. This formalism was developed for the conventional massless Skyrme model by Schroers [42], and will require two modifications to deal with the omegameson version of the model studied here: the pion field is massive, and we must introduce point sources to replicate the Skyrmion's asymptotic ω field. Although the extra forces induced by this field are subleading if the pion to ω mass ratio is given its physical value, m = 0.176, it is instructive to include them, and to consider the (unphysical) regime where m ≈ 1, so, to begin with, we keep m general.
The starting point is to observe that the 1-Skyrmion takes hedgehog form where r = |x| and the profile functions F, f satisfy the coupled ODE system subject to the boundary conditions F (0) = π, f (0) = 0, F (∞) = 0, f (∞) = 0. Of particular interest is its asymptotic form for large r.
Since F, f are small at large r, we assume they are close to solutions of the linearisation of this ODE system about (F, f ) = (0, 0), from which we deduce that where p, q are some unknown constants which can be determined by solving the nonlinear system (8.2) numerically. The factors of 4π are introduced for later convenience.
The corresponding asymptotic pion and ω fields are These coincide precisely with the solution of the linearisation of our model about the vacuum φ = (σ, π) = (1, 0), ω µ = 0, in the presence of the external point sources Viewed from afar, therefore, the 1-Skyrmion looks like a point particle in a linear field theory consisting of three uncoupled scalar boson fields of mass m (the pions) and a single vector boson of mass 1 (the ω). This point particle carries three orthogonal scalar dipole moments p a = pe a , inducing the pion fields, and a vector monopole charge q inducing the ω 0 field. It has no vector current density (j i = 0) so (or rather, because) the point Skyrmion has no ω i field. This is the point Skyrmion in standard position (located at x = 0) and orientation. We may obtain the general point Skyrmion by translation and rotation (or isorotation, since these coincide within the hedgehog Ansatz).
Since the 1-Skyrmion is asymptotically indistinguishable from a point Skyrmion inducing fields in the linearised model (8.6), and physics should be model independent, we assume that the forces between well-separated 1-Skyrmions approach those between wellseparated point Skyrmions interacting via the Lagrangian (8.6), as their separation grows. Consider the case where the 1-Skyrmions are static and located at X (1) , X (2) and have been (iso)rotated through R (1) , R (2) ∈ SO(3) respectively. Then the corresponding sources are (α = 1, 2), which induce fields The interaction Lagrangian corresponding to this configuration of fields is lin is the Lagrangian density evaluated for field and source α, and L µ ) we find, after an integration by parts, that Let us define the relative position R and orientation O of Skyrmion 2 with respect to Skyrmion 1, Then the interaction potential, according to our point source model, is which can be written explicitly as , (8.14) where R ≡ R/R.
If m < 1 (for example, m = 0.176), the leading term in V int at large R is with corrections of order me −mR /R 2 . Hence, the two-Skyrmion interaction is maximally attractive if O R = − R, that is, O represents a rotation by π about some axis orthogonal to the line joining the two Skyrmions. This is the usual prediction of an attractive channel for appropriately oriented Skyrmions, leading to the expectation that Skyrmions can coalesce and form bound states. Note, however, that if m > 1, the uniformly repulsive interaction mediated by the ω mesons dominates at large separation, so we expect no bound states in this regime. The case m = 1 is interesting. Now the (potentially) attractive scalar dipole interaction and the repulsive vector monopole interaction have exactly equal range, and which one dominates depends on the relative sizes of the dipole moment p and monopole charge q. These quantities depend on the coupling g as well as the mass m, see Figure 9. In fact, for m = 1, p 2 /4 < q 2 for all 10 ≤ g ≤ 40, so vector repulsion dominates when m = 1 and we expect no bound states. Of course, the physical pion mass, m = 0.176, is rather far from this regime. Nonetheless, the fact that the vector monopole interaction is uniformly repulsive leads one to expect that binding energies in this model may be unexpectedly small, at least for some choices of g.
We will now perform a numerical calculation of the interaction potential in the full nonlinear model by sending two 1-Skyrmions towards each other in the attractive (meaning O R = − R, one of them is rotated by 180 degrees around an axis perpendicular to the line joining them) and the maximally repulsive channels (meaning O R = R, so one is a translated copy of the other). We treat the problem adiabatically and scatter the Skyrmions at small velocity compared to that of light. This way we can calculate the static energy functional at each step, neglecting the kinetic energy contribution. The final ingredient in this calculation is to track the position of the Skyrmions. We define the position of the 1-Skyrmion to be the position of the anti-vacuum, meaning φ 0 = −1. It is numerically difficult to be precise about this point using only φ 0 , which is why our scheme is based on finding the simultaneous zero in φ 1 = φ 2 = φ 3 = 0 for φ 0 < 0. The zero can be found by determining the sign change from one lattice point to another.  34. The product channel is made by translating a copy of one Skyrmion by R in some direction. The attractive channel takes the translated Skyrmion and rotates the it by π around and axis perpendicular to the axis separating them. The repulsive channel takes instead the translated Skyrmion and rotates it by π around the axis that separates them. The mass parameter is m = 0.176. Fig. 10 shows the result of the numerical calculation of the scattering potential. We display the scattering potential for two different values of the coupling g = 34.7 and g = 14.34.
In both cases, the repulsive channel displays a growth in the energy until it becomes difficult to continue the simulation adiabatically; at the point we stop the curve, one of the two Skyrmions either strays away or rotates into a different orientation.
For the attractive channel in the case of g = 34.7, the asymptotic energy corresponds to twice the energy of the 1-Skyrmion and as the separation is shortened, the total energy drops monotonically to the level of the 2-Skyrmion, which takes the shape of a torus.
For the attractive channel in the case of g = 14.34, on the other hand, asymptotically everything is similar. However, at short distances where the asymptotic approximation breaks down, the attraction (which is very weak for this value of the coupling g) is overcome by some nonlinear repulsion and the bound state is not a torus, but two 1-Skyrmions at a distance bound extremely weakly by their soliton tails. This is reflected in the classical energy minimisers for B = 2 and B = 3 for this coupling, which resemble lightly bound clusters of spherical 1-Skyrmions, rather than fully merged bound states.

Conclusion
In this paper we have studied the omega extension of the chiral Lagrangian, which gives stable topological solitons -known as Skyrmions -without the use of the Skyrme term. The stabilisation is provided by the interaction between the omega vector meson and the baryon current, which is a topological current -whose zeroth component measures the topological degree of the field.
Although the model has been discussed in one of the seminal papers by Adkins and Nappi, numerical solutions have not been obtained from the full PDEs -until now. Our method of solving the model entails rewriting the energy functional in terms of the pion field and a scalar (the 0-th component of the omega vector meson) field. In addition to this we implement a constraint equation that is itself also a PDE, but it is linear and can readily and quickly be solved by the use of the conjugate gradients method. We check the omega field at each time step in our code and improve it iteratively once it is needed. The pion field instead is evolved by means of a second-order method which we denote arrested Newton flow. In order to settle on a minimum of the energy functional, we remove the kinetic energy once in a while and every time that the potential energy increases.
Interestingly, we find that although the model only contains 2 parameters that we can dial, it has a large parameter space which includes a line with zero classical binding energy and even negatively bound metastable classical multi-Skyrmion solutions. This happens when the mass ratio parameter m is large (but still less than one) and the coupling to the omega meson is small (g 20). Due to the possibility of extremely lightly bound Skyrmions, there is in turn an emergence of a large number of metastable solutions (local minimisers of the energy functional) and hence a large potential for nuclear clustering in the model. The model at low coupling exhibits some similarities with the lightly bound Skyrme model studied by Harland et. al. [39]. These dissociated point-like Skyrmion solutions are also found in the Witten-Sakai-Sugimoto model [43,44] at strong 't Hooft coupling [45], see Ref. [41].
The ability to accommodate very low classical binding energies is a somewhat unexpected feature of the ω-Skyrme model. Another interesting feature is that the model can reproduce, in a very elementary manner, the mass splitting between protons and neutrons [18]. It would be interesting to see what effect the isospin symmetry breaking perturbation proposed in [18] has on the Skyrmions presented here.