# Applications of Nijenhuis Geometry IV: multi-component KdV and Camassa-Holm
equations

Alexey V. Bolsinov111 School of Mathematics, Loughborough University, LE11
3TU, UK<EMAIL_ADDRESS>& Andrey Yu. Konyaev222Faculty of Mechanics and
Mathematics, Moscow State University, 119992, Moscow Russia<EMAIL_ADDRESS>& Vladimir S. Matveev333 Institut für Mathematik, Friedrich Schiller
Universität Jena, 07737 Jena Germany<EMAIL_ADDRESS>
###### Abstract

We construct a new series of multi-component integrable PDE systems that
contains as particular examples (with appropriately chosen parameters) and
generalises many famous integrable systems including KdV, coupled KdV [1],
Harry Dym, coupled Harry Dym [2], Camassa-Holm, multi-component Camassa-Holm
[14], Dullin-Gottwald-Holm and Kaup-Boussinesq systems. The series also
contains integrable systems with no low-component analogues.

###### Contents

  1. 1 Introduction
  2. 2 Explicit formulas for new integrable systems, their symmetries and conservation laws
    1. 2.1 Four types of equations
    2. 2.2 Commuting flows and conservation laws for the equations of Types I – IV
  3. 3 General construction
    1. 3.1 Parameters of the general construction
    2. 3.2 Main theorems
    3. 3.3 Parameters corresponding to known integrable systems
  4. 4 Proofs of Theorems 2 and 3
  5. 5 Proof of Theorem 4
  6. 6 Conclusion

##  1 Introduction

In our paper we consider two types of PDE systems. The first one is an
$n$-component evolutionary system of PDEs:

$u^{i}_{t}=\xi^{i}(u),\quad i=1,\dots,n,$ (1)

where $u^{i}=u^{i}(x,t)$ are unknown functions and $\xi^{i}(u)$ is a
differential polynomial in $u^{1},\dots,u^{n}$, that is a polynomial in
derivatives $u^{j}_{x},u^{j}_{xx},u^{j}_{xxx},\dots$, $j=1,\dots,n$, whose
coefficients are functions of $u^{1},\dots,u^{n}$. For $n=1$, an example of
such a system is the famous KdV equation

$u_{t}=\frac{1}{2}u_{xxx}+\frac{3}{2}uu_{x}.$

A well known two-component case is the Kaup-Boussinesq system, see e.g. [22,
eqn. (4)]:

$\displaystyle u^{1}_{t}$
$\displaystyle=u^{2}_{x}-\frac{3}{2}u^{1}u^{1}_{x}\,,$ $\displaystyle
u^{2}_{t}$
$\displaystyle=\frac{m}{2}u^{1}_{xxx}-u^{2}u^{1}_{x}-\frac{1}{2}u^{1}u^{2}_{x}\,.$

where $m$ is an arbitrary constant.

We will also consider evolutionary systems of PDEs with differential
constraints. They are given by

$\displaystyle u^{i}_{t}$ $\displaystyle=\xi^{i}(u,q),\quad i=1,\dots,n,$ (2)
$\displaystyle 0$ $\displaystyle=p(u,q).$

Here, unknown functions are $u(x,t)=(u^{1}(x,t),\dots,u^{n}(x,t))^{\top}$ and
$q(x,t)$. Each $\xi^{i}(u,q)$ is a differential polynomial in
$u^{1},...,u^{n}$ and $q$, whereas the differential constraint  $p(u,q)$ is a
differential polynomial in $q$ and a function in $u$ (but not in derivatives
of $u$). In most cases we consider just one constraint; the only exception is
Example 3.5.

In all the systems that we construct and analyse, the differential constraint
$p$ can be explicitly resolved w.r.t. at least one of $u^{i}$’s. Solving
$p(u,q)=0$ with respect to $u^{i}$ and substituting
$u^{i}=\tilde{p}(u^{1},\dots,u^{i-1},u^{i+1},\dots,u^{n},q)$ in the first
equation of (2) allows us to think of (2) as a system of $n-1$ evolutionary
PDEs and one non-evolutionary PDE on $n$ unknown functions
$u^{1},\dots,u^{i-1},u^{i+1},\dots,u^{n}$ and $q$. As an example, consider the
Dullin-Gottwald-Holm equation [11] (in this case, $n=1$):

$\displaystyle u_{t}$
$\displaystyle=\frac{\gamma}{2}q_{xxx}+uq_{x}+\frac{1}{2}qu_{x},$ (3)
$\displaystyle u$ $\displaystyle=q+\frac{m}{2}q_{xx}.$

Substituting the expression for $u$ into the first equation gives

$q_{t}+\frac{m}{2}q_{xxt}=\frac{\gamma}{2}q_{xxx}+\frac{3}{2}qq_{x}+\frac{m}{2}q_{xx}q_{x}+\frac{1}{2}qq_{xxx},$

which is a single non-evolutionary PDE on $q(x,t)$.

We will also deal with PDE systems given by (1) with the right hand side
$\xi^{i}$ being a formal differential series (i.e., infinite sum of monomials
in derivative variables $u^{j}_{x},u^{j}_{xx},u^{j}_{xxx},\dots$ with
coefficients being functions of $u^{1},...,u^{n}$). We refer to such systems
as formal evolutionary PDEs (systems of type (1) with $\xi^{i}$ being a
differential polynomial will be called non-formal).

Let us explain, using the Dullin-Gottwald-Holm equation (3) as an example, the
relation between systems of type (2) and formal evolutionary PDEs. Observe
that the second equation of (3) can be formally solved with respect to $q$:

$q=u-\frac{m}{2}u_{xx}+\frac{m^{2}}{4}u_{xxxx}-\dots$

(we discuss neither convergence of this differential series nor boundary or
other conditions). Substituting this formal expression into the first equation
of (3) leads to the formal evolutionary PDE

$u_{t}=\frac{3}{2}uu_{x}+\frac{\gamma}{2}u_{xxx}-\frac{m}{2}\Big{(}\frac{\gamma}{2}u_{xxxxx}-uu_{xxx}-\frac{1}{2}u_{xx}u_{x}\Big{)}+\dots,$

depending on $m$ as a parameter.

Recall that a (formal) evolutionary PDE-system $u^{i}_{\tau}=\eta^{i}(u)$ is a
(formal) symmetry of a (formal) evolutionary PDE-system $u^{i}_{t}=\xi^{i}(u)$
(or, equivalently, these two PDE-systems commute) if the following commutator
vanishes

$\sum\limits_{j=0}^{\infty}\sum_{\alpha=1}^{n}\Big{(}\frac{\partial\xi^{i}}{\partial
u^{\alpha}_{x^{j}}}D^{j}\big{(}\eta^{\alpha})-\frac{\partial\eta^{i}}{\partial
u^{\alpha}_{x^{j}}}D^{j}\big{(}\xi^{\alpha})\Big{)}=0.$ (4)

Here and below $u_{x^{j}}=u_{\underset{j\ \mathrm{times}}{xx...x}}$,
$D=\tfrac{d}{dx}$ is the total derivative in $x$ (for example,
$D(u^{\alpha}u^{\alpha}_{x})=(u^{\alpha}_{x})^{2}+u^{\alpha}u^{\alpha}_{xx}$)
and $D^{j}$ stands for the $j$-th power of $D$ (for example
$D^{3}(u^{\alpha})=u^{\alpha}_{xxx}=u^{\alpha}_{x^{3}}$). In the case of usual
(= non-formal) PDEs, $u^{i}_{\tau}=\eta^{i}(u)$ is a symmetry of
$u^{i}_{t}=\xi^{i}(u)$ if an only if (at least in the analytic category) both
systems can be solved simultaneously, that is, there exist $n$ functions
$u^{i}(x,t,\tau)$, which satisfy all $2n$ PDEs $u^{i}_{\tau}=\eta^{i}(u)$ and
$u^{i}_{t}=\xi^{i}(u)$ .

A (formal) differential series $\mathrm{v}$ is said to be a conservation law
density of the (formal) evolutionary equation $u^{i}_{t}=\xi^{i}(u)$ if

$\mathrm{v}_{t}=\sum\limits_{j=0}^{\infty}\sum_{\alpha=1}^{n}\frac{\partial\mathrm{v}}{\partial
u^{\alpha}_{x^{j}}}D^{j}\big{(}\xi^{\alpha}\big{)}=D\mathrm{w}$ (5)

for some formal differential series $\mathrm{w}$. Such $\mathrm{v}$ is defined
up to addition of an arbitrary total derivative $D\tilde{\mathrm{v}}$. In the
non-formal case, $\mathrm{v}$ should be understood as the density of a
conservation law: in this case, for any fast decaying or periodic solution
$u(t,x)$, the integral $\int\mathrm{v}(u,u_{x},\dots)\mathrm{d}\,x$ is
independent of $t$ for any solution $u(x,t)$ of (1). We notice that both
formulas (4) and (5) ‘respect’ the degree of differential monomials so that
they are well defined for formal differential series.

The main result of our paper is a construction of a family of multi-component
integrable PDE systems of the form (1) and (2). There are several different
notions of integrability in this context in the literature. In the case of
systems (1), we construct infinite hierarchies of (non-formal) conservation
laws and of (non-formal) pairwise commuting symmetries. The differential
degrees of the conservations laws and symmetries grow within each hierarchy.

By integrability of systems (2), we understand the existence of infinitely
many independent (possibly, formal) commuting symmetries and conservations
laws of the formal evolutionary equation obtained from (2) by the method
demonstrated in the example above; this can be done for all systems we
construct.

The paper is organised as follows. In Section 2 we describe our main examples.
They are parameterised by certain discrete and continuous parameters and are
of four types. The equations of Type II and Type IV are evolutionary, whereas
those of Types I and III are evolutionary with differential constrains. All of
them, however, can be considered as different representatives of one single
integrable system understood as a commutative algebra of evolutionary flows,
which includes both formal and non-formal flows. We give explicit recursive
formulas for common (non-formal) symmetries and conservation laws for all
these PDE systems, see Theorem 1.

All the famous integrable systems listed in the abstract correspond to certain
choices of parameters. In Section 3.3, we specify those parameters for each of
them. We also obtain other known integrable systems, e.g., the Marvan-Pavlov
system [20], which we essentially generalise. The generalisations are given by
explicit formulas and have no low-component analogues.

Theorem 1 follows directly from a more general construction described in
Section 3, see Theorems 2, 3 and 4. Namely, Theorem 2 constructs a family of
evolutionary PDEs with differential constrains, as well as formal commuting
symmetries and conservations laws for them. Next, Theorem 3 explains how one
can ‘cook’ non-formal evolutionary PDEs, non-formal symmetries and non-formal
conservation laws starting from those constructed in Theorem 2. The
constructions in Theorems 2, 3 depend on a solution of a certain, possibly
infinite, system of PDEs on a Nijenhuis manifold $(\mathrm{M}^{n},L)$. In
Theorem 4, we solve this system under the additional assumption that $L$ is
differentially non-degenerate and, hence, come to the integrable systems from
Theorem 1.

Let us now comment on the circle of ideas which led us to these results. The
construction was developed within the Nijenhuis Geomery programme [8]. Its
main ingredient is a Nijenhuis operator $L$ on a manifold $\mathrm{M}^{n}$,
that is a (1,1)-tensor field $L=L^{i}_{j}$ such that its Nijenhuis torsion
vanishes, i.e.,

$L^{2}[\nu,\eta]-L[L\nu,\eta]-L[\nu,L\eta]+[L\nu,L\eta]=0$

for arbitrary vector fields $\nu,\eta$. In our recent paper [7], we
constructed all non-degenerate pencils of compatible $\infty$-dimensional
Poisson structures of type $\mathcal{P}_{3}+\mathcal{P}_{1}$, where the
Poisson structure $\mathcal{P}_{1}$ has order 1 and $\mathcal{P}_{3}$ is a
Darboux-Poisson structure of order 3. Magri-Lenard scheme applied to these
pencils leads to certain integrable bi-Hamiltonian systems. Translating them
to the language of Nijenhuis Geometry allowed us to generalise our
construction further and obtain integrable systems which are not necessarily
Hamiltonian. We view Nijenhuis Geometry as the most natural framework for them
and expect that the systems and their properties can and should be understood
in the context of Nijenhuis operators, with no other geometric structure
involved.

##  2 Explicit formulas for new integrable systems, their symmetries and
conservation laws

###  2.1 Four types of equations

A Nijenhuis operator $L$ on $\mathrm{M}^{n}$ is called differentially non-
degenerate, if the differentials of the coefficients of its characteristic
polynomial are linearly independent at every point. Typical examples of
differentially non-degenerate Nijenhuis operators are as follows:

$L_{\mathsf{comp}}=\left(\begin{array}[]{ccccc}u^{1}&1&0&\dots&0\\\
u^{2}&0&1&\dots&0\\\ \vdots&&&\ddots&\\\ u^{n-1}&0&0&\dots&1\\\
u^{n}&0&0&\dots&0\end{array}\right)\ \ \textrm{ and }\ \
L_{\mathsf{diag}}=\left(\begin{array}[]{ccccc}x^{1}&0&0&\dots&0\\\
0&x^{2}&0&\dots&0\\\ \vdots&&\ddots&&\vdots\\\ 0&\dots&\dots&x^{n-1}&0\\\
0&0&\dots&\dots&x^{n}\end{array}\right).$ (6)

Here $u^{1},...,u^{n}$ and $x^{1},...,x^{n}$ are local coordinate charts on
$\mathrm{M}^{n}$. Moreover, in the case of $L_{\mathsf{diag}}$ we assume that
$x^{i}$ are all different. In fact, these two operators are locally
isomorphic: if we rewrite $L_{\mathsf{diag}}$ in the coordinates
$u^{1}(x),...,u^{n}(x)$ that are coefficients of its characteristic polynomial
(more precisely, we set
$\det(t\operatorname{Id}-L_{\mathsf{diag}})=t^{n}-u^{1}t^{n-1}-\dots-u^{n}$),
then it transforms into $L_{\mathsf{comp}}$. Note also that every
differentially non-degenerate Nijenhuis operator $L$ reduces to the companion
form $L_{\mathsf{comp}}$ by an appropriate coordinate transform, and to the
diagonal form $L_{\mathsf{diag}}$ near those points where $L$ has $n$ distinct
real eigenvalues, see e.g. [8].

