Quasi-integrability of deformations of the KdV equation

We investigate the quasi-integrability properties of various deformations of the Korteweg-de Vries (KdV) equation, depending on two parameters $\varepsilon_1$ and $\varepsilon_2$, which include among them the regularized long-wave (RLW) and modified regularized long-wave (mRLW) equations. We show, using analytical and numerical methods, that the charges, constructed from a deformation of the zero curvature equation for the KdV equation, are asymptotically conserved for various values of the deformation parameters. By that we mean that, despite the fact that the charges do vary in time during the scattering of solitons, they return after the scattering to the same values they had before it. That property was tested numerically for the scattering of two and three solitons, and analytically for the scattering of two solitons in the mRLW theory ($\varepsilon_2=\varepsilon_1=1$). We also show that the Hirota method leads to analytical one-soliton solutions of our deformed equation for $\varepsilon_1 = 1$, and any value of $\varepsilon_2$. We also mention some properties of soliton-radiation interactions seen in some of our simulations.


Introduction
The objective of the present paper is to study the quasi-integrability properties of deformations of the Korteweg-de Vries (KdV) equation [12,5] that include as particular cases the so-called regularized long-wave equation (RLW), proposed by Peregrine [15] and T. B. Benjamin, J. L. Bona and J. J. Mahoney [2], and also the modified regularized long-wave equation (mRLW) introduced and studied by J. D. Gibbon, J. C. Eilbeck and R. K. Dodd [11]. Concretely, the model we consider involves a scalar field u satisfying the equation In which ε 1 , ε 2 , and α are real parameters, and where Note that the integrable KdV equation [12,5], corresponds to the case where the deformation parameters vanish, i.e. ε 1 = ε 2 = 0, and it is given by 1 As is well known, this equation describes non-linear waves in shallow waters traveling in the positive direction of the x-axis only. If one considers the linearization of the KdV equation one finds that its traveling wave solutions satisfy a dispersion relation of the form ω = k − k 3 , and so the phase velocity ω/k = 1 − k 2 , and group velocity d ω d k = 1 − 3 k 2 , which become negative, and in fact unbounded, for large enough k.
Motivated by this fact, Peregine [15] and T. B. Benjamin, et.al. [2] proposed the so-called regularized long wave equation (RLW) which corresponds to (1.1) for the case ε 1 = 1 and ε 2 = 0. The advantage of the RLW equation over the KdV equation is that the RLW equation yields a dispersion relation of the form ω = k/(1 + k 2 ), and so a phase velocity that is bounded and tends to zero for short wavelengths. Its disadvantage is that the RLW equation is not integrable and that it possesses only one analytical solution, namely the one-soliton solution. The two and three-soliton solutions for RLW are only known numerically and were constructed by Eilbeck and McGuire [8,9]. 1 The KdV equation as presented here is not in its standard form, but one can perform a change of variables to get the standard form of the KdV equation (see, for instance, the notes by M. Dunajski for more information [7]). The same applies for the RLW and mRLW equations. In fact, to get the original notation used by the authors one should choose α = 12, in (1.4), and α = 8 in (1.5).
The mRLW equation, introduced by Gibbon, Eilbeck and Dodd [11], can be written in terms of the u field as and so it corresponds to (1.1) for the case ε 1 = ε 2 = 1. The linearized version of such an equation has the same dispersion relation as the RLW equation, and the same exact analytical one-soliton solution as RLW. The remarkable property of the mRLW equation, however, is that it also possesses analytic two-soliton solutions, even though it is not integrable in the sense of possessing an infinite number of conserved quantities.
The analytical one-soliton solution for the RLW equation, and the analytical one-and twosoliton solutions for the mRLW can be constructed using the Hirota direct method where the relation between the τ -function and the u-field is of the form u ∼ − (ln τ ) xt . Therefore, by changing the τ -function as τ → f (x) g (t) τ , does not change the solution for the u-field. So, integrating (1.2) one gets that w ∼ − (ln τ ) x + h (x), and v ∼ − (ln τ ) t + j (t). Consequently, as long as Hirota's solutions are concerned the integration "constants" h (x) and j (t), can be reabsorbed in the redefinition of the τ -function. Therefore, when constructing the soliton solutions, either analytically or numerically, it is convenient to work with the field q defined as Dropping the integration "constants" as explained above one can then write w x = − 8 α q xx and v t = − 8 α q tt . Replacing the u-field by the q-field into (1.1) one gets an equation for q which can in fact be written as the x-derivative of the equation q tt + q xt − 4 q 2 xt − 2 ε 2 q xx q tt + q xxxt − ε 1 (q xxtt + q xxxt ) = 0. (1.7) Therefore, any solution of (1.7) leads to a solution of (1.1), when the integration "constants" in w and v are absorbed as explained above. So in this paper we base our discussion on the study of (1.7).
First, we show in this paper that (1.7) admits an analytical Hirota one-soliton solution for the case ε 1 = 1 given by In terms of the u-field this solution takes the form (1.10) Note that (1.10) constitutes a one-parameter family of one-soliton solutions, labeled by the deformation parameter ε 2 . It interpolates between the one-soliton solution of the RLW model (for ε 2 = 0 and α = 12) constructed in [15,2], and the one-soliton solution of the mRLW model (for ε 2 = 1 and α = 8) in [11]. Note also that, as pointed out in [11], the expressions for the one-soliton solutions of the RLW and mRLW, indeed look the same for the u-field, when the rescaling parameter α is chosen as above, i.e. α = 12 for RLW and α = 8 for mRLW.
The analytical Hirota two-soliton solution, however, exists only for the case ε 1 = ε 2 = 1, and that is the solution constructed in [11] given by q = ln 1 + e Γ 1 + e Γ 2 + A 12 e Γ 1 +Γ 2 ; ε 1 = ε 2 = 1, (1.11) where and (1.13) The previous paper by two of us (FtB and WJZ) [3], looked at the scattering properties of various soliton-like configurations of the modified regularized long-wave (mRLW) equation, originally introduced and studied by J. D. Gibbon, J. C. Eilbeck and R. K. Dodd [11], and given by equation (1.7) for ε 1 = ε 2 = 1. As stated in the aforementioned paper, the three soliton-solutions are not known and, in fact, they cannot be found using Hirota's method.
The main results of the previous paper involved the observation that the scattering of two initially well-separated solitons, which, in fact, is described by an exact solution of the equations of motion (see equation (1.11)), could be approximated very well by the numerical time-evolution of a linear superposition of two one-soliton solutions. In fact the approximation was so good that it was virtually impossible to see any difference between them. Then this paper also extended this idea to the simulation of various three-soliton interactions by numerically evolving linear superpositions of three (initially well-seperated) single-soliton solutions. Such scatterings were also found to be very elastic (in the sense that the solitons preserved their shapes and velocities and there was no visible loss of radiation). Moreover, the phase shifts in the scatterings of three solitons were, again, very well approximated by the sums of the successive two two-soliton scatterings. This property holds for integrable models and the fact that it holds also for the mRLW equation, which is not an Hirota integrable system, suggests that the mRWL equation may be 'close' to being integrable.
To test this further we have decided to extend the investigations presented in the previous paper by considering this model as a perturbation of the integrable KdV equation since then we can look at the conserved quantities of the integrable model and see how they change when we perturb this model to become the mRLW model, or for that matter any other model 'nearby'. We discuss all of this in more detail in the next section.
The paper is organized as follows: In the next section we introduce our Lax-Zakharov-Shabat equation which we then use to study the quasi-integrability of the model described by equation (1.1) for various values of ε 1 and ε 2 . We perform the usual gauge transformations and finally define the charges which are truly conserved when ε 1 = ε 2 = 0 (i.e., for the KdV model).
In section 3 we introduce the parity argument for the travelling wave solutions of (1.1) and show that if the multi-soliton field configurations of the model possess this symmetry, then the charges are quasi-conserved. The next section discusses the Hirota method of finding solutions of some nonlinear equations and appling it to our equation (1.1). In it we show that this method gives a one soliton solution for any ε 2 with ε 1 = 1 being constant. The next two sections discuss the analytical quasi-integrability of the mRLW model and present arguments for the integrability based on the eveness parity properties of the multi-soliton functions. Section 7 discusses the soliton solutions of the KdV equation, obtained via the Hirota method, and the parity properties of two-and three-soliton solutions of this model. The lengthy section 8 presents some results of our numerical simulations. These simulations were performed using a specially constructed numerical program based on implicit and explicit methods of solving (1.7). The calculations involved studying the time evolution of various field configurations initially corresponding to two-or three-soliton systems and then checking whether the observed results supported quasiintegrability of the model for various values of ε 1 and ε 2 . We finish the paper with our conclusions and a few short appendices presenting more information about our numerical techniques and providing some details on the construction of conserved or quasi-conserved quantities discussed in section 2 of the paper.

