# Disorder-Induced Heating in Molecular Atmospheric Pressure Plasmas

J. LeVan  M. D. Acciarri  Department of Nuclear Engineering and Radiological
Sciences, University of Michigan, Ann Arbor, MI 48109, USA  S. D. Baalrud
Department of Nuclear Engineering and Radiological Sciences, University of
Michigan, Ann Arbor, MI 48109, USA

###### Abstract

Recent work has shown that ions are strongly coupled in atmospheric pressure
plasmas when the ionization fraction is sufficiently large, leading to a
temperature increase from disorder-induced heating that is not accounted for
in standard modelling techniques. Here, we extend this study to molecular
plasmas. A main finding is that the energy gained by ions in disorder-induced
heating gets spread over both translational and rotational degrees of freedom
on a nanosecond timescale, causing the final ion and neutral gas temperatures
to be lower in the molecular case than in the atomic case. A model is
developed for the equilibrium temperature that agrees well with molecular
dynamics simulations. The model and simulations are also applied to pressures
up to ten atmospheres. We conclude that disorder-induced heating is a
significant and predictable phenomena in molecular atmospheric pressure
plasmas.

CAPP, Atmospheric Pressure Plasmas, Strong Correlations, Strong Coupling,
Disorder Induced Heating, Molecular Dynamics, MD, LAMMPS

##  I Introduction

Cold atmospheric pressure plasmas (CAPPs) are continuously being tested and
developed for a wide range of applications, including wound healing and cancer
treatment [1], water purification [2], and food preservation [3], among
others[4, 5, 6]. This comes, in large part, from their unique capacity for
generating large quantities of reactive oxygen and nitrogen species (RONS).
Additionally, in comparison with other plasma technologies, CAPPs offer a
lower cost, complexity, and carbon footprint since they do not require the use
of a vacuum chamber [7]. This has driven increased interest in understanding
plasma dynamics, especially with concern to the transport of RONS [7]. A
significant amount of work has already been done to characterize and model
these plasmas[8, 9, 10, 11]. It is, however, also important to understand the
underlying physical mechanisms that drive CAPPs to be sure they are accounted
for in the models. In this paper, we examine and extend one such mechanism:
disorder-induced heating (DIH), which arises from strong ion-ion coupling.

Recent work showed that ions in atmospheric pressure plasmas are strongly
coupled at ionization fractions as low as $10^{-5}$, leading to rapid ion DIH
on a picosecond timescale following an instantaneous ionization[12]. This work
modelled the effect in monatomic argon plasma and found that, at high
ionization fractions ($\gtrsim 10^{-2}$), it can increase ion and neutral gas
temperatures by thousands of Kelvin. In this paper, we extend the previous
work to molecular plasmas and examine the effects of added rotational degrees
of freedom. Since many applications of CAPPs rely on N2 and O2 [13, 14, 15],
it is important to understand how disorder-induced heating impacts molecular
plasmas. We also extend the work to higher pressures (up to 10 atm) to account
for the regimes seen in plasma-assisted combustion [16].

Using molecular dynamics (MD) simulations of Nitrogen gas (N2), we find that
DIH causes the ion translational temperature to rapidly increase over
approximately one ion plasma period ($\approx$ 1 ps). This is followed by
translational-rotational and ion-neutral relaxation processes that occur in
roughly 1 ns. The latter process leads to fast neutral gas heating. A main
finding is that the translational-rotational energy exchange leads to a lower
total temperature resulting from DIH than in the atomic case. This is because
the kinetic energy gained by ions from DIH gets distributed amongst more
degrees of freedom in the molecular case. Considering N2, the kinetic energy
gain is eventually spread equally over 3 translational and 2 rotational
degrees of freedom, for a total of 5. This contrasts with the 3 degrees of
freedom in the atomic case. Furthermore, we confirm that this principle
extends to molecules with 3 rotational degrees of freedom by simulating an N4
molecular plasma.

Despite the reduction in temperature due to the presence of rotational degrees
of freedom, DIH can still cause a temperature increase on the order of
thousands of Kelvin. We derive a model for this change in temperature using
conservation of energy arguments and find that it matches well with our MD
results. This model is then generalized to account for higher pressures, where
DIH becomes more significant. Finally, we consider the common case of a
gradual ionization ramp over several nanoseconds. As in a previous study with
atomic gases, the temperature increase of ions due to DIH and the subsequent
ion-neutral energy relaxation are fast compared to the typical ionization
timescale. A consequence is that both the ion and neutral temperatures
increase at the rate of ionization. Importantly, it is also shown that
translational-rotational relaxation is fast compared to the ionization
timescale - thus, in the context of fast gas heating, one must consider the
effect of molecular rotation.

##  II Strong Coupling and Disorder-Induced Heating

An interaction between species $s$ and $s^{\prime}$ is strongly coupled if the
average potential energy is on the same order as, or greater than, the average
kinetic energy. This can be quantified by the coupling parameter exceeding
unity, $\Gamma_{ss^{\prime}}\geq 1$, where

$\Gamma_{ss^{\prime}}=\frac{\phi_{ss^{\prime}}(r=a_{ss^{\prime}})}{k_{B}T}.$
(1)

Here, $\phi_{ss^{\prime}}$ is the interaction potential between $s$ and
$s^{\prime}$, $a_{ss^{\prime}}$ is the average inter-particle spacing, and $T$
is the temperature of the two species.

Strongly coupled plasmas behave fundamentally differently than ideal (weakly
coupled) plasmas. Binary collisions characteristic of weak coupling are
replaced by many-body interactions that force the strongly coupled species
into an ordered spatial configuration. Consequently, when a gas is ionized to
the point of strong ion-ion coupling, the newly created ions move from a
highly disordered to a highly ordered state. This increased order results in a
lower potential energy which, in turn, is offset by an increase in kinetic
energy. This increase in kinetic energy is termed disorder-induced heating[17,
18]. The process is illustrated in figure 1.

Figure 1: Illustration of the disorder-induced heating process. Ionization
transforms neutral gas into a nonequilibrium plasma. Coulomb repulsion causes
the ions to separate on a picosecond timescale, causing significant
translational heating. Collisions then distribute this energy across neutral
and rotational degrees of freedom until they reach an equilibrium. The
collision relaxation timescale near atmospheric pressure is often
characteristic of nanoseconds.

Disorder-induced heating is not a new phenomenon. It has been studied in great
detail for decades in ultracold neutral plasmas[19, 18, 20]. However, only
recently has it been shown that DIH is also expected to occur in atmospheric
pressure plasmas[12]. We expect this to have significant ramifications in the
study of CAPPs[21, 22].

As an example, fast gas heating in nanosecond repetitively pulsed (NRP)
discharges has seen extensive study due to its applications in flow control,
combustion, and CO2 disassociation [23, 24, 25]. The heating is commonly
attributed in its entirety to the relaxation of electronically excited states
in atoms and molecules [26]. However, these discharges often see ionization
fractions well above $10^{-5}$ [27, 28, 29] so accordingly, DIH is expected to
occur. Since DIH occurs on a picosecond timescale, fast gas heating can no
longer be ascribed solely to the quenching of electronically excited states
and instead must also consider DIH[22].

##  III MD Simulations

Molecular dynamics simulations were run using the open-source software LAMMPS
[30]. Nitrogen (N2) molecules were modelled as rigid rotors, allowing for no
vibration, and all ions were modelled as N${}_{2}^{+}$, consisting of a
neutral atom rigidly bonded to an ion. A bond length of 1.09 Å[31] was held
constant using the “fix rigid” command in LAMMPS, which solves the forces and
torques on each rigid body by summing the forces and torques on its
constituent atoms. This approach ignores vibration, an assumption expected to
be valid because the characteristic vibrational temperature of N2 is $3374$ K
[32] and the temperature is far below this for the majority of the simulation.
Additionally, DIH and the subsequent ion-neutral energy relaxation occur at
the nanosecond timescale or less, while translational-vibrational energy
exchange in N2 occurs on a microsecond timescale [33]. A real ionization
process will certainly create excited vibrational states, but this does not
affect DIH because the ion and neutral atom temperatures are too low to couple
to the vibrational modes.

Electrons were left out of the simulation because their high temperature in
CAPPs make them weakly coupled. In strongly coupled plasmas, treating
electrons as a non-interacting background species yields accurate ion and
neutral dynamics, as commonly done in the one-component plasma (OCP) model
[34]. Chemical effects were also left out as they are beyond the scope of this
work. We aim not to simulate a CAPP as a whole, but instead examine a specific
heating mechanism. Since DIH is caused solely by a change in the spatial
orientation of ions, chemistry only affects the heating if it causes a
significant change in the number density or temperature of ion species in the
plasma on a relevant timescale. A model is provided (equation (7)) that can
account for those changes.

Interactions were modelled on a per-atom basis, with the ion-ion, ion-neutral,
and neutral-neutral interaction potentials given by the Coulomb, charge-
induced dipole, and Lennard-Jones potentials as follows

$\displaystyle\phi_{ii}(r)=\frac{q^{2}}{4\pi\epsilon_{0}}\frac{1}{r}$ (2)
$\displaystyle\phi_{in}(r)=\frac{q^{2}}{8\pi\epsilon_{0}}\frac{\alpha_{R}a_{0}^{3}}{r^{4}}\left(\frac{r_{\phi}^{8}}{r^{8}}-1\right)$
(3)
$\displaystyle\phi_{nn}(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]$
(4)

where $\epsilon=99.8k_{B}$, $\sigma=0.3667$ nm, and $\alpha_{R}=7.5$ for N2
and $a_{0}$ is the Bohr radius. The repulsive term in the charge-induced
dipole potential was added to avoid bound states, as done in previous
work[12]. Details on this addition can be found in the Appendix of Ref. 12.

Much of this paper will focus on data from the N2 simulations. However, for
comparison’s sake, simulations of N and N4 were also run. The atomic nitrogen
simulations were run in the exact same fashion as the previous work[12]. In
the N4 simulations, N4 was modelled using an open-chain geometry as described
by Glukhovtsev and Laiter [35]; a structure which, importantly, has three
rotational degrees of freedom. All N4 ions were modelled as N${}_{4}^{+}$,
with one of the two middle atoms ionized and a modified LJ coefficient
$\epsilon=29.0k_{B}$ was used in order to prevent significant neutral gas
correlations.

To start the simulations, 10,000 molecules of neutral gas was evolved under
the influence of a thermostat until equilibrium at room temperature was
reached. Then, a fraction of the molecules were instantly ionized and an NVE
(microcanonical) simulation was run with timestep $\Delta t=5\times
10^{-4}\omega_{pi}^{-1}$, a value found to be sufficiently small for
conserving energy. Molecule positions and velocities were output once every
plasma period and used to calculate translational and rotational temperatures.
The simulation domain was a 3-dimensional box with periodic boundary
conditions. Three simulations for each of N, N2, and N4 were run at every
tenth ionization fraction from $0.10$ to $0.70$. Data was averaged across the
three trials to minimize statistical deviation.

##  IV MD Results

The averaged time-evolution of the component temperatures from a particular
set of simulations, N2 at a 0.30 ionization fraction, can be seen in figure 2.
This plot shows results of the NVE stage simulation and $t=0$ refers to the
time immediately following ionization. A sharp spike can be seen in the
translational temperature of N${}_{2}^{+}$, reaching its peak after 2.22 ps
(1.52 $\omega_{pi}^{-1}$). This is disorder-induced heating. Molecules, upon
being ionized, structurally rearrange themselves until they reach their state
of minimal potential energy, increasing their translational kinetic energy in
the process. This heating is followed by collisional relaxation, during which
the N${}_{2}^{+}$ translational temperature equilibrates with the rest of the
degrees of freedom in the simulation over the course of one nanosecond,
ultimately reaching an equilibrated temperature of 732 K. For a clear physical
picture of this process, recall figure 1.

Simulations of each gas type, N, N2, and N4, followed a similar time-evolution
to figure 2. Ion translational temperatures always peaked after 1.5 plasma
periods, consistent with previous work [12], and equilibration always occurred
in roughly one nanosecond. However, they all produced differing equilibrium
temperatures. We found the equilibrium temperature scales with the ionization
fraction $x_{i}^{4/3}$, again consistent with the previous work[12]. This is
caused by an increased ion density and thus an increased chance of two ions
being generated in close enough proximity to gain an energy exceeding the
initial thermal energy as they repel.

It is also observed that the equilibrated temperature is inversely
proportional to the number of degrees of freedom present, implying
$T_{eq,N}>T_{eq,N_{2}}>T_{eq,N_{4}}$, as shown in figure 3. This comes as a
direct consequence of the equipartition theorem. Figure 4 shows that, for a
given ion density, N, N2, and N4 plasmas convert the same amount of potential
energy to kinetic energy. However, in N4, this energy must be spread across 6
degrees of freedom (3 translational and 3 rotational), while it is spread
across just 5 in N2 and 3 in N. It is clear, then, that DIH models must
consider the number of active degrees of freedom.

Figure 2: Simulated time-evolution of the temperature components of an N2
plasma at a 0.30 ionization fraction. Lines correspond to the ion
translational temperature (blue), ion rotational temperature (orange), neutral
molecule rotational temperature (green), and neutral molecule rotational
temperature (red). Figure 3: Simulated and predicted temperatures of N, N2,
and N4 plasma using the precise formulation from equation (8) with $g_{ii}(r)$
computed from MD, and from the approximate model of equation (10). Figure 4:
Change in total kinetic energy density (with respect to $t=0$) for N, N2, and
N4. This data comes from simulations at an ionization fraction of 0.30.

##  V Modelling DIH

###  V.1 Precise Formulation

Disorder-induced heating results from a change in the spatial ordering of
ions. The ion-ion radial distribution function, which outputs the ratio of the
local ion density around each ion to the bulk ion density[36]

$g_{ii}(r)=\frac{\langle\rho_{i}(r)\rangle}{\rho_{i}},$ (5)

can capture this change. Figure 5 shows an example of $g_{ii}(r)$ computed
from MD simulation data using equation (5) at various times throughout DIH.
The shifting of the curve to the right implies that ions are spreading apart
from each other and forming a highly correlated state with a reduced potential
energy.

The radial distribution function provides a convenient means to quantify the
potential energy. Integrating radially and multiplying by the ion density,
$\rho_{i}$, gives the density of ions at each distance $r+dr$ away from each
reference ion. Therefore, the potential energy of ions in terms of the ion-ion
distribution function is given by[37]

$U=\frac{N_{i}}{2}4\pi\rho_{i}\int\phi_{ii}(r)g_{ii}(r)r^{2}dr$ (6)

where $N_{i}$ is the total number of ions. For modelling the effect of
disorder-induced heating on the equilibrium temperature, we are only concerned
with how the potential energy of ions changes from immediately following
ionization to equilibrium. Additionally, to avoid simulation artifacts like
the number of ions, $N_{i}$, we are interested only in an energy density.
Therefore, we convert the previous equation to a change in energy density as
follows