Choose $N\geq 0$ and consider natural numbers $n_{0},n_{1},\dots,n_{N}$ and
$\ell_{1},...,\ell_{N}$ with conditions
$n_{0}+n_{1}+\dots+n_{N}=n=\dim\mathrm{M}$ and
$n_{0}-\ell_{1}n_{1}-\dots-\ell_{N}n_{N}=d\geq 0$, and fix a polynomial

$m(\lambda)=m_{0}+m_{1}\lambda+\dots+m_{d}\lambda^{d}$ (7)

of degree $\leq d$. These are the parameters of our construction: $2N+1$
natural numbers and $d+1$ coefficients $m_{0},\dots,m_{d}$.

Next, take the direct product $\mathrm{M}^{n}=U_{0}\times\dots\times U_{N}$ of
$N+1$ discs $U_{0},...,U_{N}$ of dimensions $n_{0},n_{1},...,n_{N}$ equipped
with differentially non-degerenerate Nijenhuis operators $L_{0},...,L_{N}$.
The operator field $L$ on $\mathrm{M}^{n}$ is defined as

$L=L_{0}\oplus\dots\oplus L_{N}.$ (8)

Consider the following family of functions $\sigma(\lambda)$ on ${M^{n}}$,
depending on $\lambda$ as a parameter (in general $\lambda\in\mathbb{C}$, so
the functions might be complex-valued):

$\sigma(\lambda)=\frac{\operatorname{det}(L_{0}-\lambda\operatorname{Id})}{\operatorname{det}(L_{1}-\lambda\operatorname{Id})^{\ell_{1}}\dots\operatorname{det}(L_{N}-\lambda\operatorname{Id})^{\ell_{N}}}$
(9)

Here in each expression
$\operatorname{det}(L_{i}-\lambda\operatorname{Id})^{\ell_{i}}$, the identity
matrix $\operatorname{Id}$ is of the same size as $L_{i}$, i.e., $n_{i}\times
n_{i}$ and $\ell_{i}$ denotes the power.

Next, consider the vector field $\zeta_{0}$ on $\mathrm{M}^{n}$ uniquely
defined by the relations

$\mathcal{L}_{\zeta_{0}}\operatorname{det}(L_{0}-\lambda\operatorname{Id})=1,\quad\mathcal{L}_{\zeta_{0}}\operatorname{det}(L_{i}-\lambda\operatorname{Id})=0,\quad
i=1,\dots,N,$

where $\mathcal{L}_{\zeta_{0}}$ denotes the Lie derivative, and define another
vector field $\zeta$ by setting

$\zeta=p(L)\zeta_{0}.$ (10)

where
$p(t)=\operatorname{det}(L_{1}-\lambda\operatorname{Id})^{\ell_{1}}\dots\operatorname{det}(L_{N}-\lambda\operatorname{Id})^{\ell_{N}}m(t)$
is a polynomial in $t$ with coefficients being functions on $\mathrm{M}^{n}$.

Based on these settings, we finally introduce four types of equations. As
unknown functions, we consider $u^{1}(x,t),\dots,u^{n}(x,t),q(x,t)$, where
$(u^{1},\dots,u^{n})$ should be understood as coordinates on $\mathrm{M}^{n}$
and $q$ as an additional function.

Type I.  For each real or complex number $\lambda$ consider the equations

$\displaystyle u_{t}$
$\displaystyle=q_{xxx}(L-\lambda\operatorname{Id})^{-1}\zeta+q(L-\lambda\operatorname{Id})^{-1}u_{x},$
(11) $\displaystyle 0$
$\displaystyle=m(\lambda)q_{xx}q-\frac{1}{2}m(\lambda)(q_{x})^{2}+\sigma(\lambda)q^{2}-1.$

This is a system of the form (2), i.e., an $n$-component system with a
differential constraint.

Type II. If $\lambda_{i}$ is a root of $m(\lambda)$, then the differential
constraint in Type I becomes degenerate and takes the form
$\sigma(\lambda_{i})q^{2}=1$. Thus, we get an evolutionary PDE

$\displaystyle u_{t}$
$\displaystyle=\Bigg{(}\frac{1}{\sqrt{\sigma(\lambda_{i})}}\Bigg{)}_{xxx}(L-\lambda_{i}\operatorname{Id})^{-1}\zeta+\Bigg{(}\frac{1}{\sqrt{\sigma(\lambda_{i})}}\Bigg{)}(L-\lambda_{i}\operatorname{Id})^{-1}u_{x}.$
(12)

Type III.  Consider the equation

$\displaystyle u_{t}$
$\displaystyle=q_{xxx}\zeta+(L+q\operatorname{Id})\,u_{x},$ (13)
$\displaystyle\frac{1}{2}\Big{(}\operatorname{tr}L_{0}$
$\displaystyle-\sum_{j=1}^{N}\ell_{j}\operatorname{tr}L_{j}\Big{)}=q+(-1)^{d}\,\frac{m_{d}}{2}\,q_{xx},$

where $m_{d}$ is the highest coefficients of the polynomial
$m(t)=m_{0}+m_{1}t+\dots+m_{d}t^{d}$. This is, again, a system of of the form
(2), i.e., an $n$-component system with a differential constraint.

Type IV.  Assume now $m_{d}=0$ (in this case we say that $m(t)$ has a root at
infinity, the terminology will be clarified later). Then (13) takes the form

$\displaystyle
u_{t}=\frac{1}{2}\Big{(}\operatorname{tr}L_{0}-\sum_{j=1}^{N}\ell_{j}\operatorname{tr}L_{j}\Big{)}_{xxx}\zeta+\Big{(}L+\frac{1}{2}\Big{(}\operatorname{tr}L_{0}-\sum_{j=1}^{N}\ell_{j}\operatorname{tr}L_{j}\Big{)}\operatorname{Id}\Big{)}u_{x}.$
(14)

###  2.2 Commuting flows and conservation laws for the equations of Types I –
IV

Here we describe an explicit procedure that generates commuting symmetries and
conservation laws for the above four types of equations.

Step 1. In the one-component case, consider the relation

$\sigma=\frac{1}{2}\mathrm{u}^{2}+\mathrm{u}_{x}.$ (15)

and its formal solution $\mathrm{u}=\mathcal{u}_{1}+\mathcal{u}_{2}+\dots$ as
a differential series in $\sigma$.

The recursion formula for the components of $\mathrm{u}$ from (15) was
essentially discovered by Kruskal and Miura. In the form we need (up to
notation) it appeared e.g. in [23, eqns. 2.16–2.19]:

$\mathcal{u}_{1}=\sqrt{2\sigma},\quad\mathcal{u}_{i+1}=-\frac{1}{\mathcal{u}_{1}}\Bigg{(}\,\frac{1}{2}\sum\limits_{j=2}^{i}\mathcal{u}_{j}\mathcal{u}_{i+2-j}+(\mathcal{u}_{i})_{x}\Bigg{)},\quad
i\geq 1.$ (16)

In the expansion $\mathrm{u}=\sum\mathcal{u}_{i}$, we are interested in the
odd terms only and introduce two formal differential series

${\mathrm{v}}(\sigma,m)=\sqrt{2}\sum_{s=0}^{\infty}(-1)^{s}m^{s}\mathcal{u}_{2s+1}$
(17)

and

$\mathrm{w}(\sigma,m)=\sqrt{2}\sum_{s=0}^{\infty}(-1)^{s}m^{s}\delta\mathcal{u}_{2s+1}.$
(18)

Here $\delta$ stands for the variational derivative w.r.t. $\sigma$ and $m$ is
considered as a formal parameter.

Important fact. The series ${\mathrm{w}}(\sigma,m)$ satisfies the following
identity

$m{\mathrm{w}}_{xxx}(\sigma,m)+2\sigma{\mathrm{w}}_{x}(\sigma,m)+\sigma_{x}{\mathrm{w}}(\sigma,m)=0.$
(19)

This identity is essentially known and comes from the theory of local
infinite-dimensional Poisson structures. It can be understood as the fact that
${\mathrm{v}}(\sigma,m)$ is a formal Casimir of the Poisson structure defined
by the operator $mD^{3}+2\sigma D+\sigma_{x}$ (for details and proof see e.g.
[23, Proposition 2.1] and also [16, Theorem 2.4] for $n$-component systems).
One can also view (19) as a way of applying the Magri-Lenard scheme to produce
commuting symmetries for the Harry Dym equation. However, we do not need such
an interpretation and will use identity (19) as it is.

We will also use the following crucial observation by Gelfand and Dikii [13]
(see also [3]). Multiplying the l.h.s. of (19) by $\mathrm{w}(\sigma,m)$ and
integrating in $x$ we get another important identity for $\mathrm{w}$ (cf. the
differential constraint from (11)):

$m\,\mathrm{w}_{xx}\mathrm{w}-\frac{1}{2}\,m\,\mathrm{w}_{x}^{2}+\sigma\,\mathrm{w}^{2}=1,$
(20)

which, in particular, allows us to reconstruct all the terms of (18) step by
step starting from the first term
$\sqrt{2}\,\delta\mathcal{u}_{1}=\frac{1}{\sqrt{\sigma}}$.

Step 2. For $m(\lambda)$ and $\sigma(\lambda)$ defined by (7) and (9)
respectively, consider the formal differential series (in any local
coordinates $u^{1},\dots,u^{n}$ on $\mathrm{M}^{n}$) with components depending
on parameter $\lambda$

$\mathrm{v}(\lambda)=\mathrm{v}(\sigma(\lambda;u),m(\lambda)),\quad\mathrm{w}(\lambda)=\mathrm{w}(\sigma(\lambda;u),m(\lambda)),$
(21)

obtained by replacing $\sigma$ with
$\sigma(\lambda)=\sigma(\lambda;u^{1},\dots,u^{n})$ and $m$ with $m(\lambda)$
in (17) and (18). Next, introduce the formal $n$-component vector field

$\xi(\lambda)=\mathrm{w}(\lambda)_{xxx}\,(L-\lambda\operatorname{Id})^{-1}\zeta+\mathrm{w}(\lambda)(L-\lambda\operatorname{Id})^{-1}u_{x}.$
(22)

This vector field is naturally related to the r.h.s. of the first equation of
system (11) of Type 1. In fact, (11) can be equivalently rewritten as
$u_{t}=\xi(\lambda)$. It follows from the fact that the differential series
$\mathrm{w}(\lambda)$ satisfies the same differential relation (20) as the
function $q$ (see second equation of (11)) and can be uniquely reconstructed
from it (see ‘Important fact’ in Step 1).

Step 3. For each root $\lambda_{i}$ of the polynomial $m(\lambda)$, expand
both $m(\lambda)$ and $\sigma(\lambda,u)$ into Taylor series in powers of
$\varepsilon=\lambda-\lambda_{i}$, i.e.,

$m(\lambda_{i}+\varepsilon)=\sum_{s=1}^{\infty}m_{s,\lambda_{i}}\varepsilon^{s},\quad\sigma(\lambda_{i}+\varepsilon,u)=\sum_{s=0}^{\infty}\sigma_{s,\lambda_{i}}(u)\varepsilon^{s}$

Substitution $m(\lambda)=\sum_{s=1}^{\infty}m_{s,\lambda_{i}}\varepsilon^{s}$
and
$\sigma(\lambda;u)=\sum_{s=0}^{\infty}\sigma_{s,\lambda_{i}}(u)\varepsilon^{s}$
“transforms” $\mathrm{v}(\lambda)=\mathrm{v}(\sigma(\lambda;u),m(\lambda))$
and $\mathrm{w}(\lambda)=\mathrm{v}(\sigma(\lambda,u),m(\lambda))$, as well as
$\xi(\lambda)$ defined by (22), into series in powers of $\varepsilon$ :

$\mathrm{v}\bigl{(}\sigma(\lambda_{i}+\varepsilon,u),m(\lambda_{i}+\varepsilon)\bigr{)}=\sum_{s=0}^{\infty}\mathcal{v}_{s,\lambda_{i}}\varepsilon^{s},\quad\mbox{with}\
\ \mathcal{v}_{0,\lambda_{i}}=2\sqrt{\sigma(\lambda_{i})}.$ (23)
$\mathrm{w}(\sigma(\lambda_{i}+\varepsilon),m(\lambda_{i}+\varepsilon))=\sum_{s=0}^{\infty}\mathcal{w}_{s,\lambda_{i}}\varepsilon^{s},\quad\mbox{with}\
\ \mathcal{w}_{0,\lambda_{i}}=\frac{1}{\sqrt{\sigma(\lambda_{i})}}.$ (24)
$\xi=\sum_{s=0}^{\infty}\xi_{s,\lambda_{i}}\varepsilon^{s},\mbox{with}\ \
\xi_{0,\lambda_{i}}=\left(\tfrac{1}{\sqrt{\sigma(\lambda_{i})}}\right)_{xxx}(L-\lambda\operatorname{Id})^{-1}\zeta+\tfrac{1}{\sqrt{\sigma(\lambda_{i})}}\,(L-\lambda\operatorname{Id})^{-1}u_{x}.$
(25)

Similarly for $\lambda=\infty$, in the above construction we replace
$\sigma(\lambda)$ and $m(\lambda)$ with
$\bar{\sigma}(\lambda)=(-\lambda)^{d}\sigma\left(\frac{1}{\lambda}\right)$,
$\bar{m}(\lambda)=(-\lambda)^{d}m\left(\frac{1}{\lambda}\right)$ and also
$L-\lambda\operatorname{Id}$ with $\operatorname{Id}-\lambda L$. In
particular, we set

