Several physically important nonlinear partial differential equations (PDEs) are completely integrable by the Inverse Scattering Transform (IST) which can be viewed as the nonlinear analog of the Fourier transform (see, [1,3,56] and references therein). Arguably, the most well-known completely integrable PDEs are the Korteweg-de Vries (KdV), modified KdV (mKdV), nonlinear Schrödinger, and sine-Gordon equations. They model wave phenomena in a wide range of applications in modern theoretical and mathematical physics, including fluid dynamics, nonlinear optics, and plasma physics, to name a few.

“Complete integrability” is an elusive term [69] but completely integrable equations have remarkable properties and a rich mathematical structure. For instance, they possess infinitely many conservation laws and high-order symmetries. They admit a Lax pair where the nonlinear PDE is replaced by a pair of linear equations whose compatibility only holds on solutions of the nonlinear PDE. They can be viewed as infinite dimensional Hamiltonian systems with two or three different Hamiltonians. By applying a change of dependent variable, completely integrable PDEs can be transformed into equations that are homogeneous of degree (second or higher) in that new variable and be recast in bilinear form in terms of Hirota operators [28,34]. Most importantly, completely integrable PDEs have solitary waves solutions that maintain their shape and speed while propagating at a constant velocity, and soliton solutions made up of such solitary waves that collide elastically. Nowadays, the terms solitary wave and soliton are used interchangeably, without reference to the elastic scattering property.

To get an initial idea about the possible complete integrability of a PDE, one should check if the equation has the Painlevé property [1,3], meaning that its solutions have no worse singularities than movable poles. To do so, one runs the so-called Painlevé test [8] which is algorithmic and can be performed with our Mathematica code PainleveTest.m [5]. If a PDE passes the Painlevé there is no guarantee that it is completely integrable but it is more likely than not to have special properties, viz., conservation laws, generalized symmetries, and so on.

In this paper we use the concept of scaling homogeneity, to symbolically compute conserved densities, fluxes, higher-order symmetries, Lax pairs, bilinear forms, and Miura-type transformations for polynomial systems of PDEs, thereby testing their complete integrability in a variety of ways. To keep the paper accessible to non-experts, we only cover PDEs in $(1+1)$ dimensions, i.e., one space variable $\left(x\right)$ in addition to time $\left(t\right)$. To illustrate the steps of the computations we focus on the Gardner equation which comprises the KdV and mKdV equations as special cases. Therefore, the results for those two important soliton equations are obtained without extra effort. This paper’s purpose is to provide a review of integrability properties and solutions of the Gardner equation and illustrate the applicability of the scaling symmetry approach for investigating the complete integrability of polynomial nonlinear PDEs.

Scaling (or dilation) homogeneity is a feature common to many nonlinear PDEs and it provides an elegant way to find, e.g., densities and higher-order symmetries. Indeed, these scalar quantities can be derived as linear combinations with undetermined (constant) coefficients of scaling homogeneous polynomial terms in the dependent variables and their derivatives. Since the defining equations for these quantities are linear, the method boils down to solving linear systems for the undetermined coefficients.

For conservation laws, the time derivative of a candidate density must be the total derivative of an unknown flux. To test this type of “exactness” we use the Euler operator from the calculus of variations. An automated computation of the corresponding fluxes requires integration by parts on the so-called jet space. In essence, one needs to integrate expressions involving arbitrary functions. Unfortunately, computer algebra systems (CAS) often fail at that task. To circumvent the shortcoming, we use the homotopy operator from differential geometry which reduces the problem to a one-dimensional integral with respect to an auxiliary parameter.

The existence of a sufficient large number of non-trivial densities suffices to establish complete integrability of the given PDE. In most cases, the corresponding fluxes are not needed. However, for some equations, e.g., the Kadomtsev-Petviashvili equation [1,3], it is necessary to swap the independent variables to be able to rewrite the PDE as a system of evolution equations. Doing so, the roles of density and flux get interchanged and that is why we also show the computation of the flux.

