Two-Field Integrable Evolutionary Systems of the Third Order and Their Differential Substitutions

A list of forty third-order exactly integrable two-field evolutionary systems is presented. Differential substitutions connecting various systems from the list are found. It is proved that all the systems can be obtained from only two of them. Examples of zero curvature representations with $4 \times 4$ matrices are presented.


Introduction
We use the term "integrability" in the meaning that a system or equation under consideration possesses a Lax representation or a zero curvature representation. Such systems can be solved by the inverse spectral transform method (IST) [1,2]. Exactly integrable evolution systems are of interest both for mathematics and applications. In particular, systems of the following form where a is a constant, excite great interest since about 1980. The paper [3] is devoted to construction of systems of the form (1.1) among others. Nine integrable systems of the form (1.1) and their Lax representations have been obtained in the paper. In particular, it contains a complete list of three integrable systems (1.1) satisfying the conditions a(a − 1) = 0 and ord(F, G) < 2.
Here ord = order, ord f < n means that f does not depend on u n , v n , u n+1 , v n+1 , . . . . Here and in what follows, the notations u n = ∂ n u/∂x n , v n = ∂ n v/∂x n are used. Two of the three mentioned systems can be written in the following form and the third system is presented below (see (3.46a)). System (1.2) was found independently in [4] and the soliton solutions were constructed there. This system is called as the Drinfeld-Sokolov-Hirota-Satsuma system. This paper contains two results: (i) a list of integrable systems of the form (1.1) with smooth functions F , G and a = −1/2; (ii) differential substitutions that allow to connect any equation from the list with (1.2) or (1.3).
There are many articles dealing with integrable systems, but some of them (see, e.g., [5,6,7]) consider multi-component systems. Other papers (see, e.g., [8,9,10,11,12]) contain twocomponent systems reducible to a triangular form. The triangular form is briefly considered below. There was possibly only one serious attempt [13] to classify integrable systems of the form (1.1) using the Painlevé test. Unfortunately, fifteen systems presented in [13] contain a large number of constants some of which can be removed by scaling and linear transformations. Note that there are two-field integrable evolutionary systems u t = Au 3 + H(u, u 1 , u 2 ) with a nondiagonal main matrix A. For example, an integrable evolutionary system with the Jordan main matrix is found in [14].
Partial solutions of the classification problem for a = 0 and ord G 1 have been obtained in [20], and in [21] for divergent systems with a = 1. A complete list of integrable systems of the form (1.1) does not exist today because the problem is too cumbersome and the set of integrable systems is very large.
Our tool is the symmetry method presented in many papers. We shall point out pioneer or review papers only. In [22] the notions of formal symmetry and canonical conserved density for a scalar evolution equation are introduced. These tools were applied to classification of the KdVtype equations in [23]. A complete theory of formal symmetries and formal conservation laws for scalar equations has been presented in [24]. A generalized theory was developed for evolutionary systems in [25]. Review paper [26] contains both general theorems of the symmetry method and classification results on integrable equations: the third and fifth order scalar equations, Schrödinger-type systems, Burgers-type equations and systems. Review paper [27] is devoted to higher symmetries, exact integrability and related problems. Peculiarities of systems (1.1) have been discussed in [21]. For the sake of completeness, the main points of the symmetry method and some results necessary for understanding of this paper are considered in the Sections 2-4.
Briefly speaking, the symmetry method deals with the so-called canonical conservation laws D t ρ n = D x θ n , D tρn = D xθn , n = 0, 1, 2, . . . , (1.4) where D t is the evolutionary derivative and D x is the total derivative with respect to x. In particular, ρ 0 = −F u 2 /3,ρ 0 = −G v 2 /(3a). The recursion relations for the canonical conserved densities ρ n andρ n are presented in Section 2. All canonical conserved densities are expressed in terms of functions F and G. That is why equations (1.4) impose great restrictions on the forms of F and G. Equations (1.4) are solvable in the jet space iff E α D t ρ n = 0, E α D tρn = 0, α = 1, 2, n = 0, 1, 2, . . . (1.5) (see [28], for example). Here is the Euler operator. Conservation law with ρ = D x χ, θ = D t χ is called trivial and the conserved density of the form ρ = D x χ is called trivial too. This can be written in the form ρ ∈ Im D x , where Im = Image. If ρ 1 − ρ 2 ∈ Im D x , then the densities ρ 1 and ρ 2 are said to be equivalent.
There are a lot of systems in the following form satisfying the integrability conditions (1.4). Such systems containing one independent equation are said to be triangular. It follows from the integrability conditions that the equation for u must be one of the known integrable equations (KdV, mKdV etc). The second equation is usually linear with respect to v, v x and v xx . Triangular systems do not possess any Lax representations and are not integrable in this sense. Therefore triangular systems and those reducible to the triangular form have been omitted as trivial.
The system of two independent equations will be called disintegrated. It is obvious that the disintegrated form is a partial case of the triangular form. Therefore the disintegrated systems and those reducible to them have been omitted. System (1.1) will be called reducible if it is triangular or can be reduced to triangular or disintegrated form. Otherwise, the system will be called irreducible.
Our computations show that for irreducible integrable systems (1.1) parameter a must belong to the following set: These values were found first in [3] and were repeated in [29]. The value of a is always defined at the end of computations when functions F and G have been found and only some coefficients are to be specified from the fifth or seventh integrability conditions (see example in Section 3.1). This means that it is enough to verify conditions (1.5) for n = 0, . . . , 7 and α = 1, 2 to obtain F , G and a. But for absolute certainty we have verified conditions (1.5) for n = 8, 9 and α = 1, 2 for each system. The presented set A consists of zero and two pairs (a, a −1 ). The transformation t ′ = at, u ′ = v, v ′ = u changes the parameter a = 0 in (1.1) into a −1 . That is why one ought to consider the values 0, − 1 2 , − 7 2 + 3 2 √ 5 of the parameter a. Integrable systems with a = 0 were mentioned above, see also [30]. This paper is devoted to investigation of the case a = −1/2 only. The case a = − 7 2 + 3 2 √ 5 will be presented in another paper. Section 2 contains recursion formulas for the canonical densities. The origin of the notion, some examples and a preliminary classification are considered.
A list of forty integrable systems and an example of computations are presented in Section 3. Section 4 contains differential substitutions that connect all systems from the list. The method of computations and an example are considered. It is shown that all systems from the list presented in Section 3 can be obtained from (1.2) and (1.3) by differential substitutions.
Section 5 is devoted to zero curvature representations. The zero curvature representations for systems (1.2) and (1.3) are obtained from the Drinfeld-Sokolov L-operators. A method of obtaining zero curvature representations for other systems is demonstrated.