$\displaystyle\bar{\mathrm{v}}(\lambda)$
$\displaystyle=\mathrm{v}(\bar{\sigma}(\lambda),\bar{m}(\lambda)),\quad\bar{\mathrm{w}}(\lambda)=\mathrm{w}(\bar{\sigma}(\lambda),\bar{m}(\lambda)),$
$\displaystyle\bar{\xi}(\lambda)$
$\displaystyle=\bar{\mathrm{w}}(\lambda)_{xxx}\,(\operatorname{Id}-\lambda
L)^{-1}\zeta+\bar{\mathrm{w}}(\lambda)(\operatorname{Id}-\lambda
L)^{-1}u_{x}.$

Then if $\bar{m}(0)=0$, we substitute $\bar{\sigma}(\varepsilon)$,
$\bar{m}(\varepsilon)$ into these relations and expand in powers of
$\varepsilon$ to get444In these power series, we shift indices of coefficients
by 1. The reason is that the first terms of these expansions are trivial and
can be ignored. This shift also allows us to keep notation consistent with the
case of $\lambda_{i}\neq\infty$.

$\bar{\mathrm{v}}(\varepsilon)=\sum_{s=0}^{\infty}\mathcal{v}_{s-1,\infty}\varepsilon^{s},\quad\mbox{with}\
\ \mathcal{v}_{-1,\infty}=2\sqrt{\bar{\sigma}(0)}=2\mbox{ and
}\mathcal{v}_{0,\infty}=-2f(u);$ (26)
$\bar{\mathrm{w}}(\varepsilon)=\sum_{s=0}^{\infty}\mathcal{w}_{s-1,\infty}\varepsilon^{s},\quad\mbox{with}\
\ \mathcal{w}_{-1,\infty}=1\mbox{ and
}\mathcal{w}_{0,\infty}=\frac{1}{2}f(u);$ (27)
$\bar{\xi}(\varepsilon)=\sum_{s=0}^{\infty}\xi_{s-1,\infty}\varepsilon^{s},\quad\mbox{with
$\xi_{-1,\infty}=u_{x}$ and
$\xi_{0,\infty}=\frac{1}{2}f_{xxx}\zeta+(L+\frac{1}{2}f\operatorname{Id})u_{x}$}$
(28)

where $f$ is defined from $\bar{\sigma}(\varepsilon)=1-\varepsilon f+\dots$,
i.e.,
$f=\frac{d}{d\varepsilon}|_{\varepsilon=0}\bar{\sigma}(\varepsilon)=\operatorname{tr}L_{0}-\sum_{j=1}^{N}\ell_{j}\operatorname{tr}L_{j}$.

The coefficients of the $\varepsilon$-expansions (23), (25), (26) and (28)
define hierarchies of common commuting symmetries and conservation laws for
the above introduced equations of types I – IV. Namely, we have

######  Theorem 1.

Let $\lambda_{1},\dots,\lambda_{r}$ be the roots of the polynomial
$m(\lambda)$ (including $\infty$ when appropriate). Then the evolutionary PDE
systems

$u_{t}=\xi_{s,\lambda_{i}},\quad i=1,\dots,r,\ \ s=0,1,\dots.$

are commuting symmetries and the differential polynomials

$\mathcal{v}_{s,\lambda_{i}},\quad i=1,\dots,r,\ \ s=0,1,\dots$

are conservation law densities for equations (11)–(14) of Types I – IV.
Moreover, the equations (12, Type II) and (14, Type IV) take the form
$u_{t}=\xi_{0,\lambda_{i}}$ for $\lambda_{i}\neq\infty$ and
$\lambda_{i}=\infty$ respectively.

Thus, for a Nijenhuis operator $L$ (decomposed into differentially non-
degenerate blocks), Theorem 1 gives a series of multi-component integrable
systems and provide, for each of them, commuting symmetries and conservation
laws that can be constructed by an explicit iterative procedure.

##  3 General construction

###  3.1 Parameters of the general construction

Let $L$ be a Nijenhuis operator on $\mathrm{M}^{n}$ and
$f:\mathrm{M}^{n}\to\mathbb{R}$ be a function such that the 1-form
$L^{*}\mathrm{d}\,f$ is closed so that locally
$L^{*}\mathrm{d}\,f=\mathrm{d}\,f_{1}$ for some function $f_{1}$. Then (see
Section 3 in [19]) there exists an infinite sequence of functions $f_{k}$,
$k=1,\dots$, such that $\mathrm{d}\,f_{k}=(L^{*})^{k}\mathrm{d}\,f$. We refer
to $f$ as a conservation law for the Nijenhuis operator $L$. The name comes
from the fact that $f$ provides a conservation law in the sense of (5) for the
quasilinear system $u_{t}=Lu_{x}$. The above observation means that $f$ is a
conservation laws for every power of $L$ or, equivalently, generates a
hierarchy of conservation laws for $L$.

In particular, this implies that
$(L^{*}-\lambda\operatorname{Id})^{-1}\mathrm{d}\,f$ is also closed for any
$\lambda$. Thus, there locally exists a function $g(\lambda;u)$ depending on
$\lambda$ as a parameter and such that
$\mathrm{d}\,g=(L^{*}-\lambda\operatorname{Id})^{-1}\mathrm{d}\,f$. Here
$u=(u^{1},\dots,u^{n})$ are coordinates on $\mathrm{M}^{n}$ and $\mathrm{d}\,$
denotes the differential of a function w.r.t. $u$. One of the main ingredients
of our construction is the function $\sigma(\lambda;u)=e^{g(\lambda;u)}$
satisfying the identity

$\big{(}L^{*}-\lambda\operatorname{Id}\big{)}\mathrm{d}\,\sigma=\sigma\,\mathrm{d}\,f,\qquad\lambda\in\mathbb{C}.$
(29)

If $L(u)-\lambda\operatorname{Id}$ is invertible, i.e.
$\lambda\notin\operatorname{Spectrum}L(u)$, then $\sigma(\lambda;u)$ is
analytic in $\lambda$, otherwise the point $(\lambda,u)$ may be singular (like
pole or zero, or branching point).

Next, assume that there exist a vector field $\zeta$ on $M^{n}$ and constant
$C\in\mathbb{R}$ such that555This condition is quite non-trivial so that the
existence of a non-zero $\zeta$ depends on $L$ and $f$. However, for $\zeta=0$
the construction still makes sense but reduces to a hydrodynamic type system,
see Corollary 3.1.

$\mathcal{L}_{\zeta}\bigl{(}\sigma(\lambda;u)\bigr{)}+C\sigma(\lambda;u)=m(\lambda),$
(30)

that is, the left hand side does not depend on $u$ and hence is a certain
function of $\lambda$ (in the examples discussed below, $m(\lambda)$ is always
a polynomial of degree $\leq n=\dim M$). This triple, i.e., Nijenhuis operator
$L$, conservation law $f$ and vector field $\zeta$ are parameters of the
construction.

Notice that $\sigma(\lambda;u)$ and $m(\lambda)$ satisfying (29), (30) are
defined up to simultaneous multiplication by an arbitrary function
$c(\lambda)$. This kind of scaling is not important for the construction below
and we will treat it as a trivial transformation.

###  3.2 Main theorems

Fix a triple $L$, $f$ and $\zeta$. Construct $\sigma(\lambda;u)$ and
$m(\lambda)$ by formulas (29) and (30). Using them construct infinite
differential series $\mathrm{v}(\lambda)$ and $\mathrm{w}(\lambda)$ by (17)
and (18).

Recall that the series $\mathrm{w}(\lambda)$ satisfies the Gelfand-Dikii
identity (20):

$m(\lambda)\left(\mathrm{w}_{xx}(\lambda)\mathrm{w}(\lambda)-\frac{1}{2}(\mathrm{w}_{x}(\lambda))^{2}\right)+\sigma(\lambda)(\mathrm{w}(\lambda))^{2}=1.$
(31)

Based on this information, we introduce an $n$-component system with a
differential constraint666The equation of Type I from Introduction is exactly
of this kind with $q=\mathrm{w}(\lambda)$ for a specific choice of parameters
$L$, $f$ and $\zeta$.

$\displaystyle u_{t_{\lambda}}$
$\displaystyle=\mathrm{w}_{xxx}(\lambda)\big{(}L-\lambda\operatorname{Id}\big{)}^{-1}\zeta+\mathrm{w}(\lambda)\big{(}L-\lambda\operatorname{Id}\big{)}^{-1}u_{x},$
(32) $\displaystyle 0$
$\displaystyle=m(\lambda)\left(\mathrm{w}_{xx}(\lambda)\mathrm{w}(\lambda)-\frac{1}{2}(\mathrm{w}_{x}(\lambda))^{2}\right)+\sigma(\lambda)(\mathrm{w}(\lambda))^{2}-1.$

In this construction one can naturally make sense of $\lambda=\infty$. Namely,
we define $\bar{\sigma}(\lambda)=c(\lambda)\sigma(\frac{1}{\lambda})$ and
$\bar{m}(\lambda)=c(\lambda)m(\frac{1}{\lambda})$, where $c(\lambda)$ is a
suitable scaling factor. The function $\bar{\sigma}(\lambda)$ does not satisfy
(29), but does satisfy a very similar relation

$(\operatorname{Id}-\lambda
L)\mathrm{d}\,\bar{\sigma}(\lambda)=-\lambda\bar{\sigma}(\lambda)\mathrm{d}\,f.$
(33)

This implies that $c(\lambda)$ can be chosen in such a way that
$\bar{\sigma}(\lambda)$ is analytic in $\lambda$ in a neighbourhood of zero
and moreover, $\bar{\sigma}(\lambda)=1-\lambda f+\dots$, where dots denote
higher order terms in $\lambda$. We will assume that $c(\lambda)$ is chosen in
this way. Then we set
$\bar{\mathrm{v}}(\lambda)=\mathrm{v}(\bar{\sigma}(\lambda;u),\bar{m}(\lambda))$
and
$\bar{\mathrm{w}}(\lambda)=\mathrm{w}(\bar{\sigma}(\lambda;u),\bar{m}(\lambda))$
and rewrite the family of equations (32) in the following equivalent form
obtained by replacing $\lambda$ with $\frac{1}{\lambda}$:

$\displaystyle u_{\bar{t}_{\lambda}}$
$\displaystyle=\bar{\mathrm{w}}_{xxx}(\lambda)\big{(}\operatorname{Id}-\lambda
L\big{)}^{-1}\zeta+\bar{\mathrm{w}}(\lambda)\big{(}\operatorname{Id}-\lambda
L\big{)}^{-1}u_{x},$ (34) $\displaystyle 0$
$\displaystyle=\bar{m}(\lambda)\left(\bar{\mathrm{w}}_{xx}(\lambda)\bar{\mathrm{w}}(\lambda)-\frac{1}{2}(\bar{\mathrm{w}}_{x}(\lambda))^{2}\right)+\bar{\sigma}(\lambda)(\bar{\mathrm{w}}(\lambda))^{2}-1.$

More precisely, we have
$\bar{\mathrm{w}}(\lambda)=\frac{1}{\sqrt{c(\lambda)}}\mathrm{w}(\frac{1}{\lambda})$
so that $u_{\bar{t}_{\lambda}}$ coincides with $u_{t_{\lambda^{-1}}}$ up to
appropriate rescaling (with a factor depending on $\lambda$).

This transformation allows us to set $\lambda=0$ in (34) which will naturally
correspond to $\lambda=\infty$ in (32). In particular, we set

$\mathrm{v}(\infty)=\bar{\mathrm{v}}(0)\quad\mbox{and}\quad\mathrm{w}(\infty)=\bar{\mathrm{w}}(0).$
(35)

However, $\bar{\sigma}(0)=1$ leading to $\bar{\mathrm{w}}(0)=1$ and hence to
the trivial evolutionary equation $u_{\bar{t}_{\lambda}}=u_{x}$. To get a non-
trivial equation ‘at infinity’, we may consider the derivative of (34) at
$\lambda=0$, namely we set

$\displaystyle u_{t_{\infty}}$ $\displaystyle=\lim_{\lambda\to
0}\frac{1}{\lambda}(u_{\bar{t}_{\lambda}}-u_{\bar{t}_{0}})=$
$\displaystyle=\frac{1}{\lambda}\left(\bar{\mathrm{w}}_{xxx}(\lambda)\big{(}\operatorname{Id}-\lambda
L\big{)}^{-1}\zeta+\bar{\mathrm{w}}(\lambda)\big{(}\operatorname{Id}-\lambda
L\big{)}^{-1}u_{x}-u_{x}\right)=$
$\displaystyle=q_{xxx}\zeta+(L+q\operatorname{Id})u_{x},$

where $\bar{\mathrm{w}}(\lambda)=1+\lambda q+\mbox{\small{(higher order terms
in $\lambda$)}}$. In other words, $q$ is the derivative of
$\bar{\mathrm{w}}(\lambda)$ w.r.t. $\lambda$ at zero. Substituting
$\bar{\mathrm{w}}(\lambda)=1+\lambda q+\dots$ into the second equation of (34)
we obtain the following constraint for $q$:

$\bar{m}(0)q_{xx}+2q-f(u)=0.$

To summarise, for $\lambda=\infty$ we consider the following evolutionary PDE
system with a constraint:

$\displaystyle u_{t_{\infty}}=q_{xxx}\zeta+(L+q\operatorname{Id})u_{x},$ (36)
$\displaystyle 0=\bar{m}(0)q_{xx}+2q-f(u).$

Equations (32) and (36) (related to $\lambda\neq\infty$ and $\lambda=\infty$
respectively) are now understood as a parametric family with
$\lambda\in\bar{\mathbb{C}}=\mathbb{C}\cup\\{\infty\\}$. The main property of
this family of PDEs is as follows.

######  Theorem 2.

Let $L$ be a Nijenhuis operator and $f$ a conservation law of $L$. Consider
$\sigma(\lambda;u)$ constructed from (29), and a vector field $\zeta$
satisfying (30) for a certain function $m(\lambda)$. Then for any
$\lambda,\mu\in\bar{\mathbb{C}}=\mathbb{C}\cup\\{\infty\\}$, the differential
series $\mathrm{v}(\mu)$ defined by (17), (35) is a conservation law density
for the evolutionary flow $u_{t_{\lambda}}$ with a differential constraint
defined by (32), (36). Moreover, if $f$ is generic in the sense that
$\mathrm{d}\,f,L^{*}\mathrm{d}\,f,\dots,(L^{n-1})^{*}\mathrm{d}\,f$ are
linearly independent, then the flows $u_{t_{\lambda}}$’s pairwise commute.