The computation of higher-order symmetries is done with the same methodology without having to rely on the relationship between symmetries and conservation laws, as expressed through Noether’s theorem. The scaling homogeneity concept and the method of undetermined coefficients also apply to recursion operators which connect the generalized symmetries. Of course, the recursion operator is hard to compute because it involves total derivatives and anti-derivatives. Fortunately, the defining equation for the recursion operator is linear, which again reduces the problem to solving a linear system. By applying the recursion operator to a low-order generalized symmetry, one can generate all higher-order symmetries one after another. Thus the existence of a recursion operator provides proof that a nonlinear PDE is completely integrable. With the recursion operator one can build an hierarchy of completely integrable equations which have the same properties as the original nonlinear PDE. For example, the well-known Lax equation [19,28] which is of fifth order in x is the first higher-order generalized symmetry of the KdV equation [20]. When encountering a allegedly “new” nonlinear PDE of high-order in the literature one should verify that it does not belong to the hierarchy of symmetries of a well-studied PDE of lower order.

The computation of Lax pairs based on scaling homogeneity is more challenging because the defining equation is no longer linear and, then the method of undetermined coefficients leads to a nonlinear algebraic system with quadratic or affine linear terms. It is still quite straightforward to solve but, if parameters are present, finding the conditions on the parameters for which solutions exist is substantially harder. The knowledge of a Lax pair is essential to apply the IST which allows one to solve the initial value problem for a nonlinear PDE and, consequently, find solitary wave solutions and solitons. Getting a Lax pair is the first step in the application of the IST and Riemann-Hilbert method to find soliton solutions. Using the Lax pair one can also derive conservation laws, and many other properties of the nonlinear PDE.

The same thing happens when applying the scaling homogeneity method to compute the Miura and Gardner transformations. Here again one has to solve nonlinear systems for the undetermined coefficients. The Miura transformation, which connects the KdV and mKdV equations, had a profound impact on the discovery of conservation laws of both equations and, in general, the early development of the notion of complete integrable of PDEs and soliton theory.

To apply Hirota’s direct method for finding solitary waves and solitons, the given PDE must be cast in bilinear form. Finding that bilinear representation is a non-trivial task which requires some insight, experience, and often ingenuity as explained in [28]. Our scaling-homogeneity approach is still helpful but has its limitations, in particular, if finding a bilinear representation requires splitting an expression into two parts in such a way that each part can be written in bilinear form. This is exactly what must be done with the mKdV and Gardner equations and, as far as we know, there is no algorithmic way to do this. Once the bilinear representation is found, special solutions, including solitary waves and solitons can be computed algorithmically.

As we will also show in Section 9, the Gardner equation can be transformed into the mKdV equation by a Lorentz transformation. Consequently, any time new solutions of the mKdV equation are discovered one has additional solutions to the Gardner equation. The search for new solutions, in particular, for the defocusing mKdV equation is still an active area of research (see, e.g., [45,59,61]).

Over the last two decades, our research team has designed and implemented fast algorithms [7,18,53], to test the integrability of nonlinear PDEs based on the scaling symmetry method [19,20,27,29,30,32,55]. In addition, we implemented a simplified version of Hirota’s method [21,28] and other direct methods [6,10] to find exact solutions of nonlinear PDEs. As a matter of fact, the computations in this paper have been largely done with our software packages which are written in Mathematica syntax but could be adapted for other computer algebra systems such as Maple and REDUCE.

The paper is organized as follows. In Section 2 we introduce the Gardner equation and mention some of its applications. As shown in Section 3, the Gardner equation passes the Painlevé test which indicates that the equation likely has many interesting properties. In Section 4 we discuss the scaling symmetry of the Gardner equation as well as the KdV and mKdV equations. Section 5 is devoted to the computation of conservation laws. Section 6 and Section 7 cover the computation of generalized symmetries and the recursion operator that connects them. In Section 8, we turn our attention to the computation of a Lax pair for the Gardner equation. The derivation of a bilinear representation is covered in Section 9. The computation of the Gardner transformation (connecting solutions of the Gardner equation and KdV equations) is shown in Section 10.

