# Theory of electron and ion holes as vortices in the phase-space of
collision-less plasmas

Allen Lobo<EMAIL_ADDRESS>Sikkim Manipal Institute of Technology,
Sikkim Manipal University, Sikkim - 737136, India.  Vinod Kumar Sayal
<EMAIL_ADDRESS>Sikkim Manipal Institute of Technology, Sikkim Manipal
University, Sikkim - 737136, India.

###### Abstract

This article studies the vortical nature and structure of phase-space holes –
nonlinear B.G.K. trapping modes found in the phase-space collision-free
plasmas. A fluid-like outlook of the particles’ phase-space is introduced,
which makes it convenient to analytically identify electron and ion holes as
vortices – similar to that of ordinary two-dimensional fluids. A fluid
velocity and vorticity field is defined for the phase-space of the electrons
and ions, Euler equations describing the flow of the phase-space fluid
representing the particle system are then developed. Using these equations,
electron holes and ion holes are analytically identified as vortices in the
phase-space of the plasma. A relation between Schamel’s trapping parameter
($\beta$), hole speed ($M$), hole phase-space depth ($-\Gamma$) and hole
potential amplitude ($\chi_{0}$) is derived. The approach introduces a new
technique to study the phase-space holes of collision-less plasmas, allowing
fluid-vortex-like treatment to these kinetic structures. Phase-space
distribution functions for electron hole regions can then be analytically
derived from this model, reproducing the schamel-df equations and thus acting
as a precursor to the pseudo-potential approach, avoiding the need to assume a
solution to the phase-space density.

##  I Introduction

The problem of electron and ion holes of collision-less plasmas is, perhaps,
one of the most complex topics that the theory of collision-free plasma has to
deal with, comprising many stochastic elements. It refers to phenomena for
which there is no linear approach. Electron holes (and their positive counter-
parts) are well-identifiable species in their respective phase-space with the
core particle density much lower than the boundary, resembling a Gaussian-like
potential in their respective configuration space, as seen in laboratory and
simulation experiments Morse and Nielson, (1969); Turikov, (1978); Saeki et
al., (1979); Lynov et al., (1979). They are the result of nonlinear trapping
mechanisms and are identified as B.G.K. modes Bernstein et al., (1957) with
velocities ($M$) comparable to thermal speeds Iizuka and Tanaca, (1987);
Hutchinson and Zhou, (2016); Saeki et al., (1979); Lynov et al., (1979). Such
trapping may be the result of one of the numerous mechanisms Schamel, 2023a .
The study of electron (and ion) holes is conducted using differential or
integral approaches, the former employing an assumed distribution function
Schamel, (1979) for the trapped and free particles in the hole region and the
latter employing an assumed potential function resembling a Gaussian or an
inverted bell-shape Turikov, (1984); Chen et al., (2005); Aravindakshan et
al., (2021). Schamel’s pseudopotential approach utilizes the schamel-df
equations Schamel, (1979, 1986) (assumed electron phase-space density
functions for free and trapped particles) and produces a nonlinear dispersion
relation expressing hole amplitude $\chi_{0}$ as a function of a particle-
trapping parameter $\beta$, wave number $k$ and hole speed $M$. However, the
value of the $\beta$ parameter can not be exactly defined – except that it
must be a negative value – allowing arbitrariness in the dispersion relation.
Interestingly, the electron and ion holes have been ubiquitously referred to
as phase-space vortices Morse and Nielson, (1969); Saeki et al., (1979);
Schamel, (1979, 1982); Turikov, (1984); Schamel, (1998); Guio et al., (2003);
Trivedi and Ganesh, (2017); Hutchinson, (2017). However, no analytical proof
of their vortical nature has been presented till now. On the other hand,
vortices and vortical behaviour of ordinary two-dimensional fluids present
themselves with an abundant history of mathematical analyses. The
correspondence between an ordinary two-dimensional fluid and a one-dimensional
(1D-1V) phase-space is well-known, with the Vlasov equation in the collision-
less plasma case being identical to the Liouville equation for ideal in-
compressible two-dimensional fluids Schamel, (1986, 2012). It should be noted
that for such an analysis, the phase-space of such systems should be studied
in analogy to ordinary fluid surfaces – electron and ion holes are vortices in
the phase-space of plasma particles, and not in the configuration space (where
they present themselves as nonlinear B.G.K. mono-polar potential waves).
Hence, a fluid-like modelling of the particle phase-space must be introduced
to analytically identify and treat electron holes and ion holes as phase-space
vortices.
In this article, the electron phase-space of collision-less plasmas is
portrayed in a fluid-like outlook. Analogous to conventional two-dimensional
fluid dynamics, the 1D-1V phase-space is transformed into a two-dimensional
plane, wherein a fluid-like dynamical model of the particle phase-space is
described. Euler equations and diffusion equations are framed to demonstrate
the fluid-like behaviour of the particle phase-space using the velocity and
vorticity fields of the phase-space fluid. On the bases of these equations and
pre-established techniques of vortex identification in fluids Tian et al.,
(2018); Lugt, (1979); Helmholtz, (1867), electron holes (and by analogy, ion
holes) are analytically identified as vortices in the phase-space of the
system. The same equations are then used to derive a steady-state structure of
the phase-space density in the hole region, reproducing the schamel-df
equations Schamel, (1979). The fluid description of the phase-space of the
plasma particles (and not the plasma itself) provides a new outlook to
understand electron and ion holes observed in the phase-space of collision-
less plasma. The equations developed using this approach analytically
identifies such holes as vortices in phase-space and provides techniques
analogous to the study of a fluid vortex to explore their structure and
behaviour. Using this approach, the phase-space density function of both free
and trapped particles can be derived, without the need of an initial
assumption.

##  II Theory of phase-space vortices

The kinetic study of many-particle systems thoroughly describes the
(stochastic) behaviour of the system – which comprises particles with a range
of energy, under influence of some interacting potential. This is represented
by the Hamiltonian of the system itself –

$T(v_{x})=H(f_{x,v_{x}})-V(x).$ (1)

Here $T(v_{x})$, $H$, $f(x,v_{x})$ and $V(x)$ represent the Kinetic Energy
function, Hamiltonian, phase-space probability density and Potential Energy
function. This analysis, however, when conducted for a dynamically evolving
system, describes a flow-like evolution of the particle phase-space –

$p_{x}\dot{v}_{x}=\int\frac{df}{dt}dE-\dot{x}\frac{\partial V(x)}{\partial
x}.$ (2)

Here, $p_{x}=\frac{\partial T}{\partial v_{x}}$ is the momentum of a particle,
$E$ is its energy and $\dot{q}$ represents an ordinary time derivative of $q$.
The terms $\dot{v}_{x}$ and $\dot{x}$ can be observed as the flow-rate of a
particle in velocity and position axes as it slides along a phase-space
trajectory, in which case the above equation (2) describes flow-like evolution
of a system in its phase-space, providing the dynamics of the many-particle
phase-space with a fluid-like outlook. The position-velocity phase-space of a
system can be described as a two-dimensional surface by multiplying a system-
specific time-constant to the velocity axis. Such a time-factor can be a
characteristic response time of the particles, characteristic to the system
itself, and we represent it as $\tau$ for our analysis. The phase-space then
describes a co-ordinate set $x-\tau v_{x}$ to locate an infinitesimal volume
$\tau dxdv_{x}$. The particle density of this phase-space is the number of
particles present in this volume –

$\frac{\partial^{2}N}{\tau\partial x\partial
v_{x}}=\frac{1}{\tau}f(x,v_{x})=\eta.$ (3)

Here, $N$ is the total number of particles in the (closed) system. The phase-
space gradient operator can similarly be defined as the derivative along each
axis –

$\nabla=\frac{\partial}{\partial x}\hat{x}+\frac{\partial}{\tau\partial
v_{x}}\hat{v}_{x}.$ (4)

One can now describe the associated position $\bm{r}$ and velocity field
$\mathcal{V}$ of the phase-space fluid111The stream function for describing
the fluid-like self-consistency of the phase-space-fluid model can be
described using a modified particle Hamiltonian ($H$). Please refer to
Appendix for the mathematics. –

$\displaystyle\bm{r}=x\hat{x}+\tau v_{x}\hat{v}_{x},$ (5a)
$\displaystyle\bm{\mathcal{V}}=\frac{d\bm{r}}{dt}=\dot{x}\hat{x}+\tau\dot{v}_{x}\hat{v}_{x}.$
(5b) The flux density of the flow in the phase-space $\bm{J}$ will then be the
multiple of the particle density and velocity field of the phase-space –
$\displaystyle\bm{J}=\eta\bm{\mathcal{V}}=\frac{1}{\tau}f\dot{x}\hat{x}+f\dot{v}_{x}\hat{v}_{x}.$
(5c)

###  II.1 The kinetic Euler equations

The velocity and flux density fields of the particle flow in phase-space,
described in equations (5a-c) can now be used to understand the nature of the
flow of this phase-space fluid. We begin by studying the time-rate of change
of the particle density in the phase-space –

$\displaystyle\frac{d\eta}{dt}=\frac{1}{\tau}\left(\frac{\partial f}{\partial
t}+\frac{\partial f}{\partial x}\frac{dx}{dt}+\frac{\partial f}{\partial
v_{x}}\frac{dv_{x}}{dt}\right).$ (6a) Also,
$\displaystyle\nabla\cdot\bm{\mathcal{V}}=\frac{\partial\dot{x}}{\partial
x}+\frac{\partial\dot{v}_{x}}{\partial v_{x}}=\frac{\partial}{\partial
x}\frac{\partial H}{\partial v_{x}}-\frac{\partial}{\partial
v_{x}}\frac{\partial H}{\partial x}=0.$ (6b)
$\displaystyle\begin{split}\nabla\cdot\bm{J}=\eta\nabla.\bm{\mathcal{V}}+\bm{\mathcal{V}}.\nabla\eta=\dot{x}\frac{\partial\eta}{\partial
x}+\tau\dot{v}_{x}\frac{\partial\eta}{\tau\partial v_{x}}\\\
=\frac{1}{\tau}\dot{x}\frac{\partial f}{\partial
x}+\frac{1}{\tau}\dot{v}_{x}\frac{\partial f}{\partial v_{x}}.\end{split}$
(6c) $\displaystyle\Rightarrow\frac{d\eta}{dt}=\frac{\partial\eta}{\partial
t}+\nabla\cdot\bm{J}=\begin{cases}0&\text{collision-less system}\\\
\left(\frac{\delta f}{\delta t}\right)_{\text{coll.}}&\text{with
collisions}\end{cases}$ (6d)

(a)

(b)

(c)

Figure 1: Velocity field $\bm{\mathcal{V}}$ plots of particle phase-space of
an ideal system (a) with Maxwellian (thermal equilibrium) distribution and (b)
in the presence of an negative bell-like attractive potential with streamlines
showing vortical motion of the phase-space fluid resembling the phase-space
trajectories of trapped (blue) and free (red) phase-space fluid elements, (c)
Unsteady flow in presence of hyperbolic secant potential waves, flow
represented by various streamlines.

Equation (6d) describes the continuity equation of the phase-space flow in the
cases of collision-less and collisional considerations of the system Gibbs,
(1902); Vlasov, (1938). The flow of the phase-space fluid representation of
the one-dimensional system is in-compressible in a collision-free condition.
The dynamics of the phase-space fluid analogue of the system can be further
explored by framing Euler phase-space momentum equations governing its flow.
We do this by determining the time-derivatives of the fluid’s velocity field.
For the same, we introduce a general internal force field $\bm{F}(x)$ and
kinetic pressure222The pressure term in the above equation is the kinetic
equivalent of fluid pressure present in the phase-space of the system and
arises (in fluid-analogue) due to particle density gradient in the phase-space
of the particle, along the position axis. This pressure term must not be
confused with the plasma fluid pressure which does not play a role in the
above equations. $\bm{P}$ into the system.

