Multisoliton complex systems with explicit superpotential interactions

We consider scalar field theories in 1+1 dimensions with N fields φ1,…φN which interact through a potential V=V(φ1,…φN) , which is defined in terms of an explicit superpotential W. We construct W for any N in terms of a known superpotential w for a single-scalar model, such as that for the sine-Gordon equation or the ϕ 4 model, leading to an expression for V which has multiple minima that supports solitons. Static solitons which minimize the total energy in each soliton sector appear as solutions of first-order Bogomolny equations, which have a gradient structure. These are identical in form to equations which arise in the context of synchronization phenomena in complex systems, with the space and time variables interchanged. The sine-Gordon superpotential, for example, leads to an explicit periodic superpotential W for N scalar fields, with associated Bogomolny equations that are equivalent to the well-known Kuramoto equations which describe the synchronization of identical phase oscillators on the unit circle. The known asymptotic properties of the Kuramoto system, for both positive and negative coupling constants, ensure that finite-energy solitons exist for any given set of intermediate values imposed at the origin. Besides the models derived from the sine-Gordon equation, we investigate ϕ 4 and ϕ 6 models with N scalar fields and show numerically that solitons again exist over a wide range of parameters. We also derive general properties of the elementary meson excitations of the system, in particular we show that meson-soliton bound states exist over a restricted range of mass parameters with respect to an exact solution of the ϕ 6 system for N = 3.


Introduction
Solitons have been widely investigated as classical configurations which also exist as particles in quantum field theory, due to the existence of topological charges which ensure that the configurations are non-dissipative with respect to both classical and quantum fluctuations.For reviews and summaries we refer to [1][2][3][4], particularly the description of single solitons, such as kinks in scalar field theory.
Here we investigate systems of solitons in scalar field theory in 1 + 1 dimensions in order to describe multiparticle states.Whereas the usual focus is on single solitons as isolated particles, multiparticle configurations can appear in a variety of ways, for example as multisolitons within a single system such as the sine-Gordon model [1] (section 5.3), with multisoliton collisions [5] and in multisoliton chains [6].In higher dimensions there exist low-charge skyrmions [1] (section 9.4) as well as multi-skyrmion configurations [2] (section 7.2.3), and multiparticle systems can also appear in the context of lattice structures, such as Skyrme crystals [2] (section 7.2.5).
A general method of constructing N-particle states is to introduce N scalar fields φ i , (i = 1, . . .N), which interact by means of a potential V = V(φ 1 , . . .φ N ).We can regard these fields as comprising a complex system consisting of N nodes at each of which there is a scalar field φ i which interacts with the other fields φ j through V, with a strength that depends on a coupling coefficient a ij .Solitons in such coupled scalar fields have long been studied, particularly for two fields, see for example [7][8][9][10][11][12] and references therein.Many two-field models have been analyzed by means of a superpotential [13][14][15][16][17][18].Models with three or more fields, including the case of N fields, have also been investigated [19,20].
A significant difficulty in the development of multiscalar models is that for a given V the superpotential W is generally not known explicitly, and so static solutions must be determined numerically by solving a second-order system.Such solutions are stationary points of the energy density, but need not correspond to a local minimum and could therefore be unstable.Furthermore, static solutions interpolate between adjacent minima of V, and so for large N there could be many such trajectories which, without the benefit of an explicit superpotential, are difficult if not impossible to classify.By contrast, if W is known explicitly, it is sufficient to solve a system of first-order equations (the Bogomolny equations, see [1] section 1.3), and the solutions are guaranteed to be stable.Explicit superpotentials are known for several multiscalar models, see for example [12][13][14][15][16][18][19][20], but not for others such as the Kuramoto-sine-Gordon system, which is the simplest case of a chirally invariant system of N interacting fields [21].Some exact solutions and other properties can be derived for N = 3, 4 in this model, but little is known for general N.
We proceed here by choosing as a starting point an explicit expression for W as a function of the scalar fields, from which we define V, and then derive properties of the first-order system.The choice of W is not arbitrary, because the derived potential V must have the appropriate vacuum structure in order to support solitons.We begin with known superpotentials for a single-scalar field theory, such as the sine-Gordon model, and from these construct W for a system of N scalar fields.The subsequent first-order equations have a gradient structure with respect to W and in some cases are well-known in a completely different context.The sine-Gordon superpotential, for example, leads to a system of N equations which are a special case of the widely-studied Kuramoto model, for which many rigorous analytic as well as numerical results are known in the context of synchronization phenomena in complex systems.Other first-order systems, such as those which follow from the φ 4 and φ 6 superpotentials, have not been studied previously yet describe interacting solitons and are also relevant to synchronization phenomena.

Outline and summary
We begin in section 2 by constructing an explicit superpotential W for N scalar fields from the known superpotential w for a single-scalar field.Because w leads to a potential V which allows solitons, the same is generally true for the model with N scalar fields determined by W. In section 3 we derive properties of the N-field system, such as the definition of V from W, the Bogomolny equations, properties special to N = 3, the linearized equations and stability, factorization of the Hamiltonian, and the elementary excitations of the system.Although some of these properties are well-known for N = 2 [14], for completeness we outline a full development for any N.
In section 4 we briefly describe what at first appears to be an unrelated topic, synchronization phenomena in complex systems with N nodes.Such systems are usually modeled by first-order equations which evolve from generic initial values such that asymptotically the variables of the system are highly correlated.Synchronization has been intensively studied for the Kuramoto model, which dates from 1975 [22], but appears here in a soliton context.There is an interplay between the time t and space variables x of the synchronization and soliton models respectively, where x is allowed to increase from an initial location at x = x 0 in the positive direction, corresponding to forward time evolution, as well as decrease in the negative direction from x = x 0 as if evolving backwards in time.By this means we can calculate not only soliton configurations numerically from random initial values at x = x 0 , but can also apply established analytic properties of synchronization models to deduce that static finite-energy solitons exist for such models, for any N.
In section 5 and following we investigate specific examples, starting with the periodic sine-Gordon superpotential and the subsequent Bogomolny equations which can be identified as a Kuramoto system.For N = 3 (sections 5.2 and 5.3) we find exact solutions which demonstrate explicit soliton properties.A remarkable property of the Kuramoto system with identical oscillators is that of partial integrability, which we discuss in section 6, where the solitons for any N have an explicit expression in terms of a single complex variable ζ.
In sections 7-9 we consider in turn the antisymmetric sine-Gordon superpotential and the φ 4 , φ 6 models with N scalar fields, each with different characteristics.In all cases solitons exist for any N and although these are mainly numerical observations, asymptotic properties of the solitons can in some cases be described exactly.The φ 6 models are easier to analyze than the φ 4 systems for several reasons, in particular we can find closed form solutions for the N = 3 case of the φ 6 model, with which we can determine explicitly the interactions between the elementary meson excitations of the system with the soliton (section 9.3).Unlike the singlescalar φ 6 model, which does not allow any bound meson-soliton states [23], the N = 3 model allows several bound states for specified meson masses.We finish with concluding remarks in section 10.

Explicit superpotentials
We construct multisoliton systems with N scalar fields φ 1 , . . .φ N as defined by the Lagrangian density: in standard relativistic notation.We define the potential V by first choosing an explicit superpotential W = W(φ 1 , . . .φ N ), then V is given by where W φi ≡ ∂W/∂φ i , and so V is non-negative with all zeroes (i.e.vacuum states) determined by the N equations W φi = 0. Stable, static kink solitons are solutions of a system of N firstorder Bogomolny equations in gradient form, namely φ ′ i = W φi , where φ ′ i ≡ dφ i /dx.In order to ensure that V has a vacuum structure that supports solitons, we construct W directly from superpotentials w = w(φ) that are well-known for single-scalar models, namely the sine-Gordon, φ 4 and φ 6 models.For any suitable potential V, however, such as a polynomial of any degree, we can in principle find w by integrating V(φ) = 1 2 w ′ (φ) 2 directly, and then follow the procedure outlined here to construct soliton models with N scalar fields.Some examples of explicit superpotentials w are given in [24].
For a given w = w(φ) we define W = W(φ 1 , . . .φ N ) according to where the coupling coefficient a ij determines the strength of interaction between distinct nodes i, j.The terms with i = j in this summation merely contribute additive constants to W, and so do not affect the dynamics, and it is convenient in some cases to specify a ii = 0.If w is an even (respectively, odd) function of φ, then a ij must be symmetric (respectively, antisymmetric) in i, j, as follows from (3) by exchanging the indices i ↔ j under the sum.Evidently, W is a function of the differences φ j − φ i only, hence the subsequent equations are invariant under the translation φ i → φ i + φ 0 for any constant φ 0 .Consequently V is a function of the differences φ j − φ i only.There are two reasons why we write W in this way: (i) For even superpotentials w the configurations with φ j = φ i for all i, j are zeroes of V and are therefore fixed points of the Bogomolny equations and so, subject to stability requirements, lead to synchronized asymptotic configurations with finite soliton energy.Such constructions are well-known for synchronization models, see for example the discussion in section V in [25].Similarly, in applications to social dynamics, such as the study of global behaviors in large groups of autonomous agents, organized global patterns can form even when individual agents interact only at a local scale, and consensus is achieved when the interactions depend only on the differences of the positions of the agents [26].Emergent phenomena have been widely studied in first-order gradient models where the potential again depends only on the differences of the agent positions [27].The properties of such models are relevant to solitons even when consensus and/or synchronization does not occur; (ii) For models defined on the unit circle, such as the Kuramoto model, the translational invariance φ i → φ i + φ 0 is equivalent to SO(2) rotational invariance which can be generalized to higher dimensional models defined on the unit sphere, as occurs for the Kuramoto model [28], although we do not develop this property here for solitons.
Having described the general procedure by which soliton systems are constructed, we now turn to four specific models.