$\Delta u_{i}=2\pi
x_{i}\rho_{i}\int\phi_{ii}(r)\left[g_{ii}(r,t_{eq})-g_{ii}(r,t_{0})\right]r^{2}dr$
(7)

where $x_{i}$ is the ionization fraction, $\phi_{ii}(r)$ is the Coulomb
potential, and $g_{ii}(r,t_{eq})$ and $g_{ii}(r,t_{0})$ are the ion-ion radial
distribution functions at equilibrium and immediately following ionization
respectively.

A change in the spatial orientation of ions also affects the ion-neutral
interaction, which can be quantified through the ion-neutral radial
distribution function ($g_{in}(r)$). However, we found the contribution to the
potential energy change from ion-neutral interactions to be much less than
from ion-ion interactions, and thus deemed the ion-neutral interaction to have
a negligible effect on the equilibrium temperature. It is worth noting that at
sufficiently high pressures, this assumption is unlikely to be valid because
the ion-neutral interaction becomes strongly coupled, becoming a source of
additional heating. However, at the pressures presented in this paper (up to
10 atm), we found this to be a valid approximation.

Hence, we use equation (7) to model the entirety of the change in potential
energy from DIH and convert this to a temperature change using traditional
thermodynamics arguments

$\Delta T=-\frac{2\Delta u_{i}}{fk_{B}},$ (8)

where $f$ is the number of degrees of freedom.

Figure 3 shows close agreement between the model and simulated temperatures.
Here, the radial distribution functions used in the model were computed from
the MD simulation data itself. We observe a slight discrepancy from the
simulated temperatures, which can be attributed to noise in the radial
distribution function. This noise can be seen in in figure 5, and comes as an
artifact of the finite number of ions and finite number of timesteps over
which we can compute this function.

Equation (8) is, by nature, a flexible formulation for the equilibrium
temperature. One can consider the effects of recombination by splitting
equation (7) into having an “initial” and “equilibrium” ionization fraction
and ion density. Additionally, it is additive, such that if multiple species
of ions exist in a plasma, one need only sum the contributions of each. To use
the model, there exist several methods for estimating radial distribution
functions including the Percus-Yevick[38] and hypernetted-chain[39]
approximations. Future work will use these methods to model DIH and be
implemented into a fluid model. This approach relies on numerically evaluating
a model for the radial distribution functions. Next, we consider a simpler
more approximate approach to estimating the potential energy change in DIH
directly.

Figure 5: Snapshots of the ion-ion radial distribution function over the
course of DIH from a simulation of N2 at a 0.30 ionization fraction. Over
time, the ion-ion radial distribution function shifts to the right as ions
spread apart from each other, reducing their potential energy in the process.

###  V.2 Approximate Model

A simpler estimate for DIH can be constructed by analyzing the change in
average ion-ion spacing. In the limiting case of transitioning from a
completely disordered to highly ordered state, we expect the average ion
spacing to change from approximately $a_{ii}/2$ to $a_{ii}$, leading to a
change in potential energy

$\Delta
u_{\max}=\frac{x_{i}e^{2}}{4\pi\epsilon_{0}}\left(\frac{1}{a_{ii}}-\frac{1}{a_{ii}/2}\right).$
(9)

We found that in the ion densities relevant to strongly coupled CAPPs, a
linear correction factor of $\kappa=2/3$ is required to fit our data, such
that $\Delta u=\kappa\Delta u_{\max}$. Converting to a change in temperature,
we find

$\Delta T=\frac{e^{2}}{2\pi\epsilon_{0}k_{B}}\frac{\kappa x_{i}}{fa_{ii}}$
(10)

Equation (10) is a slight modification of the model used in Ref. 12 to account
for the added degrees of freedom in molecular motion.

A comparison of the prediction of equation (10) and MD simulations is shown in
figure 3. The close agreement with the simple model is in part due to the use
of the $\kappa=2/3$ fit coefficient, but the model captures well the scaling
with ionization fraction. It also highlights the importance of the degrees of
freedom in distributing the energy gained by DIH. The change in temperature of
N ($f=3$), N2 ($f=5$), and N4 ($f=6$) differ only by a factor of $1/f$. While
it lacks the physical precision associated with the previous model, it can
still be useful as a fairly accurate means to estimate the extent to which DIH
will affect a given plasma.

##  VI Higher Pressures

Simulations were also run with initial gas pressures of 2.5, 5.0, 7.5, and 10
atm at constant ionization fractions to explore the pressure dependence of
DIH. This is particularly applicable to plasma-assisted combustion, where
initial gas pressures can range from 0.1-10 atm[16] and fast gas heating is
known to play a role in accelerating and intensifying the combustion
process[26].

A higher initial gas pressure leads to a larger ion density and thus, we can
expect more heating. MD results, shown in figure 6, confirm this.
Additionally, it is shown that the precise formulation from equation (8)
remains accurate, while the approximate model from equation (10) requires an
increased value for the correction factor $\kappa$ to agree with MD data.
Scaling $\kappa$ as $P^{1/20}$, where $P$ is the initial gas pressure, gives
good agreement.

Figure 6: Simulated and predicted temperatures of N2 as a function of the
initial gas pressure. The approximate model from equation (10) must be scaled
by $P^{1/20}$ to fit MD data.

##  VII Gradual Ionization

All simulations discussed thus far have involved an instantaneous ionization
process. However, in atmospheric pressure plasmas, ionization processes often
occur over the course of several nanoseconds[40, 27]. It has already been
shown that, in atomic plasma, disorder-induced heating is unaffected by
ionization dynamics[22]. Considering that translational-rotational relaxation
occurs on the same timescale as ion-neutral relaxation in these plasmas, we
expect this to be the case for molecular plasmas as well.

To confirm this, an MD simulation was run with a gradual ionization process.
The ion density was increased linearly over the course of 2.2 ns until an
ionization fraction of 0.30 was reached. The results, shown in figure 7,
confirm expectations. The simulation reaches an equilibrium temperature of 752
K, matching the simulations from an instantaneous ionization process to a 2.7%
error. It can also be observed that all degrees of freedom present in the
simulation heat at, roughly, the rate of ionization. This is because ion-
neutral and translational-rotational relaxation in these plasmas occur on a
faster timescale than the ionization process. Finally, it should be noted that
early on in the simulation, when the ionization fraction is low, there are not
enough ions to have a well-defined temperature. Consequently, the early
temperature data has a significant amount of noise.

The relaxation dynamics present in both the instantaneous and gradual
ionization process are complex. Strongly coupled ion-ion collisions play a
large role in the system’s translational-rotational energy exchange, so
attempts to model the relaxation with the Boltzmann equation[41] and Parker’s
model[42] (a model for RT relaxation in neutral gas) fail, as they do not
account for strong coupling. Thus, this work motivates the need for a theory
of translational-rotational relaxation in a strongly coupled plasma.

Figure 7: Simulated time-evolution of the temperature components of an N2
plasma subject to a gradual ionization process (2.2 ns) until a 0.30
ionization fraction is reached. Initial ion temperatures are made a lighter
color because early in the simulation, there are not enough ions to have a
well-defined temperature.

##  VIII Conclusions

This work examined how plasma and gas heating due to DIH changes with the
introduction of molecules and higher pressures. Using MD, it was shown that
energy from DIH is spread to rotational degrees of freedom on a nanosecond
timescale, leading to a lower equilibrium temperature than in the atomic case.
It was also shown that a higher initial gas pressure leads to a larger
temperature increase from DIH. Two models were derived that can predict this
behavior to a high degree of accuracy. Finally, by simulating a gradual
ionization over the course of a few nanoseconds, it was shown that the
temperature changed due to DIH in molecular plasmas is unaffected by
ionization dynamics.

This work shows that the effects of strong coupling remain significant in
molecular plasmas and establish that one must consider all active degrees of
freedom in order to accurately model the influence of DIH on the overall
heating process. It further motivates the need to develop new modelling
techniques that can account for strong coupling in CAPPs, as current models
based on the Boltzmann equation fail to do this. These models do not consider
DIH and as such, they risk vastly underestimating the temperature.

##  IX Acknowledgements

The authors thank Dr. Lucas Beving for suggesting this topic and Dr.
Christopher Moore for helpful conversations on this study. This research was
supported by the US Department of Energy under award no. DE-SC0022201. It was
also supported by computational resources and services provided by Advanced
Research Computing (ARC), a division of Information and Technology Services
(ITS) at the University of Michigan, Ann Arbor.

## References

  * [1] Thomas von Woedtke, Steffen Emmert, Hans-Robert Metelmann, Stefan Rupf, and Klaus-Dieter Weltmann.  Perspectives on cold atmospheric plasma (CAP) applications in medicine.  Physics of Plasmas, 27(7), 2020.  070601\.
  * [2] Azadeh Barjasteh, Zohreh Dehghani, Pradeep Lamichhane, Neha Kaushik, Eun Ha Choi, and Nagendra Kumar Kaushik.  Recent progress in applications of non-thermal plasma for water purification, bio-sterilization, and decontamination.  Applied Sciences, 11(8), 2021.
  * [3] Abirami R. Ganesan, Uma Tiwari, P. N. Ezhilarasi, and Gaurav Rajauria.  Application of cold plasma on food matrices: A review on current and future prospects.  Journal of Food Processing and Preservation, 45(1):e15070, 2021\.
  * [4] Noala Vicensoto Moreira Milhan, William Chiappim Junior, Aline da Graça Sampaio, Mariana Vegian, Rodrigo Pessoa, and Cristiane Koga Ito.  Applications of plasma-activated water in dentistry: A review.  International Journal of Molecular Sciences, 23:4131, 2022.
  * [5] S.P. Jeevan Kumar, Anjani Devi Chintagunta, Y. Mohan Reddy, Loïc Rajjou, Vijay Kumar Garlapati, Dinesh K. Agarwal, S. Rajendra Prasad, and Jesus Simal-Gandara.  Implications of reactive oxygen and nitrogen species in seed physiology for sustainable crop productivity under changing climate conditions.  Current Plant Biology, 26:100197, 2021.
  * [6] Som V. Thomas, Krista Dienger-Stambaugh, Michael Jordan, Yuxin Wang, Jason Hammonds, Paul Spearman, and Donglu Shi.  Inactivation of sars-cov-2 on surfaces by cold-plasma-generated reactive species.  Bioengineering, 10(3), 2023.
  * [7] Mária Domonkos, Petra Tichá, Jan Trejbal, and Pavel Demo.  Applications of cold atmospheric pressure plasma technology in medicine, agriculture and food industry.  Applied Sciences, 11(11), 2021.
  * [8] Erik C Neyts, Maksudbek Yusupov, Christof C Verlackt, and Annemie Bogaerts.  Computer simulations of plasma–biomolecule and plasma–tissue interactions for a better insight in plasma medicine.  Journal of Physics D: Applied Physics, 47(29):293001, 2014.
  * [9] Kexin Zhang, Camila A. Perussello, Vladimir Milosavljević, P. J. Cullen, Da-Wen Sun, and Brijesh K. Tiwari.  Diagnostics of plasma reactive species and induced chemistry of plasma treated foods.  Critical Reviews in Food Science and Nutrition, 59(5):812–825, 2019\.  PMID: 30676057.
  * [10] Yury Gorbanev, Angela Privat-Maldonado, and Annemie Bogaerts.  Analysis of short-lived reactive species in plasma-air-water systems: The dos and the do nots.  Analytical chemistry, 90(22):13151—13158, 2018.
  * [11] I Adamovich, S Agarwal, E Ahedo, L L Alves, S Baalrud, N Babaeva, A Bogaerts, A Bourdon, P J Bruggeman, C Canal, E H Choi, S Coulombe, Z Donkó, D B Graves, S Hamaguchi, D Hegemann, M Hori, H-H Kim, G M W Kroesen, M J Kushner, A Laricchiuta, X Li, T E Magin, S Mededovic Thagard, V Miller, A B Murphy, G S Oehrlein, N Puac, R M Sankaran, S Samukawa, M Shiratani, M Šimek, N Tarasenko, K Terashima, E Thomas Jr, J Trieschmann, S Tsikata, M M Turner, I J van der Walt, M C M van de Sanden, and T von Woedtke.  The 2022 plasma roadmap: low temperature plasma science and technology.  Journal of Physics D: Applied Physics, 55(37):373001, 2022.
  * [12] M D Acciarri, C Moore, and S D Baalrud.  Strong coulomb coupling influences ion and neutral temperatures in atmospheric pressure plasmas.  Plasma Sources Science and Technology, 31(12):125005, 2022.
  * [13] Khaled Lotfy, Nadi Awad Al-Harbi, and Hany Abd El-Raheem.  Cold atmospheric pressure nitrogen plasma jet for enhancement germination of wheat seeds.  Plasma Chemistry and Plasma Processing, 39(4):897–912, 2019.
  * [14] Katsuya Iuchi, Yukina Morisada, Yuri Yoshino, Takahiro Himuro, Yoji Saito, Tomoyuki Murakami, and Hisashi Hisatomi.  Cold atmospheric-pressure nitrogen plasma induces the production of reactive nitrogen species and cell death by increasing intracellular calcium in HEK293T cells.  Arch Biochem Biophys, 654:136–145, 2018.
  * [15] Aijun Yang, Dingxin Liu, Mingzhe Rong, Xiaohua Wang, and Michael G. Kong.  A dominant role of oxygen additive on cold atmospheric-pressure He + O2 plasmas.  Physics of Plasmas, 21(8), 2014.  083501\.
  * [16] Andrey Starikovskiy and Nickolay Aleksandrov.  Plasma-assisted ignition and combustion.  In Max Mulder, editor, Aeronautics and Astronautics, chapter 12. IntechOpen, Rijeka, 2011.
  * [17] M. S. Murillo.  Using fermi statistics to create strongly coupled ion plasmas in atom traps.  Phys. Rev. Lett., 87:115003, 2001.
  * [18] S. G. Kuzmin and T. M. O’Neil.  Numerical simulation of ultracold plasmas: How rapid intrinsic heating limits the development of correlation.  Phys. Rev. Lett., 88:065003, 2002.
  * [19] T. C. Killian, S. Kulin, S. D. Bergeson, L. A. Orozco, C. Orzel, and S. L. Rolston.  Creation of an ultracold neutral plasma.  Phys. Rev. Lett., 83:4776–4779, 1999.
  * [20] D. Murphy and B. M. Sparkes.  Disorder-induced heating of ultracold neutral plasmas created from atoms in partially filled optical lattices.  Phys. Rev. E, 94:021201, 2016.
  * [21] M D Acciarri, C Moore, and S D Baalrud.  Influence of strong coulomb coupling on diffusion in atmospheric pressure plasmas.  Plasma Sources Science and Technology, 32(11):115004, 2023.
  * [22] M D Acciarri, C Moore, and S D Baalrud.  Disorder-induced heating as a mechanism for fast neutral gas heating in atmospheric pressure plasmas.  Manuscript submitted for publication, 2023.
  * [23] N A Popov.  Fast gas heating in a nitrogen–oxygen discharge plasma: I. kinetic mechanism.  Journal of Physics D: Applied Physics, 44(28):285201, 2011.
  * [24] L. Cheng, N. Barleon, B. Cuenot, O. Vermorel, and A. Bourdon.  Plasma assisted combustion of methane-air mixtures: Validation and reduction.  Combustion and Flame, 240:111990, 2022.
  * [25] G V Pokrovskiy, N A Popov, and S M Starikovskaia.  Fast gas heating and kinetics of electronically excited states in a nanosecond capillary discharge in co2.  Plasma Sources Science and Technology, 31(3):035010, 2022.
  * [26] N.A. Popov and S.M. Starikovskaia.  Relaxation of electronic excitation in nitrogen/oxygen and fuel/air mixtures: fast gas heating in plasma-assisted ignition and flame stabilization.  Progress in Energy and Combustion Science, 91:100928, 2022.
  * [27] A Lo, A Cessou, C Lacour, B Lecordier, P Boubert, D A Xu, C O Laux, and P Vervisch.  Streamer-to-spark transition initiated by a nanosecond overvoltage pulsed discharge in air.  Plasma Sources Science and Technology, 26(4):045012, 2017.
  * [28] N Minesi, S Stepanyan, P Mariotto, G D Stancu, and C O Laux.  Fully ionized nanosecond discharges in air: the thermal spark.  Plasma Sources Science and Technology, 29(8):085003, 2020.
  * [29] R M van der Horst, T Verreycken, E M van Veldhuizen, and P J Bruggeman.  Time-resolved optical emission spectroscopy of nanosecond pulsed discharges in atmospheric-pressure n2 and n2/h2o mixtures.  Journal of Physics D: Applied Physics, 45(34):345201, 2012.
  * [30] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comp. Phys. Comm., 271:108171, 2022.
  * [31] David R Lide.  CRC handbook of chemistry and physics, volume 85.  CRC press, 2004.
  * [32] Donald Allan McQuarrie and John Douglas Simon.  Physical chemistry: a molecular approach, volume 1.  University science books Sausalito, CA, 1997.
  * [33] Wei Yang, Qianhong Zhou, Qiang Sun, Zhiwei Dong, and Eryan Yan.  Vibrational–translational relaxation in nitrogen discharge plasmas: Master equation modeling and Landau–Teller model revisited.  AIP Advances, 10(10):105311, 2020.
  * [34] Marc Baus and Jean-Pierre Hansen.  Statistical mechanics of simple coulomb systems.  Physics Reports, 59(1):1–94, 1980.
  * [35] Mikhail N. Glukhovtsev and Sergei Laiter.  Thermochemistry of tetrazete and tetraazatetrahedrane: A high-level computational study.  The Journal of Physical Chemistry, 100(5):1569–1577, 1996.
  * [36] J.P. Hansen and I.R. McDonald.  Theory of Simple Liquids.  Elsevier Science, 2006.
  * [37] Harvey Gould and Jan Tobochnik.  Statistical and thermal physics.  Princeton University Press, Princeton, NJ, 2010.
  * [38] George Stell.  The percus-yevick equation for the radial distribution function of a fluid.  Physica, 29(5):517–534, 1963.
  * [39] Tohru Morita.  Theory of Classical Fluids: Hyper-Netted Chain Approximation, I: Formulation for a One-Component System.  Progress of Theoretical Physics, 20(6):920–938, 1958.
  * [40] P Tardiveau, N Moreau, S Bentaleb, C Postel, and S Pasquiers.  Diffuse mode and diffuse-to-filamentary transition in a high pressure nanosecond scale corona discharge under high voltage.  Journal of Physics D: Applied Physics, 42(17):175202, 2009.
  * [41] J.H. Ferziger and H.G. Kaper.  Mathematical Theory of Transport Processes in Gases.  North-Holland Publishing Company, 1972.
  * [42] JG Parker.  Rotational and vibrational relaxation in diatomic gases.  The Physics of Fluids, 2(4):449–462, 1959.