Canonical densities
One of the main objects of the symmetry approach to classification of integrable equations is the infinite set of the canonical conserved densities. Let us demonstrate how canonical conserved densities can be obtained the from the asymptotic expansions for eigenfunctions of the Lax operators. The simplest Lax equations concerned with the KdV equation Here u is a solution of the KdV equation and µ is a parameter. The standard substitution ψ = exp ρ dx reduces equations (2.1) and (2.2) to the Riccati form Differentiating temporal equation (2.4) with respect to x one can rewrite it, using (2.3), as the continuity equation: To construct an asymptotic expansion one ought to set Then equation (2.3) results in the following well known recursion formula [2] and (2.5) results in infinite sequence of conservation laws: We change here ∂ t → D t and ∂ x → D x because u is a solution of the KdV equation. The obtained conservation laws are canonical. It is easy to obtain several first canonical densities: It is shown in [2] that all even canonical densities are trivial. Note that if one chooses another asymptotic expansion, for example, in powers of µ −1 instead of (2.6), then another set of canonical densities is obtained, which is equivalent to the previous set. The canonical densities that follow from (2.7) can also be obtained by using the temporal equation (2.4) only. Indeed, setting ∂ t ρ dx = θ one obtains from (2.4) −4(∂ x + ρ) 2 ρ + 6uρ + 3u x + 4µ 3 = θ. (2.9) Using the same expansions as above one can obtain from (2.9) the following recursion relation: where δ i,k is the Kronecker delta. The obtained relation provides ρ 0 = 0, ρ 1 = −u, ρ 2 = −u 1 − θ 0 /3, etc. As D t ρ 0 = D x θ 0 and ρ 0 = 0, then θ 0 = 0. The higher canonical densities ρ n , n > 2 depend on θ n−2 . The fluxes θ n must be defined now from equations (1.4). For example, The traditional method to obtain the canonical densities for an evolution system [25] u t = K(u, u x , . . . , u n ), consists, briefly, in the following. The main idea is to use the linearized equation or its adjoint as the temporal Lax equation. Here The spatial Lax operator (formal symmetry) was introduced in [25] as the infinite operator series commuting with D t − K * . R k are matrix coefficients depending on u, u x , . . . . It was shown that Tr res R (res R = R −1 ) is the conserved density for system (2.10). Canonical densities have been defined by the formulas ρ n = Tr res R n , n = 1, 2, . . . , see [26] for details. Operations with operator series (2.13) are not so simple, therefore we use an alternative method for obtaining the canonical densities. It was proposed in [32] heuristically and we present the following explanation (see also [33]).
Observation. One can obtain equation (2.9) from (2.2) by the following substitution where ρ dx+θ dt is the smooth closed 1-form, that is, D t ρ = D x θ. This implies e −ω D t e ω = D t +θ, e −ω D x e ω = D x + ρ and so (2.9) follows. Another way to obtain the same equation is to prolong the operators D t → ∂ t + θ, D x → ∂ x + ρ in (2.2) formally and to set ψ = 1. For systems, one must set ψ α = 1 for a fixed α only. We shall apply this method to system (1.1) now.
The linearized system (1.1) with prolonged operators D x → D x + ρ, D t → D t + θ takes the following form: If one sets here Ψ 1 = 1, then the first equation takes the following form It is obvious from this equation that the following forms of the asymptotic expansions are acceptable: Here µ is a complex parameter. Then, after some simple calculations, the following recursion relations are obtained (n −1): Here δ i,k is the Kronecker delta, F u 1 = ∂F/∂u 1 and so on. From the recursion relations it is obvious why the value a = 1 is singular. Some of initial elements of the sequence {ρ n , ϕ n } read others are introduced via the δ-symbols. If one sets in (2.15) Ψ 2 = 1 and a = 0, then one more pair of recursion relations for {ρ n ,φ n } is obtained. These recursion relations give us any desired number of canonical densities. As an example, we present here some more canonical densities: The tilde denotes another sequence of canonical densities. Further canonical densities are too cumbersome, therefore we do not present them here.
To simplify investigation of the integrability conditions, an additional requirement is always imposed. This is the existence of a formal conservation law [25,26]. A formal conservation law is an operator series N in powers of D −1 x . An equation for the formal conservation law can be written in the following operator form The form of this equation coincides with the form of the equation for the Noether operator [34].
That is a formal conservation law may be called a formal Noether operator. If (D t − K * , L) is the Lax pair for an equation, then (D t + K * + , L + ) is obviously the Lax pair for the same equation. Hence, canonical densities obtained from (2.11) must be equivalent to canonical densities obtained from (2.12).
It was shown in [21] that the first sequence of the canonical densities ρ n for system (1.1) obtained from (2.11) is equivalent to the first sequence of the canonical densities τ n obtained from (2.12) and the second sequence of the canonical densitiesρ n is equivalent to the second sequence of the canonical densitiesτ n . Hence, ρ n − τ n ∈ Im D x andρ n −τ n ∈ Im D x , or Equations (1.4) (or (1.5)) and (2.18) are said to be the necessary conditions of integrability. We shall refer to it simply as the integrability conditions for brevity. Our computations have shown that Other "adjoint" canonical densities τ i andτ k essentially differ from the "main" canonical densities ρ i andρ k . All canonical densities can be obtained using the Maple routines cd and acd from the package JET (see [36]). These routines generate the "main" and the "adjoint" canonical densities, correspondingly, for almost any evolutionary system (an exclusion is the case of multiple roots of the main matrix of the system under consideration). Thus, according to (2.16) and (2.19) we have F u 2 ∈ Im D x and G v 2 ∈ Im D x (a = 0). This implies the following lemma.
Indeed, one may set From higher integrability conditions one more lemma follows.