Single-scalar fields and the superpotential
Single-scalar superpotentials w and their static solitons are well-known for several models which we extend to multiscalar models, where V is related to w by V(φ) = 1 2 w ′ (φ) 2 : (i) The sine-Gordon model (even superpotential): leads to the sine-Gordon equation 2∂ µ ∂ µ θ = −K 2 sin 2θ.The static soliton is given by where x 0 is an arbitrary constant which reflects the translational invariance of the system.The solution interpolates between 0 and π, and the antisoliton is given by θ(−x).(ii) The sine-Gordon model (odd superpotential): results in the sine-Gordon equation in the form 2∂ µ ∂ µ θ = K 2 sin 2θ, which is related to (4) simply by replacing θ → θ + π 2 , however the extension to multiscalar models differs according to whether w is even or odd, and so we retain this form for w. (iii) The φ 4 model: where λ = 9K 2 /2 and w(φ) is an odd function of φ.The static soliton has the explicit expression which interpolates between the minima of V(φ) at φ = ±1.(iv) The φ 6 model: where w is even, with the static soliton given by which interpolates between the minima of V at φ = 0, 1, while −φ(x) interpolates between the minima at φ = −1, 0.
While we have chosen these examples for further development, there are many others of interest, such as the polynomial superpotential corresponding to where m, n are integers [24], for which w is generally neither even nor odd.

Soliton properties and the superpotential
We wish to solve the static second-order Euler-Lagrange equations φ ′ ′ i = V φi which follow from the Lagrangian density (1) for the N fields φ i , by means of first-order equations which generalize the well-known Bogomolny equations for single-and double-scalar systems [1,14].
Given W, we define the potential V according to (2).If V is a given function of the scalar fields, then (2) is a first-order partial differential equation for W, known as the eikonal equation, which is a special case of the Hamilton-Jacobi equation.Although solutions are known to exist for well-behaved functions V, explicit solutions are generally not available except in rare cases [29].In our approach we first define W and then calculate V, which enables us to write out explicit Bogomolny equations, an approach which has previously been followed for multiple scalar fields but in different contexts [14,30,31].

The Bogomolny equations
The Lagrangian (1) can be written where φi ≡ ∂φ i /∂t, φ ′ i ≡ ∂φ i /∂x, and the total energy E is which we minimize by choosing static configurations, φi = 0. Then we re-arrange E as follows: where η = ±1 is independent of i and where for static fields φ i (x) we denote W(x) = W(φ 1 (x), . . ., φ N (x)), and hence Finite energy requires that both φ ′ i , W φi → 0 as |x| → ∞, as follows from (12), and so W ′ (x) → 0 as |x| → ∞.Denote the values of W(x) as x → ±∞ by W ±∞ , then from (13) we have the lower bound E ⩾ |W ∞ − W −∞ |, which is attained if and only if the following Bogomolny equations hold: with the appropriate sign η, independent of i.Since the sign can be reversed simply by replacing x → −x we henceforth set η = 1.
Because any solution of the system φ ′ i = W φi minimizes the energy, the second-order equations are automatically satisfied as we verify directly.The static equations for the potential (2) are given by: where W φi,φj ≡ ∂ 2 W/∂φ i ∂φ j , and from (14): and so (15) is satisfied.In general, however, φ ′2 i − W 2 φi is not a constant of integration, meaning that there exist finite-energy solutions of the second-order equation (15) which do not satisfy (14).The following is nevertheless constant, regardless as to whether ( 14) is satisfied.

Independent fields and the reduction for N = 3
Consider now W as defined by (3), for some function w yet to be specified, then the N firstorder equation ( 14) read: where the two summations on the right-hand side can be combined into a single sum if w is either even or odd.It follows, by summing (18) over i and interchanging i, j under the double sum, that N i =1 φ ′ i = 0, and so we deduce that N i =1 φ i is a constant of integration, which is a consequence of the fact that W as defined by (3) depends only on the differences φ j − φ i .By means of the translation invariance φ i → φ i + φ 0 for any constant φ 0 we can set N i =1 φ i = 0. Of the N fields φ i , only N − 1 are therefore independent and so we could choose, for example, the N − 1 differences φ i+1 − φ i for i = 1, . . .N − 1 as independent fields.For N = 2 in particular, the independent fields are φ 2 + φ 1 (which is constant) and φ 2 − φ 1 , and so the system is equivalent to that of a single-scalar field φ = φ 2 − φ 1 .Hence we do not further consider the case N = 2.
For N = 3 there are two independent fields, which is of interest because soliton models with two scalar fields have been investigated many times, usually without the benefit of an explicit superpotential [7][8][9][10][11][12], although with some exceptions [13][14][15][16][17][18]31].The models that we consider all have explicit superpotentials by construction, but do not appear to correspond to any previously studied, such as those in [14,31].We transform to new variables u, v, σ defined by then we can write W and V, which depend only on the differences φ j − φ i , as functions of u, v only, i.e. we have W = W(u, v), V = V(u, v).The static energy density is given by showing that the σ field decouples from the system.In order to minimize the energy, we set σ ′ = 0.In terms of u, v, the Bogomolny equation ( 14) read as follows from , where we write W as a function of either φ 1 , φ 2 , φ 3 or u, v as appropriate.Then V as defined by ( 2) is given as a function of u, v by In our investigation of N = 3 models we therefore start with an explicit superpotential W = W(u, v), from which we construct V according to (22), then solve the two Bogomolny equation (21), either exactly or numerically.The fixed points are determined by W u = 0 = W v , and are zeroes of V.
For comparison with other soliton models with two scalar fields, we diagonalize the derivative terms in the energy density (20) by defining new fields ψ 1 , ψ 2 according to where Then the energy density takes the conventional form where , and similarly W can be expressed as a function of ψ 1 , ψ 2 .Polynomial superpotentials proposed in previous studies [14], with the energy density in the standard form (25), do not correspond to the N = 3 cases of the φ 4 and φ 6 models that we investigate, when transformed according to (23).

Soliton stability and the linearized equations
Because solutions of the first-order system (14) minimize the energy, they are stable against any perturbations which maintain the asymptotic values of the fields φ i .If we define the N currents J µ i = ε µν ∂ ν φ st i where φ st i (x) is a static finite-energy solution of (14), then J µ i is a conserved current.The conserved charge is the difference between the asymptotic values, namely φ st i (∞) − φ st i (−∞) which remains unchanged with respect to finite-energy perturbations.We refer to any static finite-energy configuration φ st i (x) for fixed i as a soliton if this difference is positive, otherwise as an antisoliton.Having defined W(x) = W(φ st 1 (x), . . ., φ st N (x)), it follows that when evaluated at φ i (x) = φ st i (x) and so, for solitons, W(x) is an increasing function of x.The total energy |W ∞ − W −∞ | is therefore nonzero unless the configurations φ st i (x) are zeroes of V for all i, corresponding to vacuum states.
We investigate linear stability of the static solitons φ st i (x) by writing which satisfies the Euler-Lagrange equations ∂ µ ∂ µ φ i = −V φi .These linear fluctuations ψ i also represent the elementary excitations of the system, and so determine how the elementary particles interact perturbatively with the soliton.To first order in ψ i we have: where U = (u ij ) denotes the N × N symmetric matrix of second derivatives of V evaluated at the static solution φ st i (x) for i = 1, . . .N, i.e.
Denote by ψ the column N-vector with components ψ i , then we have the linear system of equations ψ − ψ ′ ′ = −Uψ.
In order to test the stability of the soliton system we suppose that ψ has a time-dependence of the form ψ(x, t) = e Λt φ(x) for some N-vector function φ of x, where Λ is a real N × N matrix.It follows that ψ = Λ 2 ψ.We wish to establish that the eigenvalues of Λ, which can be complex, have real parts which are negative or zero, otherwise the perturbation is unstable.We write the linear equations ψ − ψ ′ ′ = −Uψ as Hψ = −Λ 2 ψ, where where I N denotes the N × N identity matrix.