Evolutionary Hirota Type (2+1)-Dimensional Bi-Hamiltonian Equations

Evolutionary Hirota Type (2+1)-Dimensional
Equations: Lax Pairs, Recursion Operators
and Bi-Hamiltonian Structures

Mikhail B. SHEFTEL † and Devrim YAZICI ‡

M.B. Sheftel and D. Yazıcı

† Department of Physics, Boğaziçi University, Bebek, 34342 Istanbul, Turkey
<EMAIL_ADDRESS>
‡ Department of Physics, Yıldız Technical University, Esenler, 34220 Istanbul,
Turkey<EMAIL_ADDRESS>
Received December 06, 2017, in final form March 02, 2018; Published online
March 07, 2018

We show that evolutionary Hirota type Euler–Lagrange equations in $(2+1)$
dimensions have a symplectic Monge–Ampère form. We consider integrable
equations of this type in the sense that they admit infinitely many
hydrodynamic reductions and determine Lax pairs for them. For two seven-
parameter families of integrable equations converted to two-component form we
have constructed Lagrangians, recursion operators and bi-Hamiltonian
representations. We have also presented a six-parameter family of tri-
Hamiltonian systems.

Lax pair; recursion operator; Hamiltonian operator; bi-Hamiltonian system

35Q75; 37K05; 37K10

##  1 Introduction

We study recursion operators, Lax pairs and bi-Hamiltonian representations for
$(2+1)$-dimensional equations of the evolutionary Hirota type

$\displaystyle u_{tt}=f(u_{t1},u_{t2},u_{11},u_{12},u_{22}).$ (1.1)

Here $u=u(t,z_{1},z_{2})$ and the subscripts denote partial derivatives of
$u$, namely, $u_{ij}=\partial^{2}u/\partial z_{i}\partial z_{j}$,
$u_{ti}=\partial^{2}u/\partial t\partial z_{i}$. Equations of this type arise
in a wide range of applications including non-linear physics, general
relativity, differential geometry and integrable systems. Some examples are
the Khokhlov–Zabolotskaya (dKP) equation in non-linear acoustics and the
theory of Einstein–Weyl structures and the Boyer–Finley equation in the theory
of self-dual gravity.

A lot of work has been done by E. Ferapontov et al. for studying integrability
of equation (1.1) which is understood as the existence of infinitely many
hydrodynamic reductions (see [2, 3, 4, 5, 6] and references therein).
Ferapontov et al. in [3] derived integrability condition which is equivalent
to the property of equation (1.1) to be linearizable by contact
transformations. This does not mean that there is a straightforward way to
obtain bi-Hamiltonian structures of the nonlinear equations (1.1) by a direct
transfer of such structures from the linear equation, because an arbitrary
transformation in the space of second partial derivatives of the unknown will
not preserve the bi-Hamiltonian structure of the equation.

Our goal here is to study bi-Hamiltonian structures of the integrable
equations of the form (1.1) together with Lax pairs and recursion operators.
We utilize the method which we used earlier [13, 17, 18] for constructing a
degenerate Lagrangian for two-component evolutionary form of the equation and
using Dirac’s theory of constraints [1] in order to obtain Hamiltonian form of
the system.

Our approach here starts with description of all equations (1.1) which have
the Euler–Lagrange form [14]. We do not consider here the equations of the
form (1.1) which become Lagrangian only after multiplication by an integrating
factor of variational calculus, postponing the appropriate generalization for
a further publication. We find that the Lagrangian evolutionary Hirota-type
equations have a symplectic Monge–Ampère form and we determine their
Lagrangians. Then we study recursion relations for symmetries and Lax pairs
for the Lagrangian equations. Here our starting point is to convert the
symmetry condition into a “skew-factorized” form. This approach extends the
method of A. Sergyeyev [16] for constructing recursion operators. According to
[16] one constructs the recursion operator from a certain Lax pair which, in
turn, is typically built from the original Lax pair for the equation under
study. On the other hand, below we construct such a special Lax pair using the
skew-factorized form of the linearized equation (symmetry condition) rather
than a previously known Lax pair, and then apply the construction from [16] to
obtain the recursion operator.

The next step is to transform Lagrangian equations converted into a two-
component form to a Hamiltonian system. Finally, composing a recursion
operator with a Hamiltonian operator we obtain the second Hamiltonian operator
and also find the corresponding Hamiltonian density. Thus, we end up with a
bi-Hamiltonian representation of an integrable equation (1.1) in a two-
component form. In this way, we obtain two seven-parameter families of bi-
Hamiltonian systems and a six-parameter family of tri-Hamiltonian systems.

The paper is organized as follows. In Section 2, we show that all equations
(1.1) of the Euler–Lagrange form have the symplectic Monge–Ampère form and we
derive a Lagrangian for such equations. In Section 3, we analyze the symmetry
condition for the Lagrangian equations. In Section 4, using an integrability
condition, we convert the symmetry condition into a “skew-factorized” form and
immediately extract Lax pair and recursion relations for symmetries from the
symmetry condition in this form. In Section 5, we convert our equation in a
two-component form and derive a degenerate Lagrangian for this system. In
Section 6, we transform the Lagrangian system into Hamiltonian system using
the Dirac’s theory of constraints. We obtain the Hamiltonian operator $J_{0}$
and corresponding Hamiltonian density $H_{1}$. In Section 7, we derive a
recursion operator $R$ in a $2\times 2$ matrix form using recursion relations
for the two-component form of the equation. In Section 8, by composing the
recursion operator with the Hamiltonian operator $J_{0}$ we obtain the second
Hamiltonian operator $J_{1}=RJ_{0}$ and the corresponding Hamiltonian density
$H_{0}$ with one additional “Hamiltonian” constraint for coefficients. Thus,
we obtain a seven-parameter family of bi-Hamiltonian systems because of the
two constraints on nine coefficients: integrability condition and Hamiltonian
condition.

We consider also an alternative skew-factorized representation of the symmetry
condition which implies different Lax pair, another recursion operator and a
different seven-parameter family of bi-Hamiltonian systems under a different
additional constraint on the coefficients. If, in addition, we require that
both additional constraints coincide and are compatible with the integrability
condition, we obtain a six-parameter family of tri-Hamiltonian systems.

In Sections 4, 7, and 8, we treat separately each of the two generic cases of
equation (2.2) when none of the coefficients $c_{1}$, $c_{2}$, $c_{3}$
vanishes, particular cases when either $c_{1}=0$ or $c_{2}=0$ which are
obtained as a specialization of one of the generic cases, and the special case
$c_{3}=0$ which cannot be obtained from the generic cases but should be
treated independently.

##  2 Lagrangian equations of evolutionary Hirota type

We start with equation (1.1) in the form

$\displaystyle F\equiv-u_{tt}+f(u_{t1},u_{t2},u_{11},u_{12},u_{22})=0.$ (2.1)

The Fréchet derivative operator (linearization) of equation (2.1) reads

$\displaystyle
D_{F}=-D_{t}^{2}+f_{u_{t1}}D_{t}D_{1}+f_{u_{t2}}D_{t}D_{2}+f_{u_{11}}D_{1}^{2}+f_{u_{12}}D_{1}D_{2}+f_{u_{22}}D_{2}^{2},$

where $D_{i}\equiv D_{z_{i}},D_{t}$ denote operators of total derivatives. The
adjoint Fréchet derivative operator has the form

$\displaystyle
D^{*}_{F}=-D_{t}^{2}+D_{t}D_{1}f_{u_{t1}}+D_{t}D_{2}f_{u_{t2}}+D_{1}^{2}f_{u_{11}}+D_{1}D_{2}f_{u_{12}}+D_{2}^{2}f_{u_{22}}.$

According to Helmholtz conditions [14], equation (2.1) is an Euler-Lagrange
equation for a variational problem iff its Fréchet derivative is self-adjoint,
$D^{*}_{F}=D_{F}$, or explicitly, equating to zero coefficients of $D_{t}$,
$D_{1}$, $D_{2}$ and the term without operators of total derivatives, we
obtain four equations on $f$

$\displaystyle D_{1}[f_{u_{t1}}]+D_{2}[f_{u_{t2}}]=0,$ $\displaystyle
D_{t}[f_{u_{t1}}]+2D_{1}[f_{u_{11}}]+D_{2}[f_{u_{12}}]=0,$ $\displaystyle
D_{t}[f_{u_{t2}}]+2D_{2}[f_{u_{22}}]+D_{1}[f_{u_{12}}]=0,$ $\displaystyle
D_{t}D_{1}[f_{u_{t1}}]+D_{t}D_{2}[f_{u_{t2}}]+D_{1}^{2}[f_{u_{11}}]+D_{1}D_{2}[f_{u_{12}}]+D_{2}^{2}[f_{u_{22}}]=0.$

The general solution of these equations for $f$ implies the Lagrangian
evolutionary Hirota equation (1.1) to have symplectic Monge–Ampère form

$\displaystyle
u_{tt}=c_{1}(u_{1t}u_{12}-u_{2t}u_{11})+c_{2}(u_{1t}u_{22}-u_{2t}u_{12})+c_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}+c_{4}u_{1t}+c_{5}u_{2t}$
$\displaystyle\hphantom{u_{tt}=}{}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22}+c_{9}.$
(2.2)

A Lagrangian for the equation (2.2) is readily obtained by applying the
homotopy formula [14] for $F=-u_{tt}+f$

$\displaystyle L[u]=\int_{0}^{1}u\cdot F[\lambda u]{\rm d}\lambda$

with the result

$\displaystyle
L=-\frac{1}{2}uu_{tt}+\frac{u}{3}\bigl{[}c_{1}(u_{1t}u_{12}-u_{2t}u_{11})+c_{2}(u_{1t}u_{22}-u_{2t}u_{12})+c_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}\bigr{]}$
$\displaystyle\hphantom{L=}{}+\frac{u}{2}(c_{4}u_{1t}+c_{5}u_{2t}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22})+c_{9}u.$
(2.3)

##  3 Symmetry condition

In the following it will be useful to introduce operator of derivative in the
direction of the vector $\vec{c}=(c_{1},c_{2})$:
$\nabla_{c}=\vec{c}\cdot\nabla=c_{1}D_{1}+c_{2}D_{2}$, so that equation (2.2)
may be written as

$\displaystyle
u_{tt}=u_{1t}\nabla_{c}(u_{2})-u_{2t}\nabla_{c}(u_{1})+c_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}$
$\displaystyle\hphantom{u_{tt}=}{}+c_{4}u_{1t}+c_{5}u_{2t}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22}+c_{9}.$
(3.1)

Symmetry condition is the differential compatibility condition of (2.2) and
the Lie equation $u_{\tau}=\varphi$, where $\varphi$ is the symmetry
characteristic and $\tau$ is the group parameter. It has the form of Fréchet
derivative (linearization) of equation (3.1)

