Complex systems of Kuramoto–sine-Gordon solitons

The 1 + 1 dimensional Kuramoto–sine-Gordon system consists of a set of N nonlinear coupled equations for N scalar fields θ i , which constitute the nodes of a complex system. These scalar fields interact by means of Kuramoto nonlinearities over a network of connections determined by N(N − 1)/2 symmetric coupling coefficients a ij . This system, regarded as a chirally invariant quantum field theory, describes a single decoupled massless field together with N − 1 scalar boson excitations of nonzero mass depending on a ij , which propagate and interact over the network. For N = 2 the equations decouple into separate sine-Gordon and wave equations. The system allows an extensive array of soliton configurations which interpolate between the various minima of the 2π-periodic potential, including sine-Gordon solitons in both static and time-dependent form, as well as double sine-Gordon solitons which can be imbedded into the system for any N. The precise form of the stable soliton depends critically on the coupling coefficients a ij . We investigate specific configurations for N = 3 by classifying all possible potentials, and use the symmetries of the system to construct static solitons in both exact and numerical form.


Introduction
Emergent behaviours of complex networks in which dynamical systems, such as phase oscillators, are located at the nodes of a network of connections have been investigated for many years. The Kuramoto model [1], in particular, in which each node behaves as a simple harmonic oscillator with nonlinear network interactions, has been well-studied as a model of synchronization in which all nodes oscillate at a common frequency [2,3].
The concept of a complex network extends to systems of classical fields interacting over a network of connections, and further extends to systems of quantum fields. We can regard, for example, the standard model of particle physics as a complex system consisting of six flavours of quarks and six flavours of leptons, each of which can be viewed as a node of a network where the interactions, mediated by the electroweak force, define the connections of the network. An example of a complex system of classical fields, which has attracted considerable recent attention, is that of 'magnetic skyrmionics' [4][5][6] in which unit 3-vectors n i in 2 + 1 space-time dimensions are located at the nodes of a lattice, with spin-spin interactions restricted to nearest-neighbour connections. The configurations of interest here are skyrmions, described as 'magnetic whirls with a unique topology' [6], which are solitons constructed from the underlying fields n i , stabilized by a topological charge. Such complex systems display a wide range of configurations, even in static form, the existence of which is due to nonlinear interactions across the connections of the network, not necessarily as a consequence of self-interactions of the underlying fields.
In this paper we develop properties of the Kuramoto-sine-Gordon system (3), which is a 1 + 1dimensional relativistic complex system which arises from chiral field theory with Kuramoto-type interactions, and includes the sine-Gordon equation as a special case. Chiral field theory has been extensively studied as a low energy description of quantum chromodynamics, and therefore as a model for nucleons and other baryons [7]. The effective Lagrangian takes the form L = Tr ∂ μ g∂ μ g −1 where g is an element of a unitary group such as SU (2), and can be expanded to include higher-derivative terms. We extend this Lagrangian so as to construct a complex system of N chiral fields g i , which are n × n unitary matrices, with dynamics defined by the relativistic, chirally invariant Lagrangian density where c ij is any set of complex coefficients. The first N terms are the usual kinetic contributions, one for each node, and the remaining terms are the simplest chirally invariant interactions that can be constructed from the fields g i . These potential terms evidently do not contribute for i = j, hence there are no explicit self-interactions of the fundamental fields g i . L is invariant under the chiral (left-right) transformation g i → g L g i g R for each i = 1, . . . , N, where g L , g R are constant unitary matrices, and the Euler-Lagrange equations are: In space dimensions 2 and 3 we also add higher derivative terms so as to stabilize skyrmions, which are finite-energy classical solutions carrying a charge of topological origin [8][9][10]. In 2 space dimensions, for example, we obtain a system of interacting planar skyrmions which can be regarded as relativistic versions of magnetic skyrmions, except that in our case the coupling coefficients are not restricted to nearest-neighbour connections.
for any two maxima there can be many trajectories of which only one is stable. The question arises then as to how to construct stable, static solutions, if indeed any exist, for two given endpoints; (b) In the sine-Gordon model static solutions interpolate only between adjacent minima of the potential, i.e. there exist exactly two trajectories which originate from a given minimum, a soliton on one side and an antisoliton on the other. For two scalar fields (the case N = 3), we find that there are as many as 12 static trajectories originating from a given minimum to a nearby minimum, although their existence and stability depends on the coupling coefficients in ways to be determined. For N scalars there are many more possibilities. For time-dependent solutions the situation is even more complicated, for example there could exist exact multisoliton time-dependent configurations, originating from sine-Gordon solitons, which interpolate between any two minima.

Outline
We first summarize the main properties of the Kuramoto-sine-Gordon system in section 2, followed by a discussion of static solutions and properties of the superpotential W even though this function, although it exists, is not known explicitly. Then in section 3 we show that for any N the system describes N − 1 scalars of nonzero mass μ i together with a decoupled massless field. In section 4 we specialize to N = 3, in particular we discuss symmetries of the static equations and classify all possible potentials V for any coefficients a ij . We show how to obtain reduced equations with exact solutions, together with a detailed discussion of the double sine-Gordon equation. A particular point of interest is that this equation allows static configurations which represent two separate solitons located at arbitrary points on the line. Whether this property extends to larger N with any number of static solitons remains to be determined. The simplest case to analyse is that of full symmetry, in which a ij takes the value either 1 or −1 for all links (i, j), as we discuss in section 5. Stability of solutions, whether in exact or numerical form, is analysed in section 6, in particular we find an exact critical value of the ratio a 13 a 12 at which the sine-Gordon soliton becomes unstable. We discuss properties of numerical solutions by means of two examples in section 7, followed by general remarks for any N in section 8, and concluding observations in section 9.