Lemma 2. Suppose system (2.20) is irreducible and satisfies the following eight integrability
Then the system must have the following form where ord(f, g, f i , g j ) 1.
A scheme of the proof has been presented in [21].

List of integrable systems
As it is shown in Section 2 the problem of the classification of integrable systems (1.1) is reduced to investigation of system (2.21). That is why it is necessary to start by investigating its symmetry properties.
and under the following permutation transformation where α, β, γ and δ are constants, h i are arbitrary smooth functions.
The classification of systems of type (2.21) has been performed by modulo of the presented transformations.
Moreover, some systems (2.21) admit invertible contact transformations. An effective tool for searching such contact transformations is investigation of the canonical conserved densities. For example, system (3.24) from the next section has the first canonical conserved density of the following form: It is obvious that the best variables for that system are This is an invertible contact transformation. In terms of U and V the system takes the following simple form: If c 1 = 0 this system can be reduced to (3.10) by scaling, otherwise the system is triangular: the equation for V will be independent single mKdV. Moreover, the equation for U becomes linear. That is why c 1 = 0 in (3.24). Canonical densities for the triangular systems contain only one highest order term in the second power as in the considered example ρ = V 2 or ρ = V 2 x + · · · , or ρ = V 2 xx + · · · etc. Triangular systems and those reducible to the triangular form have been omitted in the classification process as trivial.
To classify integrable systems (1.1) with a(a − 1) = 0 one must solve a huge number of large overdetermined partial differential systems for eight unknown functions of four variables. This work has required powerful computers and has taken about six years. All the calculations have been performed in the interactive mode of operation because automatic solving of large systems of partial differential equations is still impossible. The package pdsolve from the excellent system Maple makes errors solving some single partial differential equations. The package diffalg cannot operate with large systems because its algorithms are too cumbersome. Thus, one has to solve complicated problems in the interactive mode. Hence, to obtain a true solution one must enter true data! Under such circumstances errors are probable. The longer the computations the more probable are errors. This is the reason why we cannot state with confidence that all computations have been precise all these six years. That is why the statement on completeness of the obtained set of integrable systems is formulated as a hypothesis.
In this and in the following sections c, c i , k, k i are arbitrary constants.
Hypothesis. Suppose system (2.21) with a = −1/2 is irreducible. If the system has infinitely many canonical conservation laws, then it can be reduced by an appropriate point transformation to one of the following systems: Remark 2. Ten pairs of integrability conditions (for ρ 0 -ρ 9 andρ 0 -ρ 9 ) have been verified for each system (3.1)-(3.40), and nontrivial higher conserved densities with orders 2, 3, 4 and 5 have been found.