Factorization of H
The matrix operator H can be factorized in the form H = S † S where where M = (m ij ) is the real N × N symmetric matrix of second derivatives of W evaluated at the static solution φ st i (x) for i = 1, . . .N, explicitly: In order to show this, from (31) we have and so we wish to show that U = M ′ + M 2 .For a given potential U, this is a matrix Riccati equation which determines M, however in our case we can calculate M directly, since W is known explicitly.From the definition (2) of V we obtain from which using the fact that W φi,φj is evaluated at the static solution for which φ ′ i = W φi .Next, we evaluate (34) and (35) at φ i = φ st i to obtain, in matrix element form, and so U = M 2 + M ′ as required.
The matrix operators S, S † in (31) are hermitean conjugates in the usual L 2 inner product for N-vectors ψ, and so H = S † S is a positive semi-definite matrix operator.The equations Hψ = −Λ 2 ψ imply that −Λ 2 is a positive semi-definite matrix, and since H is also hermitean, Λ 2 is a symmetric matrix.We can express Λ = Σ + Ω as the sum of real symmetric and antisymmetric matrices Σ, Ω (respectively), then Λ 2 = Σ 2 + Ω 2 + ΣΩ + ΩΣ.On taking the transpose we obtain is required to be positive semi-definite, we must have Σ = 0, otherwise −Σ 2 can have real eigenvalues which are negative and large in magnitude, contradicting the property that −Σ 2 + Ω T Ω is positive semi-definite.On setting Σ = 0, we deduce that Λ is antisymmetric and is therefore an element of the Lie algebra of the orthogonal group O(N), and so e Λt = e Ωt is an N × N rotation matrix which preserves the length of ψ.The perturbations ψ(x, t) = e Λt φ(x) are therefore bounded as t increases, which establishes linear stability.
A state with neutral stability arises from the translational invariance of the system.Denote by ψ 0 the N-vector with components dφ st i /dx, then Sψ 0 = 0, as follows from the property φ ′ ′ i = N j =1 W φi,φj φ ′ j for any static solution φ st i .Hence Hψ 0 = 0 and so ψ 0 is a ground state, known as the zero mode, for the system.
The factorization property of Schrödinger Hamiltonians such as H in the form H = S † S has been widely discussed in many contexts, together with related concepts of supersymmetry, see [32][33][34][35] to which we refer for further development.In our case, we have the commutation relations [S, S † ] = −2M ′ from which one can develop the algebra of S and S † , regarded as ladder operators.A well-studied example is the Rosen-Morse potential which can be analyzed by algebraic means [36], and appears in our context for the φ 6 system, see section 9.3.The factorization H = S † S for single-scalar fields is well-known [2] (section 2.1.2) and has been generalized to coupled scalar fields [9,14], although usually the eigenfunctions of H are obtained directly, where possible.
For the explicit superpotential W as defined in (3), for the case where w is even for which w ′ is odd with w ′ (0) = 0 and a ij is symmetric, we have The matrix elements m ij = W φi,φj are therefore given by and M = (m ij ) has the property that the rows and columns each sum to zero.Hence M, and therefore also U, always has an eigenvector {1, 1, . . .1} T of length N with a corresponding zero eigenvalue.Generally M(x) must be evaluated numerically, however there are cases where explicit solitons are known, for which we can evaluate M and hence obtain U = M 2 + M ′ as an explicit function of x.An example is the φ 6 model for N = 3 discussed in section 9.3, where we find an explicit form for M(x) which can be diagonalized by means of a constant orthogonal matrix, enabling us to analyze the bound states of the elementary excitations with the soliton.

Elementary excitations and particle masses
The elementary excitations of the quantum system defined by the Lagrangian (1) are quantum particles with masses µ i that depend on the coefficients a ij , and propagate on the real line, interacting over the complex system by means of the potential V.These meson excitations appear in the usual way as perturbations of the vacuum states, and also determine the asymptotic form of the solitons.For nonzero masses µ i the soliton configurations approach their vacuum values exponentially quickly, and so we can construct numerical solutions for |x| < L, where L is sufficiently large to allow for the exponential fall-off to the vacuum values.
Let us restrict our investigations here to models in which w is an even function, then W φi is given by (37), where without loss of generality we set a ii = 0. Consider only the elementary excitations about the vacuum state φ i = φ 0 for all i, where φ 0 is constant, then φ j − φ i = 0 for all i, j and so W φi = 0.The state φ i = φ 0 is therefore a zero and hence an absolute minimum of V. We perturb about this vacuum according to φ i (x, t) = φ 0 + ψ i (x, t), then the linearized equations are given, as above, by (29), evaluated at φ i = φ 0 .From (36) we have U = M 2 where M is the matrix of second derivatives of W as given in (38), evaluated at the vacuum state φ i = φ 0 , for which M ′ = 0.
Since U = M 2 is real and symmetric it can be diagonalized with non-negative real eigenvalues µ 2 i , and so the equations ∂ µ ∂ µ ψ = −Uψ reduce to N independent Klein-Gordon equations.The elementary excitations about φ i = φ 0 are therefore scalar particles of mass |µ i |, where µ i is an eigenvalue of M. The explicit matrix elements of M may be deduced from (38), namely: recalling that we have set a ii = 0.As before, M has an eigenvalue of zero with the eigenvector ψ = {1, 1, . . .1} T , however these excitations merely reflect the fact that the sum N i =1 φ i decouples from the system, as shown in (20) for N = 3, and takes a constant value for any soliton solution.
As an example, for N = 3 we have For the special case a 23 = a 12 , which we consider in section 9.2 in connection with φ 6 models, the two nonzero eigenvalues of M are (µ 1 , µ 2 ) = 2w ′ ′ (0)(3a 12 , a 12 + 2a 13 ), and so the particle masses based on the vacuum state For the periodic systems of solitons discussed in section 5.1, in particular the Kuramoto system (47) for which we choose a ij = 1 for i ̸ = j, the nonzero masses are all equal to N but generally, for the Kuramoto model ( 42), solitons exist for any symmetric couplings a ij with the meson masses |µ i | given by the eigenvalues of M. We do not consider special cases where, by specific choice of coupling a ij , one or more of the masses |µ i | is zero, as these may correspond to a potential V which does not support solitons.