$\displaystyle\frac{d\bm{\mathcal{V}}}{dt}=\ddot{x}\hat{x}+\tau\ddot{v}_{x}\hat{v}_{x}.\quad\text{Introducing
advection terms,}$ (7a)
$\displaystyle\ddot{x}(x,t)=\frac{\partial\dot{x}}{\partial
t}+v_{x}\frac{\partial\dot{x}}{\partial
x},\quad\ddot{v}_{x}(x,t)=\frac{\partial\dot{v}_{x}}{\partial
t}+v_{x}\frac{\partial\dot{v}_{x}}{\partial x}$ (7b) For a system experiencing
an internal force field $\bm{F}(x)$ and pressure $\bm{P}$,
$\displaystyle\ddot{x}=\frac{\bm{F}(x)}{m}-\frac{\nabla_{x}\bm{P}}{m\eta},\qquad\ddot{v}_{x}=\frac{\dot{\bm{F}}(x)}{m}.$
(7c) Splitting into components along each axis in phase-space, we get –
$\displaystyle\frac{\partial\dot{x}}{\partial
t}+v_{x}\frac{\partial\dot{x}}{\partial
x}=\frac{\bm{F}(x)}{m}-\frac{1}{m\eta}\frac{\partial\bm{P}}{\partial x}$ (7d)
$\displaystyle\frac{\partial\dot{v}_{x}}{\partial
t}+v_{x}\frac{\partial\dot{v}_{x}}{\partial x}=\frac{\dot{\bm{F}}(x)}{m}.$
(7e)

The above equations (7d, e) describe the momentum equations for the flow of
the phase-space fluid model of the system.

Figure 2: Velocity magnitude
$|\mathcal{V}|=\sqrt{\dot{x}^{2}+\tau^{2}\dot{v}_{x}^{2}}$ of the phase-space
fluid in an electron hole region in logarithmic scale. The central, nearly
stationary, positive core is clearly visible surrounded by a trapped electron
current revolving around it.

For different particles, the equations can be generalised as –

$\begin{split}m_{i}\eta_{i}\left(\frac{\partial\dot{x}}{\partial
t}+v_{x}\frac{\partial\dot{x}}{\partial
x}\right)=\eta_{i}\bm{F_{i}}(x)-\nabla_{x}\bm{P},\\\
m_{i}\eta_{i}\left(\frac{\partial\dot{v}_{x}}{\partial
t}+v_{x}\frac{\partial\dot{v}_{x}}{\partial
x}\right)=\eta_{i}\dot{\bm{F_{i}}}(x).\end{split}$ (8)

One can introduce more terms on the right side of the equations (8) for other
interactions between particles in phase-space. Equations (6d) and (8) together
describe the kinetic Euler equations for the fluid model of the particle
phase-space. These equations are in capacity to collectively describe the
complete dynamics of the system in its phase-space, in the same nature as the
classical Euler equations describe the behaviour of an ordinary fluid. For the
collision-less condition, equation (6d) can be set to zero.

###  II.2 Collision-less diffusion in phase-space

A form of particle diffusion in this phase-space fluid model can be observed
as an outcome of the developed kinetic Euler equations. For this, we study the
kinetic momentum equation along position-space, with a time-invariant, quasi-
static flow of the fluid in phase-space, therefore equating
$\frac{\partial\dot{x}}{\partial t}$ to 0. We assume a thermal equilibrium
state, such that the pressure gradient $\nabla\bm{P}_{i}$ can be expressed in
the isothermal, barotropic form $K_{B}T_{i}\nabla\eta_{i}$, $T_{i}$ being the
temperature of the particle and $K_{B}$ the Boltzmann constant. Using
$\dot{x}=v_{x}$, the particle velocity, the equation becomes –

$\displaystyle m_{i}\eta_{i}v_{x}\frac{\partial v_{x}}{\partial
x}=\eta_{i}\bm{F_{i}}(x)-K_{B}T_{i}\nabla_{x}\eta_{i}$ (9a)
$\displaystyle\bm{J}_{x_{i}}=\eta_{i}v_{x}=\frac{\eta_{i}\bm{F}_{i}(x)}{m_{i}\frac{\partial
v_{x}}{\partial x}}-\frac{K_{B}T_{i}}{m_{i}\frac{\partial v_{x}}{\partial
x}}\nabla_{x}\eta_{i}.$ (9b)

The absence of a force field $\bm{F}(x)$ in the system produces a non-
diverging, non-rotating flow of the phase-space fluid, restricted along the
position space, since the fluid velocity along $\hat{v}_{x}$ direction becomes
zero (refer to FIG. 1a). In such a case, $\frac{\partial v_{x}}{\partial x}$
becomes equal to zero. However, in regions with non-zero force field, the
phase-space fluid velocity $\bm{\mathcal{V}}$ may retain a rotational
velocity, thus exhibiting a solenoidal flow (refer to FIG. 1b).

$\Bigl{|}\nabla\times\bm{\mathcal{V}}\Bigr{|}=\frac{\partial\tau\dot{v}_{x}}{\partial
x}-\frac{1}{\tau}\frac{\partial\dot{x}}{\partial
v_{x}}=\frac{\tau}{m_{i}}\frac{\partial\bm{F}(x)}{\partial
x}-\frac{1}{\tau}\frac{\partial v_{x}}{\partial v_{x}}.$ (10)

Here, $\dot{v}_{x}$ has been substituted by $\bm{F}/m_{i}$. The presence of a
force field also suggests that –

$\allowdisplaybreaks\begin{split}\nabla.\bm{\mathcal{V}}=\frac{\partial\dot{x}}{\partial
x}+\frac{1}{m_{i}}\frac{\partial\bm{F}}{\tau\partial
v_{x}}=0\Rightarrow\frac{\partial v_{x}}{\partial
x}=-\frac{1}{m_{i}}\frac{\partial\bm{F}(x)}{\partial x}\frac{\partial
x}{\partial v_{x}}\\\ \Rightarrow\frac{\partial v_{x}}{\partial
x}=\sqrt{\frac{1}{m_{i}}(-\nabla_{x}\bm{F})}.\end{split}$ (11)

Including equation (11) in equation (9b), we get –

$\displaystyle\begin{split}&\bm{J}_{x_{i}}=\eta_{i}v_{x}=\frac{\eta_{i}\bm{F}_{i}(x)}{\sqrt{-m_{i}\nabla_{x}\bm{F}}}-\frac{K_{B}T_{i}}{\sqrt{-m_{i}\nabla_{x}\bm{F}}}\nabla_{x}\eta_{i},\\\
&\text{or,}\quad\bm{J}_{x}=\mu\eta_{i}\bm{F}(x)-\mathcal{D}\nabla_{x}\eta_{i}\end{split}$
(12a)
$\displaystyle\text{where,}\quad\mu=\sqrt{\frac{-1}{m_{i}\nabla_{x}\bm{F}(x)}},\quad\mathcal{D}=K_{B}T_{i}\sqrt{\frac{-1}{m_{i}\nabla_{x}\bm{F}(x)}}$
(12b) $\displaystyle\text{and}\qquad K_{B}T_{i}\mu=\mathcal{D}.$ (12c)

In the above equations, $\mu$ represents particle mobility in phase-space
along position ($x$) axis. $\mathcal{D}$ represents the fluid diffusion
coefficient in phase-space, which is related to the mobility $\mu$ using
Einstein’s classical diffusion relation Einstein, (1905); Sutherland, (1905)
shown in equation (12c). Both particle mobility and diffusion coefficient are
facilitated by the presence of a negative force-field gradient in these
regions. Further understanding of this phenomena is therefore possible by
introducing field-specific equations which describe this gradient. As an
example, upon introducing an electrostatic interaction, $\bm{F}=q\bm{E}$ where
$\bm{E}(x)$ is the electric field of this interaction, the field gradient will
then become the spatial charge density, since,
$\epsilon_{0}\nabla_{x}\bm{E}=\rho(x)$, suggesting that in an electrostatic
domain, such a diffusion occurs in the presence of a charge accumulation in
space.

###  II.3 Vorticity field in phase-space

We now proceed to describe the vorticity field $\bm{\xi}(x,\tau v_{x})$
associated with the phase-space fluid analogue of a system, which is simply
the curl of the velocity field, as described in equation (10).

$\begin{split}\bm{\xi}=\nabla\times\bm{\mathcal{V}}=\left(\frac{\partial\tau\dot{v}_{x}}{\partial
x}-\frac{1}{\tau}\frac{\partial\dot{x}}{\partial v_{x}}\right)\hat{n}\\\
=\left(\frac{\tau}{m_{i}}\frac{\partial\bm{F}(x)}{\partial
x}-\frac{1}{\tau}\right)\hat{n},\\\
\hat{n}=\hat{x}\times\hat{v}_{x}.\end{split}$ (13)

Here, $\hat{n}$ is a direction normal to the phase-space plane, and is a
constant in a non-relativistic inertial frame333The relativistic
transformations of position and velocity coordinates of the particle states
and its affect on the preceding mathematics can be explored in future works on
this field of study.. It is observable that vorticity of the phase-space
fluid, which varies only along the position-space, is influenced by the
present force-field gradient, without which it remains a constant. It is hence
in the nature of this phase-space fluid model of a system to exhibit some type
of local deformation, sheering or rotational, in time. Presence of vorticities
in either direction, regulated by the direction of the force-field gradient,
can subsequently give rise to local turbulence which can then result in
formation of various characteristic structures, such as vortices (refer to
FIG. 1c).
The vorticity field of the phase-space fluid model of the system can be
further analysed using the vorticity transport equation, which can be
simplified by introducing the barotropic pressure term and restricting within
two-dimensions –

$\frac{\partial\bm{\xi}(x,t)}{\partial
t}+(\bm{\mathcal{V}}\cdot\nabla)\bm{\xi}=\mathcal{D}\nabla^{2}\bm{\xi}+\nabla\times\left(\frac{\mathcal{F}}{m\eta}\right).$
(14)

Here, $\mathcal{F}$ represents the sum-total of all body forces acting on the
system. Further simplification of this equation can be performed upon defining
the vorticity term specific to a given form of interaction between particles
in the system.
The preceding mathematics in this section describes the phase-space of a many-
particle system as a two-dimensional fluid. The Euler-like equations developed
in this section describe flow-like behaviour of the particle phase-space and
can now be used to analyse for vortical structures. We use this technique for
collision-less plasmas and identify electron holes and ion holes as vortices
in the phase-space of the plasma in the next section.

##  III Analytical identification of electron and ion holes as phase-space
vortices

Having described the phase-space as a two-dimensional dynamic fluid, we can
now analytically identify electron and ion holes as vortices in the phase-
space. For the same, we use some established techniques of vortex
identification in two-dimensional fluids Helmholtz, (1867); Lugt, (1979); Tian
et al., (2018). We analyse the phase-space vorticity fields and observe
particle trajectories in the phase-space of the plasma, and observe local
concentrations of vorticities as well as locally rotating fluids.
Upon setting particle charge $q=\delta e$ with $\delta=\pm 1$ for (singly
ionised, positive or negative) ions or electrons, $\bm{F}=qE$ with $E(x,t)$ as
electric field, $\dot{x}=v_{x}$ and $\dot{v}_{x}=\frac{\delta e}{m}E(x,t)$ in
equation (13), the vorticity associated with a plasma particle (ion or
electron) phase-space fluid flow becomes –

$\allowdisplaybreaks\begin{split}\bm{\xi}=&\Bigl{(}\tau\frac{\partial\dot{v}_{x}}{\partial
x}-\frac{1}{\tau}\Bigr{)}\hat{n}=\Bigl{(}\frac{\delta e\tau}{m}\frac{\partial
E(x,t)}{\partial x}-\frac{1}{\tau}\Bigr{)}\hat{n}\\\ =&\Bigl{(}\frac{\delta
e\tau}{\epsilon_{0}m}\rho(x,t)-\frac{1}{\tau}\Bigr{)}\hat{n}=\omega_{p}\Bigl{(}\delta\bar{\rho}(x,t)-1\Bigr{)},\\\
&\text{where,}\quad en_{0}\bar{\rho}(x,t)=\rho(x,t).\end{split}$ (15)

Here, $\tau=\sqrt{\frac{\epsilon_{0}m}{n_{0}q^{2}}}$ representing (inverse)
plasma particle frequency $\omega_{p}$, which can be set to $1$ for even
further simplification. It is a simple task to introduce multiple-ionised
particles into these equations by replacing $\delta e$ with $\delta Ze$, $Z$
being the ionisation number. As shown in equation (15), vorticity associated
with the plasma particle phase-space fluid flow is related to the local
concentration of charge in a region of phase-space, defined by the spatial
charge density $\rho(x)$.
The vorticity $\bm{\xi}$ associated with the phase-space, like any other
ordinary fluid, describes the deformations in a region of the particle phase-
space fluid. These deformations can be local sheering or rotational
deformations, the latter of which describes a local spin of fluid elements
about a common axis and is identified as a fluid vortex. These can also be
defined as a local region of concentrated vorticity distinguished by closed or
spiral streamlines Lugt, (1979), or as described by Tian et al., (2018), using
the vortex parameter $g_{z}$ (or $g_{z\theta\rvert_{max}}$ in an invariant
form) which can be defined for our system as –