As a straightforward corollary of this construction, we may consider the
‘trivial’ case when $\zeta=0$ and $m(\lambda)=0$. In this situation, the first
term in (32) disappears, but our construction still gives a non-trivial series
of integrable quasilinear systems.

######  Corollary 3.1.

Let $L(u)$ be a Nijenhuis operator, $f(u)$ a conservation law of $L$ and
$\sigma(\lambda;u)$ denote the function satisfying (29). Then the evolutionary
flows

$u_{t_{\lambda}}=\sigma(\lambda;u)\bigl{(}L(u)-\lambda\operatorname{Id}\bigr{)}^{-1}u_{x}$

pairwise commute for all $\lambda$’s. Moreover, the functions
$\dfrac{1}{\sigma(\mu,x)}$ are common conservation law densities for these
flows (for all $\lambda$ and $\mu$).

######  Remark 3.1.

In the assumptions of Theorem 2, we obtain the flows of the form
$\mathrm{w}(\lambda)(L-\lambda\operatorname{Id})^{-1}u_{x}$ with
$\mathrm{w}(\lambda)=\sigma(\lambda)^{-1/2}$ for $m=0$. This exponent
$-\frac{1}{2}$, however, is not very essential. Indeed, if $\sigma(\lambda)$
satisfies (29), then $\sigma(\lambda)^{c}$ satisfies (29) also with $f$
replaced with $\tilde{f}=c\cdot f$, so that Corollary 3.1 can be easily
obtained by an appropriate rescaling. Of course, this corollary admits a
direct proof without using Theorem 2.

We also note that Corollary 3.1 can be understood as a $\lambda$-version of
the construction by F. Magri suggested in [19] and then developed in [17]. If
$f=\mathrm{tr}\,L$ and $\lambda\to\infty$, then we obtain the system studied
by E. Ferapontov and M. Pavlov in [12] (see also [6]) and for
$f=c\cdot\mathrm{tr}\,L$, $c\in\mathbb{R}$ we obtain the so-called
$\varepsilon$-systems studied by M. Pavlov [21].

We will need another corollary from Theorem 2. Consider a formal PDE

$u_{t_{\lambda}}=\xi(\lambda)$ (37)

obtained from (32) by resolving the constraint w.r.t. $\mathrm{w}(\lambda)$,
i.e., expressing $\mathrm{w}(\lambda)$ as a formal differential series and
substituting it into the first equation of (32). As a result, the r.h.s. of
(37) becomes a formal differential series in the derivatives
$u_{x},u_{xx},\dots$ whose coefficients are functions in $\lambda$ and $u$.

We now fix $\lambda$ and, in a small neighbourhood of it, expand
$\xi(\lambda+\varepsilon)$ in powers of $\varepsilon$:

$\xi(\lambda+\varepsilon)=\sum_{s=0}^{\infty}\xi_{s,\lambda}\varepsilon^{s}.$
(38)

In the same way we defines $\varepsilon$-expansions for $\mathrm{v}(\lambda)$
and $\mathrm{w}(\lambda)$:

$\mathrm{v}(\lambda+\varepsilon)=\sum_{s=0}^{\infty}\mathrm{v}_{s,\lambda}\varepsilon^{s},\quad\mathrm{w}(\lambda+\varepsilon)=\sum_{s=0}^{\infty}\mathrm{w}_{s,\lambda}\varepsilon^{s}.$
(39)

Notice that by construction, each coefficient $\xi_{s,\lambda}$,
$\mathrm{v}_{s,\lambda}$ or $\mathrm{w}_{s,\lambda}$ is still a formal
differential series in $u_{x},u_{xx},\dots$.

######  Corollary 3.2.

In the settings of Theorem 2, assume that the conservation law $f$ is generic
in the sense that
$\mathrm{d}\,f,L^{*}\mathrm{d}\,f,\dots,(L^{n-1})^{*}\mathrm{d}\,f$ are
linearly independent. Then the (formal) evolutionary flows defined by the
(formal) vector fields $\xi_{s,\lambda}$ ($\lambda\in\bar{\mathbb{C}}$,
$s=0,1,2,\dots$) pairwise commute. Moreover, $\mathcal{v}_{r,\mu}$ are common
(formal) conservation law densities for all of them ($\mu\in\bar{\mathbb{C}}$,
$r=0,1,2,\dots$).

The next theorem is closely related to Corollary 3.2 and deals with
degeneration of the differential constraints that we observed in Section 2 for
Type I and Type III, but now in the general case.

######  Theorem 3.

In the settings of Theorem 2, let $\lambda_{i}$ be a zero of $m(\lambda)$,
i.e., $m(\lambda_{i})=0$ $(i\in\\{1,2,\dots,k\\})$. Then all the coefficients
$\xi_{s,\lambda_{i}}$, $\mathrm{v}_{s,\lambda_{i}}$ and
$\mathrm{w}_{s,\lambda_{i}}$ of $\varepsilon$-expansions (38) and (39) at the
point $\lambda_{i}$ are differential polynomials so that
$u_{t}=\xi_{s,\lambda_{i}}$ is a usual evolutionary equation as in (1). In
particular, for $s=0$ these equations have the following form

$u_{t_{\lambda_{i}}}=\xi(\lambda_{i})=\left(\frac{1}{\sqrt{\sigma(\lambda_{i})}}\right)_{xxx}(L-\lambda_{i}\operatorname{Id})^{-1}\zeta+\frac{1}{\sqrt{\sigma(\lambda_{i})}}(L-\lambda_{i}\operatorname{Id})^{-1}u_{x},\qquad\mbox{for
$\lambda_{i}\neq\infty$}$ (40)

and

$u_{t_{\infty}}=\xi(\infty)=\tfrac{1}{2}f_{xxx}\,\zeta+\left(L+\tfrac{1}{2}f\operatorname{Id}\right)u_{x},\qquad\mbox{for
$\lambda_{i}=\infty$}.$ (41)

Summarising the statements of Theorems 2, 3 and Corollary 3.2 we come to the
following conclusion. For each
$\lambda\in\bar{\mathbb{C}}=\mathbb{C}\cup\\{\infty\\}$ we define an
evolutionary multi-component PDE system (32), (36) with a differential
constraint as in (2). The corresponding (formal) evolutionary flows
$u_{t_{\lambda}}=\xi(\lambda)$ pairwise commute and admit an infinite family
of common (formal) conservation laws also parameterised by
$\lambda\in\bar{\mathbb{C}}$. For some special values of the parameter
$\lambda$, namely for the zeros $\lambda_{1},\lambda_{2},\dots$ of the
function $m(\lambda)$ ($\infty$ is also allowed when appropriate), the
corresponding PDE equations $u_{t_{\lambda_{i}}}=\xi(\lambda_{i})$ are usual
evolutionary multi-component PDEs whose r.h.s. are differential polynomials as
in (1). Each $\lambda_{i}$ generates hierarchies of commuting non-formal
symmetry fields $\xi_{s,\lambda_{i}}$ and non-formal conservation laws
$\mathrm{v}_{s,\lambda_{i}}$ for the whole family
$u_{t_{\lambda}}=\xi(\lambda)$ of formal PDE systems. Moreover, the members of
these hierarchies are defined by means of an explicit iterative procedure in
terms of the function $\sigma(\lambda;u)$ and vector field $\zeta$.

Thus, Theorems 2 and 3 give a recipe for constructing multi-component
integrable PDEs starting from a Nijenhuis operator $L$ and its conservation
law $f$ satisfying certain conditions. However, in order to construct a
specific example of such a system, we need to find a function
$\sigma(u,\lambda)$ and a vector field $\zeta(u)$ satisfying (29) and (30),
that is, to solve a (possibly, infinite) system of PDEs. It is straightforward
to check that the function $\sigma$ and vector field $\zeta$ given by (9) and
(10) in Section 2 are solutions of (29) and (30). The construction from
Theorems 2, 3 applied to these $\sigma$ and $\zeta$ gives the integrable
systems of Types I–IV from Section 2 so that Theorem 1 immediately follows.

The next theorem shows that in the differentially non-degenerate case,
$\sigma$ and $\zeta$ given by (9) and (10) provide the only non-trivial
solution of (29) and (30).

######  Theorem 4.

Let $L$ be a differentially non-degenerate Nijenhuis operator and $f$ a
conservation law of $L$ such that at a point $p\in M^{n}$ the 1-forms
$\mathrm{d}\,f,L^{*}\mathrm{d}\,f,\dots,(L^{n-1})^{*}\mathrm{d}\,f$ are
linearly independent. Assume that there exist $\sigma(\lambda;u)$,
$m(\lambda)$ and $\zeta$ satisfying (29) and (30) with $m(\lambda)\neq 0$.
Then, in a small neighborhood of $p$, the Nijenhuis operator $L$,
functions777We recall that $\sigma(\lambda;u)$ and $m(\lambda)$ are defined up
to simultaneous multiplication by an arbitrary function $c(\lambda)$ and this
freedom is assumed here. $\sigma(\lambda;u)$, $m(\lambda)$ and vector field
$\zeta$ are as in Section 2, see (8), (9), (7) and (10) respectively.

As already mentioned above, Theorem 1 follows directly from Theorems 2, 3 by
taking $\sigma(\lambda)$ and $\zeta$ given by (9) and (10). Theorems 2 and 3
are proved in Section 4 and Theorem 4 in Section 5.

###  3.3 Parameters corresponding to known integrable systems

In this section we show that for particular choice of the parameters, Type I –
IV equations from Section 2 contain many famous integrable systems so that our
approach allows one to generate a vast amount of different integrable systems
in a unifying manner.

######  Example 3.1 (KdV, Camassa-Holm, Dullin-Gottwald-Holm and their
generalisations).

In dimension $n=1$, the differentially non-degenerate Nijenhuis operator is
$L=u$. Due to Theorem 4, the only possible $\sigma(\lambda)$ is $u-\lambda$
and then $\zeta=m_{0}+m_{1}u$. Notice that

$\mathcal{L}_{\zeta}\sigma(\lambda)=m_{0}+m_{1}u=m_{0}+\lambda
m_{1}+m_{1}(u-\lambda)=m(\lambda)+m_{1}\sigma(\lambda),$

as required by (30).

The Type I equation in this case is

$\displaystyle u_{t}$
$\displaystyle=q_{xxx}\frac{m_{0}+m_{1}u}{u-\lambda}+q\frac{u_{x}}{u-\lambda},$
$\displaystyle 0$ $\displaystyle=(m_{0}+\lambda
m_{1})q_{xx}q-\frac{1}{2}(m_{0}+\lambda m_{1})(q_{x})^{2}+(u-\lambda)q^{2}-1.$

This is a three-parameter ($m_{0}$, $m_{1}$ and $\lambda$) family of
integrable evolutionary PDEs with differential constraint.

If $m_{1}\neq 0$, then taking $\lambda_{0}=-\frac{m_{0}}{m_{1}}$ we get Type
II equation

$u_{t}=m_{1}\Big{(}\frac{1}{\sqrt{u-\lambda_{0}}}\Big{)}_{xxx}+\frac{u_{x}}{(u-\lambda_{0})^{3/2}}.$

This is a two-parameter ($\lambda_{0}$ and $m_{1}\neq 0$)888The parameter
$\lambda_{0}$ is not essential unless we consider the limit as
$\lambda_{0}\to\infty$. family of equations. For $\lambda_{0}=0$ it yields
(after rescaling) the reduction of the coupled Harry Dym equation [2, eqn.
26a]. It also appeared in [9] as the first flow of the inverse Camassa-Holm
hierarchy (flow $m^{(0)}_{t}$ in Section “Bihamiltonian structure” from [9]).

The Type III equation takes the form

$\displaystyle u_{t}$ $\displaystyle=q_{xxx}(m_{0}+m_{1}u)+(u+q)u_{x},$
$\displaystyle\frac{u}{2}$ $\displaystyle=-\frac{m_{1}}{2}q_{xx}+q.$

This is a two-parameter family of the PDEs with a constraint. Differentiating
the constraint we get an expression $m_{1}q_{xxx}=2q_{x}-u_{x}$. Substituting
it into the first equation and renaming the coefficient we obtain Dullin-
Gottwald-Holm equation [11]. The case $m_{0}=0$ gives the Camassa-Holm
equation.

Finally, the Type IV equation corresponds to $m_{1}=0$, leading to the
celebrated KdV equation

$u_{t}=\frac{m_{0}}{2}u_{xxx}+\frac{3}{2}uu_{x}.$

######  Example 3.2 (Coupled KdV and Harry Dym, Kaup-Boussinesq and Ito
systems).

Take an arbitrary $n$ and consider the (differentially non-degenerate)
Nijenhuis operator $L=L_{\mathsf{comp}}$ given by the first formula of (6). In
the notation of Section 2, we take $N=0$, $n=\ell_{0}=d$ and
$m(\lambda)=m_{n}\lambda^{n}+m_{n-1}\lambda^{n-1}+\dots+m_{0}$. Then

$\displaystyle\sigma(\lambda)$
$\displaystyle=\operatorname{det}(L-\lambda\operatorname{Id})=(-1)^{n}(\lambda^{n}-u^{1}\lambda^{n-1}-\dots-u^{n}),$
$\displaystyle\zeta$
$\displaystyle=(-1)^{n+1}\left((m_{n-1}+m_{n}u^{1})\frac{\partial}{\partial
u^{1}}+\dots+(m_{0}+m_{n}u^{n})\frac{\partial}{\partial u^{n}}\right).$

It is easy to check that
$\mathcal{L}_{\zeta}\sigma(\lambda)=m(\lambda)-(-1)^{n}m_{n}\sigma(\lambda)$.
For every root $\lambda_{i}$ of $m(\lambda)$, the Type II equation is