Synchronization in complex systems
The first-order Bogomolny equations which allow stable, static finite-energy solitons in a system of scalar fields are given by (14), namely φ ′ i = W φi , where φ i satisfies the boundary conditions φ ′ i → 0 as |x| → ∞ for all i = 1, . . ., N. Such systems have long been studied in an entirely different context as evolution equations for complex systems of interacting oscillators, as models of collective behavior, in particular as models of synchronization phenomena, where x is regarded as a time variable.Here, we view a complex system as a network consisting of N discrete nodes, connected in pairs by discrete links, which can be represented by an N × N adjacency matrix (a ij ).In the simplest case we set a ij = 1 for connected nodes i, j, otherwise a ij = 0. Generally a ij = a ji if there is no preferred direction of coupling.At each node there is a dynamical system describing a phase oscillator, with a corresponding variable θ i = θ i (t) which describes the state of the system.We wish to analyze the evolving system arising from the local dynamics combined with nonlinear interactions over the network of connected nodes.
Synchronization processes, in which the dynamics of the nodes are highly correlated, are common in nature and appear in many different contexts including biology, ecology, sociology as well as physics, chemistry and engineering.We consider in particular phase and frequency synchronization in the context of the widely-studied Kuramoto model.For excellent reviews of this vast topic, including detailed discussions of network properties, nonlinear dynamical systems, stability, synchronization, with many further references, see [25,[37][38][39][40].

The Kuramoto model
A model of the transition to synchronous behavior which is mathematically tractable, and yet sufficiently complex to display a large variety of synchronization patterns, is the Kuramoto model [22] in which the dynamics are governed by the system of first-order evolution equations where the real coefficients a ij determine the strength of interaction between the ith and jth nodes and ω i is the natural frequency of oscillation at the ith node.If all interactions are switched off by setting a ij = 0, the nodes oscillate independently with frequency ω i , but if we set a ij = K/N for all i, j, where the coupling constant K > 0 is sufficiently large, then the system synchronizes in the sense that | θj − θi | → 0 as t → ∞ for all pairs i, j, i.e. all nodes oscillate at a common frequency, known as phase locking or frequency synchronization.For rigorous results we refer to [41][42][43][44].
We consider the Kuramoto model only for identical frequences ω i = ω, in which case it is known that for a ij = K/N the system completely synchronizes from generic initial values, meaning that |θ j − θ i | → 0 as t → ∞ for all pairs i, j [43,44].For ω i = ω we can transform to a rotating frame by replacing θ i → θ i + ωt in (41), which means that we can set ω = 0. We are interested in the Kuramoto model therefore in the form It is usual to quantify the overall synchrony of this system by means of a time-dependent order parameter r defined by which is the centroid of a set of N points e iθj distributed on the unit circle in the complex plane, for some time-dependent angle δ.We have 0 ⩽ r ⩽ 1, where r = 1 corresponds to a completely synchronized system in which θ i = θ j for all nodes i, j.The value of r(t), as the system evolves, conveniently indicates whether the system has attained its asymptotic value, and whether synchronization has occurred with r → 1.
We can write (42) as a gradient system by defining the potential where a ij = a ji , then (42) can be written θi = W θi , where −W is well-known as a Lyapunov function for the Kuramoto model ( 42) [45,46].Although a gradient formulation is not a necessary requirement for the modeling of synchronization phenomena, it leads to stability properties of the steady states [47], and is a common feature of the soliton systems that we construct.
Because the system ( 42) is given by θi = W θi , which is identical in form to the Bogomolny equation ( 14) with the time and spatial variables interchanged, we can use known results of the Kuramoto model to determine properties of the corresponding solitons.One difference is that models such as the Kuramoto system are usually defined on the unit circle, i.e. the phase space for the N oscillators is the N-dimensional torus, whereas for the sine-Gordon system, generalized to N scalars, the angles θ i are allowed to take values on R N .Complete synchronization, in which |θ j − θ i | → 0 as t → ∞, should therefore be re-written for solitons as the property |e iθj − e iθi | → 0 as x → ∞.

Interplay between time and spatial variables
Static soliton solutions can be viewed, by analogy with mechanical systems, as trajectories which evolve for −∞ < t < ∞ from an initial point located at a maximum of −V, to a final point at an adjoining maximum of −V (Coleman [4], chapter 6).The same analogy holds for systems of static solitons, where we regard θi = W θi as a set of time evolution equations, under which θ i (t) evolves from specified initial values at a fixed initial time, to final asymptotic values at large times.
Instead of choosing the initial values of θ i to lie exactly at a maximum of −V, however, we select an initial point t 0 ∈ (−∞, ∞), which by translational invariance can be t 0 = 0, and we specify the values θ 0 i = θ i (0) at time t = 0 through which all trajectories θ i (t) must pass.Hence θ 0 i , i = 1, . . .N denotes the set of initial values which uniquely determine the solution of θi = W θi for all t > 0. These same initial values θ 0 i also uniquely determine the solution of θi = −W θi (reversing the sign of t) where, in effect, we allow the system θi = W θi to evolve for t < 0 to final asymptotic values at large negative times.In this way we construct the unique solution θ i (t) of θi = W θi for all −∞ < t < ∞ which passes through the specified points θ 0 i at t = 0.
For numerical solutions, labeling now the spatial variable as x, we randomly generate values θ 0 i ∈ [0, 2π] and from these integrate the system θ ′ i = W θi for −L < x < L in both the positive and negative directions for x, for sufficiently large L, to obtain θ i (x) for all |x| < L. For the models that we consider, θ i (x) achieves its asymptotic values exponentially quickly due to the fact that the elementary excitations have nonzero mass, and so L can be relatively small.
The question which arises then is whether the solutions θ i (x) constructed in this way have finite energy.It is sufficient to demand that θ ′ i (x) → 0 as |x| → ∞ since, because θ ′ i = W θi is satisfied, the constant asymptotic values of θ i are then automatically zeroes of V = 1 2 i W 2  θi , and so the solution has finite energy.Numerically, we find in all cases that the condition θ ′ i (x) → 0 as |x| → ∞ is indeed satisfied, and so the above numerical procedure generates finite-energy soliton configurations for any N, for any given values θ 0 i .For the Kuramoto model ( 42) with global coupling, rigorous results have been established [41][42][43][44] which prove these numerical observations, as we discuss in the following section 5.The requirement of finite energy is generally much weaker than that of synchronization, since the Kuramoto model synchronizes only for positive (attractive) couplings when evolving for t > 0, whereas for negative (repulsive) couplings, when the system evolves in the t < 0 direction, the asymptotic values of θ i are uncorrelated.The corresponding soliton configurations nevertheless have finite energy.

Periodic systems of solitons
We construct a 2π-periodic potential V for a system of N scalar fields θ i starting with the sine-Gordon model ( 4), for which w(θ) = K cos θ.According to the formula (3) we define the superpotential W by (44), with K = 1  2 , where the coupling coefficients a ij are symmetric in i, j.We have therefore W θi = N j =1 a ij sin(θ j − θ i ), and so according to (2) the potential is given by: which expresses V as a sum of terms with both two-body and three-body interactions, where the three-body terms arise from the final summation over i, j, k.The second-order system θ ′ ′ i = V θi is solved by the first-order Bogomolny equations θ ′ i = W θi which allows only two-body interactions.
As an example, we plot −V for N = 3 in figure 1, as a function of the two independent variables u = θ 2 − θ 1 , v = θ 3 − θ 2 for the case of global coupling a ij = 1, where the explicit expression for V is given in (52).We plot −V rather than V because by mechanical analogy, solitons correspond to particle-like trajectories between adjacent maximums of −V.The zero of V at the origin is shown in green while the twelve surrounding zeroes are located at the points Figure 1.The 2π-periodic potential −V for N = 3 as defined in ( 52), plotted as a func- showing a zero of V at the origin (green) with surrounding zeroes in red, blue.
for m, n = −1, 0, 1 (marked in red, blue respectively).There are also the 2π-periodic copies of these zeroes.We expect therefore that twelve distinct solitons exist which interpolate between the origin and the adjacent zeroes of V, as we discuss further in section 5.2.

Kuramoto solitons with all-to-all coupling
Although W is defined for arbitrary symmetric coupling coefficients a ij , we now choose global (all-to-all) coupling by setting a ij = 1, then the Bogomolny equations read which we write as where r, δ are defined in (43).We see that all fixed points satisfy either sin(δ − θ i ) = 0 or r = 0, and it is known that the stable fixed points satisfy either r = 1 (for positive coupling) or r = 0 [41][42][43][44]48], where by 'stable' we refer to the asymptotic stability of θ i (x) as x → ±∞.The two cases r = 1 and r = 0 arise as follows: (a) r = 1 implies that θ i = δ for all i and so all points are co-located, corresponding to complete synchronization for the system (47), with θ i = θ j for all pairs i, j, modulo 2π.It is evident from (45) that such points are zeroes of V.There also exist fixed points of (48) for which θ i = δ + π for one or more values of i which are, however, unstable and so do not appear in any soliton configurations.Such bipolar states are discussed in [49] (section 6.1).(b) There are many stable fixed points for which r = 0, even for small values of N, all of which are zeroes of V. We refer to [50] and further references therein for the analysis and classification of such fixed points, although the application to our soliton models is for identical frequencies ω i = ω only.
We construct solitons for any N by solving (47) numerically over the interval −L < x < L, from randomly generated initial values θ 0 i at x = 0, in both the positive and negative x directions.Because the initial values are random, we do not encounter any fixed points which are the vacuum states of the system.For x > 0, it is known [41][42][43][44]48], and is easily confirmed numerically, that for generic initial values θ 0 i we always obtain r → 1, i.e. the system completely synchronizes, and so the static soliton at each node approaches a constant asymptotic value δ (modulo 2π).For rigorous statements and proofs we refer to section 4, theorem 4.1 in [43], and to theorem 4.4 in [44].For x < 0 we always obtain r → 0, and the solitons approach distinct values.Exactly which of the many such zeroes the system chooses depends on the initial values θ 0 i .For a rigorous statement and proof of the property that r → 0 we refer to [44], theorem 4.6.
The superpotential (44) for global coupling a ij = 1 is given by W = 1 2 N i,j =1 cos(θ j − θ i ) and from (43), by multiplication with its complex conjugate, we obtain: Hence W = 1 2 N 2 r 2 where we regard W, r each as functions of x.Both W, r are therefore increasing functions of x, as follows from (26) with N 2 rr ′ = 2V.Denote by W ±∞ the values of W(x) as x → ±∞, as before, then the total energy E of the system in ( 12) is 2 N 2 = E, which provides a numerical check on solutions.As an example, we have plotted the N soliton configurations θ i (x) in figure 2 for |x| < 2 for N = 15.Numerically there is no difficulty in generating configurations for much larger values of N, except that individual trajectories become less distinguishable on the plot.The randomly generated initial values θ 0 i are as shown on the vertical axis at x = 0.The system completely synchronizes from these values for x > 0 with r → 1, and for x < 0 the trajectories settle into a configuration for which the asymptotic values θ i (−L) are uncorrelated, but always such that r → 0. We make the following further observations: (i) The specific asymptotic values θ i (−L) vary according to the choice of initial values θ 0 i , however the total energy E = 1 2 N 2 remains unchanged; (ii) The N initial values θ 0 i are indirectly related to the locations of the N solitons, i.e. small perturbations in θ 0 i lead to small adjustments in the relative soliton positions; (iii) The configurations θ i represent a mixture of solitons and antisolitons, each positioned at varying locations as measured by the maximum value of the corresponding energy density (θ ′ i ) 2 .These energy densities are plotted in figure 3, showing large variations in the distribution of energy across the network; (iv) The trajectories in figure 2 never cross, and so if we plot these as points on the unit circle, then collisions never occur.This property is special to the system (47) and can be understood as a consequence of the partial integrability discussed in section 6.More directly, for N ⩾ 4 and for each set of distinct nodes i, j, k, l we define the cross ratios then it follows directly from the Kuramoto equation ( 47) that C ′ ijkl = 0, i.e.C ijkl is independent of x and so can be evaluated at x = 0 in terms of the initial values θ 0 i .Cross ratios and their properties have been investigated in detail for the Kuramoto model, see [51] section V.A, in particular there are N − 3 independent constants of integration.If the initial values θ 0 i are all distinct, then C ijkl is nonzero for all x and so any two nodes i, k cannot coincide, i.e. collisions never occur.Similarly, if any two nodes are equal at x = 0 then they remain locked together for all x.
Generally, it is straightforward to solve the Kuramoto system (42) numerically for arbitrary symmetric couplings a ij , for |x| < L. We can, for example, randomly generate the entries of the symmetrized adjacency matrix A = (a ij ) in the interval (−1, 1), and solitons are again created, but without simple asymptotic behavior such as r → 1 or r → 0, i.e. synchronization does not occur even though the configurations represent acceptable finite-energy solitons.We briefly discuss such models in section 5.3 for N = 3.

