# Temperature-reducing shocks in optically-thin radiative MHD - analytical and
numerical results

B. Snow<EMAIL_ADDRESS>University of Exeter

###### Abstract

Shocks are often invoked as heating mechanisms in astrophysical systems, with
both adiabatic compression and dissipative heating that leading to temperature
increases. Whilst shocks are reasonably well understood for ideal
magnetohydrodynamic (MHD) systems, in many astrophysical plasmas, radiation is
an important phenomena, which can allow energy to leave the system. As such,
energy becomes non-conservative which can fundamentally change the behaviour
of shocks. The energy emitted through optically-thin radiation post-shock can
exceed the thermal energy increase, resulting in shocks that reduce the
temperature of the medium, i.e., cooling shocks that have a net decrease in
temperature across the interface. In this paper, semi-analytical solutions for
radiative shocks are derived to demonstrate that both cooling (temperature
decreasing) and heating (temperature increasing) shock solutions are possible
in radiative MHD. Numerical simulations of magnetic reconnection with
optically-thin radiative losses also yield both heating and cooling shocks in
roughly equal abundances. The detected cooling shocks feature a significantly
lower pressure jump across the shock than their heating counterparts. The
compression at the shock front leads to locally-enhanced radiative losses,
resulting in significant cooling within a few grid cells in the upstream and
downstream directions. The presence of temperature-reducing (cooling) shocks
is critical in determining the thermal evolution, and heating or cooling,
across a wealth of radiative astrophysical plasmas.

††preprint: AIP/123-QED

##  I Introduction

Shocks are a universal feature of astrophysical systems, for example
supernovae Weaver (1976); Chen  _et al._ (2022), bow shocks in stellar wind
Wilkin (1996), accretion disk around black holes Das (2002), magnetic-
reconnection driven flares (Petschek, 1964), umbral flashes above solar
sunspots Beckers and Tallant (1969), to name a few. Typically, shocks are
considered to be near-discontinuous jumps in properties of the medium, and are
intrinsically linked to heating though adiabatic compression of the medium and
non-adiabatic dissipation (which is efficient due to the sharp interface
between the upstream and downstream media). As such, shocks are often studied
for their role in energy redistribution and heating of astrophysical systems.

In ideal (non-radiative) magnetohydrodynamic (MHD) systems, the role of shocks
is reasonably well understood with analytical solution that demonstrate the
necessity for temperature increases across the shock Goedbloed, Keppens, and
Poedts (2010). However, often astrophysical systems require extensions to the
ideal MHD model to accurately capture the salient physics. A common extension
of the MHD equations is radiative MHD, which uses a temperature-dependent loss
function to approximate the energy radiated from an optically-thin plasma by
photons. The optically-thin radiative MHD model has been widely applied to
study the solar atmosphere and astrophysical systems (Gudiksen  _et al._ ,
2011; Forbes and Malherbe, 1991; Avara, McKinney, and Reynolds, 2016; Bingert
and Peter, 2011; Ni  _et al._ , 2015; Rozyczka, 1985). The exact nature of the
loss function depends on the abundances and relevant processes, and even the
numerical techniques used to calculate the loss function Hitomi Collaboration
_et al._ (2018).

The simplest inclusion of optically-thin radiative losses involves adding a
temperature and density dependent source term to the energy equation. However,
the inclusion of source terms in the underlying equations allows for
significant departures from the ideal MHD for stable solutions to shocks. For
shocks in the low-temperature ($T\lessapprox 10^{4}$K) interstellar medium,
molecular disassociation is known to be an import process that can lead to
strong cooling and effectively isothermal shocks (Draine, Roberge, and
Dalgarno, 1983; Draine and McKee, 1993). Radiative cooling has been shown to
important in determining the temperature profile of stellar shocks (Fokin,
Massacrier, and Gillet, 2004), producing less heating than would be predicted
by MHD alone. Many analytical studies have been performed to determine the
structure and temperature profile of radiative shocks finding weakly heating
and iosthermal solutions (Weaver, 1976; Das, 2002; Lesaffre  _et al._ , 2004).
However, an additional, unstudied solution class exists whereby the radiative
losses within the shock exceed the adiabatic heating leading to shocks that
have net temperature reductions in the postshock region, i.e., shocks that
cool the plasma. In this paper, such shocks that have a net decrease in
temperature across the interface are referred to as ’cooling shocks’.

This paper demonstrates the existence of cooling solutions for shocks in
radiative MHD using realistic radiative loss curves for optically-thin losses
in solar and astrophysical systems. Firstly, the theory is presented, along
with semi-analytic solutions to the radiative MHD equations. Shocks are then
analysed in a numerical simulation of magnetic reconnection that is tearing
unstable. Shocks are detected using the ShockID methodology (Snow  _et al._ ,
2021) and the local properties perpendicular to the shock front analysed to
identify both heating and cooling shock solutions. The presence of cooling
shocks means that the thermal contribution of shocks in radiative plasmas
requires careful consideration.

##  II Methods

The underlying equations for a compressible, radiative magnetohydrodynamic
(MHD) plasma are given by:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\textbf{v})=0,$
(1) $\displaystyle\frac{\partial}{\partial
t}(\rho\textbf{v})+\nabla\cdot\left(\rho\textbf{v}\textbf{v}+P\textbf{I}-\textbf{B
B}+\frac{\textbf{B}^{2}}{2}\textbf{I}\right)=0,$ (2)
$\displaystyle\frac{\partial}{\partial
t}\left(e+\frac{\textbf{B}^{2}}{2}\right)+\nabla\cdot\left[\textbf{v}(e+P)-(\textbf{v}\times\textbf{B})\times\textbf{B}\right]=-\rho^{2}\Lambda(T)+\Phi_{H},$
(3) $\displaystyle\frac{\partial\textbf{B}}{\partial
t}-\nabla\times(\textbf{v}\times\textbf{B})=0,$ (4)
$\displaystyle\nabla\cdot\textbf{B}=0,$ (5) $\displaystyle
e=\frac{P}{\gamma-1}+\frac{1}{2}\rho v^{2},$ (6)

for density $\rho$, velocity $\textbf{v}=\left[v_{x},v_{y},v_{z}\right]$,
magnetic field $\textbf{B}=\left[B_{x},B_{y},B_{z}\right]$ and pressure $P$.
The temperature is assumed to obey the ideal gas law $T=\frac{\gamma
P}{\rho}$. Optically thin radiative losses are assumed to be of the form
$\rho^{2}\Lambda(T)$, where $\Lambda(T)$ is a temperature dependent function.
The exact form of $\Lambda(T)$ often depends on the abundances of
atomic/molecular species, and the relevant processes that result in photon
emission. Included in the equations is a heating term $\Phi_{H}$ which is
included as an ad-hoc balance to the radiative losses such that equilibrium
can exist initially. The heating term $\Phi_{H}$ is usually assumed to be
constant in time since the heating processes typically occur on far longer
timescales than the radiative cooling. In this paper, $\Phi_{H}$ will be
assumed to be constant and uniform.

###  II.1 Rankine-Hugoniot jump conditions (no losses)

In this section, the analytic shock jump conditions are presented for an ideal
plasma (with no radiatve losses). For the MHD equations, the shock jump
conditions can be derived by considering the equations upstream (u) and
downstream (d) of the shock in the deHoffmann-Teller shock frame, where the
electric field is identically zero. The conservative nature of the MHD
equations means that for a steady-state shock in the deHoffmann-Teller frame,
mass, momentum and energy can be equated sufficiently upstream and downstream.
The MHD equations have no source terms and hence can be integrated across the
shock to become:

$\displaystyle\left[\rho v_{x}\right]^{u}_{d}=0,$ (7) $\displaystyle\left[\rho
v_{x}^{2}+P+\frac{B_{y}^{2}}{2}\right]^{u}_{d}=0,$ (8)
$\displaystyle\left[\rho v_{x}v_{y}-B_{x}B_{y}\right]^{u}_{d}=0,$ (9)
$\displaystyle\left[v_{x}\left(\frac{\gamma}{\gamma-1}P+\frac{1}{2}\rho
v^{2}\right)\right]^{u}_{d}=0,$ (10)
$\displaystyle\left[B_{x}\right]^{u}_{d}=0,$ (11)
$\displaystyle\left[v_{x}B_{y}-v_{y}B_{x}\right]^{u}_{d}=0,$ (12)

where this notation means that

$\displaystyle\left[Q\right]^{u}_{d}\equiv Q^{u}-Q^{d},$ (13)

for any quantity $Q$ evaluated u upstream (pre-shock) or d downstream (post-
shock). Shocks are essentially 1D structures since the gradients are greatest
in one direction and hence the jump conditions can be given in terms of along
the shock direction (here this is the $x$-direction) or perpendicular to it
($y-$direction).

From these equations, the jump conditions for primitive plasma properties can
be derivedGoedbloed, Keppens, and Poedts (2010); Snow and Hillier (2021) as:

$\displaystyle\frac{\rho^{d}}{\rho^{u}}=r$ (14)
$\displaystyle\frac{v_{x}^{d}}{v_{x}^{u}}=\frac{1}{r}$ (15)
$\displaystyle\frac{v_{y}^{d}}{v_{y}^{u}}=\frac{A_{x}^{u2}-1}{A_{x}^{u2}-r}$
(16) $\displaystyle B_{x}=\mbox{const.}$ (17)
$\displaystyle\frac{B_{y}^{d}}{B_{y}^{u}}=r\frac{A_{x}^{u2}-1}{A_{x}^{u2}-r}$
(18)
$\displaystyle\frac{P^{d}}{P^{u}}=1+\frac{2}{\beta(1+\tan^{2}\theta)}\frac{A_{x}^{u2}}{r}\left[r-1+\frac{r\tan^{2}\theta}{2A_{x}^{2u}}\left(1-r^{2}\left(\frac{A_{x}^{u2}-1}{A_{x}^{u2}-r}\right)^{2}\right)\right]$
(19)

where $r$ is the compression across the shock,
$A_{x}^{u2}=\frac{\rho^{u}v_{x}^{u2}}{B_{x}^{2}}$ is the inflow (upstream)
Alfvén Mach number squared, $\theta$ is the upstream angle of the magnetic
field, and $\beta$ is the upstream plasma-beta.

These ideal MHD jump conditions can be combined with the energy jump relation
(Equation 10) to be reduced to a single equation(Hau and Sonnerup, 1989) that
determines all possible stable shock solutions for a given upstream
plasma-$\beta$ and angle of the magnetic field $\theta$:

$\displaystyle
A_{x}^{\text{u}2}=\frac{A_{x}^{\text{d}2}\left(\frac{\gamma-1}{\gamma}\left(\frac{\gamma+1}{\gamma-1}-\tan^{2}\theta\right)\left(A_{x}^{\text{d}2}-1\right)^{2}+\tan^{2}\theta\left(\frac{\gamma-1}{\gamma}A_{x}^{\text{d}2}-1\right)\left(A_{x}^{\text{d}2}-2\right)\right)-\frac{\beta}{\cos^{2}\theta}\left(A_{x}^{\text{d}2}-1\right)^{2}}{\frac{\gamma-1}{\gamma}\frac{\left(A_{x}^{\text{d}2}-1\right)^{2}}{\cos^{2}\theta^{\text{u}}}-A_{x}^{\text{d}2}\tan^{2}\theta^{\text{u}}\left(\frac{\gamma-1}{\gamma}A_{x}^{\text{d}2}-1\right)}$
(20)

These equations are valid for any conservative set of MHD equations, i.e., a
set of equations with no sources or sinks of mass/momentum/energy. However,
the shock jumps for non-conservative equations are of interest, specifically
the optically-thin MHD model, where the radiative losses (and heating term)
act as a sink (source) of energy. For such a model, a different approach is
required.

The MHD (no losses) model has neglected the role of viscosity and resistivity
in the underlying equations which can lead to heating. However, by defining
the upstream and downstream locations to be sufficiently far from the shock,
where the gradients of the magnetic and velocity fields are zero, the
contribution of these terms vanishes, and the solution reduces to the ideal
MHD solution given above. Viscosity and resistivity do however affect the
substructure of the shock (Hau and Sonnerup, 1989).

###  II.2 Radiative shock jumps

