Quasi-integrable deformations of the $SU(3)$ Affine Toda Theory

We consider deformations of the $SU(3)$ Affine Toda theory (AT) and investigate the integrability properties of the deformed theories. We find that for some special deformations all conserved quantities change to being conserved only asymptotically, {\it i.e.} in the process of the scattering of two solitons these charges do vary in time, but they return, after the scattering, to the values they had prior to the scattering. This phenomenon, which we have called quasi-integrability, is related to special properties of the two-soliton solutions under space-time parity transformations. Some properties of the AT solitons are discussed, especially those involving interesting static multi-soliton solutions. We support our analytical studies with detailed numerical ones in which the time evolution has been simulated by the 4th order Runge-Kutta method. We find that for some perturbations the solitons repel and for the others they form a quasi-bound state. When we send solitons towards each other they can repel when they come close together with or without `flipping' the fields of the model. The solitons radiate very little and appear to be stable. These results support the ideas of quasi-integrability, {\it i.e.} that many effects of integrability also approximately hold for the deformed models.


Introduction
Solitons play a very important role in the study of non-linear phenomena because often they arise in the mathematical description of the behaviour of some physical systems. Many properties of solitons are associated with the integrability of the mathematical models in which they arise. In such cases solitons are described as localised classical field configurations of the model that propagate without dissipation and dispersion. Moreover, when two such solitons are scattered they do not destroy each other but come out of their interaction region essentially unscathed. The only lasting effect of the scattering is a shift in their positions relative to the values they would have had, had they not encountered each other. The usual explanation of this behaviour involves the integrability of the model and associated with it existence of an infinite number of conserved quantities. These conservation laws dramatically constrain the soliton dynamics. The integrable theories are, however, very special as they possess highly non-trivial hidden symmetries. So, even small perturbations of these theories can destroy these symmetries and it is important to check whether any of these properties still hold when the underlying mathematical models are nonintegrable. Afterall, one would expect some 'continuity' of the properties as one introduces small (or not so small) perturbations.
We have looked at this problem and recently we have found that some non-integrable field theories in (1 + 1) dimensions, present properties similar to those of exactly integrable theories [1][2][3][4][5]. They have soliton-like field configurations that behave in a scattering process in a way which is very similar to true solitons. We have also shown that such theories possess an infinite number of quantities which are not exactly time-independent but are, however, asymptotically conserved. By that we mean that the values of these quantities change during their scattering process, and at times change a lot, but after the scattering, they return, to the values they have had before it. This is an interesting property since from the point of view of the scattering what matters are the asymptotic states, and so a theory in which solitons behave like this looks a bit as an effectively integrable theory. For these reasons we have named this phenomenon quasi-integrability. The mechanisms responsible for this behaviour are not properly understood yet, but we believe that this behaviour will play an important role in the study of many non-linear phenomena. Since integrable theories are rare and, in general, do not describe realistic physical phenomena, the quasi-integrable theories may play a significant role in the description of more realistic physical processes.
Most of the models we studied so far [1][2][3][4][5][6][7][8] involved (1+1)-dimensional theories of either one real scalar field φ subjected to a potential which is a deformation of the Sine-Gordon potential or a complex field which satisfied a modified non-linear Schrödinger equation or equation of the modified Bullough-Dodd model. The original models were integrable and the deformation of their potentials made them non-integrable.
Here we decided to extend our investigations to systems with more fields and so we have had a look at the SU(N ) Toda models and their deformations. All such undeformed models are integrable and the lowest of them (N = 2) is, in fact, the Sine-Gordon model in disguise. So, in this paper we report results of our study of the next model in this family of models, namely, of the SU(3) one.

JHEP05(2016)065
The paper is organized as follows. In section 2, we present this model and discuss some of its properties and in particular its symmetries. We also suggest a possible deformation of the model which possesses most of these symmetries. The following section discusses various properties of both the undeformed and deformed models such as their quasi-zero curvature conditions and the resulting quasi-conserved quantities. Section 4 discusses how the fields of these models change when one Lorentz transforms them and when they lead to charge conservation. We also present the explicit expression for the anomaly termswhich control the situation when the charges are only asymptotically conserved (which corresponds to our ideas of quasi-integrability). In section 5 we discuss the well known soliton solutions of the undeformed model paying particular attention to the solutions which describe static solitons. In section 6 we discuss the interplay between parity and dynamics.
The following two sections describe the numerical procedure used by us for checking some of these claims and present the results of our numerical investigations. In fact all our results were obtained using the 4th order Runge-Kutta method to simulate the time dependence of field configurations. First we performed such numerical evolutions of field configurations for which we had analytical expressions. This not only checked our numerical schemes but also demonstrated that the soliton solutions of the un-deformed SU(3) model were really stable, with respect to small numerically induced, perturbations. Then we looked at the deformed models for various values of the deformation and for solitons at rest. We followed these studies by looking at solitons moving towards each other at various speeds. In section 8 we present some of our conclusions.
In this paper we consider deformations of the integrable Affine Toda model, such that we keep the kinetic term in (2.1) unchanged, but take the potential to be of the form where φ is still given by (2.2), and v j , j = 0, 1, 2, are vectors in the root space of the SU(3) Lie algebra, which are deformations of the roots α j . The choice of the vectors v j is restricted by some conditions which we will discuss below. The Hamiltonian density and energy associated to (2.1) are given respectively by Since the fields are complex, so are the Hamiltonian density and energy. Therefore, such models do not possess vacua solutions that minimize the energy. However, in order for the energy to be conserved in time, it is necessary to require that the flows of momenta at both ends of spatial infinity are equal, i.e. that For the solutions which we consider in this paper this condition is satisfied as space and time derivatives of the fields vanish at spatial infinity. For static configurations there is a further point to take into account. It is well known that for theories of the type we are considering the quantity is independent of x for static solutions of the equations of motion, i.e. d E d x = 0. This corresponds to the energy of a mechanical problem of a particle moving in φ-space in an inverted potential with x playing the role of time. Therefore, for static solutions for which the space derivatives of the fields vanish at spatial infinity one finds that the conservation of E in x, implies that However, for the static equations (2.3) to be satisfied at spatial infinity one requires that This imposes conditions on vectors v i . Let us restrict our interest to the cases where v 1 and v 2 are linearly independent and consider the dual basis w a , such that w a · v b = δ ab , a, b = 1, 2. Then, taking the scalar product of (2.13) with w a one finds that (2.14) Next we note that we have to discard the cases where v 0 is orthogonal either to w 1 or w 2 , since (2.14) would imply that the imaginary part of φ (±) had to diverge, and so the derivatives of the fields would not vanish asymptotically at spatial infinity as we have assumed. One then concludes from (2.14) that Using (2.15) one can conclude that (2.12) implies that Thus we have two possibilities. Either However, we are really interested in theories that can be deformed away from the Affine Toda model in a continuous manner. If one takes ( w 1 + w 2 ) · v 0 = 1 then there is no way of having v j , j = 0, 1, 2, as close as possible to α j . So we shall discard the possibility (2.17). The second case (2.18) implies that the difference of the asymptotic values of the fields has to live on a dual lattice, i.e.

JHEP05(2016)065
Let us restrict our attention to deformations that preserve, as much as possible, the symmetries of the Affine Toda model. For instance, the undeformed model (2.6) is invariant under the exchange φ 1 ↔ φ 2 . In addition, for the solutions which satisfy either the condition φ 2 = −φ * 1 , or φ a = −φ * a , a = 1, 2, the energy becomes real. So, in order to keep such symmetries and the reality conditions for the energy, we consider in this paper the following deformation: with ε being a real parameter. Note that this corresponds to taking n 1 = n 2 = −1 in (2.20) and so v 0 is expressed in terms of v i like α 0 in terms of α i . It then follows that v 1 · α 1 = v 2 · α 2 = 2, and v 1 · α 2 = v 2 · α 1 = − (1 + ε). In addition, one finds that With such a choice, the potential (2.7) becomes Note that the vectors v a , a = 1, 2, correspond to the deformations of the simple roots α a of SU(3) which modify the angle between them, and rescale their lengths equally, as shown in (2.22). The dual basis associated to the choice (2.21) is given by (2.25) As we have remarked above, a given solution satisfying the condition φ 2 = −φ * 1 , has real energy. Therefore, for such static solutions one needs m 1 = −m 2 , and so At the same time we observe that a solution satisfying the condition φ a = −φ * a , a = 1, 2, also has real energy, and a static solution of this kind can only exist when m 1 = m 2 = 0. To discuss integrability of the model we introduce the Lax potentials as with v 0 being a constant, and with H n αa , a = 1, 2, and E n ±αs , s = 1, 2, 3, n ∈ Z Z, being the Chevalley basis of the SU(3) loop algebra described in appendix A.
The curvature of such potentials takes the form Here ω is a cubic root of unity other than unity itself, i.e. ω 3 = 1 and ω = 1, and so 1 + ω + ω 2 = 0. In addition, we have where Note that, as W a , a = 1, 2, are functions of ω, in the calculation of X 2 one has to interchange ω ↔ ω 2 in the expressions for W a given above. The coefficients of H 0 αa , a = 1, 2, in (3.3) are exactly the equations of motion (2.4) of the deformed models we are considering, and so they vanish when evaluated on the solutions of such models. In order for the curvature F +− to vanish one needs the anomalies X a , a = 1, 2 to vanish, and so one has to choose potentials that satisfy the four equations, W a (ω) = -7 -JHEP05(2016)065 W a ω 2 = 0, for a = 1, 2. If one takes an ansatz of the form V ∼ [exp (i γ a φ a ) − v 0 ], then these four equations become four algebraic equations for the unknowns γ 1 and γ 2 . One can check that the only possible solutions are three choices: and so any linear combination of the form , leads to the vanishing of the anomalies, and so to an exactly integrable field theory. The Affine Toda model, corresponding to all q j = 0, j = 0, 1, 2, and the so-called Conformal Toda model corresponding to q 0 = 0, are examples of such integrable models.

The quasi-conserved quantities
In order to calculate the quasi-conserved quantities for the theories (2.1) we employ a modified version of the technique widely used in integrable field theories [9][10][11][12]. This procedure is called the abelianization procedure because it consists of gauge transforming the Lax potentials into an infinite abelian sub-algebra of the SU(3) loop algebra. In our case, due to the fact that the potentials (3.1) are not really flat, we are able to gauge transform only one component of (3.1) into the abelian sub-algebra. The main ingredient of the technique relies upon the fact that the generator b −1 introduced in (3.2), is a semisimple element of the SU(3) loop algebra G. By this we mean that the kernel and image of the adjoint action of b −1 have no intersection and G splits into the vector space sum of kernel and image, i.e.
The second important ingredient of the technique is an integer gradation of the SU(3) loop algebra G, such that The relevant gradation for our case is the so-called principal gradation performed by the grading operator where H 0 αa , a = 1, 2, are the generators of the Chevalley basis of the Cartan sub-algebra of G, and λ is the so-called spectral parameter of the loop algebra (see appendix A for details).
The calculations become simpler if one uses a special basis for G, described in appendix A, where the generators of the kernel are denoted as b 3n±1 , n ∈ Z Z, and the generators of the image as F a n , n ∈ Z Z, a = 1, 2, and they have well defined grades w.r.t. D, i.e.
[ D , b 3n±1 ] = (3n ± 1) b 3n±1 ; [ D , F a n ] = n F a n . (3.11) In terms of such a basis the Lax potentials (3.1) become (see appendix A for the definition of the new basis) where we have redefined the fields as Next we perform a gauge transformation with a group element which is an exponentiation of the positive grade elements of the image of the adjoint action of b −1 , i.e.
a F a n .
(3.14) We first consider the a − -component of the transformed Lax potential, and split it into the eigensubspaces of the grading operator (3.10) as − . We then get that One can now choose the parameters ζ . In addition, note that ζ (n) a is a polynomial in x − -derivatives of the fields ϕ a , and each term of such polynomials contains precisely n x − -derivatives. Note also that in such a recursive process of canceling the image components of A (n+1) − we do not use the equations of motion. Thus we find that the a − -component of the transformed Lax potential becomes can restrict our attention to fields which satisfy the equations of motion and use them, or equivalently the quasi-zero curvature condition to determine its form. The curvature F +− , given in (3.3), gets transformed into where in the last equality we have imposed the equations of motion (2.4) (see (3.3)). Since the group element g is an exponentiation of generators of strictly positive grades, it follows that g F a 1 g −1 has also strictly positive grades only, and so we can split it into its image and kernel components as From (3.12) we observe that A + has grade one components only, and so a + has strictly positive grades only. Thus the split of a + into its image and kernel components gives us: a (n,a) + F a n . (3.20) Next we put (3.17), (3.20) and (3.19) into (3.18), and find that the kernel component leads to Note that the r.h.s. of (3.22) does not have components of zero grade but the l.h.s. does. Therefore one concludes that a  Continuing such a process one observes that the zero curvature condition implies that the a + -component of the Lax potential is also transformed into the abelian kernel generated by b M . In addition, for integrable theories X a vanish and so one gets that the r.h.s. of (3.21) also vanishes for any M . For non-integrable field theories where the anomalies X a do not vanish, none of this happens. However, the fact that a + is not -10 -

JHEP05(2016)065
transformed into the kernel does not affect (3.21), and so we can get quasi-conservation laws as we explain next.
From (2.5) we see that the x and t components of the Lax potentials are a x = a + + a − and a t = a + − a − . So we introduce the charges By imposing appropriate boundary conditions at spatial infinity on the a t component of the Lax potential one gets from (3.21) that From (3.12) and (3.20) one observes that a , and so it turns out that a (1) x is a linear combination of the energy and momentum densities. This explains the origin of the conservation of the charge Q (1) , given in (3.24), even for the non-integrable case.
In our numerical simulations we have studied the behaviour of the charge Q (2) , and so the important quantities for evaluating the anomalies are then Choosing v 0 = −1 in (3.6), it follows that quantities X a , a = 1, 2, given in (3.5), evaluated for the potential (2.23), become Using (3.25) and (3.13) we then find that

JHEP05(2016)065
We have also investigated the quasi-conservation of the second charge which satisfies (see (3.24 Thus using (3.28) we find that the total anomaly is given by Λ : with λ being the rapidity, related to the velocity v by v = tanh λ. Note that the Lax potentials (3.1), or equivalently (3.12), do not transform as vectors under the Lorentz transformation (4.1). The Lorentz group in (1 + 1)-dimensions is a non-compact one-parameter group, namely SO(1, 1). Consider also an internal one-parameter group generated by the grading operator D, defined in (3.10), and acting on the loop algebra SU(3) as an automorphism, i.e.
The structure of the Lax potentials (3.12) is such that they transform as vectors under the diagonal subgroup, i.e. (the fields φ a , or equivalently ϕ a , are scalars under the Lorentz group (4.1)) In consequence, the curvature is invariant under such a diagonal subgroup, and so is the anomalous term appearing in (3.3), i.e.
Let us now analyse how the transformed Lax potentials a ± , transform under Ω. First we consider a − and we look at the second line of (3.15) and observe that However, as A (0) − = 0, this expression has to be cancelled by the transforms of [ b −1 , F 1 ] and of b −1 and so we see that it must be that a , and so we find that Ω (F a 1 ) = e λ F a 1 . This demonstrates the validity of (4.6). Looking at the terms in the next lines of (3.15) and using (4.6) we observe that, under the action of Ω, the last three terms of the third line of (3.15) get multiplied by e −λ . Thus, in order for the term [ b −1 , F 2 ] to cancel the image part of these three terms one needs that (4.7) Continuing this process recursively, order by order in the grades, one concludes that all F n have to be invariant under Ω, and so the group element g of the gauge transformation (3.14), i.e. satisfies Ω (g) = g.
In consequence the transformed Lax potentials a ± transform as vectors under the diagonal Lorentz subgroup in the same way as A ± , i.e. they satisfy Ω (a ± ) = e ±λ a ± . (4.9) Moreover, one of the consequences of the fact that all F n 's are invariant under Ω, is that from (3.14) we see that Ω ζ a were so chosen that the a − -component of the Lax operator is gauge transformed into the kernel of the adjoint action of b −1 , it follows that it depends only on x − -derivatives of the fields, and not of their x + -derivatives. So, from its transformation under Ω, we see that each parameter ζ (n) a of the gauge transformation has to be a polynomial in the derivatives of the fields with all of its terms containing only n x − -derivatives. Moreover, from (4.8) it then follows that Ω g F a 1 g −1 = e λ g F a 1 g −1 , and so each term on the r.h.s. of (3.19) under the action of Ω gets multiplied by e λ . Since Ω (b M ) = e M λ b M , this then implies that From (3.19) we then see that α (M,a) is a function of the parameters ζ (n) a , and so depends only on the x − -derivatives of the fields. Therefore, each term in α (M,a) has to contain exactly (M − 1) x − -derivatives of the fields. Looking at (3.25) we note that α (2,a) is indeed linear in the x − -derivative. Then from (4.4) and the fact that Ω (F a 1 ) = e λ F a 1 , it follows that Ω (X a ) = e −λ X a . In consequence, we have demonstrated that the anomalies of the charges, appearing in (3.24), satisfy solution, should be time independent and so the anomalies appearing on the r.h.s. of the second equation in (3.24) should vanish. But, from (4.11) it follows that, if the anomalies vanish in one reference frame, they vanish also in in any other reference frame connected by a Lorentz transformation. Thus, we conclude that all the charges Q (M ) , for any M in the infinite set of them defined in (3.23), are exactly conserved for any traveling wave solution and, in particular, they are conserved for the one-soliton type solutions. That is a highly non-trivial result since the densities of the anomalies, namely 2 a=1 α (M,a) X a , do not vanish in general when evaluated on a traveling wave solution. It is their integral over the whole one-dimensional space that has to vanish. Note also that for finite energy solutions of the equations of motion the space and time derivatives of the fields have to vanish at spatial infinity. In consequence, the α (M,a) and X a expressions have to vanish at spatial infinity, since as we have seen above, they are polynomials in the x − -derivatives of the fields (see (3.5)). So, for any one-soliton solution the densities of the anomalies 2 a=1 α (M,a) X a , are localized in space, and their space integral vanishes. One possible reason for the vanishing of such an integral is that the densities of the anomalies are odd functions of x, in the rest frame of the traveling wave solution. We have verified that this is exactly what happens for the one-soliton solutions of the theories (2.1) with potentials given by (2.23). In section 8 we explain how the one-solitons of such theories can be constructed numerically. One can then evaluate the anomalies on such solutions numerically. In figure 1 we plot the real and imaginary parts of the density of the anomaly β (2) , given in (3.30), as functions of x, in the rest frame of the one-soliton. The value shown there is for ε = 0.0005. Note that the complex density of the anomaly is indeed an odd function of x (the imaginary part is essentially zero; its infinitesimal values are numerical artifacts).
We have not understood yet the phenomenon of the cancellation of the anomalies. However, the conservation of the infinite set of charges for traveling wave solutions is clear from the argument based on the Lorentz transformation given above. In the case of traveling wave solutions like one-solitons this argument implies that the anomalies have to vanish irrespective of their densities being odd functions of x or not. For the case of two-soliton solutions (moving with different velocities) we have found that in all examples where the anomalies cancel, there is a space-time parity transformation playing a role. It would be interesting to investigate if there is a relation between the roles of the space parity in the case of one-solitons and the space-time parity in the case of two-solitons. In the next section we discuss the role of the space-time parity in the cancellation of the anomalies.

The parity transformation and charge conservation
The properties of field configurations, specially those describing one and two soliton solutions, under space-time parity transformations do seem to play a role in the vanishing of 'total' anomalies, i.e. when the anomalies are integrated not only over space but also over time. Consider a space-time parity transformation given by P : x evaluated on them, behave as follows under this parity transformation: where c 1 and c 2 are constants. In addition, we are interested in potentials that are invariant under the parity, i.e. (4.14) Note that (4.12) and (4.13) imply that where ∂ µ stands for the space-time derivatives, and where δ stands for the functional variations of the fields. Using (4.14) and (4.16) we find from (3.6) that P : Then, (4.15) and (3.5), give us that Next we check how the quantities α (M,a) and the anomaly densities transform under this parity transformation. To determine this we need to use another automorphism of the SU(3) loop algebra which involves the following order two outer automorphism of the finite simple SU(3) Lie algebra (σ 2 = 1)

JHEP05(2016)065
One can check that (4.19) is indeed an automorphism of the algebra SU(3) given in (A.1). This automorphism is insensitive to the value of the λ parameter of the loop algebra, and so we find that (see appendix A) Next we consider the combined action of the space-time parity P and this automorphism σ From (3.13) and (4.15) we see that Thus (4.20) gives us: Then applying (1 + S) to both sides of the second equation in (3.15) we get Let us recall that the procedure in (3.15) involved choosing the group element g and so also the F n 's in such a way that the new Lax potential a − was transformed into the kernel of the adjoint action of b −1 . Hence, as a result of this procedure A (0) − belongs to the kernel. But since σ, and so S, maps kernel into kernel (see (4.20)), we note that the l.h.s. of (4.24) belongs to the kernel. However, since the r.h.s. of (4.24) is the commutator of b −1 with something, it belongs to the image of the adjoint action of b −1 . Since image and kernel do not possess common elements (see (3.8)), then both sides of (4.24) have to vanish. Also, since σ, and so S, maps image into image (see (4.20)), it follows that (1 − S) F 1 belongs to the image, and so it cannot commute with b −1 . Thus it must be that Then applying (1 + S) to both sides of the third equation in (3.15) we find Using very similar arguments to those presented above one can also conclude that Continuing this process recursively, order by order in the grade expansion of a − , one concludes that all F n 's are invariant under S and so that S (g) = g. Next, using (4.20) and (4.28) one finds that Then from (3.19), (4.20), (4.21) and (4.29) one also finds that In consequence, (4.18) allows us to conclude that Thus we have demonstrated that the anomaly densities are odd under our parity transformation. This implies that if we integrate them on a rectangle with centre at (x , t) = (x ∆ , t ∆ ), (see (4.12)), they vanish, i.e.
Finally, takingx 0 → ∞, we find from (3.24) that the charges satisfy the mirror type symmetry So, if one considers the scattering of two one-soliton fields (which make a two-soliton solution satisfying (4.13) and (4.14)), the values of the infinite number of charges Q (M ) do vary in time, but after the scattering they all return to the values they had before the scattering. Since in a scattering process what matters are the asymptotic states, we see that the properties of such scatterings resemble those of an integrable theory, and that is why we call such theories quasi-integrable.

The exact soliton solutions of the integrable Affine Toda Models
The exact soliton solutions for the Affine Toda theories (AT) can be constructed by a variety of methods, all of which are based in one way or another on the zero curvature condition or the Lax-Zakharov-Shabat equation [13,14]. Among the several methods that have been used to study such theories, we have the inverse scattering method [15], Bäcklund transformations [16], the dressing transformation method [17][18][19][20][21][22][23][24][25], the solitonic specialization [26] of the Leznov-Saveliev solution [27], the direct Hirota method [28,29], and others (see [30] for a more complete account). The soliton solutions for the SU(N ) Affine Toda field theories were first constructed by Hollowood [31] using the Hirota method. The generalization of the construction to AT models associated to other algebras were presented in [32][33][34][35][36] using the Hirota method, and in [26,[37][38][39] using the Leznov-Saveliev method and the representation theory of Kac-Moody algebras based on vertex operators. The Hirota method is perhaps the most efficient procedure for constructing explicit analytical soliton solutions. However, it does not provide a way of finding the so-called taufunctions which are crucial for the Hirota method. Such functions can however be easily -17 -

JHEP05(2016)065
found using the dressing transformation method and the representation theory of Kac-Moody algebras based on vertex operators [38,39]. Therefore, the most efficient method for constructing soliton solutions is perhaps a hybrid procedure based on the dressing transformation and the Hirota methods as explained in [40,41]. An additional advantage is that this procedure can be easily adapted to be carried out with the help a computer package for algebraic manipulations. In fact, the magic of the Hirota method, which produces exact solutions by truncations of a formal series expansion, can be understood through the nilpotency of vertex operators in highest weight representations of the Kac-Moody algebras. In such representations the central element of these algebras cannot vanish, and so the Lax potentials, like the ones given in (3.1), have to live in the full Kac-Moody algebra and not only in the loop algebra. This requires the extension of the AT models to the so-called Conformal Affine Toda models (CAT) by the introduction of one extra field (or two if one wants conformal symmetry). Such an extension explains the need for one extra tau-function for the Hirota method to work, as compared to the number of fields of the AT models (see [33] for details). Therefore, for an AT model associated to a Kac-Moody algebraĜ, affine to a finite simple Lie algebra G, of rank r, there are r + 1 taufunctions τ j , j = 0, 1, . . . r, satisfying coupled partial differential equations, the so-called Hirota's equations. These equations are quadratic, cubic or quartic, in the tau-functions, depending on the connectivity of the Cartan matrix ofĜ (see [33] for details). Then an N -soliton solution is obtained through the Hirota ansatz for the tau-functions where δ j,(k,l) , etc are obtained, recursively, from the expansion the Hirota equation in powers of κ. In the expression above the Γ function stands for where z k = η k e −α k and v k = tanh α k , with α k real and η k = ±1. So, v k is the velocity and α k is the rapidity of the soliton k. The parameters ξ k fix the positions of the solitons at t = 0, but in some cases they can even be taken to be complex. The square of the parameter m k , and the first order vectors δ (1) 's, are determined from the first order (in κ) Hirota's equations, which lead to the eigenvalue problem [33] L ij δ Here L ij = l ψ i K ij , K ij are elements of the extended Cartan matrix of the affine Kac-Moody algebraĜ, and l ψ i are positive integers appearing in the expansion of the highest co-root ψ/ψ 2 , in terms of the simple co-roots α a /α 2 a , of G, i.e. ψ/ψ 2 = r a=1 l ψ a α a /α 2 a , and l ψ 0 = 1. Moreover, the parameters m k label, together with the topological charges, the species of the soliton solutions, and they also fix the masses of the one-soliton solutions. Note that -18 -

JHEP05(2016)065
the Hirota method fixes the moduli of m k , through (5.3), but not their sign. In fact, the sign of Γ (z k ) can be changed by flipping either the sign of z k or of m k , and this changes the sign of the topological charge of the solitons. So, such a flip of the signs turns a soliton of a given species into an anti-soliton of the same species and vice-versa. The higher order vectors δ (n) 's are determined, recursively, through the expansion of the Hirota equations in powers of κ [33][34][35].
The solitons have in general short range non-trivial interactions, but there is an interesting situation, first observed in [33], where the existence multi-soliton solutions, which are at rest with respect to each other was first pointed out and which, consequently, do not have static interactions. Such solutions are more easily constructed by considering the Hirota ansatz for one-soliton solution given by j being determined by (5.3), and Γ (z) being given by (5.2). The phenomenon of the existence of static multi-soliton configurations occurs whenever a given eigenvalue of the matrix L ij is degenerate. In general such degeneracy is related to a symmetry of the Dynkin diagram of G, but it can also be an accidental degeneracy. If a given eigenvalue of L ij is degenerate, the vector δ (1) j , associated to that solution, can be taken as a generic linear combination of the degenerate eigenvectors. This situation introduces new parameters into the solutions which can make the Hirota expansion truncate at higher orders. If one takes all but one such parameters to be zero one gets a one-soliton solution. However, by taking them different from zero one gets solutions which can be interpreted as multi-soliton solutions in which solitons are at rest with respect to each other. So, there are no static interactions among them which would have set them to move. There can be, however, interactions depending on their relative velocities. The number of solitons in a given static multi-soliton solution is equal to the degree of the degeneracy of the corresponding eigenvalue m 2 k (see (5.3)). The details of such construction can be found in [33], and the results can be summarized as follows: associated to the symmetries of the Dynkin diagrams one has static two-soliton solutions for the AT models associated to the algebras SU(N ), SO(2N ) (N a positive integer) and E 6 , and static three-soliton solution for the SO(8) AT model. Associated to accidental degeneracies one has static two-soliton solutions in the AT models associated to the algebras SO(6N + 2) and SO(6N + 1) (N a positive integer).
The list however does not end there. The higher order vectors δ (n) 's are determined by algebraic equations of the form [33] where λ is an eigenvalue of L ij , and V (n−1) j is a vector made out of the vectors δ (m) 's with m < n. Therefore, if the matrix L ij has two eigenvectors λ and λ , such that λ − n 2 λ = 0, then one can add to δ (n) j a term proportional to the eigenvector associated to λ . This brings an extra parameter into the solution which makes the Hirota expansion truncate at higher orders, and so gives the solution the character of a static multi-soliton configuration.

JHEP05(2016)065
The cases where such a behaviour had occured, were first discussed in [33] through a theorem which involves Galois theory in its proof, and they corresponded to the algebras SU(6 N ) and Sp(3 N ) (N a positive integer). Therefore the AT models associated to the algebras Sp(3 N ) present static two-soliton solutions, and those associated to SU(6 N ) can be described as representing static three-soliton solutions, since two of the solitons come from the degeneracy of any SU(N ) associated to the symmetry of its Dynkin diagram.
Finally we would like to point out that static two-soliton solutions can be constructed out of solitons and anti-solitons of the same species. As we have mentioned above solitons and anti-solitons of the same species are associated to the same eigenvalue m 2 k of L ij , since they correspond to opposite choices of the signs of m k (not determined by (5.3)). Therefore one can have in (5.1) the same eigenvector δ (1) associated to two exponentials of Γ's with opposite signs, i.e. the Hirota tau-functions are given by: Since the velocity is solely determined by z, there is a rest frame where such a solution can be made static.
The phenomenon of static multi-soliton solutions which was first observed in [33], has been also explored further in some papers, in particular in those dealing with the construction of multi-soliton solutions of the AT models [34,35,42]. More recently, the behaviour of the energy density of such static multi-soliton solutions has been studied in the case of SU(N ) AT models by one of us [43].

The solitons of the SU(3) Affine Toda model
Here we discuss the exact soliton solutions of the integrable SU(3) affine Toda model, which corresponds to the theory (2.1) with potential being given by (2.6). According to (2.4) the Euler-Lagrange equations for such a theory are given by For the case of the SU(3) affine Toda model the Hirota tau-functions τ j , j = 0, 1, 2, are defined by the following field transformation When one substitutes (5.8) into (5.7) one gets two equations for three tau-functions. However, as mentioned above one needs the conformal affine extension of the model to get the Hirota's equation for the tau-functions and so these tau-functions must satisfy: One can easily check that any solution of (5.9), by substitution into (5.8), leads to a solution of (5.7).

JHEP05(2016)065
For the case of SU (3) we have that the positive integers l ψ i introduced below (5.3) are all equal to unity. Therefore, the matrix L ij is the same as the extended Cartan matrix of SU(3) and is given by (5.10) Its eigenvalues are 0 and 3, with 3 being doubly degenerate. The zero eigenvalue leads to solutions traveling with the speed of light and do not correspond to solitons. We then have two species of one-solitons associated to the degenerate eigenvalue m 2 k = 3, and they can lead to static two-soliton solutions as explained above (see [33]). Therefore, from (5.2) we have solves the Hirota equation (5.9) and corresponds, in fact, to a vacuum solution. Therefore, using the Hirota ansatz (5.4) with δ (0) j = 1 one obtains two one-soliton solutions (of two different species). The one-soliton solution of the species-1 is given by: 12) and the one-soliton solution of the species-2 is with Γ (z) given by (5.11), and where ω is a cubic root of unity, different from unity itself. So we take ω = e i 2 π/3 , 1 + ω + ω 2 = 0. (5.14) From (2.6) and (2.8) we find that the Hamiltonian for the SU(3) AT model is given by where φ is defined in (2.2). Therefore, the discrete transformations: are symmetries of the Hamiltonian, if µ · α ∈ Z Z for any root α of SU(3). The vectors µ are called co-weights of the algebra, and they form the co-weight lattice. Such a lattice -21 -

JHEP05(2016)065
describes the degenerate vacua of the theory and gives rise to topological solitons. Indeed, the topological current is defined as One can check that the topological charges of the species-1 and species-2 one-solitons, given by (5.8) and (5.12) or (5.13) are given, respectively, by and Q top. = −η where η = ±1 is the sign introduced in (5.11). Moreover, λ a , a = 1, 2 are the fundamental weights of SU(3), and we have normalized the roots as α 2 a = 2. Note that the one-soliton solutions (5.12) and (5.13) are such that Therefore from (5.8) and (2.2) we see that Thus the complex conjugation of φ amounts to a sign flip and the interchange α 1 ↔ α 2 . In consequence, the Hamiltonian (5.15) is real when evaluated on the one-soliton solutions (5.12) or (5.13). Using the Hirota ansatz (5.1) one can construct also two-soliton solutions for the SU(3) AT model by solving the Hirota equations (5.9) recursively as explained above. By combining the two species of one-solitons one gets three types of two-soliton solutions. The species-11 two-soliton solution is given by The species-22 two-soliton solution is The species-12 two-soliton solution is given by:

JHEP05(2016)065
where Γ (z k ) is given in (5.11), and the quantities ∆ 11 and ∆ 12 are given by and In these expressions α a , a = 1, 2, are the rapidities introduced in (5.2), and related to the velocities by v a = tanh α a . Note that the two-soliton solutions (5.23), (5.24) and (5.25) satisfy the conditions (5.21) and (5.22), and so the Hamiltonian (5.15) is real when evaluated on them.
As explained in [33] and mentioned above, whenever the matrix L ij has degenerate eigenvalues one can construct static multi-soliton solutions. The eigenvalue 3 of the matrix (5.10) is doubly degenerate and so we can obtain a static two-soliton solution. Such a solution is obtained using the Hirota one-soliton ansatz (5.4) and it is given by where y 1 and y 2 are the free parameters used in the expression of δ (1) j which is a linear combination of the degenerate eigenvectors of (5.10). Similarly, this solution could have been obtained from the two-soliton solution (5.25) by setting v 1 = v 2 (or equivalently α 1 = α 2 ) and η 1 = η 2 . Note that the parameters y a , a = 1, 2, can be absorbed into the exponential as y a e Γ(z) = e Γ(z)+x (a) 0 , and so they are related to the positions of each onesoliton forming the static two-soliton solution. In fact (5.28) is a particular case of the static two-soliton solution for SU(N ) AT models given in eq. (4.13) of [33].
As we have explained in (5.6) one can easily obtain static two-soliton solutions by combining soliton and anti-soliton of the same species. For the species-1 solitons we get the solution

The parity properties
In our discussions of quasi-integrability in [3,4] we have tried to relate it to the parity properties of the field configurations. So let us briefly discuss here such properties of our two-soliton configurations even though our un-deformed model is fully integrable. We will later use these results when we consider the deformed models.
To consider the parity properties we define the following quantities: with Γ (z k ) defined in (5.2) and ∆ 11 and ∆ 12 defined in (5.26) and (5.27), respectively. We then consider the following parity transformation The two-soliton solution (5.23) can be rewritten as Thus, under our parity transformation, we have P : which implies that P : The two-soliton solution (5.24) can be rewritten as In this case we see that under our parity transformation we have P : which implies that P : The most interesting, 'mixed one', two-soliton solution (5.25) can be rewritten as In this case, we have very interesting transformations properties of the fields under our parity operation as we have P : which implies that P : The two-solitons (5.23), (5.24) and (5.25) are solutions of the SU(3) Affine Toda model which is an integrable field theory possessing an infinite number of conserved quantites. However, it is worth noting that these solutions satisfy the property (4.13) (see (5.35), (5.38) and (5.41)), and that the Toda potential (2.6) satisfies (4.14). Therefore, the properties of the SU(3) Affine Toda model support our criteria for quasi-integrability. We will show in our numerical simulations that such quasi-integrability properties are preserved by some special deformations of the SU(3) Affine Toda model.

The interplay between dynamics and CPT parity
As we have seen, the parity properties of a given two-soliton solution are crucial for the vanishing of the integrated anomalies, and so also for the asymptotic conservation of the charges. In section 5.2 we have shown that the exact two-soliton solutions of the integrable SU(3) Affine Toda theory possess the desired parity properties. When we deform this theory the two-soliton solutions cannot be easily constructed analytically, even in a perturbative power series in the deformation parameter ε. Therefore, we do not have much control over what happens to the parity properties of the deformed solutions, and so we have to study these solutions numerically. The experience we have gained through all the models we have studied so far, shows that if an exact two-soliton solution of the integrable theory presents the desired parity properties, the corresponding deformed two-soliton solution also presents asymptotically conserved charges, thus indicating that it preserves the parity properties. Note that it is not easy to check the parity of the deformed solution numerically; only the conservation or not of the charges can be checked. We do not understand why the deformed solution seems to preserve the parity properties, and certainly there is much still to be understood by studying the dynamics of the quasi-integrable solutions.
Here we present an argument that gives a hint for future investigations, and is, in fact, a modified version of the argument we have used in our previous papers [1][2][3][4][5]  Let us consider a wide class of deformed theories defined by the potential (2.7) and with v 0 given by (2.20). So, this is a much more general class than we consider in detail in this paper (they are defined by the potential (2.23)). The equations of motion (2.4) for this general class of theories can be rewritten as where φ is defined in (2.2), and where we have also introduced the quantity Instead of considering just the space-time parity transformation P , introduced in (4.12), let us combine it with the complex conjugation operation C, and so consider a CPT transformationP ≡ C P . We can then split the fields into their eigen-components underP as Analogously, we can split the equations of motion (6.1) into eigen-components underP as and so where n 1 and n 2 are defined in (2.20). We can draw some important conclusions from these formulas: 1. The following transformations are symmetries of the equations of motion (6.1) with k (±) a , a = 1, 2, being integers, (k a ) being even integers, and where w a are defined below (2.13), i.e. w a · v b = δ ab , a, b = 1, 2. Note however, that after such a transformation φ (−) may cease to be an eigenstate of CPT since we are adding to it a constant real vector. Here, λ a , a = 1, 2, are the fundamental weights of SU(3) satisfying α a · λ b = δ ab , a, b = 1, 2. Therefore, the CPT parity eigen-components are given by and so, indeed the φ (+) -component is trivial for such exact two-soliton solutions.
The consequences of all these facts, for the concept of quasi-integrability, are not fully clear to us yet. However, they hint at a conclusion that, perhaps, the dynamics of the deformed and un-deformed models favours the φ (+) -component to be trivial. As we have seen in our analysis of the previous sections, the φ (−) -component is the one with the desired properties for the cancelation of the integrated anomalies, and so also for the asymptotic conservation of the charges. Our numerical simulations, which we discuss in section 7, support these views. There is certainly a lot to be explored further and understood better on the role of parities in the concept of quasi-integrability. 7 Numerical support

General comments
In this and next section we present and discuss the numerical support for our results of the previous sections.
First we concentrate our attention on the undeformed models, i.e. the integrable SU(3) AT model, and then we discuss the results for the deformed model defined by the equations (2.4) corresponding to the potential (2.23). For the numerical work and to study the time evolutions we had to solve the equations of motion which are given by

JHEP05(2016)065
Note that if we put ε = 0 we recover the equations of the undeformed model i.e. equations (5.7). As these equations involve second order time derivatives of fields φ i we treat them as a Cauchy problem and so to find their solutions we need initial values of the fields φ i and the appropriate boundary conditions that the fields have to satisfy. Of course, for ε = 0 we have the analytical forms of the full solutions (described in section 5.1) and so we can test our numerical methods and procedures by comparing the numerically determined solutions to the analytical ones.

Numerical procedures
Our numerical simulations were performed using the 4th order Runge-Kutta method of simulating time evolution. As in [3] we experimented with various grid sizes and numbers of points and most of our simulations were performed on lattices of 40001 lattice points with lattice spacing of 0.0006 (so they covered the region of (-12.0, 12.0)). The time step dt was 0.0002. At the edges of the grid (i.e. for 11.90 < |x| < 12.00) we absorbed the waves reaching this region (by decreasing progressively the time change of the magnitude of the fields there).
To perform planned numerical simulations we needed initial field configurations but unfortunately, as mentioned above, we did not have their analytical form except for ε = 0 (i.e. in the undeformed case). So we determined them numerically. Thus we did not have their exact form but our initial numerically determined configurations, we believe, were sufficiently close to the exact configurations so that we could trust all our results.
The procedure we adopted to determine these intial configurations was similar to the one used in [5]. First we constructed approximate static one soliton field configurations. To do this we used static (5.12) configurations which we multiplied by a factor µ = 3 3+ε (see (2.26)) so that they satisfied the new boundary conditions. Then, using an incredibly small time step (dt = 1.0 * 10 −7 ) we evolved these configurations using the diffusive equations, which were like the proper equations of motion in which the second order time derivatives were replaced by the first order ones. This was achieved by using the equations given by (7.1) in which ∂ + ∂ − was replaced by 1 4 (∂ 2 x − ∂ τ ) where τ is an auxiliary diffusive 'time'. This replacement had the effect of making the configuration move towards the one that solved the static equations of motion. We evolved such configurations until their energy did not change much (in practice this was the accuracy to within 0.01% and the fields were essentially τ independent). We then used such almost exact one soliton configurations to construct two soliton fields (static and non-static configurations) by exploring their symmetries and sewing the fields together at x = 0 (i.e. by putting each soliton at ±x 0 ). For the non-static fields we used Lorentz symmetry of the model to determine the time dependence of the one soliton fields by calculating ∂ t φ i from the value of the ∂ x φ i of the static fields.
To be absolutely certain that this was a good procedure we compared this way of obtaining the initial conditions of the moving solitons to their exact expressions for the un-deformed model. When we evolved configurations from the initial conditions derived both ways -we could see no difference in the properties of fields at later times. Then with the initial conditions so obtained we performed many simulations for various values of ε. In these simulations we absorbed the energy at the boundaries. In consequence, the total energy was not conserved but the only energy which was absorbed was the energy of the radiation waves which reached the boundaries. Hence the total remaining energy was effectively the energy of the field configurations which we wanted to study. In fact, in most of the simulations the energy loss was extremely small showing that our model was really almost integrable; i.e. that the ideas of quasi-integrability are quite sound.

Undeformed model
First we present our results for the un-deformed model i.e. for the model with ε = 0.
Our first set of plots shows one soliton configurations. In figure 2 we present the plots of φ 1 . The two plots show the real and imaginary parts of φ 1 . The plots of φ 2 are very similar except that its phase rotates differently. This similarity comes from the symmetry of the field configurations mentioned earlier. Note that the plots of the real parts of φ i look very similar to those for the Sine-Gordon solitons.
As we said earlier the model possesses also two different classes of two soliton solutions. They are shown in figure 3. The plot in figure 3a shows the real part of φ 1 of the first class ('the mixed' one), while figure 3b shows the configuration of the second class ('of the two of the same' one). Because of the symmetry φ 1 = −(φ 2 ) we see that in both cases Re(φ 2 ) = −Re(φ 1 ) and the imaginary parts are the same.

General comments
Our method of generating the initial conditions by reflecting one soliton fields (when solitons ended up being far apart) gave essentially the same results as the method of taking  them from the exact solutions. The results of the simulations were essentially the same; moreover, they very closely followed the analytic expressions. Hence, the method was very reliable, at least, for ε = 0 and we hope it was also reliable for ε = 0 where we do not have any analytical solutions to compare our results to.
It is interesting to observe that the energy of all our solutions was real. This, as stated before, can be checked for the exact solutions (when this reality is guaranteed by symmetries) but this was also true in all our simulations, which somehow preserved these symmetries. In fact, this reality was also true for the energy density.
Our simulations have also established the stability of the two soliton systems. Of course, numerical simulations introduce some small perturbations but these perturbations -30 -

JHEP05(2016)065
did not lead to any instabilities. In some ways, these numerical errors were extremely small and random and so canceled each other on average. In part, this was probably due to extra symmetries which the simulations preserved.

Deformed model
As we have said earlier various deformations are possible and could be considered. However, we have looked mainly at the deformation given by the potential (2.23), where ε took both positive or negative values. This deformation preserves many symmetries of the original Toda system and so is very likely to lead to quasi-integrability. Indeed, like the Toda potential (2.6), the potential (2.23) is invariant under the interchange φ 1 ↔ φ 2 . Moreover, if φ * 1 = −φ 2 the energies of the underformed and deformed models are real. As explained in section 7.2, the deformed one-soliton solution was obtained through a diffusive relaxation method using the exact one-soliton solution (5.12) of the integrable SU(3) AT model, as a seed. Note that if one had used as a seed, the exact one-soliton solution (5.13) of the other species, the result would have been the same as taking the previous result and interchanging φ 1 ↔ φ 2 . In addition, due to the boundary condition (2.26), if one has the configuration of φ 1 for a deformed one-soliton, one can obtain the configuration for φ 2 just by flipping the sign of the φ 1 -configuration. Therefore, the deformed two-soliton solutions associated to the exact two-soliton solutions of species-11 and species-22, given in (5.23) and (5.24) respectively, are related by the interchange φ 1 ↔ φ 2 , and so the numerical simulations are essentially the same. Therefore, we treat them as just one case which we refer to as two of the same type. On the other hand, the deformed two-soliton solution associated to (5.25) we call a mixed case.

Results -static cases -the 'mixed case'
Here we discuss our results corresponding to the case of two solitons of the mixed case (ie those described by φ 1 and φ 2 whose real parts are shown in figure 3a. First, we have looked at the static case. When the solitons were too far away from each other they did not interact and they did not move. In figure 4 we produce plots of energy densities obtained for ε = 0.01 at two values of time (t = 0 and t = 1000.0) The solitons were initially placed at ±6 and it is clear that at t = 1000 they are still there thus we see that the solitons were initially too far apart to move.
So we started the simulations with the solitons initially placed closer together. One soliton was placed at x = −1.5 and the other one at x = 1.5. Our results can be summarised as follows. All the plots give the trajectory of one soliton (the one placed initially at x = −1.5, the other one followed a similar trajectory -reflected in x = 0): • The two solitons for ε = 0 appear to be stable and they do not move significantly.
• For ε < 0 we observe attraction followed by repulsion resulting in interesting oscillations.
In figure 5 we present a plot of 'the motion' of our x < 0 soliton for ε = 0.0. We note essentially no motion, as to be expected from the analytical results. The small 'motion'  corresponds to the movement by only two lattice steps in t = 3000 units of time and it is very likely a numerical artifact (we did not take the exact analytical solution but a field obtained by 'sewing up' two one soliton expressions). In figure 6 we present trajectories of solitons for three simulations with negative ε and in figure 7 two simulations for ε > 0. All these figures clearly support our claims made above. Note that for negative values of ε the frequency of oscillations increases with the increase of |ε|, and in fact, as can be seen from figure 6b the oscillations gradually generate a small (numerical) instability which later destabilises the process. Moreover, in all oscillations the solitons come close together and then bounce back. Looking at the plots of the energy density of the solitons we find that in all these simulations the solitons never come closer than r min ∼ 0.5 + 0.5 = 1.0, so it would appear that they never come on top of each other (before they bounce back). This is further supported by the fact that the fields φ 1 and φ 2 look the same at all times (i.e. during the oscillations). We have tried to see what happens when we start with the fields initially further apart or for more negative values of ε. In all the cases looked by us the solitons moved down to about the same minimal distance between them and then bounced back; the only difference was the period of oscillations which increased with the decrease of the magnitude of ε and/or the increase of the initial separation between the solitons.
Can they ever come on top of each other (i.e. can r min get smaller or even become zero)? This is difficult to assess for static solitons as we would have to start with solitons much closer together but this would introduce small perturbations due to our procedure of 'sewing' two solitons together. The only way to study this would involve starting with solitons moving towards each other. This will be discussed in the next subsection. Before we do this let us say a few words about the anomalies. Of course, the ε = 0 case has no anomalies so here we present the anomalies, i.e. expressions only for β (2) (3.30) for ε = 0. In figure 8 we present the plots of the anomalies seen in two simulations for ε < 0. We clearly see that the imaginary parts of anomalies are negligible and that the real parts vary (and change when the solitons are close together) but then return to their original values. This is very much in agreement what we would expect based on the ideas of quasi-integrability. In figure 9 we present similar plots of the anomaly seen in simulation for ε = 0.001 (its trajectory is shown in figure 7a). Clearly, the anomaly is again essentially real and its (real) value is very small indeed (smaller by more than two orders of magnitude from its value for negative values of ε -this is of course, associated with the fact that solitons repel and never get very close to each other). In fact, the anomaly oscillates a little and then decreases further as the solitons move further away from each other.

Non-static cases
We have also performed many interesting simulations for various values of ε, velocity and initial positions of solitons. Here we discuss the two-soliton fields of the mixed case, and in the next section the other case.
When we sent the solitons towards each other two things could happen -solitons could reflect with or without a 'flip'. Here, by a 'flip' we denote the situation in which the two fields φ 1 and φ 2 swapped their shapes after the scattering. This 'swapping' refers only to their real parts as the imaginary parts stay the same. In figure 10 we present the plots of the real parts of fields when we had a reflection, and in figure 11 the similar plots for the case when the fields performed the 'flip'. We can try to relate this 'flipping' to the issue of the solitons coming on top of each other or not, which we alluded to in the previous subsection. In fact, all the case of the 'flipping' corresponded to the cases when the solitons got on top of each other. We have verified this in all the cases. We observed What about the anomalies? Our simulations showed that they were always very small and were essentially real. In figure 13  simulations shown in figure 10 and 11. We see that in the 'unflipped' case the anomaly does not change as much as in the 'flipped' one and so this case is more reliable.

The other class of 2 solitonic solutions
Next we present the results for the solitons from the other class, i.e. the one corresponding to figure 3b (two of the same type). In this case we always have a repulsion so below we present the results of only a few simulations.

Static case
We have performed several simulations (for several values of ε). The results are very similar so here we present 3 plots of the position of one soliton, initially placed at x = −1.5 (with the other soliton placed ast 1.5), for 3 values of ε. The results are shown in figure 15. We note that the repulsion increases with ε.

Solitons sent towards each other
We have also performed the simulations of solitons sent towards each other with various values of velocity. In each case we observed the repulsion (although with the increased velocity the solitons managed to get closer to each other). In figure 16 we present the plots of the trajectories of solitons (sent with velocity v = 0.1 towards each other) seen in simulations performed for several values of ε. As before, we plot the trajectory of one of the solitons and the other one moves symmetrically around x = 0. We do not see much difference in behaviour between all 4 plots. In figure 17 we present the plots of the anomalies seen anomalies are very small (even smaller in this case). Of course, this is due to the fact that the solitons never get very close to each other.

Further comments about our procedure
So far, in all above calculations, we have constructed the approximate (initial) two-soliton configurations by 'gluing' two one-soliton ones. However, as we have two fields φ 1 and φ 2 we have more possibilities for performing such a construction. For one soliton in the undeformed Toda model the fields φ 1 and φ 2 are related to each other by the symmetry mentioned in section 2. For two solitons we can construct the initial φ i fields by 'gluing' two one-soliton φ 1 fields into a two-soliton φ 1 field and doing the same for φ 2 fields or by taking the second one-soliton field by replacing φ 1 and φ 2 . Both resemble the undeformed exact two-soliton fields and so at first sight both procedures can be expected to give essentially the same configurations which would then be expected to evolve in the same way (whether the initial configurations started them at rest or at a velocity towards each other).
In fact, in the discussion in the previous section the initial fields were constructed using the first approach (two φ 1 's being used to construct a new φ 1 and similarly for φ 2 ).
We have performed simulations using the second method of construction (using both φ 1 and φ 2 fields to construct each of two soliton φ i fields) and the results were always the same. So our expectations were correct.

Further general numerical comments
Let us point out that the energy is very well conserved in all our simulations (and it always remains real). In figure 18a and 18b we present plots of the total energy seen in two simulations corresponding to figure 6a and figure 11. In the first simulation we see essentially no change in energy, the second one shows that after the scattering the solitons, which have already transformed themselves into new solitonic states, the new solitons radiate a little. A plot of the trajectory of one soliton seen in this scattering and the resultant transformation is also shown in this figure (18c). Note that such transformations do take place in all scattering for ε = 0. A similar scattering but for a negative epsilon is shown in figure 11. For ε = 0 the 'flip' can also take place but it does not lead to the transformation of the soliton. These transformations are very interesting and we hope to study them further. Related to these transformations are the oscillations described before and shown in figure 6. They are only seen when ε < 0 and when the solitons are initially placed close enough. For positive ε, solitons placed at rest repel, for ε = 0 remain at rest, and for ε < 0 their forces are more complicated leading to the observed oscillations. The number of oscillations depends crucially on the value of ε < 0. The larger the value of the |ε| the larger the frequency of the oscillations. Our solitons were initially placed at ±1.5 and they oscilled between this value and ±0.5. When we changed the initial value from ±1.5 to ±1.4 there was no change of the closest approach but this should be investigated further. However, this requires a lot of extra work so we leave it as one the things to do in future.
Our results suggest that for ε < 0 there may exist a bound state of two solitons but the fact that during the oscillation there does not seem to be any radiation being emitted shows that the problem may be a bit more complicated and so it also deserves further investigation.
Finally a few more words about the anomalies, which have always been real (the imaginary values are clearly numerical artifacts). In the earlier part of the text we talked about the importance of using the integrated over time anomalies but we have not their plots yet. So in figure 19a. we present the time integrated anomaly for the oscillation presented in figures 6a and 8a. We have looked in detail at the oscillations and they are not numerical artifacts, This had been checked by performing the simulations with different values of dt etc. The plot shows very clearly that the anomaly, on average does not change much. In figure 19b we have presented the integrated anomaly for the case described in figure 9. We note an initial small change followed by the stability. This, is of course, due to the fact that the initial solitons were too close to each other but in any case the change was extremely small. Incidentally, all our plots of anomaly should be multiplied by ε and in the case of integrated anomaly also multiplied by 0.1 ε. The reason for dropping these factors was due us to not wanting to have too small numbers in our plots.
The integrated anomalies for the moving soliton cases can be got from looking at plots in figures 13, 14 and 17 and summing the values in the plot. However to help the reader we present in figure 20 the corresponding plots.
All three plots show very small values (if one takes into account the factor of 0.1ε). The easiest to describe and the most reliable numerically is figure 20a. This is the case that corresponds to the no-flip scattering and our results do demonstrate that there is a change of the anomaly, though very small. This was seen in all other cases describing no-flip simulations when during the scattering the solitons never came closer than ∼ 0.5 to each other. The actual values of the 'flip' case described in figure 20b are somewhat unreliable; they are clearly small but the anomaly described already in figure 13 has very spiky behaviour and so our results can have some numerical errors. The same can be said about figure 20c (describing a typical scattering of two solitons of the same class). Though we believe in the basic features of our results the actual numerical values may not be very reliable. We feel that these aspects of the scatterings have to be studied further. This may involve calculating the conserved or quasi-conserved quantities directly and not through the additional anomalies, as when the solitons are extremely close the fields become very spiky and so may generate numerical errors. This, together with the other things mentioned above has to be postponed to further work.

Conclusions
In this paper we have discussed, in some detail, the results of our studies of the SU(3) Toda model in (1+1) dimensions and some of its deformations. First we looked at the undeformed model and studied some of its finite energy solutions. There were several of them, they all had real energy and all these solutions were stable. This was checked by performing numerical simulations and comparing the results of these simulations with explicit analytical solutions (numerical simulations introduce small perturbations and so could be used to study their stabilities). In our studies we looked at one and two soliton configurations. Amongst the solutions of the model there was one in which solitons remained at rest (i.e. the attractive and repulsive forces between them cancelled). This cancellation of forces is very reminiscent of what is seen in systems of monopoles in (2+1) dimensions and suggests the existence of a BPS condition which, so far, we have not yet been able to find.
We have also perturbed the models by introducing a small perturbation. The perturbation we have considered corresponded to the change of the angle between the root vectors of the root lattice. This changed the form of the potential V (φ 1 , φ 2 ) and it also changed the values of the vacua of the model. The perturbation made the model non-integrable and so we used it to see how its results fitted with our ideas on quasi-integrability. Of course to do this we needed our perturbations to be small. In our work we have looked only at perturbations described by ε and we varied ε between -0.1 and +0.5.
We have performed many such simulations concentrating our attention on studying the scattering behaviour of two solitons. However to do this we needed one or two soliton field configurations which we did not have. So, first of all, we determined numerically one soliton configurations. This was done, as described in sections 6 and 7, by taking one soliton configurations of the unperturbed (i.e. ε = 0) model and then perturbing them, so that the fields satisfied new boundary conditions, and then evolving them via a diffusive equation. Having determined the solutions of this equation (for various values of ε) we then constructed initial configurations for our simulations by 'tying' two one soliton configurations and, when we wanted to have moving solutions, boosting the solitons towards each other. Such a procedure was successfully used in, say, [5] and we have tested it on the undeformed model (i.e. with ε = 0). In the ε = 0 case the results of numerical evolutions of such static one soliton solutions were extremely close to the analytical solutions of the undeformed model (they were almost indistinguishable). Hence, at least for small ε, we are confident of our results.
Then we have performed many simulations with the solitons initially at rest. First we looked at the case describing two solitons of the mixed case. In this case we have found that for ε > 0 the solitons repel while for ε < 0 they attract and, of course, as knew originally, when ε = 0, the forces cancel and the configuration is static. The attractive case was found to be more interesting, as after the initial attraction, when the solitons got very close together, they started to repel and so the system oscillated. During the oscillations the field configurations always looked the same. The energy was well conserved and the anomalies were very small. Next we looked at the similar initial configurations but, this time, with the solitons initially moving towards each other with small velocities. For very small velocities nothing was very different; at larger velocities the solitons could come 'on top of each other'. In such cases, afterwards, the fields φ 1 and φ 2 'swapped' their form, and afterwards, the solitons moved away from each other (we had a genuine 'passing through each other'). For this to be the case we needed two fields, as then the rising field of one soliton in φ 1 ended up in -46 -JHEP05(2016)065 field φ 2 and vice-versa. This was observed in all cases for all values of ε (for sufficiently large velocities).
The results of our simulations bring out also an additional difference between ε = 0 and ε = 0. In the ε = 0 case the solitons after their scattering are the same as before the scattering. In the ε = 0 the solitons come out of the interaction region a little altered (in fact they oscillate and they move faster). This can be seen from figure 12a where the soliton after the scattering moves faster. This suggests to us that the ε = 0 models may have additional moving two-soliton solutions, but whether or not this is really the case, would require further studies.
We have also looked at the solitons of the second class and in all their cases the solitons always repelled. In all the scatterings, that we have looked at (even for solitons sent towards each other with some velocity), the solitons always repelled at some short distances. And this was true for all values of ε and, by this behaviour, the scattering recalled very closely the scattering of solitons in the Sine-Gordon model (unmodified or modified [1]). The anomaly also changed little. Thus, we note that our results, in addition to making some interesting observations about the properties of solitons of the unmodified SU(3) Toda model, also provide further support for the concept of quasi-integrability (as all the anomaly effects in the modified models were always very small). Moreover, our results have also indicated that the static solutions of the unmodified model changed as one introduced our perturbations. For positive values of ε the solitons repelled and for negative values of ε they got modified to interesting oscillating fields.
where ω is a cubic root of unity different from unity itself, i.e. ω 3 = 1; 1 + ω + ω 2 = 0; ω = 1. (A.5) The commutation relations of the loop algebra in such a basis can be easily obtained from their matrix construction given in (A.4).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.