With the important quantities of the Gardner equation being computed, we show in Section 11 and Section 12 how our results translate to the KdV and mKdV equations. Solitary wave solutions and solitons are covered in Section 13 and Section 14, where we show that the nature of the solutions of the Gardner and mKdV equations depends on the sign of the cubic nonlinearity. A brief discussion of our software packages used in this study is given in Section 15. Finally, in Section 16 we draw some conclusions and mention a few topics for future work.

The Gardner equation [17] is the nonlinear PDE where $u(x,t)$ is the dependent variable (or field) which is a function of space variable x and time t. The subscripts denote partial derivatives, i.e., $u_t={\displaystyle \frac{\partial u}{\partial t}},u_x={\displaystyle \frac{\partial u}{\partial x}},$ and $u_{3x}={\displaystyle \frac{{\partial }^{3}u}{\partial x^3}}.$ The parameters $\alpha$ and $\beta$ are real numbers for which the values $\pm 1$ and $\pm 6$ are frequently used in the literature.

The Gardner equation has the KdV and mKdV equations as special cases. Therefore, (1) is sometimes called the combined or mixed KdV-mKdV equation. Indeed, for $\beta =0$, the (1) reduces to the ubiquitous Korteweg-de Vries equation [37], which models, e.g., shallow water waves [26], ion-acoustic waves in plasmas [1,3,14] and many other nonlinear wave phenomena where solitons arise. When $\alpha =0$, (1) becomes the mKdV equation [1,17], which models internal ocean waves, electromagnetic surface waves, waves in plasmas, and more [1].

Eq. (1) appears for the first time in [42,44] in the context of the Miura transformation (see Section 10) which connects the mKdV and KdV equations. The Gardner equation has a long history [17,35] and many applications [70] in fluid dynamics, in particular, for modeling long internal water waves [22,41], the dynamics of undular bores [36], and waves in multi-species plasma physics [47,48,62,63]. There is a wealth of information available about the Gardner equation. The equation appears in most books about solitons and integrability (see, e.g., [26] which has a list of soliton books).

Without loss of generality we take $\alpha \ge 0$ in (1) because, if $\alpha <0$, replacing u by $-u$ would make the coefficient of $u{u}_x$ positive again. However, no discrete symmetries of $x,t,$ or u will flip the sign of the coefficient of $u^2{u}_x$ so the cases for positive and negative $\beta$ will have to be treated separately. Exact solutions of (1) with $\beta >0$, called the focusing Gardner equation, are quite different from these of the defocusing case where $\beta <0$.

We focus on (1) because it is a typical example of a scalar $(1+1)$-dimensional evolution equation, of order n in x and first order in t, with $n=3$ in (1). More importantly, many of the results for (1) lead to the corresponding results for the KdV and mKdV equations by setting $\beta =0$ or $\alpha =0$, respectively. Like the latter two equations, the Gardner equation is known to be completely integrable for both signs of $\beta$ but the solutions are quite different for $\beta >0$ and $\beta <0$. The latter case is the hardest to deal with. The solitary wave solutions and solitons for the focusing Gardner equation follow from those of a focusing mKdV equation to which the Gardner equation can be reduced. The so called “table-top” or “flat-top” soliton, corresponding to a large amplitude, only exists for the defocusing Gardner equation.

In some applications the coefficients $\alpha$ and $\beta$ are time dependent [22]. In that case (1) has only a couple of conservation laws (conservation of mass and wave action flux, see Section 5) and is non-integrable. The treatment of (1) with varying coefficients is outside the scope of this paper.

In this Section we check if (1) has the Painlevé property. Performing the Painlevé test is straightforward but usually involves lengthy computations. In essence, one investigates a Laurent series solution for (1), in which the coefficients $u_k(x,t)$ are analytic functions in a neighborhood of the singular manifold $g(x,t)$. Furthermore, $g(x,t)=0$ determines the poles and $g(x,t)$ is assumed to be non-characteristic (i.e., $g_x(x,t)\ne 0$). The negative integer $\sigma$ determines the leading order term, $u_0{g}^{\sigma }$, in (5). We summarize the main steps of the Painlevé test and refer to [8] for details.