Figure 1: Panel (a) (blue line) Normalised radiative loss curve using CHIANTI
v9 (Dere  _et al._ , 1997, 2019) with the default abundance file. The blue and
green vertical lines show the upstream and possible downstream temperatures
for a shock with $r=2$, as discussed towards the end of II.2. The orange line
is the modified version of the loss curve used in the simulation of Section
III. Panels (b)-(e) Example radiative MHD solutions for reference upstream
values of $\beta=0.1,\theta=\pi/8,T^{u}=10^{5}$K compared to the MHD solution
(black line). The radiative MHD shock solutions are coloured by the
temperature jump, with red denoting heating solution, blue denoting cooling
solutions, and effectively isothermal solutions in green. Reference upstream
temperatures are $T^{u}=1.8\times 10^{4}(b),5.0\times 10^{4}(c),2.3\times
10^{5}(d),5\times 10^{5}(e),10^{6}(f)$

Radiative losses introduce a non-conservative term in the energy equation that
fundamentally change the stable conditions across a shock. Applying the same
methodology as before (putting the underlying equations in the frame of
reference of the shock and assuming a steady-state solution) yields the
following:

$\displaystyle\nabla\cdot(\rho\textbf{v})=0,$ (21)
$\displaystyle\nabla\cdot\left(\rho\textbf{v}\textbf{v}+P\textbf{I}-\textbf{B
B}+\frac{\textbf{B}^{2}}{2}\textbf{I}\right)=0,$ (22)
$\displaystyle\nabla\cdot\left[\textbf{v}(e+P)-(\textbf{v}\times\textbf{B})\times\textbf{B}\right]=\rho^{2}\Lambda(T)+\Phi_{H},$
(23) $\displaystyle\nabla\times(\textbf{v}\times\textbf{B})=0,$ (24)
$\displaystyle\nabla\cdot\textbf{B}=0.$ (25)

The majority of these equations can be directly integrated as before, since
there are no sources of mass or momentum, however the energy equation cannot
be easily integrated due to the source term. A different approach is required
to determine the stable shock jump solutions.

Following the approach of Snow and Hillier (2021), the energy equation can be
discarded and replaced with a condition on the equilibrium, namely that for a
simple equilibrium, one requires that the heating and the loss terms balance,
i.e.,:

$\displaystyle\rho^{2}\Lambda(T)-\Phi_{H}=0$ (26)

Such condition must be true sufficiently far from an isolated shock in both
the upstream or downstream directions:

$\displaystyle\rho^{u2}\Lambda(T^{u})-\Phi_{H}^{u}=\rho^{d2}\Lambda(T^{d})-\Phi_{H}^{d}$
(27)

A common approach is to assume a temporally constant heating term to balance
the losses at time $t=0$. The heating term is also assumed to be spatially-
constant, i.e., $\Phi_{H}^{u}=\Phi_{H}^{d}$. Under these assumptions, the
equilibrium can be simplified as

$\displaystyle\frac{\rho^{d2}}{\rho^{u2}}=r^{2}=\frac{\Lambda(T^{u})}{\Lambda(T^{d})}$
(28)

i.e., a compression of the medium requires a reduction in the losses. The
compression $r$ can also be expressed in terms of the upstream and downstream
Alfvén Mach number

$\displaystyle\frac{A_{x}^{d2}}{A_{x}^{u2}}=\frac{1}{r}$ (29)

and hence the change in Mach number across the shock also depends on the
equilibrium conditions permitted by the loss function.

The shape of the radiative loss curve varies greatly depending on the medium,
composition and physics of interest. Here, as an example, a solar-like
optically-thin radiative loss curve is calculated using CHIANTI v9 (Dere  _et
al._ , 1997, 2019) with the default abundance file. As a simple example,
consider a shock that has a compression of $r=2$ with an upstream temperature
at the peak to the radiative loss curve, i.e., $T^{u}\approx 230,000$ K. From
Equation 28, for an equilibrium solution to exist, the losses must decrease by
a factor of 4 across the shock. Using the example radiative loss curve in
Figure 1a, this gives 4 potential solutions shown by the intercepts of the
horizontal dashed line with the loss curve, which give temperatures of
$T^{d}\approx 14500,25100,39800,1580000$ K. Three of these solutions are
potential cooling solutions, where the temperature reduces across the shock.
To evaluate if these are shock solutions, these potential solutions need to be
analysed using the rest of the shock jump equations.

####  II.2.1 Semi-analytic shock solution

The pressure jump across the shock (Equation 19) and the ideal gas law
$T=\gamma P/\rho$ can be combined to give a temperature jump equation, which,
when coupled with the radiative equilibrium condition, yields a pair of
equations that can be solved numerically to find the stable shock jump
conditions for a given upstream temperature, angle of magnetic field, and
plasma-beta:

$\displaystyle\frac{T^{d}}{T^{u}}=\frac{A^{d2}_{x}}{A_{x}^{u2}}\left[1+\frac{2}{\beta(1+\tan^{2}\theta)}\left(A_{x}^{u2}-A_{x}^{d2}+\frac{\tan^{2}\theta}{2}\left(1-\left(\frac{A_{x}^{u2}-1}{A_{x}^{d2}-1}\right)^{2}\right)\right)\right]$
(30)
$\displaystyle\frac{\Lambda(T^{d})}{\Lambda(T^{u})}=\left(\frac{A_{x}^{d2}}{A_{x}^{u2}}\right)^{2}$
(31)

Note that the pressure jump (Equation 19) arises without using the energy
equation, with the pressure term coming from the momentum equation. The
trivial solution where $A_{x}^{u2}=A_{x}^{d2}$, i.e., no shock, is satisfied
for no jump in temperature $T_{u}=T^{d}$, as expected. Other non-trivial
solutions exist depending on the shape of the radiative loss curve.

Numerically solving Equations 30-31 for a choice of upstream parameters as
$\beta=0.1,\theta=\pi/8$ for different upstream temperatures is shown in
Figure 1b-e for the CHIANTI cooling curve. Solutions are sought in the range
$A^{d2}<2,A^{u2}<4$ with a compression ratio in the range $1.01\leq r\leq
500$. The solution is compared to the ideal MHD solution given by Equation 20.
The shaded region of Figures 1b-f corresponds to a compression $r<1$ and thus
does not obey the definition of a shock, thus solutions are not sought in this
area. The radiative solution has both heating (red) and cooling (blue)
solutions, as well as a few effectively-isothermal shocks where the
temperature jump is defined as $0.9<T^{d}/T^{u}<1.1$.

The non-trivial ideal MHD solution is a cubic that intersects the linear
solution at the slow, Alfvén, and fast speeds. The radiative MHD simulation is
far less smooth than the ideal solution and cannot be described by a simple
(low-order) polynomial. Only solutions with a compression in the range
$1.01\leq r\leq 500$ are found by the solver, however this neglects some
highly compressible cooling solutions. The CHIANTI cooling curve has zero
losses below $T=10^{4}$ K and thus a cooling solution can be found for any
temperature above this, however the compression becomes increasingly large for
small $A^{d2}$; Equation 28 shows that a big jump in the losses requires a
large compression.

For the 5 tested upstream temperatures, all support cooling shock solutions
which follow a roughly cubic shape, see the blue line in Figures 1b-e. The
heating (i.e., temperature increasing solutions) have a far more varied shape.
For an upstream temperature of $T^{u}=1.8\times 10^{4}$K, a heating solution
is found that is close to the trivial solution where $A^{u2}=A^{d2}$, i.e., no
shock, Figure 1b. This implies that the heating solutions at this temperature
have a small increase in density across the interface. For $T^{u}=5\times
10^{4}$K, no heating solutions are found, Figure 1c. The heating solution for
reference upstream temperatures of $T=2.3\times 10^{5}$ K and $T=10^{6}$ are
reasonably similar to the MHD solution. Note that this only implies that the
compression across the shock is similar, not necessarily temperature. Overall,
Figures 1b-e show the sensitivity of the stable shock solutions relative to
the upstream temperature.

The solution in Figures 1b-e are only valid for the specific reference values
of $\beta=0.1,\theta=\pi/8,T^{u}$. The full parameter space is not presented
in this paper and different values of $\beta,\theta,T^{u}$ will produce
different sets of heating and cooling solutions. Appendix A shows a selection
of shock solutions for the CHIANTI cooling curve using a fixed temperature and
varying the plasma$-\beta$ and $\theta$.

###  II.3 Hazy Cooling Curve

Figure 2: (a) Hazy cooling curve normalised to unity at the peak. The
radiative shock solutions are for usptream temperatures of $10^{3}$K (b),
$10^{4}$K (c), $10^{5}$K (d), $10^{6}$K (e), and $10^{7}$K. The upstream
temperatures correspond to the vertical lines in panel (a).

The presence of cooling shock solutions is not unique to the CHIANTI cooling
curve used in the previous section. The CHIANTI cooling curve is specific to
the solar atmosphere, and other astrophysical systems have different
abundances and physical process, hence different temperature dependence of the
cooling curve. A common astrophysical loss function is the Hazy cooling curve
that is included in the default build of Cloudy version C17 Ferland  _et al._
(2017), shown in Figure 2a. The Hazy cooling curve assumes a gas that is in
collisional ionisation equilibrium with comic ray heating included for the low
temperature regime and solar-like abundances.

As before, Equations 30-31 can be numerically solved using the Hazy cooling
curve as the loss function, and specifying the upstream temperature,
plasma-$\beta$ and angle of the magnetic field $\theta$. Again,
plasma-$\beta=0.1$ and $\theta=\pi/8$ are held constant and the stable shock
solutions are calculated for different upstream temperatures, namely
$T^{u}=10^{3},10^{4},10^{5},10^{6},10^{7}$ K. The radiative shock solutions
coloured to indicate the temperature-increasing shocks (red), temperature-
decreasing shocks (blue), and essentially-isothermal (green) solutions, shown
in Figures 2b-e.

For upstream temperatures of less than $T^{u}\approx 10^{4}$, the radiative
loss function is lower than for any hotter temperature and thus only cooling
solutions are possible in order to satisfy Equation 31, and as such only
cooling solutions are found in Figures 2b,c. Similarly, only cooling solutions
are possible above the local loss minimum near the reference temperature of
$10^{7}$ since reductions of the radiative cooling curve only exist for cooler
post-shock temperatures. Both heating and cooling solutions are found for the
reference upstream temperatures of $T^{u}=10^{5},10^{6},10^{7}$ K, shown in
Figures 2d-f, respectively. As before, the presence of heating or cooling
solutions is very sensitive to the upstream temperature. However all tested
cases have temperature-reducing shocks as stable solutions of the optically-
thin radiative MHD equations.

###  II.4 Limitations of the analytical solutions

It should be noted that the analytical results presented are for a steady-
state, isolated shock and the jump conditions are valid for an undefined
distance in the upstream and downstream directions. Substructure may exist
between the upstream and downstream states which is not considered in the
analysis presented. Also, there is no consideration of the time taken to
obtain the equilibrium. For this, numerical simulations are useful in
determining the spatiotemporal scale on which cooling solutions can exist.

The presented shock solutions use a fixed plasma-$\beta$ and $\theta$ value,
with the upstream temperature varying. The full parameter space has not been
presented in this paper however some additional solutions are shown in
Appendix A for the CHIANTI cooling curve.

##  III Numerical Simulation of Tearing Instability

Figure 3: Numerical simulation of the tearing instability at time $t=40$
coloured by the temperature, with the identified slow-mode shocks highlighted
in red.

To investigate the heating and cooling solutions of shock, 2.5D numerical
simulations are performed using the (PIP) code (Hillier, Takasao, and
Nakamura, 2016) of a tearing-unstable reconnection current sheet with
optically-thin radiative losses. Resistivity is included in this model to
promote the instability without the reconnection relying on numerical
resistivity, specifically, the underlying equations for this simulation are:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\textbf{v})=0,$
(32) $\displaystyle\frac{\partial}{\partial
t}(\rho\textbf{v})+\nabla\cdot\left(\rho\textbf{v}\textbf{v}+P\textbf{I}-\textbf{B
B}+\frac{\textbf{B}^{2}}{2}\textbf{I}\right)=0,$ (33)
$\displaystyle\frac{\partial}{\partial
t}\left(e+\frac{\textbf{B}^{2}}{2}\right)+\nabla\cdot\left[\textbf{v}(e+P)-(\textbf{v}\times\textbf{B})\times\textbf{B}+\eta\left(\nabla\times\textbf{B}\right)\times\textbf{B}\right]=-\rho^{2}\Lambda(T),$
(34) $\displaystyle\frac{\partial\textbf{B}}{\partial
t}-\nabla\times(\textbf{v}\times\textbf{B}-\eta\nabla\times\textbf{B})=0,$
(35) $\displaystyle\nabla\cdot\textbf{B}=0,$ (36) $\displaystyle
e=\frac{P}{\gamma-1}+\frac{1}{2}\rho v^{2},$ (37)