$\displaystyle\varphi_{tt}=\nabla_{c}(u_{2})\varphi_{1t}+u_{1t}\nabla_{c}(\varphi_{2})-\nabla_{c}(u_{1})\varphi_{2t}-u_{2t}\nabla_{c}(\varphi_{1})$
$\displaystyle\hphantom{\varphi_{tt}=}{}+c_{3}(u_{22}\varphi_{11}+u_{11}\varphi_{22}-2u_{12}\varphi_{12})+c_{4}\varphi_{1t}+c_{5}\varphi_{2t}+c_{6}\varphi_{11}+c_{7}\varphi_{12}+c_{8}\varphi_{22}.$
(3.2)

It is convenient to introduce the following differential operators

$\displaystyle L_{12(1)}=u_{12}D_{1}-u_{11}D_{2},\qquad
L_{12(2)}=u_{22}D_{1}-u_{12}D_{2},$ $\displaystyle
L_{12(t)}=u_{2t}D_{1}-u_{1t}D_{2}=v_{2}D_{1}-v_{1}D_{2},\qquad v=u_{t},$

so that the symmetry condition (3.2) becomes

$\displaystyle\big{\\{}c_{1}(D_{t}L_{12(1)}-D_{1}L_{12(t)})+c_{2}(D_{t}L_{12(2)}-D_{2}L_{12(t)})+c_{3}(D_{1}L_{12(2)}-D_{2}L_{12(1)})-D_{t}^{2}$
$\displaystyle\qquad{}+c_{4}D_{1}D_{t}+c_{5}D_{2}D_{t}+c_{6}D_{1}^{2}+c_{7}D_{1}D_{2}+c_{8}D_{2}^{2}\big{\\}}\varphi=0.$
(3.3)

##  4 Recursion relations and Lax pairs

E. Ferapontov et al. in [3] have derived an integrability condition for the
symplectic Monge–Ampère equation of the following general form (equation (21)
in [3])

$\displaystyle\epsilon\det\left[\begin{matrix}u_{11}&u_{12}&u_{13}\\\
u_{12}&u_{22}&u_{23}\\\
u_{13}&u_{23}&u_{33}\end{matrix}\right]+h_{1}\big{(}u_{22}u_{33}-u_{23}^{2}\big{)}+h_{2}\big{(}u_{11}u_{33}-u_{13}^{2}\big{)}+h_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}$
$\displaystyle\qquad{}+g_{1}(u_{11}u_{23}-u_{12}u_{13})+g_{2}(u_{22}u_{13}-u_{12}u_{23})+g_{3}(u_{33}u_{12}-u_{13}u_{23})$
$\displaystyle\qquad{}+s_{1}u_{11}+s_{2}u_{22}+s_{3}u_{33}+\tau_{1}u_{23}+\tau_{2}u_{13}+\tau_{3}u_{12}+\nu=0$
(4.1)

with constant coefficients. Here integrability means that the equation (4.1)
admits infinitely many hydrodynamic reductions [5]. For our evolutionary
equation (2.2) the coefficients in (4.1) are

$\displaystyle h_{1}=h_{2}=h_{3}=g_{3}=0,\qquad g_{1}=-c_{1},\qquad
g_{2}=c_{2},\qquad s_{1}=c_{6},\qquad s_{2}=c_{8},$ $\displaystyle
s_{3}=-1,\qquad\tau_{1}=c_{5},\qquad\tau_{2}=c_{4},\qquad\tau_{3}=c_{7},\qquad\nu=c_{9}.$
(4.2)

The integrability condition given in [3, formula (22)] for the equation (4.1)
has the form

$\displaystyle
h_{1}^{2}s_{1}^{2}+h_{2}^{2}s_{2}^{2}+h_{3}^{2}s_{3}^{2}+g_{1}^{2}s_{2}s_{3}+g_{2}^{2}s_{1}s_{3}+g_{3}^{2}s_{1}s_{2}-2(h_{1}h_{2}s_{1}s_{2}+h_{1}h_{3}s_{1}s_{3}+h_{2}h_{3}s_{2}s_{3})$
$\displaystyle\qquad{}+4\epsilon s_{1}s_{2}s_{3}+4\nu
h_{1}h_{2}h_{3}+\epsilon\tau_{1}\tau_{2}\tau_{3}-\nu
g_{1}g_{2}g_{3}-\epsilon^{2}\nu^{2}-\nu\big{(}g_{1}^{2}h_{1}+g_{2}^{2}h_{2}+g_{3}^{2}h_{3}\big{)}$
$\displaystyle\qquad{}-(g_{1}\tau_{1}+g_{2}\tau_{2}+g_{3}\tau_{3}+2\epsilon\nu)(h_{1}s_{1}+h_{2}s_{2}+h_{3}s_{3}-\epsilon\nu)$
$\displaystyle\qquad{}+2(g_{1}h_{1}s_{1}\tau_{1}+g_{2}h_{2}s_{2}\tau_{2}+g_{3}h_{3}s_{3}\tau_{3})+\tau_{1}^{2}h_{2}h_{3}+\tau_{2}^{2}h_{1}h_{3}+\tau_{3}^{2}h_{1}h_{2}$
$\displaystyle\qquad{}-\epsilon(\tau_{1}^{2}s_{1}+\tau_{2}^{2}s_{2}+\tau_{3}^{2}s_{3})+s_{1}\tau_{1}g_{2}g_{3}+s_{2}\tau_{2}g_{1}g_{3}+s_{3}\tau_{3}g_{1}g_{2}$
$\displaystyle\qquad{}-(g_{1}h_{1}\tau_{2}\tau_{3}+g_{2}h_{2}\tau_{1}\tau_{3}+g_{3}h_{3}\tau_{1}\tau_{2})=0.$
(4.3)

For the equation (2.2) due to identifications (4.2), the integrability
condition (4.3) becomes

$\displaystyle
c_{2}(c_{1}c_{7}-c_{2}c_{6}+c_{3}c_{4})=c_{1}(c_{1}c_{8}+c_{3}c_{5})-c_{3}^{2}.$
(4.4)

We show explicitly that any equation of the form (2.2) satisfying (4.4) is
also integrable in the traditional sense by constructing Lax pair for such an
equation.

###  4.1 Generic case

In the integrability condition (4.4) we assume

$\displaystyle c_{1}\cdot c_{2}\cdot c_{3}\neq 0.$

The following procedure extends A. Sergyeyev’s method for constructing
recursion operators [16]. Namely, unlike [16], we start with the skew-
factorized form of the symmetry condition and extract from there a Lax pair
for symmetries instead of building it from a previously known Lax pair. After
that we construct a recursion operator from this newly found Lax pair using
Proposition 1 from [16].

The linear operator of the symmetry condition (3.3) for integrable equations
of the form (2.2) can be presented in the “skew-factorized” form

$\displaystyle(A_{1}B_{2}-A_{2}B_{1})\varphi=0.$ (4.5)

If we introduce two-dimensional vector operators $\vec{R}=(A_{1},A_{2})$ and
$\vec{S}=(B_{1},B_{2})$, then the skew-factorized form (4.5) becomes the cross
(vector) product $(\vec{R}\times\vec{S})\varphi=0$. Here differential
operators $A_{i}$ and $B_{i}$ are defined as

$\displaystyle A_{1}=c_{1}D_{t}-c_{3}D_{2},\qquad
A_{2}=-(c_{2}D_{t}+c_{3}D_{1}),$ $\displaystyle
B_{1}=c_{1}\big{(}c_{3}L_{12(2)}-c_{1}L_{12(t)}\big{)}+c_{1}c_{6}D_{1}+(c_{1}c_{7}-c_{2}c_{6}+c_{3}c_{4})D_{2},$
$\displaystyle
B_{2}=c_{1}\big{(}c_{3}L_{12(1)}+c_{2}L_{12(t)}\big{)}+(c_{3}c_{4}-c_{2}c_{6})D_{1}-c_{1}c_{8}D_{2}-c_{3}D_{t}.$
(4.6)

These operators satisfy the commutator relations

$\displaystyle[A_{1},A_{2}]=0,\qquad[A_{1},B_{2}]-[A_{2},B_{1}]=0,\qquad[B_{1},B_{2}]=0,$
(4.7)

where the last equation is satisfied on solutions of the equation (2.2).

It immediately follows that the following two operators also commute on
solutions

$\displaystyle X_{1}=\lambda A_{1}+B_{1},\qquad X_{2}=\lambda
A_{2}+B_{2},\qquad[X_{1},X_{2}]=0,$ (4.8)

and therefore constitute Lax representation for equation (2.2) with $\lambda$
being a spectral parameter.

Symmetry condition in the form (4.5) not only provides the Lax pair for
equation (2.2) but also leads directly to recursion relations for symmetries

$\displaystyle A_{1}\tilde{\varphi}=B_{1}\varphi,\qquad
A_{2}\tilde{\varphi}=B_{2}\varphi,$ (4.9)

where $\tilde{\varphi}$ is a potential for $\varphi$. This follows from a
special case of Proposition 1 of [16] which we recapitulate here together with
its proof for the readers’ convenience. Indeed, equations (4.9) together with
(4.7) imply $(A_{1}B_{2}-A_{2}B_{1})\varphi=[A_{1},A_{2}]\tilde{\varphi}=0$,
so $\varphi$ is a symmetry characteristic. Moreover, due to (4.9)

$\displaystyle(A_{1}B_{2}-A_{2}B_{1})\tilde{\varphi}=\bigl{(}[A_{1},B_{2}]-[A_{2},B_{1}]+B_{2}A_{1}-B_{1}A_{2}\bigr{)}\tilde{\varphi}=[B_{2},B_{1}]\varphi=0,$

which shows that $\tilde{\varphi}$ satisfies the symmetry condition (4.5) and
hence is also a symmetry. Thus, $\varphi$ is a symmetry whenever so is
$\tilde{\varphi}$ and vice versa. The equations (4.9) define an auto-Bäcklund
transformation between the symmetry conditions written for $\varphi$ and
$\tilde{\varphi}$. Hence, the auto-Bäcklund transformation for the symmetry
condition is nothing else than a recursion operator. Finally, we must remark
that this approach to recursion operators originated from the much older work
published in 90th [8, 12, 15].

The skew-factorized form of the symmetry condition is by no means unique. We
can derive another version by using the discrete symmetry transformation

$\displaystyle c_{1}\leftrightarrow-c_{2},\qquad c_{4}\leftrightarrow
c_{5},\qquad c_{6}\leftrightarrow c_{8},\qquad D_{1}\leftrightarrow-
D_{2},\qquad D_{t}\leftrightarrow-D_{t},\qquad v\leftrightarrow-v,$
$\displaystyle L_{12(1)}\leftrightarrow L_{12(2)},$ (4.10)

while the operator $L_{12(t)}$ is not changed. Applying (4.10) to operators
(4.6) we obtain a new set of operators

$\displaystyle A_{1}=c_{2}D_{t}+c_{3}D_{1},\qquad
A_{2}=c_{3}D_{2}-c_{1}D_{t},$ $\displaystyle
B_{1}=-c_{2}\big{(}c_{3}L_{12(1)}+c_{2}L_{12(t)}\big{)}+(c_{2}c_{7}-c_{1}c_{8}-c_{3}c_{5})D_{1}+c_{2}c_{8}D_{2},$
$\displaystyle
B_{2}=c_{2}\big{(}c_{1}L_{12(t)}-c_{3}L_{12(2)}\big{)}-c_{2}c_{6}D_{1}-(c_{1}c_{8}+c_{3}c_{5})D_{2}+c_{3}D_{t},$
(4.11)

which also satisfy the skew-factorized form (4.5) of the symmetry condition
(3.3) and the same commutator relations (4.7). Using these operators in (4.8)
and (4.9) we obtain the second Lax pair and another set of recursion relations
for symmetries, respectively.

It may be interesting to note algebraic relations between the two recursion
operators, namely, determined by the set (4.6), marked below with the
superscript (1), and by the set (4.11), marked with the superscript (2)

$\displaystyle A_{1}^{(2)}=-A_{2}^{(1)},\qquad A_{2}^{(2)}=-A_{1}^{(1)},$
$\displaystyle c_{2}B_{1}^{(1)}+c_{1}B_{2}^{(2)}=c_{3}A_{1}^{(1)},\qquad
c_{2}B_{2}^{(1)}+c_{1}B_{1}^{(2)}=c_{3}A_{2}^{(1)}.$

###  4.2 Particular cases

We now consider particular cases $c_{1}=0$, $c_{2}\neq 0$ and $c_{2}=0$,
$c_{1}\neq 0$.

In the first case integrability condition (4.4) reads

$\displaystyle c_{1}=0,\qquad c_{2}\neq 0\quad\Longrightarrow\quad
c_{2}(c_{2}c_{6}-c_{3}c_{4})=c_{3}^{2}.$ (4.12)

If we set $c_{1}=0$ in our first set of operators (4.6) from the generic case,
these operators become linearly dependent and the skew-factorized
representation (4.5) does not reproduce the symmetry condition (3.2).
Therefore, we have to put $c_{1}=0$ in the second set of operators (4.11) with
the result

$\displaystyle A_{1}=c_{2}D_{t}+c_{3}D_{1},\qquad A_{2}=c_{3}D_{2},$
$\displaystyle
B_{1}=-c_{2}\big{(}c_{3}L_{12(1)}+c_{2}L_{12(t)}\big{)}+(c_{2}c_{7}-c_{3}c_{5})D_{1}+c_{2}c_{8}D_{2},$
$\displaystyle
B_{2}=c_{3}D_{t}-c_{2}c_{3}L_{12(2)}-c_{2}c_{6}D_{1}-c_{3}c_{5}D_{2}.$ (4.13)

Operators (4.13) satisfy skew-factorized form (4.5) of the symmetry condition
(3.3) and the commutator relations (4.7). Therefore, as is shown above,
equations (4.9) yield the recursion relations for symmetries

$\displaystyle(c_{2}D_{t}+c_{3}D_{1})\tilde{\varphi}=\big{\\{}{-}c_{2}\big{(}c_{3}L_{12(1)}+c_{2}L_{12(t)}\big{)}+(c_{2}c_{7}-c_{3}c_{5})D_{1}+c_{2}c_{8}D_{2}\big{\\}}\varphi,$
$\displaystyle
c_{3}D_{2}\tilde{\varphi}=\big{(}c_{3}D_{t}-c_{2}c_{3}L_{12(2)}-c_{2}c_{6}D_{1}-c_{3}c_{5}D_{2}\big{)}\varphi,$

and the operators

$\displaystyle
X_{1}=\lambda(c_{2}D_{t}+c_{3}D_{1})-c_{2}(c_{3}L_{12(1)}+c_{2}L_{12(t)})+(c_{2}c_{7}-c_{3}c_{5})D_{1}+c_{2}c_{8}D_{2},$
$\displaystyle X_{2}=\lambda
c_{3}D_{2}+c_{3}D_{t}-c_{2}c_{3}L_{12(2)}+c_{2}c_{6}D_{1}+c_{3}c_{5}D_{2}$

commute on solutions and so constitute Lax pair for the equation (2.2) at
$c_{1}=0$.