So, at least in the neighborhood of $g(x,t)=0$, solution (9) is free of algebraic and logarithmic movable branch points. Apart from possible essential singularities (which the test is unable to detect), the movable singularities of its general solution are poles determined by $g(x,t)=0$. Note that (5) serves as a general solution because, as required by the Cauchy-Kovalevskaya theorem, (1) is of order 3 in x and that number agrees with 3 arbitrary functions $g(x,t)$, $u_3(x,t),$ and $u_4(x,t)$ in (5).

Also, notice that truncating the Laurent series at the constant level term yields an auto-Bäcklund transformation, because for an arbitrary $g(x,t)$, both $u(x,t)$ and $u_1(x,t)$ must then be solutions of (1) with $\beta >0$. Setting $u_1(x,t)=0$ in (14) motivates the Hirota transformation discussed in Section 9 where it is shown that only for $\beta >0$ one can find soliton solutions of (1).

In conclusion, the Gardner equation passes the Painlevé test and, therefore, has the Painlevé property which is a good predictor for the equation to have conservation laws, generalized symmetries, and so on. The investigation and actual computation of these quantities is based on a scaling symmetry of (1) which we will discuss next.

As stated, when $\beta =0$ the Gardner equation reduces to the KdV equation (2) which is scaling homogeneous under the scaling (dilation) symmetry where $\kappa \ne 0$ is an arbitrary scaling parameter. Indeed, replacing $(x,t,u)$ in the KdV equation in terms of $(\tilde{x},\tilde{t},\tilde{u})$ yields

A fast way to compute (15) is to introduce the notions of weight, rank, and uniformity of rank. The weight, W, of a variable is the exponent of $\kappa$ that multiplies that variable. Thus, based on (15), $W\left(x\right)=-1,W\left(t\right)=-3$, and $W\left(u\right)=2$. Equivalently, $W\left(\partial /\partial x\right)=1$ and $W\left({\displaystyle \frac{\partial }{\partial t}}\right)=3.$ The rank of a monomial is defined as the total weight of the monomial. For example, $\alpha u{u}_x$ has rank 5 since the parameter $\alpha$ has no weight. An expression (or equation) is uniform in rank if all its monomial terms have equal rank.

Now, if (15) were not known yet, it could quickly be found by requiring that the KdV equation is uniform in rank, yielding

Solving (17) gives $W\left({\displaystyle \frac{\partial }{\partial t}}\right)=3W\left({\displaystyle \frac{\partial }{\partial x}}\right)$ and $W\left(u\right)=2W\left({\displaystyle \frac{\partial }{\partial x}}\right).$ Since $\kappa$ is arbitrary, without loss of generality one can set $W\left({\displaystyle \frac{\partial }{\partial x}}\right)=1$, resulting in $W\left({\displaystyle \frac{\partial }{\partial t}}\right)=3$ and $W\left(u\right)=2.$ So, requiring uniformity in rank of a PDE allows one to compute the weights of the variables (and, hence, the scaling symmetry) with linear algebra.

Setting $\alpha =0$, the Gardner equation becomes the mKdV equation. Requiring uniformity in rank readily yields $W\left({\displaystyle \frac{\partial }{\partial t}}\right)=3$ and $W\left(u\right)=1.$ Hence, the mKdV equation is scaling homogeneous under the transformation

Because $W\left(u\right)$ is different for the KdV and mKdV equations, the Gardner equation (1) will not be uniform in rank unless we use a trick. We will let the parameter $\alpha$ scale with some power of $\kappa$. Doing so, we must solve yielding which expresses that (1) has scaling (or dilation) symmetry

Starting with conservation laws, in what follows we will use the scaling symmetry to compute important quantities related to (1) and thereby establishing its complete integrability.