$\displaystyle\begin{split}&g_{z}=\frac{\tau\partial\dot{v}_{x}}{\partial
x}\frac{\partial\dot{x}}{\tau\partial v_{x}}=\frac{\delta e}{m}\frac{\partial
E(x)}{\partial x}\frac{\partial v_{x}}{\partial v_{x}}\\\ &=\frac{\delta
e}{m\epsilon_{0}}\rho(x)=\delta\omega_{p}^{2}\bar{\rho}(x)\end{split}$ (16a)
$\displaystyle g_{z\theta}\lvert_{\max}=\alpha^{2}-\beta^{2}$ (16b)
$\displaystyle\begin{split}=\frac{1}{4}\Biggl{[}\left(\frac{\tau\partial\dot{v}_{x}}{\tau\partial
v_{x}}-\frac{\partial\dot{x}}{\partial
x}\right)^{2}+\left(\frac{\tau\partial\dot{v}_{x}}{\partial
x}+\frac{\partial\dot{x}}{\tau\partial v_{x}}\right)^{2}\\\
-\left(\frac{\tau\partial\dot{v}_{x}}{\partial
x}-\frac{\partial\dot{x}}{\tau\partial v_{x}}\right)^{2}\Biggr{]}\end{split}$
(16c)
$\displaystyle=\frac{1}{4}\Biggl{[}\Bigl{(}\frac{\partial\dot{x}}{\partial
x}+\frac{\partial\dot{v}_{x}}{\partial
v_{x}}\Bigr{)}^{2}-4\frac{\partial\dot{x}}{\partial
x}\frac{\partial\dot{v}_{x}}{\partial
v_{x}}+4\frac{\partial\dot{v}_{x}}{\partial x}\frac{\partial\dot{x}}{\partial
v_{x}}\Biggr{]}$ (16d)
$\displaystyle\begin{split}=\Biggl{(}\frac{\partial\dot{v}_{x}}{\partial
x}\frac{\partial\dot{x}}{\partial v_{x}}-\frac{\partial\dot{x}}{\partial
x}\frac{\partial\dot{v}_{x}}{\partial v_{x}}\Biggr{)}=\\\
\Biggl{(}\frac{\delta e}{m}\frac{\partial E(x)}{\partial x}\frac{\partial
v_{x}}{\partial v_{x}}-\frac{\partial v_{x}}{\partial x}\frac{\delta
e}{m}\frac{\partial E(x)}{\partial v_{x}}\Biggr{)}\\\ =\frac{\delta
e}{m}\frac{\partial E(x)}{\partial x}=\frac{\delta
e}{m\epsilon_{0}}\rho(x)=\delta\omega_{p}^{2}\bar{\rho}(x)=g_{z}.\end{split}$
(16e) From ref. Tian et al., (2018), for a vortex region, $\displaystyle
g_{z\theta}\lvert_{\max}<0\implies\delta\omega_{p}^{2}\bar{\rho}(x)<0\implies\delta\bar{\rho}(x)<0.$
(16f)

Figure 3: Phase-space portrait of the Q-machine Plasma showing electron hole
formed using potential pulse with amplitude $A=3.0\frac{K_{B}T_{e}}{e}$,
showing variation of phase-space density. Figure generated using kinetic
simulation of the Q-machine experiment as performed by Saeki et al., (1979).

The results described in equation (16a,e,f) state the relation between the
local accumulation of charge in position-space and the rotational deformation
caused by this accumulation in the phase-space fluid. It is well-known that
phase-space electron and ion holes produce a well-defined local depression of
particle charge in their regions. This would be in complete agreement with
equation (16f) since a vortical presence will be detected in any region of
reduced particle charge, i.e.,

$g_{z\theta\rvert_{max}}=\begin{cases}\bar{\rho}(x)>0&\text{for a -ve particle
phase-space vortex}\\\ \bar{\rho}(x)<0&\text{for a +ve particle phase-space
vortex }\end{cases}$ (17)

A key phenomena of this vortex region in the phase-space fluid model and its
correspondence to the trapping condition must be included in the preceeding
analysis – the presence of fluid elements with both low and high velocities in
the phase-space plane. Like any ordinary fluid, in the region of a vortex,
fluid particles with low velocity produce complete rotational deformation,
indicating the ‘trapping’ of particles in the vortex. However, fluid elements
with higher fluid velocity enter and exit the vortex region showing an arc-
like trajectory. Clearly, this trapping will be determined by the kinetic
energy of the fluid element in comparison to the potential $\Phi$ developed
due to the local accumulation of charge in the vortex, as shown in equation
(16a,e), i.e.,

$\displaystyle\frac{1}{2}mv_{x}^{2}+\delta e\Phi\geq 0\implies\text{particles
escape the vortex.}$ (18a) $\displaystyle\frac{1}{2}mv_{x}^{2}+\delta
e\Phi<0\implies\text{particle trapping.}$ (18b)

The conditions defined by equations (16f), (17) and (18) can be used to
identify vortices in the phase-space fluid (see Fig. 4).
We therefore conclude in this section that electron holes and ion holes are
in-fact vortices in the phase-space of the collision-free plasma. Having found
this analytical proof of their vortical nature, we can now approach the
problem of electron phase-space holes using the fluid-like kinetic approach
described in sec. II of this article. We do this in the succeeding section
wherein we apply the phase-space-fluid model to study the structure of an
electron phase-space hole.

(a)

(b)

(c)

Figure 4: (a) Vorticity of phase-space electron fluid, showing concentrated
vorticity at the cite of the hole (see FIG. 3). (b) Streamline plots of
electron phase-space fluid elements showing rotational streamline (red) in the
region of the electron hole. (c) Map showing regions of rotational deformation
(electron hole encircled in red) with $g_{z\theta\rvert_{max}}<0$. Note that
other regions also identified as vortices due to formed turbulence produced by
the travelling compression and rarefaction waves.

##  IV Application of phase-space fluid model to derive the phase-space
density equations for an electron hole and the validity of schamel-df
equations

As an application of the phase-space-fluid model towards the study of electron
phase-space vortices, we determine a steady-state equation for the phase-space
density structure of this vortical region in the fluid phase-space. For the
same, we use the fluid-momentum equation derived in equation (8) along-with
the isothermal pressure gradient introduced in equation (9a). In the region of
a vortex, the temperature $T=T_{e}^{\text{tr.}}$, which is the temperature of
the trapped electrons.

$\eta m\left(\frac{\partial v_{x}}{\partial
t}+v_{x}\cdot\nabla_{x}v_{x}\right)=\eta
qE(x)-K_{B}T_{e}^{\text{tr.}}\nabla_{x}\eta$ (19)

We aim to find the steady state solution444Here, we have used the momentum
equation along position-space (19) which seems to be a more descriptive
relation between the phase-space density function $f(x,v_{x})$ and its
coordinates, and not the Vlasov equation (6), which is the continuity equation
of the phase-space fluid. of this equation in the region of the electron
vortex, by setting the material derivative of the phase-space fluid velocity
(left side of equation (19)) equal to zero, as follows –

$\displaystyle\eta
m_{e}\frac{D(v_{x})}{Dt}=0\Rightarrow\frac{q}{m_{e}\tau}fE(x)=\frac{1}{\tau}\nabla_{x}f\frac{K_{B}T_{e}^{\text{tr.}}}{m_{e}}$
(20a) $\displaystyle\Rightarrow-fq\frac{\partial\phi(x)}{\partial
x}=K_{B}T_{e}^{\text{tr.}}\frac{\partial f}{\partial x}$ (20b) $\displaystyle-
fq\frac{\partial\phi(x)}{\partial x}dx=K_{B}T_{e}^{\text{tr.}}\frac{\partial
f}{\partial x}dx\Rightarrow-fqd\phi=K_{B}T_{e}^{\text{tr.}}df$ (20c)
$\displaystyle\Rightarrow\int\frac{df}{f}=-\frac{q}{K_{B}T_{e}^{\text{tr.}}}\int\frac{\partial\phi}{\partial
x}dx=-\frac{q}{K_{B}T_{e}^{\text{tr.}}}\int d\phi.$ (20d)

For the above integration, we set the limits from within the vortex region (
$0\leq v_{x}<\sqrt{2e\phi/m_{e}}$) till its boundary
($v_{x}=\sqrt{2e\phi/m_{e}}$). For a constant Hamiltonian $H\leq 0$, change in
potential can be expressed in terms of particle kinetic energy ($\mathcal{K}$)
as follows –

$\displaystyle-
qd\phi=d\mathcal{K}\Rightarrow-q\Delta\phi=\Delta\mathcal{K}=q\phi-\frac{1}{2}mv_{x}^{2}$
(21a) Inserting in equation (20b) and normalising
$\phi=\chi(K_{B}T_{e}q^{-1})$, $\mathcal{K}=\bar{v}_{x}^{2}K_{B}T_{e}$, we get
– $\displaystyle\begin{split}\ln f_{\text{bd.}}-\ln
f_{\text{tr.}}=-\frac{1}{K_{B}T_{e}^{\text{tr}}}\left(\frac{1}{2}mv_{x}^{2}-q\phi\right)\\\
=-\frac{K_{B}T_{e}}{K_{B}T_{e}^{\text{tr}}}\left(\bar{v}_{x}^{2}-\chi\right)\end{split}$
(21b) $\displaystyle\Rightarrow\ln f_{\text{tr}}=-\left[(-\ln
f_{\text{bd.}})+\Bigl{(}-\frac{T_{e}}{T_{e}^{\text{tr}}}\Bigr{)}(\bar{v}_{x}^{2}-\chi)\right],\text{
or,}$ (21c) $\displaystyle\begin{split}\ln
f_{\text{tr}}=-\left(M^{2}+\beta(\bar{v}_{x}^{2}-\chi)\right),\quad\text{where,}\quad\\\
M^{2}=-\ln
f_{\text{bd.}},\qquad\beta=-\frac{T_{e}}{T_{e}^{\text{tr}}}.\end{split}$ (21d)
$\displaystyle\Rightarrow\bar{f}_{\text{tr}}=\frac{1}{\sqrt{\pi}}\exp\left[-\bigl{\\{}M^{2}+\beta(\bar{v}_{x}^{2}-\chi)\bigr{\\}}\right].$
(21e)

Using the same approach around the hole, it can be shown that the phase-space
density remains a shifted Maxwellian in the un-trapped region –

$\bar{f}_{\text{free}}=\frac{1}{\sqrt{\pi}}\exp\Bigl{[}-\left(\bar{M}^{2}+\bar{v}_{x}^{2}-\chi\right)\Bigr{]},\qquad\bar{v}_{x}^{2}\geq\chi.$
(22)

The mathematics presented in equations (21a-e, 22) analytically derive the
schamel-dfSchamel, (1979) equations for an electron hole. These equations,
which describe the phase-space density $f$ for both trapped and free electrons
and were initially assumed B.G.K. solutions to the Vlasov-Poisson system, can
be very-well derived upon treating the electron phase-space of collision-less
plasmas as a fluid.
The equations (21a-e) also provide a relation between the normalised electron
hole speed $M$ in the reference frame of the perturbation and its structure in
phase-space (refer to FIG. 5 ). The normalised electron hole speed $M$ is, as
described in equation (21d), equal to the square root of the phase-space
density of the electron hole at its boundary, i.e.,

$M=\sqrt{-\ln(f^{\text{tr.}}_{\text{bd.}})}=\sqrt{-\ln\Bigl{[}\sqrt{\pi}\bar{f}^{\text{tr.}}\bigl{(}{x,v_{x}=\sqrt{{2e\phi}/{m_{e}}}}\bigr{)}\Bigr{]}}.$
(23)

Third, the definition of the particle trapping parameter $\beta$, as presented
by Schamel in his works on particle trapping solitary B.G.K. waves Schamel,
(1971, 1972, 1975, 1979) can be extracted as a ratio of the plasma electron
temperature to trapped electron temperature, as shown in equation (21d). This
is also in well-agreement to the initially assumed relation Hutchinson,
(2017). Upon further analysis of the parabolic structure of this phase-space
vortex, as clearly depicted from equation (21c), a relation between this
$\beta$ parameter and the depth of this vortex in the phase-space logarithmic
surface can be deduced –