The anomalous Lax-Zakharov-Shabat equation
Our motivation is to study equation (1.1) in the context of quasi-integrability as a deformation of the integrable KdV equation expressed by equation (1.3). To this end, we start by constructing the quasi-zero curvature equation, and subsequently produce the quasi-conserved charges from it. To do this we introduce the Lax potentials A x and A t , as where the generators b 2n+1 and F n , are defined in the appendix B, and where 3) The curvature of the Lax connection is given by where Y = 0 corresponds to the equation of motion of our deformed model (1.1) and where X ≡ α 6 Thus, if we now assume that u is a solution of the equations of motion (i.e., Y = 0), then X represents the anomaly of the zero curvature equation. In the case of the KdV equation (i.e., ε 2 = ε 1 = 0), the anomaly vanishes resulting in the well-known infinite number of truly conserved charges.

Quasi-conserved charges
In order to construct the quasi-conserved charges we follow the procedure discussed in [10], which is an adaptation to quasi-integrable theories, of the so-called abelianization procedure for exactly integrable theories [6,14,1]. The key ingredient is that the generator b 1 appearing in A x in (2.1), is a semi-simple element of the loop algebra and so splits this algebra into its kernel and image under its adjoint action as explained in the appendix B. Then, we perform a gauge transformation to rotate the Lax potentials into an infinite abelian subalgebra of the sl(2) algebra, generated by b 2n+1 . The gauge transformation transforms the potentials as where the group element g is chosen to be Performing the expansion we find that a x can be expressed as Using the algebra defined by equations (B.25) to (B.30), we can write a x as where we have written down all the terms of the same grade on the same line and the grading is defined by the grading operator d given in (B.31). Note that the parameter ζ k , in the expansion (2.11), multiplying the commutator [b 1 , F −k ], and so first appears among the terms of grade −k + 1. Thus, one can choose the parameters ζ k , recursively, to cancel the image part of a x at the grade −k + 1, i.e. its component in the direction of F −k+1 . The expressions for the parameters ζ k , obtained this way, are given in the appendix C, and they are polynomials in u and its x-derivatives. Having chosen the parameters ζ k this way, the expression for a x becomes where the first values of a (−2n−1) x are given by Note that in our gauge transformation we have not used the equations of motion Y = 0, with Y given in (2.5). In the transformation (2.8) of the A t component of the Lax connection, the group element g is already fixed, but we still can use the equations of motion Y = 0 to simplify it, i.e. we are performing an on-shell gauge transformation. The 'on-shell result is then given by (2.17) In the appendix C we give explicit expressions for the first few quantities a . Note that, due to the anomaly X, the quantities c (−n) t do not vanish, and so the potential a t is not really rotated to the abelian sub-algebra generated by b 2n+1 . That is a difference with respect to the case of exactly integrable theories, but it will not be a concern for us as we explain below.
The 'on-shell gauge-transformed curvature then becomes and so it takes the form where we have assumed that the equations of motion are satisfied (i.e., Y = 0), and where In addition we have that β 0 = 1, and β −m = 0 for m = 1, 2, 3, 4. The first nonvanishing β i is β −5 , and it is given by From the commutation relations equations (B.25) to (B.30) one can deduce that the commutators of a given element b 2n+1 with any element of the algebra never produce an element of the abelian sub-algebra generated by b 2n+1 . Therefore, the commutator [a x , a t ] does not contain any terms in the direction of the b −2n−1 generators, since a x lies in this abelian sub-algebra. Thus, if we now consider only the terms in the direction of the b −2n−1 generators of the gauge transformed curvature, we find that which can be rewritten as Since all the terms of the parameters ζ n depend on u and its x-derivatives, and u → 0 when x → ±∞, we see that g → 1l as x → ±∞. This implies that and so from (2.2), we get that a and so the quantities Q (−2n−1) are candidates for our quasi-conserved charges.
It is important to point out that the Lax potential A x , given in (2.1), does not depend upon the deformation parameters ε 1 and ε 2 , and so it is the same as the Lax potential for the integrable KdV equation. Therefore, the charge densities a (−2n−1) x and, consequently, the charges Q (−2n−1) , are the same as those for the KdV theory. The dependence upon the deformation parameters in (2.27), comes only from the anomaly X. Furthermore, we note that the lowest of these charges is exactly conserved, i.e.
Indeed, from (2.20) one sees that γ (−1) = 0, and this implies that Q (−1) is conserved and so u is a density of a conserved quantity for any value of ε 1 , ε 2 , and α.
The remaining charges are not properly conserved. However, as we show below, for some special soliton solutions such charges are quasi-conserved. By that we mean that in the process of a soliton scattering these charges do vary in time but, after the scattering, they all return to the same values they had prior to this scattering (i.e. the values before and after the scattering are the same). This property is not well understood yet, but we have found that it is accompanied by a space-time parity symmetry of the soliton solutions, and this is useful in trying to gain an understanding why the anomalies of the charges vanish when integrated over time during the whole scattering process. In section 3 we explain how this works for the case under consideration here.