In the second case $c_{2}=0$, integrability condition (4.4) has the form

$\displaystyle c_{2}=0,\qquad c_{1}\neq 0\quad\Longrightarrow\quad
c_{1}(c_{1}c_{8}+c_{3}c_{5})=c_{3}^{2}.$ (4.14)

We set $c_{2}=0$ in our first set of operators (4.6) in the generic case to
obtain

$\displaystyle A_{1}=c_{1}D_{t}-c_{3}D_{2},\qquad A_{2}=-c_{3}D_{1},$
$\displaystyle
B_{1}=c_{1}\big{(}c_{3}L_{12(2)}-c_{1}L_{12(t)}\big{)}+c_{1}c_{6}D_{1}+(c_{1}c_{7}+c_{3}c_{4})D_{2},$
$\displaystyle
B_{2}=c_{1}c_{3}L_{12(1)}+c_{3}c_{4}D_{1}-c_{1}c_{8}D_{2}-c_{3}D_{t}.$ (4.15)

Operators (4.15) also satisfy skew-factorized form (4.5) of the symmetry
condition (3.3) and the commutator relations (4.7). Therefore, as is shown
above, using these operators in (4.8) and (4.9) we obtain the Lax pair and
recursion relations for symmetries in this case.

We could not use the second family of operators (4.11) at $c_{2}=0$ because
skew-factorized form (4.5) would give identical zero instead of reproducing
the symmetry condition.

We note again that the skew-factorized form of the invariance condition (4.5)
is still not unique. One could obtain such a form at $c_{2}=0$ with a
different choice of the operators $A_{i}$ and $B_{i}$

$\displaystyle A_{1}=c_{1}D_{1},\qquad A_{2}=c_{3}D_{2}-c_{1}D_{t},$
$\displaystyle B_{1}=c_{1}^{2}L_{12(1)}+(c_{1}c_{5}-c_{3})D_{2}-c_{1}D_{t},$
$\displaystyle
B_{2}=c_{1}\big{\\{}c_{3}L_{12(2)}-c_{1}L_{12(t)}+c_{4}D_{t}+c_{6}D_{1}+c_{7}D_{2}\big{\\}}.$
(4.16)

Recursion relations and the Lax pair are still valid with the new definitions
(4.16). Since the latter choice leads to more complicated recursion operator
and second Hamiltonian operator, we stick to our previous definitions (4.15).

Finally, the case $c_{1}=0$ together with $c_{2}=0$ implies $c_{3}=0$ from the
integrability condition (4.4) and hence corresponds to linear equation (2.2).

###  4.3 Special case $\boldsymbol{c_{3}=0}$

In the integrability condition (4.4) we assume

$\displaystyle c_{1}\cdot c_{2}\neq 0,\qquad c_{3}=0\quad\Longrightarrow\quad
c_{1}c_{2}c_{7}=c_{1}^{2}c_{8}+c_{2}^{2}c_{6}.$

Then the linear operator of the symmetry condition (3.3) can be presented in
the skew-factorized form (4.5) with the following operators $A_{i}$ and
$B_{i}$

$\displaystyle A_{1}=c_{1}D_{t},\qquad
A_{2}=c_{1}D_{1}+c_{2}D_{2}=\nabla_{c},$ $\displaystyle
B_{1}=c_{1}c_{2}L_{12(t)}-c_{2}c_{6}D_{1}-c_{1}c_{8}D_{2},$ $\displaystyle
B_{2}=c_{2}\big{\\{}c_{1}L_{12(1)}+c_{2}L_{12(2)}+c_{4}D_{1}+c_{5}D_{2}-D_{t}\big{\\}}.$
(4.17)

Operators (4.17) satisfy the commutator relations (4.7) and hence, as is shown
above, the equations (4.9) produce the recursion relations for symmetries

$\displaystyle
c_{1}D_{t}\tilde{\varphi}=\\{c_{1}c_{2}L_{12(t)}-c_{2}c_{6}D_{1}-c_{1}c_{8}D_{2}\\}\varphi,$
$\displaystyle\nabla_{c}(\tilde{\varphi})=c_{2}\\{c_{1}L_{12(1)}+c_{2}L_{12(2)})+c_{4}D_{1}+c_{5}D_{2}-D_{t}\\}\varphi,$
(4.18)

and the operators

$\displaystyle X_{1}=\lambda
c_{1}D_{t}+c_{1}c_{2}L_{12(t)}-c_{2}c_{6}D_{1}-c_{1}c_{8}D_{2},$
$\displaystyle
X_{2}=\lambda\nabla_{c}+c_{2}\big{\\{}c_{1}L_{12(1)}+c_{2}L_{12(2)}+c_{4}D_{1}+c_{5}D_{2}-D_{t}\big{\\}}.$

commute on solutions and so constitute Lax representation for the equation
(2.2) at $c_{3}=0$.

We note that if we considered this case as a particular case at $c_{3}=0$ of
the operators (4.6) or (4.11) of the generic case, then operators $A_{i}$,
$B_{i}$ would be linearly dependent and the skew-factorized form (4.5) would
not yield the symmetry condition (3.3).

##  5 Two-component form

Introducing the additional dependent variable $v=u_{t}$, we convert equation
(3.1) into the evolutionary system

$\displaystyle u_{t}=v,$ $\displaystyle
v_{t}=v_{1}\nabla_{c}(u_{2})-v_{2}\nabla_{c}(u_{1})+c_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}$
$\displaystyle\hphantom{v_{t}=}{}+c_{4}v_{1}+c_{5}v_{2}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22}+c_{9}.$
(5.1)

Lie equations become $u_{\tau}=\varphi$, $v_{\tau}=\psi$, so that $u_{t}=v$
implies the first invariance condition $\varphi_{t}=\psi$.

The second invariance condition is obtained by differentiating the second
equation in (5.1) with respect to the group parameter $\tau$

$\displaystyle\psi_{t}=\nabla_{c}(u_{2})\psi_{1}+v_{1}\nabla_{c}(\varphi_{2})-\nabla_{c}(u_{1})\psi_{2}-v_{2}\nabla_{c}(\varphi_{1})$
$\displaystyle\hphantom{\psi_{t}=}{}+c_{3}(u_{22}\varphi_{11}+u_{11}\varphi_{22}-2u_{12}\varphi_{12})+c_{4}\psi_{1}+c_{5}\psi_{2}+c_{6}\varphi_{11}+c_{7}\varphi_{12}+c_{8}\varphi_{22}.$

The Lagrangian for system (5.1) is obtained by a suitable modification of the
Lagrangian (2.3) of the one-component equation (2.2), skipping some total
derivative terms

$\displaystyle
L=u_{t}v-\frac{v^{2}}{2}+\frac{u_{t}}{3}\\{c_{1}(u_{2}u_{11}-u_{1}u_{12})+c_{2}(u_{2}u_{12}-u_{1}u_{22})\\}$
$\displaystyle\hphantom{L=}{}-\frac{u_{t}}{2}(c_{4}u_{1}+c_{5}u_{2})+\frac{c_{3}}{3}u\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}+\frac{u}{2}(c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22})+c_{9}u.$
(5.2)

##  6 Hamiltonian representation

To transform from Lagrangian to Hamiltonian description, we define canonical
momenta

$\displaystyle\pi_{u}=\frac{\partial L}{\partial
u_{t}}=v+\frac{1}{3}\bigl{\\{}c_{1}(u_{2}u_{11}-u_{1}u_{12})+c_{2}(u_{2}u_{12}-u_{1}u_{22})\bigr{\\}}-\frac{1}{2}(c_{4}u_{1}+c_{5}u_{2}),$
$\displaystyle\pi_{v}=\frac{\partial L}{\partial v_{t}}=0,$ (6.1)

which satisfy canonical Poisson brackets

$\displaystyle\big{[}\pi_{i}(z),u^{k}(z^{\prime})\big{]}=\delta_{i}^{k}\delta(z-z^{\prime}),$

where $u^{1}=u$, $u^{2}=v$, $z=(z_{1},z_{2})$, the only nonzero Poisson
bracket being
$[\pi_{u},u]=\delta(z_{1}-z_{1}^{\prime})\delta(z_{2}-z_{2}^{\prime})$. The
Lagrangian (5.2) is degenerate because the momenta cannot be inverted for the
velocities. Therefore, following the theory of Dirac’s constraints [1], we
impose (6.1) as constraints

$\displaystyle\Phi_{u}=\pi_{u}-v-\frac{1}{3}\bigl{[}c_{1}(u_{2}u_{11}-u_{1}u_{12})+c_{2}(u_{2}u_{12}-u_{1}u_{22})\bigr{]}+\frac{1}{2}(c_{4}u_{1}+c_{5}u_{2}),$
$\displaystyle\Phi_{v}=\pi_{v},$

and calculate Poisson bracket for the constraints

$\displaystyle
K_{11}=[\Phi_{u}(z_{1},z_{2}),\Phi_{u^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})],\qquad
K_{12}=[\Phi_{u}(z_{1},z_{2}),\Phi_{v^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})],$
$\displaystyle
K_{21}=[\Phi_{v}(z_{1},z_{2}),\Phi_{u^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})],\qquad
K_{22}=[\Phi_{v}(z_{1},z_{2}),\Phi_{v^{\prime}}(z^{\prime}_{1},z^{\prime}_{2})].$

We obtain the following matrix of Poisson brackets, which for convenience we
multiply by the overall factor $(-1)$

$\displaystyle K=\left(\begin{matrix}-K_{11}&-1\\\ 1&0\end{matrix}\right),$
(6.2)

where

$\displaystyle
K_{11}=c_{1}(u_{11}D_{2}-u_{12}D_{1})+c_{2}(u_{12}D_{2}-u_{22}D_{1})-c_{4}D_{1}-c_{5}D_{2}.$

The Hamiltonian operator is the inverse to the symplectic operator
$J_{0}=K^{-1}$

$\displaystyle J_{0}=\left(\begin{matrix}0&1\\\
-1&-K_{11}\end{matrix}\right)=\left(\begin{matrix}0&1\\\
-1&c_{1}L_{12(1)}+c_{2}L_{12(2)}+c_{4}D_{1}+c_{5}D_{2}\end{matrix}\right).$
(6.3)

Operator $J_{0}$ is Hamiltonian if and only if its inverse $K$ is symplectic
[7], which means that the volume integral $\Omega=\iiint_{V}\omega{\rm d}V$ of
$\omega=(1/2){\rm d}u^{i}\wedge K_{ij}{\rm d}u^{j}$ should be a symplectic
form, i.e., at appropriate boundary conditions ${\rm d}\omega=0$ modulo total
divergence. Another way of formulation is to say that the vertical
differential of $\omega$ should vanish [9]. In $\omega$ summations over $i$,
$j$ run from 1 to 2 and $u^{1}=u$, $u^{2}=v$. Using (6.2), we obtain

$\displaystyle\omega=\frac{1}{2}\bigl{[}c_{1}(u_{12}{\rm d}u\wedge{\rm
d}u_{1}-u_{11}{\rm d}u\wedge{\rm d}u_{2})+c_{2}(u_{22}{\rm d}u\wedge{\rm
d}u_{1}-u_{12}{\rm d}u\wedge{\rm d}u_{2})$
$\displaystyle\hphantom{\omega=}{}-c_{4}{\rm d}u\wedge{\rm d}u_{1}-c_{5}{\rm
d}u\wedge{\rm d}u_{2}-2{\rm d}u\wedge{\rm d}v\bigr{]}.$ (6.4)

Taking exterior derivative of (6.4) and skipping total divergence terms, we
have checked that ${\rm d}\omega=0$ which proves that operator $K$ is
symplectic and hence $J_{0}$ defined in (6.3) is indeed a Hamiltonian
operator.

The first Hamiltonian form of this system is

$\displaystyle\left(\begin{matrix}u_{t}\\\
v_{t}\end{matrix}\right)=J_{0}\left(\begin{matrix}\delta_{u}H_{1}\\\
\delta_{v}H_{1}\end{matrix}\right),$

where we still need to determine the corresponding Hamiltonian density
$H_{1}$. We convert $L$ from (5.2) to the form

$\displaystyle
L=u_{t}\pi_{u}-\frac{v^{2}}{2}+\frac{c_{3}}{3}u\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}+\frac{u}{2}(c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22})+c_{9}u,$

and apply the formula $H_{1}=\pi_{u}u_{t}+\pi_{v}v_{t}-L$, where $\pi_{v}=0$,
with the final result

$\displaystyle
H_{1}=\frac{v^{2}}{2}-\frac{c_{3}}{3}u\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}-\frac{1}{2}u(c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22})-c_{9}u.$

##  7 Recursion operators in $\boldsymbol{2\times 2}$ matrix form

###  7.1 Generic case: $\boldsymbol{c_{1}\cdot c_{2}\cdot c_{3}\neq 0}$

We define two-component symmetry characteristic $(\varphi,\psi)^{\rm T}$
(where T means transposed matrix) with $\psi=\varphi_{t}$ and
$(\tilde{\varphi},\tilde{\psi})^{\rm T}$ with
$\tilde{\psi}=\tilde{\varphi}_{t}$ for the original and transformed
symmetries, respectively. The recursion relations (4.9) with the use of (4.6)
for the operators $A_{i}$, $B_{i}$ take the form

$\displaystyle
c_{1}\tilde{\psi}-c_{3}D_{2}\tilde{\varphi}=B_{1}\varphi\equiv\big{\\{}c_{1}(c_{3}L_{12(2)}-c_{1}L_{12(t)})+c_{1}c_{6}D_{1}+(c_{3}c_{4}+c_{1}c_{7}-c_{2}c_{6})D_{2}\big{\\}}\varphi,\\!\\!\\!$
(7.1) $\displaystyle-
c_{2}\tilde{\psi}-c_{3}D_{1}\tilde{\varphi}=B_{2}\varphi\equiv\big{\\{}c_{1}(c_{3}L_{12(1)}+c_{2}L_{12(t)})+(c_{3}c_{4}-c_{2}c_{6})D_{1}-c_{1}c_{8}D_{2}\big{\\}}\varphi-
c_{3}\psi.$

Combining these two equations to eliminate first $\tilde{\psi}$ and then
$\tilde{\varphi}$, we obtain

$\displaystyle\nabla_{c}(\tilde{\varphi})=c_{1}\big{\\{}\nabla_{c}(u_{1})\varphi_{2}-\nabla_{c}(u_{2})\varphi_{1}\big{\\}}-c_{1}c_{4}\varphi_{1}-(c_{1}c_{5}-c_{3})\varphi_{2}+c_{1}\psi,$
$\displaystyle\nabla_{c}(\tilde{\psi})=c_{1}\big{\\{}\nabla_{c}(v_{1}\varphi_{2}-v_{2}\varphi_{1})+c_{3}(u_{22}\varphi_{11}+u_{11}\varphi_{22}-2u_{12}\varphi_{12})$
$\displaystyle\hphantom{\nabla_{c}(\tilde{\psi})=}{}+c_{6}\varphi_{11}+c_{7}\varphi_{12}+c_{8}\varphi_{22}\big{\\}}+c_{3}\psi_{2},$

