Hierarchies of Manakov–Santini Type by Means of Rota–Baxter and Other
Identities

Hierarchies of Manakov–Santini Type
by Means of Rota–Baxter and Other Identities

Błażej M. SZABLIKOWSKI

B.M. Szablikowski

Faculty of Physics, Adam Mickiewicz University, Umultowska 85, 61-614 Poznań,
Poland<EMAIL_ADDRESS>
Received January 11, 2016, in final form February 22, 2016; Published online
February 27, 2016

The Lax–Sato approach to the hierarchies of Manakov–Santini type is formalized
in order to extend it to a more general class of integrable systems. For this
purpose some linear operators are introduced, which must satisfy some
integrability conditions, one of them is the Rota–Baxter identity. The theory
is illustrated by means of the algebra of Laurent series, the related
hierarchies are classified and examples, also new, of Manakov–Santini type
systems are constructed, including those that are related to the
dispersionless modified Kadomtsev–Petviashvili equation and so called
dispersionless $r$-th systems.

Manakov–Santini hierarchy; Rota–Baxter identity; classical $r$-matrix
formalism; generalized Lax hierarchies; integrable $(2+1)$-dimensional systems

37K10; 37K30

##  1 Introduction

In recent years, one of the significant achievements in the theory of
integrable systems was the construction of formal solutions of the Cauchy
problems for a wide class of $(2+1)$-dimensional dispersionless systems by
means of the inverse scattering transform [14, 15, 16, 17, 18]. In this
process one of the crucial steps was the introduction by Manakov and Santini
of a two-field system that generalizes the dispersionless
Kadomtsev–Petviashvili (KP) equation. The Manakov–Santini system possesses a
non-Hamiltonian Lax pair and the construction of related hierarchy [15] within
the Lax–Sato formalism [6, 8] unifies two original approaches based on
different underlying structures: first by Takasaki and Takebe [33, 34] and the
second one by Martínez Alonso and Shabat [21, 22, 23]. The Manakov–Santini
hierarchy and its generalizations were further studied in several works, see
for instance [6, 7, 8, 9, 10, 25].

The aim of this work is an extension of the Lax–Sato formalism of
Manakov–Santini hierarchy to a more general class of integrable systems, in
particular such as the dispersionless modified KP equation or the so-called
$r$-th systems [2, 3, 4]. Influenced by the papers [7, 8, 9, 10] we generalize
the Lax–Sato formalism of Manakov–Santini hierarchy by means of the Lax
hierarchy (2.1), where two linear operators $\mathcal{P}$ and $\mathcal{R}$
are introduced. In Theorem 2.3 we find the conditions, on the operators
$\mathcal{P}$ and $\mathcal{R}$, for the mutual commutativity of equations
from the hierarchy (2.1). One of the conditions turns out to be the well-known
Rota–Baxter identity [1, 26, 27]. The general source of the relations in
Theorem 2.3 is explained in Section 3. In fact, these relations are directly
connected with those that are used in the work [31] for the construction of
Frobenius manifolds in the cotangent bundles. In Section 4 we illustrate the
above construction by means of the algebra of Laurent series. The related
hierarchies (2.1) are classified and there are presented examples of
integrable systems of Manakov–Santini type, including new ones.

##  2 Generalized hierarchy

Let $\mathbb{A}$ be a commutative associative unital algebra111We assume, for
simplicity, that all structures in this work are defined over the field of
complex numbers.. We define the generalized Lax hierarchy of evolution
equations by

$\displaystyle\Psi_{t_{n}}=A_{n}\Psi_{x}-B_{n}\partial\Psi,\qquad\Psi=\begin{pmatrix}L\\\
M\end{pmatrix},\qquad n\in\mathbb{N},$ (2.1)

where $L,M\in\mathbb{A}$ are Lax and Orlov operators (functions),
respectively. The independent variables are evolution parameters (times)
$t_{n}$ and the spatial variable $x$. We assume that $\partial$ is some
(auxiliary) derivation in the algebra $\mathbb{A}$ invariant with respect to
all $t_{n}$ and $x$. The hierarchy is generated by the functions

$\displaystyle A_{n}:=\mathcal{P}\bigl{(}J^{-1}\partial
X_{n}\bigr{)}\qquad\text{and}\qquad
B_{n}:=\mathcal{R}\bigl{(}J^{-1}(X_{n})_{x}\bigr{)},\qquad X_{n}:=L^{n},$
(2.2)

where $\mathcal{P}$ and $\mathcal{R}$ are some linear maps
$\mathbb{A}\rightarrow\mathbb{A}$ and

$\displaystyle J:=\\{L,M\\}\equiv\partial LM_{x}-L_{x}\partial M.$ (2.3)

We assume that the endomorphisms $\mathcal{P}$ and $\mathcal{R}$ are invariant
with respect to times $t_{n}$ and the spatial variable $x$, that is
$\mathcal{P}$ and $\mathcal{R}$ commute with derivatives related to $t_{n}$
and $x$.

In all the following proofs we will skip most of straightforward computations,
however we will exhibit all the crucial intermediate steps.

######  Proposition 2.1.

The evolution equations from the hierarchy (2.1) pairwise commute if the
following pair of zero-curvature type equations is satisfied by the generating
functions $A_{n}$ and $B_{n}$:

$\displaystyle(A_{n})_{t_{m}}-(A_{m})_{t_{n}}+\langle
A_{n},A_{m}\rangle_{x}+B_{m}\partial A_{n}-B_{n}\partial A_{m}=0$ (2.4a)

and

$\displaystyle(B_{n})_{t_{m}}-(B_{m})_{t_{n}}+A_{n}(B_{m})_{x}-A_{m}(B_{n})_{x}-\langle
B_{n},B_{m}\rangle_{\partial}=0,$ (2.4b)

where

$\displaystyle\langle a,b\rangle_{x}:=ab_{x}-ba_{x},\qquad\langle
a,b\rangle_{\partial}:=a\partial b-b\partial a.$

######  Proof 2.2.

The commuting of the respective flows means that
$(\Psi_{t_{n}})_{t_{m}}=(\Psi_{t_{m}})_{t_{n}}$. Comparing the coefficients of
both sides with respect to the independent variables $\Psi_{x}$ and
$\partial\Psi$ we obtain the required pair of zero-curvature conditions.

######  Theorem 2.3.

The following set of conditions on the endomorphisms $\mathcal{P}$ and
$\mathcal{R}$:

$\displaystyle\mathcal{P}\bigl{(}\mathcal{P}(a)b+a\mathcal{P}(b)\bigr{)}-\mathcal{P}(a)\mathcal{P}(b)=\kappa_{1}\mathcal{P}(ab),$
(2.5a) $\displaystyle\mathcal{P}\bigl{(}\mathcal{R}(a)\partial
b+a\partial\mathcal{P}(b)\bigr{)}-\mathcal{R}(a)\partial\mathcal{P}(b)=\kappa_{1}\mathcal{P}(a\partial
b),$ (2.5b)
$\displaystyle\mathcal{R}\bigl{(}\mathcal{P}(a)b+a\mathcal{R}(b)\bigr{)}-\mathcal{P}(a)\mathcal{R}(b)=\kappa_{2}\mathcal{R}(ab),$
(2.5c) $\displaystyle\mathcal{R}\bigl{(}\mathcal{R}(a)\partial
b+a\partial\mathcal{R}(b)\bigr{)}-\mathcal{R}(a)\partial\mathcal{R}(b)=\kappa_{2}\mathcal{R}(a\partial
b),$ (2.5d)

where $a,b\in\mathbb{A}$ and $\kappa_{1},\kappa_{2}\in\mathbb{C}$, is
sufficient for the zero-curvature equations (2.4) to be identically fulfilled
by the generating functions $A_{n}$ and $B_{n}$. This means that the above
equations are sufficient conditions for mutual commutation of the flows from
the hierarchy (2.1).

######  Proof 2.4.

Let’s define

$\displaystyle a_{n}:=J^{-1}\partial X_{n}\qquad\text{and}\qquad
b_{n}:=J^{-1}(X_{n})_{x},$

so that $A_{n}=\mathcal{P}(a_{n})$ and $B_{n}=\mathcal{R}(b_{n})$. Using (2.1)
one finds that

$\displaystyle(X_{n})_{t_{m}}=A_{m}(X_{n})_{x}-B_{m}\partial X_{n}$

and

$\displaystyle
J_{t_{n}}=(A_{n}J)_{x}-\partial(B_{n}J)\quad\Rightarrow\quad\big{(}J^{-1}\big{)}_{t_{n}}=\langle
A_{n},J^{-1}\rangle_{x}-\langle B_{n},J^{-1}\rangle_{\partial}.$

Hence,

$\displaystyle(a_{n})_{t_{m}}=\langle A_{m},a_{n}\rangle_{x}+b_{n}\partial
A_{m}-B_{m}\partial
a_{n},\qquad(b_{n})_{t_{m}}=A_{m}(b_{n})_{x}-a_{n}(B_{m})_{x}-\langle
B_{m},b_{n}\rangle_{\partial}.$

Now substituting $(A_{n})_{t_{m}}\equiv\mathcal{R}((a_{n})_{t_{m}})$ and
$(B_{n})_{t_{m}}\equiv\mathcal{R}((a_{n})_{t_{m}})$ into the pair of zero-
curvature equations (2.4) and taking advantage of the relations (2.5) one
finds that

$\displaystyle(A_{n})_{t_{m}}+A_{n}(A_{m})_{x}+B_{m}\partial
A_{n}-n\leftrightarrow m$
$\displaystyle\qquad{}=\mathcal{P}\bigl{(}A_{m}(a_{n})_{x}+a_{m}(A_{n})_{x}\bigr{)}-A_{m}(A_{n})_{x}-\mathcal{P}\bigl{(}B_{m}\partial
a_{n}+b_{m}\partial A_{n}\bigr{)}+B_{m}\partial A_{n}-n\leftrightarrow m$
$\displaystyle\qquad{}=\kappa_{1}\mathcal{P}\bigl{(}a_{m}(a_{n})_{x}-a_{n}(a_{m})_{x}-b_{m}\partial
a_{n}+b_{n}\partial a_{m}\bigr{)}$
$\displaystyle\qquad{}=\kappa_{1}\mathcal{P}\bigl{(}J^{-1}\partial(J^{-1}\\{X_{m},X_{n}\\})\bigr{)}=0$