The parity argument
Here we explain in detail how the quasi-integrability concept is related to the existence of some parity symmetries of the soliton solutions. Before we show how it works for some specific soliton solutions, let us give the general argument for the deformations of the KdV theory given by the equation (1.1).

Parity argument for the traveling wave solutions
Consider a class of traveling wave solutions u ≡ u x − ω k t + δ of the equation (1.1) which, of course, includes in it the one-soliton solutions as its particular cases. For a fixed value of the time t, we define the shifted space coordinate as and introduce the space parity transformation The only hypothesis that we are making here is that the traveling wave solution is invariant under parity, i.e. that Px (u) = u.
In consequence, we see that u xt = −ω uxx and so Px (u xt ) = u xt . Thus, the anomaly X given in (2.6), for such type of solutions satisfies Px (X) = X. (3.5) Note that the anomalies coefficients, γ's given in (2.20), always involve an odd number of xderivatives of the field u, and so Thus we see that and so the charges (2.24) are exactly conserved for all traveling wave solutions, and so also for the one-soliton solutions.

Parity argument for multi-soliton fields
Now we consider the space-time parity transformation around a given point (x ∆ , t ∆ ), of spacetime P : Again we make the hypothesis that the soliton solution we are considering is even under such a parity transformation, i.e. P (u) = u, (3.9) where u in (3.9) is evaluated on that particular soliton solution. In addition, since w = dt u, and v = dx u (see (1.2)), we see that Next we consider the following order 2 automorphism (i.e. σ 2 = 1) of the sl(2) loop algebra where d is the grading operator defined in (B.31). Then, We also consider the combination of these two operations, the parity and the automorphism: (3.14) One can check that the Lax operators (2.1) and (2.2), when evaluated on soliton solutions satisfying (3.9) and (3.10), satisfy and so the curvature (2.4) satisfies Let us rewrite (2.12) as Since where F −n is defined in (2.9), and so x lies in the kernel and Ω maps the kernel into the kernel, it follows that the l.h.s. of (3.20) also lies in the kernel. But the r.h.s. of the same equation clearly lies in the Image. Therefore, we conclude that both sides vanish. Since (1 − Ω) (F −1 ) cannot lie in the kernel we conclude that or Using (3.21), we get from the third line of (2.11) that lies in the kernel, we get that the l.h.s. of (3.22) lies in the kernel, and its r.h.s. lies in the image. Therefore both side have to vanish and so similarly to the previous case we conclude that Continuing this process we conclude that Ω (F −n ) = F −n , for any positive integer n, and so the group element g defined in (2.9) satisfies Indeed from (C.35) we observed that Thus we see that and so from (2.18) and (2.19) we get that Thus we see that all the anomalies appearing in (2.27) satisfy (see (3.11)) In consequence, integrating the r.h.s of this expression over a rectangle with centre at the origin of the (x ,t) coordinates (see (3.8)) we get Sendingx to infinity we find that the charges (2.24) satisfy the mirror type symmetry (see (2.27)) Thus, for any soliton solution satisfying the property (3.9), all the charges Q (−2n−1) are quasiconserved, and the sector defined by such soliton solutions constitutes a quasi-integrable sector of the theory.