N = 3 solitons for the Kuramoto model
Solutions to the system (47) can be determined numerically for very large N, however we are also interested in small N for which some exact solutions can be determined.Solitons for N = 2 are equivalent to the sine-Gordon soliton (5), due to the fact that θ 1 + θ 2 is constant.For N = 3 we wish to verify general properties of the system, as well as those specific to N = 3.Because solitons such as those depicted in figure 2 are found numerically by choosing random initial values θ 0 i , we do not encounter special cases such as sine-Gordon solitons, which arise when two or more of the initial values coincide.Previous work on the N = 3 Kuramoto model includes [40], section 4.1.2and [52], for oscillators with arbitrary natural frequencies, whereas we have identical oscillators with ω = 0. Define , and the superpotential is given up to additive constants by and V has the explicit form: which is plotted in figure 1. V satisfies symmetries which follow from the cyclic properties of θ 1 , θ 2 , θ 3 such as: V is non-negative by construction, with zeroes determined by W u = 0 = W v , as given explicitly in (46) where m, n are integers.The Kuramoto equation ( 47) reduce to as follows from (21).Having solved for u, v we can recover the three angles θ 1 , θ 2 , θ 3 by using the fact that σ = θ 1 + θ 2 + θ 3 is constant.
Figure 4 shows a contour plot of V where the zeroes are marked corresponding to those in figure 1. Solitons interpolate between one of these zeroes and the origin, due to the fact that as x → ∞ the system completely synchronizes, i.e. we have (u, v) → (0, 0).(We assume here that any such solution is translated by a multiple of 2π so that the trajectory terminates at the origin).As x → −∞, however, (u, v) takes an asymptotic value corresponding to one of the twelve zeroes surrounding the origin.For any solution (u, v) of ( 54) we can verify directly that the following are also solutions: together with combinations of these symmetric transformations.Again, this follows from the cyclic properties of θ 1 , θ 2 , θ 3 .Hence, for any soliton interpolating between (π, 0) and (0, 0), there are further solitons interpolating between any one of (−π, 0), (0, π), (0, −π), (π, −π), (−π, π) and the origin.A similar remark applies to solitons interpolating between points such as ± 2π 3 , − 4π 3 , i.e. those marked in red in figure 4, and the origin.
The order parameter r in ( 49) is given by and takes the value r = 1 at the origin, and r = 0 at the points marked in red in figure 4, consistent with the asymptotic values for soliton configurations θ i for general N.For the zeroes of V marked in blue in figure 4, however, we obtain r = 1 3 , which is a case not covered by the previous analysis, and occurs only for exactly specified initial values, as described in the next section.Hence there are essentially two types of soliton for N = 3, those interpolating between the origin and either the red or blue zeroes of V, as shown in figure 4.

Reduced Kuramoto model solutions for N = 3
There are special cases for which we can find exact solutions to (54) which interpolate between (0, 0) and one of the surrounding twelve zeroes of V depicted in figure 4, which demonstrates that such solitons do indeed exist.The two types are given by: (i) If v = 0 then u ′ = −3 sin u, which has the solution u = 2 arctan e −3(x−x0) , for which u → 0 as x → ∞ and u → π as x → −∞.This is a sine-Gordon soliton interpolating between the points (0, 0) and (π, 0) in the (u, v) plane, whereas −u interpolates between (0, 0) and (−π, 0).By symmetry, there are similar sine-Gordon solitons for u = 0 and u + v = 0.By setting u = 0, for example, we ensure that the variables θ 1 , θ 2 are locked together, and so the N = 3 equations reduce to the N = 2 equations.These solutions do not appear for the system (54) for random (generic) values (u 0 , v 0 ) unless we specify exact initial values, such as v 0 = 0 which implies v = 0. (ii) For v = u (54) is satisfied by u ′ = − sin u − sin 2u which may be integrated to obtain where x 0 is the constant of integration.As x → ∞ we have either sin u = 0 or cos u = −1, which again implies sin u = 0, and so u = nπ for integers n.As x → −∞ we have cos u = − 1 2 which implies u = ± 2π 3 + 2mπ.Hence this soliton interpolates between (0, 0) and one of the points ±( 2π 3 , 2π 3 ), shown in figure 4 as red dots on the line u = v.There are similar solitons obtained by setting 2u + v = 0 and u + 2v = 0.
We have so far chosen all-to all coupling a ij = 1, however solitons exist for general a ij , although properties such as synchronization are much less understood.For N = 3 some exact solutions may again be obtained.The superpotential which generalizes ( 51) is: and the Kuramoto equation ( 42) read: As the properties of this system, such as the fixed points, are rather complicated for general a ij , let us choose a 23 = a 12 then we can solve (59) by setting v = u, provided that The fixed points satisfy either sin u = 0, or cos u = −a 12 /(2a 13 ), which exists only if |a 12 | ⩽ 2|a 13 |.The solution of ( 60) is given by: and reduces to (57) on setting a 12 = 1 = a 13 .From this expression we can determine the asymptotic behavior of cos u, depending on the relative values of a 12 , a 13 .We observe that there are special cases, such as a 12 = 6a 13 , for which an explicit expression for cos u can be obtained.In this case the linearized fluctuations about the soliton can be calculated explicitly, as explained in sections 3.3 and 3.4 by evaluating the 3 × 3 matrix M(x) in (32).We do not pursue this further here, but perform such an evaluation for the φ 6 system in section 9.3.

Partial integrability of the Kuramoto system
The Kuramoto system (47) of N equations with global coupling a ij = 1 has the remarkable property that it can be reduced to just two independent equations by means of a partial integration, implemented by the Watanabe-Strogatz (WS) transform.This occurs because the model has N − 2 constants of motion, which become identically constant by transforming to a system consisting of a single complex variable ζ.Complete information about the solutions and properties of the system for any N ⩾ 3 can therefore be deduced from the properties of ζ.The WS transform was first outlined by Watanabe and Strogatz in [53,54] in trigonometrical form, but was later seen to be equivalent to the Möbius group of transformations which map the unit disk {z ∈ C : |z| < 1} onto itself [51].This property of partial integrability has been extensively investigated for the Kuramoto model and its generalizations to higher dimensions [28,51], and has applications to neuroscience and biological systems amongst others [55], because it allows a substantial simplication for large N systems.Here we outline the application to solitons, in particular how the asymptotic behaviour of ζ leads to finite-energy solitons.

The Möbius map and partial integration
Let z i = e iθi then we can write the Kuramoto equation ( 47) in the Riccati form where It follows from ( 62) and its complex conjugate that z i z i is constant in x and so the N equation ( 62) are equivalent to the Kuramoto system (47) provided we specify that z i z i = 1 at some point x = x 0 .We implement the Möbius map by substituting for z i according to where we have replaced the N variables θ i by N variables ξ i plus one complex variable z (not to be confused with z i ), and so the new variables z, ξ i are not uniquely determined by θ i .After substituting (64) into (62), where we regard Z, Z each as fixed functions of x, and clearing the fractions on both sides, we obtain an expression which is quadratic in e iξi .We set the coefficients of e iξi and e 2iξi each to zero, to obtain the conjugate equations Next, we back-substitute for z ′ and z ′ into the remaining equation, which we find is satisfied provided that i ξ ′ i = −z Z + z Z. Since the right-hand side of this equation is independent of i, we must have ξ i (x) = α(x) + ξ 0 i for real constants ξ 0 i and some real function α(x) satisfying i α ′ = −z Z + z Z.It is here that the partial integration occurs, where ξ 0 i are the N constants of integration.Since we can arbitrarily shift constants between α and ξ 0 i , we choose the initial value α(0) = 0. We also set z(0) = 0, which corresponds to a choice of base point as discussed in [56], then by evaluating (64) at x = 0 we obtain ξ 0 i = θ 0 i , i.e. the constants ξ 0 i are precisely the given initial values for the Kuramoto system (47).
Let us rewrite (64) by defining ζ = z e −iα , then the WS transform reads and (65) becomes together with the complex conjugate equation, where We also have i α ′ = −ζ F + ζ F, however we can obtain a closed form expression for α using the fact that N j =1 θ j is a constant of integration, i.e.
N j =1 θ j = N j =1 θ 0 j , which we exponentiate.From (66) we obtain: which is equal to its value at x = 0, hence In summary, in order to solve the Kuramoto system (47) given θ i (0) = θ 0 i , for any N ⩾ 3, we first solve (67) for ζ(x) with ζ(0) = 0, then calculate e iα from (70), and so determine e iθi for all i = 1, . . .N from (66).The number of real equations to be solved has therefore been reduced from N to only 2, which is essentially due to the fact that there are N − 3 independent constants of integration ( 50), together with the additional constant N i =1 θ i .Since solitons are constructed in terms of the differences θ i − θ j , we require only the formula which follows from ( 66): independent of α.We refer to [57] for a summary of the Watanabe-Strogatz theory, and for developments such as the connection with the Ott-Antonsen theory in the N → ∞ limit, and to [51] for a detailed description of the action of the Möbius group on the unit circle, and to [56,58] for derivations of the properties of ζ.