A conservation law of (4) is an equation of the form where the dot means that the equality should only hold on solutions $u(x,t)$ of (4). $\rho$ is called a conserved density (or charge) and J is the corresponding flux (or current).

The scalar functions $\rho$ and J are functions of u and its partial derivatives with respect to x. All subsequent computations are done on the jet space which means that u and its partial derivatives with respect to x are treated as independent, and also all monomials such as $u^3,u_x^2,$ etc. The density and flux could also explicitly depend on x and t but we will not cover such exceptional cases.

If J vanishes at infinity (because u and its x-derivatives decay to zero), then is constant in time often referred to as a conserved quantity or constant of motion. Depending on the physical setting, the first few constants of motion express conservation of mass, momentum, and energy. Preserving these types of quantities plays in important role in testing the accuracy of numerical integrators.

In (22), ${\mathrm{D}}_t$ is the total derivative with respect to $t,$ defined by where ${\rho }^{\prime }\left[u_t\right]$ is the Fréchet derivative of $\rho$ in the direction of $u_t$ and M is the highest order of $\rho$ in x. In practice, one simply applies the chain rule for differentiation with respect to t treating $u,u_x,u_{xx}$, etc., as independent functions.

Likewise, ${\mathrm{D}}_x$ is the total derivative with respect to x, where N is the order of J.

Since (22) is linear in the densities (and fluxes) a linear combination of densities with constant coefficients is still a density. The matching flux would, of course, be a linear combination of the corresponding fluxes with the same constant coefficients.

Returning to (1), one can readily verify that

Indeed, (26) is (1) written as a conservation law. Next, (27) straightforwardly follows after multiplication of (1) by $2u$ and a bit of integration by parts. Clearly, (28) is far less obvious and will require a computational strategy.

Notice that the above densities are uniform in rank. Indeed, ${\rho }^{\left(1\right)}=u$ has rank 1, ${\rho }^{\left(2\right)}=u^2$ has rank 2, and ${\rho }^{\left(3\right)}$ is of rank 4. The corresponding fluxes are also uniform with ranks $3,4,$ and $6.$ As a matter of fact, the entire conservation laws themselves are uniform in rank with ranks $4,5,$ and 7, respectively. This comes as no surprise because the defining equation (22) is only non-trivial if evaluated on solutions of the PDE and therefore the densities, fluxes, and conservation laws “inherit” (or adopt) the scaling homogeneity of the given PDE (and all its other continuous and discrete symmetries for that matter).

It turns out that the list of conservation laws of (1) continues ad infinitum. The Gardner equation has infinitely many conservation laws which is a clear indication that the PDE is completely integrable.

Using the scaling symmetry (21), we will now show how to compute (28) which is the shortest possible density of rank 4 since it is free of any terms that could be moved into the flux.

The homotopy operator for variable $u\left(x\right),$ acting on an exact expression E of order $K,$ is given [31,54,55] by with integrand

In (36), $\left(I_{u}E\right)\left[\lambda u\right]$ means that once $I_{u}E$ is computed one must replace u by $\lambda u$, $u_x$ by $\lambda u_x$, $u_{xx}$ by $\lambda u_{xx}$, etc. Use (35), to get which already has the terms of $J^{\left(3\right)}$ but still with incorrect coefficients (and the opposite sign). Finally, using (36), which is exactly the flux in (28).

A second criterion to establish the complete integrability of (1) is to show that the PDE has infinitely many generalized (higher-order) symmetries. As we will see, the computation of a hierarchy of such symmetries [20] is algorithmic and easier than for conservation laws.