$\displaystyle\log_{e}\sqrt{\pi}\bar{f}_{tr}=-M^{2}+\beta\chi-\beta\bar{v}_{x}^{2}$
(24a) The minima of this parabolic surface is equal to $\ln
f_{\text{min.}}^{\text{tr.}}$ at the centre of the hole, for which
$\chi=\chi_{0}$ and $\bar{v}_{x}^{2}=0$, with $\chi_{0}$ being the amplitude
of the hole potential. Therefore, $\displaystyle\beta=\frac{\ln
f_{\text{min.}}^{\text{tr.}}+M^{2}}{\chi_{0}}=\frac{1}{\chi_{0}}\ln\left(\frac{f_{\text{min.}}^{\text{tr.}}}{f_{\text{bd.}}^{\text{tr.}}}\right)=\frac{\Gamma}{\chi_{0}},\quad\Gamma<0.$
(24b)

The term
$\Gamma=\ln{f_{\text{min.}}^{\text{tr.}}}-ln{f_{\text{bd.}}^{\text{tr.}}}$ in
the above expression has been used to represent the (negative) depth of this
phase-space vortex in the $\log_{e}f(x,v_{x})$ surface. Thus,

$|\beta\chi_{0}|=|\Gamma|$ (25)

It is clear from this relation as to why the electron trapping parameter must
be negative – the value of $f_{\text{min.}}^{\text{tr.}}$ lies between $0$ and
$f_{\text{bd.}}^{\text{tr.}}$, implying that $\Gamma\in(-\infty,0]$. Hence,
$\beta$ must always be a non-positive value, its minima representing a hole
with an empty core and maxima describing a flat surface electron hole. It is
also clear from equation (25) that $\beta$ does not represent the hole depth
in the phase-space plane ($-\Gamma$)Schamel, 2023b , but is actually the ratio
of hole phase-space depth to the hole potential amplitude.
Upon moving across an electron phase-space hole along position-space, one can
observe that the hole phase-space depth ($-\Gamma$) must be directly
proportional to the hole potential, keeping a constant $\beta$ for the hole.
Thus, for a single hole,

$\chi_{0}\propto-\Gamma.$ (26)

The same is shown in the numerical study of a solitary hole formed using the
Q-machine plasma kinetic simulationSaeki et al., (1979) with a pulse amplitude
of $4.3K_{B}T_{e}e^{-1}$ (refer to FIG. 6).

Figure 5: Hole speeds measured using kinetic simulation study of electron-hole
formation in Q-machine plasma Saeki et al., (1979); Lynov et al., (1979,
1980). Relation between hole speeds recorded in simulation and the same
calculated using relation (23) (red points), fitted curve using quadratic
fitting (red line) and its comparison to the $x=y$ curve(black line). Figure
6: Variation of hole phase-space depth with hole potential across position-
space (red curve) for a single hole. $1/|\beta|$ is the slope of the curve
(approximated by the black line).

##  V Discussions and Conclusion

In this article, a fluid outlook of the particle phase-space has been
presented. This approach, as has been shown, analytically confirms the
vortical nature of phase-space holes. This implies the validity not just of
their vortical behaviour, but also of the applicability of a fluid-like
treatment in the study of these phase-space vortices using a different
mathematical framework – the kinetic-Euler equations – which use the fluid-
like nature of the phase-space to describe the structure of the holes. It
confirms that, when the particle phase-space is treated as a two-dimensional
fluid, a velocity field and a vorticity field associated with the flow-like
behaviour of the particles (electrons and ions) in their respective phase-
space can be defined. Various contortions can then be analysed using these
fields and analogous to the analytical techniques of two-dimensional fluid-
dynamics, the phase-space characteristics of such kinetic structures can then
be studied.
The study of electron holes using this approach provides an exact derivation
of the schamel-df equations Schamel, (1979) which depict the phase-space
structure of electron holes. It reproduces the solutions both for free and
trapped particle densities and derives the definition of the particle trapping
$\beta$ parameter as it was originally defined (as the inverse, normalised
trapped electron temperature Hutchinson, (2017)). It also provides an
analytical relation between electron hole speed $M$, potential amplitude
$\chi_{0}$, the trapping parameter $\beta$, and phase-space hole depth
($-\Gamma$), expressing $\beta$ as a ratio of $\Gamma$ and $\chi_{0}$. This is
a measurable definition of the electron trapping parameter and
simplifies555The measurement of hole trapping parameter can be directly
calculated from simulation data of phase-space particle density minima and
reduces both the analytical as well as computational burden of measuring the
ratio of trapped and free electron temperatures. the dispersion relation of
the hole to

$\beta=\beta(\chi_{0},-\Gamma),\quad\Gamma=\Gamma(M,\ln
f_{\text{min.}}^{\text{tr.}}).$ (27)

Hence, this approach provides a precursor to the B.G.K. differential approach
– allowing one to not assume but derive a well-defined particle distribution
function for the study of the phase-space vortices using the pseudo-potential
method Schamel, (1979, 1986, 2012); Schamel, 2023b .
It is also noteworthy to observe that the electron hole B.G.K. wave is
observed as a vorticity wave travelling in the phase-space fluid, and
therefore can be studied using the vorticity transport relation presented in
equation (14), with the diffusion coefficient $\mathcal{D}$ defined in
equation (12b), which can be further calculated using the exact solution to
the phase-space structure in the trapped region of electron holes which has
been already derived. Hence, an alternate approach to the pseudo-potential
method is presented as an outcome of the fluid-phase-space model.
One must note that the stochastic nature of the studied kinetic phenomena and
the abundance in the trapping mechanisms can also be accounted for in this
model upon further exploration and modifications of the kinetic-Euler
equations. One such modification is the rejection of the thermal equilibria of
the trapped particles and adoption of a non-equilibrium state, allowing a
time-evolution and stabilisation of the electron-hole structure. Such
phenomena can also be explored as a part of further development of this model,
which we leave for future study in this domain.

## References

  * Aravindakshan et al., (2021) Aravindakshan, H., Kakad, A., Kakad, B., and Yoon, P. H. (2021).  Structural Characteristics of Ion Holes in Plasma.  Plasma 2021, Vol. 4, Pages 435-449, 4(3):435–449.
  * Bernstein et al., (1957) Bernstein, I. B., Greene, J. M., and Kruskal, M. D. (1957).  Exact nonlinear plasma oscillations.  Physical Review, 108(3):546.
  * Chen et al., (2005) Chen, L. J., Pickett, J., Kintner, P., Franz, J., and Gurnett, D. (2005).  On the width-amplitude inequality of electron phase space holes.  Journal of Geophysical Research: Space Physics, 110(A9):9211.
  * Einstein, (1905) Einstein, A. (1905).  Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen.  Annalen der Physik, vol. 4, t. 17.
  * Gibbs, (1902) Gibbs, J. W. (1902).  Elementary principles in statistical mechanics : developed with especial reference to the rational foundation of thermodynamics : Gibbs, J. Willard (Josiah Willard), 1839-1903 : Free Download & Streaming : Internet Archive.
  * Guio et al., (2003) Guio, P., Børve, S., Daldorff, L. K. S., Lynov, J. P., Michelsen, P., Pécseli, H. L., Rasmussen, J. J., Saeki, K., and Trulsen, J. (2003).  Nonlinear Processes in Geophysics Phase space vortices in collisionless plasmas.  Nonlinear Processes in Geophysics, 10:75–86.
  * Helmholtz, (1867) Helmholtz, H. (1867).  About integrals of hydrodynamic equations related with vortical motions.  Phil. Mag., 33(226):485–512.
  * Hutchinson, (2017) Hutchinson, I. H. (2017).  Electron holes in phase space: What they are and why they matter.  Physics of Plasmas, 24(5):055601.
  * Hutchinson and Zhou, (2016) Hutchinson, I. H. and Zhou, C. (2016).  Plasma electron hole kinematics. I. Momentum conservation.  Physics of Plasmas, 23(8):82101.
  * Iizuka and Tanaca, (1987) Iizuka, S. and Tanaca, H. (1987).  Mach number limit of electron holes in bounded beam-plasma system.  Plasma physics and controlled fusion, 29(2):133.
  * Lugt, (1979) Lugt, H. J. (1979).  The dilemma of defining a vortex.  Theoretical and Experimental Fluid Mechanics, pages 309–321.
  * Lynov et al., (1980) Lynov, J.-P., Michelsen, P., Pécseli, H. L., and Rasmussen, J. J. (1980).  Interaction between electron holes in a strongly magnetized plasma.  Physics Letters A, 80(1):23–25.
  * Lynov et al., (1979) Lynov, J.-P., Michelsen, P., Pécseli, H. L., Rasmussen, J. J., Saeki, K., and Turikov, V. A. (1979).  Observations of solitary structures in a magnetized, plasma loaded waveguide.  Physica Scripta, 20(3-4):328.
  * Morse and Nielson, (1969) Morse, R. L. and Nielson, C. W. (1969).  One-, two-, and three-dimensional numerical simulation of two-Beam plasmas.  Physical Review Letters, 23(19):1087–1090.
  * Saeki et al., (1979) Saeki, K., Michelsen, P., Pécseli, H. L., and Rasmussen, J. J. (1979).  Formation and coalescence of electron solitary holes.  Physical Review Letters, 42(8):501.
  * Schamel, (1971) Schamel, H. (1971).  Stationary solutions of the electrostatic Vlasov equation.  Plasma Physics, 13(6):491.
  * Schamel, (1972) Schamel, H. (1972).  Stationary solitary, snoidal and sinusoidal ion acoustic waves.  Plasma Physics, 14(10):905.
  * Schamel, (1975) Schamel, H. (1975).  Analytic BGK Modes and Their Modulational Instability.  Journal of Plasma Physics, 13(1):139–145.
  * Schamel, (1979) Schamel, H. (1979).  Theory of electron holes.  Physica Scripta, 20(3-4):336.
  * Schamel, (1982) Schamel, H. (1982).  Kinetic theory of phase space vortices and double layers.  Physica Scripta, 1982(T2A):228.
  * Schamel, (1986) Schamel, H. (1986).  Electron holes, ion holes and double layers: Electrostatic phase space structures in theory and experiment.  Physics reports, 140(3):161–191.
  * Schamel, (1998) Schamel, H. (1998).  Phase space vortices in unbunched beams in synchrotrons.  Physica Scripta T, 75(T75):23–27.
  * Schamel, (2012) Schamel, H. (2012).  Cnoidal electron hole propagation: Trapping, the forgotten nonlinearity in plasma and fluid dynamics.  Physics of Plasmas, 19(2):20501.
  * (24) Schamel, H. (2023a).  Pattern formation in Vlasov–Poisson plasmas beyond Landau caused by the continuous spectra of electron and ion hole equilibria.  Reviews of Modern Plasma Physics, 7(1):1–41.
  * (25) Schamel, H. (2023b).  Pattern formation in Vlasov–Poisson plasmas beyond Landau caused by the continuous spectra of electron and ion hole equilibria.  Reviews of Modern Plasma Physics, 7(1):1–41.
  * Sutherland, (1905) Sutherland, W. (1905).  LXXV. A dynamical theory of diffusion for non-electrolytes and the molecular mass of albumin.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 9(54):781–785.
  * Tian et al., (2018) Tian, S., Gao, Y., Dong, X., and Liu, C. (2018).  Definitions of vortex vector and vortex.  Journal of Fluid Mechanics, 849:312–339.
  * Trivedi and Ganesh, (2017) Trivedi, P. and Ganesh, R. (2017).  Driven phase space vortices in plasmas with nonextensive velocity distribution.  Physics of Plasmas, 24(3).
  * Turikov, (1978) Turikov, V. A. (1978).  Computer simulation of the formation of langmuir solitons and holes in a cylindrical magnetized plasma column.  Technical report, Risoe national lab.
  * Turikov, (1984) Turikov, V. A. (1984).  Electron phase space holes as localized BGK solutions.  Physica Scripta, 30(1):73–77.
  * Vlasov, (1938) Vlasov, A. A. (1938).  On vibration properties of electron gas.  Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki, 8(3):291.

## Appendix

The self-consistency of a fluid requires solving a Laplace equation of the
stream function $\mathcal{H}(x,y)$, such thatSchamel, (2012) –

$\mathcal{V}(x,y)=\nabla\mathcal{H}\times\hat{n},$ (28)