Solutions ζ of the reduced system
One consequence of ( 67) is that ζ remains inside the unit circle on the complex ζ-plane for all x, since it follows from (67) and so ζ(x) cannot cross the unit circle.Complete synchronization as x → ∞ corresponds to |ζ| → 1, for which ζ → e iβ , then from (66) we obtain e iθi → e i(α+β) which is independent of i, i.e. all trajectories completely synchronize, as the example in figure 2 shows.As x → −∞, however, ζ approaches a fixed point of ( 67) which is a zero of F as defined in (68), see the discussion in [58].For any such configuration we have r = 0, as follows from the definition (43).
In figure 5   Exact solutions for ζ(x) and hence θ i generally do not exist for arbitrary N, except for certain restricted initial values θ 0 i for which two or more nodes are locked together.If θ 0 i = θ 0 j for some i, j then from (71) we have e i(θi−θj) = 1, hence trajectories which are co-located at x = 0 remain so for all x.If θ 0 i = θ 0 j for all i, j then θ i = θ j for all x, which represents the vacuum state θ i (x) = θ 0 for some constant angle θ 0 , corresponding to a completely synchronized configuration.The solution ζ of (67) in this case is ζ(x) = e iθ 0 tanh(Nx/2), which takes the asymptotic values ζ ±∞ = ±e iθ 0 .
For even N with N = 2M, and initial values θ i (0) = θ 0 for i ∈ S + = {1, 2 . . .M} and θ i (0) = −θ 0 for i ∈ S − = {M + 1, . . .2M}, for some constant angle θ 0 , all nodes θ i take one of two possible trajectories, either a soliton or an antisoliton, according to whether i ∈ S − or S + .We can find an exact expression for ζ, which in this case is real, by writing F as the sum of two terms which can be combined to obtain (67) in the form: The solution with which has the asymptotic values (for cos θ 0 > 0): which are fixed points of (73).From (71) we obtain: where i ∈ S + , j ∈ S − , and directly from (66), with α = 0, we obtain the soliton or antisoliton configurations where the ± sign corresponds to i ∈ S ± .We have the asymptotic values θ i (x) → 0 as x → ∞ and θ i (x) → ± π 2 as x → −∞, consistent with previous observations.These solutions are essentially the sine-Gordon soliton or antisoliton (5), due to the restricted initial values which we have imposed at x = 0, and although they can be derived directly from the Kuramoto system (47), we have proceeded as above as an example of the properties of the ζ-function.
Finally, we consider properties of ζ for the N = 3 Kuramoto equations discussed in section 5.3.The formulas ( 66) and (70) for ζ are valid for N = 3 except that the conserved crossratios (50) do not exist, and so there are two independent real equations, whether expressed in terms of ζ or (u, v).We regain the two equation ( 54) by evaluating e iu and e iv explicitly in terms of ζ and ζ and the initial values θ 0 i using (71), where Then, by expressing ζ, ζ in terms of u, v, we can derive the u, v equations ( 54) from ( 67) for ζ, which is best accomplished using a computer algebra system.
Since the resulting formulas are rather lengthy, let us consider only the restricted initial values u 0 = v 0 , which imply that v = u, in which case (54) then the ζ-equations ( 67) and ( 72) each reduce correctly to u ′ = − sin u − sin 2u, which has the solutions (57).Evidently, ζ(0) = 0 follows from u(0) = u 0 .If u 0 = 0 then u(x) ≡ 0 for all x, which is a vacuum state, but otherwise u 0 and u are of the same sign, and so sin u 0 sin u > 0, which from (79) implies correctly that |ζ| < 1.For x → ∞ we have |ζ| → 1 which corresponds to sin u → 0, as noted for the solution (57) but as x → −∞, ζ approaches a fixed point ζ −∞ of F, i.e. cos u → − 1 2 as (80) shows correctly.The corresponding value ζ −∞ is given exactly by (78) in which u → 2π 3 .In summary, we have outlined how the WS transformation can be applied to the Kuramoto soliton system (47) to derive soliton configurations in which properties such as finite energy follow from the asymptotic behavior of ζ.

The antisymmetric Kuramoto model
The superpotential for the sine-Gordon model can be expressed as either an even or odd function of θ, as shown in ( 4) and ( 5), each of which extends to interacting systems of N scalar fields in different ways.Whereas the even extension leads to the well-known Kuramoto model, the odd extension leads to Kuramoto-type equations but with antisymmetric couplings, which nevertheless allow conventional solitons.We consider such systems here only briefly in order to point out the main differences compared to those of the usual Kuramoto model, essentially that solitons with antisymmetric couplings repel each other, whereas for symmetric couplings the configurations are attractive for x > 0.
We choose therefore as in (6), w(θ) = −K sin θ and define the superpotential according to (3): where a ij = −a ji , and so (with K = 1 2 ) we obtain the system of N first-order equations These equations have been previously studied [59] (section 7), as well as a similar system with symmetric all-to-all coupling [54] (section 4).The 2π-periodic potential is defined by (2) with the zeroes, which are absolute minima, given by W θi = 0 for all i.The system (82) differs from the Kuramoto model ( 41) in several ways: (i) Complete synchronization cannot occur because the configuration with θ j = θ i for all i, j is not a fixed point; (ii) (82) is invariant under θ i (x) ↔ −θ i (−x), i.e. properties of the system in the negative x direction, such as its asymptotic values, are essentially the same as those for x > 0; (iii) the coupling coefficients a ij are antisymmetric, as follows from the definition of W in (81), and so properties of the system depend on the allocation of signs to a ij , even for the simplest case |a ij | = 1 for all i ̸ = j.A convenient coupling which we choose henceforth is i.e. the adjacency matrix A = (a ij ) has values equal to 1 above the diagonal, −1 below the diagonal, and zero on the diagonal.
Upon solving (82) numerically for a ij = ε ij , for randomly generated initial values θ 0 i ∈ (−π, π), we find that synchronization always occurs in the sense that θ i → θ 0 − i π N for all i as x → ∞, for some constant angle θ 0 , and similarly The asymptotic values of the nodes are therefore equally spaced on the half-circle, a configuration which is sometimes referred to as a splay state, see the discussion of properties of splay states and further references in [60].The fact that the states θ i = θ 0 ± i π N are fixed points is proved in [59], section 7.1, using the property The asymptotic value of the order parameter r, defined in (43), for these fixed points is given by r = (N sin π 2N ) −1 , and the asymptotic value of W is as follows from [59], equation (54), and we also have Solitons are generated numerically as before, by choosing random initial values θ 0 i at x = 0, and then integrating in both the positive and negative directions to obtain θ i (x) for |x| < L, for sufficiently large L. Since θ i always attains constant asymptotic values, these solutions correspond to finite-energy solitons.
As an example we have plotted θ i (x) for N = 15 in figure 6, showing explicit kink solitons for the choice of initial values θ 0 i as shown on the vertical axis at x = 0.The asymptotic values as x → ∞ have equal spacings of π N and similarly for x → −∞, provided that we take the values modulo 2π.Configurations such as those in figure 6 depend strongly on the selected values θ 0 i at x = 0, i.e. a small change in initial values can lead to completely different configurations, but always with equally spaced asymptotic values (modulo 2π).Unlike the solitons for the Kuramoto model as shown in figure 2, soliton trajectories always cross before attaining their asymptotic values, a property that holds for general N.

Exact solutions for the antisymmetric Kuramoto model for N = 3
We investigate the N = 3 case of (82) in order to determine exact solutions and their properties, in particular we find solitons which appear only for special initial values and so do not appear in general numerical computations.Define u = θ 2 − θ 1 , v = θ 3 − θ 2 as before, then the superpotential (81) with a ij = ε ij and K = 1 2 is given by from which we determine V using (22): which differs from (52) only in the sign of various terms.We have symmetries such as From (82) we obtain the Bogomolny equations: and the fixed points, which satisfy cos u = cos v, are given by together with 2π translations in either u or v.A contour plot of V in figure 7 shows the zeroes at ± π 3 , π 3 and their 2π translations (in red), as well as the fixed points at ±(π, π) (in blue).
For any solution (u, v) of (88) further solutions are given by (v, u) and (u, −u − v + π), and so we can satisfy (88) by setting either v = u or v = −2u + π, provided that This can be integrated directly, see Gradshteyn and Ryzhik [61], section 2.553, to obtain: There are three possibilities, depending on the value of and so u interpolates between the zeroes of V at (u, v) = ± π 3 , π 3 , indicated by the red dots in figure 7 lying on the line u = v either side of the origin.This corresponds to the solitons previously found numerically, such as those in figure 6, for which the spacing between asymptotic values is π N = π 3 , corresponding to the asymptotic value of u = v.
The other two possibilities arise when √ 3 tan u0 2 > 1, for which the solution satisfies either By appropriate choice of sign in the expression (91) we can deduce the asymptotic values of u, and find that u interpolates between either π 3 and π, or between − π 3 and −π.With reference to figure 7, these solitons join the zeroes of V, marked in red or blue, at either (u, v) = −(π, π) and − π 3 , π 3 , or at (u, v) = (π, π) and π 3 , π 3 .In addition to the solitons on the line v = u in the (u, v) plane, there are also solitons with v = −2u + π, for which (90) is again satisfied, which join points such as (π, −π) and π 3 , π 3 .By symmetry, there are also solitons which lie on the line u + 2v = π, where v satisfies (90), with similar properties.Such solitons do not appear for general N for random numerical values θ 0 i , since restrictions such as u = v are not satisfied unless we specifically choose u 0 = v 0 .
In this section we have briefly described solitons constructed using the odd superpotential for the sine-Gordon model, in order to indicate how their properties depend on those of the coupling coefficients a ij , in particular how the antisymmetry of a ij leads to equally-spaced asymptotic values for θ i , for both x > 0 and x < 0.

