Non-holonomic constraints and their impact on discretizations of Klein-Gordon lattice dynamical models

We explore a new type of discretizations of lattice dynamical models of the Klein-Gordon type relevant to the existence and long-term mobility of nonlinear waves. The discretization is based on non-holonomic constraints and is shown to retrieve the"proper"continuum limit of the model. Such discretizations are useful in exactly preserving a discrete analogue of the momentum. It is also shown that for generic initial data, the momentum and energy conservation laws cannot be achieved concurrently. Finally, direct numerical simulations illustrate that our models yield considerably higher mobility of strongly nonlinear solutions than the well-known"standard"discretizations, even in the case of highly discrete systems when the coupling between the adjacent nodes is weak. Thus, our approach is better suited for cases where an accurate description of mobility for nonlinear traveling waves is important.


1.
Introduction. Discrete solitons such as kinks and breathers exist in many physical systems, including Josephson juctions, electrical circuits, and granular crystals [1,2,3]. In the models that describe these systems, the mobility of soliton solutions, and specifically kinks -which will be the focus of the present work-, is an important and well-studied property [4,5]. For instance, kinks in extemely discrete systems need to overcome the celebrated Peierls-Nabarro barrier, while they undergo energy radiation which leads to their trapping in the lattice [9]. More recent work revealed the existence of kinks that do not experience the action of the Peierls-Nabarro potential [10].
For a class of Klein-Gordon discretizations [5] it has been shown that energy and linear momentum cannot be both conserved. This is problematic in cases when the objective is the discretization of the continuum version of equation, which conserves both momentum and energy. Going beyond the properties of isolated kinks, the above result has consequences for the studies of soliton collisions in Klein-Gordon lattices. In general, the mobility of kinks in discrete manifestations of the Klein-Gordon equation may vary from that in the continuous version as well as among different discretizations.
In this work, we present a class of discretizations based on a non-holonomic constraint that preserve the momentum and alleviate the problem of decreased mobility in kink solutions. Specifically, as we enforce the conservation of momentum it in turn implies an affine constraint in kink velocities, which is of the non-holonomic type. The extra term added in the discretization due to the non-holonomic constraint, vanishes in the continuum limit, hence demonstrating that the discretized version of the equation converges to its continuous counterpart. Our numerical results validate our theoretical considerations and we present multiple instances where the additional constraint that we impose increases the mobility of the kinks.
We shall note that for generalized non-holonomic constraints that are not linear in velocities, one has to be careful in applying the variational principle. This was demonstrated in, for example, [6,7] for the case of stabilizing servomechanisms, where an alternative formulation through Poisson structures has been derived. The discrete equivalent of a conservation law for continuous system (1) below, enforced by the constraint, cannot be readily interpreted as a consequence of a continuous symmetry through a Noether theorem. Thus, applications of alternative methods of non-holonomic mechanics [8], in addition to the Lagrange-d'Alembert's principle developed here, are of particular interest and will be considered in future work.
This article is structured as follows. In section 2 we describe the Klein-Gordon model that will be studied and present the theoretical framework of our method. In section 3 we present numerical results that corroborate our theoretical predictions. Finally, in section 4 we discuss our findings and offer insight for future directions.
2. Model Theoretical Setup and Proposed Discretizations. The models that we will consider herein will involve discretizations of continuum partial differential equations of the nonlinear wave type: There are numerous key model examples of this class of equations, such as the sine-Gordon [1] or the φ 4 model arising from classical field theory and high-energy physics [2], as well as variants thereof. The corresponding discrete models are of considerable interest in their own right as models of coupled torsion pendula, dislocation dynamics, Josephson junctions, electrical circuits, among numerous other settings [1,2,11,12]. The generic discrete analogue of these models is derived by using a standard lattice discetization of the second space derivative as follows: The field is defined at the integer nodes n and ǫ denotes the coupling constant between adjacent sites, which could be interpreted as ǫ = 1/dx 2 , with dx being the spacing between the adjacent sites (see below). Note that lattice model (2) is Hamiltonian with where H n is the energy density.