where, $\hat{n}$ is a direction normal to the two-dimensional fluid. The same,
in the case of the fluid-analogue of the particle phase-space, can be
described as a modified Hamiltonian $H$ of the particles –

$\mathcal{H}=\frac{\tau}{m}H=\frac{\tau}{m}\left(\frac{1}{2}mv_{x}^{2}+q\phi(x)\right).$
(29)

Here, $\phi(x)$ is the potential field experienced by the particle of charge
$q$. Inserting equation (29) in equation (28), we see that666Since
$\frac{d}{dx}\neq\partial_{x}$, the same can be implied to show that
$\partial_{x}(\frac{1}{2}mv_{x}^{2})=0$. –

$\displaystyle\nabla\mathcal{H}$
$\displaystyle=\frac{\tau}{m}\left(q\frac{\partial\phi}{\partial
x}\hat{x}+\frac{m}{2\tau}\frac{\partial v_{x}^{2}}{\partial
v_{x}}\hat{v}_{x}\right)$ (30a)
$\displaystyle=\frac{q\tau}{m}\frac{\partial\phi}{\partial
x}\hat{x}+v_{x}\hat{v}_{x}.$ (30b)
$\displaystyle\nabla\mathcal{H}\times\hat{n}$
$\displaystyle=\frac{q\tau}{m}\frac{\partial\phi}{\partial
x}(\hat{x}\times\hat{n})+v_{x}(\hat{v}_{x}\times\hat{n})$ (30c)
$\displaystyle=\frac{q\tau}{m}\left(-\frac{\partial\phi}{\partial
x}\right)\hat{v}_{x}+v_{x}\hat{x}=\mathcal{V}.$ (30d)

Therefore, upon using the relation between $\hat{n}$, $\hat{x}$ and
$\hat{v}_{x}$ defined in equation (13), equation (28) is satisfied. The stream
function (29) can then be acted upon by the Laplacian operator in phase-space
to reproduce the vorticity expression defined in equation (13) –

$\displaystyle\nabla\cdot\nabla\mathcal{H}$
$\displaystyle=\frac{q\tau}{m}\frac{\partial^{2}\phi}{\partial
x^{2}}+\frac{1}{\tau}\frac{\partial v_{x}}{\partial v_{x}}$ (31a)
$\displaystyle=-\left(\frac{\tau}{m}\frac{\partial\bm{F}(x)}{\partial
x}-\frac{1}{\tau}\right)=-\bm{\xi}(x),$ (31b) where,
$\displaystyle\quad\bm{F}(x)=-q\frac{\partial\phi}{\partial x}\hat{x}.$ (31c)

This shows the self-consistency of the fluid-analogue of the particle phase-
space.


# Recurrent Neural Network from Adder’s Perspective:
Carry-lookahead RNN

Haowei Jianga, Feiwei Qina, Jin Caob,∗, Yong Penga, Yanli Shaoa
aSchool of Computer Science and Technology, Hangzhou Dianzi University, China
bWhiting School of Engineering, Johns Hopkins University, USA

###### Abstract

The recurrent network architecture is a widely used model in sequence
modeling, but its serial dependency hinders the computation parallelization,
which makes the operation inefficient. The same problem was encountered in
serial adder at the early stage of digital electronics. In this paper, we
discuss the similarities between recurrent neural network (RNN) and serial
adder. Inspired by carry-lookahead adder, we introduce carry-lookahead module
to RNN, which makes it possible for RNN to run in parallel. Then, we design
the method of parallel RNN computation, and finally Carry-lookahead RNN (CL-
RNN) is proposed. CL-RNN takes advantages in parallelism and flexible
receptive field. Through a comprehensive set of tests, we verify that CL-RNN
can perform better than existing typical RNNs in sequence modeling tasks which
are specially designed for RNNs. Code and models are available at:
https://github.com/WinnieJiangHW/Carry-lookahead_RNN.

###### keywords:

Deep learning; Carry-lookahead; Parallel computation; Sequence modeling

††journal: Neural Networks

##  1 Introduction

Recurrent network architecture is a common model used to handle sequence
tasks. The title of sequence modeling chapter in the canonical textbook _Deep
Learning_ [1] is ”Sequence Modeling: Recurrent and Recursive Nets”. Recurrent
neural network (RNN) proposed by Rumelhart et al. [2] and its improved models
such as long short-term memory (LSTM) [3], gated recurrent neural networks
[4], are firmly established to deal with language modeling, machine
translation, etc. [5, 6, 7] Numerous researches on recurrent network
architecture push its boundaries and application in natural language
processing [8, 9].

However, recurrent neural networks rely on the preorder state $h_{t-1}$ to
generate the output and postorder state $h_{t}$. This serial dependency
hinders the parallelization of computation, which makes the model run
inefficient.

In the history of digital electronics, the adder as the most basic computing
unit has faced the same problem. The traditional serial adder has to wait for
the completion of the preorder computation in order to obtain the carry bit.
Then it generates next carry bit and delivers. Considering the serial adder is
analogous to the recurrent network architecture in various aspects, we suppose
that there is an analogy between them.

Digital electronics has developed into a mature subject nowadays. Meanwhile,
parallelization improvement has always been a research hotspot in the
discipline [10, 11, 12, 13, 14]. As one of the most fundamental components,
adder has been deeply studied. Given the analogy between adder and recurrent
network architecture, we suppose that it is valid to apply the proven
optimization experience on adder to improving recurrent neural network.

By extension, we suppose that the knowledge of digital electronics could be
informative to the design of neural network architectures. Neural network
architectures have commonalities with electronic circuits in terms of
structure and other characteristics. Based on that, the proven experience on
digital electronics can be applied to improving neural networks. For example,
shortcut connection used in ResNet [15] is similar to short circuit, which is
a common method used in electronic circuits; branch structure adopted by
GoogLeNet [16] is also a basic structure in electronic circuits.

Carry-lookahead adder [17] is today’s mainstream adder. It innovatively
proposes the carry-lookahead module and applies independent formulas to enable
the parallel calculation of multi-bit. Compared with the traditional serial
adder, it greatly improves the computation efficiency.

Inspired by carry-lookahead adder, we introduce carry-lookahead module to RNN.
Furthermore, the method of parallel RNN computation is designed and finally
Carry-lookahead RNN (CL-RNN) is proposed.

Through a comprehensive set of tests, CL-RNN not only achieves parallelization
compared to recurrent neural networks, but also shows stronger representation
capability in sequence modeling.

In summary, the main contributions of this paper are as follows.

  1. 1)

Consider the similarities between RNN and serial adder in various aspects.
Taking their analogy as the core point, we propose the idea of transferring
the optimization experience on adder to RNN.

  2. 2)

Inspired by carry-lookahead adder, we propose Carry-lookahead RNN, which
consists of carry-lookahead module and parallel RNN module. The carry-
lookahead module provides support for model’s parallelization, and the
parallel RNN module improves the representation capability and eases the model
optimization difficulty.

  3. 3)

Compare CL-RNN with classic RNNs as well as TCN. Through ablation experiment,
we verify the effectiveness of the model improvement. After a series of stress
tests, it’s found that CL-RNN can perform better than classic TCN, RNN and
LSTM in sequence modeling, which recurrent network architecture excel in.

##  2 Related Work

Since the classic recurrent neural network has been proposed by Rumelhart in
1986 [2], a lot of related researches are conducted and many variants are
proposed. Recurrent neural networks now have been a large family, and
extensively applied in language modeling, machine translation, etc. [5, 6, 7]

Some of the classic variants greatly improves the performance of RNN. Long
short-term memory [3] introduces the structure of ”gates” to RNN, which allows
RNN to forget irrelevant information and keep important information in long-
term; Gate Recurrent Unit [4] simplify LSTM to some extent, which accelerate
the training efficiency.

In recent years, more and more improved models are proposed. krueger et al.
regularize the penalty term of LSTM and propose Norm-stabilized LSTMs [18];
Vohra et al. superimpose Deep Belief Net [19] on LSTM and propose DBN LSTM
[20], which can perform very well on certain datasets. Other recent works
combine RNN with CNN, such as Convolutional LSTM [21], which replaces the
fully-connected layer in LSTM with CNN; Quasi-RNN [22], which interleaves
convolutional layers with recurrent layers; Dilated RNN [23], which introduces
dilation into recurrent network architecture.

The various variants proposed make the performance of RNN improve a lot.
However, RNNs are still limited by its serial dependency, which seriously
lowers the computation efficiency and hinders its large-scale applications.

Although some recent works have somewhat optimized the efficiency through
factorization techniques [24], conditional operations [25], etc., the
fundamental constraints of sequential computation still exist and limit its
further development.

On the other hand, some researches try to use convolutional neural network
(CNN) [26, 27] to deal with sequence tasks, because CNN is easy to run in
parallel. Temporal convolutional network (TCN), which mainly refers to the
Dilated TCN proposed in [28], is a paradigm of using CNN for sequence
modeling. Compared with classic CNN, temporal convolutional network has two
distinctive features:

  1. 1)

By causal convolution [29], each neuron can only receive the input sequence
before its corresponding word, which means there is no possibility of
information ”leakage” from the future to the past.

  2. 2)

It can receive input sequence of any length and map them to output sequence of
the same length, which is consistent with RNN.

While classic TCN is proposed to process video data, new researches have found
that TCN can also provide effective support for language modeling, with
results comparable to mature sequence models such as LSTM, Gate Recurrent Unit
[4, 30, 29, 31]. Among the works, Bai et al. proposes a variant of TCN [30],
which introduces residual structure to TCN, and achieves great performance in
sequence modeling tasks.

##  3 Model

Figure 1: Overview of Carry-lookahead RNN. The model receives sequence data as
input. Embedding layer is optional depending on the type of input sequence.
Carry-lookahead module generates hidden states $h$ in advance. Parallel RNN
module receives both original input $x$ (residual connection marked in red)
and hidden states $h$. A linear layer is added in the end to generate the
output.

In this section, we first discuss the similarities between serial adder and
RNN. Then, inspired by carry-lookahead adder, the carry-lookahead module is
introduced to RNN, which makes it possible for RNN to run in parallel. Also,
we design the method of parallel RNN computation, and finally Carry-lookahead
RNN (CL-RNN) is proposed.

The structure of our proposed model is shown in Fig. 1. The core part of CL-
RNN consists of two modules: Carry-lookahead Module and Parallel RNN Module.
Sequence input is firstly sent to carry-lookahead module to calculate hidden
states in advance. Then parallel RNN module receives both hidden states and
original sequecne input to generate the output. Meanwhile, parallel RNN module
plays a role as residual component for the overall model, which eases the
model optimization difficulty.

###  3.1 Motivation

(a)

(b)

Figure 2: The structure of (a) serial adder and (b) RNN. They are typical
recurrent structure in their respective fields and similar in various aspects.

Serial adder consists of several full adders connected in series, and the
structure is shown in Fig. 2(a). Each full adder $FA_{i}$ waits and receives
the carry bit $C_{i}$ generated by the previous full adder $FA_{i-1}$, and
generates the next carry bit $C_{i+1}$ recursively.

Recurrent neural network unfolds with a similar structure to the serial adder,
as shown in Fig. 2(b). Each neuron receives the hidden state $h_{t-1}$
generated by the previous neuron. Then it produces the result $o_{t}$ and the
new hidden state $h_{t}$. Noting that the internal parameters of each neuron
in RNN are shared, which means the internal function is same. This is
consistent with the fact that each full adder in serial adder has the same
operation function.

Their similar structures lead to the same problem, i.e., serial dependency
hinders parallel computation and lowers efficiency in large-scale
applications. Due to their similarities in various aspects, we consider them
analogous. Since the optimization for adder is a mature topic, the
optimization experience on adder can be applied to improving recurrent neural
network.

###  3.2 Carry-lookahead Module

(a) (b)

(c) (d)

Figure 3: (a) An adder-style edition of Carry-lookahead RNN, which hides the
internal implementation. (b) The structure of carry-lookahead adder. (c) The
structure of Carry-lookahead Module with dilation factor $d$ = 1, 2, 4. (d)
The implementation of the hidden layer in (c).

Carry-lookahead adder is an improved version of serial adder, which is the
mainstream adder nowadays. It consists of a carry-lookahead module and a
number of full adders. The structure is shown in Fig. 3(b), and the basic
principle of carry-lookahead module is shown in Equation (1).

$\begin{cases}G_{i}={A_{i}\cdot B_{i}}\\\ P_{i}={A_{i}+B_{i}}\\\
C_{i+1}={G_{i}+P_{i}\cdot C_{i}}\\\ \end{cases}$ (1)