$u_{t}=\Bigg{(}\frac{1}{\sqrt{\det(L-\lambda_{i}\operatorname{Id})}}\Bigg{)}_{xxx}(L-\lambda_{i}\operatorname{Id})^{-1}\zeta+\frac{1}{{\sqrt{\det(L-\lambda_{i}\operatorname{Id})}}\big{)}}(L-\lambda_{i}\operatorname{Id})^{-1}u_{x}.$

This is a $(n+1)$-parameter family of integrable equations (with $m_{i}$’s as
parameters involved in the formula for $\zeta$ above). For the rather special
case $m_{n}=1$ and $m_{0}=\dots=m_{n-1}=0$, we get $\lambda_{i}=0$ and taking
$L=L_{\mathsf{comp}}$ as in (6) we come to coupled Harry Dym equations

$u_{t}=\left(\frac{1}{\sqrt{u^{n}}}\right)_{xxx}e_{1}\,+\,\left(\frac{1}{\sqrt{u^{n}}}\right)L_{\mathsf{comp}}^{-1}\,u_{x},\quad
e_{1}=\begin{pmatrix}1\\\ 0\\\ \vdots\\\ 0\end{pmatrix},\
u=\begin{pmatrix}u^{1}\\\ u^{2}\\\ \vdots\\\ u^{n}\end{pmatrix},\ $

introduced in [2] by M. Antonowicz and A. Fordy.

If $m_{n}=0$, then $m(\lambda)$ has a root at infinity and Type IV equation is

$u_{t}=\tfrac{1}{2}\big{(}\operatorname{tr}L\big{)}_{xxx}\zeta+\big{(}L+\tfrac{1}{2}\operatorname{tr}L\operatorname{Id}\big{)}u_{x}.$

This is a family of integrable multi-component PDE systems with $n$ parameters
$m_{0},\dots,m_{n-1}$. More specifically, for $L=L_{\mathsf{comp}}$ given by
(6) we get

$u_{t}=\tfrac{1}{2}\,u^{1}_{xxx}\,\zeta+\big{(}L_{\mathrm{comp}}+\tfrac{1}{2}\,u^{1}\operatorname{Id}\big{)}u_{x},\quad\zeta=\begin{pmatrix}m_{n-1}\\\
\vdots\\\ m_{1}\\\ m_{0}\end{pmatrix}=\sum_{i=1}^{n}m_{n-i}e_{i},\
m_{i}\in\mathbb{R}.$

For $\zeta=e_{i}$, $i=1,\dots,n$, we get $n$ different systems known as
coupled KdV systems and introduced by Antonowicz and Fordy in [1].

The latter have two important examples for $n=2$. For $m_{2}=m_{1}=0,m_{0}\neq
0$ after coordinate change $u^{1}\to-u^{1},u^{2}\to-u^{2}$ we get the Kaup-
Boussinesq system [22, eqn. (4)]:

$\displaystyle u^{1}_{t}$ $\displaystyle=u^{2}_{x}-\frac{3}{2}u^{1}u^{1}_{x},$
$\displaystyle u^{2}_{t}$
$\displaystyle=\frac{m_{0}}{2}u^{1}_{xxx}-u^{2}u^{1}_{x}-\frac{1}{2}u^{1}u^{2}_{x}.$

For $m_{0}=m_{2}=0,m_{1}\neq 0$ after coordinate change
$u^{1}\to-u^{1},u^{2}\to-u^{2}$ the same formula yields Ito system [22, eqn.
(25)]:

$\displaystyle u^{1}_{t}$
$\displaystyle=\frac{m_{1}}{2}u^{1}_{xxx}-\frac{3}{2}u^{1}u^{1}_{x}+u^{2}_{x},$
$\displaystyle u^{2}_{t}$
$\displaystyle=-u^{2}u^{1}_{x}-\frac{1}{2}u^{1}u^{2}_{x}.$

######  Example 3.3 (Marvan-Pavlov system).

Now consider a pair of differentially non-degenerate Nijenhuis operators
$L_{0},L_{1}$ in dimensions $n_{0}$ and $n_{1}$. Assume $n_{0}-n_{1}=d\geq 0$
and consider coordinates $u^{1},\dots,u^{n_{0}}$ and $v^{1},\dots,v^{n_{1}}$
in which $L_{0}$ and $L_{1}$ are given by the first formula of (6). In the
notations of Section 2, take $N=1$ and $\ell_{1}=1$. We get
$m(\lambda)=m_{d}\lambda^{d}+m_{d-1}\lambda^{d-1}+\dots+m_{0}$. In these
coordinates

$\sigma(\lambda)=(-1)^{d}\frac{\lambda^{n_{0}}-u^{1}\lambda^{n_{0}-1}-\dots-u^{n_{0}}}{\lambda^{n_{1}}-v^{1}\lambda^{n_{1}-1}-\dots-v^{n_{1}}}$

and

$\begin{aligned}
\zeta&=(-1)^{d}m_{d}\Bigg{(}\sum_{r=1}^{n_{1}}v^{r}\frac{\partial}{\partial
u^{r}}-\sum_{j=1}^{n_{0}}u^{j}\frac{\partial}{\partial
u^{j}}\Bigg{)}+\sum_{s=1}^{d}(-1)^{d}m_{d-s}\Bigg{(}-\frac{\partial}{\partial
u^{s}}+\sum_{j=1}^{n_{1}}v^{j}\frac{\partial}{\partial
u^{j+s}}\Bigg{)}\end{aligned}.$

By direct computation we have
$\zeta(\sigma(\lambda))=m(\lambda)-(-1)^{d}\sigma(\lambda)$. For every root
$\lambda_{i}$ of $m(\lambda)$, the Type II equation is

$\displaystyle
u_{t}=\Bigg{(}\sqrt{(-1)^{d}\frac{\lambda_{i}^{n_{1}}-\sum_{r=1}^{n_{1}}v^{r}\lambda_{i}^{n_{1}-r}}{\lambda_{i}^{n_{0}}-\sum_{j=1}^{n_{0}}u^{j}\lambda_{i}^{n_{0}-j}}}\Bigg{)}_{xxx}$
$\displaystyle(L-\lambda_{i}\operatorname{Id})^{-1}\zeta+$
$\displaystyle+\sqrt{(-1)^{d}\frac{\lambda_{i}^{n_{1}}-\sum_{r=1}^{n_{1}}v^{r}\lambda_{i}^{n_{1}-r}}{\lambda_{i}^{n_{0}}-\sum_{j=1}^{n_{0}}u^{j}\lambda_{i}^{n_{0}-j}}}(L-\lambda_{i}\operatorname{Id})^{-1}u_{x}.$

This is a $d$-parameter family of integrable equations (with the coefficients
of $m(\lambda)$ as parameters involved into the formula for $\zeta$). If
$m_{d}=0$ and the infinity is a root of $m(\lambda)$, we get the Type IV
equation

$u_{t}=\frac{1}{2}\big{(}\operatorname{tr}L_{1}-\operatorname{tr}L_{2}\big{)}_{xxx}\zeta+\Big{(}L+\frac{1}{2}\big{(}\operatorname{tr}L_{1}-\operatorname{tr}L_{2}\big{)}\operatorname{Id}\Big{)}u_{x}.$

This is a $d$-parameter family of integrable equations. Taking $m_{0}\neq 0$
and all $m_{i}=0,i\geq 1$ yields example by M. Marvan and M. Pavlov [20, 22].

######  Example 3.4 (Two-component Camassa-Holm and Dullin-Gottwald-Holm
systems).

Fix $n=2$ and consider the Nijenhuis operator of the from

$L=\left(\begin{array}[]{cc}2u^{1}&u^{2}\\\ u^{2}&0\end{array}\right).$

This operator is related to left-symmetric algebras and plays an important
role in the linearization problem (see [15] for details). We take

$\sigma(\lambda)=\operatorname{det}(L-\lambda\operatorname{Id})=\lambda^{2}-2u^{1}\lambda-(u^{2})^{2},\quad
m(\lambda)=m_{2}\lambda^{2}+m_{1}\lambda+m_{0}.$

The vector field $\zeta$ is

$\zeta=-\Big{(}\frac{m_{1}}{2}+m_{2}u^{1}\Big{)}\frac{\partial}{\partial
u^{1}}-\Big{(}\frac{m_{0}}{2u^{2}}+\frac{m_{2}}{2}u^{2}\Big{)}\frac{\partial}{\partial
u^{2}}.$

We get $\mathcal{L}_{\zeta}(\sigma(\lambda))=m(\lambda)-m_{2}\sigma(\lambda)$.
For $m_{2}\neq 0$, the Type III system in this setting is

$\displaystyle u^{1}_{t}$
$\displaystyle=q_{xxx}\Big{(}-\frac{m_{1}}{2}-m_{0}u^{1}\Big{)}+2u^{1}u^{1}_{x}+u^{2}u^{2}_{x}+qu^{1}_{x},$
$\displaystyle u^{2}_{t}$
$\displaystyle=q_{xxx}\Big{(}-\frac{m_{2}}{2u^{2}}-\frac{m_{0}}{2}u^{2}\Big{)}+u^{2}u^{1}_{x}+qu^{2}_{x},$
$\displaystyle u^{1}$ $\displaystyle=\frac{m_{2}}{2}q_{xx}+q.$

Differentiating the last equation and rearranging terms, we get
$-\frac{m_{2}}{2}q_{xxx}=q_{x}-u^{1}_{x}$. Substituting it into the first, we
get the equivalent form of the previous PDE with constraint

$\displaystyle u^{1}_{t}$
$\displaystyle=\frac{m_{1}}{2}q_{xxx}+2u^{1}q_{x}+qu^{1}_{x}+u^{2}u^{2}_{x},$
$\displaystyle u^{2}_{t}$
$\displaystyle=-\frac{m_{0}}{4}\frac{q_{xxx}}{u^{2}}+(qu^{2})_{x},$
$\displaystyle u^{1}$ $\displaystyle=q+\frac{m_{2}}{2}q_{xx}.$

This is a 3-parameter family of integrable systems. For $m_{2}=0$ we obtain
the two-component Dullin-Gottwald-Holm equation [24, (3)]. If, in addition,
$m_{1}=0$, then we get two-component Camassa-Holm equation [10, (3) and (4)].

######  Example 3.5.

We actually can generalise the equations of type I and III to the case of $k$
constraints, $k>1$. Consider the expansion (27) up to $\varepsilon^{2}$ :

$\displaystyle\bar{\mathrm{w}}(\varepsilon)$
$\displaystyle=1+\varepsilon\mathcal{w}_{0,\epsilon}+\varepsilon^{2}\mathcal{w}_{1,\infty}+\dots.$

If $\bar{m}(0)\neq 0$, then $\mathcal{w}_{0,\infty}$ and
$\mathcal{w}_{1,\infty}$ are formal differential series. Now substitute the
decompositions for $\bar{\mathrm{w}}(\varepsilon),\bar{\sigma}(\varepsilon)$
and $\bar{m}(\varepsilon)$ into the Gelfand-Dikii identity (31). Renaming
$q^{1}=\mathcal{w}_{0,\infty}$, $q^{2}=\mathcal{w}_{1,\infty}$ we get

$\displaystyle 0$
$\displaystyle=\bar{m}(\varepsilon)\left(\bar{\mathrm{w}}_{xx}(\varepsilon)\bar{\mathrm{w}}(\varepsilon)-\frac{1}{2}(\bar{\mathrm{w}}_{x}(\varepsilon))^{2}\right)+\bar{\sigma}(\varepsilon)(\bar{\mathrm{w}}(\varepsilon))^{2}-1=\varepsilon\Big{(}\bar{\sigma}_{1}+2q^{1}+\bar{m_{0}}q^{1}_{xx}\Big{)}+$
$\displaystyle+\varepsilon^{2}\Big{(}\bar{\sigma}_{2}+2q^{2}+\bar{m}_{0}q^{2}_{xx}+\bar{m}_{1}q^{1}_{xx}+\bar{m}_{0}q^{1}q^{1}_{xx}+2\bar{\sigma}_{1}q^{1}+(q_{1})^{2}-\frac{1}{2}\bar{m}_{0}(q^{1}_{x})^{2}\Big{)}+\dots.$

This yields differential relations for $q^{1},q^{2}$. Differential operator
$\operatorname{Id}+D^{2}$ is formally invertible, so these constraints imply
that $q^{1}$ is a differential series in $\bar{\sigma}_{1}$ and its
derivatives and $q^{2}$ is a differential series in
$\bar{\sigma}_{1},\bar{\sigma}_{2}$ and their derivatives.

Now consider the expansion (28) up to $\varepsilon^{2}$:

$\displaystyle\bar{\xi}(\varepsilon)$
$\displaystyle=\xi_{-1,\infty}+\varepsilon\xi_{0,\infty}+\varepsilon^{2}\xi_{1,\infty}+\dots=$
$\displaystyle=u_{x}+\varepsilon\Big{(}q^{1}_{xxx}\zeta+(L+q^{1}\operatorname{Id})u_{x}\Big{)}+\varepsilon^{2}\Big{(}q^{2}_{xxx}\zeta+q^{1}_{xxx}L\zeta+(L^{2}+q^{1}L+q^{2}\operatorname{Id})u_{x}\Big{)}+\dots.$

Taking $\xi_{1,\infty}$ we get a PDE with two differential constraints

$\displaystyle u_{t}$
$\displaystyle=q^{2}_{xxx}\zeta+q^{1}_{xxx}L\zeta+(L^{2}+q^{1}L+q^{2}\operatorname{Id})u_{x},$
$\displaystyle 0$
$\displaystyle=\bar{\sigma}_{1}+2q^{1}+\bar{m}_{0}q^{1}_{xx},$ $\displaystyle
0$
$\displaystyle=\bar{\sigma}_{2}+2q^{2}+\bar{m}_{0}q^{2}_{xx}+\bar{m}_{1}q^{1}_{xx}+\bar{m}_{0}q^{1}q^{1}_{xx}+2\bar{\sigma}_{1}q^{1}+(q_{1})^{2}-\frac{1}{2}\bar{m}_{0}(q^{1}_{x})^{2}.$