The Kuramoto-sine-Gordon system
We specialize now to the system (2) for the case g i ∈ U(1) in 1 + 1 dimensions, and refer to [8] for properties of chiral systems in higher dimensions. With the parametrization g i = e iθ i and with real, symmetric c ij = a ij , (2) reduces to the Kuramoto-sine-Gordon system: and the Lagrangian (1) can be written in terms of the N scalar fields θ i as where for a ij 0 the potential V is given by which depends only on the differences θ j − θ i . The diagonal values a ii do not appear in these equations, meaning that there are no self-couplings, and so we have N(N − 1)/2 independent coefficients a ij . We can summarize the main properties of (3) as follows: (a) We have ∂ μ ∂ μ i θ i = 0 by antisymmetry of the summand, and the sum i θ i also decouples from the remaining independent fields in the Lagrangian (4). Since (3) is invariant under θ i → θ i + θ wave , for all i = 1, . . . , N, where θ wave satisfies the wave equation, ∂ μ ∂ μ θ wave = 0, we can without loss of generality choose i θ i = 0. The independent unknowns are constructed from the differences θ j − θ i , for example we can define the N − 1 independent fields ϕ i = θ i+1 − θ i . A special case of the transformation θ i → θ i + θ wave is chiral invariance under global U(1) transformations, for which θ i → θ i + θ 0 for any constant angle θ 0 . (b) (3) is invariant under 2π translations of the angles θ i due to the periodicity of V. Permutational symmetry arises from the action of the symmetric group S N on the N nodes, for example (3) is invariant under the (12) transposition θ 1 ↔ θ 2 provided that a 1i ↔ a 2i for all i = 3, . . . , N, and similarly for any two nodes i, j. (c) Because V is bounded above and below, we can always choose an additive constant such that V 0 for any combination of signs for a ij . For a ij 0 the vacuum states (zeroes of V) are constant fields θ i such that θ j − θ i = 2πn ij for integers n ij , for any i, j. Finite energy solutions of (3) require thatθ i , θ i → 0, together with cos(θ j − θ i ) → 1 as |x| → ∞ for all i, j. (d) For N = 2 the system (3) reduces to ∂ μ ∂ μ θ 1 = −∂ μ ∂ μ θ 2 = a 12 sin(θ 2 − θ 1 ), with the general solution where θ sG satisfies the sine-Gordon equation ∂ μ ∂ μ θ sG = −2a 12 sin θ sG . This equation has been wellstudied, see for example [22][23][24], also [9] (section 5.3) and [10] (chapter 1), in particular sine-Gordon solitons are known to correspond to particles in the quantum theory. Hence, for N = 2 there is a rich array of exact solutions consisting of single and multiple soliton and antisoliton configurations, together with a periodic breather solution. The mass (energy) of a single soliton is nonzero and proportional to √ 2a 12 , even though (4) contains no explicit mass terms. The soliton is stabilized by the network connections in the Lagrangian, not by self-interactions of the underlying fields. The well-known integrability properties of the sine-Gordon equation [23] do not extend to the system (3) for N 3.

Static solutions
Non-trivial time-dependent solutions exist for the sine-Gordon equation and have also been studied numerically for the double sine-Gordon equation, both of which appear as special cases of (3), as we discuss explicitly for N = 3 in section 4. However, let us restrict our considerations now to static solutions of (3), which by means of a Lorentz boost lead to travelling wave solutions, and satisfy from which we obtain the first integral where we have used the finite-energy boundary conditions to fix the constant of integration. The energy density for static solutions is given by and so the total energy is equipartitioned between the kinetic and potential contributions.

The superpotential
We can bound the total energy below in the usual way by means of a superpotential W to obtain first-order Bogomolny equations, as explained in [9] (sections 1.3 and 5.1), however this is of limited use because an explicit superpotential is not known for general coefficients a ij . Define W = W(θ 1 , . . . , θ N ) according to the first-order partial differential equation where W θ i ≡ ∂W/∂θ i and V is given by (5), then we can arrange the total energy E in the usual way as follows: For any solution θ i (x) of (7), denote W(x) = W(θ 1 (x), . . . , θ N (x)) then W (x) = i θ i W θ i and so we obtain the bound E |W(∞) − W(−∞)|, which is attained provided that the N first-order equations θ i = ± W θ i are satisfied. Solutions satisfy the static equation (7) and are stable against all perturbations which maintain the finite-energy boundary conditions. Explicit solutions of (10) are not known for N 3, for general coefficients a ij . For N = 2 we have W(θ 1 , θ 2 ) = √ 8a 12 cos 1 2 (θ 2 − θ 1 ), which satisfies W θ 1 + W θ 2 = 0 as well as (10) and, in general, W is a function of the differences θ j − θ i . Equation (10) is well-known as the eikonal equation which appears in geometric optics, and can be written as |∇u| 2 = f where the unknown u is defined on R N , or on an open set Ω ⊂ R N together with boundary conditions on ∂Ω, where f is a given function. The eikonal equation is a special case of the Hamilton-Jacobi equation which has been extensively investigated [25]. It is known that u exists for continuous functions f and, with extra conditions, is unique. Various existence theorems are proved in [26], see in particular section 4.1 for results with Ω = R N , also Theorem 5.1 and Remark 5.5 for the eikonal equation on R N . A smooth solution W exists therefore for the potential V in (5), with corresponding stable static solutions satisfying θ i = ± W θ i , which solve (7) and satisfy the finite-energy boundary conditions.
Multi-scalar soliton models are known for which it is possible to find an explicit superpotential, see for example [17] (section 2), and there are also methods for generating models with explicit superpotentials [27], but in our case, lacking an explicit superpotential, we proceed by solving the second-order equation (7) directly, either exactly or numerically. It is then necessary to determine whether this solution is stable against linear perturbations. Generally there exist multiple static solutions of (7) for any fixed boundary conditions, of which only that with the lowest total energy is stable, this being the configuration that would be derived from the superpotential. We consider stability of static solutions in section 6.