Taking 4-bit adder as an example, Equation (2) can be derived based on
Equation (1), so that the carry bit of each full adder can be calculated in
advance. Carry-lookahead module breaks the serial dependency between full
adders, realizes the parallelization of multi-bit addition, and greatly
improves the efficiency.

$\begin{cases}C_{1}={G_{0}+P_{0}\cdot C_{0}}\\\ C_{2}={G_{1}+P_{1}\cdot
G_{0}+P_{1}\cdot P_{0}\cdot C_{0}}\\\ C_{3}={G_{2}+P_{2}\cdot G_{1}+P_{2}\cdot
P_{1}\cdot G_{0}+P_{2}\cdot P_{1}\cdot P_{0}\cdot C_{0}}\\\
C_{4}={G_{3}+P_{3}\cdot G_{2}+P_{3}\cdot P_{2}\cdot G_{1}+P_{3}\cdot
P_{2}\cdot P_{1}\cdot G_{0}+P_{3}\cdot P_{2}\cdot P_{1}\cdot P_{0}\cdot
C_{0}}\\\ \end{cases}$

(2)

Through Equation (2), the core idea of carry-lookahead module can be concluded
as to calculate each carry bit with respective preorder input, which can be
expressed as:

$C_{i+1}={f_{i}(C_{0},A_{0},A_{1},\cdots,A_{i},B_{0},B_{1},\cdots,B_{i})}$ (3)

Based on the principle, we introduce carry-lookahead module to RNN (Fig.
3(a)). Carry-lookahead module can be expressed as:

$h_{i+1}=1D\underline{~{}~{}}CNN_{i}(h_{0},x_{0},\cdots,x_{i})$ (4)

which is consistent with Equation (3) in form. It’s implemented with multiple
dilated causal convolutional blocks, and the structure is shown in Fig. 3(c).

Dilated causal convolution combines the techniques of causal convolution and
dilated convolution. They have their own unique purpose. Causal convolution
avoids the information ”leakage” from the future to the past, which enables
carry-lookahead module to process sequential data effectively. Dilated
convolution exponentially enlarges the receptive field. As a result, a very
large receptive field can be achieved with a few layers, which prevents the
model from over-fitting the training data. The receptive field on each layer
can be determined as:

$Receptive\underline{~{}~{}}Field(l)=\frac{k-1}{m-1}\times(m^{l}-1)+1$ (5)

where $l$ is level of layer, $k$ is kernel size, $m$ is multiple of dilation
factor. According to the formula, receptive field can be flexibly adjusted by
modifying dilation factor or the number of layers.

In a word, carry-lookahead module is the key component to achieve parallelism.
It generates hidden states in advance to replace the “recurrent” computation
in RNN. Besides, it has the good characteristics that there is no maximum
length limitation of input, and the output sequence of the module is the same
length as the input sequence, which provides convenience for subsequent
operation.

###  3.3 Parallel RNN Module

Figure 4: Part Left shows the operations of parallel RNN computation. By
carry-lookahead module, hidden states $h$ are calculated in advance. The input
sequence $x$ and hidden states $h$ are split to word vectors. Word vector pair
$(x_{i},h_{i})$ is fed to a _RNNCell_ as input and all _RNNCell_ s run in
parallel. The outputs are merged in the end. Part Right shows the structure of
parallel RNN module. _RNNCell_ s are integrated to one _RNNCell_. The whole
length of sequence is calculated as an entirety.

The classic RNN (Fig. 2(b)) can be viewed as consisting of a series of
parameter-shared _RNNCell_ s. Each _RNNCell_ receives the hidden state $h_{i}$
with the corresponding word vector $x_{i}$ to produce the output $o_{i}$. The
formula is shown as:

$o_{i}=tanh(x_{i}W_{ih}+b_{ih}+h_{i}W_{hh}+b_{hh})$ (6)

where $W_{ih}$, $W_{hh}$, $b_{ih}$ and $b_{hh}$ are internal trainable
parameters.

Limited by serial dependency, RNNs can only process on a word-by-word mode,
because each _RNNCell_ must wait for the completion of the previous one to
obtain the hidden state, then to generate the next one.

In the last section, through carry-lookahead module, hidden states $h$ are
calculated in advance, which makes it possible for RNN to run in parallel. In
this section, we transform the classic RNN into a parallel version and give
the detailed method.

Parallel RNN computation is implemented vis three operations - _Split_ ,
_Parallel computation_ and _Merge_ , as illustrated in Fig. 4 (Part Left). RNN
deals with sequence in word-level, so the input sequence $x$ is split to word
vectors $x_{0},x_{1},\cdots,x_{n}$, and $h$ is also split correspondingly.
Each _RNNCell_ takes the respective $(x_{i},h_{i})$ as input and runs in
parallel. The outputs are merged in the end.

Taking the input sequence ”I love the bird” as an example.

_Input preprocess:_ Each word in the sequence is first encoded to numeric
vector (word vector) by embedding layer [32]. In embedding layer, the word is
converted to an unique value based on dictionary order, then encoded to one-
hot vector and mapped to word vector by linear projection. The internal
parameters in embedding layer can be optimized with stochastic gradient
descent (SGD) to better represent the words. When embedding size is set to
100, each word vector has 100 elements.

$x=Embedding\left(\begin{array}[]{l}I\\\ love\\\ the\\\
bird\end{array}\right)=\left(\begin{array}[]{l}{\rm{0}}{\rm{.014,0}}{\rm{.000,}}\cdots{\rm{,-0}}{\rm{.011}}\\\
{\rm{-0}}{\rm{.008,-0}}{\rm{.000,}}\cdots{\rm{,0}}{\rm{.000}}\\\
{\rm{-0}}{\rm{.006,-0}}{\rm{.002,}}\cdots{\rm{,0}}{\rm{.006}}\\\
{\rm{0}}{\rm{.002,0}}{\rm{.018,}}\cdots{\rm{,-0}}{\rm{.003}}\end{array}\right)$

Here $x$ is formed by the four word vectors and the shape is $(4,100)$.

Then, $x$ is sent to carry-lookahead module (CL Module) to generate the hidden
states in advance.

$h=CL\underline{~{}~{}}Module(h_{0},x)=\left(\begin{array}[]{l}{\rm{-0}}{\rm{.010,0}}{\rm{.009,}}\cdots{\rm{,0}}{\rm{.067}}\\\
{\rm{-0}}{\rm{.024,0}}{\rm{.024,}}\cdots{\rm{,0}}{\rm{.064}}\\\
{\rm{-0}}{\rm{.018,0}}{\rm{.009,}}\cdots{\rm{,0}}{\rm{.084}}\\\
{\rm{-0}}{\rm{.017,0}}{\rm{.028,}}\cdots{\rm{,0}}{\rm{.071}}\end{array}\right)$

_Split:_ The input sequecne $x$ and hidden states $h$ are split to separate
word vectors.

$\left\\{\begin{array}[]{l}x=({x_{0}},{x_{1}},{x_{2}},{x_{3}})\\\
h=({h_{0}}^{\prime},{h_{1}},{h_{2}},{h_{3}})\end{array}\right.$

_Parallel computation:_ Corresponding word vectors form the input pair
$(x_{i},h_{i})$ and sent to a _RNNCell_. There is no dependency between
_RNNCell_ s, so they can calculate in parallel.

$\left\\{\begin{array}[]{l}{o_{0}}=RNNCell({x_{0}},{h_{0}}^{\prime})=({\rm{-0}}.{\rm{210}},{\rm{0}}.{\rm{035}},\cdots,{\rm{0}}.{\rm{036)}}\\\
{o_{1}}=RNNCell({x_{1}},{h_{1}})=({\rm{-0}}.{\rm{064}},{\rm{0}}.{\rm{034}},\cdots,{\rm{0}}.{\rm{044}})\\\
{o_{2}}=RNNCell({x_{2}},{h_{2}})=({\rm{0}}.{\rm{035}},{\rm{-0}}.{\rm{048}},\cdots,{\rm{0}}.{\rm{078}})\\\
{o_{3}}=RNNCell({x_{3}},{h_{3}})=({\rm{0}}.{\rm{036}},{\rm{0}}.{\rm{067}},\cdots,{\rm{-0}}.{\rm{041}})\end{array}\right.$

_Merge:_ Finally, the outputs of _RNNCell_ s
$\\{{o_{0}},{o_{1}},{o_{2}},{o_{3}}\\}$, which are vectors with 100 elements
respectively, are merged in the end. The shape of $o$ is $(4,100)$, which has
the same length as $x$.

$({o_{0}},{o_{1}},{o_{2}},{o_{3}})=o$

_RNNCell_ in the classic RNN, limited by its serial dependency, has to deal
with the input $x$ word by word. While, with the proposed carry-lookahead
module, hidden states $h$ are calculated in advance, so the input $x$ can be
used to calculate as an entirety. Thus, _RNNCell_ in parallel RNN module,
shown as Fig. 4 (Part Right), considers the sequence $x$ and $h$ as
entireties. The formula is expressed as:

$o=RNNCell(x,h)=tanh(xW_{ih}+b_{ih}+hW_{hh}+b_{hh})$ (7)

where $W_{ih}$, $W_{hh}$, $b_{ih}$ and $b_{hh}$ refer to the same things as in
Equation (6).

In Equation (6), which shows the principle of RNNCell in the classic RNN, $x$
is separately dealt as $\\{{x_{0}},{x_{1}},\cdots,{x_{n}}\\}$. In Equation
(7), $x$ and $h$ are considered as entireties, which is further integrated.

Meanwhile, from a global perspective, parallel RNN module plays a role as
residual operation to carry-lookahead module (show as the red line in Fig. 1).
Residual operation is introduced to neural network by He et al. [15] to ease
the model optimization difficulty and it now has been widely applied to image
classification, image dehazing, object detection, etc. [33, 34, 35, 36, 37]
Compared with the classic residual operation, using parallel RNN module as
residual component has following characteristics:

  1. 1)

Compound residual operation: parallel RNN module is a generalization of the
classic residual operation, which adds linear projection and bias to inputs.
When all internal parameters are initialized to identity matrices, it
degenerates to the classic residual operation.

  2. 2)

Trainability of residual operation: the parameters in parallel RNN module are
optimized through iteration. The optimized residual operation makes the model
have stronger representation capability.

  3. 3)

Merging of any dimension: the input sequence will have the same shape as the
hidden state after the linear projection and can be added directly. There is
no need to design a special padding layer or convolution layer for specific
dataset to increase or reduce the dimension, which improves its
generalizability.

Due to the introduction of linear projection, parallel RNN module takes more
parameters than the classic residual operation. However, since CL-RNN is not
large in model size, it is acceptable to appropriately increase the number of
parameters to improve the model performance, which is verified by ablation
experiment in Sec. 4.1.

###  3.4 Evaluation of the Proposed Model

Using CL-RNN for sequence modeling has the following characteristics.

  * 1.

Parallelism. CL-RNN break the serial dependency of RNN by introducing carry-
lookahead module to calculate hidden states in advance, which makes it
possible for RNN to run in parallel. Based on that, parallel RNN module
specifically achieves parallel RNN computation.

  * 2.

Flexible receptive field. The receptive field of neurons in carry-lookahead
module is exponential expanded with the number of layers and can be further
enhanced by dilated convolution. As a result, a very large receptive field can
be achieved with a few layers.

  * 3.

Compound residual operation. Parallel RNN module plays a role as residual
operation in the global model. It improves the representation capability and
eases the model optimization difficulty, which avoids the model being
restricted to local optimal solutions.

##  4 Experiment

In this section, we validate the performance of CL-RNN in sequence modeling on
5 representative datasets. The models are implemented with PyTorch and trained
in Google Colab with NVIDIA Tesla T4 GPU.

The MNIST database is one of the most common databases used for model testing.
Sequential MNIST is one of its varietas, which is specially designed to test
sequence models. To validate the effectiveness of modules in CL-RNN, ablation
experiment is executed on Sequential MNIST first. Then, the model is tested in
character-level language modeling on PTB and text8 datasets, which is a
classic problem in natural language processing. Finally, to test the
generalization ability of the model, polyphonic music modeling is executed on
JSB Chorales and Nottingham datasets, where CL-RNN shows leading advantage.

If not specified, all models are trained through 100 epochs with learning rate
decay [38]. Learning rate decay is a dynamic learning rate strategy, which can
help model to converge better. When the test loss is greater than the maximum
of the last three times, the learning rate will be divided by 10.