where the subscripts of $\varphi$ and $\psi$ denote partial derivatives.
Applying the inverse operator $\nabla_{c}^{-1}$, which is defined to satisfy
the relations $\nabla_{c}^{-1}\nabla_{c}=1$, we obtain the explicit form of
recursion relations

$\displaystyle\tilde{\varphi}=\nabla_{c}^{-1}\big{\langle}c_{1}\big{\\{}\nabla_{c}(u_{1})\varphi_{2}-\nabla_{c}(u_{2})\varphi_{1}-c_{4}\varphi_{1}-c_{5}\varphi_{2}\big{\\}}+c_{3}\varphi_{2}+c_{1}\psi\big{\rangle},$
$\displaystyle\tilde{\psi}=c_{1}(v_{1}\varphi_{2}-v_{2}\varphi_{1})+\nabla_{c}^{-1}\big{\langle}c_{1}\big{\\{}c_{3}(u_{22}\varphi_{11}+u_{11}\varphi_{22}-2u_{12}\varphi_{12})$
$\displaystyle\hphantom{\tilde{\psi}=}{}+c_{6}\varphi_{11}+c_{7}\varphi_{12}+c_{8}\varphi_{22}\big{\\}}+c_{3}\psi_{2}\big{\rangle}.$
(7.2)

Here an important remark is due. The operator $\nabla_{c}^{-1}$ can make sense
merely as a formal inverse of $\nabla_{c}$. Thus, the relations (7.2) are
formal as well. The proper interpretation of the quantities like
$\nabla_{c}^{-1}$ and of (7.2) requires the language of differential
coverings, see the original papers [8, 12] and the recent survey [9].

In a two component form, the recursion relations (7.1) read

$\displaystyle\left(\begin{matrix}\tilde{\varphi}\\\
\tilde{\psi}\end{matrix}\right)=R\left(\begin{matrix}\varphi\\\
\psi\end{matrix}\right)$ (7.3)

with the recursion operator $R$ in the $2\times 2$ matrix form

$\displaystyle R=\left(\begin{matrix}R_{11}&c_{1}\nabla_{c}^{-1}\\\
R_{21}&c_{3}\nabla_{c}^{-1}D_{2}\end{matrix}\right)$ (7.4)

with the matrix elements

$\displaystyle
R_{11}=\nabla_{c}^{-1}\big{\langle}c_{1}\big{\\{}\nabla_{c}(u_{1})D_{2}-\nabla_{c}(u_{2})D_{1}-c_{4}D_{1}-c_{5}D_{2}\big{\\}}+c_{3}D_{2}\big{\rangle},$
$\displaystyle
R_{21}=c_{1}\big{\langle}v_{1}D_{2}-v_{2}D_{1}+\nabla_{c}^{-1}\big{\\{}c_{3}\big{(}u_{22}D_{1}^{2}+u_{11}D_{2}^{2}-2u_{12}D_{1}D_{2}\big{)}$
$\displaystyle\hphantom{R_{21}=}{}+c_{6}D_{1}^{2}+c_{7}D_{1}D_{2}+c_{8}D_{2}^{2}\big{\\}}\big{\rangle}.$
(7.5)

Next, we use the alternative set (4.11) of operators $A_{i}$, $B_{i}$ in the
recursion relations (4.9) presented in a two-component form

$\displaystyle
c_{2}\tilde{\psi}+c_{3}D_{1}\tilde{\varphi}=B_{1}\varphi\equiv\big{\\{}{-}c_{2}(c_{3}L_{12(1)}+c_{2}L_{12(t)})+(c_{2}c_{7}-c_{1}c_{8}-c_{3}c_{5})D_{1}+c_{2}c_{8}D_{2}\big{\\}}\varphi,$
$\displaystyle
c_{3}D_{2}\tilde{\varphi}-c_{1}\tilde{\psi}=B_{2}\varphi\equiv\big{\\{}c_{2}(c_{1}L_{12(t)}-c_{3}L_{12(2)})-c_{2}c_{6}D_{1}-(c_{1}c_{8}+c_{3}c_{5})D_{2}\big{\\}}\varphi+c_{3}\psi.$

Combining these two equations to eliminate first $\tilde{\psi}$ and then
$\tilde{\varphi}$, we obtain the explicit form of recursion relations

$\displaystyle\tilde{\varphi}=\nabla_{c}^{-1}\big{\\{}{-}c_{2}(c_{1}L_{12(1)}+c_{2}L_{12(2)})-(c_{2}c_{4}+c_{3})D_{1}-c_{2}c_{5}D_{2}\big{\\}}\varphi+c_{2}\nabla_{c}^{-1}\psi,$
$\displaystyle\tilde{\psi}=\left\langle
c_{3}\nabla_{c}^{-1}D_{1}\left\\{c_{1}L_{12(1)}+c_{2}L_{12(2)}+\left(c_{4}+\frac{c_{3}}{c_{2}}\right)D_{1}+c_{5}D_{2}\right\\}\right.$
$\displaystyle\left.\hphantom{\tilde{\psi}=}{}+\left\\{-(c_{3}L_{12(1)}+c_{2}L_{12(t)})+\frac{1}{c_{2}}(c_{2}c_{7}-c_{1}c_{8}-c_{3}c_{5})D_{1}+c_{8}D_{2}\right\\}\right\rangle\varphi-
c_{3}\nabla_{c}^{-1}D_{1}\psi,$

and we immediately extract the second recursion operator $R^{\prime}$ in the
$2\times 2$ matrix form

$\displaystyle
R^{\prime}=\left(\begin{matrix}R^{\prime}_{11}&c_{2}\nabla_{c}^{-1}\\\
R^{\prime}_{21}&-c_{3}\nabla_{c}^{-1}D_{1}\end{matrix}\right),$ (7.6)

where

$\displaystyle
R^{\prime}_{11}=-c_{2}\nabla_{c}^{-1}\left\\{c_{1}L_{12(1)}+c_{2}L_{12(2)}+\left(c_{4}+\frac{c_{3}}{c_{2}}\right)D_{1}+c_{5}D_{2}\right\\},$
$\displaystyle
R^{\prime}_{21}=c_{3}\nabla_{c}^{-1}D_{1}\left\\{c_{1}L_{12(1)}+c_{2}L_{12(2)}+\left(c_{4}+\frac{c_{3}}{c_{2}}\right)D_{1}+c_{5}D_{2}\right\\}$
$\displaystyle\hphantom{R^{\prime}_{21}=}{}-(c_{3}L_{12(1)}+c_{2}L_{12(t)})+\frac{1}{c_{2}}(c_{2}c_{7}-c_{1}c_{8}-c_{3}c_{5})D_{1}+c_{8}D_{2}.$
(7.7)

###  7.2 Particular cases

We consider again particular cases $c_{1}=0$, $c_{2}\neq 0$ and $c_{2}=0$,
$c_{1}\neq 0$.

As we know from Section 4, the first case, $c_{1}=0$, $c_{2}\neq 0$ should be
considered as a particular case of the second recursion operator $R^{\prime}$
from (7.6) and (7.7) with the result

$\displaystyle R^{\prime}=\left(\begin{matrix}\displaystyle-
c_{2}D_{2}^{-1}\\!\left(L_{12(2)}+\frac{c_{6}}{c_{3}}D_{1}\right)-c_{5},&D_{2}^{-1}\vspace{1mm}\\\
R^{\prime}_{21},&-\displaystyle\frac{c_{3}}{c_{2}}D_{2}^{-1}D_{1}\end{matrix}\right),$

where

$\displaystyle
R^{\prime}_{21}=c_{3}D_{2}^{-1}D_{1}\left(L_{12(2)}+\frac{c_{6}}{c_{3}}D_{1}\right)-\big{(}c_{3}L_{12(1)}+c_{2}L_{12(t)}\big{)}+c_{7}D_{1}+c_{8}D_{2}$

and the integrability condition (4.12) has been used.

The second case, $c_{2}=0$, $c_{1}\neq 0$ should be considered as a particular
case of the first recursion operator $R$ from (7.4) and (7.5) with the result

$\displaystyle R=\left(\begin{matrix}R_{11}&D_{1}^{-1}\vspace{1mm}\\\
R_{21}&\displaystyle\frac{c_{3}}{c_{1}}D_{1}^{-1}D_{2}\end{matrix}\right)$

with the matrix elements

$\displaystyle
R_{11}=-D_{1}^{-1}\big{\\{}c_{1}L_{12(1)}+(c_{1}c_{5}-c_{3})D_{2}\big{\\}}-c_{1}c_{4},$
$\displaystyle
R_{21}=-c_{1}L_{12(t)}+c_{6}D_{1}+c_{7}D_{2}+D_{1}^{-1}\big{\\{}c_{3}\big{(}u_{22}D_{1}^{2}+u_{11}D_{2}^{2}-2u_{12}D_{1}D_{2}\big{)}+c_{8}D_{2}^{2}\big{\\}}.$

###  7.3 Special case: $\boldsymbol{c_{1}\cdot c_{2}\neq 0}$,
$\boldsymbol{c_{3}=0}$

Recursion relations (4.18) in a two-component form become

$\displaystyle
c_{1}\tilde{\psi}=c_{1}c_{2}(v_{2}\varphi_{1}-v_{1}\varphi_{2})-c_{2}c_{6}\varphi_{1}-c_{1}c_{8}\varphi_{2},$
$\displaystyle\nabla_{c}(\tilde{\varphi})=c_{2}\big{\\{}\nabla_{c}(u_{2})\varphi_{1}-\nabla_{c}(u_{1})\varphi_{2}+c_{4}\varphi_{1}+c_{5}\varphi_{2}-\psi\big{\\}}.$
(7.8)

The explicit two-component form of the recursion relations (7.8) is

$\displaystyle\tilde{\varphi}=\nabla_{c}^{-1}c_{2}\big{\\{}\nabla_{c}(u_{2})\varphi_{1}-\nabla_{c}(u_{1})\varphi_{2}+c_{4}\varphi_{1}+c_{5}\varphi_{2}-\psi\big{\\}},$
$\displaystyle\tilde{\psi}=c_{2}(v_{2}\varphi_{1}-v_{1}\varphi_{2})-\frac{c_{2}c_{6}}{c_{1}}\varphi_{1}-c_{8}\varphi_{2}.$
(7.9)

In the matrix form (7.3), the recursion operator arising from (7.9) reads

$\displaystyle
R=\left(\begin{matrix}c_{2}\nabla_{c}^{-1}\\{\nabla_{c}(u_{2})D_{1}-\nabla_{c}(u_{1})D_{2}+c_{4}D_{1}+c_{5}D_{2}\\},&-c_{2}\nabla_{c}^{-1}\vspace{1mm}\\\
c_{2}(v_{2}D_{1}-v_{1}D_{2})-\displaystyle\frac{c_{2}c_{6}}{c_{1}}D_{1}-c_{8}D_{2},&0\end{matrix}\right).$
(7.10)

##  8 Bi-Hamiltonian systems

###  8.1 Generic case

####  8.1.1 First family of bi-Hamiltonian systems

The second Hamiltonian operator $J_{1}$ is obtained by composing the recursion
operator (7.4) with the first Hamiltonian operator $J_{1}=RJ_{0}$

$\displaystyle\left(\begin{matrix}J_{1}^{11}&J_{1}^{12}\\\
J_{1}^{21}&J_{1}^{22}\end{matrix}\right)=\left(\begin{matrix}R_{11}&c_{1}\nabla_{c}^{-1}\\\
R_{21}&c_{3}\nabla_{c}^{-1}D_{2}\end{matrix}\right)\left(\begin{matrix}0&1\\\
-1&\nabla_{c}(u_{2})D_{1}-\nabla_{c}(u_{1})D_{2}+c_{4}D_{1}+c_{5}D_{2}\end{matrix}\right),$

where we have used an alternative equivalent expression for the matrix element
$J_{0}^{22}$, with the final result

$\displaystyle
J_{1}=\left(\begin{matrix}-c_{1}\nabla_{c}^{-1}&c_{3}\nabla_{c}^{-1}D_{2}\\\
-c_{3}\nabla_{c}^{-1}D_{2}&J_{1}^{22}\end{matrix}\right),$ (8.1)

where

$\displaystyle
J_{1}^{22}=c_{3}L_{12(2)}-c_{1}L_{12(t)}+\nabla_{c}^{-1}\big{\\{}c_{1}\big{(}c_{6}D_{1}^{2}+c_{7}D_{1}D_{2}+c_{8}D_{2}^{2}\big{)}+c_{3}D_{2}(c_{4}D_{1}+c_{5}D_{2})\big{\\}}.$

Here operator $J_{1}$ is manifestly skew symmetric. A check of the Jacobi
identities and compatibility of the two Hamiltonian structures $J_{0}$ and
$J_{1}$ is straightforward but too lengthy to be presented here. The method of
the functional multi-vectors for checking the Jacobi identity and the
compatibility of the Hamiltonian operators is developed by P. Olver in [14,
Chapter 7] and has been applied recently for checking bi-Hamiltonian structure
of the general heavenly equation [18] and the first heavenly equation of
Plebański [17] under the well-founded conjecture that this method is
applicable for nonlocal Hamiltonian operators as well.

The next problem is to derive the Hamiltonian density $H_{0}$ corresponding to
the second Hamiltonian operator $J_{1}$ such that implies the bi-Hamiltonian
representation of the system (5.1)

$\displaystyle\left(\begin{matrix}u_{t}\\\
v_{t}\end{matrix}\right)=J_{0}\left(\begin{matrix}\delta_{u}H_{1}\\\
\delta_{v}H_{1}\end{matrix}\right)=J_{1}\left(\begin{matrix}\delta_{u}H_{0}\\\
\delta_{v}H_{0}\end{matrix}\right)=\left(\begin{matrix}v\\\
v_{t}\end{matrix}\right),$ (8.2)

where $v_{t}$ should be replaced by the right-hand side of the second equation
in (5.1). Then we could conclude that our system (5.1) is also integrable in
the sense of Magri [10, 11].

######  Proposition 8.1.

Bi-Hamiltonian representation (8.2) of the system (5.1) is valid under the
constraint

$\displaystyle
c_{9}=\frac{1}{c_{1}^{2}}\big{\\{}c_{6}(c_{3}-c_{1}c_{5})+c_{4}(c_{1}c_{7}-c_{2}c_{6}+c_{3}c_{4})\big{\\}}$
(8.3)

with the following Hamiltonian density

$\displaystyle
H_{0}=v\left\\{\frac{1}{c_{1}}\nabla_{c}(u)+\frac{c_{4}}{c_{1}}z_{2}+\frac{(c_{3}-c_{1}c_{5})}{c_{1}^{2}}z_{1}+s_{0}\right\\}-\frac{c_{3}}{2c_{1}^{2}}u_{2}\nabla_{c}(u)+\frac{c_{3}c_{4}}{c_{1}^{2}}u.$
(8.4)