where the resistivity $\eta=0.001$ is constant and uniform in the entire
domain. A modified form of the radiative loss curve is used, where the loss
function is smoothly reduced to zero for high temperatures, see Figure 1a. The
gradual decrease in losses was implemented to have a zero-gradient region at
extreme temperatures, however this is not a requirement of the code.

The initial condition is a pair of force-free Harris current sheets with a
superimposed perturbation:

$\displaystyle
B_{x}=B_{0}+B_{0}\tanh\left(\frac{y-\hat{y}}{w}\right)-B_{0}\tanh\left(\frac{y+\hat{y}}{w}\right)+B_{x,\mbox{pert}}$
(38) $\displaystyle B_{y}=B_{y,\mbox{pert}}$ (39) $\displaystyle
B_{z}=\frac{B_{0}}{\cosh\left(\frac{y-\hat{y}}{w}\right)}-\frac{B_{0}}{\cosh\left(\frac{y+\hat{y}}{w}\right)}$
(40)

where $\hat{y}=5$ is the $y-$location of the current sheet centres, $B_{0}=1$
is the reference magnetic field strength, and $w=0.1$ defines the width of the
current sheet. Since the initial conditions are force-free, the density and
pressure (and hence temperature) can be specified as constant, specifically as
$\rho=1$, $P=B_{0}^{2}\beta/2$, where the plasma-$\beta=0.01$. The pair of
current sheets allows periodic boundaries to be specified in all directions,
ensuring energy remains within the numerical box, with the exception of the
energy lost due to radiation.

The system is normalised according to an initial temperature of $T=10^{4}$K
across the domain. Since the radiative loss curve evaluated at this reference
temperature is identically zero, no heating term is required to balance the
initial conditions and hence $\Phi_{H}=0$. The radiative losses are normlised
such that the peak value is $\max(\Lambda(T))=0.01$, which, for the reference
density of $\rho=1$, means that thermal energy is removed from the system on
timescales of 100, which is longer than the simulation time. This is
calculated throughout the simulation so that dense regions, such as shocks,
feature enhanced losses depending on their local density and temperature.

The reconnection is triggered by specifying a localised pinch in the magnetic
field

$\displaystyle
B_{x,\mbox{pert}}=B_{0,\mbox{pert}}\frac{x}{2}e^{-\frac{(x-\hat{x})^{2}+(y-\hat{y})^{2}}{w/16}}-B_{0,\mbox{pert}}\frac{x}{2}e^{-\frac{(x+\hat{x})^{2}+(y+\hat{y})^{2}}{w/16}}$
(41) $\displaystyle
B_{y,\mbox{pert}}=B_{0,\mbox{pert}}\frac{y}{2}e^{-\frac{(x-\hat{x})^{2}+(y-\hat{y})^{2}}{w/16}}-B_{0,\mbox{pert}}\frac{y}{2}e^{-\frac{(x+\hat{x})^{2}+(y+\hat{y})^{2}}{w/16}}$
(42)

with reference amplitude $B_{0,\mbox{pert}}=0.1$ which is a small perturbation
centred on $\pm[\hat{x},\hat{y}]$ that decays with distance. The distance from
the current sheets is chosen to limit the interaction between the reconnection
regions on the timescale of the simulation. The simulation is resolved using
$4096\times 2048$ grid cells spanning $-20\leq x\leq 20$, $-10\leq y\leq 10$
($dx=dy\approx 0.01$) with a 1st order HLLD solver.

As the initial conditions evolve, the current sheet reconnects and becomes
tearing-unstable leading to a wealth of smaller plasmoids being generated, as
shown by the snapshot at time $t=40$ in Figure 3. The evolution of the tearing
instability is well discussed in the literature (Del Zanna  _et al._ , 2016;
Loureiro and Uzdensky, 2016). Here, the purpose is not to study the
intricacies of the tearing instability, but rather use is as a test bed for
identifying and studying shocks. Plasmoids generated by fast magnetic
reconnection are known to host a wealth of shocks (Zenitani, 2015) making this
simulation ideal for identifying the presence of heating and cooling shocks.

The system features a broad distribution of temperature that spans the full
extent of the cooling curve, as shown in Figure 3. The background temperature
remains stable at $T=10^{4}$, however the reconnection region features strong
increases to a few million Kelvin. The high temperatures produced are a result
of the low plasma beta used in the simulation, which provides a large
reservoir of magnetic energy that can be converted to thermal and kinetic
energy during magnetic reconnection.

Figure 4: Histogram of the temperature jump of the detected slow-mode shocks
in the tearing instability coloured.

###  III.1 Shock Detection

Shocks are detected and classified using the open-source ShockID routine (Snow
_et al._ , 2021) which has recently been ported to Python with additional
features such as parallisation. The shock detection method works by
identifying candidate shocks from locations where $\nabla\cdot\textbf{v}<0$,
i.e., converging flow, calculating the shock direction from the maximum
density gradient, interpolating the velocity, magnetic field, density, and
pressure along the shock direction, and finally using these values to classify
the shock transition. The result is a list of grid cell locations that satisfy
shock jump transitions. Figure 3 shows the detected slow-mode shocks at an
instance in time. It should be noted that the shock detection method is
approximate and only detects whether the candidate cell satisfies the shock
jump equations.

The majority of the detected shocks are categorised as slow-mode shocks and
centred around the plasmoids, which is to be expected from the known structure
of shocks around plasmoids (Zenitani, 2015). Here only the slow-mode shocks
are analysed in detail since these are the most numerous in the simulation.

Figure 5: Scatter plot of the compression vs the temperature jump for the
detected slow-mode shocks, coloured by the pressure jump Figure 6: Scatter
plot of the upstream temperature vs the downstream temperature for the
detected slow-mode shocks, coloured by the total losses along the 1D shock
width. The black line represents an isothermal temperature across the shock.

A histogram of the temperature jumps across the shocks is shown in Figure 4,
where ’across the shock’ is a distance of fewer than 5 grid cells in the
upstream and downstream direction (the exact length of the shock depends on
the density profile along the shock line and the maximum line length checked
is 5 cells either side of the shock for a total of 11 grid cells). Three broad
categories of shocks exist: heating shocks, cooling shocks, and near-
isothermal shocks. The near-isothermal shocks are determined as having a
temperature jump $0.9<=T^{d}/T^{u}<=1.1$. The temperature jump across the
shock is calculated local to the shock, as in within 5 grid cells in the
upstream and downstream direction. There are 2670 heating solutions, 1722
cooling solutions, and 357 near-isothermal shocks detected in the simulation.
Heating and cooling solutions are therefore abundant in the simulation. This
can be seen by cool ‘halos’ that form around the plasmoids in Figure 3.

The detected compression across the shocks is in the range $1<r<40$, with a
large range of compression in both the heating and cooling shocks, see Figure
5. In MHD, it is possible to determine a maximum compression ratio as a
function of the adiabatic index by considering an infinitely fast fast-mode
shock, which yeilds:

$\displaystyle r_{max}=\frac{\gamma+1}{\gamma-1}$ (43)

which, for $\gamma=5/3$ which is used in this paper, yields a maximum
compression of $r_{max}=4$. In the radiative simulation, slow-mode shocks are
detected with an order of magnitude more compression than is possible in ideal
MHD. The enhanced compression is seen in other simulations of shocks with non-
conservative energy equations to account for losses of macroscopic fluid
energy (Snow and Hillier, 2021; Snow, Druett, and Hillier, 2023).

The detected heating shocks are generally associated with a increase in
pressure across the shock, as shown in Figure 5. This is similar to MHD, where
a pressure jump is necessary for a shock. However, for the cooling shocks,
there are cases where the pressure decreases and increases. Cooling shocks can
feature a decrease in temperature despite an increase in pressure, as long as
the pressure increase is small compared to the density increase. Similarly
isothermal shock solutions can be obtained when the density and pressure jumps
are roughly equal.

A cooling shock that has a pressure jump less than unity presents an odd
situation where the pressure and density gradients act in opposite direction.
As such, shocks of this nature may become baroclinically unstable, however,
this is not detected in the simulation. Further study is needed on the long-
term stability of the cooling shocks with pressure reductions.

Figure 6 shows a scatter-plot of the upstream and downstream temperatures
coloured by the total losses along the 1D shock width. Both heating and
cooling shocks exist across a wide range of temperatures in the upstream and
downstream regions. This was expected from Section II.2 where a single
upstream temperature allows stable solution to exist for multiple downstream
temperatures, which can be either heating or cooling. The total losses along
the shock, calculated as

$\displaystyle L=\sum^{u}_{d}\rho^{2}\Lambda(T)$ (44)

gives the losses per unit area over the finite width of the shock. A trend
exist whereby the total losses are stronger for cooling solutions. This is
reasonably intuitive since for a cooling shock to exist, the energy lost due
to radiation must exceed the adiabatic temperature increase due to
compression. The timescale for the losses becomes significantly enhanced in
all the shocks due to the large increases in density across the interface. For
the initial reference density, losses occur on timescales of $10^{2}$ at the
peak of the cooling curve. Here, losses are seen to have timescales of around
$10^{-1}$, three orders of magnitude quicker than the reference values. As
such, radiative losses are very efficient in the shocks.

##  IV Discussion

This paper has proven both analytically and numerically that the optically-
thin radiative MHD equations can support both heating and cooling shock
solutions. The presence of cooling solutions depends on the radiative cooling
curve and the upstream temperature.

###  IV.1 Model assumptions

The underlying model used for this analysis is widely used across a number of
astrophysical plasmas. The results of this paper are only as valid as the
assumptions of the underlying model, as is true of any analytical or numerical
study.

A strong assumption here is that the heating function is constant. In Snow and
Hillier (2021), a related set of equations was studied where the heating
function was constant and energy loss was due to the ionisation energy
required to change the state of the atom. It was concluded that their
underlying model necessarily resulted in shocks that reduce the temperature of
the plasma. However, in a subsequent paper, the underlying model was extended
to include more realistic treatment of hydrogen and a time dependent heating
function revealed that the same type of shock (switch-off slow-mode) now
featured an increase in temperature postshock Snow, Druett, and Hillier
(2023).

The heating function should likely be dependent on the local properties of the
plasma, e.g., for photon heating, the energy absorbed by the plasma is
dependent on the optical depth and hence the density. At shock fronts, the
density increases greatly (up to a factor of 40 was found in the simulation)
and hence one may expect that the energy absorbed by photons would increase,
leading to locally enhanced heating in the post-shock region. Further work is
required to study the consequences of non-constant heating functions and how
this affects the post-shock temperature.

The analytical solutions assume steady-state shocks and only consider the jump
conditions an arbitrary distance from the shock. In dynamic systems, the time
taken to obtain a steady-state shock may be long and thus shocks in
observations and simulations may not strictly adhere to the developed
solutions. However, the simulation presented in this paper shows that the
cooling solution can be obtained over a relatively short distance and hence
provides evidence of cooling shocks forming on small lengthscales.

###  IV.2 Interstellar medium

Shocks in the interstellar medium are often studied using the isothermal
assumption since the losses within the shock are assumed to be strong enough
to balance the compressional temperature increase post-shock. This has been
used to explain the observed presence of molecules in the post-shock regions
that should have been disassociated by the MHD temperature predictions (Draine
and McKee, 1993). As such, a modelling approximation is to the set the
adiabatic index to unity such that the system is isothermal.

For typical temperatures of the interstellar medium ($T<10^{4}$ K), the stable
shock solutions using the Hazy cooling curve result in only cooling solutions,
see Figure 2. Cooling shocks could also explain the observed lack of molecular
disassociation. However, the cooling shock solutions are not possible to
obtain by changing the adiabative index alone, since this would require
$\gamma<1$ local to the shock. Setting $\gamma<1$ is problematic from a
thermodynamic standpoint. The reason why such departures from classical
thermodynamic laws occur for this model is that the problem becomes non-local,
with energy available to leave the system through the radiative losses.

###  IV.3 Relation to the thermal stability