NON-HOLONOMIC CONSTRAINTS IN KLEIN-GORDON LATTICES 3
In that context, an alternative limit, namely ǫ = 0, the so-called anti-continuum limit becomes particularly important and provides a valuable tool for assessing the potential stationary or time-periodic states of the model and their respective spectral stability properties [13,14,15].
One of the "pathologies" of such discrete models is that as the coupling parameter ǫ decreases, the mobility of the nonlinear waves accordingly weakens [4]. The decreased mobility stems, to a considerable degree, from the fact that the momentum P = u t u x dx is no longer a conservation law in the discrete case of Eq. (2), even though it is in the corresponding continuum limit of ǫ → ∞. This curious fact has spurred a considerable volume of recent work (for a partial summary, see [16]) focused on providing alternative discretizations that would enable higher mobility of the waves, even at small ǫ. Our aim in the present work is to propose a novel class of such discretizations, based upon the imposition of a non-holonomic constraint providing an exact conservation of the discrete momentum.
Indeed, the continuum equation (1) provides a conservation of the momentum quantity P c = u t u x dx. On the other hand, after defining the analogue of discrete momentum for (2) as a quick check shows that such discrete momentum is not conserved [17]. The goal of this paper is to enforce the conservation of (4) using the methods of non-holonomic mechanics. The equation (2) can be obtained as the Euler-Lagrange equation with the Lagrangian through the critical action principle δ Ldt = 0 on arbitrary variations δu n (t) with appropriately vanishing boundary conditions at the end time points. Next, enforcing the conservation law (4) leads to an affine constraint in velocities of the type n a n (u)u n = b(u) with a n (u) = u n+1 − u n−1 , b(u) = P (6) using the notation u = (u 1 , u 2 , . . . , u n ). Since this constraint involves velocities in a non-trivial way, i.e., that constraint cannot be integrated to contain coordinates only, the constraint is non-holonomic, and linear in velocities. The usual way of dealing with this type of constraints is through the Lagrange-d'Alembert's principle [18]. For the constraints of the type (6), the variations δu n are no longer arbitrary, but must satisfy the linear relations n a n (u)δu n = n (u n+1 − u n−1 ) δu n = 0 .
Notice that we do not enforce (6) directly, as this would correspond to the socalled vakonomic approach, relevant for e.g. the control theory but normally not leading to correct expressions for mechanical systems, as described in [19]. On the contrary, Lagrange-d'Alembert's approach (7) was consistently shown to be physically appropriate. For more in-depth discussion of various approaches to nonholonomic systems, we refer the reader to [18,19].
If r is the Largange multiplier enforcing the constraint (7), the critical action principle coupled with the constraint gives provided δu n vanish at t = t 0,1 so no boundary integration terms appear in the equations. Since δu n is arbitrary, the equations of motion conserving momentum areü Note the new term proportional to r and enforcing the constraint on the right hand side of (8), as compared to (2). The physical meaning of r is the (generalized) constraint force. In order to find r, we compute the derivative dP/dt = 0 and set it equal to zero. We find that while the terms associated with ǫ vanish telescopically [17], the nonlinear term does not. This leads to the explicit condition for r as Making this selection for r, multiplying both sides of Eq. (8) byu n , and summing, we observe that the derivative of the (former) Hamiltonian of Eq. (3) is not zero anymore, but is given by the following condition: Importantly, this identity suggests that for generic initial data, such that P = 0, the Hamiltonian will not be conserved. This is due to the constraint (6) being non-homogeneous in velocities for P = 0. It is relevant to point out that such an "incompatibility" of the momentum and energy conservation laws for other broad classes of discretizations of the Klein-Gordon dynamical lattices have also been reported previously [5].
An additional comment about this non-holonomic discretization is that the extra term introduced by the non-holonomic constraint vanishes at the continuum limit ǫ → ∞, ensuring that Eq. (8) converges to the proper continuum limit of Eq. (1). To see this, we consider the continuum analogue of the new r(u n+1 −u n−1 ) term, which is u x ( V ′ (u)u x dx/ u 2 x dx). The numerator of the term in the parenthesis can be integrated to yield V (u 2 ) − V (u 1 ) where u 2 and u 1 are, respectively the asymptotic values of the field at ±∞. Assuming that we examine the dynamics of homoclinic orbit or of a heteroclinic orbit in a symmetric double well (as is typically the case in the examples of interest mentioned above), this quantity vanishes identically and the discretization reverts to its proper continuum limit. We will return to this point and corroborate it through our numerical computations below. A final remark about the nature of our discretization procedure, different from all other schemes that we are aware of, is that our method is inherently nonlocal since the computation of r in (9) entails a summation over the entire lattice.
3. Numerical Computations. To provide a comparison between the properties of the different discretizations, we examine the dynamics of Eq. (2) vs. that of Eq. (8). We use an initial condition in the form of a prototypical nonlinear excitation such as the kink. In the continuum limit, the latter stationary solution (which can be boosted to any sub-luminal speed via Lorentz transformation) is of the form u(x) = tanh(x). This is the case for the φ 4 model with a potential V (u) = (1/2)(u 2 − 1) 2 . For the discrete system, we have used the initial condition u n (0) = tanh(x) = tanh(n × dx), where dx is used to signal the spatial discretization step (see also below), approximating the continuous solution on a discrete set of points.  In addition, we enforce an initial momentum of P f ix = 0.5 via the definition of Eq. (4). The simplest way that we envisioned for achieving this was to assume all sites as having vanishing velocity except for the central one of the chain, say n 0 , and then usingu n0 = P f ix /(u n0+1 − u n0−1 ). This prescription provides the given amount of momentum to the central portion of the kink and we expect that it should set the kink in motion.
This is exactly what we observe in Fig. 1, which is close to the continuum limit of ǫ = 1/dx 2 → ∞. Here, as usual dx is effectively the corresponding spacing between the lattice sites, leading as a finite-difference three-point approximation to the second derivative. In particular, in this figure dx = 1/4 or ǫ = 16 is used. The top left panel illustrates the traveling kink of the standard discretization of Eq. (2). The proximity of the model to its continuum analogue coupled with the imparting of momentum via the initial condition leads the kink to travel through the chain with an apparently constant speed (evident in the constant slope of the two contour plots of the top panel, the left one illustrating the energy density of Eq. (3) and the right one the actual field u n (t)). Given that the model of Eq. (2) is Hamiltonian, our explicit fourth-order time integration (Runge-Kutta) scheme conserves the initial energy (of O(1)) up to a few parts in 10 −14 verifying its conservation for all practical purposes. On the other hand, the momentum fluctuates in the third decimal digit indicating its non-preservation by the scheme; instead, it is only approximately conserved because of the proximity to the continuum limit where the continuum analogue of Eq. (4) can be shown to be preserved by Eq. (1).
On the other hand, we can see that the fluctuations in this quantity (P ) are far more dramatic in the considerably more discrete case ǫ = 1, corresponding to  Figure 2. Same as the previous plot, only now we consider a highly discrete case dx = ǫ = 1. Here, we observe the detrimental effect of the energy-conserving discretization with regards to the mobility of the kink initial condition. The strong discreteness/weak coupling precludes the kink from moving through the lattice and the initial momentum (which as can be seen is far from being conserved now) is expended in a jerky motion with a range of a few discrete sites (trapping of the kink).
dx = 1, shown on Fig. 2. There, the fluctuations in P are of O(1) while the energy is still conserved with the same precision as before. In addition, the motion of the kink is fundamentally different, featuring a jerky back and forth trajectory between a few sites around the origin. Thus, in this case, the kink appears trapped due to the weak inter-site coupling which does not allow the momentum to translate itself into mobility as it would in the corresponding continuum limit. We now turn to the corresponding dynamics of the newly proposed model of the form of Eq. (8). The case of dx = 1/4 is again shown in Fig. 3. The fundamental difference in this case is that the momentum P of Eq. (4) in the bottom left panel is conserved to a similar accuracy as the energy in the previous model, as expected by our construction. On the contrary, the energy of the bottom right panel varies at the third significant digit remaining roughly constant solely due to the proximity to the continuum limit, where the analogue of Eq. (3) is conserved for Eq. (1).
The main advantage of our discretization method is evident in Fig. 4. More specifically, as the example of dx = ǫ = 1 shown in Fig. 3 illustrates, the discretization of Eq. (8) presented here can accurately describe traveling waves in the case of strong coupling. Fig. 3 shows that the momentum preserved up to numerical accuracy, as expected, while the energy fluctuations become more significant, revealing the absence of energy conservation already in the second significant digit. It should also be mentioned that we have also separately verified Eq. (10) numerically, although we do not present it here. Hence, our principal numerical conclusion is that the models of Eq. (8) ensure a higher mobility of their corresponding coherent structures compared with the standard discretization of Eq. (2). This comes at the price of the incorporation of the nonlocal effect through the summation that ensures the conservation of the linear momentum.
Lastly, we have confirmed that the relevant discretizations are indeed properly converging to the continuum limit as dx → 0. A typical example of this approach is shown in Fig. 5. In particular, the figure showcases the evolution of the quantity r(t), which constitutes the prefactor of the newly proposed nonlocal term arising through the non-holonomic constraint. It can be seen that as we go from the weaklycoupled strongly discrete case of the left panel (for dx = 1) to the right panel, closer to continuum case of dx = 0.25, the value of r(t) decreases by more than an order of magnitude 4. Conclusions and Future Challenges. In the present work, we have revisited the theme of discretizations of continuum nonlinear partial differential equations of the wave type. We have focused on the example of Klein-Gordon equations and offered a novel class of schemes which can be summarized as follows: (a) they can be straightforwardly designed to conserve linear momentum; (b) when they do so, they do not in general conserve the energy, (c) they revert to the continuum PDE in the appropriate limit of the lattice spacing tending to zero; and (d) they are based on the imposition of a non-holonomic constraint and thus are fundamentally different than previously proposed schemes due to their nonlocal nature. We have argued that these models have intriguing characteristics, such as the ability to impart high mobility to localized coherent states, a feature corroborated via numerical computations in the case of kinks of the φ 4 model. This feature is preserved even in the regime of high discreteness/weak coupling where the standard discretizations  Figure 4. Same as Fig. 4 but now for ǫ = 1 (corresponding to dx = 1). Notice that despite the high discreteness of the model and the weakness of the coupling, the mobility of the state is well preserved. under the same initial data produce kinks that are trapped within the Peierls-Nabarro energy barrier and cannot move through the chain. This study opens an interesting new direction for the investigation of the relevance of such schemes from a numerical perspective and also for a more detailed understanding of their properties. While in the present initial communication we have focused on the dynamics of kinks, it would be relevant to explore corresponding structural features for discrete breather type solutions [20]. Furthermore, while in the present work we have focused on the one-dimensional domains only, it would be particularly interesting to consider generalizing such models in higher dimensions. The latter have progressively become a theme of considerable attention over the past few years [21,22] and hence it would again be particularly relevant to compare the standard higher dimensional discretizations, with the ones produced by the multi-dimensional generalization of the non-holonomic constraint concepts presented herein. Such studies are presently in progress and will be reported in future publications.