For $L$ differentially non-degenerate, $N=0$ and $m(t)=t^{n}$ this yields the
general form of Camassa-Holm equation CH($n$,$2$) from [14]. Taking the
expansion up to a higher order, one obtains a greater number of differential
constraints.

######  Example 3.6.

In the previous examples, we have shown that many notable integrable systems
are special cases of the systems from Section 2. We now describe one of the
simplest new examples. By construction, integrable systems we deal with are
written in invariant form that is independent on the choice of a local
coordinate chart. In particular, in order to make our system more symmetric,
we may choose local coordinates related to the roots of the polynomial
$m(\lambda)$.

The next 3-component example is build starting with
$\sigma(\lambda)=\det(L-\lambda\operatorname{Id})$ and
$m(\lambda)=(\lambda_{1}-\lambda)(\lambda_{2}-\lambda)(\lambda_{3}-\lambda)$
and choosing local coordinates $(u^{1},u^{2},u^{3})$ to be

$u^{i}=\det(\lambda_{i}\operatorname{Id}-L)\prod_{s\neq
i}\frac{1}{\lambda_{i}-\lambda_{s}}.$

In particular, $u^{i}$ is proportional to $\sigma(\lambda_{i})$ with some
constant factor which is not essential. Moreover,
$\zeta=-(u^{1},u^{2},u^{3})^{\top}$ and
$(L-\lambda_{i}\operatorname{Id})^{-1}\zeta=e_{i}$. Now if we take an
arbitrary linear combination of the (commuting) evolutionary PDEs (12) of Type
II, we get the following integrable system:

$\begin{pmatrix}u_{t}^{1}\\\ u_{t}^{2}\\\
u_{t}^{3}\end{pmatrix}=\begin{pmatrix}c_{1}\left({1}/{\sqrt{u^{1}}}\right)_{xxx}\\\
c_{2}\left({1}/{\sqrt{u^{2}}}\right)_{xxx}\\\
c_{3}\left({1}/{\sqrt{u^{3}}}\right)_{xxx}\end{pmatrix}+A(u)\begin{pmatrix}u^{1}_{x}\\\
u^{2}_{x}\\\ u^{3}_{x}\end{pmatrix}.$ (42)

where $c_{1},c_{2},c_{3}$ are arbitrary constants and $A(u)$ is the $3\times
3$ matrix with the components

$A^{j}_{i}=\frac{u^{j}}{\lambda_{i}-\lambda_{j}}\left(\frac{c_{i}}{(u^{i})^{3/2}}-\frac{c_{j}}{(u^{j})^{3/2}}\right),\quad
j\neq i,\quad\mbox{and}\quad
A^{i}_{i}=-\frac{c_{i}}{(u^{i})^{3/2}}-\sum_{j\neq i}A^{j}_{i}.$

One can also write it as follows:

$u^{j}_{t}=c_{j}\left(\left(\frac{1}{\sqrt{u^{j}}}\right)_{xxx}-\frac{u^{j}_{x}}{(u^{j})^{3/2}}\right)+\sum_{i\neq
j}\frac{u^{j}u^{i}_{x}-u^{i}u^{j}_{x}}{\lambda_{i}-\lambda_{j}}\left(\frac{c_{i}}{u_{i}^{3/2}}-\frac{c_{j}}{u_{j}^{3/2}}\right).$

Here the first term represents the system of three uncoupled Harry Dym type
equations (see Example 3.1), but the second term mixes all the variables.

The recursion formula for the conservation laws gives the following explicit
formula for the first six of them. The first two corresponding to the root
$\lambda_{1}$ are

$\sqrt{u^{1}}\quad\mbox{and}\quad\frac{\left({u^{1}}u^{1}_{xx}-\frac{5}{4}\left(u^{1}_{x}\right)^{2}\right)-2{\left(u^{1}\right)}^{2}\left(\frac{u^{1}+u^{2}}{\lambda_{1}-\lambda_{2}}+\frac{u^{1}+u^{3}}{\lambda_{1}-\lambda_{3}}-1\right)}{(u^{1})^{3/2}}.$
(43)

The other four correspond to the roots $\lambda_{2}$, $\lambda_{3}$ and can be
obtained from the above formulas by cyclic permutation of indices 1,2,3.

One can also find, using the procedure described in Section 2, the formulas
for commuting flows. Actually, the commuting flows of the lowest order are
(42) with arbitrarily chosen $c_{1},c_{2},c_{3}$. Notice that this example can
be naturally generalised to the case of an arbitrary number of components.

##  4 Proofs of Theorems 2 and 3

We start with the following Lemma.

######  Lemma 4.1.

Under the assumptions of Theorem 2, we have

$(L^{*}_{\lambda})^{-1}\mathrm{d}\,\sigma(\mu)=\frac{1}{\mu-\lambda}\Bigg{(}\mathrm{d}\,\sigma(\mu)-\frac{\sigma(\mu)}{\sigma(\lambda)}\mathrm{d}\,\sigma(\lambda)\Bigg{)}.$

###### Proof.

Condition (29) reads

$L^{*}\mathrm{d}\,\sigma(\mu)=\sigma(\mu)\mathrm{d}\,f+\mu\mathrm{d}\,\sigma(\mu)=\sigma(\mu)\mathrm{d}\,f+(\mu-\lambda)\mathrm{d}\,\sigma(\mu)+\lambda\mathrm{d}\,\sigma(\mu).$
(44)

Recall that by construction
$(L_{\lambda}^{*})^{-1}\mathrm{d}\,f=\frac{1}{\sigma(\lambda)}\mathrm{d}\,\sigma(\lambda).$
Rearranging the terms and multiplying both sides of (44) by
$(L^{*}_{\lambda})^{-1}$, we get the statement of Lemma. ∎

Now let us recall some basic formulas and introduce some notations. We denote
the derivative coordinates of order $j$ by $u^{\alpha}_{x^{j}}$ and set
$u^{\alpha}_{x^{0}}=u^{\alpha}$. Consider a formal evolutionary vector field
$\xi$ with components $\xi^{i}$. The (Lie) derivative of a formal differential
series $\mathrm{w}$ along $\xi$ is defined by

$\mathcal{L}_{\xi}\mathrm{w}=\sum\limits_{j=0}^{\infty}\frac{\partial\mathrm{w}}{\partial
u^{\alpha}_{x^{j}}}D^{j}(\xi^{\alpha})$

with summation over $\alpha$ assumed, $\alpha=1,\dots,n$.

Let $\xi=\xi(\lambda)$ be the formal vector field defined by (37) and
associated with the PDEs from Theorem 2. We have:

$\displaystyle\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=\frac{\partial\sigma(\mu)}{\partial
u^{\alpha}}\xi^{\alpha}(\lambda)=\frac{\partial\sigma(\mu)}{\partial
u^{\alpha}}\Big{(}\mathrm{w}_{xxx}(\lambda)\big{(}L-\lambda\operatorname{Id}\big{)}^{-1}\zeta+\mathrm{w}(\lambda)\big{(}L-\lambda\operatorname{Id}\big{)}^{-1}u_{x}\Big{)}^{\alpha}=$
$\displaystyle=\mathrm{w}_{xxx}(\lambda)\frac{\partial\sigma^{\mu}}{\partial
u^{\alpha}}\Big{(}L^{-1}_{\lambda}\Big{)}^{\alpha}_{q}\zeta^{q}+\mathrm{w}(\lambda)\frac{\partial\sigma^{\mu}}{\partial
u^{\alpha}}\Big{(}L^{-1}_{\lambda}\Big{)}^{\alpha}_{q}u^{q}_{x}=\frac{1}{\mu-\lambda}\mathrm{w}_{xxx}(\lambda)\Big{(}m(\mu)-\frac{\sigma(\mu)}{\sigma(\lambda)}m(\lambda)\Big{)}+$
$\displaystyle+\frac{1}{\mu-\lambda}\mathrm{w}(\lambda)\Big{(}\sigma_{x}(\mu)-\frac{\sigma(\mu)}{\sigma(\lambda)}\sigma_{x}(\lambda)\Big{)}.$

The last step follows from Lemma 4.1. Further rearranging terms and using
(19), we get

$\displaystyle\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=$ (45)
$\displaystyle=\frac{1}{\mu-\lambda}\Bigg{(}m(\mu)\mathrm{w}_{xxx}(\lambda)+w(\lambda)\sigma_{x}(\mu)-\frac{\sigma(\mu)}{\sigma(\lambda)}\Big{(}m(\lambda)\mathrm{w}_{xxx}(\lambda)+\mathrm{w}(\lambda)\sigma_{x}(\lambda)\Big{)}\Bigg{)}=$
$\displaystyle=\frac{1}{\mu-\lambda}\Bigg{(}m(\mu)\mathrm{w}_{xxx}(\lambda)+w(\lambda)\sigma_{x}(\mu)-\frac{\sigma(\mu)}{\sigma(\lambda)}\Big{(}-2\sigma(\lambda)\mathrm{w}_{x}(\lambda)\Big{)}\Bigg{)}=$
$\displaystyle=\frac{1}{\mu-\lambda}\Bigg{(}m(\mu)\mathrm{w}_{xxx}(\lambda)+2\sigma(\mu)\mathrm{w}_{x}(\lambda)+\sigma_{x}(\mu)\mathrm{w}(\lambda)\Bigg{)}.$

Now let us proceed with the proof. Consider a pair of formal differential
series $\mathrm{w}_{1},\mathrm{w}_{2}$. We use notation
$\mathrm{w}_{1}\sim\mathrm{w}_{2}$, if there exists a formal differential
series $\mathrm{u}$, such that $\mathrm{w}_{1}-\mathrm{w}_{2}=D\mathrm{u}$. In
particular, the Leibnitz rule for $D$ implies, that
$\mathrm{w}_{1}D(\mathrm{w}_{2})\sim-D(\mathrm{w}_{1})\mathrm{w}_{2}$. More
generally, it yields the formula

$\mathrm{w}_{1}D^{j}(\mathrm{w}_{2})\sim(-1)^{j}D^{j}(\mathrm{w}_{1})\mathrm{w}_{2}.$

Using the chain rule $\frac{\delta\mathrm{v}(\mu)}{\delta
u^{\alpha}}=\mathrm{w}(\mu)\frac{\partial\sigma(\mu)}{\partial u^{\alpha}}$,
we get the following sequence of relations

$\displaystyle\mathcal{L}_{\xi(\lambda)}\mathrm{v}(\mu)=\sum\limits_{j=0}^{\infty}\frac{\partial\mathrm{v}(\mu)}{\partial
u^{\alpha}_{x^{j}}}D^{j}(\xi^{\alpha}(\lambda))\sim\sum\limits_{j=0}^{\infty}(-1)^{j}D^{j}\Big{(}\frac{\partial\mathrm{v}(\mu)}{\partial
u^{\alpha}_{x^{j}}}\Big{)}\xi^{\alpha}(\lambda)=\frac{\delta\mathrm{v}(\mu)}{\delta
u^{\alpha}}\xi^{\alpha}(\lambda)=$
$\displaystyle=\mathrm{w}(\mu)\frac{\partial\sigma(\mu)}{\partial
u^{\alpha}}\xi^{\alpha}(\lambda)=\mathrm{w}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=\frac{1}{\mu-\lambda}m(\mu)\mathrm{w}(\mu)\mathrm{w}_{xxx}(\lambda)+\frac{2}{\mu-\lambda}\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\sigma(\mu)+$
$\displaystyle+\frac{1}{\mu-\lambda}\mathrm{w}(\mu)\mathrm{w}(\lambda)\sigma_{x}(\mu)\sim-\frac{1}{\mu-\lambda}m(\mu)\mathrm{w}_{xxx}(\mu)\mathrm{w}(\lambda)-\frac{2}{\mu-\lambda}\big{(}\mathrm{w}(\mu)\sigma(\mu)\big{)}_{x}\mathrm{w}(\lambda)+$
$\displaystyle+\frac{1}{\mu-\lambda}\mathrm{w}(\mu)\sigma_{x}(\mu)\mathrm{w}(\lambda)=-\frac{1}{\mu-\lambda}\Bigg{(}m(\mu)\mathrm{w}_{xxx}(\lambda)+2\sigma(\mu)\mathrm{w}_{x}(\lambda)+\sigma_{x}(\mu)\mathrm{w}(\lambda)\Bigg{)}\mathrm{w}(\lambda)=0$

Here we used formula (45) and, again, identity (19). Thus, we get
$\mathcal{L}_{\xi(\lambda)}\mathrm{v}(\mu)\sim 0$, meaning that
$\mathrm{v}(\mu)$ is a formal conservation law for the flow $\xi(\lambda)$.
For $\lambda=\infty$ the proof is essentially the same.

Now let us proceed to the commuting flows. We will need the following Lemma.

######  Lemma 4.2.

Assume that $\sigma(\mu)$ and $\mathrm{w}(\mu)$ are related by Gelfand-Dikii
identity (31) and $\mathcal{L}_{\xi(\lambda)}\sigma(\mu)$ is given by (45).
Then

$\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)=\frac{1}{\mu-\lambda}\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}.$
(46)

###### Proof.

We start with applying $\xi(\lambda)$ to Gelfand-Dikii identity (31) and
multiplying the result by $\mathrm{w}(\mu)$ (we also take into account the
fact that $\mathcal{L}_{\xi}$ commute with $D$):

$\displaystyle 0$
$\displaystyle=\mathrm{w}(\mu)\mathcal{L}_{\xi(\lambda)}\Bigg{(}m(\mu)\left(\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)-\frac{1}{2}(\mathrm{w}_{x}(\mu))^{2}\right)+\sigma(\mu)\mathrm{w}^{2}(\mu)\Bigg{)}=$
(47)
$\displaystyle=m(\mu)\Big{(}\mathrm{w}^{2}(\mu)\mathcal{L}_{\xi(\lambda)}\mathrm{w}_{xx}(\mu)+\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)-\mathrm{w}_{x}(\mu)\mathrm{w}(\mu)\mathcal{L}_{\xi(\lambda)}\mathrm{w}_{x}(\mu)\Big{)}+$
$\displaystyle+2\sigma(\mu)\mathrm{w}^{2}(\mu)\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)+\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=$
$\displaystyle=\left(m(\mu)\Big{(}\mathrm{w}^{2}(\mu)D^{2}-\mathrm{w}_{x}(\mu)\mathrm{w}(\mu)D+\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)\Big{)}+2\sigma(\mu)\mathrm{w}^{2}(\mu)\right)\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)+$
$\displaystyle+\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=\cal
R\big{(}\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)\big{)}+\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)$