The thermal instability occurs when the radiative losses cause a reduction in
the pressure, driving an inflow which enhances the temperature, allowing for
additional cooling. The thermal instability is critical for the formation of
cool condensations in the solar atmosphereField (1965); Sen and Keppens
(2022), planetary nebulaeHunter (1970), and the intergalactic medium Sánchez-
Salcedo, Vázquez-Semadeni, and Gazol (2002). It is possible that non-steady
cooling shocks could result in a post-shock condition that is unstable to the
thermal instability, potentially enhancing the formations of condensations.

The idea of shock-induced coronal condensations has been proposed in other
papers that use very similar underlying model equations, i.e., MHD with
optically-thin radiative losses (Forbes and Malherbe, 1991). The post-shock
temperature calculated in Forbes and Malherbe (1991) is based on the MHD
predicted model and thus shocks were assumed to heat, then gradually cool
post-shock over relatively long timescales. In this paper it has been
demonstrated that cooling solutions are possible and can occur over short
length-scales. As such, cooling shock solutions may be an efficient way to
induce condensations, which requires further study.

##  V Conclusion

In this paper, it has been shown both analytically and numerically that the
radiative MHD equations allow both heating and cooling solutions of shocks to
exist, where a cooling shock has a temperature decrease across the interface.
Introducing optically-thin radiative losses results in a non-conservative
energy equation, however the shock jump conditions can be closed through
consideration of an equilibrium state either side of the shock, resulting in a
pair of equations that can be numerically solved to find the steady-state
shock solutions for given upstream angle of magnetic field, plasma beta, and
temperature.

The analytical solution shows that radiative shocks have a strong departure
from the ideal MHD solution and multiple shock solutions exist that can be
heating or cooling. However, this is under the assumption that the shocks are
steady state and the jump conditions have reached an equilibrium over some
unknown distance. In the simulation data, the same behaviour is found, i.e.,
both heating and cooling solutions, when analysed over a finite and small
distance specified as 5 grid cells in the upstream and downstream direction.
This shows that cooling shocks can be obtained locally. Since the optically-
thin radiative losses considered in this paper are dependent on the density
squared, compression at the shock front can drastically increase rate at which
radiative losses remove energy.

The shape of the cooling curve is important in determining the presence of
heating and cooling solutions. Equation 28 relates the loss rate either side
of the shock to the compression:

$\displaystyle\frac{\rho^{d2}}{\rho^{u2}}=r^{2}=\frac{\Lambda(T^{u})}{\Lambda(T^{d})}$
(45)

For a hypothetical radiative loss curve that monotonically decreases with
temperature, a compression of the medium requires a temperature increase
across the interface. Similarly, a monotonically increasing loss curve would
only allow cooling shock solutions. Realistic radiative loss curves often have
a complex shape with multiple turning points and thus allow both heating and
cooling solutions to exist.

Only the optically-thin regime is studied here. For optically-thick plasmas,
the emitted photons can become reabsorbed and potentially lead to heating
processes. This process may become important for shocks, where the density
increases greatly, and thus the opacity may increase. As such, including
radiative heating processes in the underlying equations may provide extra
heating or even inhibit cooling shocks from forming. Further work is required
to study this effect.

Often, shocks are studied for their role in heating. However, for the
optically-thin raditive MHD model presented here, shocks can also provide
substantial cooling due to the locally enhanced radiative losses in the
finite-width of the shock. As such, the thermal contribution of shocks is far
less definitive when radiative losses are included than ideal MHD results
would imply.

###### Acknowledgements.

BS is supported by STFC research grant ST/V000659/1. I would thank Prof.
Andrew Hillier (University of Exeter) and Dr. Malcolm Druett (KU Leuven) for
support and feedback with this research.

## Data Availability Statement

The data that support the findings of this study are available from the BS
upon reasonable request. The (PIP) solver is freely available on GitHub:
https://github.com/AstroSnow/PIP; as is the ShockID code:
https://github.com/AstroSnow/shockid. The post-processing scripts, including
the codes required to find shock solutions for given upstream properties with
a user-specified loss function, are available here:
https://github.com/AstroSnow/RadiativeShockCool.

##  Appendix A Shock solutions for different plasma$-\beta$ and $\theta$

Figures 7 shows shock solutions for the CHIANTI cooling curve using different
plasma-$\beta$ and $\theta$ values for a constant upstream temperature of
$T^{u}=10^{5}$K. Solutions are sought in the compression range $1.01\leq r\leq
50$ for downstream Alfvén Mach numbers of $0<A^{d2}<2$.

Figure 7: Shock solutions for the CHIANTI cooling curve with a constant
upstream temperature of $T^{u}=10^{5}$K and varying plasma-$\beta$ and
$\theta$ values. Cooling solutions are shown in blue, and heating solutions in
red.

## References

  * Weaver (1976) T. A. Weaver, “The structure of supernova shock waves.” ApJS 32, 233–282 (1976).
  * Chen  _et al._ (2022) W. Chen, P. L. Kelly, M. Oguri, T. J. Broadhurst, J. M. Diego, N. Emami, A. V. Filippenko, T. L. Treu, and A. Zitrin, “Shock cooling of a red-supergiant supernova at redshift 3 in lensed images,” Nature (London) 611, 256–259 (2022), arXiv:2306.12985 [astro-ph.GA] .
  * Wilkin (1996) F. P. Wilkin, “Exact Analytic Solutions for Stellar Wind Bow Shocks,” ApJ 459, L31 (1996).
  * Das (2002) T. K. Das, “Generalized Shock Solutions for Hydrodynamic Black Hole Accretion,” Astrophys. J. 577, 880–892 (2002), arXiv:astro-ph/0212119 [astro-ph] .
  * Petschek (1964) H. E. Petschek, “Magnetic Field Annihilation,” in  _NASA Special Publication_ , Vol. 50 (1964) p. 425.
  * Beckers and Tallant (1969) J. M. Beckers and P. E. Tallant, “Chromospheric Inhomogeneities in Sunspot Umbrae,” Sol Phys 7, 351–365 (1969).
  * Goedbloed, Keppens, and Poedts (2010) J. P. Goedbloed, R. Keppens, and S. Poedts,  _Advanced Magnetohydrodynamics_ (2010).
  * Gudiksen  _et al._ (2011) B. V. Gudiksen, M. Carlsson, V. H. Hansteen, W. Hayek, J. Leenaarts, and J. Martínez-Sykora, “The stellar atmosphere simulation code Bifrost. Code description and validation,” A&A 531, A154 (2011), arXiv:1105.6306 [astro-ph.SR] .
  * Forbes and Malherbe (1991) T. G. Forbes and J. M. Malherbe, “A numerical simulation of magnetic reconnection and radiative cooling in line-tied current sheets,” Sol Phys 135, 361–391 (1991).
  * Avara, McKinney, and Reynolds (2016) M. J. Avara, J. C. McKinney, and C. S. Reynolds, “Efficiency of thin magnetically arrested discs around black holes,” MNRAS 462, 636–648 (2016), arXiv:1508.05323 [astro-ph.HE] .
  * Bingert and Peter (2011) S. Bingert and H. Peter, “Intermittent heating in the solar corona employing a 3D MHD model,” A&A 530, A112 (2011), arXiv:1103.6042 [physics.space-ph] .
  * Ni  _et al._ (2015) L. Ni, B. Kliem, J. Lin, and N. Wu, “Fast Magnetic Reconnection in the Solar Chromosphere Mediated by the Plasmoid Instability,” Astrophys. J. 799, 79 (2015), arXiv:1509.06895 [astro-ph.SR] .
  * Rozyczka (1985) M. Rozyczka, “Two-dimensional models of stellar wind bubbles. I - Numerical methods and their application to the investigation of outer shell instabilities,” A&A 143, 59–71 (1985).
  * Hitomi Collaboration  _et al._ (2018) Hitomi Collaboration, F. Aharonian, H. Akamatsu, F. Akimoto, S. W. Allen, L. Angelini, M. Audard, H. Awaki, M. Axelsson, A. Bamba, M. W. Bautz, R. Blandford, L. W. Brenneman, G. V. Brown, E. Bulbul, E. M. Cackett, M. Chernyakova, M. P. Chiao, P. S. Coppi, E. Costantini, J. de Plaa, C. P. de Vries, J.-W. den Herder, C. Done, T. Dotani, K. Ebisawa, M. E. Eckart, T. Enoto, Y. Ezoe, A. C. Fabian, C. Ferrigno, A. R. Foster, R. Fujimoto, Y. Fukazawa, A. Furuzawa, M. Galeazzi, L. C. Gallo, P. Gandhi, M. Giustini, A. Goldwurm, L. Gu, M. Guainazzi, Y. Haba, K. Hagino, K. Hamaguchi, I. M. Harrus, I. Hatsukade, K. Hayashi, T. Hayashi, K. Hayashida, N. Hell, J. S. Hiraga, A. Hornschemeier, A. Hoshino, J. P. Hughes, Y. Ichinohe, R. Iizuka, H. Inoue, Y. Inoue, M. Ishida, K. Ishikawa, Y. Ishisaki, M. Iwai, J. Kaastra, T. Kallman, T. Kamae, J. Kataoka, S. Katsuda, N. Kawai, R. L. Kelley, C. A. Kilbourne, T. Kitaguchi, S. Kitamoto, T. Kitayama, T. Kohmura, M. Kokubun, K. Koyama, S. Koyama, P. Kretschmar, H. A. Krimm, A. Kubota, H. Kunieda, P. Laurent, S.-H. Lee, M. A. Leutenegger, O. Limousin, M. Loewenstein, K. S. Long, D. Lumb, G. Madejski, Y. Maeda, D. Maier, K. Makishima, M. Markevitch, H. Matsumoto, K. Matsushita, D. McCammon, B. R. McNamara, M. Mehdipour, E. D. Miller, J. M. Miller, S. Mineshige, K. Mitsuda, I. Mitsuishi, T. Miyazawa, T. Mizuno, H. Mori, K. Mori, K. Mukai, H. Murakami, R. F. Mushotzky, T. Nakagawa, H. Nakajima, T. Nakamori, S. Nakashima, K. Nakazawa, K. K. Nobukawa, M. Nobukawa, H. Noda, H. Odaka, T. Ohashi, M. Ohno, T. Okajima, N. Ota, M. Ozaki, F. Paerels, S. Paltani, R. Petre, C. Pinto, F. S. Porter, K. Pottschmidt, C. S. Reynolds, S. Safi-Harb, S. Saito, K. Sakai, T. Sasaki, G. Sato, K. Sato, R. Sato, M. Sawada, N. Schartel, P. J. Serlemtsos, H. Seta, M. Shidatsu, A. Simionescu, R. K. Smith, Y. Soong, Ł. Stawarz, Y. Sugawara, S. Sugita, A. Szymkowiak, H. Tajima, H. Takahashi, T. Takahashi, S. Takeda, Y. Takei, T. Tamagawa, T. Tamura, T. Tanaka, Y. Tanaka, Y. T. Tanaka, M. S. Tashiro, Y. Tawara, Y. Terada, Y. Terashima, F. Tombesi, H. Tomida, Y. Tsuboi, M. Tsujimoto, H. Tsunemi, T. G. Tsuru, H. Uchida, H. Uchiyama, Y. Uchiyama, S. Ueda, Y. Ueda, S. Uno, C. M. Urry, E. Ursino, S. Watanabe, N. Werner, D. R. Wilkins, B. J. Williams, S. Yamada, H. Yamaguchi, K. Yamaoka, N. Y. Yamasaki, M. Yamauchi, S. Yamauchi, T. Yaqoob, Y. Yatsu, D. Yonetoku, I. Zhuravleva, A. Zoghbi, and A. J. J. Raassen, “Atomic data and spectral modeling constraints from high-resolution X-ray observations of the Perseus cluster with Hitomi,” PASJ 70, 12 (2018), arXiv:1712.05407 [astro-ph.HE] .
  * Draine, Roberge, and Dalgarno (1983) B. T. Draine, W. G. Roberge, and A. Dalgarno, “Magnetohydrodynamic shock waves in molecular clouds.” Astrophys. J. 264, 485–507 (1983).
  * Draine and McKee (1993) B. T. Draine and C. F. McKee, “Theory of interstellar shocks.” Ann Rev A Ap 31, 373–432 (1993).
  * Fokin, Massacrier, and Gillet (2004) A. B. Fokin, G. Massacrier, and D. Gillet, “Radiative cooling of shocked gas in stellar atmospheres. II. Self-consistent LTE shock wake model with Fe lines and H continua contributions,” A&A 420, 1047–1059 (2004).
  * Lesaffre  _et al._ (2004) P. Lesaffre, J. P. Chièze, S. Cabrit, and G. Pineau des Forêts, “Temporal evolution of magnetic molecular shocks. II. Analytics of the steady state and semi-analytical construction of intermediate ages,” A&A 427, 157–167 (2004), arXiv:astro-ph/0407588 [astro-ph] .
  * Snow  _et al._ (2021) B. Snow, A. Hillier, G. Murtas, and G. J. J. Botha, “Shock identification and classification in 2d magnetohydrodynamiccompressible turbulence—orszag–tang vortex,” Experimental Results 2, e35 (2021).
  * Snow and Hillier (2021) B. Snow and A. Hillier, “Collisional ionisation, recombination, and ionisation potential in two-fluid slow-mode shocks: Analytical and numerical results,” A&A 645, A81 (2021), arXiv:2010.06303 [astro-ph.SR] .
  * Hau and Sonnerup (1989) L. N. Hau and B. U. O. Sonnerup, “On the structure of resistive MHD intermediate shocks,” J. Geophys. Res 94, 6539–6551 (1989).
  * Dere  _et al._ (1997) K. P. Dere, E. Landi, H. E. Mason, B. C. Monsignori Fossi, and P. R. Young, “CHIANTI - an atomic database for emission lines,” A&AS 125, 149–173 (1997).
  * Dere  _et al._ (2019) K. P. Dere, G. Del Zanna, P. R. Young, E. Landi, and R. S. Sutherland, “CHIANTI—An Atomic Database for Emission Lines. XV. Version 9, Improvements for the X-Ray Satellite Lines,” ApJS 241, 22 (2019), arXiv:1902.05019 [astro-ph.SR] .
  * Ferland  _et al._ (2017) G. J. Ferland, M. Chatzikos, F. Guzmán, M. L. Lykins, P. A. M. van Hoof, R. J. R. Williams, N. P. Abel, N. R. Badnell, F. P. Keenan, R. L. Porter, and P. C. Stancil, “The 2017 Release Cloudy,” Revista Mexicana de Astronomía y Astrofísica 53, 385–438 (2017), arXiv:1705.10877 [astro-ph.GA] .
  * Hillier, Takasao, and Nakamura (2016) A. Hillier, S. Takasao, and N. Nakamura, “The formation and evolution of reconnection-driven, slow-mode shocks in a partially ionised plasma,” A&A 591, A112 (2016), arXiv:1602.01112 [astro-ph.SR] .
  * Del Zanna  _et al._ (2016) L. Del Zanna, E. Papini, S. Landi, M. Bugli, and N. Bucciantini, “Fast reconnection in relativistic plasmas: the magnetohydrodynamics tearing instability revisited,” MNRAS 460, 3753–3765 (2016), arXiv:1605.06331 [astro-ph.HE] .
  * Loureiro and Uzdensky (2016) N. F. Loureiro and D. A. Uzdensky, “Magnetic reconnection: from the Sweet-Parker model to stochastic plasmoid chains,” Plasma Physics and Controlled Fusion 58, 014021 (2016), arXiv:1507.07756 [physics.plasm-ph] .
  * Zenitani (2015) S. Zenitani, “Magnetohydrodynamic structure of a plasmoid in fast reconnection in low-beta plasmas: Shock-shock interactions,” Physics of Plasmas 22, 032114 (2015), arXiv:1503.01920 [astro-ph.SR] .
  * Snow, Druett, and Hillier (2023) B. Snow, M. K. Druett, and A. Hillier, “Partially ionized two-fluid shocks with collisional and radiative ionization and recombination - multilevel hydrogen model,” MNRAS 525, 4717–4734 (2023), arXiv:2308.12802 [astro-ph.SR] .
  * Field (1965) G. B. Field, “Thermal Instability.” Astrophys. J. 142, 531 (1965).
  * Sen and Keppens (2022) S. Sen and R. Keppens, “Thermally enhanced tearing in solar current sheets: Explosive reconnection with plasmoid-trapped condensations,” A&A 666, A28 (2022), arXiv:2208.04355 [astro-ph.SR] .
  * Hunter (1970) J. Hunter, James H., “Generalized Thermal Stability and its Application to the Interstellar Gas,” Astrophys. J. 161, 451 (1970).
  * Sánchez-Salcedo, Vázquez-Semadeni, and Gazol (2002) F. J. Sánchez-Salcedo, E. Vázquez-Semadeni, and A. Gazol, “The Nonlinear Development of the Thermal Instability in the Atomic Interstellar Medium and Its Interaction with Random Fluctuations,” Astrophys. J. 577, 768–788 (2002), arXiv:astro-ph/0203067 [astro-ph] .