Example of computations
Let us consider the simplest case of system (1.1): where a(a − 1) = 0. Formulas (2.16) are reduced now to the following The further canonical densities read where indices after commas denote derivatives.
The first integrability condition (1.5) for ρ 1 can be split with respect to u 3 , v 3 , u 2 and v 2 . This provides the following equations Analogously, condition (1.5) forρ 1 implies It is obvious from (3.43) that the second integrability conditions (2.18) are τ 2 ∈ Im D x andτ 2 ∈ Im D x . These conditions provide f 2,uu = g 1,vv = 0 or f 2 (u, v) = uf 3 (v) + f 4 (v), g 1 (u, v) = vg 3 (u) + g 4 (u). Thus, system (3.41) takes the following form: in an explicit form. The expressions D t ρ 1 and D tρ1 are not the total derivatives yet: where q 1 andq 1 are unknown functions and D x q 1 = R 1 , D xq1 =R 1 . This trick allows us to evaluate ρ 4 , τ 4 ,ρ 4 ,τ 4 and verify the fourth integrability conditions (2.18). These conditions imply f ′′ 3 = g ′′ 3 = 0, hence To simplify the further analysis one must list all irreducible cases of f 1 (or g 2 ). Let us take Lemma 4. Using complex dilatations of u and v, translations u → u + λ 1 , v → v + λ 2 and the Galilei transformation u t → u t + αu x , v t → v t + αv x one can reduce f 1 to one of the following forms: where α is any constant. Moreover, in the cases 4-7 the function g 2 must be linear (b 1 = b 3 = 0) because otherwise the permutation u ↔ v gives one of the cases 1-3.
There are many branches in case 7 but all of them provide linear or triangular systems only. As one can see, classification of integrable systems of the form (3.41) is a sufficiently laborious task. System (2.21) contains eight unknown functions depending on four variables, therefore classification of these systems is much more difficult.