and

$\displaystyle(B_{n})_{t_{m}}-A_{m}(B_{n})_{x}-B_{n}\partial
B_{m}-n\leftrightarrow m$
$\displaystyle\qquad{}=\mathcal{R}\bigl{(}A_{m}(b_{n})_{x}+a_{m}(B_{n})_{x}\bigr{)}-A_{m}(B_{n})_{x}-\mathcal{R}\bigl{(}B_{m}\partial
b_{n}+b_{m}\partial B_{n}\bigr{)}+B_{m}\partial B_{n}-n\leftrightarrow m$
$\displaystyle\qquad{}=\kappa_{2}\mathcal{R}\bigl{(}a_{m}(b_{n})_{x}-a_{n}(b_{m})_{x}-b_{m}\partial
b_{n}+b_{n}\partial b_{m}\bigr{)}$
$\displaystyle\qquad{}=\kappa_{2}\mathcal{R}\bigl{(}J^{-1}(J^{-1}\\{X_{m},X_{n}\\})_{x}\bigr{)}=0.$

In the above $n\leftrightarrow m$ stands for all the remaining terms arising
from the permutation of the indices $m$ and $n$ in the preceding terms. Since
$\\{X_{m},X_{n}\\}=0$ we see that the zero-curvature equations (2.4) are
indeed satisfied identically.

######  Lemma 2.5.

Assume that the endomorphisms $\mathcal{P}$ and $\mathcal{R}$ satisfy the
following relation

$\displaystyle\mathcal{P}(\partial a)=\partial\mathcal{R}(a)$ (2.6)

for arbitrary $a\in\mathbb{A}$. Then, under the constraint

$\displaystyle J\equiv\\{L,M\\}=L^{s},$ (2.7)

where $s$ is some fixed integer, the generalized hierarchy (2.1) reduces to
the standard hierarchy of the form

$\displaystyle L_{t_{n}}=\frac{n}{n-s}\\{\mathcal{R}(L^{n-s}),L\\},$ (2.8)

where the Poisson bracket is defined through the formula (2.3), that is
$\\{,\\}:=\partial\wedge\partial_{x}$.

On the other hand, if we assume that $L$ is invariant with respect to times
$t_{n}$ and the variable $x$, that is $L_{t_{n}}=L_{x}=0$, then the hierarchy
reduces to

$\displaystyle M_{t_{n}}=n\mathcal{P}\bigl{(}M_{x}^{-1}L^{n-1}\bigr{)}M_{x}.$
(2.9)

######  Proof 2.6.

Taking advantage of the constraint (2.7) and the relation (2.6) we have

$\displaystyle A_{n}=\mathcal{P}\bigl{(}J^{-1}\partial
X_{n}\bigr{)}=\frac{n}{n-r}\mathcal{P}\bigl{(}\partial
L^{n-r}\bigr{)}=\frac{n}{n-r}\partial\mathcal{R}(L^{n-r})$

and

$\displaystyle
B_{n}=\mathcal{R}\bigl{(}J^{-1}(X_{n})_{x}\bigr{)}=\frac{n}{n-r}\mathcal{R}\bigl{(}(L^{n-r})_{x}\bigr{)}=\frac{n}{n-r}\bigl{(}\mathcal{R}(L^{n-r})\bigr{)}_{x}.$

The last relation follows from the assumption that $\mathcal{R}$ commutes with
$\partial_{x}$. As result the hierarchy (2.1) reduces into

$\displaystyle\Psi_{t_{n}}=\frac{n}{n-r}\bigl{[}\partial\mathcal{R}(L^{n-r})\Psi_{x}-\bigl{(}\mathcal{R}(L^{n-r})\bigr{)}_{x}\partial\Psi\bigr{]}\equiv\frac{n}{n-r}\\{\mathcal{R}(L^{n-r}),\Psi\\},$

where $\Psi=(L,M)^{\rm T}$. The first equation on $L$ coincides with (2.8).
The second equation gives evolution of $M$ consistent with the constraint
(2.7).

If $L_{t_{n}}=L_{x}=0$, then $B_{n}=0$ and $J=\partial LM_{x}$. Hence

$\displaystyle A_{n}=\mathcal{P}\bigl{(}J^{-1}\partial
X_{n}\bigr{)}=n\mathcal{P}\bigl{(}M_{x}^{-1}L^{n-1}\bigr{)}.$

In this case, the first equation in (2.1) for $\Psi=L$ is satisfied
identically and the second equation for $\Psi=M$ takes the form (2.9).

##  3 Rota–Baxter and other identities

Consider some algebra $(\mathbb{A},\cdot)$. The Rota–Baxter identity [1, 26,
27] for some linear operator
$\mathcal{P}\colon\mathbb{A}\rightarrow\mathbb{A}$ has the form

$\displaystyle\mathcal{P}(\mathcal{P}(a)\cdot
b+a\cdot\mathcal{P}(b))-\mathcal{P}(a)\cdot\mathcal{P}(b)=\kappa\mathcal{P}(a\cdot
b),$ (3.1)

where $a,b,c\in\mathbb{A}$ and $\kappa$ is some fixed scalar weight.
Alternatively we can write the identity as

$\displaystyle\ell(\ell(a)\cdot
b+a\cdot\ell(b))-\ell(a)\cdot\ell(b)=\frac{1}{4}\kappa^{2}a\cdot b,$

where $\ell:=\mathcal{P}-\frac{1}{2}\kappa$. For an operator $\mathcal{P}$
satisfying (3.1) there is always associated operator
$\mathcal{P}^{\prime}:=\kappa-\mathcal{P}$, which satisfies the identity (3.1)
for the same weight $\kappa$.

There is a source of simple solutions to the Rota–Baxter identity (3.1).
Assume that the algebra $\mathbb{A}$ can be decomposed into direct sum of
subalgebras, that is

$\displaystyle\mathbb{A}=\mathbb{A}_{+}\oplus\mathbb{A}_{-},\qquad\mathbb{A}_{\pm}\cdot\mathbb{A}_{\pm}\subset\mathbb{A}_{\pm},\qquad\mathbb{A}_{+}\cap\mathbb{A}_{-}=\\{0\\}.$

Then, the projections $P_{+}$ and $P_{-}$ on the subalgebras $\mathbb{A}_{+}$
and $\mathbb{A}_{-}$ satisfy the identity (3.1) for the weight $\kappa=1$.
Notice that $P_{+}+P_{-}=1$.

The main feature of the Rota–Baxter identity is that in the case of
associative algebra $\mathbb{A}$ the identity (3.1) is a sufficient condition
for associativity of another multiplication in $\mathbb{A}$ given by

$\displaystyle a*_{\mathcal{P}}b:=\mathcal{P}(a)\cdot
b+a\cdot\mathcal{P}(b)-\kappa a\cdot b\equiv\ell(a)\cdot b+a\cdot\ell(b).$

For more information on the Rota–Baxter algebras see [12, 13].

The special case of the Rota–Baxter identity (3.1), for a Lie algebra
$(\mathbb{A},[\cdot,\cdot])$, is the modified Yang–Baxter equation

$\displaystyle\mathcal{R}\bigl{(}[\mathcal{R}(a),b]+[a,\mathcal{R}(b)]\bigr{)}-[\mathcal{R}(a),\mathcal{R}(b)]=\kappa\mathcal{R}\bigl{(}[a,b]\bigr{)}.$
(3.2)

The equation (3.2) is a sufficient condition for the bracket

$\displaystyle[a,b]_{\mathcal{R}}:=[\mathcal{R}(a),b]+[a,\mathcal{R}(b)]-\kappa[a,b]$

to define second Lie bracket in $\mathbb{A}$, in this case the linear map $R$
is called the classical $r$-matrix. The more standard convention is to
consider the endomorphism $r:=\mathcal{R}-\frac{1}{2}\kappa$ instead of
$\mathcal{R}$. The classical $r$-matrix formalism [28, 29] is known to be very
useful in the construction of very broad classes of integrable systems, see
also the survey [5] and references therein.

Let’s now consider commutative and associative algebra $(\mathbb{A},\cdot)$.
We define Poisson bracket in $\mathbb{A}$ be means of two commuting
derivations $\partial,\partial_{x}\in\operatorname{Der}\mathbb{A}$:

$\displaystyle\\{a,b\\}:=\partial
a\cdot\partial_{x}b-\partial_{x}a\cdot\partial b.$ (3.3)

Here, the derivation $\partial_{x}$ is a counterpart of the derivative with
respect to a spatial variable $x$ and the derivation $\partial$ is a
counterpart of the derivative with respect to some auxiliary variable. The
following identity on the endomorphism $\mathcal{R}$ turns out to be
important:

$\displaystyle\mathcal{R}\bigl{(}\mathcal{R}(a)\partial
b\bigr{)}+\mathcal{R}\bigl{(}a\partial\mathcal{R}(b)\bigr{)}-\mathcal{R}(a)\partial\mathcal{R}(b)=\kappa\mathcal{R}(a\partial
b),$ (3.4)

where $\kappa$ is some constant. Analogically as before, $R^{\prime}=\kappa-R$
solves (3.4) for the same weight $\kappa$.

######  Proposition 3.1.

Let us assume that the endomorphism $\mathcal{R}$ commutes with the derivation
$\partial_{x}$, that is $\mathcal{R}\partial_{x}=\partial_{x}\mathcal{R}$.
Then, the identity (3.4) is a sufficient condition for $\mathcal{R}$ to
satisfy the modified Yang–Baxter equation

$\displaystyle\mathcal{R}\bigl{(}\\{\mathcal{R}(a),b\\}+\\{a,\mathcal{R}(b)\\}\bigr{)}-\\{\mathcal{R}(a),\mathcal{R}(b)\\}=\kappa\mathcal{R}\bigl{(}\\{a,b\\}\bigr{)},$
(3.5)