Elementary excitations of the Kuramoto-sine-Gordon system
For N > 2 explicit solutions of the Kuramoto-sine-Gordon system (3) exist only in special cases for restricted values of the coefficients a ij . We can, however, determine asymptotic solutions which correspond to the elementary excitations of the system, which in the quantum theory describe particles of mass μ i depending on the coefficients a ij , and which propagate and interact over the complex system by means of the potential V. These excitations appear in the usual way as perturbations of the vacuum states.
These same perturbative solutions are also useful in order to solve the static equation (7) for large |x|. Numerically we can solve (7) only over a restricted interval x ∈ [−L, L] for some fixed, large L, and we match these asymptotic solutions to the numerical solutions at x = L in order to generate solutions over R. We find that the angles θ i approach their asymptotic values with an exponential fall-off, from which we can identify the masses μ i of the fundamental excitations as functions of the coupling coefficients a ij .
If we consider only the case a ij > 0 for all nodes i, j, then the potential V is given by (5), and finite energy static solutions satisfy cos(θ j − θ i ) → 1 as |x| → ∞. The vacuum states are given by θ i = 2πn i for integers n i and so we linearize (3) by setting θ i = 2πn i + ε i , to obtain In matrix form we have ∂ μ ∂ μ ψ = −Mψ, where ψ = (ε 1 , ε 2 , . . . , ε N ) T and the N × N symmetric matrix M = (m ij ) is given by Since M is symmetric it can be diagonalized and has real eigenvalues. M has a zero eigenvalue corresponding to the eigenfunction (1, 1, . . . , 1) T which arises because the sum of the matrix elements of M in any row is zero, i.e. j m ij = 0. It follows from (12) that ∂ μ ∂ μ ( i ε i ) = 0, corresponding to the property that i θ i describes massless excitations of the system. All other eigenvalues of M are strictly positive, however, since M is a positive-semidefinite matrix. Let u = (u 1 , u 2 , . . . , u N ) T be any vector in R N then: The summand (u i − u j ) 2 a ij is non-negative, since we have assumed that a ij > 0, and so i,j (u i − u j ) 2 a ij is zero only if u i = u j for all i, j, corresponding to the eigenfunction (1, 1, . . . , 1) T of M with a zero eigenvalue. Otherwise u T Mu > 0 and so all other eigenvalues of M are strictly positive. Denote these eigenvalues by μ 2 i , then we can solve ∂ μ ∂ μ ψ = −Mψ by diagonalizing M to obtain N − 1 independent Klein-Gordon equations ∂ μ ∂ μ ψ i = −μ 2 i ψ i for functions ψ i . Hence the elementary excitations of (3) consist of N − 1 scalar fields of mass μ i > 0, together with a single decoupled massless field.
The two nonzero eigenvalues μ 2 1 , μ 2 2 of M are given by corresponding to two scalar particles of mass μ 1 , μ 2 . Properties of the elementary excitations enable us to determine asymptotic solutions of the static equation (7) for the differences θ i − θ j . These functions attain their asymptotic values with an exponential fall-off determined by the masses μ i , i.e. the asymptotic solutions are linear combinations of the functions e −μ i x . Hence we can deduce the appropriate boundary conditions to be imposed on the differences θ i − θ j at x = L, for large L, when finding numerical solutions on [0, L].

Static solutions for N = 3
We now consider the static equation (7) for N = 3 for arbitrary coefficients a 12 , a 23 , a 13 . The sum θ 1 + θ 2 + θ 3 decouples from the system, and we define the differences (7) reads which we write as where the potential U is given by If a 12 , a 23 , a 13 0, for example, we set c = a 12 + a 23 + a 13 and the zeroes of U occur at cos ϕ 1 = 1 = cos ϕ 2 . The first integral of (18), corresponding to (8), is and the energy density E is Although we can write E in diagonal form in terms of ψ ± = 1 2 (ϕ 1 ± ϕ 2 ), for example, we leave E as shown so that the permutational symmetries are more evident.
Static solutions of (18) are determined firstly by specifying the endpoints, i.e. the asymptotic values of ϕ 1 , ϕ 2 , which must be adjacent zeroes of U, where the meaning of 'adjacent' has yet to be determined. Each solution can be viewed as a trajectory traced out by the 3-vector (ϕ 1 (x), ϕ 2 (x), −U[ϕ 1 (x), ϕ 2 (x)]) on the surface of −U as x varies. This generalizes a property of solitons for a single scalar field, see for example Coleman [21], chapter 6, where solutions are described as particle motions in the potential −U, with x regarded as the time variable. In our case, the trajectory traverses the ridges and valleys of −U as it interpolates between the maxima of −U, to form a stationary point of the energy.
As an example, we have plotted the surface U defined by (19) for the symmetric case a 12 = a 13 = a 23 = 1 (with c = 3) in figure 1, where the zeroes of U are the points 2π(m, n) for integers m, n in the (ϕ 1 , ϕ 2 ) plane. U takes its maximum value at the points marked in red, which are also the minima of U for a 12 = a 13 = a 23 = −1, since for this case U simply changes sign, up to an additive constant. Solitons for a 12 = a 13 = a 23 = −1 interpolate between these maximum values, for example there is a soliton (shown in blue) which interpolates between (− 2π 3 , − 2π 3 ) and ( 2π 3 , 2π 3 ), passing through the origin, and a second (in green) interpolating between (− 2π 3 , 4π 3 ) and ( 2π 3 , 2π 3 ), of which only the second is stable, see section 6. We discuss solitons for the fully symmetric case in section 5, where we conclude that a static soliton joining the points (− 2π 3 , − 2π 3 ) and (− 2π 3 , 4π 3 ), for example, does not exist.

Symmetries of the static equations
The unbroken symmetries of the Kuramoto-sine-Gordon model for the static system (18) include translational invariance in which x → x − x 0 for any constant x 0 , and parity inversion x → −x, as well as the periodic and permutational symmetries, although not in manifest form, since we have defined preferred fields ϕ 1 , ϕ 2 . We can summarize the main invariant transformations of (18) as follows, where the first four relate to periodicity properties of ϕ 1 , ϕ 2 and the next five describe permutational symmetries: Many of these symmetries follow from combinations of others, for example (d) follows from (b) and (c), similarly for the periodic property (a), and the permutations (h) and (i) follow from the transpositions. We have stated symmetries such as (k) and (l) separately for future reference.
It is therefore not necessary to solve (18) for all possible values of a ij since, having found a solution for certain coefficients, others follow by symmetry. We can also use these symmetries to reduce (18) to a single equation, for example (k) is satisfied by setting ϕ 1 = ϕ 2 , provided that a 23 = a 12 , as discussed in section 4.3.