The analytical Hirota soliton solutions
Next we investigate the existence of analytical soliton solutions for the deformed KdV equation (1.7), or equivalently (1.1). To do this we introduce the Hirota τ -function as for some parameter β to be appropriately chosen later. Putting (4.1) into (1.7) one gets the following Hirota's equation Next we take the one-soliton ansatz and insert it into (4.2). We consider the expansion of the equation in powers of η. In the lowest order (η 0 ) the equation is automatically satisfied. In the next order (η 1 ) we find that the equation is satisfied if either ω = 0 or The order η 2 terms show that the equation is satisfied if If we do not want k to depend on ε's then we have to take With that choice one can check that all the higher order terms, in powers of η, vanish automatically, showing that the truncation of the series leads to an exact solution. Thus, we see that we have found a family of one-soliton solutions, parameterized by ε 2 , and given by where we have absorbed the η parameter into δ, by writing η = eδ and shifting δ +δ → δ. This is the solution given in (1.10). Note that this solution, written in terms of the q field, takes the form: Note also that, as we have stated in the introduction, for ε 2 = 0 we get the RLW one-soliton solution, and for ε 2 = 1, we get the mRLW one-soliton. In between we get a whole new family of one-soliton solutions. The one-soliton solution for ε 1 = 1, would need k to depend on the deformation parameters ε 1 and ε 2 (see (4.5)), and so it appears to be unphysical.
We have also checked that applying the above Hirota-type procedure for a two-soliton type ansatz, i.e. one that expands τ up to a second order in η, leads to a solution only for the case ε 2 = ε 1 = 1 and β = 1. In such a case, the Hirota's equation (4.2) to Its two-soliton solution is and This is precisely the two-soliton solution given in (1.11)-(1.13), and first found in [11]. One can check that by replacing ω i given in (1.12) and putting it into (1.13), one finds that A 12 given in (4.13) is the same as the one in (1.13).

The analytical quasi-integrability of the mRLW theory
In subsection 3.1 we have shown that if the field u evaluated on a one-soliton is even under the space parity (3.2) (see (3.3)), then the anomalies vanish (see (3.7)), and so all the charges Q (−2n−1) , introduced in (2.24), are exactly conserved. Note that the family of one-soliton solutions (4.7) are even under the parity (3.2), and so this infinity of charges are exactly conserved for such one-soliton solutions.
Let us now analyse the two-soliton solution (4.11)-(4.13). Denoting we can write the two-soliton tau-function (4.11) as where we have defined Note that if k 1 = k 2 we see that A 12 = 0 and so the two-soliton solution reduces to a onesoliton solution. Therefore, for truly two-soliton solutions, i.e. when k 1 = k 2 , z + and z − can be considered independent space-time variables, i.e. they are linearly independent combinations of x and t. Then we have that From (1.6) and (4.1), for β = 1, we find that the two-soliton solution (4.11)-(4.13) for the u-field can be rewritten as (see (5.2)) To find the point of space-time around which we perform our parity transformation, we consider the linear relation among z ± and x and t, i.e.
Thus, the space-time parity transformation P : is of the same form as (3.8), and from (5.5) we see that the field u evaluated on the two-soliton solution of the mRLW model is invariant under this parity, and so it satisfies the hypothesis made in (3.9), i.e. that P (u) = u. So, as shown in subsection 3.2, all the charges Q (−2n−1) , defined in (2.24), satisfy the mirror symmetry described in (3.30), and they are asymptotically conserved in the scattering of two solitons.
We have thus presented an analytical proof of the quasi-integrability of the mRLW theory. It is worth adding that this is the first analytical proof of the quasi-integrability of a (nonintegrable) field theory in 1 + 1 dimensions (as in this case we have an analytical form of a two-soliton solution).