and so to be a classical $r$-matrix for $\mathcal{R}$ with respect to the
Poisson bracket (3.3). This means that

$\displaystyle\\{a,b\\}_{\mathcal{R}}:=\\{\mathcal{R}(a),b\\}+\\{a,\mathcal{R}(b)\\}-\kappa\\{a,b\\}$

is a Lie bracket when the relation (3.4) holds.

The proof is straightforward expanding the formula (3.4) by means of (3.3). In
fact, when $\mathcal{R}$ also commutes with the derivation $\partial$, the
Rota–Baxter identity (3.1) is sufficient for $\mathcal{R}$ to solve the
equation (3.4). However, in general (3.1) is more restrictive than (3.4).

######  Proposition 3.2.

Assume that, on a commutative associative algebra $(\mathbb{A},\cdot)$, there
exists non-degenerate bilinear product
$(\cdot,\cdot)_{\mathbb{A}}\colon\mathbb{A}\times\mathbb{A}\rightarrow\mathbb{C}$,
such that the Frobenius condition holds:

$\displaystyle(a\cdot b,c)_{\mathbb{A}}=(a,b\cdot c)_{\mathbb{A}}$

and the product is invariant with respect to the derivation
$\partial\in\operatorname{Der}\mathbb{A}$:

$\displaystyle(\partial a,b)_{\mathbb{A}}+(a,\partial b)_{\mathbb{A}}=0.$

Then, the Rota–Baxter identity (3.1) is equivalent to the ‘dual’ relation:

$\displaystyle\mathcal{R}\bigl{(}\mathcal{P}(a)b+a\mathcal{R}(b)\bigr{)}-\mathcal{P}(a)\mathcal{R}(b)=\kappa\mathcal{R}(ab),\qquad\mathcal{R}:=\kappa-\mathcal{P}^{*},$

where $(\mathcal{P}^{*}a,b)_{\mathbb{A}}:=(a,\mathcal{P}b)_{\mathbb{A}}$.

Let additionally assume that the following relation is valid:

$\displaystyle\mathcal{P}(\partial a)=\partial\mathcal{R}(a).$ (3.6)

Then, from the Rota–Baxter identity (3.1) we obtain the relation (3.4) and

$\displaystyle\mathcal{P}\bigl{(}\mathcal{R}(a)\partial
b+a\partial\mathcal{P}(b)\bigr{)}-\mathcal{R}(a)\partial\mathcal{P}(b)=\kappa\mathcal{P}(a\partial
b).$

For proof see Proposition 3.3 and Lemma 5.3 in [31]222Notice the misprint in
formula (5.9) in [31], there is missing minus sign on the right-hand side of
the equality.. In the work [31] it is shown that the Rota–Baxter identity
(3.1) and the equation (3.4) can be significant in the construction of
Frobenius manifolds on cotangent bundles inherent to the integrable systems of
hydrodynamic type.

Summarizing, the Rota–Baxter identity (3.1) coincides with the condition
(2.5a) from Theorem 2.3, the sufficient condition (3.4) on $\mathcal{R}$ to
satisfy the modified Yang–Baxter equation (3.5) coincides with (2.5d), and if
the relation (3.6) or (2.6) and some additional natural assumptions are
satisfied, then the remaining conditions (2.5b) and (2.5c) hold automatically.

##  4 Application to the algebra of Laurent series

Consider the algebra of Laurent polynomials
$\mathbb{A}:=\mathbb{C}[p,p^{-1}]$, it is commutative and associative. When
necessary the algebra $\mathbb{A}$ can be extended to the algebra of formal
Laurent series at $\infty$: $\mathbb{A}^{\infty}:=\mathbb{C}(\\!(p^{-1})\\!)$
or the algebra of formal Laurent series at $0$:
$\mathbb{A}^{0}:=\mathbb{C}(\\!(p)\\!)$.

Consider decomposition of $\mathbb{A}$ in the form

$\displaystyle\mathbb{A}=\mathbb{A}_{\geqslant l}\oplus\mathbb{A}_{<l},\qquad
A_{\geqslant l}:=p^{l}\mathbb{C}[p],\qquad
A_{<l}:=p^{l-1}\mathbb{C}\big{[}p^{-1}\big{]}.$

The related projections are defined by

$\displaystyle\left[\sum_{i}a_{i}p^{i}\right]_{\geqslant l}:=\sum_{i\geqslant
l}a_{i}p^{i}\qquad\text{and}\qquad\left[\sum_{i}a_{i}p^{i}\right]_{<l}:=\sum_{i<l}a_{i}p^{i}.$

The subsets $\mathbb{A}_{\geqslant l}$ and $\mathbb{A}_{<l}$ are subalgebras
only for $l=0$ or $l=1$. As result, the projections on these subalgebras solve
the Rota–Baxter identity (3.1) or (2.5a) with the weight $\kappa=1$, that is
for

$\displaystyle\mathcal{P}=[\cdot]_{\geqslant
l}\qquad\text{or}\qquad\mathcal{P}=[\cdot]_{<l}\quad\text{if}\quad l=0,1.$
(4.1)

Remember that $[\cdot]_{\geqslant l}+[\cdot]_{<l}=1$.

We will look now for solutions of the identity (3.4), where we take the
derivative:

$\displaystyle\partial:=p^{r}\partial_{p},\qquad r\in\mathbb{Z}.$

Then, projections

$\displaystyle\mathcal{R}=[\cdot]_{\geqslant
k-r}\qquad\text{or}\qquad\mathcal{R}=[\cdot]_{<k-r}$ (4.2)

solve the identity (3.4) or (2.5d) with $\kappa=1$ if

  1. 1)

$r=0$ and $k=0$;

  2. 2)

$r\in\mathbb{Z}$ and $k=1$ or $k=2$;

  3. 3)

$r=2$ and $k=3$.

Notice that the above solutions coincide with the $r$-matrices from [3, 32]
with respect to the Poisson bracket defined by

$\displaystyle\\{,\\}_{r}:=p^{r}\partial_{p}\wedge\partial_{x}.$ (4.3)

######  Proposition 4.1.

  * •

All combinations of the above operators $\mathcal{P}$ and $\mathcal{R}$, (4.1)
and (4.2), that satisfy the identities (2.5a) and (2.5c) also satisfy the
remaining identities (2.5b) and (2.5d).

  * •

For

$\displaystyle\mathcal{P}=[\cdot]_{\geqslant
l},\qquad\mathcal{R}=[\cdot]_{\geqslant
k-r}\qquad\text{and}\qquad\partial=p^{r}\partial_{p}$ (4.4a) or
$\displaystyle\mathcal{P}=[\cdot]_{<l},\qquad\mathcal{R}=[\cdot]_{<k-r}\qquad\text{and}\qquad\partial=p^{r}\partial_{p}$
(4.4b)

in addition to the identities (2.5) the constraint (2.6) is also satisfied if

    1. $1)$

$k=0$: $l=0$ and $r=0$;

    2. $2)$

$k=1$: $l=0$ and $r\in\mathbb{Z}$ or $l=1$ and $r=1$;

    3. $3)$

$k=2$: $l=1$ and $r\in\mathbb{Z}$ or $l=0$ and $r=1$;

    4. $4)$

$k=3$: $l=1$ and $r=2$.

  * •

The respective hierarchies (2.1) for the solutions (4.4) turn out to be
independent of $r$. That is, for (4.4a) we get

$\displaystyle\Psi_{t^{\infty}_{n}}=\left[\frac{(L^{n})_{p}}{\\{L,M\\}_{0}}\right]_{\geqslant
l}\Psi_{x}-\left[\frac{(L^{n})_{x}}{\\{L,M\\}_{0}}\right]_{\geqslant
k}\Psi_{p},$ (4.5a) where $\Psi=(L,M)^{\rm T}$ for the Lax function
$L\in\mathbb{A}^{\infty}$ and $M$ being associated Orlov function. For (4.4b)
we get
$\displaystyle\Psi_{t^{0}_{n}}=\left[\frac{(\tilde{L}^{n})_{p}}{\\{\tilde{L},\tilde{M}\\}_{0}}\right]_{<l}\Psi_{x}-\left[\frac{(\tilde{L}^{n})_{x}}{\\{\tilde{L},\tilde{M}\\}_{0}}\right]_{<k}\Psi_{p},$
(4.5b)

where $\Psi=(\tilde{L},\tilde{M})^{\rm T}$ for the Lax function
$\tilde{L}\in\mathbb{A}^{0}$ and $\tilde{M}$ the associated Orlov function.

  * •

The solutions (4.4a) and (4.4b) as well as the Lax hierarchies (4.5a) and
(4.5b) are mutually equivalent through the transformation:

$\displaystyle p^{\prime}=p^{-1}\qquad\text{with}\quad l^{\prime}=1-l,\quad
k^{\prime}=3-k,\quad r^{\prime}=2-r.$ (4.6)

######  Proof 4.2.

The first two points are consequence of straightforward verification. To show
the third point consider the solutions (4.4a) and let
$L\in\mathbb{A}^{\infty}$. Then, the generating functions (2.2) are

$\displaystyle
A_{n}=\left[\frac{p^{r}(L^{n})_{p}}{\\{L,M\\}_{r}}\right]_{\geqslant
l}=\left[\frac{(L^{n})_{p}}{\\{L,M\\}_{0}}\right]_{\geqslant l}$

and

$\displaystyle B_{n}=\left[\frac{(L^{n})_{x}}{\\{L,M\\}_{r}}\right]_{\geqslant
k-r}=p^{-r}\left[\frac{(L^{n})_{x}}{\\{L,M\\}_{0}}\right]_{\geqslant k},$

where the Poisson bracket is defined by (4.3). Hence, the related hierarchy
(2.1) takes the form (4.5a), since the factors $p^{-r}$ and $p^{r}$ cancel
out. For the solutions (4.4b) and $\tilde{L}\in\mathbb{A}^{0}$ the reasoning
is similar. The last point follows from simple verification.