Minima of the potential
In order to determine the vacuum states and hence the boundary values for static solutions, we need to find the minimum values of U for general coefficients a ij . Stationary points of U occur when the right-hand side of (18) is zero, which also determines the fixed points, then: We have either sin where this fixed point exists only for coefficients a 12 , a 13 , a 23 such that the right-hand side in absolute value is less than or equal to unity. It follows from (24) that a 2 12 sin 2 ϕ 1 = a 2 23 sin 2 ϕ 2 , and so the relative signs of sin ϕ 1 , sin ϕ 2 must agree so as to satisfy a 12 sin ϕ 1 = a 23 sin ϕ 2 .
There are four possibilities for (23), namely cos ϕ 1 = ±1, cos ϕ 2 = ±1, however it is sufficient to consider only the case cos ϕ 1 = 1 = cos ϕ 2 , since the others follow from the symmetries (b), (c) and (d). In order to determine whether U attains a minimum at cos ϕ 1 = 1 = cos ϕ 2 we first calculate the matrix U (2) = (U ϕ i ,ϕ j ) of second derivatives, given by: then at cos ϕ 1 = 1 = cos ϕ 2 we have det U (2) = s 2 and U ϕ 1 ,ϕ 2 = a 12 + a 13 , together with U = c − s 1 , where s 1 , s 2 are defined in (16). The point cos ϕ 1 = 1 = cos ϕ 2 is therefore a minimum of U provided that s 2 > 0, a 12 + a 13 > 0, in which case we set c = s 1 in (19). Provided that s 2 > 0, it is sufficient that any one of the three parameters a 12 + a 13 , a 13 + a 23 , a 12 + a 23 be positive, then the others are also positive, and cos ϕ 1 = 1 = cos ϕ 2 is a minimum of U for all such cases. We can include here the special case in which one of a 12 , a 13 , a 23 is zero.
If the fixed point (24) exists, then at this point we have U ϕ 1 ,ϕ 1 = −a 12 a 13 /a 23 with det U (2) = a 12 a 23 sin ϕ 1 sin ϕ 2 = a 2 12 sin 2 ϕ 1 > 0, hence (24) is a minimum of U provided that s 3 < 0, where s 3 is defined in (16). In this case we choose then U 0 with U = 0 only at the points (24). We have found therefore five classes of potentials U in (19), corresponding to the four stationary points (23) which are related by symmetries, and (24), at which U attains a minimum. The finite-energy requirement of solitons means that static solutions must attain one of these minima as |x| → ∞.
As an example, let us consider the fully symmetric case in which the coefficients a ij are all equal and so by rescaling take the value either +1 or −1. For a 12 = a 13 = a 23 = 1 we have from (19) U(ϕ 1 , ϕ 2 ) = 3 − cos ϕ 1 − cos ϕ 2 − cos(ϕ 1 + ϕ 2 ), which has minima at the points 2π(n, m) in the (ϕ 1 , ϕ 2 ) plane, for integers n, m, at which cos ϕ 1 = 1 = cos ϕ 2 . Evidently these minima form a square lattice array in the (ϕ 1 , ϕ 2 ) plane. Figure 1 plots U as a function of ϕ 1 , ϕ 2 , where the periodicity of U is evident. The points at which U attains its maximum value of 9 2 are marked in red. A contour plot of U is also shown in figure 3 (left), showing the square lattice array of minima.
If a 12 = a 13 = a 23 = −1 then U(ϕ 1 , ϕ 2 ) = 3/2 + cos ϕ 1 + cos ϕ 2 + cos(ϕ 1 + ϕ 2 ), where from (26) c = 3 2 , which has minima at the points (24) at which cos ϕ 1 = − 1 2 = cos ϕ 2 . This potential is the negative of U as plotted in figure 1, up to an additive constant, and so the points in red in figure 1 mark the zeroes of U near the origin. In this case the minima form a hexagonal lattice array in the (ϕ 1 , ϕ 2 ) plane, which is also evident in the contour plot of U in figure 3 (right). We consider solitons for the fully symmetric case in section 5.
The exact solutions obtained by this reduction are not necessarily stable against perturbations, however, for example in (1) below the solution satisfies ϕ 1 + ϕ 2 = 0, which is a constraint that is not necessarily maintained with respect to perturbations. We return to the question of stability in section 6.
Exact solutions are significant because they furnish explicit details of solution properties with respect to underlying parameters. For the double sine-Gordon equation there are exact static solutions which describe double solitons located at arbitrary points on the line, a property which extends to larger values of N, and can therefore be analysed precisely. Exact solutions also indicate the precise conditions for which static configurations are possible. For models with a single scalar field, such as the sine-Gordon equation, static solitons interpolate only between adjacent minima of the potential, but for two or more fields it is not immediately clear whether all adjacent minima are associated with static solitons. Taking the example a 12 = a 13 = a 23 = −1, we find that there exist exact soliton configurations which join the maxima of −U as shown in figure 1, together with their symmetric counterparts. Even if these solutions are unstable, there nevertheless exist stable static configurations which interpolate between these maxima. On the other hand, there is no exact solution interpolating between (− 2π 3 , − 2π 3 ) and (− 2π 3 , 4π 3 ), suggesting that a static soliton does not exist for this case. We have plotted −U for the symmetric case in figure 3 (right), showing the interpolating trajectories.
Hence there exist six trajectories which interpolate between the origin in the (ϕ 1 , ϕ 2 ) plane and a nearby minimum of U, for restricted coefficients a ij > 0 as stated, which follow from (a) and (b) by cyclic symmetry. For each such solution there also exists an anti-soliton which interpolates between (0, 0) and the negative of the soliton endpoint. Altogether then, there exist 12 distinct trajectories, six solitons and six antisolitons, from (0, 0) to a neighbouring minimum of U. Numerical investigations indicate that these trajectories exist also for general coefficients a ij > 0.
The other cases of the double sine-Gordon equation, for a 13 = a 12 and a 23 = a 13 as in (d) and (f) in section 4.3, are related to (28) by cyclic symmetry. We can find exact solutions of (28) by means of the first integral ϕ 2 = U(ϕ) which follows from (20), where from (19) we have for some constant c . Of the known exact solutions, whether in static form or as travelling waves, as derived for example in [32,33], we are interested only in those with finite energy on the real line, of which there are three types depending on the relative values of a 12 , a 13 . We wish to determine the dependence of properties on a 12 , a 13 , and we show in particular that static double solitons exist for certain values of a 12 , a 13 . This accords with the picture that the underlying constituent fields θ 1 , θ 2 , θ 3 can give rise to distinct multi-soliton configurations, a property which has received little attention except in [45].
Corresponding to the symmetries discussed in section 4.1, (28) is invariant under ϕ → −ϕ, hence to each soliton ϕ there is an associated antisoliton −ϕ, and (28) is also invariant under which relates solitons with opposite signs of a 12 . There are three cases to be considered: which has zeroes at cos ϕ = −a 12 /(2a 13 ), corresponding to the stationary point (24), which exists because |a 12 | 2|a 13 |.
These three cases account for all generic values of the coefficients in the a 12 , a 13 plane.
For the first case, in which a 12 > 0, a 12 + 2a 13 > 0, the explicit solution can be expressed in the form where α = √ a 12 + 2a 13 , and where the arbitrary point x 0 can be regarded as the centre of the soliton. Evidently lim |x|→∞ cos ϕ(x) = 1 and cos ϕ(x 0 ) = −1, hence ϕ is a kink or anti-kink soliton which interpolates smoothly between adjacent maxima of the cosine function, for example from 0 to 2π passing through ϕ(x 0 ) = π. The energy density is given by E = (ϕ ) 2 + U(ϕ) = 2(ϕ ) 2 , explicitly: and the total energy for a 13 > 0 is and for − a 12 2 < a 13 < 0 we obtain These expressions are useful in order to determine whether a numerical solution with the same endpoints is of higher energy and therefore unstable. For a 12 > 2a 13 , E(x) has a turning point with a corresponding maximum at x = x 0 . This is verified by evaluating E (x 0 ) = −16a 12 (a 12 − 2a 13 ), hence E (x 0 ) < 0 for a 12 > 2a 13 , recalling that a 12 > 0. In this case the soliton is located at x = x 0 , since this is the point at which the energy density is peaked. For a 12 < 2a 13 , (32) represents a double soliton since E(x) has turning points at x = x ± , at which cos ϕ = −a 12 /(2a 13 ), where At each of these points E(x) attains a local maximum, and E(x 0 ) is a local minimum.
As an example of such a double soliton we plot ϕ(x) in figure 2 (left) as a function of x (in blue) for |x| < 8, for the parameters a 23 = a 12 = 10 −4 , a 13 = 1 with x 0 = 0, together with E(x) (in green). Evidently E(x) attains its maximum at the distinct points x = x ± . These locations can be spaced arbitrarily far apart by suitable choice of a 12 , a 13 , as follows from (34), since the difference x + − x − can be made arbitrarily large by choosing a 12 to be sufficiently small, for fixed a 13 . We also plot E(x) in figure 2 (right) as a function of x for 0.01 < a 12 < 1 with a 13 = 1, showing the distinct double peaks of E moving together as a 12 increases; eventually these peaks coincide as a 12 increases to the value a 12 = 2a 13 (not shown), and then for a 12 > 2a 13 the two solitons remain co-located at x = x 0 .
The second case of interest, item (b) for a 12 < 0, a 12 < 2a 13 , is obtained from the first by application of the transformation (30). This leads to an explicit solution similar to (32), and includes double solitons located at distinct positions.
The third and final case in (c) occurs for |a 12 | 2|a 13 |, a 13 < 0, which corresponds to the fixed point (24), where U attains a minimum at cos ϕ = −a 12 /(2a 13 ). The condition a 13 < 0 follows from the requirement that s 3 < 0, with U as in (31). Two distinct solitons exist in this case, with the first given explicitly by where β 2 = (4a 2 13 − a 2 12 )/(−2a 13 ). This solution is antisymmetric about x = x 0 , with ϕ(x 0 ) = 0, and interpolates between adjacent minima of U located at ϕ = ±2 arctan β √ 2|a 13 | a 12 −2a 13 , at each of which cos ϕ = −a 12 /(2a 13 ), as is evident from the alternative expression The energy density is peaked at x = x 0 , and so double solitons do not exist for this case. The second soliton is given by which differs slightly from (35). We can also write from which we see that this soliton interpolates between adjacent minima of U at ϕ = π ± 2 arctan β √ 2|a 13 | a 12 +2a 13 , at each of which again cos ϕ = −a 12 /(2a 13 ). The energy density for the solitons (35) and (37) is given by where the minus sign is for (35) and the plus sign for (37). The total energy of (35) is E = 4β − 4a 12 τ/ √ −2a 13 , where τ = arccos − a 12 2a 13 , and for (37) we have E = 4β + 4β(π − τ ) cot τ . For a 23 = a 12 we have therefore obtained explicit solutions of (18) as single and double sine-Gordon solitons, which by cyclic symmetry extend to the cases a 23 = a 13 and a 13 = a 12 . The stability of these solutions is investigated in section 6, but in any case the existence of stable solutions which interpolate between the corresponding endpoints is assured because of the existence of these exact solutions.