# Two-component integrable extension
of general heavenly equation

Wojciech Kryński  Institute of Mathematics, Polish Academy of Sciences,
Śniadeckich 8, 00-656 Warsaw, Poland<EMAIL_ADDRESS>and  Artur Sergyeyev
Mathematical Institute, Silesian University in Opava, Na Rybníčku 1, 746 01
Opava, Czech Republic<EMAIL_ADDRESS>
###### Abstract.

We introduce an integrable two-component extension of the general heavenly
equation and prove that the solutions of this extension are in one-to-one
correspondence with 4-dimensional hyper-para-Hermitian metrics. Furthermore,
we demonstrate that if the metrics in question are hyper-para-Kähler, then our
system reduces to the general heavenly equation. We also present an infinite
hierarchy of nonlocal symmetries and a recursion operator for the system under
study.

##  1\. Introduction

The partial differential equation in four independent variables derived by
Schief in [24], known as the general heavenly equation, serves as a natural
framework for anti-self-dual gravity and has a number of important
applications in physics and integrable systems [6, 11, 15, 16, 23]. In the
present paper, we introduce and study its two-component integrable extension,
given below as system (1). We demonstrate that solutions of this new system
descend to a special class of metrics $[g]$ of neutral signature. These
metrics admit a triple of anti-commuting vector bundle endomorphisms
$I,J,K\colon TM\to TM$ such that $J$ is a complex structure, $I$ and $K$ are
para-complex structures,

$IK=-KI=J$

and

$g(X,Y)=-g(JX,JY)=g(IX,IY)=g(KX,KY),$

for any $X,Y\in T_{x}M$, $x\in M$. Following [10, Section 6.1], we shall refer
to the structures as hyper-para-Hermitian (see also [4, 5]). They are also
known as pseudo-hyperhermitian in [3, 9] or hyper-Hermitian of neutral
signature in [13]. The underlying algebra spanned by $(I,J,K)$ is called
split-quaternions [10] or para-quaternions [3], and it can be easily shown
that the triple itself uniquely determines the conformal class $[g]$, see [4]
and [13] for higher-dimensional generalizations. The hyper-para-Hermitian
structures are necessarily anti-self-dual and are characterized in the class
of anti-self-dual metrics by the condition that the corresponding twistor
space, defined as the space of $\alpha$-surfaces, fibers over a projective
line [9].

In this paper we work in a local coordinate system $(x^{1},x^{2},x^{3},x^{4})$
on manifold $M$. For any function $u\in C^{\infty}(M)$ we denote
$u_{i}=\partial_{x^{i}}u$ and $u_{ij}=\partial_{x^{i}}\partial_{x^{j}}u$. Our
main result is the following theorem giving a characterization of hyper-para-
Hermitian metrics on 4-manifolds.

######  Theorem 1.

Any hyper-para-Hermitian metric $[g]$ a 4-dimensional manifold can be locally
put in the form

$g=\omega^{1}\omega^{4}-\omega^{2}\omega^{3}$

where 1-forms $\omega^{i}$, $i=1,\ldots,4$ are defined as

$\displaystyle\omega^{1}=dx^{4}-\lambda_{4}\left(\frac{1}{\lambda_{2}-\lambda_{4}}\frac{v_{1}}{v_{4}}dx^{1}+\frac{1}{\lambda_{1}-\lambda_{4}}\frac{u_{2}}{u_{4}}dx^{2}\right)$
$\displaystyle\omega^{2}=\frac{1}{\lambda_{4}}dx^{4}-\frac{1}{\lambda_{4}}\left(\frac{\lambda_{2}}{\lambda_{2}-\lambda_{4}}\frac{v_{1}}{v_{4}}dx^{1}+\frac{\lambda_{1}}{\lambda_{1}-\lambda_{4}}\frac{u_{2}}{u_{4}}dx^{2}\right)$
$\displaystyle\omega^{3}=dx^{3}-\lambda_{3}\left(\frac{1}{\lambda_{2}-\lambda_{3}}\frac{v_{1}}{v_{3}}dx^{1}+\frac{1}{\lambda_{1}-\lambda_{3}}\frac{u_{2}}{u_{3}}dx^{2}\right)$
$\displaystyle\omega^{4}=\frac{1}{\lambda_{3}}dx^{3}-\frac{1}{\lambda_{3}}\left(\frac{\lambda_{2}}{\lambda_{2}-\lambda_{3}}\frac{v_{1}}{v_{3}}dx^{1}+\frac{\lambda_{1}}{\lambda_{1}-\lambda_{3}}\frac{u_{2}}{u_{3}}dx^{2}\right)$

and functions $u$ and $v$ satisfy the following system

(1)
$\displaystyle(\lambda_{3}-\lambda_{4})(\lambda_{1}-\lambda_{2})v_{1}u_{2}u_{34}+(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{4})(v_{1}u_{4}u_{23}-v_{3}u_{4}u_{12}+v_{3}u_{2}u_{14})$
$\displaystyle\quad-(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda_{3})(v_{1}u_{3}u_{24}-v_{4}u_{3}u_{12}+v_{4}u_{2}u_{13})=0,$
$\displaystyle(\lambda_{3}-\lambda_{4})(\lambda_{1}-\lambda_{2})v_{1}u_{2}v_{34}+(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{4})(v_{1}u_{4}v_{23}-v_{3}u_{4}v_{12}+v_{3}u_{2}v_{14})$
$\displaystyle\quad-(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda_{3})(v_{1}u_{3}v_{24}-v_{4}u_{3}v_{12}+v_{4}u_{2}v_{13})=0,$

where $\lambda_{i}\in\mathbb{R}$, $i=1,\ldots,4$, are fixed, arbitrary
pairwise distinct constants.

Moreover, the above system is integrable and admits an isospectral Lax pair
with the Lax operators

$\displaystyle
L_{0}=(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda)u_{2}v_{4}\partial_{x^{1}}-(\lambda_{1}-\lambda_{4})(\lambda_{2}-\lambda)u_{4}v_{1}\partial_{x^{2}}+(\lambda_{1}-\lambda_{2})(\lambda_{4}-\lambda)u_{2}v_{1}\partial_{x^{4}}$
$\displaystyle
L_{1}=(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda)u_{2}v_{3}\partial_{x^{1}}-(\lambda_{1}-\lambda_{3})(\lambda_{2}-\lambda)u_{3}v_{1}\partial_{x^{2}}+(\lambda_{1}-\lambda_{2})(\lambda_{3}-\lambda)u_{2}v_{1}\partial_{x^{3}}$

Alternative descriptions of hyper-para-Hermitian structures through integrable
systems have been studied in [7, 9, 13]. Notice that both equations in (1)
share the same symbol, which, in turn, recovers the conformal class by means
of the characteristic variety coinciding with the null cone of $[g]$. The
relation of system (1) to the aforementioned general heavenly equation is
explained in the following result.

######  Theorem 2.

If $w$ satisfies the general heavenly equation

(2)
$(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda_{3})w_{13}w_{24}-(\lambda_{3}-\lambda_{4})(\lambda_{1}-\lambda_{2})w_{12}w_{34}-(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{4})w_{14}w_{23}=0$

then

(3) $u=w_{1},\quad v=w_{2}$

satisfy (1).

Moreover, a hyper-para-Hermitian structure is hyper-para-Kähler if and only if
it can be locally put in the form from Theorem 1 with $(u,v)$ given by (3) for
some function $w$ satisfying (2).

The hyper-para-Kähler structures appearing in Theorem 2 can be introduced as
Ricci-flat hyper-para-Hermitian structures (see [9, 18]). According to the
Mason–Newman formalism [18], they are described by an integrable system with a
Lax pair whose Lax operators are vector fields that are divergence-free with
respect to some volume form. One realization of this formalism is provided by
the general heavenly equation (2) (we refer to [11, 16, 24] for details).
Other approaches include the Plebański equations and the Husain–Park equation,
among others, cf. [11, 23].