The parity versus dynamics argument
We shall now check whether the dynamics of the deformed model (1.1) favours or not the correct parity property of the field u (3.9), so as to make the infinite set of charges Q (−2n−1) , defined in (2.24), quasi-conserved. The plan of our approach is to write a given solution of (1.1) as a perturbative expansion around an exact solution of the integrable KdV equation. In addition, if that exact solution does satisfy the desired property (3.9) under the parity, we want to understand how the higher terms in the expansion behave under this parity. Our studies will show that the dynamics of the deformed model (1.1) does indeed favour the parity property (3.9) in a quite non-trivial and interesting way.
To make the perturbative expansion simpler we parameterise our two deformation parameters as Then (1.1) becomes Next we expand the field u as where each u (i) , in general, depends upon ξ. Thus, we find that Next we split all fields into their eigen-function parts under the parity (3.8). Thus and the same for w and v. Of course, the ε 0 part of (6.2) is the KdV equation, i.e.
The next orderε 1 equation is then given by Next we split this equation into its even and odd parts under the parity. The odd part is given by and the even part takes the form Let us now suppose that the zeroth order solution is even under the parity, i.e. that In this case we find that In consequence, we see that (6.9) becomes The order ε 2 part of (6.2) takes the form Splitting (6.15) into its even and odd parts one gets that the odd part becomes while the even part takes the form Note that if (6.11), (6.12) and (6.19) are all satisfied, then these equations become Continuing that process it seems that we can always have a solution in which the u field is even under our parity. Such a solution of the perturbed model fits into the scheme presented in subsection 3.2 and so all the charges Q (−2n−1) (from their infinite set)are quasi-conserved, i.e. they satisfy the property (3.30). In this sense, the dynamics of the perturbed model favours the even u field, since as we saw there cannot exist pure odd u field solution. Of course, there can exist mixed solutions, i.e. us with an even plus a (perhaps small) odd parts.
We do not understand yet how one could fine tune the perturbed solution to be purely even under the parity. Note that the condition of making the odd part of the solution vanish at each order of perturbation, does not work like an initial condition. So, we cannot prepare the solution at a given initial time t 0 and guarantee that it will evolve in time keeping its evenness. The conditions discussed above are imposed at all times. However, as our analysis has shown that a purely odd solution cannot exist, perhaps there might be a mechanism making the odd part of the solution small (maybe its appearance is not energetically favorable). We have not found such a mechanism yet. However, our numerical simulations, which we describe below, show that if we start with a seed configuration which is even under the parity, the numerical evolution of the configuration, under the perturbed dynamics, essentially keeps the configuration even. This is a very intriguing property of our quasi-integrable field theories, and we shall describe it, in more detail, in the next sections.

The exact Hirota's soliton solutions for the KdV equation
In section 6 we have shown how the dynamics of the perturbed equation (1.1) favours the even property of the u field under the parity. Since our discussion involved a perturbative expansion around an exact solution of the KdV equation (1.3), in this section we discuss the properties of the exact multi-soliton solutions under the space-time parity transformation (3.8). In order to do that, we construct these solutions using the Hirota's method.
We introduce the Hirota's τ -function for the KdV equation (1.3) as Putting this expression into (1.3) we get the Hirota's equation for KdV in the form The one-soliton solution of this equation is given by and so So we see that the KdV one-soliton is invariant under the space parity Px, defined in (3.2).
The two-soliton solution of the KdV equation corresponds to To understand its parity properties we introduce and so we find that where k ± ≡ (k 1 ± k 2 ) /2. Thus, the KdV two-soliton solution is even under the parity transformation P : (z + , z − ) → (−z + , −z − ), which is of the same form as (3.8).
Next we look at the three-soliton solution which corresponds to To simplify the expressions we define ∆ ij by and so where we have defined and so z + = z − . We can put where . From (7.12) we find that ln τ three−sol = ln 2 + z + + ln F Therefore, from (7.1), we get For this to happen, the determinant of this 3 × 3 matrix should vanish, and this implies Note that if we choose any pair of k i 's equal, we reduce the three-soliton solution to a twosoliton solution. So, we need k 1 = k 2 = k 3 . One way of satisfying this is to choose δ i , i = 1, 2, 3, in such a way that Note that in such a case the z (i) − 's, given in (7.14), (and so z + ) become homogeneous in x and t, i.e. z (i) Therefore, the parity transformation − , and so z + → −z + . The condition (7.21) leads to the linear system and so Denoting we get from (7.12) that One can check that all these cases correspond to the three soliton colliding simultaneously at the origin x = 0 and t = 0. The conclusion is therefore that for the KdV three-soliton solution in which the three solitons collide at the same point in space-time, the field u is even under the parity (7.23), when evaluated on such a solution.

Numerical analysis of quasi-conserved charges
In this section we discuss results of some of oyr numerical simulations of equation (1.7) using various values of ε 1 and ε 2 . We particularly focus on investigating the corresponding first nontrivial quasi-conserved charge Q (−3) , defined in equation (2.24). The numerical scheme used to approximate equation (1.7) is very similar to the algorithm used to solve the mRLW equation [3]. We have adjusted this scheme appropriately, which is discussed in appendix A, for the simulations presented in this section.
We have used several different initial conditions depending on the values of ε 1 and ε 2 used for each simulation, and we discuss these initial conditions in more detail in the subsections below. We have performed the numerical experiments for two-soliton simulations for a range of values for Γ 1 and Γ 2 but, for consistency, all the plots presented in this section regarding two-soliton simulations are obtained using the same values regardless of the initial conditions used, and similarly for all the three-soliton simulations (Γ 1 , Γ 2 and Γ 3 ) presented in this section. For all the simulations presented in this section we have used a grid spacing h = 0.1 and time level τ = 0.001. More details on the variables used for the simulations discussed in this section can be found in appendix A.2.

Quasi-conserved charges of the mRLW equation
In this subsection we discuss the results of our investigations of the quasi-conserved charges of the mRLW equation. To perform the simulations of the mRLW equation, we have set ε 1 = ε 2 = 1 in equation (1.7).

Two-soliton solutions of the mRLW equation
In figure 1 we have plotted the scattering of two solitons described by equation (1.11) at various values of t, and in figure 2a we have plotted the time dependence of ∂ t Q (−3) for this interaction. Since the mRLW equation possesses an analytical form of the two-soliton solution, in figure 2 we have in fact plotted the values ∂ t Q (−3) for both the analytical and numerical simulation. To better illustrate how small the discrepancy is, we have added an insert of the region near the global maximum on a much smaller scale. We note that there is hardly any discrepancy between the analytical and numerical values. Now, comparing figures 1 and 2a, we see that as t → ±∞ the charge is conserved except for during the collision, where it first decreases and then increases back to the same value as before the collision. Figure 2b  analytically in section 5. The agreement of the analytical and numerical results is a good test of the numerical code that we are using.