with

$\displaystyle\mathcal{R}$
$\displaystyle=m(\mu)\Big{(}\mathrm{w}^{2}(\mu)D^{2}-\mathrm{w}_{x}(\mu)\mathrm{w}(\mu)D+\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)\operatorname{Id}\Big{)}+2\sigma(\mu)\mathrm{w}^{2}(\mu)\operatorname{Id}=$
$\displaystyle=\big{(}2-m(\mu)\mathrm{w}(\mu)\mathrm{w}_{xx}(\mu)+m(\mu)\mathrm{w}^{2}_{x}(\mu)\big{)}\operatorname{Id}-m(\mu)\mathrm{w}(\mu)\mathrm{w}_{x}(\mu)D+m(\mu)\mathrm{w}^{2}(\mu)D^{2},$

where in the latter relation we substitute
$\sigma(\mu)\mathrm{w}^{2}(\mu)=1-m(\mu)\left(\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)-\frac{1}{2}(\mathrm{w}_{x}(\mu))^{2}\right)$
from (31). Thus, we have the identity

$\mathcal{R}\big{(}\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)\big{)}=-\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu).$
(48)

Note that $\cal R$ is a (formally) invertible differential operator.
Therefore, it suffices to verify that
$\mathcal{L}_{\xi(\lambda)}\mathrm{w}(\mu)$ defined by (46) satisfies (48) or,
equivalently,

$\cal
R\big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\big{)}+(\mu-\lambda)\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=0.$

Computing the l.h.s. of this relation gives:

$\displaystyle\ \ \
\,m(\mu)\mathrm{w}^{2}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}_{xx}+m(\mu)\mathrm{w}_{xx}(\mu)\mathrm{w}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}-$
$\displaystyle-m(\mu)\mathrm{w}_{x}(\mu)\mathrm{w}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}_{x}+2\sigma(\mu)\mathrm{w}^{2}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}+$
$\displaystyle+(\mu-\lambda)\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=$
$\displaystyle=m(\mu)\mathrm{w}^{2}(\mu)\Big{(}\mathrm{w}_{xxx}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{xxx}(\lambda)\Big{)}+2\sigma(\mu)\mathrm{w}^{2}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{x}(\lambda)\Big{)}$
$\displaystyle+(\mu-\lambda)\mathrm{w}^{3}(\mu)\mathcal{L}_{\xi(\lambda)}\sigma(\mu).$

Adding and subtracting $\sigma_{x}(\mu)\mathrm{w}(\lambda)\mathrm{w}^{3}(\mu)$
we arrive to the identity

$\displaystyle\Big{(}(\mu-\lambda)\mathcal{L}_{\xi(\lambda)}\sigma(\mu)-$
$\displaystyle
m(\mu)\mathrm{w}_{xxx}(\lambda)-2\sigma(\mu)\mathrm{w}_{x}(\lambda)-\sigma_{x}(\mu)\mathrm{w}(\lambda)\Big{)}\mathrm{w}^{3}(\mu)+$
$\displaystyle+$
$\displaystyle\Big{(}m(\mu)\mathrm{w}_{xxx}(\mu)+2\sigma(\mu)\mathrm{w}_{x}(\mu)+\sigma_{x}(\mu)\mathrm{w}(\mu)\Big{)}\mathrm{w}^{2}(\mu)\mathrm{w}(\lambda)=0,$

where the first term vanishes due to (45) and the second due to (19),
completing the proof. ∎

As we deal with evolutionary vector fields, it is enough to check that
$\xi(\lambda)$ and $\xi(\nu)$ commute, acting on coordinate functions. Fix
three pairwise distinct $\lambda,\mu,\nu$. From Lemma 4.2, we get

$\displaystyle\mathcal{L}_{\xi(\nu)}\mathcal{L}_{\xi(\lambda)}\sigma(\mu)=\frac{1}{\mu-\lambda}\mathcal{L}_{\xi(\nu)}\Big{(}m(\mu)\mathrm{w}_{xxx}(\lambda)+2\sigma(\mu)\mathrm{w}_{x}(\lambda)+\sigma_{x}(\mu)\mathrm{w}(\lambda)\Big{)}=$
$\displaystyle=\frac{m(\mu)}{(\mu-\lambda)(\lambda-\nu)}\Big{(}\mathrm{w}_{xxxx}(\mu)\mathrm{w}(\lambda)+2\mathrm{w}_{xxx}(\mu)\mathrm{w}_{x}(\lambda)-\mathrm{w}(\mu)\mathrm{w}_{xxxx}(\lambda)-2\mathrm{w}_{x}(\mu)\mathrm{w}_{xxx}(\lambda)\Big{)}+$
$\displaystyle+\frac{2}{(\mu-\lambda)(\mu-\nu)}\Big{(}m(\mu)\mathrm{w}_{xxx}(\nu)+2\sigma(\mu)\mathrm{w}_{x}(\nu)+\sigma_{x}(\mu)\mathrm{w}(\nu)\Big{)}\mathrm{w}_{x}(\lambda)+$
$\displaystyle+\frac{1}{(\mu-\lambda)(\lambda-\nu)}\Big{(}2\sigma(\mu)\Big{(}\mathrm{w}_{xx}(\mu)\mathrm{w}(\nu)-\mathrm{w}(\mu)\mathrm{w}_{xx}(\nu)\Big{)}+\sigma_{x}(\mu)\Big{(}\mathrm{w}_{x}(\mu)\mathrm{w}(\nu)-\mathrm{w}(\mu)\mathrm{w}_{x}(\nu)\Big{)}\Big{)}+$
$\displaystyle+\frac{1}{(\mu-\lambda)(\mu-\nu)}\Big{(}m(\mu)\mathrm{w}_{xxxx}(\nu)+2\sigma_{x}(\mu)\mathrm{w}_{x}(\nu)+2\sigma(\mu)\mathrm{w}_{xx}(\nu)+\sigma_{xx}(\mu)\mathrm{w}(\nu)+\sigma_{x}(\mu)\mathrm{w}_{x}(\nu)\Big{)}\mathrm{w}(\lambda)$

The identity

$\frac{1}{(\mu-\lambda)(\mu-\nu)}-\frac{1}{(\mu-\lambda)(\lambda-\nu)}+\frac{1}{(\mu-\nu)(\lambda-\nu)}=0.$

implies that the r.h.s. of formula for
$\mathcal{L}_{\xi(\nu)}\mathcal{L}_{\xi(\lambda)}\sigma(\mu)$ is symmetric in
$\lambda,\nu$. Thus, $\mathcal{L}_{\xi(\nu)}$ and $\mathcal{L}_{\xi(\lambda)}$
commute on $\sigma(\mu)$.

Now recall that $\sigma(\mu)$ is constructed from a generic conservation law
$f$. This implies that in the expansion

$\sigma(\mu)=\sigma_{0}+\mu\sigma_{1}+\dots$

the differentials of $\sigma_{0},\dots,\sigma_{n-1}$ are linearly independent
almost everywhere. Thus, one can take them as coordinates $u^{i}=\sigma_{i-1}$
and in these coordinates
$(\mathcal{L}_{\xi(\nu)}\mathcal{L}_{\xi(\lambda)}-\mathcal{L}_{\xi(\lambda)}\mathcal{L}_{\xi(\nu)})u^{i}=0$,
as required. This completes the proof of Theorem 2.

To verify the statement of Theorem 3, we first need to show that the
coefficients $\xi_{s,\lambda_{i}}$, $\mathrm{v}_{s,\lambda_{i}}$ and
$\mathrm{w}_{s,\lambda_{i}}$ of the $\varepsilon$-expansions (38) and (39) are
well defined and are differential polynomials in $u^{1},\dots,u^{n}$. Indeed,
by definition,

$\mathrm{v}(\lambda)=\mathrm{v}(\sigma(\lambda),m(\lambda))=\sqrt{2}\sum_{s=0}^{\infty}\bigl{(}m(\lambda)\bigr{)}^{s}\mathcal{v}_{s}(\lambda),$
(49)

where $\mathcal{v}_{s}(\lambda)$ is a differential polynomial obtained from
the homogeneous differential polynomial
$\mathcal{u}_{2s+1}(\sigma,\sigma_{x},\sigma_{xx},\dots)$ of degree $2s$ by
substitution $\sigma=\sigma(\lambda;u)$.

We are interested in the expansion of
$\mathrm{v}(\lambda_{i}+\varepsilon)=\sum\mathrm{v}_{s,\lambda_{i}}\varepsilon^{s}$
under the assumption that $m(\lambda_{i})=0$. Since
$m(\lambda_{i}+\varepsilon)=a_{1}\varepsilon+a_{2}\varepsilon^{2}+\dots$ and,
therefore,
$\bigl{(}m(\lambda_{i}+\varepsilon)\bigr{)}^{s}=a_{1}^{s}\varepsilon^{s}+\dots$,
we see from expansion (49) that $\mathrm{v}_{s,\lambda_{i}}$ is defined from
the first $s+1$ coefficients $\mathcal{v}_{0},\dots,\mathcal{v}_{s}$. Hence,
$v_{s,\lambda_{i}}$ is a non-homogeneous differential polynomial of degree at
most $2s$.

The proof for $\mathrm{w}_{s,\lambda_{i}}$ is literally the same. The
conclusion for $\mathrm{w}_{s,\lambda_{i}}$ immediately follows from the
explicit formula of $\xi(\lambda)$ in terms of $\mathrm{w}(\lambda)$, see (32)
and (37).

The explicit form (40) of the flows $u_{t_{\lambda_{i}}}=\xi(\lambda_{i})$ for
a root $\lambda_{i}\in\mathbb{C}$ of $m(\lambda)$ is straightforward. Indeed,
setting $m(\lambda_{i})=0$ in (32) gives
$\sigma(\lambda_{i})\bigr{(}\mathrm{w}(\lambda_{i})\bigr{)}^{2}-1=0$, or
equivalently, $\mathrm{w}(\lambda_{i})=\frac{1}{\sqrt{\sigma(\lambda_{i})}}$.
Substituting this expression into the first equation of (32) gives (40), as
required. Similarly, for $\lambda_{i}=\infty$ we set $\bar{m}(0)=0$ in (36) to
get $q=\frac{1}{2}f$, which after substitution into the first equation of (36)
immediately gives (41), completing the proof of Theorem 3.

##  5 Proof of Theorem 4

Since $L$ is differentially non-degenerate, this operator is diagonalisable
almost everywhere. At “non-diagonalisable” points, the conclusion of Theorem 4
can be derived by continuity arguments. So w.l.o.g. we assume that
$L=\operatorname{diag}(x^{1},\dots,x^{n})$. Recall that $f$ from equation (29)
is a conservation law for $L$. Then $f$ is a sum of $n$ functions such that
the $i$th function depends on $x^{i}$ only. Therefore, for every $i$ the
function $f_{i}:=\tfrac{\partial f}{\partial x^{i}}$ depends on $x^{i}$ only.

Next, we consider relation (29). In coordinates, it reads:

$\displaystyle\frac{\partial\ln(\sigma)}{\partial x^{1}}$ $\displaystyle=$
$\displaystyle\frac{f_{1}(x^{1})}{x^{1}-\lambda}$ $\displaystyle\vdots$
$\displaystyle\frac{\partial\ln(\sigma)}{\partial x^{n}}$ $\displaystyle=$
$\displaystyle\frac{f_{n}(x^{n})}{x^{n}-\lambda}.$

Hence, the system (29) of $n$ PDEs is actually a system of $n$ ODEs in
different variables. Its solution must be of the form
$\sigma=c(\lambda)\cdot\sigma_{1}\dots\sigma_{n}$ with

$\sigma_{i}(\lambda,x^{i})=\exp\left(\int_{s_{i}}^{x^{i}}\frac{f_{i}(s)}{s-\lambda}\mathrm{d}\,s\right)$
(50)

where $c(\lambda)$ is an arbitrary function and $(s_{1},...,s_{n})$ is an
arbitrary point; we assume that all $s_{i}\neq 0$.

Next, consider relation (30). For our
$\sigma(\lambda)=c(\lambda)\cdot\sigma_{1}\dots\sigma_{n}$ it reads

$\left(C+\sum_{i=1}^{n}\zeta^{i}\frac{f_{i}}{x^{i}-\lambda}\right)\,\sigma_{1}\dots\sigma_{n}=\frac{m(\lambda)}{c(\lambda)}:=\widehat{m}(\lambda).$
(51)

In the left hand side of this relation, $f_{i}$ and $\zeta_{i}$ are smooth
functions in $x$ which are independent on $\lambda$, whereas the r.h.s. is a
function independent of $x$. The following statement shows that under these
conditions, $f_{i}$’s have to be constants and, moreover, very special.

######  Lemma 5.1.

The functions $f_{i}$ are integer constants different from zero and no greater
than $1$. Moreover, for every $i$ such that $f_{i}\neq 1$ we have
$\zeta_{i}=0$.

###### Proof.

Integration by parts gives

$\sigma_{i}(\lambda,x^{i})=\exp\left(\int_{s_{i}}^{x^{i}}\frac{f_{i}(s)}{s-\lambda}\mathrm{d}\,s\right)=\exp\left(f_{i}(x^{i})\ln(x^{i}-\lambda)-f_{i}(s_{i})\ln(s_{i}-\lambda)-\int_{s_{i}}^{x^{i}}f^{\prime}_{i}(s)\ln(s-\lambda)\mathrm{d}\,s\right)$