The rest of the paper is organized as follows.

Section 2 gives the proofs of Theorems 1 and 2. We apply the approach from
[13], where the hyper-para-Hermitian structures are studied through the lens
of Kronecker webs. Recall that the webs are special families of foliations,
originally introduced as reductions of certain bi-Hamiltonian systems [26].
Similar concepts have also been recently exploited in [23] and earlier in [8,
14] in the context of 3-dimensional Einstein–Weyl geometry.

Section 3 explores the properties of (1) as an integrable system. In
particular, we present an infinite hierarchy of nonlocal symmetries and a
recursion operator for (1).

##  2\. Hyper-para-hermitian structures

The aim of this section is to give proofs of Theorems 1 and 2. As mentioned in
the Introduction, we shall use a description of hyper-para-Hermitian
structures via Kronecker webs, which was recently exploited in [13] and [23].
In the 4-dimensional case studied in this paper, the webs are defined as
follows.

######  Definition.

A _Kronecker web_ on a 4-dimensional manifold is a 1-parameter family of
foliations $\\{\mathcal{F}_{\lambda}\\}_{\lambda\in\mathbb{R}}$ such that,
locally, there exist point-wise independent 1-forms
$\alpha^{1},\ldots,\alpha^{4}$ such that

$T\mathcal{F}_{\lambda}=\ker\left(\alpha^{1}+\lambda\alpha^{2},\alpha^{3}+\lambda\alpha^{4}\right).$

Notice that we can equivalently write

$T\mathcal{F}_{\lambda}=\operatorname{span}\left\\{X_{0}(\lambda),X_{1}(\lambda)\right\\}$

for $X_{0}(\lambda)=Y_{2}-\lambda Y_{1}$, $X_{1}(\lambda)=Y_{4}-\lambda
Y_{3}$, where $(Y_{i})_{i=1,\ldots,4}$ is a frame on $M$ dual to the co-frame
$(\alpha^{i})_{i=1,\ldots,4}$.

In [13], a one-to-one correspondence between hyper-para-Hermitian structures
and Kronecker webs on 4-dimensional manifolds is established. Here, the
endomorphisms $I$ and $K$ are defined such that the corresponding eigenspaces
$\mathcal{D}_{I}^{+}$, $\mathcal{D}_{I}^{-}$, $\mathcal{D}_{K}^{+}$ and
$\mathcal{D}_{K}^{-}$, associated with the eigenvalues $\pm 1$, respectively,
are given by $T\mathcal{F}_{\lambda_{i}}$ for certain fixed values of
$\lambda_{i}$, $i=1,2,3,4$ (see Corollaries 2.3, 2.5 and the discussion on
page 461 in [13]). Note that the complex structure $J$ is induced from $I$ and
$K$ by means of the composition $J=IK$. Additionally, from the viewpoint of
the anti-self-dual metrics, the leaves of foliations
$\\{\mathcal{F}_{\lambda}\\}_{\lambda\in\mathbb{R}}$ are $\alpha$-submanifolds
of $[g]$ and the collection of all leaves is usually referred to as the
twistor space, see [9]. Proof of Theorem 1. Let
$\lambda_{1},\ldots,\lambda_{4}\in\mathbb{R}$ be four distinct values of
$\lambda$. There are functions $f^{i},g^{i}$, $i=1,\ldots,4$, defined locally
in a neighborhood of a point $x\in M$, such that the corresponding foliations
$\mathcal{F}_{\lambda_{i}}$ are defined as $f^{i}=\mathrm{const}$,
$g^{i}=\mathrm{const}$, respectively. Without lost of generality, since the
foliations $\mathcal{F}_{\lambda}$ corresponding to different values of
$\lambda$ are transversal, one can define local coordinates in the
neighborhood of $x$ as $x^{i}=f^{i}$. Moreover, let $u=g^{1}$ and $v=g^{2}$.

We shall look for the $\lambda$-dependent vector fields $X_{0}(\lambda)$ and
$X_{1}(\lambda)$, spanning $T\mathcal{F}_{\lambda}$ for any
$\lambda\in\mathbb{R}$. They can be written in the most general form as

$X_{j}=(a^{1}_{j}+\lambda b^{1}_{j})\partial_{x^{1}}+(a^{2}_{j}+\lambda
b^{2}_{j})\partial_{x^{2}}+(a^{3}_{j}+\lambda
b^{3}_{j})\partial_{x^{3}}+(a^{4}_{j}+\lambda
b^{4}_{j})\partial_{x^{4}},\qquad j=0,1,$

for some functions $a^{i}_{j}$ and $b^{i}_{j}$, $i=1,\ldots,4$, $j=0,1$, in a
neighborhood of $x$. Now, since $dx^{i}$ annihilates both $X_{0}(\lambda_{i})$
and $X_{1}(\lambda_{i})$ we get $a^{i}_{j}=-\lambda_{i}b^{i}_{j}$ and
consequently

$X_{j}=(\lambda-\lambda_{1})b^{1}_{j}\partial_{x^{1}}+(\lambda-\lambda_{2})b^{2}_{j}\partial_{x^{2}}+(\lambda-\lambda_{3})b^{3}_{j}\partial_{x^{3}}+(\lambda-\lambda_{4})b^{4}_{j}\partial_{x^{4}},\qquad
j=0,1.$

Further on, $du$ annihilates $X_{0}(\lambda_{1})$ and $X_{1}(\lambda_{1})$ and
$dv$ annihilates $X_{0}(\lambda_{2})$ and $X_{1}(\lambda_{2})$. These
conditions give equations

(4)
$\displaystyle(\lambda_{1}-\lambda_{2})u_{2}b^{2}_{j}+(\lambda_{1}-\lambda_{3})u_{3}b^{3}_{j}+(\lambda_{1}-\lambda_{4})u_{4}b^{4}_{j}=0,$
$\displaystyle\qquad j=0,1$
$\displaystyle(\lambda_{2}-\lambda_{1})v_{1}b^{1}_{j}+(\lambda_{2}-\lambda_{3})v_{3}b^{3}_{j}+(\lambda_{2}-\lambda_{4})v_{4}b^{4}_{j}=0,$
$\displaystyle\qquad j=0,1.$

Moreover, $X_{1}$ and $X_{2}$ are defined up to a $GL(2)$-action. Hence, it
can be assumed that $b^{3}_{0}=0$, $b^{4}_{0}=1$, $b^{3}_{1}=1$ and
$b^{4}_{1}=0$. There are four $b^{i}_{j}$ left, with $i=1,2$ and $j=0,1$.
However, system (4) fixes them uniquely:

$b^{1}_{0}=\frac{(\lambda_{2}-\lambda_{4})}{(\lambda_{1}-\lambda_{2})}\frac{v_{4}}{v_{1}},\quad
b^{1}_{1}=\frac{(\lambda_{2}-\lambda_{3})}{(\lambda_{1}-\lambda_{2})}\frac{v_{3}}{v_{1}},\quad
b^{2}_{0}=-\frac{(\lambda_{1}-\lambda_{4})}{(\lambda_{1}-\lambda_{2})}\frac{u_{4}}{u_{2}},\quad
b^{2}_{1}=-\frac{(\lambda_{1}-\lambda_{3})}{(\lambda_{1}-\lambda_{2})}\frac{u_{3}}{u_{2}}.$

We get

(5) $\displaystyle
X_{0}=(\lambda-\lambda_{1})\frac{(\lambda_{2}-\lambda_{4})}{(\lambda_{1}-\lambda_{2})}\frac{v_{4}}{v_{1}}\partial_{x^{1}}-(\lambda-\lambda_{2})\frac{(\lambda_{1}-\lambda_{4})}{(\lambda_{1}-\lambda_{2})}\frac{u_{4}}{u_{2}}\partial_{x^{2}}+(\lambda-\lambda_{4})\partial_{x^{4}},$
$\displaystyle
X_{1}=(\lambda-\lambda_{1})\frac{(\lambda_{2}-\lambda_{3})}{(\lambda_{1}-\lambda_{2})}\frac{v_{3}}{v_{1}}\partial_{x^{1}}-(\lambda-\lambda_{2})\frac{(\lambda_{1}-\lambda_{3})}{(\lambda_{1}-\lambda_{2})}\frac{u_{3}}{u_{2}}\partial_{x^{2}}+(\lambda-\lambda_{3})\partial_{x^{3}}.$

Note that setting $L_{i}=(\lambda_{1}-\lambda_{2})u_{2}v_{1}X_{i}$ recovers
the Lax operators from Theorem 1. By assumption, $X_{0}(\lambda)$ and
$X_{1}(\lambda)$ span an integrable distribution for any $\lambda$. Hence,
$X_{0}(\lambda)$ and $X_{1}(\lambda)$ necessarily commute. One readily checks
that vanishing of the commutator of $X_{i}(\lambda)$ is equivalent to system
(1). Likewise, it is easily seen that vanishing of the commutator of $L_{i}$
modulo a certain linear combination of $L_{i}$ is equivalent to system (1).
Thus an overdetermined linear system for $\psi$, $L_{i}\psi=0$, $i=0,1$, whose
compatibility condition is nothing but (1), provides a Lax pair for (1).
$\hfill\Box$

######  Lemma 1.

Suppose that $(u,v)$ is a solution to (1). Then, for any two nowhere vanishing
functions $a,b\colon\mathbb{R}^{2}\to\mathbb{R}$, the pair
$(\tilde{u},\tilde{v})$ defined as

$\displaystyle\tilde{u}(x^{1},x^{2},x^{3},x^{4})=a(x^{1},u(x^{1},x^{2},x^{3},x^{4}))u(x^{1},x^{2},x^{3},x^{4}),$
$\displaystyle\tilde{v}(x^{1},x^{2},x^{3},x^{4})=b(x^{2},v(x^{1},x^{2},x^{3},x^{4}))v(x^{1},x^{2},x^{3},x^{4}),$

is a solution to (1) descending to the same hyper-para-Hermitian metric.

###### Proof.

Observe that, in view of the proof of Theorem 1, the functions $\tilde{u}$ and
$\tilde{v}$ belong to the rings of functions constant on leaves of foliations
$\mathcal{F}_{\lambda_{1}}$ and $\mathcal{F}_{\lambda_{2}}$, respectively.
Hence $\tilde{u}$ and $\tilde{v}$ can replace $u$ and $v$ in the derivation of
the Lax pair in the proof of Theorem 1, and therefore in system (1) itself,
and the result follows. ∎

Proof of Theorem 2. The first part of the theorem can be verified by direct
computations. In order to prove the second part, we shall work using the
coordinates from Theorem 1. The Lax distribution
$\operatorname{span}\\{L_{0}(\lambda),L_{1}(\lambda)\\}$ is annihilated by two
$\lambda$-dependent 1-forms that can be taken in the following form

$\displaystyle\alpha_{0}(\lambda)=$
$\displaystyle(\lambda_{1}-\lambda_{2})(\lambda_{4}-\lambda)\left(\frac{u_{3}}{u_{2}}dx^{1}+\frac{v_{3}}{v_{1}}dx^{2}+\frac{u_{3}v_{3}}{u_{2}v_{1}}dx^{3}\right)$
$\displaystyle+\left((\lambda_{1}-\lambda_{4})(\lambda_{2}-\lambda)\frac{u_{4}v_{3}}{u_{2}v_{1}}-(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda)\frac{u_{3}v_{4}}{u_{2}v_{1}}\right)dx^{4}$
$\displaystyle\alpha_{1}(\lambda)=$
$\displaystyle(\lambda_{1}-\lambda_{2})(\lambda_{3}-\lambda)\left(\frac{u_{4}}{u_{2}}dx^{1}+\frac{v_{4}}{v_{1}}dx^{2}+\frac{u_{4}v_{4}}{u_{2}v_{1}}dx^{3}\right)$
$\displaystyle+\left((\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda)\frac{u_{4}v_{3}}{u_{2}v_{1}}-(\lambda_{1}-\lambda_{3})(\lambda_{2}-\lambda)\frac{u_{3}v_{4}}{u_{2}v_{1}}\right)dx^{4}$

Let $\beta=\alpha_{0}\wedge\alpha_{1}$. According to [23, Corollary 2.4] a
hyper-para-Hermitian structure is divergence-free (which in the terminology of
the present paper means that the structure is hyper-para-Kähler) if and only
if there is a nowhere vanishing function $c\colon M\to\mathbb{R}$ such that
$c\beta$ is closed. Since $L_{0}$ and $L_{1}$ commute modulo a linear
combination thereof with nonconstant coefficients, see the proof of Theorem 1,
we have

