About the solution of the numerical instability for topological solitons with long range interaction

The computations of solutions of the field equations in the Model of Topological Particles, formulated with a scalar SU(2)-field, have shown instabilities leading to discrepancies between the numerical and analytical solutions. We identify the origin of these deviations in misalignments of the rotational axes corresponding to the SU(2) elements. The system of a single soliton we use as an example to show that a constraint suppressing the wave-like disturbances is able to lead to excellent agreement between the result of the numerical minimisation procedure and the analytical solution.


Introduction
Topological solitons are interesting objects. Their masses are given by integrals over the energy density, the particle number is a topological quantum number and their interaction is a consequence of the topology. Well-known models of this type are the Sine-Gordon model and the Skyrme model. Both obey the laws of special relativity. The Sine-Gordon model [1] is a 1+1D model with 1 degree of freedom and two types of solitons differing by their chirality acting like a charge. They behave as expected for extended charged particles. Particles of equal charge repel and of opposite charge attract each other and annihilate. Bound states of particle pairs oscillate and are therefore dubbed breathing modes. Skyrme's model [2,3,4] is formulated in 3+1D with the three degrees of freedom of an SU(2) field with the interpretation of the meson field in nuclei. It was intended as a model for particles with the strong, short range interaction of nucleons. Skyrme's model can not model the Coulomb interaction by topological properties.
In [5] a model in 3+1D with the three degrees of freedom of an SO(3) field was suggested, modelling a long range Coulomb interaction for topologically stable solitons. It is of similar spirit as the Skyrme model, but due to the different Lagrangian it allows for three topological quantum numbers π 2 (S 2 ),π 3 (S 2 ) and π 3 (S 3 ). The solitons of this model can be interpreted as Dirac monopoles [6,7] without any singularity, without Dirac string and with a soft core. There are four types of stable solitons differing in the two topological quantum numbers π 3 (S 3 ) and π 2 (S 2 ) which can be interpreted as spin up and down, positive and negative charge. The equations of motion can be solved analytically for the one soliton systems. In addition, the model describes two types of Goldstone bosons, i.e. massless excitations propagating with the speed of light [8]. In Ref. [5] the model was published under the title "model of topological fermions" but possibly it should rather be dubbed model of topological particles since it turned out that also the Goldstone bosons of the model are characterised by a topological quantum number [9]. Due to the non-linearity of the model more complicated systems have to be solved numerically. But the numerics has suffered from numerical instabilities leading to large uncertainties [10,11,12]. In this article we report about a successful method to avoid these instabilities [13]. For the one soliton system we are able to compare in a careful analysis the numerical outcome with the exact analytical results. In this article we discuss the sources and the size of the errors of the numerical evaluations. In Sect. 2 we give a short overview of the model and the analytical solution for the one soliton system in Sect. 3. After discussing the cylindrical and the lattice formulation in Sects. 4 and 5 we demonstrate the failure of the calculation in Sect. 6. Further, we suggest in this section, to improve the calculations by a constraint. With the improved numerics we then get good agreement with the analytical results and give characteristic numbers for the achieved accuracy.
2 Formulation of the model SU(2) is the double covering group of SO (3). With the Rodriguez formula, SU(2) matrices can be expanded in a sin and a cos term, whereas the real 3 × 3 SO(3) rotational matrices in 3D need three terms. It is therefore simpler to do the calculations in SU(2) than in SO (3). For smooth field configurations in space-time there is the only difference that for every configuration of the SO(3)-field of the model we get two SU(2) configurations, differing by a multiplication with the non-trivial centre element of SU (2). Taking this into account, the three degrees of freedom of the model are formulated with a scalar field of SU(2) matrices of unit quaternions, where q σ = q i σ i is an element of the su(2) algebra, with the usual Pauli matrices σ i and Einstein's summation convention applied.
The Lagrangian of the model reads [5] where α f and r 0 are in principle arbitrary constants. Choosing for α f the value of Sommerfeld's fine-structure constant the force field of solitons can be compared with the Coulomb field and the size parameter r 0 can be adjusted to the mass of the lightest fundamental, charged particles existing in nature, to electrons. The dynamical term proportional to R µν R µν is equivalent to the Skyrme term of the Skyrme model. It can be formulated with the connection field on the SU(2) manifold Γ µ defined by R µν is an area density on the SU(2) manifold is satisfied. After a rotation of these local coordinates systems on the SU(2) manifold one recognizes in R µν the well-known form of the field strength tensor in QCD and the curvature tensor in general relativity The potential term proportional to q 2m 0 defines a two-fold degenerate vacuum at q 0 = 0. It fixes the size and mass of solitons.
We can relate the geometry to physics introducing a dual vector potential C µ and a dual field strength tensor * F µν by Since we are considering static cases, we have no magnetic fields, B i = 0 and get for the energy density, the 00 component of the energy-momentum tensor