For fixed parameters $k,r$ and $l$ we can consider the Lax functions $L$ and
$\tilde{L}$ to be analytic extensions of some function outside and inside,
respectively, unit circle on the complex plane. Thus $L$ and $\tilde{L}$ are
defined near $\infty$ and $0$, respectively, where they have poles. We can
choose these poles to be simple. As result we can consider two families of
related hierarchies together and postulate
$\Psi=(L,M,\tilde{L},\tilde{M})^{\rm T}$ in (4.5). This is, inter alia,
consequence of the fact that the solutions (4.4a) and (4.4b) are mutually
associated through the decomposition of unity:

$\displaystyle 1=[\cdot]_{\geqslant l}+[\cdot]_{<l}=[\cdot]_{\geqslant
k-r}+[\cdot]_{<k-r}.$

The compatibility equations (2.4) for the hierarchies (4.5) takes the
following, independent of $r$, form:

$\displaystyle\bigl{(}A_{n}^{\lambda}\bigr{)}_{t_{m}^{\mu}}-\bigl{(}A_{m}^{\mu}\bigr{)}_{t_{n}^{\lambda}}+\bigl{\langle}A_{n}^{\lambda},A_{m}^{\mu}\bigr{\rangle}_{x}+B_{m}^{\mu}\bigl{(}A_{n}^{\lambda}\bigr{)}_{p}-B_{n}^{\lambda}\bigl{(}A_{m}^{\mu}\bigr{)}_{p}=0$
(4.7a) and
$\displaystyle\bigl{(}B_{n}^{\lambda}\bigr{)}_{t_{m}^{\mu}}-\bigl{(}B_{m}^{\mu}\bigr{)}_{t_{n}^{\lambda}}+A_{n}^{\lambda}\bigl{(}B_{m}^{\mu}\bigr{)}_{x}-A_{m}^{\mu}\bigl{(}B_{n}^{\lambda}\bigr{)}_{x}-\bigl{\langle}B_{n}^{\lambda},B_{m}^{\mu}\bigr{\rangle}_{p}=0,$
(4.7b)

where $m,n\in\mathbb{N}$; $\mu,\lambda=\infty,0$ and $\langle
f,g\rangle_{\xi}:=fg_{\xi}-gf_{\xi}$ for $\xi=x,p$. The related generating
functions are defined as

$\displaystyle
A_{n}^{\infty}:=\left[\frac{(L^{n})_{p}}{\\{L,M\\}_{0}}\right]_{\geqslant
l},\qquad
B_{n}^{\infty}:=\left[\frac{(L^{n})_{x}}{\\{L,M\\}_{0}}\right]_{\geqslant k},$

and

$\displaystyle
A_{n}^{0}:=\left[\frac{(\tilde{L}^{n})_{p}}{\\{\tilde{L},\tilde{M}\\}_{0}}\right]_{<l},\qquad
B_{n}^{0}:=\left[\frac{(\tilde{L}^{n})_{x}}{\\{\tilde{L},\tilde{M}\\}_{0}}\right]_{<k}.$

Notice that if we treat $\Psi$ as a common eigenfunction, then the pairs of
equations from (4.5) give the Lax pairs for the respective systems from (4.7).

Let us consider reductions from Lemma 2.5. First for constraints of the type
(2.7) we take

$\displaystyle\\{L,M\\}_{r}=L^{r}\qquad\text{and}\qquad\\{\tilde{L},\tilde{M}\\}_{r}=\tilde{L}^{2-r}.$
(4.8)

We assume that they are consistent with respect to the appropriate choice of
Lax and Orlov functions, see the forthcoming examples. Then, in accordance
with equations (2.7) and (2.8) the hierarchies (4.5) reduce to

$\displaystyle L_{t^{\infty}_{n}}=\frac{n}{n-r}\bigl{\\{}[L^{n-r}]_{\geqslant
k-r},L\bigr{\\}}_{r}\qquad\text{and}\qquad\tilde{L}_{t^{0}_{n}}=\frac{n}{n-2+r}\bigl{\\{}[\tilde{L}^{n+2-r}]_{<k-r},\tilde{L}\bigr{\\}}_{r}.$

So we obtain the standard Poisson hierarchies of dispersionless systems, see
[3, 32] and references therein. For the next reduction from Lemma 2.5 we
postulate that

$\displaystyle L=p\qquad\text{and}\qquad\tilde{L}=p^{-1}.$

Then, the hierarchies (4.5) reduce to

$\displaystyle
M_{t^{\infty}_{n}}=n\bigl{[}p^{n-1}M_{x}^{-1}\bigr{]}_{\geqslant
l}M_{x}\qquad\text{and}\qquad\tilde{M}_{t^{0}_{n}}=n\bigl{[}p^{1-n}\tilde{M}_{x}^{-1}\bigr{]}_{<l}\tilde{M}_{x}.$

Taking $G:=M_{x}^{-1}$ and $\tilde{G}:=\tilde{M}_{x}^{-1}$ the above
hierarchies can be rewritten in the form

$\displaystyle
G_{t^{\infty}_{n}}=n\bigl{\langle}\bigl{[}p^{n-1}G\bigr{]}_{\geqslant
l},G\bigr{\rangle}_{x}\qquad\text{and}\qquad\tilde{G}_{t^{0}_{n}}=n\bigl{\langle}\bigl{[}p^{1-n}\tilde{G}_{x}\bigr{]}_{<l},\tilde{G}\bigr{\rangle}_{x}.$

These are the universal hierarchies considered in [21, 22, 23].

For appropriate choice of Lax functions and associated Orlov functions the Lax
hierarchies (4.5), in principle, yield construction of $(1+1)$-dimensional
integrable infinite-field (chain) systems, while the compatibility conditions
(4.7) provide finite-field systems that include $(2+1)$-dimensional integrable
equations of Manakov–Santini type, which are of our interest. Although it
would be fairly easy, we will not consider in this work finite-field
reductions of Lax hierarchies (4.5).

###  4.1 The case of $\boldsymbol{k=l=r=0}$

The Lax function [2, 33, 34] has the form

$\displaystyle
L=p+u(x)p^{-1}+u_{2}(x)p^{-2}+u_{3}(x)p^{-3}+\cdots\in\mathbb{A}^{\infty}$

and the associated Orlov function [33, 34] is

$\displaystyle M=M_{0}+x+v(x)L^{-1}+v_{2}(x)L^{-2}+\cdots,$

where $M_{0}$ is the part of $M$ that explicitly depends on times
$t^{\infty}_{n}$ and commutes with $L$, that is $\\{M_{0},L\\}_{0}\equiv 0$.
This means that the choice of $M_{0}$ does not influence the construction of
related systems from (4.7). In this case, consistent Lax hierarchy for
$\tilde{L}\in\mathbb{A}^{0}$ does not exist.

Then, we calculate

$\displaystyle\\{L,M\\}_{0}=1+v_{x}p^{-1}+\bigl{(}(v_{2})_{x}-u\bigr{)}p^{-2}+\cdots$

and the first generating functions

$\displaystyle A^{\infty}_{1}=1,\qquad$ $\displaystyle
A^{\infty}_{2}=p-v_{x},\qquad$ $\displaystyle
A^{\infty}_{3}=p^{2}-v_{x}p+2u+v_{x}^{2}-(v_{2})_{x},$ $\displaystyle
B^{\infty}_{1}=0,\qquad$ $\displaystyle B^{\infty}_{2}=u_{x},\qquad$
$\displaystyle B^{\infty}_{3}=u_{x}p+(u_{2})_{x}-u_{x}v_{x}.$

From the generalized zero-curvature equations (4.7) for $n=1$, $m=2$ and
$\lambda=\mu=\infty$ we get

$\displaystyle u_{t_{1}}=u_{x},\qquad v_{t_{1}}=v_{x},$

where $t_{1}\equiv t^{\infty}_{1}$. For $n=1$, $m=3$ we get

$\displaystyle(u_{1})_{t_{1}}=(u_{1})_{x},\qquad(v_{2})_{t_{1}}=(v_{2})_{x}.$

Which in fact means that we can identify the time $t_{1}$ with $x$. First
nontrivial equations are for $n=2$ and $m=3$. Let $t\equiv t_{3}$ and $y\equiv
t_{2}$. Hence, we obtain two compatibility conditions:

$\displaystyle(v_{2})_{x}=u+v_{y}+v_{x}^{2},\qquad(u_{2})_{x}=u_{y}+u_{x}v_{x}$

and the celebrated Manakov–Santini system [15]

$\displaystyle\begin{split}&v_{xt}=v_{yy}+v_{x}v_{xy}+uv_{xx}-v_{y}v_{xx},\\\
&u_{xt}=u_{yy}+u_{x}^{2}+u_{xy}v_{x}+uu_{xx}-u_{xx}v_{y}.\end{split}$ (4.9)

The related Lax pair which follows from (4.5a) is

$\displaystyle\partial_{y}\Psi=\left[(p-v_{x})\partial_{x}-u_{x}\partial_{p}\right]\Psi,$
$\displaystyle\partial_{t}\Psi=\left[\bigl{(}p^{2}-v_{x}p+u-v_{y}\bigr{)}\partial_{x}-\bigl{(}u_{x}p+u_{y}\bigr{)}\partial_{p}\right]\Psi.$

The reduction given by the condition $\\{L,M\\}_{0}=1$ means that $v=0$. Thus,
for $v=0$ from (4.9) we obtain the dispersionless KP equation

$\displaystyle u_{xt}=u_{yy}+u_{x}^{2}+uu_{xx}.$

The second possible reduction is given by the constraint: $L=p$, from which it
follows that $u=0$ and the system (4.9) reduces to the Pavlov equation [24]
(see also [11, 22, 30])

$\displaystyle v_{xt}=v_{yy}+v_{x}v_{xy}-v_{y}v_{xx}.$

###  4.2 The case: $\boldsymbol{k=1}$, $\boldsymbol{l=0}$ and
$\boldsymbol{r\in\mathbb{Z}}$