Soliton systems with polynomial φ 4 interactions
The well-known solitons of the φ 4 model in 1 + 1 dimensions have been studied from multiple perspectives and display many remarkable properties as discussed and summarized in the contributions in [62].We describe here a generalization of the φ 4 potential to N scalar fields, choosing as a starting point the superpotential w in (7) and hence define W = W(φ 1 , . . ., φ N ) as in (3): where the coupling coefficients a ij = −a ji are necessarily antisymmetric because w is odd.The quartic potential V is defined by (2): where we have set K = −1/6, and the Bogomolny equations read: Because a ij is antisymmetric, configurations for which φ i = φ j for all i, j are not fixed points, unlike the Kuramoto system (47) for which complete synchronization is possible.Nevertheless, finite-energy solitons exist with asymptotic values that are fixed points of (94).Since these fixed points are roots of multinomials, they have no simple closed-form expressions such as those for the antisymmetric model discussed in section 7.In common with this antisymmetric model, we have the invariance φ i → −φ i for all i, together with x → −x, and so the asymptotic values of the fields for x → ∞ are the negative of those for x → −∞, i.e. the system evolves for x < 0 with essentially the same properties as those for x > 0, unlike the standard Kuramoto model.
As before, while there are many ways to choose global couplings such that |a ij | = 1, a convenient choice is a ij = ε ij as defined in (83).We again generate numerical solutions by choosing random values φ 0 i ∈ (−1, 1) at x = 0, then integrating in both the forward and backward directions in x to obtain φ i (x) for all |x| < L. Unlike the periodic systems previously discussed, the initial values for the differences φ i − φ j must lie within a restricted interval, such as (−1, 1) or a smaller subset, otherwise solutions do not exist.The φ 4 kink soliton shown in (8), for example, satisfies |φ| < 1 and so, in order to generate this solution numerically, the initial value φ 0 , which corresponds to a choice for x 0 , must lie within (−1, 1).If a solution does not exist for any specified set of values φ 0 i , we simply repeat the process, perhaps with a smaller interval containing the values φ 0 i , until we obtain a valid numerical solution.We find that in all cases for which solutions exist for all |x| < L, they correspond to finite-energy solitons because they attain constant asymptotic values that are zeroes of V.
As an example, in figure 8 we plot soliton trajectories for N = 10 with initial values as shown on the vertical axis at x = 0.The asymptotic values for x > 0 correspond to those for x < 0, except that they appear in reverse order, but are not quite equally spaced.Unlike the Kuramoto model for which complete synchronization occurs, corresponding to an attractive force, for the φ 4 model the solitons repel each other and so asymptotically are spaced apart.If N is odd, one of the trajectories has identical values as x → ±∞, and so describes a solitonantisoliton pair.Different initial values φ 0 i lead to similar trajectories, but with the same value for |W ∞ − W −∞ |, for fixed N, where W ∞ = −W −∞ .Solitons can be generated numerically without difficulty for much larger values of N, but it remains to investigate their properties for different choices of coupling coefficients a ij .

Properties of the φ 4 model for N = 3
The φ 4 model has been extended many times to two or more independent fields with quartic potentials, see for example [7,9,12,14,16,19], however these models differ from those investigated here.For N = 3 the system (94) has two independent scalar fields, but for a direct comparison with previous models the derivative terms in the energy density need to be diagonalized as discussed in section 3.2.
As before, define u = φ 2 − φ 1 , v = φ 3 − φ 2 , then the superpotential (92) reduces to and depends on three parameters a 12 , a 23 , a 13 , which can be reduced to two by rescaling x in the subsequent equations.The quartic potential V (93) has the explicit form: where and satisfies V(u, v) = V(−u, −v).The Bogomolny equation (94) read: and have various symmetries, such as invariance under u ←→ v provided that a 12 ←→ a 23 .
As the fixed points of (98) have rather complicated expressions let us consider for simplicity the case a 23 = a 12 , then there are four fixed points (u, v) = ±α(1, −1) and (u, v) = ±β (1,1), where  which interpolates between the zeroes ±β(1, 1) of V, where β is defined in (99).Figure 9 shows a contour plot of V for a 12 = a 23 = a 13 = 1 for which showing four zeroes at ± 2 5 (1, 1) in red, and at ± √ 2 (1, −1) in green.The soliton (101) interpolates between the red points, however a soliton which joins the other two fixed points, in green, does not appear to exist.
As previously noted for general N, in solving (98) the initial values u 0 , v 0 must be suitably restricted in order that solutions exist, for example (101) exists only for initial values u 0 = v 0 such that |u 0 | < β.The coefficients a ij in W in (92) must also be selected such that V has a vacuum structure which supports solitons.For example, for a 12 = −a 13 the kink soliton (101) does not exist, since then α = β = 0.If we set a 12 = a 23 = −a 13 = 1 then V = (u 2 + v 2 + uv) 2  which is zero only at u = v = 0, and so solitons do not exist.
The generalization to N scalar fields is of interest because although the potential V has a more complicated vacuum structure than the φ 4 models, it is simpler in some respects due partly to the fact that the superpotential w in ( 9) is an even function of φ, whereas for the φ 4 model w is odd.Hence, the coupling coefficients a ij are symmetric, like those of the Kuramoto model but unlike the φ 4 systems.One consequence is that complete synchronization can occur, meaning that |φ j − φ i | → 0 as x → ∞ for all pairs i, j.Unlike the Kuramoto system, however, for which the nodes are uncorrelated as x → −∞, we find that for the φ 6 system the differences |φ j − φ i | approach either 0 or 1, which is a form of bipolar synchronization.In all cases we find numerically that the fields take constant values as |x| → ∞, and so can be identified as finite-energy solitons in the usual way.