Differential substitutions
A differential substitution is a pair of equations where f and g are some smooth functions.
In all cases that we know, the new systems (Σ) belong to the same class (1.1) with some smooth functions P and Q. There exist some group-theoretical explanation of this fact for KdV type equations [35]. Our attempts to introduce another parameter a ′ = a in (S) had no success. Substituting (4.1) into (1.1) one obtains the following equations It is obvious that transition to the manifold (S) in (4.2) is equivalent to a replacement of ∂ t by the evolutionary differentiation D t performed in accordance with (S): where Another way to obtain (4.3) is to differentiate equations (4.1) with respect to t in accordance with (1.1) and (S) and exclude u and v by using (4.1). This algorithm and many others are coded in Maple, see for example [36].
To find the admissible functions f , g, P , Q from (4.3) one can use the following easily provable formula: and the analogous formula for ∂/∂V k . Differentiating (4.3) with respect to U n+3 and V n+3 , one obtains Other corollaries of (4.3) are too cumbersome to consider them in the general form. Let us consider, as an example, the first order differential substitutions for system (1.2) According to (4.4) one has f = f (U, V, U 1 ), g = g(U, V, V 1 ), hence equations (4.3) now read Differentiating (4.6) with respect to U 3 and V 3 one can obtain four equations: Let us consider some corollaries of these equations.
2. If ∂f /∂U 1 = 0, then u = f (U ) and one can set f (U ) = U by modulo of the point transformation. In this case P = gU 1 from the first of equations (4.6).
3. If ∂g/∂V 1 = 0, then v = g(V ) and one can set g(V ) = V by modulo of the point transformation. In this case Q = f D x f − V V 1 from the second of equations (4.6).
Investigation of cases 2-4 provides seven nontrivial solutions of equations (4.6) (see below (3.1) → (3.2), . . . , (3.1) → (3.17)). Note that integrable system (3.6) admits strange differential substitutions that generate nonintegrable systems. For example, system (3.6) admits the following differential substitution: so that the functions U and V satisfy the following system: with arbitrary function f . This system does not satisfy the integrability conditions (1.4). To comprehend this unusual phenomenon we evaluate V 2 , V 3 and V 4 from the first equation (4.8) and substitute them into the second one. The result is It is easily verified that the obtained constraint is a reduction of system (3.6) into the single KdV equation u t = −1/2u 3 − uu 1 . This means that using a substitution like (4.8) we are trying to construct an integrable system from the single KdV equation. There are other such examples for system (3.6). Note that the reduction obtained above follows from the reduction u = const for system (3. To organize the presented list of systems we have computed admissible differential substitutions for each system and present the results in this section. The formula will denote that if u ′ and v ′ are substituted into system (A), then system (B) follows for u and v.
We say in this case that system (B) is obtained from system (A) by the differential substitution. Substitution (A) → (B) establishes an interrelation between the sets of solutions of systems (A) and (B): (u, v) → (u ′ , v ′ ) is a single valued map. And conversely, if for some solution (u ′ , v ′ ) of system (A) one solves the system of two ordinary differential equations (A) → (B) for u and v, then one or more solutions of system (B) are obtained. Of course, explicit solutions can be obtained very rarely when the substitution is linear or invertible (see below).
Let us consider the following simple example: This substitution is possible for any divergent system This is a well known fact, that is why the substitutions (u, v) → (u 1 , v 1 ) are not written for the divergent systems below. In some cases analogous substitutions are not so obvious and we present them likewise (3.1) → (3.2), for example. The proof can be obtained by a direct verification. List of the substitutions: either. Systems (3.4) and (3.6) admit the substitutions from the first till fourth orders. We do not present higher order substitutions for (3.4) because simpler substitutions exist for (3.6). Fifth and higher order substitutions for systems (3.4) and (3.6) have not been computed because the computations are extremely cumbersome.

Examples of zero curvature representations
The IST method for nonlinear equations with two independent variables is based on investigation of a linear overdetermined system where L and A are ordinary linear operators, u is a smooth (vector) function satisfying a nonlinear partial differential equation E(u) = 0 and λ is a parameter. The operators L, A may be both scalar and matrix. The operators L, A and E may be also pseudodifferential or integrodifferential. The compatibility condition of system (5.1) reads There are two ways to satisfy this condition. The first operator condition was introduced by P.D. Lax [37]: The second more general operator condition was introduced in [38], see also [39]: In all cases, operator L must essentially depend on the parameter λ. This parameter cannot be removed by a gauge transformation L → f −1 Lf with some smooth function f , in particular.
If system (5.1) is differential, then the standard substitution ψ = Ψ 1 , ψ x = Ψ 2 and so on, provides the following first order system: where U and V are square matrices depending on u and λ. The compatibility condition of linear system (5.4) reads if E(u) = 0. Usually a stronger condition is required: (5.5) is valid iff E(u) = 0. In this case equation (5.5) is said to be the zero curvature representation. For an evolutionary system u t = K(u, u x , . . . , u n ), u = {u α } the matrix U usually depend on u only, but it may depend on u, u x , u xx , and so on. Let us consider the general case. If some smooth functions F = F (u, u x , . . . , u r ) and Φ = Φ(u, u x , . . . , u p ) satisfy the condition where the summation over i = 0, . . . , r and α = 1, . . . , m is implied. This implies according to (5.6). Using this result one obtains the following identity: for any F and Φ satisfying (5.6). Let us apply identity (5.7) to equation (5.5). If the matrix U depends on u only then It is obvious now that equation (5.5) is equivalent to u t = K iff the matrices ∂U/∂u α , α = 1, . . . , m are linearly independent. Suppose now that the matrix U depends on u and u x then one obtains If the matrices ∂U/∂u α , ∂U/∂u β x , α, β = 1, . . . , m are linearly independent, then equation (5.5) is equivalent to u t = K again. Otherwise, equation (5.5) would be equivalent to some differential consequence of the system u t = K that is a more general system than the original one. In this case we call equation (5.5) the generalized zero curvature representation.
It is well known that equations (5.4) and (5.5) are covariant under the following transforma-tionΨ where S is any non-degenerate matrix. This transformation is called a gauge one. Any gauge transformation is invertible and preserves compatibility of system (5.4). Two (L, A)-pairs were proposed for system (1.2) in [3]. One of these (L, A)-pairs coincides with the (L, A)-pair that was presented in [4]. The L-operator of the common (L, A)-pair takes the form L = (∂ 2 The temporal Lax equation reads ψ t = Aψ, where A is a fractional degree of L. The spatial Lax equation Lψ = λ 2 ψ can be transformed into the system (∂ 2 x − g)ψ = λϕ, (∂ 2 x + f )ϕ = λψ and then into the normal form (5.4), were Here f 1 and f 2 take the following form: Matrices (5.11) realize the zero curvature representation of system (1.2). System (1.3) also has two Lax representations (see [3]). Using the simpler L-operator, we have found, similarly to the previous case, the following matrices that realize the zero curvature representation of system (1.3): where Let us consider an admissible differential substitution u = f (ũ i ,ṽ j ), v = g(ũ i ,ṽ j ) of system (1.1). Substituting u and v in the matrices U (u, v) and V (u i , v j ) one obtainŝ AsÛ depends on derivatives ofũ orṽ, then one has a generalized zero curvature representation.
To obtain an ordinary zero curvature representation one can try to remove higher order derivatives from the matrixÛ using the gauge transformation (5.10). But this is not always possible (see example B below).
A. Performing the substitution (3.1) → (3.9) into the matrix U from (5.11) one obtainŝ One can easily verify that matrices A =Û u 2 and B =Û v 2 are commutative, hence the system S x = (Au xx + Bv xx )S has the following solution S = exp(Au x + Bv x ). The matrix U 1 =Ū evaluated according to (5.10) takes the following form Now another gauge transformation is possible with the following diagonal matrix: where diag (U 1 ) is the main diagonal of U 1 . This gauge transformation provides the following matrix A corresponding V -matrix can be obtained by solving equation (5.5) directly or by the previous twofold gauge transformation. This matrix takes the following form The matrices U 2 and V 2 realize the zero curvature representation of system (3.9). B. Substitution (3.1) → (3.17) reduces matrix U from (5.11) to the following form where R = 2(u 1 + v 2 )/6. It is obvious that one cannot remove u 1 fromÛ by a gauge transformation. It is clear from the structure of the matrixÛ that where A and B are some linearly independent matrices. Thus, this zero curvature representation for system (3.17) is generalized. Of course, one may introduce here the new variable u ′ = u 1 to obtain an ordinary zero curvature representation. But we do not know if it is always possible. C. Performing the substitution (3.3) → (3.6) into the matrix U from (5.12) one obtainŝ The first gauge transformation is performed using S 1 = exp(u 1 (∂Û/∂u 2 )): The transformed U -matrix is The second gauge transformation is performed using S 2 = exp(u (∂Ũ 1 /∂u 1 )): MatricesŨ 2 andṼ 2 realize the zero curvature representation of system (3.6).

Conclusion
The examples in Section 5 illustrate the fact that some systems possess ordinary zero curvature representation while others possess generalized zero curvature representation. All these representations are obtained from the Drinfeld-Sokolov L, A operators by using corresponding differential substitutions listed in Section 4. Matrices U and V that realize all zero curvature representations have the size 4 × 4. Thus, the two-field evolutionary systems presented above are integrable in principle by the inverse spectral transform method. But the fact is that the inverse scattering problem for differential equations with order more than two is extremely difficult. That is why other methods for solution of equations may be useful [40]. They may be Bäcklund transformations [41], Darboux transformations [42,43], Hirota method [44] or numeric simulating (see [45], for example).