Solitonic solution
For static monopoles at the origin [5] the model reduces to spherical symmetry in 3D with spherical coordinates r, ϑ, ϕ 1 where r denotes a vector in space and n a unit vector in the su(2) algebra. The Euler-Lagrange equation, a non-linear differential equation Due to the simplicity of this solution we will use further on the case m = 3. Corresponding to the three terms in Eq. (3.0.2) the radial energy density has three contributions a radial field from the contributions of R ϑϕ approaching at large distances the Coulomb field of a point charge, a tangential field from R rϑ and R rϕ and a potential contribution. They are depicted in Fig. 1. Integrating this radial energy density over ρ we get the energy for the monopole. Comparing this energy to the rest energy of the lightest fundamental monopole existing in nature 2 we can fix the radius r 0 , defining the size of the soliton [16] (3.0.5)

Electrodynamic limit
As long as we do not have analytical solutions for static systems with several charges or scattering problems, we are able to treat the fields outside the numerical integration region within classical electrodynamics only, a scenario which we dubbed in Ref. [8] electrodynamic limit.
In this limit we assume q 0 = 0 and neglect therefore tangential and potential energy contributions. To get estimates of the accuracy of the numerical calculations, it is sufficient to determine the size of these neglected contributions in spherical coordinates in the region ρ = r r0 > ρ > = r> r0 . Their contributions to the radial energy density are given by the integrals over the second and third term in Eq.
Dividing by the integral over the whole r-axis, α f c r0 π 8 , results in the relative error of the sum of these contributions In the further calculations we are using r > = 10 r 0 with an error of 0.167 %. From  see that for increasing r > the relative error decreases nicely with (r > /r 0 ) −3 . For comparison, in the same region the first term in Eq. (3.0.4) contributes with 12.6 %, which should not be neglected and is taken care of in the electrodynamic limit. These contributions we are going to determine in cylindrical coordinates.

Cylindrical formulation
The algorithm presented here [10,12,13], uses cylindrical coordinates to work with, because it should be capable of computing dipoles as well. Therefore, we will also do monopole calculations in cylindrical coordinates r, ϕ, z in order to get insight in the accuracy of the numerical calculations. The general soliton field (2.0.1) for any configuration in cylindrical coordinates reads The components of the curvature tensor (2.0.4) read and its squares (4.0.8) With the proper conversion factor to SI units and the length scales for cylindrical coordinates, l ϕ = r and l r = l z = 1 we adjust relation (2.0.7) between curvature tensor and electric field strength to cylindrical coordinates

Lattice computation
We consider the case of a monopole at the centre of a cylinder. Inside this cylinder we introduce the dimensionless coordinatesr,φ,z ∈ Z, defined by the relations r = ar, ϕ =φ, z = az,r ∈ {0, 1, . . . , n r },z ∈ {−n z , −n z + 1, . . . , n z } (5.0.1) who just number the lattice points. Due to the rotational symmetry around the z-axis we perform the ϕ-integrations analytically and are left with a two dimensional lattice with spacing a, often referred to as box in the following. It is characterized with the number n r of points in r direction and the number n z of points in ±z direction. Such a finite lattice suffers from boundary effects.
In the example of lattice QCD the boundary problems are often diminished by periodic boundary conditions. In the present model with solitons with long range Coulomb interaction periodic boundary conditions are in general not useful. Since we know the field configurations of charged particles from classical electrodynamics analytically, we use them for the boundary conditions.
We have to integrate | E| 2 over the whole volume except the cylinder, characterised by its radius R and its half length Z.
Due to the cylindrical symmetry, we get an integration factor of 2π and we are able to restrict our consideration to the rzplane, where the cylinder gets projected to a rectangular box. We split the remaining two dimensional area into three smaller regions, which are shown in figure 3. The energy outside the box may be written as (5.0.3)