Three-soliton solutions of the mRLW equation
In figure 3 we present a selection of plots (at different values of t) of the u field obtained in a typical numerical simulation of a superposition of three one-soliton solutions. This was done in this way since we do not have analytic three-soliton solutions of the mRLW equation. However, as was discussed in [3], such a field is a very good approximation to a solution and the solitons of the numerical three-soliton configuration behave as integrable (or quasi-integrable) solitons. So here we have tested this by looking at the behaviour of corresponding charges. Figure 4 shows the corresponding time dependence of ∂ t Q (−3) and t 0 dt ∂ t Q (−3) . As Q (−3) changes only during the collision, and has the same value in the asymptotic regions (see figure 4b), we conclude that our charge is also quasi-conserved for the three-soliton configuration.
In section 7 we have shown that the u-field for the exact Hirota three-soliton of the KdV equation is even under the parity transformation (7.23), when the three solitons collide all at the   same point. As we have discussed in section 6, one should therefore expect that the deformed three-soliton solution should keep that parity property, and so to have an infinity of quasiconserved charges given by (2.24) (see section 3.1). The numerical results presented in figure 4 show that the charge Q (−3) is indeed quasi-conserved for the mRLW equation, and so it confirms this expectation.

Quasi-conserved charges of the RLW equation
In this subsection we discuss similar topics for the RLW equation. This model, is described by equation (1.4), in which the conventional re-scaling of the u-field used in the literature sets α = 12. The model possesses the same exact one-soliton solution with the same dispersion relation as the mRLW equation, but it does not possess an analytical two-soliton solutions [4]. However, since we are simulating equation (1.7) with ε 1 = 1 and ε 2 = 0, we must perform the appropriate rescaling of q (i.e., q → 3 2 q) in order to obtain the analytical one-soliton solution. Thus, to construct the two-and three numerical systems we have used, as initial conditions for the numerical simulation, a linear superposition of the following analytical single-soliton solutions where n = 2 corresponds to the two-soliton simulation and n = 3 to the three-soliton simulation.

Two-soliton solutions of the RLW equation
The results of the two-soliton simulation for the RLW equation are shown in figure 5. Clearly, the figure shows some radiation (see the inserts in figures 5e and 5f), and this agrees with the results presented in [13]. This is expected since the initial conditions used for this simulation do not solve the RLW equation analytically. Note that the amplitudes of these solitons are bigger compared to the solitons of the mRLW equation due to the aforementioned rescaling of q. Furthermore, in figure 6 we have presented the plots of the time dependence of ∂ t Q (−3) and t 0 dt ∂ t Q (−3) seen in this simulation. The shapes of the curves in these plots are very similar to the ones found in  . We note that the small observed emission of radiation does not seem to have a visible effect on the quasi-conservation of the charge.

Three-soliton solutions of the RLW equation
Next we have looked at the three solitons systems of the RLW equation. In figure 7 we present the plots of the fields u seen at various times in their evolution, and in figure 8 the plots of the corresponding ∂ t Q (−3) and t 0 dt ∂ t Q (−3) . Note that the amplitudes of the solitons are again larger when compared with the solitons of the mRLW equation (see figure 3). This is due to the factor of 3/2, mentioned before. The three solitons behave in a very similar way to the solitons shown in figure 3 except that, around the time they collide, a small bit of radiation is emitted (see the inserts in figures 7d, 7e and 7f). Just as for the two-soliton simulation, this emission of radiation does not seem to have a visible effect on the quasi-conservation of the charge.
where n = 2 corresponds to the two-soliton simulation and n = 3 to the three-soliton simulation. We present the results of such two-and three-soliton simulations of equation (1.7), for various values of ε 2 = 1 while keeping ε 1 = 1 constant.

Two-soliton configurations
As before, first we present our results for two solitons. Figure 9 shows the time-evolution of the u field for a simulation with ε 2 = 0.5 and a simulation with ε 2 = 1.5. Figure 10a and 11a show the time dependence of ∂ t Q (−3) and t 0 dt ∂ t Q (−3) seen these simulations. The system for ε 2 = 0.5 emits a small amount of (visible) radiation during the interaction of the soliton fields (similarly to the two-soliton simulation for the RLW equation (see subsubsection 8.2.1)). This radiation is difficult to spot on the scale shown in figure 9, and so we have added inserts in plots 9e and 9f to show more clearly this radiation emitted during the interaction. Looking at figure 10a, we see that Q (−3) changes significantly only during the interaction of the two soliton