###  4.1 Ablation Experiment on Sequential MNIST

Table 1: The hyperparameters on Sequential MNIST. Hyperparameter | Value
---|---
Batch size | 64
Epochs | 12
Initial learning rate | 2e-3
Kernel size | 7
Number of layers | 8
Number of kernels (per layer) | 1
Optimizer | Adam
Dropout | 0.05
Table 2: Evaluation of models on Sequential MNIST on test set. Model | Top-1 Accuracy | Model Size ($\approx$)
---|---|---
CL-RNN | 98.02% | 1M
TCN [28] | 11.35% | 8K
TCN variant [30] | 94.43% | 10K
RNN [2] | 92.83% | 8K
LSTM [3] | 93.68% | 8K

The MNIST database [39] is a large database of handwritten digits, which is
the most common database used in image classification. Sequential MNIST is
based on the MNIST, but each image (size: $28\times 28$) is represented as a
$784\times 1$ sequence. It is adopted by a great deal of researches to test
the ability of recurrent neural networks to retain past information [40, 41,
42, 43, 44, 45].

All the neural network models are not pre-trained, and the hyperparameters are
set as shown in Table 1, where epochs are set to 12, due to the overfit of
models when epochs exceed 12.

It is worth noting that in order to show the differences between models,
control variable method, instead of grid search method, is adopted for
hyperparameters. Otherwise, if the model size is violently enlarged, such as
setting the number of convolution kernels in the hidden layer to 25, all CL-
RNNs are able to reach 99% accuracy, which makes models lose their comparative
meaning. So, we adopt control variable method to highlight the impact of the
adjustment to the models.

The results of ablation experiment on Sequential MNIST are shown in Table 3
and the architectures of the involved models are shown in Table 4. Through the
comparison between CL-RNN(1) shortcut and CL-RNN(1), it can be seen that the
performance of CL-RNN is greatly improved after applying parallel RNN module
as residual component while the model size is similar, which proves the
effectiveness of parallel RNN module.

Table 3: Ablations of modules on Sequential MNIST. CL-RNN(1) refers to the linear layer (at the end of the model) uses only the last item of parallel RNN module’s output sequence; while CL-RNN(784) uses the whole output sequence instead; CL-RNN(1) shortcut refers to the network that parallel RNN module is replaced by addition operation, which is adopted in classic residual operation; addCNN/addRNN refers to CL-RNN with 1D-CNN or RNN layer added behide parallel RNN module. Model | Top-1 Accuracy | Model Size ($\approx$)
---|---|---
CL-RNN(1) shortcut | 92.58% | 8K
CL-RNN(1) | 94.72% | 8K
CL-RNN(784) | 98.02% | 1M
CL-RNN(784) ReLU | 97.71%) | 1M
CL-RNN(784) Dropout(0.2) | 98.05% | 1M
CL-RNN(784) ReLU+Dropout | 97.72% | 1M
CL-RNN(784) ReLU+Dropout addCNN | 96.51% | 1M
CL-RNN(784) addRNN | 97.87% | 1M
Table 4: Architectures of the models for Sequential MNIST. Inside the brackets
are the general shape of dilated causal convolutional blocks, including filter
sizes and feature dimensionalities. $d$ denotes dilation factor, $p$ denotes
padding, $s$ denotes stride, and $i$, $h$ are the feature dimensionalities of
input and output respectively.

Output size | CL-RNN(1) | CL-RNN(784) | CL-RNN(784) addCNN | CL-RNN(784) addRNN
---|---|---|---|---
$784\times 1$ | $\left[\begin{array}[]{l}7\times 7,1,d=1,p=6\\\ 7\times 7,1,d=2,p=12\\\ 7\times 7,1,d=4,p=24\\\ 7\times 7,1,d=8,p=48\\\ 7\times 7,1,d=16,p=96\\\ 7\times 7,1,d=32,p=192\\\ 7\times 7,1,d=64,p=384\\\ 7\times 7,1,d=128,p=768\end{array}\right]$
$784\times 1$ | $RNNCell(i=1,h=1)$ | $RNNCell(i=784,h=784)$
$784(23)\times 1$ | Identity | $\left[\begin{array}[]{l}5\times 5,1,s=5,p=1\\\ 7\times 7,1,s=7,p=2\end{array}\right]$ | $RNN(i=1,h=1)$
$10\times 1$ | 784-d $fc$, log softmax | 23-d $fc$, log softmax | 784-d $fc$, log softmax

Figure 5: Comparison of model test accuracy on Sequential MNIST in the 10
classes. Class size denotes the number of items in this class on test set.
Total size of the 10 classes is 10,000.

(a) (b)

Figure 6: Comparison of model bpc on (a) PTB and (b) text8 on test set.

The last item of output sequence is generated by the neuron with the widest
receptive field. So CL-RNN(1) uses only the last item can achieve such good
results. While CL-RNN(784) uses the entire output sequence for classification,
which means the result is calculated based on multiple scales of input
information. Although it brings much more parameters to linear layer making
the model size bigger, the accuracy further increases to 98%, which is
significantly improved. To balance the performance and efficiency, CL-RNN(1)
is the comprehensive optimal choice; to aim at the higher accuracy, CL-
RNN(784) is the best choice.

Fig. 5 shows the distribution of the test set and the test accuracy of the
models in the 10 classes. All models perform well in the classes of digit with
simple stroke, such as 0 and 1. But when the stroke are relatively complex,
such as 2 and 8, the accuracy fluctuates considerably, which can well reflect
the representation capability of the models. CRN(784) not only perform best in
each class, but also has the smoothest curve, which indicates its
representation capability is sufficient for the complex digits.

Some other modules are also tried to be added on the model. Since the output
sequence of parallel RNN module is as same length as the input sequence. It’s
convenient to add 1D-CNN or RNN layer behind it. The results show no
significant improvement in accuracy, which indicates that CL-RNN(1) and CL-
RNN(784) are suitable enough on Sequential MNIST. Therefore, in the subsequent
experiments, CL-RNN(784) is used as the main model.

In comparison experiment, CL-RNN is compared with TCN, RNN, LSTM. The results
are shown in Table 2. Considering the models are trained from scratch on the
same dataset, the results can objectively verify the performance and
feasibility of CL-RNN. According to Table 2, CL-RNN performs much better than
the classic sequence models.

Besides, it’s found that TCN performs quite bad in the tasks. We suppose the
reason is that the network can’t be well initialized in this task, which makes
the model restricted to a local optimal solution and can’t get further
optimized. Although several methods (e.g. change the initialization methods,
change the random seed, change the initial learning rate) have been tried, the
results remain similar. Compared with that, Parallel RNN Module in our
proposed CL-RNN plays a role as the residual component in the network, which
eases the optimization difficulty and avoids such problem. As a result, the
performance gets significant improvement.

As TCN is unable to normally perform, CL-RNN is also compared with the variant
of TCN. According to the results, CL-RNN has better performance, which further
verifies the effectiveness of the proposed model.

###  4.2 Character-level Language Modeling

Table 5: The hyperparameters on character-level language modeling. Hyperparameter | Value
---|---
Batch size | 32
Epochs | 100
Kernel size | 3
Number of layers | 3
Number of kernels (per layer) | 450
Optimizer | SGD
Dropout | 0.1
Gradient Clip | 0.15
Initial Learning rate | 4
Embedding size | 100
Valid sequence length | 320
Sequence length | 400
Table 6: Evaluation of models on character-level language modeling on test set. CL-RNN LSTM refers to CL-RNN with LSTM added behind parallel RNN module. Model | bpc (PTB) | bpc (text8) | Model Size ($\approx$)
---|---|---|---
CL-RNN | 1.394 | 1.769 | 0.9M
CL-RNN LSTM | 1.382 | 1.741 | 1M
TCN [28] | 1.391 | 1.771 | 0.9M
TCN variant [30] | 1.376 | 1.730 | 3M
LSTM [3] | 1.587 | 1.975 | 85K

(a) (b)

Figure 7: Comparison of model loss on (a) JSB and (b) Nott on test set.

Character-level language modeling treats each character in text as a word,
which greatly compresses the encoding space. Compared with word-level language
modeling, it ignores the semantics of words in the text but requires the model
to capture the semantic information. The idea of character-level language
modeling comes from signal processing and it’s a classic problem in natural
language processing. This paper makes a comparison test of character-level
language modeling on two datasets.

  * 1.

PennTreebank (PTB)
PennTreebank [46] consists of 49 characters and contains 5059K characters for
training, 396K characters for validation, and 446K characters for test. It’s a
relatively small but highly studied dataset [47, 44, 48].

  * 2.

text8
text8 [49] consists of 27 characters and has about 20 times the amount of data
as PennTreebank. It contains 100M characters from Wikipedia, 90M of which are
used for training, 5M for validation, and 5M for test.

In implementation, _Sequence length_ characters are intercepted from the
dataset as a sentence. For forward propagation, the whole sentence is used as
the input sequence. For backward propagation, only the last _Valid sequence
length_ characters of the sentence are used, which ensures that the training
data are also provided with sufficient history.

We use cross-entropy as the loss function and bits-per-character (bpc) [50] to
measure the performance of the language model, which is an usual measure index
in language modelling. Bpc is calculated with the following formula, and it
can be seen that the lower bpc is better.

$bpc=\frac{{\overline{loss}}}{{\ln(2)}}$ (8)

The hyperparameters of the two tests are set as shown in Table 5, and the test
results are shown in Table 6. Fig. 6(a) and Fig. 6(b) show the loss of model
versus iteration on the two datasets.

According to Table 6, it’s found that the performance of LSTM is much worse
than other models, probably due to the small number of parameters which limits
the representation capability of the model. Besides, the performance of TCN
and CL-RNN is similar.

Although the performance of CL-RNN and LSTM itself is mediocre, CL-RNN LSTM,
which is formed by the superposition of the two, is much better. It achieves
the lowest bpc in both two tests. Moreover, Fig. 6(a) and Fig. 6(b) shows that
CL-RNN LSTM is not fully converged after 100 epochs of training, which
indicates that there is still room for further optimization.

With three times the size than other models, TCN variant achieves better
performance in character-level language modeling. However, big model takes
long time to train. Compared with that, CL-RNN LSTM well balances the
relationship between the model size and performance.

###  4.3 Polyphonic Music Modeling

Table 7: The hyperparameters on polyphonic music modeling. Hyperparameter | Value
---|---
Epochs | 100
Kernel size | 5
Number of layers | 4
Number of kernels (per layer) | 150
Optimizer | Adam
Dropout | 0.25
Gradient Clip | 0.2
Initial Learning rate (JSB) | 1e-2
Initial Learning rate (Nott) | 1e-3
Table 8: Evaluation of models on polyphonic music modeling on test set.

Model | Loss (JSB) | Loss (Nott) | Model Size ($\approx$)
---|---|---|---
CL-RNN | 8.565 | 2.776 | 455K
TCN [28] | 9.158 | 3.305 | 418K
TCN variant [30] | 8.589 | 2.969 | 883K
RNN [2] | 8.707 | 3.506 | 50K
LSTM [3] | 8.442 | 3.567 | 157K

A piece of music consists of a sequence of notes, and many researches have
applied recurrent network architecture to music modeling tasks [4, 51, 52,
53]. In this paper, CL-RNN is tested on two polyphonic music datasets.

  * 1.

JSB Chorales (JSB)
JSB Chorale [54] is a polyphonic music dataset consisting of 382 four-part
harmonic works composed by J. S. Bach. Each element in the input sequence
consists of 88 bits, corresponding to the 88 keys on the piano, and value of 1
indicates that the key is pressed at that moment.

  * 2.

Nottingham (Nott)
Nottingham is a dataset consisting of 1200 English and American folk tunes,
which is a larger dataset compared to JSB Chorales.

The hyperparameters of the two tests are set as shown in Table 7, except that
the initial learning rates for TCN on Nott dataset and Residual TCN on JSB and
Nott datasets are set to 1e-4 (gradient explosion occurs at 1e-3). Cross-
entropy is used as the loss function and the value of loss on test set is used
to measure the performance of the models. We compare CL-RNN with TCN, RNN, and
LSTM, and the results are shown in Table 8. Fig. 7(a) and Fig. 7(b) show the
loss of models versus iteration.