Outside potential energy
After we did the appraisal for the involved energies in-and outside a sphere in Sect. 3, we want go on and find an analytical expression for the potential energy outside the cylindrical box, which we neglect in the electrodynamic limit. We integrate over H pot in Eq. (2.0.8) where we were splitting up the integral in the exact same way as done in (5.0.3) for the electrical energy outside the box. If we consider a lattice with the specifications of Z = 10 r 0 and R = 10 r 0 we get a relative error of 0.122 %. This value is a little smaller than the error 0.167 %, which we got for a spherical volume in Eq.

Discretization
Having gathered all energy contributions, we are now ready to go on with the discretization, which is necessary for the numerical calculations. The soliton field Q stated in (4.0.1), with its components q 0 , q r and q z is only defined on the sites of the rz-lattice.

Derivatives
To evaluate the curvature energy density H cur (4.0.10), we need the soliton field's derivatives of the form ∂ i q j in r and z directions only. Labelling the points in each of the two directions with integers as indicated in Eq. (5.0.1) we are using the five-point method for the first derivative This energy is computed on the lattice and minimised by a conjugate gradient descent algorithm.
Here, we want to emphasize again, that quantities with a bar over them are always dimensionless as can be seen in (5.1.5), because α f has no dimension, c ≈ 200 MeV fm and a measures the distance between two neighbouring points in fm. To enhance the accuracy of our calculations, we decided to do a cubic interpolation of the curvature-and potential energy density between the lattice points before integrating them.
Last, but not least we express the electrical energy outside the box H out el in terms of dimensionless units Summing up the various energy contributions, the total energy H tot , which will be minimised, reads

Results and accuracy
The presented algorithm essentially consists of two different parts, both carrying a numerical error. The first part is to calculate the various energy components and the second one is the energy minimisation to find the associated configuration. We realized, that we have to be very careful to avoid numerical instabilities, therefore we will take a look at both steps separately in the following.

Precision of the energy computation
The computation of the energy for a monopole with radiusr 0 in lattice units a results inH tot , the total energy H tot in units of  Fig. 4, where the two measures of accuracy are plotted versus different grid sizes. Blue lines indicate the total energy H tot and brown lines the ratio H tot /H pot . The analytical values are marked with dashed lines. As expected, for bigger lattices the calculations become more accurate. Below we will discuss the results after minimisation.

Precision of the conjugate gradient minimisation
Up to now, all we did was to place the centre of the analytical monopole configuration at the origin of the lattice and calculate its total energy. Carrying on, we want to test the minimising procedure 3 . If there were no numerical errors, the algorithm would stop after a few iterations. Every deviation from the original configuration may be seen as an error. In our case, the function to be minimised is the total energyH tot which depends on the field componentsH tot =H tot q 0 (r,z), q r (r,z), q z (r,z) with q 2 0 + q 2 r + q 2 z = 1.
The conjugate gradient method starts at a given point on the hypersurface ofH tot , which in our case is associated with the analytical configuration and iterates its way to a local minimum. After the difference of energies of two consecutive iterations falls below a certain value, the algorithm stops. It turns out that the procedure fails.
In Fig. 5 we compare the two-component vectors (q r , q z ) in the box before and after the minimisation and find clear differences between the numerical and analytical minima. Further, it can be clearly seen in Fig. 6, that after the minimisation the q 0 component showed up wave-like discontinuities. These deviations from the analytical solution are not physical and need to be avoided. From the analytical solution, we know that the field should be smooth. Figure 5: Shown are the vector components q = (q r , q z ), to allow for a comparison of the monopole configurations before (exact, blue) and after (with error, brown) the minimisation procedure. The magnification on the left illustrates the differences more clearly. Figure 6: q 0 -field in the rz-plane for the minimised monopole configuration. Observe the wave-like discontinuities along the radial direction.
The contribution of the discontinuities to the overall energy is very tiny and for bigger lattices they get negligible small. Nevertheless they are always present in the computations and are dangerous for the determination of static dipole configurations. Even fixing the monopole centres these small fluctuations lead to tedious collapses of the monopoles as was shown in several diploma theses [11,12,13].