Three-soliton configurations
Next, we also looked at some three solitons systems which were constructed using equation (8.3), with n = 3, as the initial conditions. Figure 12 shows the time evolution seen in two simulations for the values ε 1 = 1, ε 2 = 0.5 and ε 1 = 1, ε 2 = 1.5. These plots show that the three solitons interact with each other in a way very similar to the mRLW and RLW equation (see subsections 8.1.2 and 8.2.2). Unlike for the two-soliton simulation, we now see slightly more radiation emitted in the ε 2 = 1.5 case than in the ε 2 = 0.5 case (see the insert in figure 12f). We also would like to consider other values of ε 1 and ε 2 . However, we have no analytical expressions of one-soliton solutions when ε 1 = 1. Therefore, we have used the analytical twosoliton solutions of the mRLW equation as initial conditions for the two-soliton configurations, and the linear superposition of three single-soliton solutions of the mRLW equation for the threesoliton simulations. In the next subsection we use the same values as those used in subsection 8.3 (for ε 1 = 1 and ε 2 = 1) to simulate the two-soliton configurations. We have repeated these simulations to study the effect of using the 'wrong' initial conditions. Then in the subsequent subsubsections we investigate the quasi-integrability for various values of ε 1 = 1 while keeping ε 2 = 1.    Figure 15a shows the initial condition at t = 0. Figure 15b shows these two soliton fields after they have evolved for 5 units of time. Looking at that figure, we note that both solitons have moved to the right and have left some radiation behind immediately after the start of the simulation, while their amplitudes have increased slightly. This is because the true soliton of such a system is probably significantly different the one described by the initial conditions used in this simulation. This implies that the field u changes from the form given by the initial conditions and settles to a true solution, and in doing so the two solitons leave behind some radiation. Then figure 15c shows the largest soliton field after it has interacted with the radiation emitted by the smallest soliton field. Figure 15d shows the solution at t = 15 when the largest soliton starts to interact with the smallest soliton. Figure 15e shows the actual soliton fields interacting with each other. Finally, figure 15f shows the field when the solitons have moved away from each other again. These changes are all seen in figure 16a; from the beginning (at t = 0 to t 2) the initial conditions settle down to the true solution of the equation leaving behind some radiation, from about t 6 to t 10 the largest soliton interacts with the radiation left behind by the smallest soliton. Finally, from t 16 to t 23 the two soliton fields interact, after which the solitons keep moving unhindered. Notice also that figure 16b shows that, after the the initial conditions settle down to the true solution, all the lumps cancel each other out.
For other values of ε 1 = 1 and ε 2 = 1 we have observed a very similar behaviour to the one described above. To illustrate this better, figure 17 presents three plots comparing ∂ t Q (−3) for ε 2 = 0.5 and ε 2 = 1.5 (see figure 17a), ε 2 = 0.7 and ε 2 = 1.3 (see figure 17b), and ε 2 = 0.9 and ε 2 = 1.1 (see figure 17c). We see that the larger |1 − ε 2 | becomes, the more intense are the changes of Q (−3) when the solitons readjust their form from the one described by the initial conditions and when the larger soliton interacts with the radiation left behind by the smaller soliton. This is expected since the larger the value of |1 − ε 2 | is, the larger is the difference between the initial conditions and the true solutions of the equation, which results in a faster rate of change of the charge Q (−3) . Furthermore, this implies that the radiation left behind by the soliton is larger, and so the interaction between the larger soliton and the radiation due to the smaller soliton is more intense. However, when the two actual soliton fields collide, we see that as |1 − ε 2 | becomes smaller, the peak of ∂ t Q (−3) converges closer to the peak shown in figure 2a, both in terms of amplitude as well as in terms of its location on the horizontal axis. The horizontal movement is due to the fact that when the solitons emit radiation, their amplitudes increase/decrease (depending on whether ε 2 < 1 or ε 2 > 1) which causes their velocities to increase/decrease as well.
To illustrate better these observed properties, figure 18 presents again the various plots shown in figure 15 corresponding to ε 1 = 1 and ε 2 = 1.5 but it also shows the results of the simulation corresponding to ε 1 = 1 and ε 2 = 0.5, in order to compare the two simulations. This figure also illustrates the fact that if ε 2 > 1, then the solitons' amplitudes and velocities increase as a result of the emission of this radiation. On the other hand, if ε 2 < 1, the solitons' amplitudes and velocities decrease. Figure 19 presents the corresponding values of Armed with the observation of the last subsection we present here results of some simulations for ε 1 = 1 and ε 2 = 1 where, just as in the previous subsubsection, the initial conditions were constructed from the analytical two-soliton solution of the mRLW equation. Figure 20 shows the u fields seen in the two-soliton simulation corresponding to ε 1 = 0.5 (respresented by the red curve) and ε 1 = 1.5 (represented by the green curve). We observe that when ε 1 < 1 the radiation emitted by the solitons has initially a positive amplitude, and the solitons' amplitudes and velocities increases as a result of emission of this radiation. On the