The Lax and associated Orlov functions [2, 19, 20, 34] for
$p\rightarrow\infty$ are given by

$\displaystyle
L=p+u(x)+u_{1}(x)p^{-1}+u_{2}(x)p^{-2}+\cdots\in\mathbb{A}^{\infty}$

and

$\displaystyle
M=M_{0}+x+\bigl{(}\partial_{x}^{-1}u(x)+w(x)\bigr{)}L^{-1}+w_{2}(x)L^{-2}+\cdots,$

where as before $\\{M_{0},L\\}_{0}\equiv 0$ and we made some modification for
convenience of further calculations. Then

$\displaystyle\\{L,M\\}_{0}=1+(u+w_{x})p^{-1}+\bigl{(}(w_{2})_{x}-u_{1}-u^{2}-uw_{x}\bigr{)}p^{-2}+\cdots$

and

$\displaystyle A^{\infty}_{1}=1,\qquad$ $\displaystyle
A^{\infty}_{2}=p-w_{x},\qquad$ $\displaystyle
A^{\infty}_{3}=p^{2}+(u-w_{x})p+u^{2}+uw_{x}+w_{x}^{2}+2u_{1}-(w_{2})_{x},$
$\displaystyle B^{\infty}_{1}=0,\qquad$ $\displaystyle
B^{\infty}_{2}=u_{x}p,\qquad$ $\displaystyle
B^{\infty}_{3}=u_{x}p^{2}+\bigl{(}uu_{x}-u_{x}w_{x}+(u_{1})_{x}\bigr{)}p.$

From the generalized zero-curvature equations (4.7) for $n=1$, $m=2$ and
$\lambda=\mu=\infty$ we get

$\displaystyle u_{t_{1}}=u_{x},\qquad w_{t_{1}}=w_{x},$

where $t_{1}\equiv t^{\infty}_{1}$. For $n=1$, $m=3$ we get

$\displaystyle(u_{1})_{t_{1}}=(u_{1})_{x},\qquad(w_{2})_{t_{1}}=(w_{2})_{x}.$

Which means that we can identify the time $t_{1}$ with $x$. First nontrivial
equations are for $n=2$, $m=3$. Let $t\equiv t^{\infty}_{3}$ and $y\equiv
t^{\infty}_{2}$. After some simplifications, we obtain two compatibility
conditions

$\displaystyle(u_{1})_{x}=u_{y}+u_{x}w_{x},\qquad(w_{2})_{xx}=2u_{y}+w_{xy}+2uu_{x}+3u_{x}w_{x}+2uw_{xx}+2w_{x}w_{xx}$

and the new $(2+1)$-dimensional integrable system

$\displaystyle\begin{split}&u_{xt}=u_{yy}+u_{x}u_{y}+u_{x}^{2}w_{x}+uu_{xy}+u_{xy}w_{x}+u_{xx}a,\\\
&w_{xt}=uw_{xy}+u_{y}w_{x}+w_{x}w_{xy}+aw_{xx}-a_{y},\end{split}$ (4.10)

where

$\displaystyle a_{x}=u_{x}w_{x}-w_{xy}.$

The Lax pair for the system (4.10) is given by

$\displaystyle\partial_{y}\Psi=\bigl{[}(p-w_{x})\partial_{x}-u_{x}p\partial_{p}\bigr{]}\Psi,$
$\displaystyle\partial_{t}\Psi=\bigl{[}\bigl{(}p^{2}+(u-w_{x})p-w_{y}-\partial_{x}^{-1}uw_{xx}\bigr{)}\partial_{x}-\bigl{(}u_{x}p^{2}+(uu_{x}+u_{y})p\bigr{)}\partial_{p}\bigr{]}\Psi.$

Consider first reduction (4.8) given by the condition

$\displaystyle\\{L,M\\}_{r}=L^{r},\qquad r\in\mathbb{Z},$ (4.11)

which is consistent since the order at $\infty$ of both sides of the equality
is the same. Hence,

$\displaystyle\\{L,M\\}_{0}=p^{-r}L^{r}=1+rup^{-1}+\left(ru_{1}-\frac{1}{2}r(1-r)u^{2}\right)p^{-2}+\cdots,$

and we obtain the constraint

$\displaystyle w_{x}=(r-1)u$

by means of which (4.10) gives the $r$-th dispersionless modified KP equation
[2, 4]:

$\displaystyle
u_{t}=\frac{1}{2}(r-1)u^{2}u_{x}+ruu_{y}+\partial_{x}^{-1}u_{yy}+(1-r)u_{x}\partial_{x}^{-1}u_{y},$

which for $r=0$ gives the standard dispersionless modified KP equation. The
second reduction is given by the constraint: $L=p$, from which it follows that
$u=0$ and in this case the system (4.10) reduces again to the Pavlov equation

$\displaystyle w_{xt}=w_{yy}+w_{x}w_{xy}-w_{y}w_{xx}.$

The Lax and the associated Orlov functions for $p\rightarrow 0$ are

$\displaystyle\tilde{L}=v(x)p^{-1}+v_{0}(x)+v_{1}(x)p+v_{2}(x)p^{2}+\cdots\in\mathbb{A}^{0}$

and

$\displaystyle\tilde{M}=\tilde{M}_{0}+m(x)+m_{1}(x)\tilde{L}^{-1}+m_{2}(x)\tilde{L}^{-2}+\cdots,$

where $\\{\tilde{M}_{0},\tilde{L}\\}_{0}\equiv 0$. We have

$\displaystyle\\{\tilde{L},\tilde{M}\\}_{0}=-vm_{x}p^{-2}-(m_{1})_{x}p^{-1}+\cdots$

and

$\displaystyle A^{0}_{1}=0,\qquad$ $\displaystyle
A^{0}_{2}=\frac{v}{m_{x}}p^{-1},\qquad$ $\displaystyle
A^{0}_{3}=\frac{v^{2}}{m_{x}}p^{-2}+\frac{2vv_{0}m_{x}-v(m_{1})_{x}}{m_{x}^{2}}p^{-1},$
$\displaystyle B^{0}_{1}=0,\qquad$ $\displaystyle
B^{0}_{2}=-\frac{v_{x}}{m_{x}},\qquad$ $\displaystyle
B^{0}_{3}=-\frac{vv_{x}}{m_{x}}p^{-1}+\left(\frac{(m_{1})_{x}}{m_{x}^{2}}-\frac{2v_{0}}{m_{x}}\right)v_{x}-\frac{v(v_{0})_{x}}{m_{x}}.$

From the generalized zero-curvature equations (4.7) for $n=2$, $m=3$ and
$\lambda=\mu=0$ we get the constraints

$\displaystyle(v_{0})_{x}=\frac{v_{z}m_{x}}{v},\qquad\left(\frac{(m_{1})_{x}}{m_{x}}\right)_{x}=m_{xz}+\frac{v_{z}}{v}m_{x},$

where $z\equiv t^{0}_{2}$, and the system

$\displaystyle\left(\frac{v}{m_{x}}\right)_{\tau}=-\left[\frac{v}{m_{x}}\partial_{x}^{-1}v\left(\frac{m_{x}}{v}\right)_{z}\right]_{z},\qquad\left(\frac{v_{x}}{m_{x}}\right)_{\tau}=v_{zz}-\left[\frac{v_{x}}{m_{x}}\partial_{x}^{-1}v\left(\frac{m_{x}}{v}\right)_{z}\right]_{z},$
(4.12)

where $\tau\equiv t^{0}_{3}$. The Lax pair for (4.12) is given by

$\displaystyle\partial_{\tau}\Psi=\left[\frac{v}{m_{x}}p^{-1}\partial_{x}+\frac{v_{x}}{m_{x}}\partial_{p}\right]\Psi,$
$\displaystyle\partial_{z}\Psi=\left[\left(\frac{v^{2}}{m_{x}}p^{-2}-\frac{v}{m_{x}}\partial_{x}^{-1}v\left(\frac{m_{x}}{v}\right)_{\\!z}p^{-1}\right)\partial_{x}+\left(\frac{vv_{x}}{m_{x}}p^{-1}+v_{z}-\frac{v_{x}}{m_{x}}\partial_{x}^{-1}v\left(\frac{m_{x}}{v}\right)_{\\!z}\right)\partial_{p}\right]\Psi.$

Consider first reduction given by the condition

$\displaystyle\\{\tilde{L},\tilde{M}\\}_{r}=\tilde{L}^{2-r},\qquad
r\in\mathbb{Z},$ (4.13)

which is consistent since order at $0$ of both sides is the same. Thus,

$\displaystyle\\{\tilde{L},\tilde{M}\\}_{0}=p^{-r}\tilde{L}^{2-r}=v^{2-r}p^{-2}+(2-r)v^{1-r}v_{0}p^{-1}+\cdots$

and we obtain the constraint

$\displaystyle m_{x}=-v^{1-r},$

by means of which (4.12) gives the $r$-th dispersionless Harry–Dym equation
[2, 4]:

$\displaystyle
v_{\tau}=-v^{1-r}\bigl{[}v^{r}\partial_{x}^{-1}v^{-r}v_{z}\bigr{]}_{z}.$

The second reduction is given by the constraint: $\tilde{L}=p^{-1}$, from
which it follows that $v=1$ and the system (4.12) reduces to [22, 30]

$\displaystyle m_{x\tau}=m_{x}m_{zz}-m_{xz}m_{z}.$

Now, we will consider the mixed case. For $\lambda=\infty$, $n=1$ and $\mu=0$
with $m=2$ or $m=3$ in (4.7) we obtain

$\displaystyle v_{t_{1}}=v_{x},\qquad m_{t_{1}}=m_{x},$

which holds automatically as $t_{1}\equiv x$. From the zero-curvature
equations (4.7), for $\lambda=\infty$, $n=2$ and $\mu=0$, $m=2$, we obtain the
following compatibility equations

$\displaystyle u_{z}=\frac{v_{x}}{m_{x}},\qquad
w_{xz}=-v\frac{m_{xx}}{m_{x}^{2}}$