N = 3 solitons for the fully symmetric case
We have discussed properties of the Kuramoto-sine-Gordon system (3), such as the elementary excitations in section 3, for arbitrary a ij > 0 for general N. Since the system depends on N(N − 1)/2 independent coefficients a ij , which can be of any sign, there are many possible soliton configurations which generally require caseby-case consideration. Of particular interest, however, is the fully symmetric case in which all coefficients a ij are all equal, where by rescaling we can set a ij = ±1. This describes a complex system of identical nodes with global couplings of equal strength which is invariant under the transposition of any two nodes. For this case the elementary excitations have equal mass, since for a ij = 1 the matrix M defined by (13) has N − 1 eigenvalues μ 2 i = N, plus a single zero eigenvalue. For N = 3, for example, we have from (16) s 1 = 3 = s 2 and so μ 2 1 = μ 2 2 = 3. We now consider solitons for the fully symmetric case with N = 3 as special cases of the exact solitons already discussed.

N = 3 symmetric solutions for a ij = 1
For N = 3 with a ij = 1, the potential is given by U(ϕ 1 , ϕ 2 ) = 3 − cos ϕ 1 − cos ϕ 2 − cos(ϕ 1 + ϕ 2 ) and has minima at the points 2π(m, n) on the square lattice in the (ϕ 1 , ϕ 2 ) plane, for integers m, n, at which cos ϕ 1 = 1 = cos ϕ 2 . We have plotted U in figure 1, and a contour plot is shown in figure 3 (left) showing the square lattice array of zeroes of U. There exist a 12 static solitons emanating from the origin (0, 0) to adjacent zeroes marked in in red in figure 3 (left), consisting of six static solitons and six antisolitons, as follows from the solutions discussed for the six cases (a)-(f) in section 4.3. Of these, (a), (c) and (e) are sine-Gordon solitons satisfying ϕ 1 + ϕ 2 = 0, ϕ 2 = 0, ϕ 1 = 0, respectively, with the linear trajectories marked in cyan in figure 3 (left). The exact solutions are given by (27) with α = √ 3, with a total energy of 16/ √ 3 units. The remaining three trajectories correspond to (b), (d) and (f) which are double sine-Gordon solitons satisfying ϕ 1 = ϕ 2 , 2ϕ 1 + ϕ 2 = 0, ϕ 1 + 2ϕ 2 = 0, respectively, marked in blue in figure 3 (left). The exact solutions are given by (32) with α = √ 3, and have a total energy E = 8 √ 3 + 4 √ 2 tanh −1 2/3 . We show in section 6.3 that the exact sine-Gordon solitons are stable for a ij = 1, but that the double sine-Gordon solitons are unstable. There nevertheless exist stable trajectories which interpolate between the zeroes of U which are not linear in ϕ 1 , ϕ 2 and remain to be found. By periodicity, similar trajectories exist originating from all other zeroes of U.
There exists a second double sine-Gordon soliton given by (37), equivalently (38), which also solves ϕ = −sin ϕ − sin 2ϕ and interpolates between 2π 3 and 4π 3 , together with the 2π-periodic translations of these points. There are three possibilities, corresponding to ϕ 1 = ϕ 2 = ϕ, 2ϕ 1 = −ϕ 2 = 2ϕ, ϕ 1 = 2ϕ 2 = −2ϕ which are shown in cyan in figure 3 (right), where we have included 2π-periodic copies in order to form a ring of six soliton trajectories encompassing the origin. Figure 1 also shows the trajectory (in green) interpolating between (− 2π 3 , 4π 3 ) and ( 2π 3 , 2π 3 ) along the ridge between these maxima, and is stable. The energy of this soliton is 2 √ 2 3 √ 3 − π /3 units. Note that figure 3 does not show all possible trajectories, only those passing through the origin in figure 3 (left), and those immediately surrounding the origin in figure 3 (right); to these must be added the 2π-periodic copies. Figure 3 also shows that static solitons do not exist for all adjacent zeroes of U, for example there is no soliton connecting the zeroes (0, 0) and (2π, 4π) (left) or the zeroes (− 2π 3 , − 2π 3 ) and (− 2π 3 , 4π 3 ) (right) of U.