Conclusions
We have performed a detailed study of the quasi-integrability properties of various deformations of the KdV equation which include among them the RLW and mRLW equations. The charges were constructed by introducing an anomalous zero-curvature condition as explained in section 2. Note that only the time component A t of the Lax pair was deformed away from the KdV potentials. The space component A x was not deformed and so the charges Q (−2n−1) , that we presented in (2.24), are the same as the exact conserved charges of the KdV equation. The difference lies in the anomalies α (−2n−1) , also introduced in (2.24), which vanish for KdV but not for its deformations. However, we have shown that the deformations of the KdV theory considered in this paper present a very interesting property. Even though the charges are not conserved and do vary during the scatterings of solitons, they all return, after the scattering, to the values they had before the scattering. It is this property, the asymptotic conservation of the charges, that defines what we call a quasi-integrability of the theory. The mechanisms underlying such asymptotic conservation of the charges is still not fully understood. What we have learned so far is that the quasi-integrability goes hand in hand with some special properties of the soliton solutions under a specific space-time parity transformation, as explained in section 3, and that these properties lead to the asymptotic vanishing of the anomalies. In addition, another intriguing property that correlates with quasi-integrability is that if a soliton solution of the exactly integrable theory possesses the correct behaviour under this parity, then the dynamics of the deformed theories seem to preserve this behaviour under a perturbative expansion around the integrable theory. In section 6 we have presented the details of such an expansion and showed its compatibility with the expected parity behaviour of the solutions.
Our analysis of deformations of the KdV theory present two novel ideas in the study of the quasi-integrability. First, we believe it is the first time that the quasi-integrability ideas were tested for the scattering of three solitons, where the implementation of the parity transformation is much more involved. When the three solitons scatter in separated pairs, it is expected that the quasi-integrability of two-soliton scattering would be preserved, since the third soliton is away and not interacting with the other two. However, when all three solitons collide together at the same time and in a given position in space, the three body interactions come into play, and one cannot rely on the two body interactions to analyze the quasi-integrability. We have shown in section 7 that the exact Hirota three-soliton solutions of the KdV equation do have the correct parity properties when their solitons all collide together at the same point in space. Using the results of section 6 one could therefore expect that deformed three-soliton solutions would possess the same property. That is exactly what our numerical simulations have shown: the charges are indeed asymptotically conserved when the three solitons of various deformations of the KdV theory collide together. This provides a strong support for our working definition of quasi-integrability.
The other key result of the present paper is that it presents the first example of an analytical demonstration, and not only numerical, of the quasi-integrability of a truly non-integrable theory. This was discussed in section 5 where we showed that the two-soliton sector of the mRLW theory was analytically quasi-integrable. The proof of this result is valid for all charges from their infinite set. Our numerical simulations have confirmed it for the case of the first non-trivial charge . For next charges the results were also supportive; however, as the next charges involved expressions which contained terms with more derivatives, our numerical results had more numerical errors and so were less reliable. Clearly, better numerical techniques have to be developed to prove our claims more decisively (i.e. with the same accuracy as we had for Q (−3) ). This work lies beyond the present paper.
We must emphasise, however, that the strongest support for our working definition of quasiintegrability comes from the results of our numerical simulations. In all simulations reported in section 8 it has been shown that the charge Q (−3) is quasi-conserved, i.e. that it returns, after the scattering of two or three solitons, to the value it had before the scattering. This quasiconservation has been tested to a very good precision. In addition, our simulations have shown what had been observed before in other theories; namely, that the charge only varies in time when the solitons come together and interact. The charge remains constant when the solitons are far apart. Another important and intriguing result of our simulations is related to the observation that these simulations have shown that the interaction of radiation with the solitons also seems to respect quasi -integrability. In the cases, in which we did not have an analytical configuration to be used as a seed for these simulations, there was a reasonable amount of radiation sent out while the initial configuration gradually settled to a proper solution of the deformed theory. The radiation emitted by each soliton interacted non-linearly with the other solitons, and we observed that such a radiation-soliton interaction made the charge Q (−3) to vary quite reasonably in time. However, when the radiation and the soliton have got separated, the charge returned to its original value. So, we have observed that the radiation-soliton interaction also seems to respect quasi-integrability, but at this stage we do not have any clear idea why this property holds too. Our analytical parity arguments do not apply to such an interaction.
The numerical methods used in this paper involved the specially constructed numerical programs which were based on the original program of [8] and have been appropriately adapted to numerically approximate equation (1.7). More details are given in appendix A. The results of our numerical simulations supported very well our expectations. This was particularily true for the lowest charge but it was also true to higher charges. However, expressions for higher charges involve more derivatives of the fields and so our expressions are more prone to numerical errors. Thus, although their behaviour supported our claims, the reliability of these results is not as sound as for the ones we have included in this paper.

A.1 Numerical approximation of the deformations of the KdV equation
In order to perform numerical simulations of equation (1.7) we follow the techniques discussed by J. C. Eilbeck and G. R. McGuire [8]. The equation necessitates the introduction of implicit methods. Hence first of all we introduce a new field p(x, t) by so that equation (1.7) can be written as We see that the main difference in approximating this equation compared to the numerical scheme used in [3] is an extra term proportional to p xxx .
In order to adjust our scheme to this extra term, we first discretize in both x and t by taking a finite set of points x 0 , x 1 , . . . , x N and t 0 , t 1 , . . . , t K , where h and τ denote the step size in space and time. Furthermore, we denote the grid points as (ih, mτ ) ≡ (i, m), where i = 0, 1, 2, . . . , N and m = 0, 1, 2, . . . , K, and we employ the notation p i m ≡ p(ih, mτ ) and q i m ≡ q(ih, mτ ).
Finally, let v i m denote our approximation to p i m and let w i m denote our approximation to q i m .
Next we introduce the following central finite difference operators by their actions on v i m : and similarly on w i m . Applying these operators to equation (A.7) in a straightforward manner we get which can be rewritten as (A.12) Next we introduce the following matrices and vectors   where Note that we need 2 elements per boundary, which is a consequence of the term proportional to p xxx . Having introduced this notation, our problem reduces to solving the following matrix equation where we need to solve for vector B. In our algorithm we have used the well-known LU decomposition method for this problem, see for instance the book [16] for more information. Repeating this procedure for many time levels allows us to determine the numerical time evolution of a system.
It is not too difficult to verify that this scheme is both second-order accurate in τ and in h. Furthermore, this scheme is an extension of the scheme discussed in [3], and so that gives us an indication that we can trust the results of our simulations. Tables 1, 2 and 3 shows all the parameters used to produce the figures shown in this paper.    where m ∈ Z, and λ is the so-called spectral parameter. The algebra is then given by The important ingredient in such procedure is that the generator b 1 in A x given in (2.1), is a semi-simple element of the algebra in the sense that its adjoint action on the loop algebra, splits it into kernel and image without an intersection, i.e. where we have used the notation ∂ m x ∂ n t u ≡ u (m,n) (C.36)