and

$\displaystyle\left(\frac{v}{m_{x}}\right)_{y}=\frac{v}{m_{x}}(u_{x}+w_{xx})-w_{x}\left(\frac{v}{m_{x}}\right)_{x},\qquad\left(\frac{v_{x}}{m_{x}}\right)_{y}=\frac{(u_{x}v)_{x}}{m_{x}}-w_{x}\left(\frac{v_{x}}{m_{x}}\right)_{x}.$

The above equations can be equivalently rewritten in the following form, for
which we have the conditions

$\displaystyle u_{x}=(\log v)_{y}-(\log v)_{x}\frac{m_{y}}{m_{x}},\qquad
w_{x}=-\frac{m_{y}}{m_{x}}$

and the system

$\displaystyle m_{xx}v=m_{x}m_{yz}-m_{y}m_{xz},\qquad v_{xx}=(\log
v)_{yz}m_{x}-(\log v)_{xz}m_{y},$ (4.14)

which is a system of Manakov–Santini type recently obtained in [7]. The
related Lax pair is

$\displaystyle\partial_{y}\Psi=\left[(p+w)\partial_{x}-u_{x}p\partial_{p}\right]\Psi,\qquad\partial_{z}\Psi=\left[\frac{v}{m_{x}}p^{-1}\partial_{x}+\frac{v_{x}}{m_{x}}\partial_{p}\right]\Psi.$

Consider now the reduction by means of (4.11) and (4.13), that is we have
constraints in the form

$\displaystyle w_{x}=(r-1)u,\qquad m_{x}=-v^{1-r},\qquad r\in\mathbb{Z}.$

Hence, from (4.14) we obtain the $r$-th dispersionless Toda system [4, 19]:

$\displaystyle u_{z}+v^{r-1}v_{x}=0,\qquad v_{y}=u_{x}v+(1-r)uv_{x},$

which for $r=1$ gives the Boyer–Finley equation being $(2+1)$-dimensional
version of dispersionless Toda equation. On the other hand from the reduction:
$L=p$, $\tilde{L}=p^{-1}$ we have the constraints: $u=0$, $v=1$, from which we
get the equation [22, 24]

$\displaystyle m_{xx}=m_{x}m_{yz}-m_{y}m_{xz}.$

The remaining not considered cases from Proposition 4.1 are equivalent to the
above examples through the transformation (4.6) or there is no consistent
hierarchies for them on the level of equations (4.5).

### Acknowledgement

The very preliminary results related to this work were obtained in 2010 during
the visit of the author to Institute of Mathematics in Academia Sinica in
Taipei. The author wishes to express his thanks to Professor Jen-Hsu Chang for
valuable discussions and Professor Jyh-Hao Lee for warm hospitality.

## References

  * [1] Baxter G., An analytic problem whose solution follows from a simple algebraic identity, Pacific J. Math. 10 (1960), 731–742.
  * [2] Błaszak M., Classical $R$-matrices on Poisson algebras and related dispersionless systems, Phys. Lett A 297 (2002), 191–195.
  * [3] Błaszak M., Szablikowski B.M., Classical $R$-matrix theory of dispersionless systems. I. $(1+1)$-dimensional theory, J. Phys. A: Math. Gen. 35 (2002), 10325–10344, nlin.SI/0211008.
  * [4] Błaszak M., Szablikowski B.M., Classical $R$-matrix theory of dispersionless systems. II. $(2+1)$ dimension theory, J. Phys. A: Math. Gen. 35 (2002), 10345–10364, nlin.SI/0211018.
  * [5] Błaszak M., Szablikowski B.M., Classical $R$-matrix theory for bi-Hamiltonian field systems, J. Phys. A: Math. Theor. 42 (2009), 404002, 35 pages, arXiv:0902.1511.
  * [6] Bogdanov L.V., On a class of multidimensional integrable hierarchies and their reductions, Theoret. and Math. Phys. 160 (2009), 887–893, arXiv:0810.2397.
  * [7] Bogdanov L.V., Non-Hamiltonian generalizations of the dispersionless 2DTL hierarchy, J. Phys. A: Math. Theor. 43 (2010), 434008, 8 pages, arXiv:1003.0287.
  * [8] Bogdanov L.V., On a class of reductions of the Manakov–Santini hierarchy connected with the interpolating system, J. Phys. A: Math. Theor. 43 (2010), 115206, 11 pages, arXiv:0910.4004.
  * [9] Bogdanov L.V., Chang J.-H., Chen Y.-T., Generalized dKP: Manakov–Santini hierarchy and its waterbag reduction, arXiv:0810.0556.
  * [10] Chang J.-H., Chen Y.-T., Hodograph solutions for the Manakov–Santini equation, J. Math. Phys. 51 (2010), 042701, 18 pages, arXiv:0904.4595.
  * [11] Dunajski M., A class of Einstein–Weyl spaces associated to an integrable system of hydrodynamic type, J. Geom. Phys. 51 (2004), 126–137, nlin.SI/0311024.
  * [12] Guo L., What is $\ldots$ a Rota–Baxter algebra?, Notices Amer. Math. Soc. 56 (2009), 1436–1437.
  * [13] Guo L., An introduction to Rota–Baxter algebra, Surveys of Modern Mathematics, Vol. 4, International Press, Somerville, MA, 2012.
  * [14] Manakov S.V., Santini P.M., Inverse scattering problem for vector fields and the Cauchy problem for the heavenly equation, Phys. Lett. A 359 (2006), 613–619, nlin.SI/0604024.
  * [15] Manakov S.V., Santini P.M., A hierarchy of integrable partial differential equations in dimension $2+1$, associated with one-parameter families of vector fields, Theoret. and Math. Phys. 152 (2007), 1004–1011.
  * [16] Manakov S.V., Santini P.M., On the solutions of the dKP equation: the nonlinear Riemann Hilbert problem, longtime behaviour, implicit solutions and wave breaking, J. Phys. A: Math. Theor. 41 (2008), 055204, 23 pages, arXiv:0707.1802.
  * [17] Manakov S.V., Santini P.M., The dispersionless 2D Toda equation: dressing, Cauchy problem, longtime behaviour, implicit solutions and wave breaking, J. Phys. A: Math. Theor. 42 (2009), 095203, 16 pages, arXiv:0810.4676.
  * [18] Manakov S.V., Santini P.M., On the solutions of the second heavenly and Pavlov equations, J. Phys. A: Math. Theor. 42 (2009), 404013, 11 pages, arXiv:0812.3323.
  * [19] Mañas M., On the $r$th dispersionless Toda hierarchy: factorization problem, additional symmetries and some solutions, J. Phys. A: Math. Gen. 37 (2004), 9195–9224, nlin.SI/0404022.
  * [20] Mañas M., $S$-functions, reductions and hodograph solutions of the $r$th dispersionless modified KP and Dym hierarchies, J. Phys. A: Math. Gen. 37 (2004), 11191–11221, nlin.SI/0405028.
  * [21] Martínez Alonso L., Shabat A.B., Energy-dependent potentials revisited: a universal hierarchy of hydrodynamic type, Phys. Lett. A 299 (2002), 359–365, nlin.SI/0202008.
  * [22] Martínez Alonso L., Shabat A.B., Towards a theory of differential constraints of a hydrodynamic hierarchy, J. Nonlinear Math. Phys. 10 (2003), 229–242, nlin.SI/0310036.
  * [23] Martínez Alonso L., Shabat A.B., Hydrodynamic reductions and solutions of the universal hierarchy, Theoret. and Math. Phys. 140 (2004), 1073–1085, nlin.SI/0312043.
  * [24] Pavlov M.V., Integrable hydrodynamic chains, J. Math. Phys. 44 (2003), 4134–4156, nlin.SI/0301010.
  * [25] Pavlov M.V., Chang J.-H., Chen Y.-T., Integrability of the Manakov–Santini hierarchy, arXiv:0910.2400.
  * [26] Rota G.C., Baxter algebras and combinatorial identities. I, Bull. Amer. Math. Soc. 75 (1969), 325–329.
  * [27] Rota G.C., Baxter algebras and combinatorial identities. II, Bull. Amer. Math. Soc. 75 (1969), 330–334.
  * [28] Semenov-Tian-Shansky M.A., Integrable systems and factorization problems, in Factorization and Integrable Systems (Faro, 2000), Oper. Theory Adv. Appl., Vol. 141, Birkhäuser, Basel, 2003, 155–218, nlin.SI/0209057.
  * [29] Semenov-Tyan-Shanskii M.A., What is a classical $r$-matrix?, Funct. Anal. Appl. 17 (1983), 259–272.
  * [30] Sergyeyev A., Szablikowski B.M., Central extensions of cotangent universal hierarchy: $(2+1)$-dimensional bi-Hamiltonian systems, Phys. Lett. A 372 (2008), 7016–7023, arXiv:0807.1294.
  * [31] Szablikowski B.M., Classical $r$-matrix like approach to Frobenius manifolds, WDVV equations and flat metrics, J. Phys. A: Math. Theor. 48 (2015), 315203, 47 pages, arXiv:1304.2075.
  * [32] Szablikowski B.M., Błaszak M., Meromorphic Lax representations of $(1+1)$-dimensional multi-Hamiltonian dispersionless systems, J. Math. Phys. 47 (2006), 092701, 23 pages, nlin.SI/0510068.
  * [33] Takasaki K., Takebe T., ${\rm SDif\/f}(2)$ KP hierarchy, Internat. J. Modern Phys. A 7 (1992), 889–922, hep-th/9112046.
  * [34] Takasaki K., Takebe T., Integrable hierarchies and dispersionless limit, Rev. Math. Phys. 7 (1995), 743–808, hep-th/9405096.


# Meta-modeling strategy for data-driven forecasting

Dominic J. Skinner
Department of Mathematics
Massachusetts Institute of Technology
Cambridge, MA 02139-4307
<EMAIL_ADDRESS>&Romit Maulik
Argonne Leadership Computing Facility
Argonne National Laboratory
Lemont, IL 60439
<EMAIL_ADDRESS>
###### Abstract