Stability of static solutions
In this section we investigate linear stability of static solutions of (7), firstly for general N, then for N = 3 for the specific equation (18), in order to determine if a given solution, whether in numerical or exact form, is stable. As discussed in the introduction, since a superpotential does not exist in closed form we cannot determine by direct integration an explicit configuration which minimizes the energy in any sector. Exact or numerical static solutions, although they are stationary points of the energy functional, are not necessarily the solutions of lowest energy, and so we need to test such solutions for stability.

The perturbed equations for general N
and so, writing the perturbed fields as θ i (x, t) = θ s i (x) + ε i (x) cos ωt, we obtain the N linearized equations Denote ψ = (ε 1 , ε 2 , . . . , ε N ) T then we have the matrix equation −ψ + M(x)ψ = ω 2 ψ, where M = (m ij ) is the symmetric N × N matrix of second derivatives of V, as defined in (5), evaluated at the static solution, i.e.
This generalizes the expression (13), which is the special case of M for vacuum solutions satisfying cos(θ s j − θ s i ) = 1. We write (40) as the matrix Schrödinger equation Hψ = ω 2 ψ, where where I N denotes the N × N identity matrix. For the solution θ s i (x) to be stable, all eigenvalues ω 2 of H must be non-negative. The operator H is Hermitian in the space of square-integrable functions ψ, for which meaning that H is a symmetric operator which we assume to be self-adjoint. H has two zero eigenvalues, the first corresponding to the zero mode of oscillation obtained by differentiating (39), with the eigenfunction ψ 0 = (θ s 1 , θ s 2 , . . . , θ s N ) T , as follows from the translational invariance of the system (3). The second zero eigenvalue corresponds to the constant eigenfunction (1, 1, . . . , 1) T which lies in the null space of M, and arises from the massless excitations of the fields θ i .
It follows from (14) that for any N-vector u = (u 1 , u 2 , . . . , u N ) T we have and so M(x) has non-negative eigenvalues for sufficiently large |x| that cos(θ s j − θ s i ) > 0 for all i, j, assuming that a ij > 0. For static solitons, however, M(x) always has negative eigenvalues for small x, since cos(θ s j − θ s i ) takes negative values. It follows from Hψ = ω 2 ψ that ψ , ψ + ψ, Mψ = ω 2 ψ, ψ but due to the negative eigenvalues of M we cannot conclude that ω 2 > 0.
For static solutions determined directly from the superpotential W, however, the eigenvalues of H are always non-negative. Such solutions satisfy θ i = W θ i where W is defined by (10), as explained in section 2.2, and in this case we can factorize H according to  (44) is indeed satisfied. It follows in the usual way that H has nonnegative eigenvalues, since from (43) we can write H = S † S, where S † is the adjoint of S, then from Hψ = ω 2 ψ we obtain ψ, Hψ = Sψ, Sψ = ω 2 ψ, ψ , and hence ω 2 0. The factorization of H provides a method to demonstrate stability of a numerical solution, even though an explicit superpotential is not known. For a given static solution and hence for a given matrix M(x), we solve the matrix Riccati equation M = A 2 + A directly for A, assuming that a regular solution exists, with boundary conditions to be determined. If there exists an acceptable solution, then we have in effect found the matrix (W θ i ,θ j ) of second derivatives of W evaluated at that particular solution, and stability of this solution then follows from the factorization property of H. The matrix Riccati equation M = A 2 + A , for given M, can be converted to the second-order linear equation R = MR in the usual way by means of the substitution A = R R −1 , where we require the N × N matrix R(x) to be nonsingular at all points x.
We proceed more directly in section 6.4, however, by estimating the lowest eigenvalue of H using the Rayleigh-Ritz variational method. With suitable trial states this method is able to detect negative eigenvalues of H, and hence demonstrate instability of the soliton configuration. If no such negative eigenvalues are found, and if the trial function successfully detects the zero eigenvalue corresponding to zero mode oscillations then, although strictly we are unable to draw any conclusion, we state that the soliton 'appears to be stable'.
We first consider sine-Gordon solitons for N = 3 for which we can find exact stability results. Generally, however, −ψ + M(x)ψ = ω 2 ψ can be solved exactly only in special cases, such as soliton models with a single scalar field with a polynomial potential [46], or some two-component scalar field theories [16].
is Hermitian with respect to the inner product ψ, φ A , for twocomponent functions ψ, φ, defined by where we note that A is positive-definite. H has a zero eigenvalue with the zero mode eigenfunction ψ 0 = (ϕ s 1 , ϕ s 2 ) and we wish to determine whether H has any negative eigenvalues which render the soliton unstable.