Solution to the numerical instability
As shown in Fig. 5, the problem discussed in the previous section is accompanied by misalignments of the q-vectors. This fact suggests to diminish the variations in the directions of the q-field. Therefore, we have added a third term in Eq. (6.3.1) to the energy functional on the lattice in order to minimise the fluctuations of n := q/| q| and forcing discontinuities to vanish in the minimisation procedure. The energy functional which we minimise is thereforē where the order of magnitude is controlled by λ. This extra termH λ reminds us of the mass term for the pion in the Skyrme model. Likewise, the physical interpretation of our new term is, that it gives mass to photons and hence suppresses photonic excitations. Since we are computing static snapshots, the application of this method is justified. We are not interested in the size ofH λ , we want to use it as a constraint and to minimise among the set of lattice configurations with minimal H λ . In Fig. 7, we show a plot of H tot for simulations with different values for λ for 10 × 2 · 10 and 30 × 2 · 30 lattices. As one may expect, we see that the results do not change significantly any more, if the parameter λ exceeds a certain value. With this method the two-component vectors before and after the minimisation show no visible difference as Fig. 8 demonstrates at λ = 100. The correction of the minimisation procedure by the additional term leads to a nicely smooth distribution of q 0 -values around the centre of the monopole, see Fig. 9. Figure 7: Plot of the total energy H tot versus λ for lattices of size n r × 2n z with n r = n z = 10 and 30. Increasing λ, discontinuities get smoothed out in both cases very soon. By choosing a high enough value, they vanish completely and rising it further has no effect. The difference in the left and right scales is caused by the much smaller accuracy for n r = 10. Figure 8: Same simulation and parameters as in Fig. 5, but using in the minimization procedure the additional termH λ of Eq. (6.3.1) at λ = 100. As can be seen, the vector part of the configuration is nearly indistinguishable before and after the minimisation. Figure 9: Same simulation and parameters as in Fig. 6, but with the additional termH λ at λ = 100. Observe, that the waves in radial direction have vanished.
In Fig. 10 we compare the results of the different numerical evaluations of the total energy of the monopole for Z = 10 r 0 and R = 10 r 0 and various lattice sizes n r × 2n z and n r = n z . We can clearly see that deviations from the analytical results due to the lattice discretization. As discussed in Sect. 6.2, the minimisation without constraint leads to a strange shape of the field distribution but to a small effect on the energy only, as can be seen comparing the results before and after minimisation. With the additional constraint (6.3.1) the minimal configurations shows the expected behaviour and the expected energy. The results nicely agree before and after the minimisation with the constraint. Remaining differences to the energy of the analytical configuration are due to the chosen discretization of the energy functional and the missing terms for the potential and tangential energy outside the box. Adding the analytical expressions 2H out pot (5.0.5) = 0.312 keV for the sum of both terms would lead to a small overshooting of the analytical result. Figure 10: Total energy H tot of a monopole configuration with Z = 10 r 0 and R = 10 r 0 and different lattice sizes n r × 2n z and n r = n z , determined with different methods as explained in the text. The inset shows that the energy scale was chosen in order to show the tiny differences in the results clearly.
A similar comparison we show for the ratio H tot /H pot in Fig. 11. Due to the Hobart-Derrick theorem this value should approach 4. Figure 11: Same situations as Fig. 10. The ratio H tot /H pot is shown, which is 4 for the analytical minimum configuration.

Conclusions
Only the three field variables of a scalar SO(3)-field and an appropriate Lagrangian with two terms are necessary to describe the behaviour of topological solitons with long-range interaction. SU(2), the double covering group of SO (3), is well suited for investigations of the topology of field configurations. The known solutions of the field equations gives a nice opportunity to investigate the problems appearing in the numerical computations. Several attempts [11,12,13] to determine field distributions of static charges in two [12,13] and three [11] dimensions had failed due to numerical instabilities. We have shown in this article that the origin of these instabilities are misalignments of the rotational axes corresponding to the scalar SU(2)-field. Suppressing these wave-like disturbances with a corresponding constraint we were able to get results of the numerical minimisation procedure in excellent agreement with the analytical solutions. The application of this method to two soliton systems will in future be of special importance. It will allow to get accurate results for the interaction at short distances. Deviations from the Coulomb behaviour are expected and should be compared with the the running coupling, well-known from field theory.