Accurately forecasting the weather is a key requirement for climate change
mitigation. Data-driven methods offer the ability to make more accurate
forecasts, but lack interpretability and can be expensive to train and deploy
if models are not carefully developed. Here, we make use of two historical
climate data sets and tools from machine learning, to accurately predict
temperature fields. Furthermore, we are able to use low fidelity models that
are cheap to train and evaluate, to selectively avoid expensive high fidelity
function evaluations, as well as uncover seasonal variations in predictive
power. This allows for an adaptive training strategy for computationally
efficient geophysical emulation.

##  1 Introduction

While numerical weather forecasting dates back over a century [10], climate
change is expected to significantly alter the predictability of the
atmosphere, increasing the error of weather forecasts [18]. Simultaneously,
climate mitigation requires more accurate forecasts of weather events for
power grid optimization [1], and extreme event prediction [20], from flooding
[9], to heatwaves [6]. Even the largest climate simulations have
$O(10\,\text{km})$ separation between grid points, meaning dynamics of smaller
scales can not be explicitly resolved [14, 17, 15], and often ad-hoc closure
conditions are postulated to account for the unresolved dynamics. Machine
learning has the potential to find data-driven closure conditions [16] and
parameterize sub-grid scale modeling [3, 7]. Convolutional neural networks can
forecast weather [4], deep learning can predict extreme weather events [5],
and neural network architectures can be optimized automatically to enhance the
quality of forecasts [11, 2]. However, these models can become expensive to
train and deploy, especially when ensemble predictions are required [8], and
often their predictions lack interpretability. Here, we use a combination of
proper orthogonal decomposition (POD), and long-short term memory (LSTM)
reccurent neural networks to create forecasts for two real-world data sets. In
addition to the LSTM network, which is our high-fidelity (HF) model, we
introduce simpler or low-fidelity (LF) models, which do not have the same
predictive power, but are faster to deploy. By studying the prediction
differences of the LF models, we are able to selectively avoid HF function
evaluations, as well as uncovering seasonal variations in prediction accuracy.
With the proliferation of data driven methods for climate forecasting, the
approach introduced here, could be used to reduce computational resources,
which is of particular relevance for ensemble forecasting [8], as well as to
identify uncertainties in neural network predictions.

##  2 Proper Orthogonal Decomposition

To first reduce the dimensionality of the problem, we project spatially
resolved fields, such as temperature, onto a set of principal modes which
capture the salient features and then track the evolution of these modes in
time. We use the technique of Proper Orthogonal Decomposition (POD) in order
to find the dominant modes [21, 12, 13]. In short, suppose we are given
spatial snapshots of the data at various times,
$\bm{\theta}_{1},\dots,\bm{\theta}_{T}$. We take a truncated set of
orthonormal basis vectors $\bm{v_{1}},\dots,\bm{v}_{M}$, that approximates the
spatial snapshots optimally with respect to the $L_{2}$ norm, i.e. minimizing
$\sum_{i}^{T}||\bm{\theta_{i}}-\bm{\phi}_{i}||^{2}$ over
$\bm{\phi}_{i}\in\text{Span}\\{\bm{v}_{1},\dots,\bm{v}_{M}\\}$. Defining the
snapshot matrix, $\bm{S}=[\bm{\theta}_{1}|\cdots|\bm{\theta}_{T}]$, the
optimal basis is given by the $M$ eigenvectors of $\bm{SS}^{T}$, with largest
eigenvalues [21], which is readily found numerically. See appendix A, for
further details.

Figure 1:  Proper orthogonal decomposition of sea surface temperature. (a)
Historical climate data measuring ocean surface temperature weekly from 1981
to 2000. (b) Evolution of the first 3 modal coefficients, together with their
respective basis elements. (c) Reconstructed temperature field from first the
6 modal coefficients against measured temperature field.

##  3 NOAA sea surface temperature

Our first data set is the NOAA Optimum Interpolation SST V2 data set,
containing the sea surface temperature weekly on a 1 degree grid across the
period 1981-2018 111Available at https://www.esrl.noaa.gov/psd/. The first 20
years were used as training data, the rest was reserved for testing. From
this, we build data driven forecasts, predicting the next 5 weeks of
temperature evolution from historical data. We project the system onto the
first 6 POD modes, which is sufficient to approximate the temperature field
(Fig. 1), and examine the evolution of these modal coefficients (Fig. 1a-b).
From predictions of the modal coefficients, we can reconstruct the temperature
field and test the temperature predictions at a sensor located at a specific
coordinate, Fig. 2. The baseline to improve on is the climatology prediction;
the historical average temperature at that specific location for the time of
year.

For our purposes, we will use a bi-directional LSTM network [19] as our high-
fidelity model, taking 5 weeks of historical data as input and making a 5 week
forecast. This model is proto-typical for more complex machine learning models
that could be deployed on data like this [11]. In short, a bi-directional LSTM
is a recurrent neural network, that instead of acting purely sequentially on
the data, has an additional pass backwards in time, Fig. 3c, in practice
improving predictive power. In all cases, we take the first 20 years as
training data, and use the remainder for testing. We compare the LSTM
predictions against the recorded temperature and climatology baseline, Fig. 2,
finding that while both the climatology and LSTM forecast captures the
seasonal trends, the LSTM predictions have lower $L_{2}$ error and cosine norm
closer to 1, as compared to the recorded temperature, than the climatology
prediction for all sensors tested.

Figure 2: Sea surface temperature at 4 randomly chosen sensors (orange points)
for the first 4 years of testing data (2002-2006), with POD-LSTM prediction
and climatology baseline. After training, the LSTM model was applied non-
autoregressively, making 5 weeks of predictions at a time. The spatial
predictions were then reconstructed from the POD basis. Also shown is the
climatology baseline.

Our aim now is to investigate when and where the model breaks down, as well as
identifying where expensive HF function evaluations can be avoided, and
substituted for low-fidelity alternatives. To do so, we introduce two LF
models. The first is linear regression, applied autoregressively to each mode
individually, with an input window of 5 weeks, and making a single prediction.
The best linear model is fitted to the training data for a particular mode,
and then applied to the testing data. The second is a random forest regressor
[22], which takes in all 5 weeks of input, and trains a model to predict each
modal coefficient across the forecast. This works by training 100 decision
trees on the input sequence, and then taking an average of the results as the
output, Fig. 3d.

Figure 3:  Comparison of high and low-fidelity models. (a) Input and
predictions for a representative test example, true modal coefficients are
solid lines, predictions are dashed. (b) Distributions of forecast errors for
the testing data. (c) Architecture of the bi-directional LSTM, in our case
$n=10$. The input information travels forwards through time, and then
backwards again to allow for a more accurate global prediction. (d)
Architecture of the random forest; the input is sent to several decision
trees, and their result is aggregated at the end for a prediction. (e) Total
error of adaptive meta-modeling strategy against the optimal and random
strategy

Unsurprisingly, the LF models perform worse on testing data than the HF one,
Fig. 3(a,b). However, the low-fidelity models are extremely cheap to train and
evaluate, we therefore seek to make use of them whenever possible. Suppose now
that we have a limited computational budget for evaluating a number of
forecasts, and so are forced to use a mixture of HF and LF models. If we made
the choice randomly with probability $p$ for the HF model, we would expect the
average error to be $p\times(\text{HF average error})+(1-p)\times(\text{LF
average error})$. Similarly, we expect the computational cost to be
$p\times(\text{HF cost})+(1-p)\times(\text{LF cost})$. We will use the two LF
models to create a meta-modeling strategy that achieves lower error than the
random strategy for the same computational cost. The strategy is as follows:

  1. 1.

Take input data and perform the random forest and linear forecasts.

  2. 2.

Calculate the difference between the predictions.

  3. 3.

If this difference exceeds some threshold, evaluate the HF model. If not,
stick with the LF model.

This simple strategy using two LF models significantly outperforms the random
strategy, Fig. 3e. Moreover, it performs similarly to the theoretical optimum
strategy, defined as the best possible choice of HF and LF evaluations, given
their values and the true value in advance. This optimum strategy is clearly
not viable in practice, yet the meta-model strategy performs similarly.

##  4 DayMet North America Daily Surface Temperature

We now validate the ideas we have introduced on an alternative data set. We
consider the forecasting problem for the Daymet data set222available at
https://daymet.ornl.gov/, containing the maximum daily temperature field
across North America for 2000-2015, Fig. 4. The first 11 years are used as
training data, the rest is reserved for testing. The temperature over land is
more variable than the sea temperature, so here we forecast 7 days from 7 days
of historical data, and take only the first 4 POD modes.

Figure 4:  (a) Maximum daily temperature field across North America in 2011,
representative seasonal snapshots shown. (b) Temporal evolution of the first
four POD modes of the temperature field across one year, together with their
representative basis elements. (c) Reconstruction of the temperature field
from 4 modes against the true value.

The bi-directional LSTM has the lowest total prediction error, Fig. 5, and
outperforms the climatology baseline for all sensors tested. The linear model
performs almost as well, and both perform better than the random forest, Fig.
5a-b. Beyond overall error analysis, we aim to interrogate where predictions
break down. Averaging the prediction error across 5 years of testing data, we
find that error varies with season, with the worst predictions in winter, Fig.
5c. Examining the difference between models, we also find this seasonal
variation, meaning we need not know the error to deduce the regions of low
predictive power, Fig. 5d. On this data, the meta-model approach is almost as
good as the theoretical optimum strategy, Fig. 5d. Under this strategy, we see
that in the summer when the LF models agree, a LF model is used for
prediction, whereas in the winter the HF model is used, Fig. 5e-f.

Figure 5:  (a) Forecasting error distribution for HF and LF models. (b)
Seasonal prediction error, averaged across 5 years of testing data and over a
window of 5 days. (c) Seasonal prediction difference between models. Model
predictions were made, and the norm of their difference was calculated and
averaged as in (b). Difference between LSTM and linear models recovers
seasonal trend in error. (d) Total error of meta-modeling strategy against the
theoretical optimum and random strategy. (e-f) Seasonal distribution for where
the LF model is evaluated for $25\%$ and $50\%$ LF model evaluations.