###### Proof.

We will need the following simple relation between $J_{1}^{22}$ and the
operator $B_{1}$ from (4.6), which is involved in the recursion relations
(4.9) and Lax pair (4.8)

$\displaystyle
J_{1}^{22}=\frac{1}{c_{1}}\big{(}B_{1}+\nabla_{c}^{-1}c_{3}^{2}D_{2}^{2}\big{)}.$
(8.5)

Acting by the first row of $J_{1}$ on the column of variational derivatives of
$H_{0}$ in (8.2) and applying $\nabla_{c}$ we obtain

$\displaystyle\nabla_{c}\big{(}J_{1}^{11}\delta_{u}H_{0}+J_{1}^{12}\delta_{v}H_{0}\big{)}\equiv
c_{3}D_{2}\delta_{v}H_{0}-c_{1}\delta_{u}H_{0}=\nabla_{c}(v)=c_{1}v_{1}+c_{2}v_{2}.$
(8.6)

On account of the relation (8.5), the second row of the last equation in (8.2)
reads

$\displaystyle
J_{1}^{21}\delta_{u}H_{0}+J_{1}^{22}\delta_{v}H_{0}\equiv\frac{1}{c_{1}}B_{1}\delta_{v}H_{0}+\frac{c_{3}}{c_{1}}D_{2}\nabla_{c}^{-1}(c_{3}D_{2}\delta_{v}H_{0}-c_{1}\delta_{u}H_{0})=v_{t},$

which with the use of (8.6) becomes

$\displaystyle\frac{1}{c_{1}}B_{1}(\delta_{v}H_{0})+\frac{c_{3}}{c_{1}}v_{2}=v_{t}\quad\iff\quad
B_{1}(\delta_{v}H_{0})=c_{1}v_{t}-c_{3}v_{2}.$ (8.7)

We assume a linear dependence of $H_{0}$ on $v$

$\displaystyle
H_{0}=b[u]v+c[u]\quad\Longrightarrow\quad\delta_{v}H_{0}=\frac{\partial
H_{0}}{\partial v}=b[u],$ (8.8)

where $b$ and $c$ depend only on $u$ and its derivatives.

We note that adding to $H_{0}$ the term $av^{2}$ with constant $a$ does not
contribute to the Hamiltonian flow

$\displaystyle J_{1}\left(\begin{matrix}\delta_{u}H_{0}\\\
\delta_{v}H_{0}\end{matrix}\right)$

and hence this is unnecessary.

Plugging (8.8) into equation (8.7) and using the definition of $B_{1}$ from
(4.6), we obtain

$\displaystyle
c_{1}\big{\\{}c_{3}(u_{22}D_{1}[b]-u_{12}D_{2}[b])+c_{1}(v_{1}D_{2}[b]-v_{2}D_{1}[b])\big{\\}}+c_{1}c_{6}D_{1}[b]$
$\displaystyle\qquad{}+(c_{1}c_{7}-c_{2}c_{6}+c_{3}c_{4})D_{2}[b]=c_{1}\big{\\{}v_{1}\nabla_{c}(u_{2})-v_{2}\nabla_{c}(u_{1})+c_{3}\big{(}u_{11}u_{22}-u_{12}^{2}\big{)}$
$\displaystyle\qquad{}+c_{4}v_{1}+c_{5}v_{2}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22}+c_{9}\big{\\}}-c_{3}v_{2}.$
(8.9)

Splitting equation (8.9) with respect to $v_{1}$ and $v_{2}$ and collecting
separately terms with $v_{1}$ and $v_{2}$ implies the following two equations

$\displaystyle
D_{1}[b]=\frac{1}{c_{1}}\left\\{\nabla_{c}(u_{1})+\frac{(c_{3}-c_{1}c_{5})}{c_{1}}\right\\},\qquad
D_{2}[b]=\frac{1}{c_{1}}(\nabla_{c}(u_{2})+c_{4}).$ (8.10)

Integrating these two equations we obtain

$\displaystyle
b=\frac{1}{c_{1}}\nabla_{c}(u)+\frac{c_{4}}{c_{1}}z_{2}+\frac{(c_{3}-c_{1}c_{5})}{c_{1}^{2}}z_{1}+s_{0},$
(8.11)

where $s_{0}$ is a constant of integration. Plugging (8.10) into the remaining
terms in (8.9), we note that the terms $c_{1}c_{3}(u_{11}u_{22}-u_{12}^{2})$
and all terms proportional to $u_{11}$ and $u_{12}$ are canceled identically
while the terms proportional to $u_{22}$ cancel due to integrability condition
(4.4). The remaining constant terms in (8.9) imply the relation (8.3) which is
an additional (“Hamiltonian”) condition for integrable system (5.1) to have
bi-Hamiltonian form.

Using the result (8.11) for $b$ in our ansatz (8.8) for $H_{0}$ yields

$\displaystyle
H_{0}=v\left\\{\frac{1}{c_{1}}\nabla_{c}(u)+\frac{c_{4}}{c_{1}}z_{2}+\frac{(c_{3}-c_{1}c_{5})}{c_{1}^{2}}z_{1}+s_{0}\right\\}+c[u].$
(8.12)

From (8.6) we have

$\displaystyle\delta_{u}H_{0}\equiv\delta_{u}(b[u]v)+\delta_{u}(c[u])=\frac{1}{c_{1}}\big{\\{}c_{3}D_{2}[b]-\nabla_{c}(v)\big{\\}},$
(8.13)

where (8.11) implies that

$\displaystyle\delta_{u}(bv)=-\frac{1}{c_{1}}\nabla_{c}(v).$

After cancelations, equation (8.13) yields

$\displaystyle\delta_{u}(c[u])=\frac{c_{3}}{c_{1}}D_{2}[b]\equiv\frac{c_{3}}{c_{1}^{2}}\big{(}\nabla_{c}(u_{2})+c_{4}\big{)},$
(8.14)

where $D_{2}[b]$ from (8.10) is used. This result obviously suggests $c[u]$ to
be of the form

$\displaystyle c[u]=\alpha u_{1}u_{2}+\beta u_{2}^{2}+\gamma(u)$ (8.15)

with constant $\alpha$ and $\beta$, so that $\delta_{u}(c)=-D_{1}(\alpha
u_{2})-D_{2}(\alpha u_{1}+2\beta u_{2})+\gamma^{\prime}(u)$. Plugging the
latter expression in the equation (8.14) we find the coefficients in (8.15) to
be

$\displaystyle\alpha=-\frac{c_{3}}{2c_{1}},\qquad\beta=-\frac{c_{2}c_{3}}{2c_{1}^{2}},\qquad\gamma(u)=\frac{c_{3}c_{4}}{c_{1}^{2}}u+\gamma_{0}$

with the constant of integration $\gamma_{0}$, so that $c[u]$ in (8.15)
becomes

$\displaystyle
c[u]=-\frac{c_{3}}{2c_{1}^{2}}u_{2}\nabla_{c}(u)+\frac{c_{3}c_{4}}{c_{1}^{2}}u+\gamma_{0}.$

We have now completely determined the Hamiltonian density $H_{0}$ in (8.12) to
be (8.4), where we have skipped the additive constant $\gamma_{0}$. ∎

With the Hamiltonian density $H_{0}$ from (8.4), corresponding to the second
Hamiltonian operator $J_{1}$, the system (5.1) admits bi-Hamiltonian
representation (8.2), provided that the integrability condition (4.4) and bi-
Hamiltonian constraint (8.3) are satisfied in the generic case. Thus, we
obtain the first seven-parameter family of bi-Hamiltonian systems.

####  8.1.2 Second family of bi-Hamiltonian systems

Composing an alternative recursion operator $R^{\prime}$, determined by (7.6)
and (7.7), with the first Hamiltonian operator $J_{0}$ we obtain an
alternative second Hamiltonian operator $J^{\prime}_{1}=R^{\prime}J_{0}$ with
the resulting expression

$\displaystyle
J_{1}^{\prime}=\left(\begin{matrix}-c_{2}\nabla_{c}^{-1}&-c_{3}\nabla_{c}^{-1}D_{1}\vspace{1mm}\\\
c_{3}\nabla_{c}^{-1}D_{1}&\displaystyle\frac{1}{c_{2}}\left(c_{3}^{2}\nabla_{c}^{-1}D_{1}^{2}+B_{1})\right)\end{matrix}\right),$
(8.16)

where $B_{1}$ is defined in (4.11). Next task is to determine the Hamiltonian
density $H_{0}^{\prime}$ corresponding to $J_{1}^{\prime}$ in the bi-
Hamiltonian representation

$\displaystyle J_{0}\left(\begin{matrix}\delta_{u}H_{1}\\\
\delta_{v}H_{1}\end{matrix}\right)=J^{\prime}_{1}\left(\begin{matrix}\delta_{u}H^{\prime}_{0}\\\
\delta_{v}H^{\prime}_{0}\end{matrix}\right)=\left(\begin{matrix}v\\\
v_{t}\end{matrix}\right).$ (8.17)

A check of the Jacobi identities for $J^{\prime}_{1}$ and compatibility of the
two Hamiltonian structures $J_{0}$ and $J^{\prime}_{1}$ is straightforward and
too lengthy to be presented here.

######  Proposition 8.2.

Bi-Hamiltonian representation (8.17) of the system (5.1) is valid under the
constraint

$\displaystyle
c_{9}=\frac{1}{c_{2}^{2}}\big{\\{}c_{8}(c_{2}c_{4}+c_{3})+c_{5}(c_{1}c_{8}-c_{2}c_{7}+c_{3}c_{5})\big{\\}}$
(8.18)

with the following Hamiltonian density

$\displaystyle
H^{\prime}_{0}=\frac{v}{c_{2}}\left\\{\nabla_{c}(u)-c_{5}z_{1}+\left(c_{4}+\frac{c_{3}}{c_{2}}\right)z_{2}\right\\}+s_{0}v+\frac{c_{3}}{2c_{2}^{2}}\big{(}u_{1}\nabla_{c}(u)+2c_{5}u\big{)}.$
(8.19)

###### Proof.

Acting by the first row of $J^{\prime}_{1}$ on the column of variational
derivatives of $H^{\prime}_{0}$ in (8.17) and applying $\nabla_{c}$ we obtain

$\displaystyle\delta_{u}H^{\prime}_{0}=-\frac{1}{c_{2}}\big{(}c_{3}D_{1}\delta_{v}H^{\prime}_{0}+\nabla_{c}v\big{)}.$
(8.20)

The second row of equation (8.17) becomes

$\displaystyle
c_{3}\nabla_{c}^{-1}D_{1}\delta_{u}H^{\prime}_{0}+\frac{1}{c_{2}}\big{(}c_{3}^{2}\nabla_{c}^{-1}D_{1}^{2}+B_{1}\big{)}\delta_{v}H^{\prime}_{0}=v_{t},$

which with the use of (8.20) becomes

$\displaystyle B_{1}(\delta_{v}H^{\prime}_{0})=c_{2}v_{t}+c_{3}v_{1}.$ (8.21)

We again assume $H^{\prime}_{0}$ to be linear in $v$

$\displaystyle
H^{\prime}_{0}=b[u]v+c[u]\quad\Longrightarrow\quad\delta_{v}H^{\prime}_{0}=\frac{\partial
H^{\prime}_{0}}{\partial v}=b[u],$ (8.22)

where $b$ and $c$ depend only on $u$ and its derivatives. Plugging (8.22) into
equation (8.21), using the expression for $B_{1}$ from (4.11) and $v_{t}$ from
our equation (5.1), we obtain the equation linear in $v_{1}$ and $v_{2}$,
similar to (8.9). Equating the coefficients of $v_{1}$ and $v_{2}$ in this
linear equation on both sides, we obtain

$\displaystyle D_{1}[b]=\frac{1}{c_{2}}\\{\nabla_{c}(u_{1})-c_{5}\\},\qquad
D_{2}[b]=\frac{1}{c_{2}}\left\\{\nabla_{c}(u_{2})+c_{4}+\frac{c_{3}}{c_{2}}\right\\}.$
(8.23)

Integrating these equations we obtain

$\displaystyle
b=\frac{1}{c_{2}}\left\\{\nabla_{c}(u)-c_{5}z_{1}+\left(c_{4}+\frac{c_{3}}{c_{2}}\right)z_{2}\right\\}+s_{0}.$
(8.24)

Plugging the expressions (8.23) in remaining terms of the equation linear in
$v_{1}$ and $v_{2}$, we find that all variable terms cancel and we end up with
the relation (8.18) among constant terms. This is again an additional
condition for integrable system (5.1) to have bi-Hamiltonian form which is in
general different from (8.3) obtained for our first choice of the operators
$A_{i}$, $B_{i}$. Thus, we obtain two different families of seven-parameter
bi-Hamiltonian systems.

Plugging the expression (8.24) for $b$ into the formula (8.22) for
$H^{\prime}_{0}$, calculating $\delta_{u}H^{\prime}_{0}$ and using equation
(8.20), we obtain

$\displaystyle\delta_{u}c[u]=-\frac{c_{3}}{c_{2}^{2}}(c_{1}u_{11}+c_{2}u_{12}-c_{5}),$

which implies

$\displaystyle
c[u]=\frac{c_{3}}{2c_{2}^{2}}(u_{1}\nabla_{c}(u)+2c_{5}u)+\gamma_{0}.$ (8.25)

Using our results (8.24) for $b$ and (8.25) for $c[u]$ in the formula (8.22),
we obtain the Hamiltonian density $H^{\prime}_{0}$ to be (8.19) for the
alternative bi-Hamiltonian system (8.17) in the generic case. Here again we
have skipped the additive constant $\gamma_{0}$. ∎

####  8.1.3 A family of tri-Hamiltonian systems

If we impose an additional condition on the coefficients such that both
constraints (8.3) and (8.18) coincide, we eliminate $c_{9}$ and consider the
resulting condition together with the integrability condition (4.4) and the
constraint (8.3). With these three conditions, we obtain a six-parameter
family of tri-Hamiltonian systems

$\displaystyle\left(\begin{matrix}u_{t}\\\
v_{t}\end{matrix}\right)=J_{0}\left(\begin{matrix}\delta_{u}H_{1}\\\
\delta_{v}H_{1}\end{matrix}\right)=J_{1}\left(\begin{matrix}\delta_{u}H_{0}\\\
\delta_{v}H_{0}\end{matrix}\right)=J^{\prime}_{1}\left(\begin{matrix}\delta_{u}H^{\prime}_{0}\\\
\delta_{v}H^{\prime}_{0}\end{matrix}\right)$

with compatibility of the three Hamiltonian operators $J_{0}$, $J_{1}$,
$J^{\prime}_{1}$ been checked. An explicit solution of the two constraints
together with the integrability condition (4.4) for the tri-Hamiltonian system
has the form