implying
$\sigma(x,\lambda)=(x^{1}-\lambda)^{f_{1}(x^{1})}(x^{2}-\lambda)^{f_{2}(x^{2})}\dots(x^{n}-\lambda)^{f_{n}(x^{n})}\tilde{\sigma}(\lambda,x)$
where the function $\tilde{\sigma}$ has neither zeros nor poles. The equation
(51) reads then

$\left(C+\sum_{s=1}^{n}\frac{\zeta^{s}f_{s}(x^{s})}{x^{s}-\lambda}\right)(x^{1}-\lambda)^{f_{1}(x^{1})}(x^{2}-\lambda)^{f_{2}(x^{2})}\dots(x^{n}-\lambda)^{f_{n}(x^{n})}\tilde{\sigma}(\lambda,x)=\widehat{m}(\lambda).$
(52)

Note that the function
$\int_{s_{i}}^{x^{i}}f^{\prime}_{i}(s)\ln(s-\lambda)\mathrm{d}\,s$ is locally
bounded, so the function $\tilde{\sigma}(\lambda,x)$ is bounded for small
$x^{s}-\lambda$ and is not zero.

Assume for a certain $i$ that $f_{i}$ is not integer at a point
$(x^{1},...,x^{n})$. Substituting $\lambda=x^{i}$ (for generic $x^{i}$) makes
$\widehat{m}(\lambda)=0$ or $\widehat{m}(\lambda)=\infty$ leading to a
contradiction. Indeed, the first factor of (52) has integer order of zeros and
poles, so to compensate it $f_{i}$ must be integer for each $x^{i}$ and
therefore constant. Note that the case $\widehat{m}(\lambda)\equiv\infty$ is
not allowed since the left hand side is finite almost everywhere.

Thus, all $f_{i}$ are integer. Then the function $\tilde{\sigma}(\lambda,x)$
depends on $\lambda$ only. Further we assume that it is equal to $1$ since we
can “hide” it in $\widehat{m}(\lambda)$ (we keep the same notation). Our
equation (52) then becomes:

$\left(C+\sum_{s=1}^{n}\frac{\zeta^{s}f_{s}}{x^{s}-\lambda}\right)(x^{1}-\lambda)^{f_{1}}(x^{2}-\lambda)^{f_{2}}\dots(x^{n}-\lambda)^{f_{n}}=\widehat{m}(\lambda),$
(53)

where $f_{i}$ are integer constants. If $\zeta^{i}\neq 0$ then the first
factor of (53) has a pole of order $1$ implying $f_{i}=1$. If $\zeta^{i}=0$
and $f_{i}>0$, then (53) has zero for $\lambda=x^{i}$ which again leads to
contradiction. It remains to notice that $f_{i}\neq 0$ since
$\mathrm{d}\,f,L^{*}\mathrm{d}\,f,\dots,L^{n-1}\mathrm{d}\,f$ are linearly
independent by our assumption. ∎

W.l.o.g. we assume that $f_{i}=1$ for $i=1,...,k$ and the other $f_{i}$ are
negative integers; we set $\ell_{k+1}=-f_{k+1},\dots,\ell_{n}=-f_{n}$. The
equation (53) reads then

$\left(C+\sum_{s=1}^{k}\frac{\zeta^{s}}{x^{s}-\lambda}\right)=\frac{(x^{k+1}-\lambda)^{\ell_{k+1}}\dots(x^{n}-\lambda)^{\ell_{n}}}{(x^{1}-\lambda)(x^{2}-\lambda)\dots(x^{k}-\lambda)}\,\widehat{m}(\lambda).$
(54)

Notice that the expression in the l.h.s. can be written as a rational function
in $\lambda$ of the form $\frac{P(\lambda)}{\prod(x^{s}-\lambda)}$, where
$P(\lambda)=C(-\lambda)^{k}+\dots$ is a polynomial of degree at most $k$.
Similarly, the r.h.s. is $\frac{F(\lambda)}{\prod(x^{s}-\lambda)}$ where
$F(\lambda)=(x^{k+1}-\lambda)^{\ell_{k+1}}\dots(x^{n}-\lambda)^{\ell_{n}}\widehat{m}(\lambda)$.
Since $P(\lambda)=F(\lambda)$ we conclude that
$\widehat{m}(\lambda)=m_{d}\lambda^{d}+\dots$ is a polynomial of degree at
most $d=k-\ell_{k+1}-\dots-\ell_{n}$. In particular, $d\geq 0$ and
$m_{d}=(-1)^{d}C$.

Finally, it remains to notice that $\zeta^{i}$ can be found from (54) by using
the partial fraction decomposition theorem which gives:

$\zeta^{i}=\frac{(x^{k+1}-x^{i})^{\ell_{k+1}}\dots(x^{n}-x^{i})^{\ell_{n}}\widehat{m}(x^{i})}{\prod_{s=1,s\neq
i}^{k}(x^{s}-x^{i})},\qquad i=1,\dots,k.$ (55)

Summarising this discussion we conclude that for
$L=\operatorname{diag}(x^{1},\dots,x^{n})$, we have (up to scaling with a
factor $c(\lambda)$)

$\displaystyle f(x)$
$\displaystyle=\sum_{i=1}^{k}x^{i}-\sum_{i=k+1}^{n}\ell_{i}x^{i}$
$\displaystyle\sigma(\lambda,x)$
$\displaystyle=\frac{\prod_{i=1}^{k}(x_{i}-\lambda)}{\prod_{i=k+1}^{n}(x_{i}-\lambda)^{\ell_{i}}}$

and $\zeta$ given by (55).

To get the conclusion of Theorem 4 in its final form, we only need to combine
the coordinates $x^{i}$’s, $i=k+1,\dots,n$ into groups depending on the values
of the exponents $\ell_{i}$ ($x^{i}$ and $x^{j}$ go to the same group iff
$\ell_{i}=\ell_{j}$). After appropriate renumbering, we come to the desired
description of $L$, $\sigma$, $m$ and $\zeta$.

##  6 Conclusion

In our paper, we constructed new explicit families of integrable multi-
component evolutionary equations with and without differential constraints,
see Section 2. The equations
$u_{t_{\lambda}}=\xi(\lambda;u,u_{x},u_{xx},\dots)$ within each family are
parametrised by $\lambda\in\bar{\mathbb{C}}=\mathbb{C}\cup\\{\infty\\}$. The
corresponding (formal) evolutionary flows pairwise commute for all values of
parameters and admit a family of common (formal) conservation laws
$\mathrm{v}(\mu;u,u_{x},u_{xx},\dots)$ also parametrised by
$\mu\in\bar{\mathbb{C}}$. For specific values of $\lambda$ (namely, roots
$\lambda_{1},\dots,\lambda_{d}$ of a certain polynomial) the above equation
generates an hierarchy of usual (i.e. non-formal) commuting flows
$\xi_{s,\lambda_{i}}$, $s=0,1,\dots$, defined by means of a differential
polynomial of degree $2s+3$. The equation
$u_{t_{\lambda_{i}}}=\xi(\lambda_{i};u,u_{x},u_{xx},u_{xxx})$, initial term of
this hierarchy, is a third order PDE system given by an elegant explicit
formula. All the subsequent terms can be found by means of explicit recurrent
formulas. Similar for conservation laws: $\mathrm{v}(\lambda_{i})$ generates
an hierarchy of common polynomial conservation laws
$\mathrm{v}_{s,\lambda_{i}}$, $s=0,1,\dots$, for all the flows, where
$\mathrm{v}_{s,\lambda_{i}}$ is a differential polynomial of degree $2s$ that
can be found explicitly by an iterative procedure. These families of
integrable equations, for a simple choice of parameters, include and
generalise many known examples of integrable systems. Some of multi-component
evolutionary equations we constructed are essentially new and they have no
low-component analogues.

The construction is based on a new approach, which is rather differential-
geometric than algebraic (in contrast to many other constructions of
integrable systems which are often based on algebraic or algebraic-geometric
concepts). Our results have been naturally obtained within the Nijenhuis
Geometry programme initiated in [8]. This suggests that Nijenhuis Geometry
might be a convenient framework for studying further properties of the
constructed systems and generalizing them. Because of its differential-
geometric nature, our constructions are invariant with respect to the choice
of coordinates on $\mathrm{M}^{n}$; that is, the systems behave covariantly if
we change unknown functions $u$ by a diffeomorphism
$u_{\mathsf{new}}=u_{\mathsf{new}}(u_{\mathsf{old}})$. One can use this fact
in the search of applications of our systems in natural sciences. Examples
discussed in Section 3.3 actually suggest that ‘physically relevant’ variables
correspond to those coordinates on $\mathrm{M}^{n}$ in which the Nijenhuis
operator has a nice form, e.g., the ‘companion’ and ‘diagonal’ forms from (6),
or the form in which the components of $L$ are linear in coordinates.

The famous integrable systems that we generalise to an arbitrary number of
components (such as KdV, Camassa-Holm, Dullin-Gottwald-Holm, Harry Dym, Kaup-
Boussinesq) have been intensively studied for decades; for these studies, a
number of non-trivial geometric, algebraic and analytical methods were
invented and successfully applied. The next natural step would be to figure
out how to adapt these methods to new systems. In particular, it would be
interesting to construct Lax representations for new systems, to find explicit
solutions by the inverse scattering method, to construct a recursion operator
and, of course, to find physically relevant models that are described by new
systems. We invite our fellow mathematicians and physicists to join this
research.

## References

  * [1] M. Antonowicz, A. Fordy, Coupled KdV equations with multi-Hamiltonian structures, Physica D: Nonlinear Phenomena, 28 (1987), 3, 345–357.
  * [2] M. Antonowicz, A. Fordy, Coupled Harry Dym equations with multi-Hamiltonian structures, Journal of Physics A: Mathematical and General, 21(1988), 5, L269–75.
  * [3] M. Antonowicz, A. Fordy, Factorisation of energy dependent Schrödinger operators: Miura maps and modified systems, Communication in Mathematical Physics, 124 (1989), 465–486.
  * [4] A. V. Bolsinov, A. Yu. Konyaev, V. S. Matveev, Nijenhuis Geometry III: gl-regular Nijenhuis operators, arXiv:2007.09506.
  * [5] A. V. Bolsinov, A. Yu. Konyaev, V. S. Matveev, Applications of Nijenhuis geometry: non-degenerate singular points of Poisson–Nijenhuis structures, European Journal of Mathematics, (2021), arXiv: 2001.04851.
  * [6] A. V. Bolsinov, A. Yu. Konyaev, V. S. Matveev, Applications of Nijenhuis geometry II: maximal pencils of multi-Hamiltonian structures of hydrodynamic type, Nonlinearity, 34 (2021), 8, 5136–5162.
  * [7] A. V. Bolsinov, A. Yu. Konyaev, V. S. Matveev, Applications of Nijenhuis geometry III: Frobenius pencils and compatible non-homogeneous Poisson structures, arXiv: 2112.09471v2.
  * [8] A. V. Bolsinov, A. Yu. Konyaev, V. S. Matveev, Nijenhuis geometry, Advances in Mathematics, 394(2022), 22, 108001.
  * [9] R. Camassa, D. Holm, An integrable shallow water equation with peaked solitons, Phys. Rev. Lett. 71(1993), 1661–1664.
  * [10] A. Constantin, R. Ivanov, On an integrable two-component Camassa-Holm shallow water system. Physics Letters A, 372(2008), 7129–7132.
  * [11] H. R. Dullin, G. A. Gottwald, D. D. Holm, An integrable shallow water equation with linear and nonlinear dispersion, Physical Review Letters, 87 (2001), 194501.
  * [12] E. V. Ferapontov, M. V. Pavlov, Quasiclassical limit of coupled KdV equations. Riemann invariants and multi-Hamiltonian structure, Physica D: Nonlinear Phenomena, 52 (1991), 2–3, 211–219, https://doi.org/10.1063/1.2982232
  * [13] I. M. Gelfand, L. A. Dikii, Asymptotic behavior of the resolvent of Sturm-Liouville equations and the algebra of the Korteweg-De Vries equations, Russ. Math. Surveys 30 (1975), 5, 77–113, Usp. Mat. Nauk 30(5) (1975), 67–10.
  * [14] D. Holm, R. Ivanov, Multi-component generalizations of the CH equation: geometrical aspects, peakons and numerical examples, J. Phys. A: Math. Theor. 43(2010), 492001.
  * [15] A. Yu. Konyaev, Nijenhuis geometry II: Left-symmetric algebras and linearization problem for Nijenhuis operators, Differential Geometry and its Applications, 74(2021), 101706.
  * [16] A. Yu. Konyaev, Geometry of inhomogeneous Poisson brackets, multi-component Harry Dym hierarchies and multi-component Hunter-Saxton equations, arXiv:2112.05635.
  * [17] P. Lorenzoni, F. Magri, A cohomological construction of integrable hierarchies of hydrodynamic type, International Mathematics Research Notices, 34 (2005), 2087–2100.
  * [18] F. Magri, A simple model of the integrable Hamiltonian equation, Journal of Mathematical Physics, 19 (1978), 1156–1162.
  * [19] F. Magri, Lenard chains for classical integrable systems, Theoretical and Mathematical Physics, 137 (2003), 424–432 (Russian); English transl.: Theoret. and Math. Phys., 137 (2003), 1716–1722.
  * [20] M. Marvan, M. V. Pavlov, Integrable dispersive chains and their multi-phase solutions, Letters in Mathematical Physics, 109 (2019), 1219–1245.
  * [21] M. Pavlov, Integrable hydrodynamic chains, Journal of Math. Physics, 44 (2003), 9, 4134–4156.
  * [22] M. V. Pavlov, Integrable dispersive chains and energy dependent Schrödinger operator, J. Phys. A: Math. Theor. 47 (2014), 295204.
  * [23] M. Pedroni, V. Sciacca, J. Zubelli, The Bi-Hamiltonian Theory of the Harry Dym Equation, Theoretical and Mathematical Physics, 133 (2002), 1585–1597.
  * [24] Yong Chen, Hongjun Gao, Yue Liu, On the Cauchy problem for the two-component Dullin-Gottwald-Holm system. Discrete and Continuous Dynamical Systems, 33 (2013) 8, 3407–3441.