##  5 Summary

In this work, we considered the forecasting problem for two real-world
geophysical data sets, the NOAA sea surface temperature and the NASA Daymet
land surface temperature data for North America. We confirmed that the
combination of dimensionality reduction through proper orthogonal
decomposition and recurrent neural network based predictions outperforms a
climatology baseline. While the neural network prediction outperforms simple
random forest and linear predictions, we were able to extract information from
these LF models nevertheless. Specifically, by using the two LF models, we
were able to devise a simple rule for deciding whether to use a HF or LF
model, which allowed a near optimal decision strategy. When used as a strategy
for complicated architectures, this could allow substantial computational
savings, especially when ensemble predictions are required. We were further
able to use the LF models to understand seasonal variation in predictive
power, finding that our models predictive power is highest in the summer for
land surface temperature, a valuable insight for climate modeling. In future,
it will be of interest to apply this to more advanced network architectures.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy
(DOE), Office of Science, Office of Advanced Scientific Computing Research,
under Contract DE-AC02-06CH11357. This research was funded in part and used
resources of the Argonne Leadership Computing Facility, which is a DOE Office
of Science User Facility supported under Contract DE-AC02-06CH11357. DS
acknowledges support from the NSF-MSGI fellowship. This paper describes
objective technical results and analysis. Any subjective views or opinions
that might be expressed in the paper do not necessarily represent the views of
the U.S. DOE or the United States Government.

##  Appendix A Proper Orthogonal Decomposition

Proper orthogonal decomposition (POD) provides a systematic method to project
dynamics of a high dimensional system onto a lower dimensional subspace. We
suppose that a single snapshot of the full system is a vector in
$\mathbb{R}^{N}$, where $N$ could be the number of grid points at which the
field is resolved. Observing the system across a number of time points gives
us the snapshots $\bm{\theta}_{1},\dots,\bm{\theta}_{T}$, with mean subtracted
by convention. The aim of POD, is to find a small set of orthonormal basis
vectors $\bm{v_{1}},\dots,\bm{v}_{M}$, with $M\ll N$, which approximates the
spatial snapshots,

$\bm{\theta}_{t}\approx\sum_{j=1}^{M}a_{j}(t)\bm{v}_{j},\quad t=1,\dots,T,$
(1)

and so allows us to approximate the evolution of the full $N$ dimensional
system, by considering only the evolution of the $M$ coefficients $a_{j}(t)$.
POD chooses the basis, $\bm{v}_{j}$, to minimize the residual with respect to
the $L_{2}$ norm,

$R=\sum_{t=1}^{T}||\bm{\theta}_{t}-\sum_{j=1}^{M}a_{j}(t)\bm{v}_{j}||^{2}.$
(2)

Defining the snapshot matrix,
$\bm{S}=[\bm{\theta}_{1}|\cdots|\bm{\theta}_{T}]$, the optimal basis is given
by the $M$ eigenvectors of $\bm{SS}^{T}$, with largest eigenvalues, after
which, the coefficients are found by orthogonal projection,
$a_{j}(t)=\langle\bm{\theta}_{t},\bm{v}_{j}\rangle$ [21].

For both of our data sets, we take only the training data snapshots, say
$\bm{D}_{1},\dots,\bm{D}_{T}$, from which we calculate the mean
$\bar{\bm{D}}=(1/T)\sum_{t}\bm{D}_{t}$, hence defining the mean subtracted
snapshots $\bm{\theta}_{t}=\bm{D}_{t}-\bar{\bm{D}}$. We then create the
snapshot matrix, $\bm{S}$, and find numerically the $M$ eigenvectors of
$\bm{S}\bm{S}^{T}$ with largest eigenvalues. From this, we train models,
$\mathcal{N}$, to forecast the coefficients

$\bm{a}(t+1)\approx\hat{\bm{a}}(t+1)=\mathcal{N}(\bm{a}(t),\bm{a}(t-1),\dots).$
(3)

making predictions of future coefficients given previous ones.

To test the predictions on unseen data, $\bm{E}_{1},\dots,\bm{E}_{k}$, we take
the mean $\bar{\bm{D}}$, and vectors $\bm{v}_{j}$ calculated from the training
data to get that

$a_{j}(t)=\langle\bm{E}_{t}-\bar{\bm{D}},\bm{v}_{j}\rangle,\quad j=1,\dots,M,$
(4)

which will be used by the model $\mathcal{N}$ to make a prediction for future
coefficients. The prediction for the coefficients $\hat{\bm{a}}$, can be
converted into predictions in the physical space by taking
$\bar{\bm{D}}+\sum_{j}\hat{a}_{j}\bm{v}_{j}$. This procedure only makes use of
testing data to pass into the model, not to train the model in any way.
Crucially, to make a forecast of $\bm{E}_{t+1}$, only previous measurements
$\bm{E}_{t},\bm{E}_{t-1},\dots$ are needed.

## References

  * [1] M. Auffhammer, P. Baylis, and C. H. Hausman.  Climate change is projected to have severe impacts on the frequency and intensity of peak electricity demand across the united states.  Proceedings of the National Academy of Sciences, 2017.
  * [2] P. Balaprakash, R. Egele, M. Salim, S. Wild, V. Vishwanath, F. Xia, T. Brettin, and R. Stevens.  Scalable reinforcement-learning-based neural architecture search for cancer deep learning research.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–33, 2019.
  * [3] N. D. Brenowitz and C. S. Bretherton.  Prognostic validation of a neural network unified physics parameterization.  Geophysical Research Letters, 45(12):6289–6298, 2018.
  * [4] A. Chattopadhyay, P. Hassanzadeh, and S. Pasha.  Predicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data.  Scientific Reports, 10(1):1317, Jan 2020.
  * [5] A. Chattopadhyay, E. Nabizadeh, and P. Hassanzadeh.  Analog forecasting of extreme-causing weather patterns using deep learning.  Journal of Advances in Modeling Earth Systems, 12(2):e2019MS001958, 2020.  e2019MS001958 10.1029/2019MS001958.
  * [6] T. W. Ford, P. A. Dirmeyer, and D. O. Benson.  Evaluation of heat wave forecasts seamlessly across subseasonal timescales.  npj Climate and Atmospheric Science, 1(1):20, Oct 2018.
  * [7] P. Gentine, M. Pritchard, S. Rasp, G. Reinaudi, and G. Yacalis.  Could machine learning break the convection parameterization deadlock?  Geophysical Research Letters, 45(11):5742–5751, 2018.
  * [8] T. Gneiting and A. E. Raftery.  Weather forecasting with ensemble methods.  Science, 310(5746):248–249, 2005.
  * [9] E. J. Kendon, N. M. Roberts, H. J. Fowler, M. J. Roberts, S. C. Chan, and C. A. Senior.  Heavier summer downpours with climate change revealed by weather forecast resolution model.  Nature Climate Change, 4(7):570–576, Jul 2014.
  * [10] P. Lynch.  The origins of computer weather prediction and climate modeling.  Journal of Computational Physics, 227(7):3431 – 3444, 2008.  Predicting weather, climate and extreme events.
  * [11] R. Maulik, R. Egele, B. Lusch, and P. Balaprakash.  Recurrent Neural Network Architecture Search for Geophysical Emulation.  arXiv e-prints, Apr. 2020.
  * [12] R. Maulik, A. Mohan, B. Lusch, S. Madireddy, P. Balaprakash, and D. Livescu.  Time-series learning of latent-space dynamics for reduced-order model closure.  Physica D: Nonlinear Phenomena, 405:132368, 2020.
  * [13] M. Milano and P. Koumoutsakos.  Neural network modeling for near wall turbulent flow.  Journal of Computational Physics, 182(1):1 – 26, 2002.
  * [14] R. Mizuta, H. Yoshimura, H. Murakami, M. Matsueda, H. Endo, T. Ose, K. Kamiguchi, M. Hosaka, M. Sugi, S. Yukimoto, et al.  Climate simulations using mri-agcm3. 2 with 20-km grid.  Journal of the Meteorological Society of Japan. Ser. II, 90:233–258, 2012.
  * [15] W. Ohfuchi, H. Sasaki, Y. Masumoto, and H. Nakamura.  Mesoscale resolving simulations of the global atmosphere and ocean on the earth simulator.  Eos, Transactions American Geophysical Union, 86(5):45–46, 2005\.
  * [16] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, and A. Ramadhan.  Universal differential equations for scientific machine learning.  arXiv preprint arXiv:2001.04385, 2020.
  * [17] M. Satoh, T. Matsuno, H. Tomita, H. Miura, T. Nasuno, and S. Iga.  Nonhydrostatic icosahedral atmospheric model (nicam) for global cloud resolving simulations.  Journal of Computational Physics, 227(7):3486 – 3514, 2008.  Predicting weather, climate and extreme events.
  * [18] S. Scher and G. Messori.  How global warming changes the difficulty of synoptic weather forecasting.  Geophysical Research Letters, 46(5):2931–2939, 2019.
  * [19] M. Schuster and K. K. Paliwal.  Bidirectional recurrent neural networks.  IEEE Transactions on Signal Processing, 45(11):2673–2681, 1997\.
  * [20] J. Sillmann, T. Thorarinsdottir, N. Keenlyside, N. Schaller, L. V. Alexander, G. Hegerl, S. I. Seneviratne, R. Vautard, X. Zhang, and F. W. Zwiers.  Understanding, modeling and predicting weather and climate extremes: Challenges and opportunities.  Weather and Climate Extremes, 18:65 – 74, 2017.
  * [21] K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. T. M. Dawson, and C.-A. Yeh.  Modal analysis of fluid flows: Applications and outlook.  AIAA Journal, 58(3):998–1022, 2020.
  * [22] Tin Kam Ho.  Random decision forests.  In Proceedings of 3rd International Conference on Document Analysis and Recognition, volume 1, pages 278–282 vol.1, 1995.