Stability properties of sine-Gordon solitons for N = 3
In section 4.3 we listed six cases for which the N = 3 static equations can be solved exactly by either the sine-Gordon or double sine-Gordon equations. Since these follow from the first two cases (a) and (b) by symmetry, we consider stability for (a) and (b) only, for which a 23 = a 12 . We suppose that U as defined in (19) has minima at cos ϕ 1 = 1 = cos ϕ 2 which is the case provided that s 2 > 0, a 12 + a 13 > 0, which for a 23 = a 12 reduces to a 12 > 0, 2a 13 + a 12 > 0. For the sine-Gordon soliton (27) we can solve the corresponding perturbation equations exactly, and find that the soliton is unstable if the inequality (55) holds, but is otherwise stable. For the double sine-Gordon system we resort to numerical observations and observe that the soliton is unstable for a 13 a 12 9. Consider firstly case (a) for the sine-Gordon soliton (27), then we wish to solve (45) for ψ = (ψ 1 , ψ 2 ), where M r (x) is given by (46) in which ϕ s 1 = −ϕ s 2 = ϕ, where cos ϕ(x) is given by (27). We have M r = 2a 12 cos ϕ + a 13 −a 12 cos ϕ + a 13 −a 12 cos ϕ + a 13 2a 12 cos ϕ + a 13 .
Define ψ ± = ψ 1 ± ψ 2 , then (45) decouples into two independent equations: It follows from ψ, ψ A < ∞ that ψ ± are each square-integrable. The Schrödinger equation (51), with a sech 2 potential, is well-known as that which determines linear perturbations about static sine-Gordon solitons, and can be solved exactly in terms of hypergeometric functions, see Morse and Feshbach [47], equation (12.3.22), page 1651, also [48]. The parameters in [47] correspond to ours according to Bound states occur only if ε < ve −2μ , i.e. if ω 2 < 3a 12 , and for our case there is only one such bound state at b = 1 which is the zero mode state with ω 2 = 0. The corresponding wavefunction ψ is given by ψ − (x) = sech √ 3a 12 (x − x 0 ) together with ψ + = 0, which solves (51) and (52) with ω 2 = 0. There is a second bound state with ψ − = 0 which we find by solving (52) for ψ + in the same way. In this case we identify the parameters in Morse and Feshbach [47] as in (53), except for Bound states occur for non-negative integers n such that b n = v + 1 4 − n − 1 2 is positive and, again, only the value n = 0 is allowed, hence b = b 0 = 11 12 − 1 2 , for which with a corresponding wavefunction ψ + (x) = sech b √ 3a 12 (x − x 0 ). The sine-Gordon soliton (27), which exists for − 1 2 < a 13 a 12 , is therefore unstable for but is otherwise linearly stable, since (51) and (52) have bound state solutions with either ψ − = 0 or ψ + = 0, and these are the only such solutions. The case ψ + = 0 corresponds to the zero mode solution, but ψ − = 0 leads to a nontrivial solution for ψ + representing a meson-soliton bound state which does not occur in the simple sine-Gordon model for N = 2. The value of the ratio a 13 a 12 given by the right-hand side of (55) can be regarded as a critical coupling, at which point the stable soliton changes its form. We investigate the numerical form of the solution when (55) holds in section 7.1, where we find a stable solution for a 12 = 1 = a 23 and a 13 = − 2 5 . For the double sine-Gordon equation, which appears in (2) by setting a 23 = a 12 together with ϕ 1 = ϕ 2 = ϕ, the exact soliton is given explicitly by (32), and the reduced matrix M r is given by M r = 2a 12 cos ϕ + a 13 cos 2ϕ −a 12 cos ϕ + a 13 cos 2ϕ −a 12 cos ϕ + a 13 cos 2ϕ 2a 12 cos ϕ + a 13 cos 2ϕ .
Define, as before, ψ ± = ψ 1 ± ψ 2 then (45) decouples into the independent equations In this case an exact solution for general a 12 , a 13 is not available, however the zero mode eigenfunction with An exact solution is known for a 13 = 0 in which case we find that the soliton is unstable. Numerically, we can use trial functions as described in the following section 6.4 to probe the stability of the soliton (32), which we find is unstable for a 13 a 12 9, but otherwise appears to be stable.

Eigenvalue upper bounds for N = 3 solitons
We can obtain an upper bound λ on the lowest eigenvalue of H by the Rayleigh-Ritz variational method, in which we choose a normalizable trial function ψ tr (x), depending on several parameters, and evaluate λ = ψ tr , ψ tr A + ψ tr , M r ψ tr A ψ tr , ψ tr A .
We minimize this expression with respect to the parameters in ψ tr (x) to find the smallest value of λ, then if λ < 0 the static solution ϕ s 1 , ϕ s 2 is unstable, since ω 2 λ. In order to demonstrate that a given solution is unstable, therefore, it is sufficient to show that the numerator of (59) is negative for some choice of parameters for ψ tr (x), i.e. we show that where A, U (2) are defined in (47) and (48). In our calculations we chose two-component Gaussian trial functions of the form which depend on parameters a 1 , a 2 , b 1 , b 2 , x 1 , x 2 , although one can of course refine this by choosing more general trial functions which are orthogonal to the zero mode eigenfunction ψ 0 = (ϕ s 1 , ϕ s 2 ). In some cases it is sufficient to choose with only two parameters b, x 0 in order to prove instability. In the expression (60) the first term can be evaluated exactly for these trial functions by means of a computer algebra system. The second term can be evaluated numerically, for any exact or numerical solution, by numerical integration over an interval [−L, L] for sufficiently large L < ∞, due to the exponential fall-off of the trial functions. By means of this variational method we determine, for example, that the soliton (32) is unstable for a 13 a 12 9.  the energy of the sine-Gordon soliton of 16 a 12 /3 ≈ 9.24 units. In figure 4 (right) we plot the computed solution (green) as a trajectory on the surface of the potential −U(ϕ 1 , ϕ 2 ), showing the soliton interpolating between the endpoints (0, 0) and (2π, −2π). By comparison, the exact sine-Gordon soliton (red) descends into the valley along the line ϕ 1 + ϕ 2 = 0, leading to a higher energy.