$d\alpha_{i}=0\mod\alpha_{0},\alpha_{1},$

and consequently

$d\beta=\varphi\wedge\beta$

where $\varphi$ is explicitly written as

$\displaystyle\varphi(\lambda)=\frac{1}{(\lambda_{1}-\lambda_{2})(\lambda-\lambda_{3})u_{2}v_{1}}$
$\displaystyle\left((\lambda_{2}-\lambda_{3})(\lambda-\lambda_{1})(u_{12}(v_{3}-u_{23}v_{1})\right.$
$\displaystyle\left.+(\lambda_{1}-\lambda_{3})(\lambda-\lambda_{2})(u_{2}v_{13}-u_{3}v_{12})\right)dx^{3}$
$\displaystyle+\frac{1}{(\lambda_{1}-\lambda_{2})(\lambda-\lambda_{4})u_{2}v_{1}}$
$\displaystyle\left((\lambda_{2}-\lambda_{4})(\lambda-\lambda_{1})(u_{12}v_{4}-u_{24}v_{1})\right.$
$\displaystyle\left.+(\lambda_{1}-\lambda_{4})(\lambda-\lambda_{2})(u_{2}v_{14}-u_{4}v_{12})\right)dx^{4}$

Notice that $\varphi$ is specified only up to transformations
$\varphi\mapsto\varphi_{A^{0},A^{1}}=\varphi+A^{0}\alpha_{0}+A^{1}\alpha_{1}$
for some functions $A^{0}$ and $A^{1}$. Since the result is local, we can
apply the Poincaré lemma for differential forms, and get that function $c$
exists if and only if $\varphi_{A^{0},A^{1}}$ is exact for an appropriate
choice of $A^{0}$ and $A^{1}$. Moreover $\varphi_{A^{0},A^{1}}$ has to be a
function on $M$ (i.e. it cannot depend on $\lambda$). We now consider the
equation

$\varphi_{A^{0},A^{1}}=\varphi+A^{0}\alpha_{0}+A^{1}\alpha_{1}=df$

where $A^{0},A^{1}$ and $f$ are unknown. Examining the coefficients at
$dx^{1}$ and $dx^{2}$, we obtain an algebraic system for $A^{i}$’s. This
system can be effortlessly solved, resulting in

$A^{0}=\frac{1}{(\lambda_{1}-\lambda_{2})(\lambda-\lambda_{4})}\left(\frac{u_{4}v_{1}f_{2}-u_{2}v_{4}f_{1}}{u_{3}v_{4}-u_{4}v_{3}}\right),\
\
A^{1}=\frac{1}{(\lambda_{1}-\lambda_{2})(\lambda-\lambda_{3})}\left(\frac{u_{2}v_{3}f_{1}-u_{3}v_{1}f_{2}}{u_{3}v_{4}-u_{4}v_{3}}\right).$

Then, the coefficients at $dx^{3}$ and $dx^{4}$ reduce to differential
equations for $f$, which, in turn, can be written in the form

(6) $L_{0}(f)=(\lambda-\lambda_{4})\varphi_{4},\qquad
L_{1}(f)=(\lambda-\lambda_{3})\varphi_{3},$

where $\varphi_{3}$ and $\varphi_{4}$ stand for the coefficients defined as
$\varphi=\varphi_{3}dx^{3}+\varphi_{4}dx^{4}$. The system always satisfies the
compatibility condition, which can be verified directly. Hence, it yields a
solution $f$. However, in general, $f$ depends on $\lambda$. Imposing the
condition that $\partial_{\lambda}f=0$ and substituting $\lambda=\lambda_{1}$
and $\lambda=\lambda_{2}$ in (6) gives the following system

$\displaystyle u_{2}v_{1}f_{3}-u_{3}v_{1}f_{2}=u_{2}v_{13}-u_{3}v_{12},$
$\displaystyle u_{2}v_{1}f_{3}-u_{2}v_{3}f_{1}=u_{23}v_{1}-u_{12}v_{3},$
$\displaystyle u_{2}v_{1}f_{4}-u_{4}v_{1}f_{2}=u_{2}v_{14}-u_{4}v_{12},$
$\displaystyle u_{2}v_{1}f_{4}-u_{2}v_{4}f_{1}=u_{24}v_{1}-u_{12}v_{4}.$

It turns out that this system can be rewritten as

$\displaystyle\left(\partial_{x^{3}}-\frac{u_{3}}{u_{2}}\partial_{x^{2}}\right)(f-\ln(v_{1}))=0,\qquad\left(\partial_{x^{3}}-\frac{v_{3}}{v_{1}}\partial_{x^{1}}\right)(f-\ln(u_{2}))=0,$
$\displaystyle\left(\partial_{x^{4}}-\frac{u_{4}}{u_{2}}\partial_{x^{2}}\right)(f-\ln(v_{1}))=0,\qquad\left(\partial_{x^{4}}-\frac{v_{4}}{v_{1}}\partial_{x^{1}}\right)(f-\ln(u_{2}))=0.$

Denoting $g=f-\ln(v_{1})$ and $h=f-\ln(u_{2})$ and using the fact that the
vector fields $\partial_{x^{3}}-\frac{u_{3}}{u_{2}}\partial_{x^{2}}$ and
$\partial_{x^{4}}-\frac{u_{4}}{u_{2}}\partial_{x^{2}}$ commute, and so do
$\partial_{x^{3}}-\frac{v_{3}}{v_{1}}\partial_{x^{1}}$ and
$\partial_{x^{4}}-\frac{v_{4}}{v_{1}}\partial_{x^{1}}$, we get that $g$ and
$h$ can be written as $g=g(x^{1},u)$ and $h=h(x^{2},v)$, because $x^{1}$ and
$u$ are constant along the first pair of vector fields, and $x^{2}$ and $v$
are constant along the second pair of vector fields. Further, expressing $f$
in terms of $g$ and $h$ and comparing the expressions gives

$\ln(u_{2})-\ln(v_{1})=h-g.$

By Lemma 1 there is a different solution to (1) such that
$\ln(\tilde{u}_{2})-\ln(\tilde{v}_{1})=0$. Indeed, one can adjust functions
$a$ and $b$ from Lemma 1, such that $h$ and $g$ cancel out. Consequently,

$\ln\left(\frac{\tilde{u}_{2}}{\tilde{v}_{1}}\right)=0.$

Hence, we have $\tilde{u}_{2}=\tilde{v}_{1}$ and it implies that the 1-form

$\tilde{u}dx^{1}+\tilde{v}dx^{2}$

is closed with respect to coordinates $(x^{1},x^{2})$, meaning that locally
$\tilde{u}=w_{1}$ and $\tilde{v}=w_{2}$, for some function $w$, i.e. (3)
holds. Furthermore, the Lax pair in Theorem 1, under assumption (3), reduces
to the Lax pair for (2), as can be found, for instance, in [11, formula
(4.11)]. It follows that $w$ is a solution to (2).$\hfill\Box$

######  Remark.

The second assertion of Theorem 2 can be expressed in terms of a curvature.
Namely, in order for a hyper-para-Hermitian structure to be hyper-para-Kähler
the Obata connection associated with the former has to be Ricci-flat. Recall
that the Obata connection is the unique connection $\nabla$ such that $\nabla
I=\nabla J=\nabla K=0$. It was initially introduced within the framework of
hyper-Hermitian structures [21], but its definition extends seamlessly to the
neutral signature, see [1]. In this context it coincides with the Chern
connection of webs [20], defined by the following expression (see also [13,
formula (2.1)])

$\nabla_{X}Y=\pi_{H}(j[\pi_{H}(X),j\pi_{H}(Y)]+[\pi_{V}(X),\pi_{H}(Y)])+\pi_{V}(j[\pi_{V}(X),j\pi_{V}(Y)]+[\pi_{H}(X),\pi_{V}(Y)])$

where $[.,.]$ is the Lie bracket of vector fields, and $\pi_{V}$ and $\pi_{H}$
are projections to the factors of the decomposition $TM=V\oplus H$, where

$V=\operatorname{span}\\{L_{0}(\lambda_{1}),L_{1}(\lambda_{1})\\},\quad
H=\operatorname{span}\\{L_{0}(\lambda_{2}),L_{1}(\lambda_{2})\\},$

and $j\colon TM\to TM$ is a mapping satisfying $j^{2}=1$ and uniquely
determined by properties

$j\colon V\to H,\quad j\colon H\to V,\quad j(\pi_{H}(X))-\pi_{H}(X)\in T,\quad
j(\pi_{V}(X))-\pi_{V}(X)\in T,$

for any $X\in TM$, where
$T=\operatorname{span}\\{L_{0}(\lambda_{3}),L_{1}(\lambda_{3})\\}$. In the
present context

$V=\operatorname{span}\\{u_{4}\partial_{2}-u_{2}\partial_{4},u_{3}\partial_{2}-u_{2}\partial_{3}\\},\quad
H=\operatorname{span}\\{v_{4}\partial_{1}-v_{1}\partial_{4},v_{3}\partial_{1}-v_{1}\partial_{3}\\}$

and

$\displaystyle
j(u_{3}\partial_{2}-u_{2}\partial_{3})=\frac{u_{2}}{v_{1}}(v_{3}\partial_{1}-v_{1}\partial_{3}),\quad
j(u_{4}\partial_{2}-u_{2}\partial_{4})=C\frac{u_{2}}{v_{1}}(v_{4}\partial_{1}-v_{1}\partial_{4}),$
$\displaystyle
j(v_{3}\partial_{1}-v_{1}\partial_{3})=\frac{v_{1}}{u_{2}}(u_{3}\partial_{2}-u_{2}\partial_{3}),\quad
j(v_{4}\partial_{1}-v_{1}\partial_{4})=C^{-1}\frac{v_{1}}{u_{2}}(u_{4}\partial_{2}-u_{2}\partial_{4}).$

where
$C=\frac{(\lambda_{3}-\lambda_{1})(\lambda_{4}-\lambda_{2})}{(\lambda_{4}-\lambda_{1})(\lambda_{3}-\lambda_{2})}$.

##  3\. Symmetries and Recursion Operator

Now turn to the study of symmetries for (1). To this end we first note that
the operators $X_{i}$ given by (5) from the proof of Theorem 1 are linear in
$\lambda$ and thus can be written as

(7) $X_{i}=X_{i}^{(0)}-\lambda X_{i}^{(1)},\quad i=0,1$

with obvious expressions for $X_{i}^{(j)}$.

Consider now the following ‘adjoint Lax pair’ for (1):

(8) $[Q,X_{i}]=0,\quad i=0,1,$

where $[\cdot,\cdot]$ is again the usual Lie bracket of vector fields and

(9) $Q=\sum\limits_{j=1}^{2}\xi^{j}\partial_{x^{j}}$

(in this connection recall that $[X_{0},X_{1}]=0$ modulo (1)).

The formal expansions

(10) $\xi^{j}=\sum\limits_{r=0}^{\infty}\xi_{r}^{j}\lambda^{r},\quad j=1,2$

give rise to an infinite hierarchy of nonlocal variables $\xi_{r}^{j}$
associated with (1) as follows.

Let $Q_{r}=\sum\limits_{j=1}^{2}\xi_{r}^{j}\partial_{x^{j}}$. Then
substituting (9) and (10) into (8) yields

(11) $[Q_{0},X_{i}^{(0)}]=0,\quad i=0,1$

and a set of recursion relations

(12) $[Q_{r},X_{i}^{(0)}]=[Q_{r-1},X_{i}^{(1)}],\quad i=0,1,\quad r=1,2,\dots$

relating $Q_{r}$ and $Q_{r-1}$.

Equations (11) and (12) can be solved with respect to
$\partial\xi_{r}^{j}/\partial x^{m}$, $m=3,4$, to yield relations of general
form

$\partial\xi_{r}^{j}/\partial x^{m}=A_{mr}^{j},\quad m=3,4,\quad j=1,2,\quad
r=0,1,2,\dots$

which recursively define $\xi_{r}^{j}$ starting from $r=0$; here $A_{mr}^{j}$
are certain functions, a bit too cumbersome to spell out here in full, of
$\xi_{s}^{j}$, $s=0,\dots,r$ and $x$\- and $y$-derivatives of those, and of a
number of first- and second-order derivatives of $u$ and $v$.

In other words, we have here an infinite-dimensional differential covering
over (1) with the nonlocal variables $\xi_{r}^{j}$, $j=1,2$, $r=0,1,2,\dots$;
the said covering is defined via (11) and (12). For generalities on nonlocal
variables, differential coverings and nonlocal symmetries the reader is
referred to [12, 17, 25] and references therein.111Note, however, that the
authors of [12] refer to what we call below nonlocal symmetries as to the
shadows of nonlocal symmetries.

With this in mind, we readily arrive at the following result

######  Proposition 1.

The flows

(13) $u_{\tau_{r}}=\xi_{r}^{2}u_{2},\quad v_{\tau_{r}}=-\xi_{r}^{1}v_{1},\quad
r=0,1,2,\dots,$

with $\xi_{r}^{k}$ defined above, are compatible with (1), and thus define an
infinite hierarchy of nonlocal symmetries for (1) with the characteristics
$(\xi_{r}^{2}u_{2},-\xi_{r}^{1}v_{1})$, $r=0,1,2,\dots$.

Note that such infinite hierarchies of symmetries are a common occurrence for
integrable partial differential systems, and for systems in more than two
independent variables the symmetries in question are usually nonlocal, cf.
e.g. [12, 17, 19, 22, 25] and references therein.

Sketch of proof. It is straightforward to verify that $U=\xi^{2}u_{2}$ and
$V=-\xi^{1}v_{1}$ satisfy the linearized version of (1), that is,

(14)
$\displaystyle(\lambda_{3}-\lambda_{4})(\lambda_{1}-\lambda_{2})(V_{1}u_{2}u_{34}+v_{1}U_{2}u_{34}+v_{1}u_{2}U_{34})$
$\displaystyle\quad+(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{4})(V_{1}u_{4}u_{23}+v_{1}U_{4}u_{23}+v_{1}u_{4}U_{23}$
$\displaystyle\quad-
V_{3}u_{4}u_{12}-v_{3}U_{4}u_{12}-v_{3}u_{4}U_{12}+V_{3}u_{2}u_{14}+v_{3}U_{2}u_{14}+v_{3}u_{2}U_{14})$
$\displaystyle\quad-(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda_{3})(V_{1}u_{3}u_{24}+v_{1}U_{3}u_{24}+v_{1}u_{3}U_{24}$
$\displaystyle\quad-
V_{4}u_{3}u_{12}-v_{4}U_{3}u_{12}-v_{4}u_{3}U_{12}+V_{4}u_{2}u_{13}+v_{4}U_{2}u_{13}+v_{4}u_{2}U_{13})=0,$
$\displaystyle(\lambda_{3}-\lambda_{4})(\lambda_{1}-\lambda_{2})(V_{1}u_{2}v_{34}+v_{1}U_{2}v_{34}+v_{1}u_{2}V_{34})$
$\displaystyle\quad+(\lambda_{2}-\lambda_{3})(\lambda_{1}-\lambda_{4})(V_{1}u_{4}v_{23}+v_{1}U_{4}v_{23}+v_{1}u_{4}V_{23}$
$\displaystyle\quad-
V_{3}u_{4}v_{12}-v_{3}U_{4}v_{12}-v_{3}u_{4}V_{12}+V_{3}u_{2}v_{14}+v_{3}U_{2}v_{14}+v_{3}u_{2}V_{14})$
$\displaystyle\quad-(\lambda_{2}-\lambda_{4})(\lambda_{1}-\lambda_{3})(V_{1}u_{3}v_{24}+v_{1}U_{3}v_{24}+v_{1}u_{3}V_{24}$
$\displaystyle\quad-
V_{4}u_{3}v_{12}-v_{4}U_{3}u_{12}-v_{4}u_{3}V_{12}+V_{4}u_{2}v_{13}+v_{4}U_{2}v_{13}+v_{4}u_{2}V_{13})=0,$

modulo (1), (8) and differential consequences thereof, i.e., (1) admits a
nonlocal symmetry with the characteristic $(\xi^{2}u_{2},-\xi^{1}v_{1})$ .
Upon taking into account linearity of (14) in $U$ and $V$ and substituting the
expansions (10) into (14) it is readily seen that $U=\xi_{s}^{2}u_{2}$ and
$V=-\xi_{s}^{1}v_{1}$ for $s=0,1,2,\dots$ also satisfy (14) modulo (1), (11),
and (12) for $r=1,\dots,s$, and differential consequences thereof, so (1)
indeed admits an infinite hierarchy of nonlocal symmetries with the
characteristics $(\xi_{s}^{2}u_{2},-\xi_{s}^{1}v_{1})$, $s=0,1,2,\dots$.
Compatibility of (13) with (1) immediately follows from this, cf. e.g. [12,
22] and references therein. $\hfill\Box$

In closing note that, using the technique from [25], we have arrived at the
following result, which is readily verified by straightforward albeit fairly
tedious computation (for background on total derivatives and recursion
operators seen essentially as Bäcklund auto-transformations for the linearized
version of the system under study see e.g. [12, 17, 19, 25]):

######  Proposition 2.

Let the flow

(15) $u_{\tau}=U,\quad v_{\tau}=V,$

where $U$ and $V$ are assumed to be functions of independent variables, jet
variables intrinsic to (1), and some nonlocal variables over (1), be
compatible with (1) and thus define a nonlocal symmetry with the
characteristic $(U,V)$ for (1).

Consider nonlocal variables $\zeta^{i}$ defined by the relations

(16) $[X_{i}^{(1)},S]=\left([X_{i}^{(0)},R]\right)^{TD},\quad i=0,1$

where $X_{i}^{(j)}$ are defined via (7), and

(17) $R=-(V/v_{1})\partial_{x^{1}}+(U/u_{2})\partial_{x^{2}},\quad
S=\sum\limits_{j=1}^{2}\zeta^{j}\partial_{x^{j}}$

and the superscript $TD$ means that expressions involving the derivatives of
$U$ and $V$ like $U_{x^{i}}$ and $V_{x^{j}}$ should be replaced by
$D_{x^{i}}U$ and $D_{x^{j}}V$, where $D_{x^{i}}$ are total derivatives.

Then the new flow

(18) $u_{\sigma}=\tilde{U},\quad v_{\sigma}=\tilde{V},$

where

(19) $\tilde{U}=u_{2}\zeta^{2},\quad\tilde{V}=-v_{1}\zeta^{1},$

is also compatible with (1) and thus defines another nonlocal symmetry with
the characteristic $(\tilde{U},\tilde{V})$ for (1).

In other words, the relations (16), (17), and (19) define a recursion operator
for (1) associating to any nonlocal symmetry with the characteristic $(U,V)$
for (1) a new (again in general nonlocal) symmetry with the characteristic
$(\tilde{U},\tilde{V})$ for (1).

The above is a natural generalization of the recursion operator for general
heavenly equation found in [25].

Notice that one now can construct (in general nonlocal) symmetries for system
(1) beyond those from Proposition 1 by applying the recursion operator from
Proposition 2 e.g. to ‘simpler’ symmetries like, say, those with the
characteristics $(u_{i},v_{i})$, $i=1,\dots,4$, that result from the
translation invariance of (1).

### Acknowledgments

The research of WK, and the visit of AS to Warsaw, were supported in part by
the grant 2019/34/E/ST1/00188 from the National Science Centre (NCN), Poland.
The research of AS was also supported in part through institutional funding
for the development of research organizations (RVO) for IČ 47813059. The
authors gratefully acknowledge the support from the above sources.

WK and AS thank respectively Silesian University in Opava and Institute of
Mathematics of Polish Academy of Sciences in Warsaw for warm hospitality
extended to them in the course of their visits to the respective institutions.

Some of the computations in the paper were performed employing the package
Jets [2] for Maple® whose use is hereby gratefully acknowledged.

## References

  * [1] D. Alekseevsky, C. Medori, A. Tomassini, Homogeneous para-Kähler Einstein manifolds, Russ. Math. Surv. 64, No. 1 (2009), 1–43.
  * [2] H. Baran, M. Marvan, Jets. A software for differential calculus on jet spaces and diffieties. Available online at http://jets.math.slu.cz/
  * [3] A. Dancer, H. Jorgensen, A. Swann, Metric geometries over the split quaternions, Rend. Sem. Mat. Univ. Politec. Torino 63 (2005), 119–139.
  * [4] J. Davidov, G. Grantcharov, O. Mushkarov, M. Yotov, Para-hyperhermitian surfaces, Bull. Math. Soc. Sci. Math. Roumanie (N.S.) 52(100) (2009), no. 3, 281–289.
  * [5] J. Davidov, G. Grantcharov, O. Mushkarov, Complex surfaces and null conformal Killing vector fields, J. Geom. Anal. 33 (2023), no. 7, Paper No. 224, 18 pp.
  * [6] B. Doubrov, E.V. Ferapontov, On the integrability of symplectic Monge–Ampère equations, J. Geom. Phys. 60 (2010), 1604–1616
  * [7] M. Dunajski, The twisted photon associated to hyper-hermitian four-manifolds, J. Geom. Phys. 30 (1999), 266–281
  * [8] M. Dunajski and W. Kryński, Einstein–Weyl geometry, dispersionless Hirota equation and Veronese webs, Math. Proc. Camb. Phil. Soc. 157 (2014), no. 1, 139–150.
  * [9] M. Dunajski, S. West, Anti-self-dual conformal structures in neutral signature, in Recent developments in pseudo-Riemannian geometry, 113–148. ESI Lect. Math. Phys. European Mathematical Society (EMS), Zürich, 2008
  * [10] S. Ivanov, S. Zamkovoy, ParaHermitian and paraquaternionic manifolds, Diff. Geom. App. 23(2) (2005) 205–234
  * [11] B. Konopelchenko, W. Schief, A. Szereszewski, Self-dual Einstein spaces and the general heavenly equation. Eigenfunctions as coordinates, Class. Quant. Grav. 38 (2021), no. 4, Paper No. 045007, 30 pp.
  * [12] J. Krasil’shchik, A.M. Verbovetsky, R. Vitolo, The Symbolic Computation of Integrability Structures for Partial Differential Equations, Springer, Cham, 2017.
  * [13] W. Kryński, Webs and the Plebański equations, Math. Proc. Camb. Phil. Soc. 161(3) (2016) 455–468.
  * [14] W. Kryński, On deformations of the dispersionless Hirota equation, J. Geom. Phys. 127 (2018) 46–54.
  * [15] W. Kryński, Deformations of dispersionless Lax systems, Class. Quant. Grav. 40(23) (2023), Paper No. 235013, 12 pp.
  * [16] A.A. Malykh, M.B. Sheftel, General heavenly equation governs anti-self-dual gravity, J. Phys. A: Math. Theor. 44 (2011), Paper No. 155201, 11 pp.
  * [17] M. Marvan, A. Sergyeyev, Recursion operators for dispersionless integrable systems in any dimension, Inverse Problems 28 (2012), no. 2, Paper No. 025011, 12 pp.
  * [18] L. Mason, E.T. Newman, A connection between the Einstein and Yang-Mills equations, Comm. Math. Phys. 121 (1989), 659–668.
  * [19] O.I. Morozov, A. Sergyeyev, The four-dimensional Martínez Alonso–Shabat equation: reductions and nonlocal symmetries, J. Geom. Phys. 85 (2014), 40–45.
  * [20] P. Nagy, Webs and curvature, In Web Theory and Related Topics. (Toulouse, December, 1996), World Scientific, River Edge, 2001, 48–91.
  * [21] M. Obata, Affine connections on manifolds with almost complex, quaternion or Hermitian structure, Jap. J. Math. 26 (1956), 43–77.
  * [22] P.J. Olver, Applications of Lie groups to differential equations, 2nd ed. Springer, New York, 1993.
  * [23] A. Panasyuk, A. Szereszewski, Webs, Nijenhuis operators, and heavenly PDEs, Class. Quantum Grav. 40 (2023), Paper No. 235003, 33 pp.
  * [24] W. Schief, Self-dual Einstein spaces via a permutability theorem for the Tzitzeica equation, Phys. Lett. A 223 (1996), 55–62.
  * [25] A. Sergyeyev, A simple construction of recursion operators for multidimensional dispersionless integrable systems, J. Math. Anal. Appl. 454 (2017), 468–480.
  * [26] I. Zakharevich, Kronecker webs, bihamiltonian structures and the method of argument translation, Transform. Groups. 6, No. 3 (2001), 267–300.