Fixed points for the φ 6 model
We begin with the superpotential w in ( 9) and so from (3) define W according to where a ij is symmetric, leading to the Bogomolny equations for i = 1, . . .N, where we have set K = −1/8.Evidently this system is invariant under φ i ↔ −φ i .The potential is a multinomial in φ i of degree 6, specifically: Besides the fixed point corresponding to complete synchronization, for which φ i = φ 0 for constant φ 0 , there are further fixed points in which φ i = φ 0 for some nodes, and φ i = 1 + φ 0 for the remaining nodes.For any pair i, j, therefore, there are fixed points for which either φ j − φ i = 0 or |φ j − φ i | = 1.To be precise, denote N N = {1, 2, . . .N} and define the disjoint sets of indices S 0 , S 1 by for constant φ 0 , which can be transformed to zero by translation invariance.If both i, j ∈ S 0 , or if both i, j ∈ S 1 , we have φ j − φ i = 0, but if i ∈ S 0 and j ∈ S 1 then In all cases we have (φ j − φ i ) 1 − (φ j − φ i ) 2 = 0. Hence the system (104) has many fixed points, depending on how the nodes are distributed between S 0 , S 1 .
Numerically, for x → ∞ we always obtain φ j − φ i → 0 for all i, j, but for x → −∞ we have either φ j − φ i → 0 or |φ j − φ i | → 1, for which we observe that the nodes are almost equally distributed in number between S 0 , S 1 , depending on the initial values φ 0 i .If N is even, for example, there are configurations with N/2 nodes in each of S 0 , S 1 , but there are also others with an unequal distribution.These numerical observations apply when N ⩾ 4 for the couplings a ij = 1 for all i, j, and for randomly generated initial values φ 0 i .For N = 3 there exist further nontrivial fixed points and corresponding solitons, which we discuss in the next section.In figure 10 we plot the trajectories for N = 15 with a ij = 1, where the initial values φ 0 i are shown on the vertical axis at x = 0.As for the φ 4 model, the initial values for the differences φ i − φ j must lie within a restricted interval, such as (−1, 1), for solutions to exist.Complete synchronization as x → ∞ and bipolar synchronization as x → −∞ are evident in figure 10.We have W ∞ = 0 for any N, due to complete synchronization as x → ∞, as follows from (103), however the value of W −∞ depends on the distribution of nodes between S 0 , S 1 .One noticeable property for a ij = 1 is that trajectories never cross, as can be seen in figure 10.This also occurs for the Kuramoto system, where it follows from properties of the conserved cross ratios (50), but for the φ 6 system it can be proved directly from (104), see appendix, although there is no evidence that any conserved quantities exist which would prevent such crossings.
There is no numerical difficulty in calculating solitons for randomly-generated coefficients a ij for any N, provided that the initial values φ 0 i are such that solutions exist for all |x| < L. If a ij > 0 for most pairs i, j, then complete synchronization still occurs as x → ∞, however in such cases there exist many more nontrivial fixed points which the trajectories approach as x → −∞.9.2.The φ 6 model for N = 3 The system (104) has special properties for N = 3 which merit separate investigation.There are again two independent fields and so with u where a 12 , a 13 , a 23 are arbitrary parameters.From (105) V is given by:  V has zeroes at (0, 0), ±(1, 0), ±(0, 1), ±(1, −1), at each of which P = Q = R = 0, as well as nontrivial zeroes at two other points ±(α, β).
For simplicity let us choose a 23 = a 12 , then these two nontrivial zeroes lie at ±α(1, 1) where α = a 12 + 2a 13 a 12 + 8a 13 , (110) which exists only when the argument of the square root is non-negative.The equations for a 23 = a 12 are: and have nine distinct fixed points as given above.
Figure 11 shows a contour plot of V for a 13 = 2, a 12 = 1 = a 23 , with six zeroes at ±(1, 0), ±(0, 1) and ±(1, −1) marked in white, the zero at the origin in green, and the remaining two zeroes at ±α(1, 1) in red.We expect to find solitons which interpolate between the origin and any one of the eight surrounding zeroes.We show that such solitons exist by finding exact solutions for special initial values (u 0 , v 0 ).
One such solution is obtained by choosing v 0 = u 0 then v = u and both equations (111) and (112) are satisfied provided that which has the solution u (x) = ± a 12 + 2a 13 a 12 + 8a 13 + e 2(a12+2a13)(x−x0) , ( where we can choose either sign, due to the invariance of (113) under u → −u.This describes solitons which are similar in form to the φ 6 soliton (10), and interpolate between (u, v) = (0, 0) and one of the zeroes ±α(1, 1) of V, as shown by the green and red points in figure 11, where α is defined in (110).
A second exact solution is obtained by setting v = −u, then again both equations ( 111) and ( 112) are satisfied provided that u ′ = 3a 12 u (u 2 − 1), which has the solution The parameter a 13 does not appear in this expression, as is evident from the formula (107) for W, since u + v = 0.This soliton interpolates between one of the zeroes at ±(1, −1) (in white) and the origin (in green), depending on the sign in (115).Further exact solutions are obtained by setting v = 0, then both equations ( 111) and ( 112) are satisfied provided that a 13 = a 12 , and u ′ = 3a 12 u (u 2 − 1) which again has the solution (115).Similarly we can set u = 0, and so for these two cases we obtain solitons which interpolate between the zeroes of V at one of the points ±(1, 0), ±(0, 1) (in white) and the origin.These are essentially the same solitons as for the single-scalar φ 6 model.
These exact solutions can also be obtained numerically, but only for initial values which satisfy the stated restrictions.For randomly generated values (u 0 , v 0 ) we obtain configurations as described for general N, i.e. u, v → 0 as x → ∞, and for x → −∞ we have one of |u|, |v| → 1 and the other → 0. The exponents in these solutions, namely a 12 + 2a 13 in (114), and 3a 12 in (115), correspond to the particle masses µ 1 , µ 2 derived in section 3.5 for N = 3.The soliton (114) does not appear for randomly generated values (u 0 , v 0 ), unless we specify v 0 = u 0 .
We do not consider identical couplings with a 12 = 1 = a 23 = a 13 , since then V factorizes according to for which a continuum of zeroes lies on the ellipse u 2 + v 2 + uv = 1.

Bound states for the φ 6 model for N = 3
We apply here the properties of the linear stability equations, derived in section 3.3, to the φ 6 system with respect to the exact solution (114).The perturbative excitations ψ i in (27), when regarded as quantum particles, interact with the soliton (114) with the possibility of forming bound states, but more generally are reflected or transmitted through the soliton as determined by the matrix Schrödinger equation Hψ = Ω T Ωψ, where H is defined in (30).It is known for the single-scalar φ 6 theory that bound meson-soliton states do not occur [23], however for the N = 3 system we find that several such bound states can appear with respect to the soliton (114), as we now describe.
Without loss of generality we set a 12 = 1 in the solution (114), which is equivalent to rescaling the x variable in (113) by a 12 > 0, and define the mass parameter µ = 2a 13 + 1, with µ > 1.By suitable choice of x 0 in (114) we can write the soliton as then the underlying fields φ i are given by where c is an arbitrary constant.We observe that u, as defined here, differs from the singlescalar φ 6 soliton (10), because it interpolates between 0 (for x → ∞) and a nontrivial zero of V at µ/(4µ − 3) (for x → −∞), which appears only for N = 3. Next, we calculate the 3 × 3 symmetric matrix M = (m ij ) as an explicit function of x using the formula (38) with φ i as given in (118).Surprisingly, M can be diagonalized with a constant orthogonal matrix.Define then PP T = P T P = I 3 , and we find that PM(x)P T = D(x) is diagonal with D = diag(0, d 1 , d 2 ) where Hence the 3 × 3 potential matrix U, defined in (29) and satisfying U = M 2 + M ′ , can also be diagonalized according to PUP T = diag(0, U 1 , U 2 ), where We ignore the zero diagonal element of U, which arises because the fields φ i in (118) sum to a constant value.
We wish to evaluate therefore the spectra of the two operators H 1 , H 2 given by where we write each potential U 1 , U 2 in the form sometimes referred to as a Rosen-Morse potential, for constants c 0 , c 1 , c 2 to be identified.The corresponding Schrödinger equation can be solved exactly in terms of hypergeometric functions using the formulas in Morse and Feshbach [67], page 1651.Alternatively, one can use the factorization method together with further algebraic properties, as derived in [36], to analyze H 1 , H 2 , however we use the formulas in [67] directly.Let us firstly dispose of the second case with H 2 , then the Schrödinger equation reads −ψ ′ ′ + U 2 (x)ψ = ω 2 ψ where U 2 is given by ( 122), where we identify c 0 = 5µ 2 2 , c 1 = − 3µ 2 2 , c 2 = − 15µ 2 4 .There exists only one bound state for which ω = 0, however this is merely the zero mode corresponding to the wavefunction ψ(x) = u ′ (x), where u is defined in (117), which satisfies −ψ ′ ′ + U 2 ψ = 0 together with −ψ ′ + d 2 ψ = 0, where d 2 is given in (120).There are no other bound states, and so this case is similar to that of the single-scalar φ 6 theory [23].
We do not further develop here the properties of the meson-soliton interactions determined by the matrix Schrödinger equation, as our aim is to show that the N = 3 system allows much more intricate interactions than does the single-scalar system, as well as providing a tractable example of the properties of the matrix operator H in (30).Similar properties hold for the Kuramoto system discussed in section 5.2.

Conclusion
We have shown that models of synchronization in complex systems, in which ensembles of phase oscillators display coherent behaviour and other emergent phenomena as they evolve from generic initial values, are related to scalar field theories with N fields φ 1 , . . .φ N through corresponding superpotentials or Lyapunov functions W, sometimes referred to as disagreement functions in the context of consensus seeking multi-agent systems [47].The first-order Bogomolny equations for static fields have a gradient form with respect to W, leading to solitons which are defined on the real line and have finite energy.The same equations, regarded as time evolution equations, can be viewed as models of synchronization phenomena in complex systems of oscillators, in which solutions asymptotically approach fixed points, which correspond to absolute minima of the potential V. The established properties of phase oscillators, such as their asymptotic values, lead to corresponding properties for solitons such as finite energy.
The relation between synchronization phenomena and static solitons does not apply in all circumstances, however, since there are many models of synchronization which do not have a gradient structure and so do not correspond to any Bogomolny equations.In any case, synchronization is relevant to solitons only for identical natural frequencies of oscillation, for which the phase oscillators asymptotically approach constant values, as is required for finite energy solitons.Generally, the finite energy requirement is weaker than the property of synchronization, since some of the models that we analyze do not show any discernable synchronization behavior for general coupling coefficients, yet still allow acceptable soliton solutions.
We have investigated in detail the Bogomolny equations for the sine-Gordon system, generalized to N scalar fields, since they are identical in form to the well-known equations of the Kuramoto model.We have also taken as examples the soliton systems that follow from the well-studied sine-Gordon, φ 4 and φ 6 models, even though the polynomial models no longer correspond to phase oscillators, and have highlighted their differences as determined by the vacuum structure of the potential V. The polynomial models also suggest new dynamical models of complex systems which allow a much broader range of synchronization phenomena beyond the well-known phase or frequency synchronization properties.We have also discussed the elementary excitations of the classical fields and their soliton interactions, which allows us to investigate bound soliton-meson states.
It is clear that multisoliton configurations can be constructed starting with any suitable superpotential w, for example w(φ) = Kφ 3 (5 − 3φ 2 ) leads to a φ 8 system where W is defined as in (2) and V is constructed according to (3), then solitons again exist for any number of scalar fields φ i .Many other generalizations can also be considered, for example the definition of W in terms of w in (3), which allows pairwise interactions only, can be generalized to include higher-order interactions.

Figure 2 .
Figure 2.Soliton trajectories θ i (x) for N = 15, showing complete synchonization for large x in which θ i (x) → δ for some constant angle δ, for which the order parameter r → 1, whereas for negative x the asymptotic trajectories settle into one of the many zeroes of V such that r → 0.

Figure 3 .
Figure 3.The soliton energy densities (θ ′ i ) 2 for each of the N = 15 trajectories in figure 1, showing distributed soliton locations and energies.

Figure 4 .
Figure 4. Contour plot of V for N = 3 as given in (52), showing the zeroes surrounding the origin (in green) at the points (m, n) π (in blue), and 2π 3 (m + 2n, m − n) (in red) for integers m, n.
we plot ζ(x) in the complex ζ-plane for |x| < 2 for the same solitons as in figure 2 for N = 15, showing the asymptotic endpoints ζ −∞ and ζ ∞ .Having calculated ζ, we can construct the soliton trajectories shown in figure 2 by means of (66) and (70), modulo 2π.

Figure 5 .
Figure 5.The trajectory ζ(x) (in red) in the complex ζ-plane for |x| < 2, passing through the initial value ζ(0) = 0, with the asymptotic limit |ζ| → 1 for large x, for the case N = 15 with initial values as for the solitons in figure 2.

Figure 8 .
Figure 8. Soliton trajectories φ i (x) for N = 10 for the φ 4 equation (94), showing asymptotic values which are almost equally spaced, but with the order reversed from left to right.

Figure 10 .
Figure 10.Soliton trajectories φ i (x) for N = 15 for the φ 6 system (104) with a ij = 1, showing complete synchronization as x → ∞, and bipolar synchronization as x → −∞, with a difference of unity between the two asymptotic branches.The trajectories never intersect.
a 12 a 23 PQ + a 12 a 13 PR + a 13 a 23 QR, ( Evidently these exist only if the arguments of the square root functions are non-negative.Although (98) does not in general have exact solutions, we can satisfy both equations for a 23 = a 12 by setting v = u, provided that u ′ = − (a 12 + a 13 ) + (a 12 + 4a 13 ) u 2 .