According to Table 8, CL-RNN performs well in both tests. Although it is just
slightly inferior to LSTM on JSB dataset, it is far ahead compared to other
models on Nott dataset, which shows the strong representation and
generalization capability.

Furthermore, as shown in Fig. 7(a) and especially Fig. 7(b), CL-RNN converges
significantly faster than the other models and achieves good result at the
early stage of training, which indicates the model is well and easily
optimized on this dataset.

##  5 Conclusion and Future Work

Considering the analogy between serial adder and RNN, this paper puts forward
the idea of applying mature optimization experience on adders to improving
RNN. Carry-lookahead adder is the mainstream adder at present. Based on its
core idea, CL-RNN is proposed, which consists of carry-lookahead module and
parallel RNN module. It overcomes the problem of serial dependency of
recurrent architecture and has the advantages of parallelism, flexible
receptive field, compound residual operation, etc. The proposed CL-RNN
specifically targets a comprehensive set of sequence modeling tasks that have
been repeatedly used to compare the effectiveness of different recurrent
neural networks. Experimental results prove that CL-RNN can perform better
than classic TCN, RNN and LSTM in this domain.

The core idea of this paper is the transfer of adder optimization experience
and the proposal of CL-RNN, but deeper stack of CL-RNN is not explored.
Referring to the architecture of ResNet [15], the core part of CL-RNN can be
viewed as an integrated RNN residual block to build deeper networks in the
future. Besides, more variants of TCN have been proposed as TCN gradually gets
popular in sequence modeling tasks. Meanwhile, Carry-lookahead RNN based on
various variants of TCN can be further explored.

## Acknowledgement

This work was supported by National Natural Science Foundation of China (Nos.
61972121, 61971173, U1909210), Zhejiang Provincial Natural Science Foundation
of China (Nos. LY21F020015, LY20F020015, LY21F030005), Fundamental Research
Funds for the Provincial Universities of Zhejiang (No. GK209907299001-008),
National Social Science Foundation of China (No. 19ZDA348). The authors would
like to thank the reviewers for their comments and suggestions in advance.

## References

  * [1] I. Goodfellow, Y. Bengio, A. Courville, Y. Bengio, Deep learning, Vol. 1, MIT Press Cambridge, 2016.
  * [2] D. E. Rumelhart, G. E. Hinton, R. J. Williams, Learning representations by back-propagating errors, Nature 323 (6088) (1986) 533–536.
  * [3] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural Computation 9 (8) (1997) 1735–1780.
  * [4] J. Chung, C. Gulcehre, K. Cho, Y. Bengio, Empirical evaluation of gated recurrent neural networks on sequence modeling, arXiv preprint arXiv:1412.3555 (2014).
  * [5] I. Sutskever, O. Vinyals, Q. V. Le, Sequence to sequence learning with neural networks, arXiv preprint arXiv:1409.3215 (2014).
  * [6] D. Bahdanau, K. Cho, Y. Bengio, Neural machine translation by jointly learning to align and translate, arXiv preprint arXiv:1409.0473 (2014).
  * [7] K. Cho, B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, Y. Bengio, Learning phrase representations using RNN encoder-decoder for statistical machine translation, arXiv preprint arXiv:1406.1078 (2014).
  * [8] Y. Wu, M. Schuster, Z. Chen, Q. V. Le, M. Norouzi, W. Macherey, M. Krikun, Y. Cao, Q. Gao, K. Macherey, et al., Google’s neural machine translation system: Bridging the gap between human and machine translation, arXiv preprint arXiv:1609.08144 (2016).
  * [9] R. Jozefowicz, O. Vinyals, M. Schuster, N. Shazeer, Y. Wu, Exploring the limits of language modeling, arXiv preprint arXiv:1602.02410 (2016).
  * [10] A. F. Kirichenko, I. V. Vernik, M. Y. Kamkar, J. Walter, M. Miller, L. R. Albu, O. A. Mukhanov, ERSFQ 8-bit parallel arithmetic logic unit, IEEE Transactions on Applied Superconductivity 29 (5) (2019) 1–7.
  * [11] H. Thapliyal, M. Srinivas, A novel reversible TSG gate and its application for designing reversible carry look-ahead and other adder architectures, in: Asia-Pacific Conference on Advances in Computer Systems Architecture, Springer, 2005, pp. 805–817.
  * [12] Y. Zhou, F. He, N. Hou, Y. Qiu, Parallel ant colony optimization on multi-core SIMD CPUs, Future Generation Computer Systems 79 (2018) 473–487.
  * [13] N. Hou, X. Yan, F. He, A survey on partitioning models, solution algorithms and algorithm parallelization for hardware/software co-design, Design Automation for Embedded Systems 23 (1) (2019) 57–77.
  * [14] Y. Zhou, F. He, Y. Qiu, Accelerating image convolution filtering algorithms on integrated CPU–GPU architectures, Journal of Electronic Imaging 27 (3) (2018) 033002.
  * [15] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2016, pp. 770–778.
  * [16] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, A. Rabinovich, Going deeper with convolutions, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 1–9.
  * [17] G. B. Rosenberger, Simultaneous carry adder, US Patent 2,966,305 (Dec. 27 1960).
  * [18] D. Krueger, R. Memisevic, Regularizing RNNs by stabilizing activations, arXiv preprint arXiv:1511.08400 (2015).
  * [19] G. E. Hinton, Deep belief networks, Scholarpedia 4 (5) (2009) 5947.
  * [20] R. Vohra, K. Goel, J. K. Sahoo, Modeling temporal dependencies in data using a DBN-LSTM, in: IEEE International Conference on Data Science and Advanced Analytics, IEEE, 2015, pp. 1–4.
  * [21] X. Shi, Z. Chen, H. Wang, D. Y. Yeung, W. K. Wong, W. C. Woo, Convolutional LSTM network: A machine learning approach for precipitation nowcasting, Advances in Neural Information Processing Systems 2015 (2015) 802–810.
  * [22] J. Bradbury, S. Merity, C. Xiong, R. Socher, Quasi-recurrent neural networks, arXiv preprint arXiv:1611.01576 (2016).
  * [23] S. Chang, Y. Zhang, W. Han, M. Yu, X. Guo, W. Tan, X. Cui, M. Witbrock, M. A. Hasegawa-Johnson, T. S. Huang, Dilated recurrent neural networks, Advances in Neural Information Processing Systems 30 (2017) 77–87.
  * [24] O. Kuchaiev, B. Ginsburg, Factorization tricks for LSTM networks, arXiv preprint arXiv:1703.10722 (2017).
  * [25] N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, J. Dean, Outrageously large neural networks: The sparsely-gated mixture-of-experts layer, arXiv preprint arXiv:1701.06538 (2017).
  * [26] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, L. D. Jackel, Backpropagation applied to handwritten zip code recognition, Neural Computation 1 (4) (1989) 541–551.
  * [27] A. Krizhevsky, I. Sutskever, G. E. Hinton, ImageNet classification with deep convolutional neural networks, Advances in Neural Information Processing Systems 25 (2012) 1097–1105.
  * [28] C. Lea, M. D. Flynn, R. Vidal, A. Reiter, G. D. Hager, Temporal convolutional networks for action segmentation and detection, in: Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2017, pp. 156–165.
  * [29] A. v. d. Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, K. Kavukcuoglu, WaveNet: A generative model for raw audio, arXiv preprint arXiv:1609.03499 (2016).
  * [30] S. Bai, J. Z. Kolter, V. Koltun, An empirical evaluation of generic convolutional and recurrent networks for sequence modeling, arXiv preprint arXiv:1803.01271 (2018).
  * [31] N. Kalchbrenner, L. Espeholt, K. Simonyan, A. v. d. Oord, A. Graves, K. Kavukcuoglu, Neural machine translation in linear time, arXiv preprint arXiv:1610.10099 (2016).
  * [32] Q. Le, T. Mikolov, Distributed representations of sentences and documents, in: International Conference on Machine Learning, PMLR, 2014, pp. 1188–1196.
  * [33] F. Wang, M. Jiang, C. Qian, S. Yang, C. Li, H. Zhang, X. Wang, X. Tang, Residual attention network for image classification, in: Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2017, pp. 3156–3164.
  * [34] Z. Zhong, J. Li, Z. Luo, M. Chapman, Spectral–spatial residual network for hyperspectral image classification: A 3-D deep learning framework, IEEE Transactions on Geoscience and Remote Sensing 56 (2) (2017) 847–858.
  * [35] S. Zhang, F. He, DRCDN: learning deep residual convolutional dehazing networks, The Visual Computer 36 (9) (2020) 1797–1808.
  * [36] S. Chen, X. Tan, B. Wang, H. Lu, X. Hu, Y. Fu, Reverse attention-based residual network for salient object detection, IEEE Transactions on Image Processing 29 (2020) 3763–3776.
  * [37] D. Zhang, L. He, M. Luo, Z. Xu, F. He, Weight asynchronous update: Improving the diversity of filters in a deep convolutional network, Computational Visual Media 6 (4) (2020) 455–466.
  * [38] Y. LeCun, I. Kanter, S. Solla, Second order properties of error surfaces: Learning time and generalization, in: R. P. Lippmann, J. Moody, D. Touretzky (Eds.), Advances in Neural Information Processing Systems, Vol. 3, Morgan-Kaufmann, 1991, pp. 918–924.
  * [39] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, Gradient-based learning applied to document recognition, Proceedings of the IEEE 86 (11) (1998) 2278–2324.
  * [40] Q. V. Le, N. Jaitly, G. E. Hinton, A simple way to initialize recurrent networks of rectified linear units, arXiv preprint arXiv:1504.00941 (2015).
  * [41] S. Zhang, Y. Wu, T. Che, Z. Lin, R. Memisevic, R. Salakhutdinov, Y. Bengio, Architectural complexity measures of recurrent neural networks, arXiv preprint arXiv:1602.08210 (2016).
  * [42] S. Wisdom, T. Powers, J. R. Hershey, J. L. Roux, L. Atlas, Full-capacity unitary recurrent neural networks, arXiv preprint arXiv:1611.00035 (2016).
  * [43] N. Kalchbrenner, E. Grefenstette, P. Blunsom, A convolutional neural network for modelling sentences, arXiv preprint arXiv:1404.2188 (2014).
  * [44] D. Krueger, T. Maharaj, J. Kramár, M. Pezeshki, N. Ballas, N. R. Ke, A. Goyal, Y. Bengio, A. Courville, C. Pal, Zoneout: Regularizing RNNs by randomly preserving hidden activations, arXiv preprint arXiv:1606.01305 (2016).
  * [45] L. Jing, Y. Shen, T. Dubcek, J. Peurifoy, S. Skirlo, Y. LeCun, M. Tegmark, M. Soljačić, Tunable efficient unitary neural networks (EUNN) and their application to RNNs, in: International Conference on Machine Learning, PMLR, 2017, pp. 1733–1741.
  * [46] M. Marcus, B. Santorini, M. A. Marcinkiewicz, Building a large annotated corpus of English: The Penn Treebank, Computational Linguistics 19 (2) (1993) 313–330.
  * [47] Y. Miyamoto, K. Cho, Gated word-character recurrent language model, arXiv preprint arXiv:1606.01700 (2016).
  * [48] S. Merity, N. S. Keskar, R. Socher, Regularizing and optimizing LSTM language models, arXiv preprint arXiv:1708.02182 (2017).
  * [49] T. Mikolov, I. Sutskever, A. Deoras, H.-S. Le, S. Kombrink, J. Cernocky, Subword language modeling with neural networks, preprint (http://www.fit.vutbr.cz/imikolov/rnnlm/char.pdf) 8 (2012) 67.
  * [50] A. Graves, Generating sequences with recurrent neural networks, arXiv preprint arXiv:1308.0850 (2013).
  * [51] R. Pascanu, C. Gulcehre, K. Cho, Y. Bengio, How to construct deep recurrent neural networks, arXiv preprint arXiv:1312.6026 (2013).
  * [52] R. Jozefowicz, W. Zaremba, I. Sutskever, An empirical exploration of recurrent network architectures, in: International Conference on Machine Learning, PMLR, 2015, pp. 2342–2350.
  * [53] K. Greff, R. K. Srivastava, J. Koutník, B. R. Steunebrink, J. Schmidhuber, LSTM: A search space odyssey, IEEE Transactions on Neural Networks and Learning Systems 28 (10) (2016) 2222–2232.
  * [54] M. Allan, C. K. Williams, Harmonising chorales by probabilistic inference, Advances in Neural Information Processing Systems 17 (2005) 25–32.