A general numerical example
As a second, more general, example we choose a 12 = 1/10, a 13 = 1, a 23 = 3/10 with a trajectory interpolating between (0, 2π) to (2π, 0) in the (ϕ 1 , ϕ 2 ) plane. This is related to the previous example, upon replacing ϕ 2 → ϕ 2 − 2π, the main difference being that the equality a 23 = a 12 is now broken, and so the symmetries previously used are no longer valid. We can assume, however, that both ϕ 1 (−x) = 2π − ϕ 1 (x) and ϕ 2 (−x) = 2π − ϕ 2 (x) hold, i.e. that π − ϕ(x) is odd in x for ϕ = ϕ 1 , ϕ 2 , consistent with (18). Hence ϕ 1 (0) = ϕ 2 (0) = π and we can choose ϕ 2 (0) as a starting value for the shooting method. We have found two distinct solutions for ϕ 1 , ϕ 2 as plotted in figure 5. The left trajectory satisfies the second order equation (18) and the consequent integrated equation (20) to high accuracy but is unstable. This trajectory is also shown in figure 6 (in red) and evidently follows a highly convoluted path between the specified endpoints. The other much simpler trajectory in figure 5 (right), for which ϕ 1 + ϕ 2 is approximately equal to 2π, appears to be stable according to the tests described in section 6.4, and joins the same endpoints as shown in figure 6 (green). The energies of these configurations are 28.56 and 4.12 units for the unstable and stable trajectories, respectively.
It remains to determine stable trajectories interpolating between other zeroes of U for these same values of a 12 , a 13 , a 23 , for example between (0, 0) and (0, 2π), as well as for those cases in which the exact double sine-Gordon solitons are unstable.

The Kuramoto-sine-Gordon model for N 4
We have analysed the N = 3 case of the Kuramoto-sine-Gordon system (3) in detail in order to derive underlying properties of the static equations, in particular we have systematically discussed the symmetries of these equations in section 4.1, and how they can be used to reduce the static equations to lower order systems with exact soliton solutions. These properties extend to any N, but there are now many more possible reductions.
The system (63) has symmetries which follow from those of (3), such as permutational symmetry obtained by transposing any two nodes. For example, (3) is invariant under the (12) transposition θ 1 ↔ θ 2 provided that a 1i ↔ a 2i for all i = 3, . . . , N, which for (63) is equivalent to ϕ 1 → −ϕ 1 , ϕ 2 → ϕ 1 + ϕ 2 . Such symmetries can also be combined with invariances such as sign reversal, ϕ i → −ϕ i for all i. We can use these to reduce the system to one of lower dimension, as for N = 3, for example we can set θ 1 = θ 2 (i.e. ϕ 1 = 0) provided that a 1i = a 2i for all i = 3, . . . , N, to obtain a system of N − 2 equations, which can then be further reduced. In general, we can reduce the system to the sine-Gordon equation by setting θ i = θ j for all pairs i, j but one.
Then, by setting ϕ 2 = −ϕ together with a 14 = a 23 we obtain the sine-Gordon equation, or by setting a 13 = a 12 with either ϕ 2 = 0 or ϕ 2 = −2ϕ, we obtain the double sine-Gordon equation. A different reduction, given in [8], is obtained by setting ϕ 1 = −ϕ 2 = −ϕ 3 = ϕ in (64) together with a 14 = a 12 , a 34 = a 23 , then again we obtain the double sine-Gordon equation, this time in the form ϕ = (a 12 + a 23 ) sin ϕ + a 24 sin 2ϕ. Equivalent reductions follow by cyclic symmetry of the nodes. We then use known exact solutions to construct solitons with boundary values that are appropriate for the chosen values of a ij . In particular, these reductions apply to the fully symmetric case for which the coefficients a ij are all equal.
For general N, reductions to the double sine-Gordon equation are less evident, although we have verified explicitly that this is possible for N 8. These reductions apply also to the time-dependent case, hence solutions of the sine-Gordon and double sine-Gordon equations can be used to construct time-dependent solutions for any N, subject to stability considerations.

Conclusion
The Kuramoto-sine-Gordon system (3) can be regarded as a 1 + 1-dimensional version of 3 + 1-dimensional skyrmion systems which describe pion interactions and nucleons by means of chirally invariant Lagrangians [8]. Even for a single pion system (the case N = 1) it is known that there are many possible skyrmion configurations [9,10], and so there is merit in studying the simpler 1 + 1-dimensional Kuramoto-sine-Gordon system and its possible soliton configurations, as a means of determining the effect of different network connections. A similar remark holds for 2 + 1-dimensional skyrmions in systems of relativistic chirally invariant field theories. The wide variety of skyrmion configurations which arise in 'magnetic skyrmionics' [6] suggests that similar configurations can occur in relativistic planar skyrmion systems, which again motivates the study of the simpler Kuramoto-sine-Gordon equations.
We have developed properties of the Kuramoto-sine-Gordon system for general N and general coupling coefficients a ij , including the particle content expressed as elementary excitations of mass μ i , where the massless excitations conveniently decouple from the system and can therefore be ignored. Since there are many possible soliton configurations, depending on the N(N − 1)/2 coefficients a ij , we have focussed in particular on N = 3, but even then the classification of stable static solitons for any a ij is not complete. We have shown by means of exact solutions that a wide variety of static solitons exist, some of which are double solitons located at arbitrary points on the line. Whether this property generalizes for larger N, in which any number of static solitons exist at arbitrary locations on the line, remains to be determined. Properties of time-dependent solitons and how they imbed into the system for large N also remain to be investigated, for example whether there exist stable periodic solutions which generalize the exact sine-Gordon breather solution. Further generalizations of the system (3) are also possible since one could include, for example, arbitrary coefficients for the kinetic terms in (1), as well as additional antisymmetric couplings for the potential terms [8].

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.