$\displaystyle
c_{4}=\frac{1}{c_{2}c_{3}}\big{(}c_{2}^{2}c_{6}-c_{1}^{2}c_{8}-c_{1}c_{3}c_{5}\big{)},\qquad
c_{7}=\frac{1}{c_{1}c_{2}}\big{(}2c_{1}^{2}c_{8}+2c_{1}c_{3}c_{5}-c_{3}^{2}\big{)},$
$\displaystyle
c_{9}=\frac{(c_{3}c_{5}+c_{1}c_{8})}{c_{1}c_{2}^{2}c_{3}}\big{\\{}c_{3}(c_{3}-c_{1}c_{5})-c_{1}^{2}c_{8}\big{\\}}+\frac{c_{6}c_{8}}{c_{3}}$

with the six free parameters $c_{1}$, $c_{2}$, $c_{3}$, $c_{5}$, $c_{6}$,
$c_{8}$.

###  8.2 Particular cases

The first case $c_{1}=0$, $c_{2}\neq 0$ we should treat as a particular case
of the second bi-Hamiltonian family. Therefore, we put $c_{1}=0$ in the second
Hamiltonian operator $J^{\prime}_{1}$ from (8.16)

$\displaystyle
J^{\prime}_{1}=\left(\begin{matrix}-D_{2}^{-1}&-\displaystyle\frac{c_{3}}{c_{2}}D_{2}^{-1}D_{1}\vspace{1mm}\\\
\displaystyle\frac{c_{3}}{c_{2}}D_{2}^{-1}D_{1}&\displaystyle\frac{c_{3}^{2}}{c_{2}^{2}}D_{2}^{-1}D_{1}^{2}+\frac{1}{c_{2}}B_{1}\end{matrix}\right),$

where $B_{1}$ is defined in (4.13), so that explicitly we have

$\displaystyle(J^{\prime}_{1})^{22}=-(c_{3}L_{12(1)}+c_{2}L_{12(t)})+\left(c_{7}-\frac{c_{3}c_{5}}{c_{2}}\right)D_{1}+c_{8}D_{2}.$

The corresponding Hamiltonian density is a particular case of the density
$H^{\prime}_{0}$ from (8.19) at $c_{1}=0$

$\displaystyle
H^{\prime}_{0}=v\left\\{u_{2}-\frac{c_{5}}{c_{2}}z_{1}+\frac{c_{6}}{c_{3}}z_{2}\right\\}+s_{0}v+\frac{c_{3}}{2c_{2}}u_{1}u_{2}+\frac{c_{3}c_{5}}{c_{2}^{2}}u,$

where the integrability condition $c_{3}(c_{2}c_{4}+c_{3})=c_{2}^{2}c_{6}$ has
been used and an additional constraint has the form

$\displaystyle
c_{9}=\frac{1}{c_{2}^{2}}\\{c_{8}(c_{2}c_{4}+c_{3})+c_{5}(c_{3}c_{5}-c_{2}c_{7})\\}.$

The second case $c_{2}=0$, $c_{1}\neq 0$ is treated as a particular case of
the first bi-Hamiltonian family, so we put $c_{2}=0$ in the second Hamiltonian
operator $J_{1}$ from (8.1)

$\displaystyle
J_{1}=\left(\begin{matrix}-D_{1}^{-1}&\displaystyle\frac{c_{3}}{c_{1}}D_{1}^{-1}D_{2}\\\
-\displaystyle\frac{c_{3}}{c_{1}}D_{1}^{-1}D_{2}&\displaystyle\frac{c_{3}^{2}}{c_{1}^{2}}D_{1}^{-1}D_{2}^{2}+\frac{1}{c_{1}}B_{1}\end{matrix}\right),$

where $B_{1}$ is defined in (4.15), so that we have explicitly

$\displaystyle(J_{1})^{22}=\frac{c_{3}^{2}}{c_{1}^{2}}D_{1}^{-1}D_{2}^{2}+c_{3}L_{12(2)}-c_{1}L_{12(t)})+c_{6}D_{1}+\left(c_{7}+\frac{c_{3}c_{4}}{c_{1}}\right)D_{2}.$

The corresponding Hamiltonian density is a particular case of the density
$H_{0}$ from (8.4) at $c_{2}=0$

$\displaystyle
H_{0}=v\left(u_{1}+\frac{c_{4}}{c_{1}}z_{2}+\frac{c_{8}}{c_{3}}z_{1}+s_{0}\right)-\frac{c_{3}}{2c_{1}}u_{1}u_{2}+\frac{c_{3}c_{4}}{c_{1}^{2}}u,$

where integrability condition (4.14) has been used, with the additional
constraint

$\displaystyle
c_{9}=\frac{1}{c_{1}^{2}}\big{\\{}c_{6}(c_{3}-c_{1}c_{5})+c_{4}(c_{1}c_{7}+c_{3}c_{4})\big{\\}}.$

###  8.3 Special case $\boldsymbol{c_{3}=0}$

The second Hamiltonian operator is a composition $J_{1}=RJ_{0}$ of the
recursion operator (7.10) for this special case with the first Hamiltonian
operator $J_{0}$. Multiplication of these two matrices yields

$\displaystyle J_{1}=\left(\begin{matrix}c_{2}\nabla_{c}^{-1}&0\\\
0&\displaystyle\frac{1}{c_{1}}B_{1}\end{matrix}\right),$

where

$\displaystyle
J_{1}^{22}=-c_{2}\left(v_{1}D_{2}-v_{2}D_{1}+\frac{c_{6}}{c_{1}}D_{1}+\frac{c_{8}}{c_{2}}D_{2}\right)=\frac{1}{c_{1}}B_{1}$

with $B_{1}$ defined in (4.17).

The remaining problem is to obtain the Hamiltonian density $H_{0}$
corresponding to the second Hamiltonian operator $J_{1}$ according to (5.1)

$\displaystyle J_{1}\left(\begin{matrix}\delta_{u}H_{0}\\\
\delta_{v}H_{0}\end{matrix}\right)=\left(\begin{matrix}v\\\
v_{t}\end{matrix}\right).$ (8.26)

######  Proposition 8.3.

Bi-Hamiltonian representation (8.26) of the system (5.1) is valid under the
constraint

$\displaystyle c_{9}=\frac{(c_{1}c_{4}c_{8}-c_{2}c_{5}c_{6})}{c_{1}c_{2}}$
(8.27)

with the following Hamiltonian density

$\displaystyle
H_{0}=\frac{v}{c_{2}}\\{-\nabla_{c}(u)+c_{5}z_{1}-c_{4}z_{2}+s_{0}\\}.$ (8.28)

###### Proof.

The equation following from the first row of the matrix equation (8.26) is
$c_{2}\nabla_{c}^{-1}\delta_{u}H_{0}=v$ or equivalently

$\displaystyle\delta_{u}H_{0}=\frac{1}{c_{2}}\nabla_{c}(v)=\frac{c_{1}}{c_{2}}v_{1}+v_{2}.$
(8.29)

The second row of (8.26) yields

$\displaystyle\frac{1}{c_{1}}B_{1}(\delta_{v}H_{0})=v_{t}.$ (8.30)

We assume again a linear dependence of $H_{0}$ on $v$

$\displaystyle H_{0}=b[u]v+c[u]\quad\Longrightarrow\quad\delta_{v}H_{0}=b[u],$
(8.31)

where $b$ and $c$ may depend on $u$ and its partial derivatives.

Using (8.31) and also using the definition of $B_{1}$ and $v_{t}$, given by
(5.1), on the left-hand side and right-hand side of (8.30), respectively, we
obtain

$\displaystyle
c_{2}(v_{1}D_{2}[b]-v_{2}D_{1}[b])+\frac{c_{2}c_{6}}{c_{1}}D_{1}[b]+c_{8}D_{2}[b]$
$\displaystyle\qquad{}=-\big{\\{}v_{1}\nabla_{c}(u_{2})-v_{2}\nabla_{c}(u_{1})+c_{4}v_{1}+c_{5}v_{2}+c_{6}u_{11}+c_{7}u_{12}+c_{8}u_{22}+c_{9}\big{\\}}.$
(8.32)

Splitting this equation with respect to $v_{1}$ and $v_{2}$ we obtain two
equations

$\displaystyle
D_{1}[b]=-\frac{1}{c_{2}}\nabla_{c}(u_{1})+\frac{c_{5}}{c_{2}},\qquad
D_{2}[b]=-\frac{1}{c_{2}}\nabla_{c}(u_{2})-\frac{c_{4}}{c_{2}},$

which are integrated in the form

$\displaystyle
b=\frac{1}{c_{2}}\big{\\{}{-}\nabla_{c}(u)+c_{5}z_{1}-c_{4}z_{2}+s_{0}\big{\\}}.$
(8.33)

Using (8.33), we check that all the remaining nonconstant terms in (8.32)
cancel and we obtain an additional relation (8.27) between coefficients.
Plugging the expression (8.33) for $b$ into (8.31) for $H_{0}$ we obtain

$\displaystyle
H_{0}=\frac{v}{c_{2}}\big{\\{}{-}\nabla_{c}(u)+c_{5}z_{1}-c_{4}z_{2}+s_{0}\big{\\}}+c[u].$
(8.34)

Computing $\delta_{u}H_{0}$ from (8.34) and plugging the result into (8.29),
we conclude that $\delta_{u}c[u]=0$ and hence we may choose $c[u]=0$. Thus, we
end up with the formula (8.28). ∎

##  9 Conclusion

We have shown that all the Euler–Lagrange equations of the evolutionary Hirota
type in $(2+1)$ dimensions have the symplectic Monge–Ampère form. The symmetry
condition for such an equation is converted to a skew-factorized form. Then
recursion relations and Lax pairs are obtained as immediate consequences of
this representation. We have converted the equation into a two-component
evolutionary form and obtained Lagrangian and recursion operator for this two-
component system. The Lagrangian is degenerate because the momenta cannot be
inverted for the velocities. Applying to this degenerate Lagrangian the
Dirac’s theory of constraints, we have obtained a symplectic operator and its
inverse, the latter being a Hamiltonian operator $J_{0}$. We have found the
corresponding Hamiltonian density $H_{1}$, thus presenting our system in a
Hamiltonian form. Composing the recursion operator $R$ with $J_{0}$, we have
obtained the second Hamiltonian operator $J_{1}=RJ_{0}$. We have found the
Hamiltonian density $H_{0}$ corresponding to $J_{1}$ and thereby obtained bi-
Hamiltonian representation of our system under one constraint on the
coefficients, which is additional to the integrability condition. Thus, we end
up with a seven-parameter class of bi-Hamiltonian systems in $(2+1)$
dimensions.

Representation of symmetry condition in an alternative skew-factorized form
produces another seven-parameter family of bi-Hamiltonian systems which
satisfy another constraint, additional to the integrability condition. This
constraint is different from the additional constraint for the first family of
bi-Hamiltonian systems. If we require both constraints to coincide, then,
accounting also for the integrability condition, we end up with a six-
parameter family of tri-Hamiltonian systems.

We apply similar methods to Hirota equations in $3+1$ dimensions with similar
results. This work is in progress and will be published elsewhere.

### Acknowledgments

The authors are grateful to an anonymous referee for important remarks which
contributed to the improvement of our paper. The research of M.B. Sheftel is
partly supported by the research grant from Boğaziçi University Scientific
Research Fund (BAP), research project No. 11643.

## References

  * [1] Dirac P.A.M., Lectures on quantum mechanics, Belfer Graduate School of Science Monographs Series, Vol. 2, Belfer Graduate School of Science, New York, 1967.
  * [2] Doubrov B., Ferapontov E.V., On the integrability of symplectic Monge–Ampère equations, J. Geom. Phys. 60 (2010), 1604–1616, arXiv:0910.3407.
  * [3] Ferapontov E.V., Hadjikos L., Khusnutdinova K.R., Integrable equations of the dispersionless Hirota type and hypersurfaces in the Lagrangian Grassmannian, Int. Math. Res. Not. 2010 (2010), 496–535, arXiv:0705.1774.
  * [4] Ferapontov E.V., Khusnutdinova K.R., Hydrodynamic reductions of multidimensional dispersionless PDEs: the test for integrability, J. Math. Phys. 45 (2004), 2365–2377, nlin.SI/0312015.
  * [5] Ferapontov E.V., Khusnutdinova K.R., On the integrability of $(2+1)$-dimensional quasilinear systems, Comm. Math. Phys. 248 (2004), 187–206, nlin.SI/0305044.
  * [6] Ferapontov E.V., Kruglikov B., Novikov V., Integrability of dispersionless Hirota type equations in 4D and the symplectic Monge–Ampère property, arXiv:1707.08070.
  * [7] Fuchssteiner B., Fokas A.S., Symplectic structures, their Bäcklund transformations and hereditary symmetries, Phys. D 4 (1981), 47–66.
  * [8] Guthrie G.A., Recursion operators and non-local symmetries, Proc. Roy. Soc. London Ser. A 446 (1994), 107–114.
  * [9] Krasil’shchik J., Verbovetsky A., Geometry of jet spaces and integrable systems, J. Geom. Phys. 61 (2011), 1633–1674, arXiv:1002.0077.
  * [10] Magri F., A simple model of the integrable Hamiltonian equation, J. Math. Phys. 19 (1978), 1156–1162.
  * [11] Magri F., A geometrical approach to the nonlinear solvable equations, in Nonlinear Evolution Equations and Dynamical Systems (Proc. Meeting, Univ. Lecce, Lecce, 1979), Lecture Notes in Phys., Vol. 120, Springer, Berlin-New York, 1980, 233–263.
  * [12] Marvan M., Another look on recursion operators, in Differential Geometry and Applications (Brno, 1995), Masaryk University, Brno, 1996, 393–402.
  * [13] Neyzi F., Nutku Y., Sheftel M.B., Multi-Hamiltonian structure of Plebanski’s second heavenly equation, J. Phys. A: Math. Gen. 38 (2005), 8473–8485, nlin.SI/0505030.
  * [14] Olver P.J., Applications of Lie groups to differential equations, Graduate Texts in Mathematics, Vol. 107, Springer-Verlag, New York, 1986\.
  * [15] Papachristou C.J., Potential symmetries for self-dual gauge fields, Phys. Lett. A 145 (1990), 250–254.
  * [16] Sergyeyev A., A simple construction of recursion operators for multidimensional dispersionless integrable systems, J. Math. Anal. Appl. 454 (2017), 468–480, arXiv:1501.01955.
  * [17] Sheftel M.B., Yazıcı D., Recursion operators and tri-Hamiltonian structure of the first heavenly equation of Plebański, SIGMA 12 (2016), 091, 17 pages, arXiv:1605.07770.
  * [18] Sheftel M.B., Yazıcı D., Malykh A.A., Recursion operators and bi-Hamiltonian structure of the general heavenly equation, J. Geom. Phys. 116 (2017), 124–139, arXiv:1510.03666.