A scalar function $G(x,t,u,u_x,u_{xx},\dots ,u_{mx})$ is called a generalized symmetry of (4) if and only if it leaves (4) invariant under the replacement $u\to u+ϵG$ within order $ϵ$, that is, must hold up to order $ϵ$ on the solutions of (4). Consequently, G must satisfy the linearized equation where $F^{\prime }$ is the Fréchet derivative of F in the direction of G: where n is the (highest) order of F. In (40) and (42), one must not only replace u by $u+ϵG,$ but also $u_x$ by $u_x+ϵ{\mathrm{D}}_{x}G$,   $u_{xx}$ by $u_{xx}+ϵ{\mathrm{D}}_x^{2}G$, etc. As before, $\mathrm{I}$ is the identity operator. The total derivative operators, ${\mathrm{D}}_t$ and ${\mathrm{D}}_x$, were defined in (24) and (25), respectively. Since the defining equation (41) is evaluated on solutions of the PDE, the higher-order symmetries inherit the scaling homogeneity of (4).

As we will show, (1) has infinitely many generalized symmetries starting with

With regard to the weights in (20) the above symmetries are of ranks $2,4$, and 6, respectively.

Except for $G^{\left(1\right)}$, each of these symmetries leads to a completely integrable nonlinear PDE, $u_t+G^{\left(j\right)}=0,j=2,3,\dots$, where the weight of t increases as j increases. Notice that $u_t+G^{\left(2\right)}=0$ corresponds to (1) itself.

Although the computation of symmetries is algorithmic and technically simpler than for densities, it still requires the computation of a plethora of terms and is best handled with symbolic software [20]. As an example, we show how to compute $G^{\left(3\right)}$ of rank 6.

To prove that there are infinitely many higher-order symmetries we will compute the recursion operator which generates those symmetries sequentially starting from the lowest order symmetry $G^{\left(1\right)}$ in (44).

As expected, the recursion operator for the Gardner equation is a combination of the well-known recursion operators for the KdV and mKdV equations ([49,50], p. 312), where ${\mathrm{D}}_x^{-1}$ denotes the anti-derivative (or integral) operator. $\mathcal{R}$ connects the symmetries sequentially (without gaps in the simplest cases)

For example, which is $G^{\left(2\right)}$, and which is $G^{\left(3\right)}$. One can use computer algebra to verify that $\mathcal{R}$ is a hereditary operator [15].

If $\mathcal{R}$ is a recursion operator for (4), then the Lie derivative [27,49,66,67] of $\mathcal{R}$ is zero, yielding where $[,]$ and ∘ denote the commutator and composition of operators, respectively. $F^{\prime }$ is the Fréchet derivative operator defined as where n is the order of F in (4). ${\mathcal{R}}^{\prime }\left[F\right]$ is the Fréchet derivative of $\mathcal{R}$ in the direction of F. For recursion operators of the form where T is the total number of terms in $\mathcal{R}$, and $m_1$ and $m_2$ are the orders of U and V, respectively, one has

With regard to (52), $\mathcal{R}$ is a linear integro-differential operator [57,58] which naturally splits into two pieces [9,11,67], where ${\mathcal{R}}_0$ is a local differential operator and ${\mathcal{R}}_1$ is a non-local integral operator. ${\mathcal{R}}_0$ is a linear combination of monomials involving ${\mathrm{D}}_x$, u, and parameters with weight (if applicable) so that each monomial has the correct rank. Also, note that ${\mathrm{D}}_x$ will always be “propagated” to the right in ${\mathcal{R}}_0$. For example,

Based on the theory of recursion operators, ${\mathcal{R}}_1$ is a linear combination with constant coefficients of terms of the form where $G^{\left(i\right)}$ is a symmetry and ${\mathcal{L}}_u\left({\rho }^{\left(j\right)}\right)$ is a cosymmetry (Euler operator applied to a density ${\rho }^{\left(j\right)})$, selected such that ${\mathcal{R}}_1$ has the correct rank [12,67]. To standardize the form of ${\mathcal{R}}_1$, one propagates ${\mathrm{D}}_x$ to the left. For example, ${\mathrm{D}}_x^{-1}{u}_x{\mathrm{D}}_x=u_x\mathrm{I}-{\mathrm{D}}_x^{-1}{u}_{xx}\mathrm{I}$. The local and non-local operators are then added to obtain a candidate recursion operator.