figurec

# A numerical study of gravity driven instability in strongly coupled dusty
plasmas. Part I: Rayleigh-Taylor instability and Buoyancy-driven instability

Vikram S. Dharodi<EMAIL_ADDRESS>Mechanical Engineering, Michigan State
University, East Lansing, Michigan 48824, USA.  Amita Das<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology, New Delhi, 110016,
Delhi, India

###### Abstract

Rayleigh-Taylor and Buoyancy-driven instabilities are very common
instabilities for an inhomogeneous medium. We examine here how these
instabilities grow for incompressible viscoelastic fluids like a strongly
coupled dusty plasma by using incompressible generalized hydrodynamic (i-GHD)
fluid model. Since the dust particles are pretty massive, the gravitational
attraction of earth has a significant role in its dynamics. In this paper the
acceleration due to gravity g and the role of strong coupling in the context
of gravitationally stratified dusty plasma fluid have been considered. We find
that the appearance of elasticity speed up the growth of viscoelastic RT
instability by reducing the effect of viscosity at timescales shorter than the
Maxwell relaxation time $\tau_{m}$. The buoyancy driven situation with
spatially localized (in both dimensions) of low/high density regions placed in
a background medium has also been evolved in the presence of gravity. Our
studies show that both the instabilities get suppressed with increasing
coupling strength of the medium. This suppression has been illustrated
analytically as well as by carrying out a two-dimensional fully nonlinear
simulation.

##  I Introduction

The Rayleigh-Taylor (RT) instability is one of the most prominent fluid
instabilities. This instability, for a fluid in a gravitational field, was
first observed by Rayleigh Rayleigh (1900) and later implemented to all
accelerated fluids by Taylor Taylor (1950). Since then numerous experiments
and numerical simulations have been conducted in order to better understand
it. For a hydrodynamic (HD) fluid under gravity, a detailed theoretical work
has been done by S. Chandrasekhar Chandrasekhar . A lot of research studies
has been pursued and published on this topic so far Youngs (1984); Waddell,
Niederhaus, and Jacobs (2001); Livescu (2004); Guo and Tice (2010); Baldwin,
Scase, and Hill (2015); Lyubimova, Vorobev, and Prokopev (2019). Furthermore,
in the context of plasma medium the charged species respond to electromagnetic
forces and various scenarios arise where an inverse stratification against a
force (and/or a pseudo force due to the choice of the frame) leads to this
instability. Studies for different types of plasmas under various conditions
have been carried out Dickinson  _et al._ (1962); Ariel and Bhatia (1970);
Takabe  _et al._ (1985); Mikhailenko, Mikhailenko, and Weiland (2002); Ma  _et
al._ (2006); Sen, Fukuyama, and Honary (2010); Weber  _et al._ (2014);
Kessler, Dahlburg, and Ganguli (2014); Haines  _et al._ (2014); Hoshoudy
(2014); Khiar  _et al._ (2019). This instability has been observed in diverse
situations such as a supernova explosion Norman  _et al._ (1981); Arnett  _et
al._ (1989), geophysics Plag and Jüttner (1995), astrophysics Allen and Hughes
(1984); Zingale  _et al._ (2005), atmospheric physicsKeskinen  _et al._
(1981), liquid atomization Beale and Reitz (1999); Kong, Senecal, and Reitz
(1999), fusion physics Chertkov (2003); Kadau  _et al._ (2004); Hoshoudy
(2011); Srinivasan and Tang (2013), oceanography Debnath (1994), turbulent
mixing Abarzhi (2010, 2008); Abarzhi and Rosner (2010), etc. Another form of
RT instability is the buoyancy-driven (BD) instability, here a buoyant force
acts in the direction opposite to gravity for low density regions. The BD
instability decides whether an object will float (if the density of an object
is less than background fluid) or sink (if the density of an object is greater
than background fluid) in fluid. The floating of boats and ships while sinking
of small objects like rocks in water, or the pouring cream (heavy fluid) into
coffee (light fluid) and petroleum wells are some examples where buoyant force
plays a major role.

Our objective here is to understand the both RT and BD instabilities when a
medium is in strongly coupled state. For this purpose we consider the case of
a dusty plasma. A dusty plasma can be easily prepared/found in the strong
coupling state and also the gravity plays an important role for a micron-sized
dust particle has a mass of about $10^{-15}-10^{-10}$Kg. Thus the role of
strong coupling in behaviour of these instabilities can be readily
investigated in the context of strongly coupled dusty plasmas (SCDPs). The
SCDPs are modelled, here, by using the incompressible generalized hydrodynamic
(i-GHD) fluid model. This model captures the strong coupling behaviour in
terms of a viscoelastic fluid characteristics through two crucial parameters,
the shear viscosity $\eta$ and elasticity $\tau_{m}$ (representing memory
relaxation time) Frenkel (1955). For those phenomena which are faster compared
to $\tau_{m}$, the GHD system retains the memory of the past configurations
and the elasticity effects dominate. This elastic behaviour of the system is
known to produce transverse shear (TS) waves in the medium. However, for times
longer than $\tau_{m}$ the memory effects are insignificant and the usual
viscous characteristics of the fluid phase dominate.

Here, we demonstrate that elasticity leads to a speed up of the RT instability
by reducing the coupling strength ($\eta$/$\tau_{m}$, viscous fraction). Such
reducing effect of elasticity on RT instability has also been observed for
geodynamic settings in Kaus and Becker (2007) and for the viscoelastic fluids
within the framework of the Oldroyd-B model in Guido  _et al._ (2010). In
order to evident this effect, for the same value of $\eta$, we have simulated
the fully viscous fluid and viscoelastic fluids for different values of
$\tau_{m}$. To develop a better physical insight into the dynamics of each
phenomena, the inviscid limit of fluids are also simulated. The BD instability
has been discussed by considering the evolution of localized low/high density
bubbles/droplets in the medium. For both instabilities, the simulations show a
gradual suppression with increasing coupling strength, which is consistent
with the linear theory. The coupling strength has been understood in terms of
shear wave whose phase velocity is proportional to it.

The paper has been organized as follows. Section II contains the governing
equations of i-GHD model for viscoelastic medium in the presence of gravity.
In Subsections II.1 and II.2, we obtain the analytical dispersion relations
for the gradually and sharply increasing density gradients against
gravitational force, respectively. Section III contains the description of
numerical scheme. In Section IV, we describe the numerical evolution of
various density profiles with time for different values of coupling
parameters. The suppression of the both instabilities i.e. RT and BD has been
clearly depicted by numerical simulation, as one moves from weakly coupled to
strongly coupled regime. In Section VI, we describe the results and
observations.

##  II Analytical Description

Generalized hydrodynamic (GHD) fluid model is a phenomenological model which
is used to study the strongly coupled dusty plasma system Kaw and Sen (1998);
Kaw (2001); Tiwari  _et al._ (2012); S. Dharodi, K. Tiwari, and Das (2014);
Tiwari  _et al._ (2014); Dharodi  _et al._ (2016). GHD model treats dusty
plasma as a viscoelastic fluid in which the coupling strength depends on the
ratio of the shear viscosity coefficient $\eta$ and relaxation time parameter
$\tau_{m}$ Frenkel (1955). This model supports both the existence of
incompressible transverse shear and compressible longitudinal modes. Here, we
consider only the incompressible limit of the GHD and term is as the i-GHD
fluid model for which the set of equations is obtained. In the incompressible
limit the Poisson equation is replaced by the quasi-neutrality condition and
charge density fluctuations are ignored. The coupled set of continuity and
momentum equations for the dust fluid under gravity acceleration $\vec{g}$ can
be written as:

$\frac{\partial\rho_{d}}{\partial
t}+\nabla\cdot\left(\rho_{d}\vec{v}_{d}\right)=0{,}$ (1)
$\displaystyle\left[1+{\tau_{m}}\left(\frac{\partial}{\partial{t}}+{\vec{v}_{d}}\cdot\nabla\right)\right]$
$\displaystyle\left[{{\rho_{d}}\left(\frac{\partial{\vec{v}_{d}}}{\partial{t}}+{\vec{v}_{d}}{\cdot}\nabla{\vec{v}_{d}}\right)}+{\rho_{d}}\vec{g}+{\rho_{c}}\nabla\phi_{d}\right]$
$\displaystyle=\eta\nabla^{2}\vec{v}_{d}{,}$ (2)

respectively and the incompressible condition is given as

${\nabla}{\cdot}{\vec{v_{d}}}=0{.}$ (3)

The derivation of these reduced equations has been discussed in detail in our
earlier papers S. Dharodi, K. Tiwari, and Das (2014); Dharodi  _et al._ (2016)
along with the procedure of its numerical implementation and validation. Here,
the variables $\rho_{d}$ and $\rho_{c}$, $\vec{v}_{d}$ and $\phi_{d}$ are the
dust mass density, dust charge density, dust fluid velocity and dust charge
potential respectively. The time, length, velocity and potential are
normalised by inverse of dust plasma frequency
$\omega^{-1}_{pd}=\left({4\pi(Z_{d}e)^{2}n_{d0}}/{m_{d0}}\right)^{-1/2}$ and
plasma Debye length
$\lambda_{d}=\left({K_{B}T_{i}}/{4{\pi}{Z_{d}}{n_{d0}}{e^{2}}}\right)^{1/2}$,
${\lambda_{d}}{\omega_{pd}}$ and ${{Z_{d}}e}/{{K_{B}}{T_{i}}}$ respectively.
The parameters $m_{d}$, $T_{i}$ and $K_{B}$ are the dust grain mass, ion
temperature and Boltzmann constant respectively. $Z_{d}$ is the charge on each
dust grain with no consideration of charge fluctuation. The number density
$n_{d}$ is normalised by the equilibrium value $n_{d0}$.

###  II.1 Gradual density gradient

Here, we consider two-dimensional (x-y plane) incompressible system where
density gradient and potential gradient are chosen along the y-axis i.e.
${{\partial{\rho_{d}}}/{\partial y}}$, ${{\partial\phi_{d0}}/{\partial y}}$
respectively, and acceleration $\vec{g}$ applied in opposite direction of
fluid density gradient $-g\hat{y}$. Initially, we consider no initial flow
i.e. $\vec{v}_{d0}=0$ at t=0, the equilibrium condition Eq. (II) becomes

${\rho_{d0}}{g}=-{\rho_{c}}{\frac{\partial\phi_{d0}}{\partial y}}{.}$ (4)

A small perturbation in the various fields, e.g. density, scalar potential and
dust velocity can be written as

${\rho_{d}}(x,y,t)={\rho_{d0}}(y,t=0)+{\rho_{d1}}(x,y,t){,}$ (5)
${\phi_{d}}(x,y,t)={\phi_{d0}}(y,t=0)+{{\phi_{d1}}}(x,y,t){,}$ (6)
${\vec{v}_{d}}(x,y,t)=0+{\vec{v}_{d1}}(x,y,t){,}$ (7)

respectively. Retaining only linear terms in the perturbed fields we obtain
the following equations for the linearized instability analysis

$\frac{\partial\rho_{d1}}{\partial
t}+\left(\vec{v_{d1}}\cdot\nabla\right)\rho_{d0}=0{,}$ (8)
$\left[1+\tau_{m}\frac{\partial}{\partial
t}\right]\left[{{\rho_{d0}}{\frac{\partial{v}_{d1y}}{\partial
t}}}+{\rho_{d1}}g{\hat{y}}+{\rho_{c}}{\frac{\partial\phi_{d1y}}{\partial
y}}\right]=\eta\nabla^{2}{v}_{d1y}{,}$ (9)
$\left[1+\tau_{m}\frac{\partial}{\partial
t}\right]\left[{{\rho_{d0}}{\frac{\partial{v}_{d1x}}{\partial
t}}}+{\rho_{c}}{\frac{\partial\phi_{d1x}}{\partial
x}}\right]=\eta\nabla^{2}{v}_{d1x}{,}$ (10) $\nabla\cdot\vec{v}_{d1}=0{.}$
(11)

Since equilibrium fields vary along $\hat{y}$, the above set of equations can
only be Fourier analysed in time and spatial coordinate $x$. However, assuming
the perturbations to be of much smaller scale compared to the equilibrium
variation along $y$ we invoke the local approximation and proceed with the
Fourier decomposition along y also. This leads to

${\rho_{d1}}=-\frac{i}{{\omega}}\frac{\partial\rho_{d0}}{\partial
y}{v_{1y}}{.}$ (12)

The $\hat{y}$ component

${\left[1-i{\omega}\tau_{m}\right]}{\left[{-i{\omega}{\rho_{d0}}{{v}_{1y}}}+{\rho_{d1}}g+i{\mathrm{k}_{y}}{\rho_{c}}{\phi_{1}}\right]}=-{\eta{\mathrm{k}}^{2}{v}_{1y}}{.}$
(13)

The $\hat{x}$ component

${\left[1-i{\omega}\tau_{m}\right]}{\left[{-i{\omega}{\rho_{d0}}{{v}_{1x}}}+i{\mathrm{k}_{x}}{\rho_{c}}{\phi_{1}}\right]}=-{\eta{\mathrm{k}}^{2}{v}_{1x}}{,}$
(14)
$i{\mathrm{k}_{x}}{v_{1x}}=-i{\mathrm{k}_{y}}{v_{1y}}\Rightarrow{v_{1y}}={-\frac{\mathrm{k}_{x}}{\mathrm{k}_{y}}}{v_{1x}}{.}$
(15)

Obtain $\phi_{d1}$, using above relation and Eqs. (12) (13) (14), and (15), we
get

${{\phi_{d1}}}=-{\frac{{g}{\mathrm{k}_{x}}{v_{1x}}}{{\rho_{c}}{\omega}{\mathrm{k}^{2}}}}\frac{\partial\rho_{d0}}{\partial
y}{,}$ (16)

using the above relation Eq. (16) in Eq. (14), we get the dispersion relation
as

${\left(1-i{\omega}\tau_{m}\right)}{\left[{\omega^{2}}{\mathrm{k}^{2}}+{\frac{g\mathrm{k}^{2}_{x}}{{{\rho_{d0}}}}}\frac{\partial\rho_{d0}}{\partial
y}\right]}={-i{\omega}{\eta^{*}}{\mathrm{k}}^{4}}{,}$ (17)

where kinematic viscosity $\eta^{*}$=${\eta}/{\rho_{d0}}$, the ratio of
absolute viscosity ${\eta}$ to density ${\rho_{d0}}$. A similar equation for
the electron ion plasma in the context of inertial fusion has been predicted
by A. Das et. al. Das and Kaw (2014). When $\tau_{m}$ the memory relaxation
parameter is zero, the RT growth rate shows a viscous damping due to
$\eta^{*}$. In strongly coupled regime i.e. ${\omega}{\tau_{m}}\gg 1$, the
above dispersion Eq. (17) becomes

${\omega^{2}}={\frac{\eta^{*}}{{\tau_{m}}}{\mathrm{k}}^{2}}-{\frac{g}{{{\rho_{d0}}}}}\frac{\partial\rho_{d0}}{\partial
y}{\frac{{\mathrm{k^{2}_{x}}}}{{\mathrm{k}}^{2}}}{.}$ (18)

Here $\sqrt{\eta^{*}/\tau_{m}}=V_{s}$ is the phase speed of the transverse
shear wave. For hydrodynamic fluid case where ($\eta^{*}/\tau_{m}$=0), the
above dispersion relation reduces to the usual expression of the RT growth
rate is

${\omega^{2}}=-{\frac{g}{{{\rho_{d0}}}}}\frac{\partial\rho_{d0}}{\partial
y}{\frac{{\mathrm{k^{2}_{x}}}}{{\mathrm{k}}^{2}}}{.}$ (19)

###  II.2 Sharp interface

When the density interface is sharp the above dispersion relation becomes,

${\omega^{2}}={\frac{\eta^{*}}{{\tau_{m}}}{{\mathrm{k}}^{2}}}-g{\mathrm{k}}{A_{T}}{,}$
(20)

here $A_{T}={(\rho_{d02}-\rho_{d01})}/{(\rho_{d01}+\rho_{d02})}$ is the
Atwood’s number. $\rho_{d02}$ is heavier fluid density is supported by lighter
fluid density $\rho_{d01}$. Here, too the hydrodynamic limit is recovered when
the TS wave velocity $\eta^{*}/\tau_{m}=0$ is chosen.

${\omega^{2}}=-g{\mathrm{k}}{A_{T}}{.}$ (21)

For $g=0$, Eqs. (18) and (20) imply the existence of TS wave moving with phase
velocity ${\sqrt{{\eta^{*}}/{\tau_{m}}}}$. Detailed analytical and numerical
studies of TS waves for the homogeneous S. Dharodi, K. Tiwari, and Das (2014);
Dharodi  _et al._ (2016) and inhomogeneous Dharodi (2020) density media have
been carried out by us. From these two equations one can immediately infer
that in the strong coupling limit as the value of the first term
${{\eta^{*}}/{\tau_{m}}}$ (coupling strength parameter) increases, it
suppresses the RT instability growth rate. Alternatively, this suppressing
effect of strong coupling on RT instability can also be understood in terms of
shear waves whose phase velocity (${\sqrt{{\eta}/{{\rho_{d}}{\tau_{m}}}}}$) is
proportional to it. These equations also imply that elasticity, in the
viscoelastic fluids, enhances the growth of RT instability by reducing the
coupling strength . The effects of elasticity on the RT instability have been
given in a number of studies Biot (1965); Biot and Odé (1965); Odé (1966);
Poliakov  _et al._ (1993); Naimark and Ismail-Zadeh (1994); Kaus and Becker
(2007); Guido  _et al._ (2010). Most of these studies showed that elasticity
speedup the growth of the RT instability. It should be noted that the
analytical dispersion relations for dusty plasma medium has also been obtained
in 2015 by Avinash et. al.Avinash and Sen (2015). However, we had presented
this result earlier in 2014 in a conferenceDharodi and Das (2014).

##  III Numerical simulation

For the purpose of numerical simulation the generalized momentum Eq. (II) has
been expressed as a set of following two coupled convective equations,

$\displaystyle{{\rho_{d}}\left(\frac{\partial{\vec{v}_{d}}}{\partial{t}}+{\vec{v}_{d}}{\cdot}\nabla{\vec{v}_{d}}\right)}+{\rho_{d}}\vec{g}+{\rho_{c}}\nabla\phi_{d}={\vec{\psi}}$
(22) $\frac{\partial{\vec{\psi}}}{\partial
t}+\vec{v}_{d}\cdot\nabla{\vec{\psi}}={\frac{\eta}{\tau_{m}}}{\nabla^{2}}{\vec{v}_{d}}-{\frac{\vec{\psi}}{\tau_{m}}}{.}$
(23)

For two-dimensional (2D) studies all the variables are functions of $x$ and
$y$ only, $i.e$ ${\vec{\psi}}(x,y)$, ${\vec{v}_{d}}(x,y)$, ${\rho_{d}}(x,y)$
and $\vec{g}(y)$. From Eq. (22) it is clear that the quantity
${\vec{\psi}}(x,y)$ is the strain created in the elastic medium by the time-
varying velocity fields. We re-write the Eq. (22) using the equilibrium
condition of ${\rho_{d0}}{g}=-{\rho_{c}}{\nabla\phi_{d0}}$ in terms of
perturbed density Eq. (5).

$\displaystyle\frac{\partial{\vec{v}_{d}}}{\partial{t}}+{\vec{v}_{d}}{\cdot}\nabla{\vec{v}_{d}}+{\frac{\rho_{d1}}{\rho_{d}}}\vec{g}+{\frac{\rho_{c}}{\rho_{d}}}{\nabla{\phi_{d1}}}=\frac{\vec{\psi}}{\rho_{d}}{.}$
(24)

Taking the curl of Eq. (24) and using the Boussinesq approximation, we get

$\frac{\partial{\xi_{z}}}{\partial
t}+\left(\vec{v}_{d}\cdot\vec{\nabla}\right){\xi_{z}}={\frac{1}{\rho_{d0}}}{{\nabla}{\times}{\rho_{d1}}{\vec{g}}}+{\nabla}{\times}{\frac{\vec{\psi}}{\rho_{d}}}{.}$
(25)

Here, ${\xi_{z}}(x,y)=\left[\vec{\nabla}{\times}{\vec{v}_{d}}(x,y)\right]_{z}$
is the $z$ (and the only finite) component of vorticity for the 2D case
considered here. It is normalised by dust plasma frequency. The acceleration
$\vec{g}$ is applied opposite to the fluid density gradient $-g\hat{y}$. The
coupled set of model equations for numerical simulations are the following:

$\frac{\partial\rho_{d}}{\partial
t}+\left(\vec{v}_{d}\cdot\nabla\right)\rho_{d}=0{,}$ (26)
$\frac{\partial{\vec{\psi}}}{\partial
t}+\left(\vec{v}_{d}\cdot\vec{\nabla}\right){\vec{\psi}}={\frac{\eta}{\tau_{m}}}{\nabla^{2}}{\vec{v}_{d}}-{\frac{\vec{\psi}}{\tau_{m}}}{,}$
(27) $\frac{\partial{\xi}_{z}}{\partial
t}+\left(\vec{v}_{d}\cdot\vec{\nabla}\right){{\xi}_{z}}=-{\frac{g}{\rho_{d0}}}{\frac{\partial{\rho_{d1}}}{\partial
x}}+{\frac{\partial}{\partial
x}}\left({\frac{\psi_{y}}{\rho_{d}}}\right)-{\frac{\partial}{\partial
y}}\left({\frac{\psi_{x}}{\rho_{d}}}\right){.}$ (28)

We have used the LCPFCT method (Boris et al. Boris  _et al._ (1993)) to evolve
the coupled set of Eqs. (26), (27) and (25) for various kinds of density
profiles. This method is based on a finite difference scheme associated with
the Flux-Corrected Algorithm. The velocity at each time step is updated by
using the Poisson’s equation
${\nabla^{2}}{\vec{v}_{d}}=-{\vec{\nabla}}{\times}{\vec{\xi}}$, which has been
solved by using the FISPACK package Swarztrauber, Sweet, and Adams (1999). In
the subsequent sections we will present the results of 2D numerical
simulations of gravitational and buoyancy-driven instabilities and discuss how
the growth rate of these instabilities gets suppressed with increasing
coupling strength as expected analytically in Eqs. (18)-(20). Boundary
conditions are periodic in the horizontal direction (x-axis) and absorbing
(non-periodic a) along the vertical (y-axis) direction throughout for all
simulation studies.

##  IV Rayleigh Taylor instability

We consider two types of inhomogeneous density profiles, namely (A) with sharp
interface between two different densities where heavier fluid $\rho_{d02}$
(upper) is supported by lighter one $\rho_{d01}$ (lower) as shown in Fig. 1(a)
and (B) with gradually increasing density gradient along the y-axis as shown
in Fig. 1(b). In colorbar, letter H is the acronym of the heavy density region
and L stands for lighter density region.

Figure 1: Initial density profiles for Rayleigh-Taylor at time t=0. (a) sharp
density interface profile and (b) gradually increasing density gradient along
the vertical y-axis; where H stands for heavy density and L for lighter one.

For both cases A and B, the system length is $lx=ly=2{\pi}$ units and density
gradient is chosen along the vertical y-axis opposite to gravitational
acceleration $g(=10)$. A sufficient electric field is applied for equilibrium
balance against gravity to levitate the dust particles, we had a discussion
about it in Section II.

###  IV.1 Sharp interface

For case A (Fig. 1(a)), we consider a system which consists of two
incompressible fluids of constant densities $\rho_{d01}$=1 for $-{\pi}\leq
y\leq 0$ (lower half) and $\rho_{d02}$=2 for $0=y\leq{\pi}$ (upper half). The
denser fluid $\rho_{d02}$ is positioned above the lighter fluid of density
$\rho_{d01}$. This equilibrium configuration remains stable in absence of
gravity. To hasten the evolution of our equilibrated system under gravity, we
impose a small amplitude sinusoidal perturbation on the density profile at the
interface ($y=0$) $i.e.$

$\rho_{d1}={\phi_{0}}{cos({\mathrm{k}}_{x}x)}exp(-y^{2}/{\epsilon^{2}}){.}$
(29)

Here ${\mathrm{k}}_{x}$ is the wavenumber associated with the perturbation.
Also, the parameter $\epsilon$ defines the width of the perturbation along
$y$. For our simulations we have chosen these parameters as ${\phi_{0}}=0.01$,
$\epsilon=0.1$ and ${\mathrm{k}}_{x}=1$. The total density is
$\rho_{d}=\rho_{d0s}+\rho_{d1}$, where the suffix $s$ corresponds to $1$ and
$2$ which are labeled for the low and high density of the dust fluid,
respectively. The evolution of RT instability is reproduced for the inviscid
HD fluid from our code (see Fig. 2).

Figure 2: Time evolution of RT instability for inviscid HD fluid at the sharp
interface of two different densities. The heavy fluid (upper half) penetrate
into the lighter fluid (lower half) and the lighter fluid rises above into the
domain of heavy fluid. The RT instability perturbations first grows
exponentially then saturate to form the bubbles and spikes types of
structures.

The initial perturbation on the interface was given at ${\mathrm{k}}_{x}=1$
which is observed. The heavy fluid is observed to penetrate into the lighter
fluid and the lighter fluid rises above into the domain of heavy fluid.
Subsequently, the falling heavy fluid develops rolls at the edges acquiring a
mushroom like profile. These rolls grow up with time as the heavy fluid keeps
falling with time. Such an evolution of RT instability is a well known result
of HD fluid Hosseini  _et al._ (2017); Talat  _et al._ (2018), thereby
substantiating our numerical simulation. In order to see the pure viscous
damping effect on the growth of RT instability, we have simulated the viscous
HD fluid of $\eta=0.1$ with the value of $\tau_{m}=0$ shown in Fig. 3.

Figure 3: Time evolution of RT instability for pure viscous HD fluid ($\eta$=
0.1; $\tau_{m}$=0) at the sharp interface of two different densities. Here,
due to the viscous damping force, the RT instability grows at a considerably
slower pace compared to the inviscid HD fluid (Fig. 2).

Here, the instability grows at a considerably slower pace compared to the
inviscid HD fluid (Fig. 2) because of the viscous damping force. In the
interest of elastic effects on the growth of instability we have also
considered two viscoelastic fluids with two different values of $\tau_{m}$ (20
and 5) for the same value of $\eta=0.1$ (viscous damping is similar) shown in
Fig. 4.

Figure 4: The growth of RT instability at the sharp interface of two
viscoelastic fluids of different densities. (a) $\eta$= 0.1; $\tau_{m}$=20,
and (b) $\eta$= 0.1; $\tau_{m}$=5. RT instability gets weaker with increasing
coupling strength i.e. ${\eta}/{\tau_{m}}$ as expected analytically from Eq.
(20). Moreover, a comparison between viscous HD fluid (Fig. 3) and
viscoelastic fluids (Fig. 4) evidents that the appearance of elasticity
$\tau_{m}$ speed up the growth of RT instability by reducing the viscous
fraction ($\eta$/$\tau_{m}$).

A comparison between viscous HD fluid (Fig. 3), viscoelastic fluids (Fig. 4)
evidents that the RT instability in the dust fluid phenomenologically similar
to the HD case as observed in dusty plasma experiment by Pacha et al. Pacha
_et al._ (2012). It can also be noticed that the appearance of elasticity
$\tau_{m}$ speed up the growth of RT instability by reducing the viscous
fraction ($\eta$/$\tau_{m}$). Such reducing effect of elasticity on RT
instability has also been observed in many studies Biot (1965); Biot and Odé
(1965); Odé (1966); Poliakov  _et al._ (1993); Naimark and Ismail-Zadeh
(1994); Kaus and Becker (2007); Guido  _et al._ (2010). Or the RT instability
gets weaker with increasing the value of coupling strength ($\eta$/$\tau_{m}$)
as expected analytically from Eq. (20). In Fig. 5, for the quantitative
comparison, the density interface growth for all these has also been evaluated
and found to similar with the above observations.

Figure 5: The density interface growth (distance from the initial interface)
in vertical downward direction with time. The RT instability gets weaker with
increasing the ratio ${\eta}/{\tau_{m}}$ as expected analytically from Eq.
(20).

###  IV.2 Gradual variation of density profile

Next, the density of the dust fluid is stratified against gravity and has a
gradually changing profile as shown in Fig. 1(b).

$\rho_{d}={\rho_{d0}}+{\phi_{0}}exp({\sigma}y)+{\rho_{d1}}{.}$ (30)

Here $\rho_{d0}$ is constant background density. The value of parameters taken
for the present case are $\rho_{d0}$=5, amplitude of inhomogeneity
$\phi_{0}$=0.1 and $\sigma$=0.025 decide the ramp of inhomogeneity in
background density. The sinusoidal perturbation ${\rho_{d1}}$ (see Eq. 29) is
imposed at y = 0, which can also be clearly seen in Fig. 1(b). In Fig. 6, we
have compared two different viscoelastic fluids with same $\eta$
(${\eta}=0.1,{\tau_{m}}=20,{\eta}/{\tau_{m}}$=0.005 for Fig. 6(a),
(${\eta}=0.1,{\tau_{m}}=5,{\eta}/{\tau_{m}}$=0.02 for Fig. 6(b)). From the
comparison of Fig. 6(a) and Fig. 6(b), it is clear that at each time step, the
growth of RT instability gets weaker with increasing ${\eta}/{\tau_{m}}$ as
suggested by Eq. (18) or the reducing effect of elasticity to the viscous
fraction ($\eta$/$\tau_{m}$) enhances the RT growth. The elastic behaviour of
the system is known to produce TS waves in the medium. Similar observations
have been made for the sharp density gradient profiles (Subsection IV.1). It
is clear that in all these cases the evolution is identical suggesting that
the shear wave associated with strong coupling is operational.

Figure 6: The time evolution of RT instability for two viscoelastic fluids
having gradually increasing densities against the gravity. (a) Viscoelastic
fluid with $\eta$= 0.1; $\tau_{m}$=20, and (b) viscoelastic fluid with $\eta$=
0.1; $\tau_{m}$=5. RT instability gets weaker with increasing coupling
strength i.e. ${\eta}/{\tau_{m}}$ as analytically suggested by Eq. (18).

##  V Buoyancy-driven evolution

We consider the two types of density inhomogeneity profiles, namely (A) an
initially static and circular bubble whose density is less than that of the
surrounding quiescent Newtonian/viscoelastic fluid as shown in Fig. 7(a) and
(B) initially static and a circular droplet whose density is higher than that
of the surrounding medium as shown in Fig. 7(b). In both the cases the
variation in density is symmetric about an axis which is perpendicular to
simulation plane.

Figure 7: Initial densities’ profiles at time t=0. (a) Bubble of less density
inside the heavier fluid and (b) droplet of high density inside the lighter
fluid. In colorbar, letter H is the acronym of the heavy density region and L
stands for lighter density region.

For both cases, we have considered a system of length $lx=ly=12{\pi}$ units,
approximating the gravitational acceleration as $g=10$. Under the influence of
gravity, the less dense regions (bubble) have a tendency to rise upwards,
whereas the higher density (droplet) sinks against the background low density
fluid. Now let us discuss the equations of model to develop the some simple
qualitative understanding of simulations results. For an inviscid flow, the
vorticity Eq. (25)/Eq. (28) in hydrodynamic limit, becomes

$\frac{\partial{\xi_{z}}}{\partial
t}+\left(\vec{v}_{d}\cdot\vec{\nabla}\right){\xi_{z}}={\frac{1}{\rho_{d0}}}{{\nabla}{\times}{\rho_{d1}}{\vec{g}}}=-{\frac{g}{\rho_{d0}}}{\frac{\partial{\rho_{d1}}}{\partial
x}}{.}$ (31)

The spatial variation in density $i.e.$
$-{({g}/{\rho_{d0}})}{({\partial{\rho_{d1}}}/{\partial x})}$ acts as a buoyant
force. In order to simulate the case of rising/falling bubble/droplet, we have
considered a Gaussian density profile given by,

$\rho_{d1}={\rho^{\prime}_{0}}exp\left(\frac{-\left({\left(x-x_{c}\right)^{2}+(y-y_{c})^{2}}\right)}{a^{2}_{c}}\right){,}$
(32)

where $a_{c}$ is the bubble/droplet core radius. As the simulation begins, the
buoyant force comes into a picture which in turn induces the dipolar vorticity
(a pair of counter rotating vortices) for our considered Gaussian density
profile $i.e.$

$-{\frac{\partial{\rho_{d1}}}{\partial
x}}=2{\left(x-x_{c}\right)}{\rho_{d1}}{.}$ (33)

The net resultant of counter rotating pair of vortices will decide the net
flow direction of the medium Horstmann  _et al._ (2014) $i.e$ vertical upward
(bubble) or vertical downward (droplet). Figure 8, corresponds to the Fig. 7,
shows the schematic diagram of vorticities (using Eq. (33)), it is clear that
two similar counter-rotating lobes cause the net propagation of bubble/droplet
in the vertical upward/downward direction (indicated by arrows).

Figure 8: Schematic vorticity diagram (a pair of counter counter-rotating
lobes) of the bubble and droplet. (a) Dipolar vorticity due to the buoyant
force on bubble w.r.t. Fig. 7(a). Here, the net rotation causes the motion in
vertical upward direction, indicated by vertical upward arrow. Similarly,
vorticity in (b) w.r.t. Fig. 7(b) causes the net downward flow for the
droplet.

While in viscoelastic fluids, in addition to the buoyant force term, the
vorticity Eq. (25)/Eq. (28) includes
${\nabla}{\times}{({\vec{\psi}}/{\rho_{d}})}$ as second term (RHS) which
incorporate the shear waves into the medium moving with phase velocity
$\sqrt{{\eta}/{{\rho_{d}}{\tau_{m}}}}$.

###  V.1 Rising bubble

In order to simulate the case of rising bubble the numerical simulation has
been carried out for ${a_{c}}$=2.0, ${\rho^{\prime}_{0}}=-0.5$ and
${x_{c}}=0$, ${y_{c}}=-4\pi$ in above Eq. 32. The homogeneous background
density with ${\rho_{d0}}=5$ has been taken. Thus, the total density is
$\rho_{d}=\rho_{d0}+\rho_{d1}$. Figure 9 displays the time evolution of the
bubble density profile of inviscid fluid. Initially at time $t=0$, the bubble
is axisymmetric (see Fig. 7(a)) and as the system evolves, the bubble simply
rises without any significant change in shape for a short period of time. At a
later stage, the initially circular profile assumes a crescent shape as
evident from the first subplot of Fig. 9. Further, as time progresses, a
mushroom-like structure which is characteristic of RT instability begins to
appear along with rolling and mixing at the edge of the bubble. This mushroom-
like structure then breaks into two distinct elliptical lobes as evident from
the figure. These lobes propagate forward as a single entity leaving behind
wake-like structure in background fluid. Next, with the understanding of
bubble evolution for HD fluids, the time evolution of bubble for viscoelastic
fluids has been performed.

Figure 9: Time evolution of bubble density (see Fig. 12 for the respective
vorticity) profile for inviscid HD fluid.

To envisage the effect of coupling strength, we compare one viscoelastic fluid
(Fig. 10(a)) having ${\eta}=2.5,{\tau_{m}}=20$ i.e. $\eta/\tau_{m}$=0.125 with
another viscoelastic fluid (Fig. 10(b)) having higher coupling strength i.e.
${\eta}$=5, ${\tau_{m}}$=10, $\eta/\tau_{m}$=0.5.

Figure 10: Evolution of bubble density profile (see Fig. 11 for the respective
vorticity evolution) in time for two viscoelastic fluids (a) $\eta$= 2.5;
$\tau_{m}$=20, and (b) $\eta$= 5; $\tau_{m}$=10. The growth of bubble gets
reduced (rising rate) with increasing coupling strength.

We observe a decrease in rising rate of bubble with increasing coupling
strength. Also, the elliptical lobes start moving apart in horizontal
direction earlier in stronger coupled medium. At the same time, the rolling
rate of lobes increases causing the expansion of lobes and disappearance of
the dragging tail. Thus, we observe that similar to the RT instability the
growth of bubble gets reduced (rising rate) with increasing coupling strength
and this reducing effect can be easily visualised through the respective
vorticity evolution given in Fig. 11.

Figure 11: Time evolution of bubble vorticity profile (see Fig. 10 for the
respective density evolution) for two viscoelastic fluids (a) $\eta$= 2.5;
$\tau_{m}$=20, and (b) $\eta$= 5; $\tau_{m}$=10. There is emission of TS wave
surrounding the vorticity lobes whose phase velocity is proportional to the
coupling strength. It is clear that the increase of phase velocity of TS waves
with coupling strength suppresses the BD instability ($i.e$ reduction in
rising rate of lobes).

In Fig. 11, the vorticity subplots in first and second row correspond to the
density profiles shown in Fig. 10(a) and Fig. 10(b), respectively. It is
evident from Fig. 11 that there is emission of TS waves surrounding the
vorticity lobes for visco-elastic fluids and no such TS waves exist for
inviscid HD fluid (Fig. 12). The emission of waves from either lobe pushes the
other lobe in the direction perpendicular to the direction of propagation of
the entire structure and as a result the lobes get well separated with time
and proportional to the coupling strength. Besides this lobe separation, the
emission of TS wave significantly reduced the strength of dipole thereby
impeding the dipole propagation and this reduction in dipole propagation is
also proportional to the coupling strength. The relative observations of Fig.
11(a) and Fig. 11(b) clearly reflect it through the TS enclosure which is
larger for the $\eta$= 5; $\tau_{m}$=10 than $\eta$= 2.5; $\tau_{m}$=20.

Figure 12: Evolution of bubble vorticity profile in time for inviscid HD fluid
(see Fig. 9 for the respective density).

Thus, we conclude that the increase of phase velocity of TS waves with
coupling strength suppresses the BD instability ($i.e$ reduction in vertical
propagation direction and enhancement in horizontal separation between lobes).
To understand the separate effects of viscosity and elasticity on BD growth
the snapshots of density/vorticity snapshots at the same timestep
$({\approx}20)$ are shown in Fig. 13/Fig. 14. From Fig. 13, the damping effect
of viscosity on the BD instability can be easily noticed, and elasticity
reduces the viscous effects by reducing the viscous fraction.

Figure 13: Viscosity $\eta$ vs elasticity $\tau_{m}$. The snapshots of bubble
density profiles at a particular timestep (see Fig. 14 for the respective
vorticity). The BD instability gets weaker with increasing viscosity (in the
vertical direction with a fixed value of $\tau_{m}$), and elasticity only
reduces the viscous effects (along the horizontal direction with a fixed value
of $\eta$).

The role of TS waves can be visualised from the vorticity snapshots given in
Fig. 14.

Figure 14: Viscosity $\eta$ vs elasticity $\tau_{m}$. The snapshots of bubble
vorticity profiles at a particular timestep (see Fig. 13 for the respective
density). The BD instability gets weaker with increasing viscosity (in the
vertical direction with a fixed value of $\tau_{m}$), and elasticity only
reduces the viscous effects (along the horizontal direction with a fixed value
of $\eta$).

From viscoelastic theory it well known that the strong coupling increases the
material properties like its rigidity (giving rise to shear waves). So this
extra rigidity should give rise to an extra torque thereby increasing the
rotation. This would imply that in strong coupling limit rotation should be
faster. This effect can be observed from Fig. 13/Fig. 14, in strong coupling
limit the expansion of lobes is higher, and disappearance of the dragging tail
in density evolution. Also, the vortex evolution is similar for the same value
of coupling strength (${\eta}/{\tau_{m}}$).

###  V.2 Falling droplet

For simulating the falling droplet $i.e$ case B (Fig. 7(b)), the values of
parameters in Eq. 32 for the droplet density profile are ${a_{c}}$=2.0,
${\rho^{\prime}_{0}}=0.5$ and ${x_{c}}=0$, ${y_{c}}=4\pi$. We considered the
homogeneous background density as ${\rho_{d0}}=5$. The total density is
$\rho_{d}=\rho_{d0}+\rho_{d1}$. Figure 15 shows the dynamics of the droplet,
falling in tranquil inviscid HD fluid. The suspended drop starts to fall. As
time passes, the drop breaks up forming first a semicircular structure then a
two lobes similar to the case of bubble.

Figure 15: Time evolution of droplet density (see Fig. 18 for the respective
vorticity evolution) profile for inviscidHD fluid.

Figure 16 reveals the different stages of droplet for viscoelsatic fluids. As
done for rising bubble, we shall now compare Fig. 16(a) having
${\eta}=3.5,{\tau_{m}}=10$ $i.e.$ $\eta/\tau_{m}$=0.35) to Fig. 16(b) having
${\eta}=5,{\tau_{m}}=10$ $i.e.$ $\eta/\tau_{m}$=0.5 .

Figure 16: Evolution of droplet density profile (see Fig. 17 for the
respective vorticity evolution) in time for two viscoelastic fluids (a)
$\eta$= 2.5; $\tau_{m}$=20, and (b) $\eta$= 5; $\tau_{m}$=10. The growth of
droplet gets reduced (falling rate) with increasing coupling strength.

For a higher coupling strength shown in Fig. 16(b), the downward motion is
slower. The transverse dimension is even larger and the dragging tail does not
seem to appear at all. Again, similar to the case of bubbles, the vorticity
plots are provided for observing a clear role of TS waves i.e. pushing the two
lobes apart, as seen in Fig. 17.

Figure 17: Time evolution of droplet vorticity profile (see Fig. 16 for the
respective density evolution) for two viscoelastic fluids (a) $\eta$= 2.5;
$\tau_{m}$=20, and (b) $\eta$= 5; $\tau_{m}$=10. There is emission of TS wave
surrounding the vorticity lobes whose phase velocity is proportional to the
coupling strength. It is clear that the increase of phase velocity of TS waves
with coupling strength suppresses the BD instability ($i.e$ reduction in
falling rate of lobes).

For the inviscid HD fluid the corresponding vorticity subplots (Fig. 18) to
the mentioned density profile (Fig. 15), no such wave are observed which can
constraint the falling rate. Thus, we conclude that the increase of phase
velocity of TS waves with coupling strength suppresses the buoyant nature of a
droplet.

Figure 18: Time evolution of droplet vorticity (see Fig. 15 for the respective
density evolution) profile for inviscid HD fluid.

To understand the separate effects of viscosity and elasticity on BD growth
the snapshots of density/vorticity snapshots at the same timestep
$({\approx}20)$ are shown in Fig. 19/Fig. 20.

Figure 19: Viscosity $\eta$ vs elasticity $\tau_{m}$. The snapshots of droplet
density profiles at a particular timestep (see Fig. 20 for the respective
vorticity). The BD instability gets weaker with increasing viscosity (in the
vertical direction with a fixed value of $\tau_{m}$), and elasticity only
reduces the viscous effects (along the horizontal direction with a fixed value
of $\eta$).

The role of TS waves can easily be understood from the vorticity snapshots
given in Fig. 20.

Figure 20: Viscosity $\eta$ vs elasticity $\tau_{m}$. The snapshots of droplet
vorticity profiles at a particular timestep (see Fig. 19 for the respective
density). The BD instability gets weaker with increasing viscosity (in the
vertical direction with a fixed value of $\tau_{m}$), and elasticity only
reduces the viscous effects (along the horizontal direction with a fixed value
of $\eta$).

Similar to the bubble evolution, it can be easily noticed that the BD
instability gets weaker with increasing viscosity, and elasticity reduces the
viscous effects. Also, in strong coupling limit the expansion of lobes is
higher, and disappearance of the dragging tail in density evolution. Also, the
vortex evolution is similar for the same value of coupling strength
(${\eta}/{\tau_{m}}$). This means the spatio-temporal evolution of
vortices/instabilities depends on the coupling strength.

##  VI Observations and conclusion

We present a detailed discussion on the evolution of Rayleigh-Taylor
instability and Buoyancy-Driven instability in two-dimensional strongly
coupled dusty plasmas (SCDPs) in the presence of gravity. The SCDPs have been
treated as viscoelastic fluids in the framework of incompressible generalized
hydrodynamics (i- GHD) fluid model. The density stratification against gravity
leading to Rayleigh-Taylor instability. For the buoyancy-driven evolution, we
consider a spatially localized high/low (droplet/bubble) density region placed
in a constant background density medium in the presence of gravity. Some main
observations are as follows:

  * •

RT instability in the dust fluid phenomenologically similar to the
hydrodynamic case as observed in dusty plasma experiment by Pacha et al. Pacha
_et al._ (2012).

  * •

Our studies, analytically as well as numerical simulations, show that both the
instabilities get suppressed with increasing coupling strength of the medium.

  * •

The appearance of elasticity $\tau_{m}$ speed up the growth of RT instability
by reducing the effect of viscosity.

  * •

The spatio-temporal evolution of vortices/instabilities depends on the value
of coupling parameter.

The addition of third dimension (3D) will make these results more general.
These observations would be presented in a subsequent publication. For the
hydrodynamic fluids Youngs et. al. reported the difference between 2D and 3D
simulation of RT mixing Youngs (1991) and observed that dissipation of density
fluctuations is less in 2D than in 3D. To investigate three-dimensional
effects on the Rayleigh–Taylor instability, numerically observe the two-layer
roll-up phenomenon of the heavy fluid, which does not occur in the two-
dimensional case. Lee and Kim (2013)

## References

  * Rayleigh (1900) L. Rayleigh, “Investigation of the character of the equilibrium of an incompressible heavy,” (1900).
  * Taylor (1950) G. I. Taylor, “The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. i,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 201, 192–196 (1950).
  * (3) S. Chandrasekhar, “Hydrodynamic and hydromagnetic stability, dover, new york, 1981,” .
  * Youngs (1984) D. L. Youngs, “Numerical simulation of turbulent mixing by rayleigh-taylor instability,” Physica D: Nonlinear Phenomena 12, 32–44 (1984).
  * Waddell, Niederhaus, and Jacobs (2001) J. Waddell, C. Niederhaus, and J. W. Jacobs, “Experimental study of rayleigh–taylor instability: low atwood number liquid systems with single-mode initial perturbations,” Physics of Fluids 13, 1263–1273 (2001).
  * Livescu (2004) D. Livescu, “Compressibility effects on the rayleigh–taylor instability growth between immiscible fluids,” Physics of fluids 16, 118–127 (2004).
  * Guo and Tice (2010) Y. Guo and I. Tice, “Linear rayleigh–taylor instability for viscous, compressible fluids,” SIAM Journal on Mathematical Analysis 42, 1688–1720 (2010).
  * Baldwin, Scase, and Hill (2015) K. A. Baldwin, M. M. Scase, and R. J. Hill, “The inhibition of the rayleigh-taylor instability by rotation,” Scientific reports 5, 11706 (2015).
  * Lyubimova, Vorobev, and Prokopev (2019) T. Lyubimova, A. Vorobev, and S. Prokopev, “Rayleigh-taylor instability of a miscible interface in a confined domain,” Physics of Fluids 31, 014104 (2019).
  * Dickinson  _et al._ (1962) H. Dickinson, W. Bostick, J. DiMarco, and S. Koslov, “Experimental study of rayleigh-taylor instability in plasma,” The Physics of Fluids 5, 1048–1056 (1962).
  * Ariel and Bhatia (1970) P. Ariel and P. Bhatia, “Rayleigh-taylor instability of a rotating plasma,” Nuclear Fusion 10, 141 (1970).
  * Takabe  _et al._ (1985) H. Takabe, K. Mima, L. Montierth, and R. Morse, “Self-consistent growth rate of the rayleigh–taylor instability in an ablatively accelerating plasma,” The Physics of fluids 28, 3676–3682 (1985).
  * Mikhailenko, Mikhailenko, and Weiland (2002) V. S. Mikhailenko, V. V. Mikhailenko, and J. Weiland, “Rayleigh–taylor instability in plasmas with shear flow,” Physics of Plasmas 9, 2891–2895 (2002).
  * Ma  _et al._ (2006) J. Ma, Y.-h. Chen, B.-x. Gan, and M. Yu, “Behavior of the rayleigh–taylor mode in a dusty plasma with rotational and shear flows,” Planetary and Space Science 54, 719–725 (2006).
  * Sen, Fukuyama, and Honary (2010) S. Sen, A. Fukuyama, and F. Honary, “Rayleigh taylor instability in a dusty plasma,” Journal of Atmospheric and Solar-Terrestrial Physics 72, 938–942 (2010).
  * Weber  _et al._ (2014) C. Weber, D. Clark, A. Cook, L. Busby, and H. Robey, “Inhibition of turbulence in inertial-confinement-fusion hot spots by viscous dissipation,” Physical Review E 89, 053106 (2014).
  * Kessler, Dahlburg, and Ganguli (2014) D. A. Kessler, R. B. Dahlburg, and G. Ganguli, “Gravitational instability and shear stabilization in a dusty plasma layer,” in  _2014 United States National Committee of URSI National Radio Science Meeting (USNC-URSI NRSM)_ (IEEE, 2014) pp. 1–1.
  * Haines  _et al._ (2014) B. Haines, E. Vold, K. Molvig, R. Rauenzahn, and C. Aldrich, “Plasma transport in rayleigh-taylor and kelvin-helmholtz instabilities,” Phys. Plasmas 21, 092306 (2014).
  * Hoshoudy (2014) G. Hoshoudy, “Rayleigh-taylor instability in magnetized plasma,” World Journal of Mechanics 4, 260 (2014).
  * Khiar  _et al._ (2019) B. Khiar, G. Revet, A. Ciardi, K. Burdonov, E. Filippov, J. Béard, M. Cerchez, S. Chen, T. Gangolf, S. Makarov, _et al._ , “Laser-produced magnetic-rayleigh-taylor unstable plasma slabs in a 20 t magnetic field,” Physical Review Letters 123, 205001 (2019).
  * Norman  _et al._ (1981) M. L. Norman, L. Smarr, M. D. Smith, and J. R. Wilson, “Hydrodynamic formation of twin-exhaust jets,” Astrophysical 247, 52–58 (1981).
  * Arnett  _et al._ (1989) W. D. Arnett, J. N. Bahcall, R. P. Kirshner, and S. E. Woosley, “Supernova 1987a,” Annual Review of Astronomy and Astrophysics 27, 629–700 (1989).
  * Plag and Jüttner (1995) H.-P. Plag and H.-U. Jüttner, “Rayleigh-taylor instabilities of a self-gravitating earth,” Journal of geodynamics 20, 267–288 (1995).
  * Allen and Hughes (1984) A. J. Allen and P. A. Hughes, “The Rayleigh-Taylor instability in astrophysical fluids,” mnras 208, 609–621 (1984).
  * Zingale  _et al._ (2005) M. Zingale, S. E. Woosley, C. A. Rendleman, M. S. Day, and J. B. Bell, “Three-dimensional numerical simulations of rayleigh-taylor unstable flames in type ia supernovae,” The Astrophysical Journal 632, 1021–1034 (2005).
  * Keskinen  _et al._ (1981) M. Keskinen, E. Szuszczewicz, S. Ossakow, and J. Holmes, “Nonlinear theory and experimental observations of the local collisional rayleigh-taylor instability in a descending equatorial spread f ionosphere,” Journal of Geophysical Research: Space Physics 86, 5785–5792 (1981).
  * Beale and Reitz (1999) J. C. Beale and R. D. Reitz, “Modeling spray atomization with the kelvin-helmholtz/rayleigh-taylor hybrid model,” Atomization and sprays 9 (1999).
  * Kong, Senecal, and Reitz (1999) S. C. Kong, P. K. Senecal, and R. D. Reitz, “Developments in spray modeling in diesel and direct-injection gasoline engines,” Oil & Gas Science and Technology 54, 197–204 (1999).
  * Chertkov (2003) M. Chertkov, “Phenomenology of rayleigh-taylor turbulence,” Phys. Rev. Lett. 91, 115001 (2003).
  * Kadau  _et al._ (2004) K. Kadau, T. C. Germann, N. G. Hadjiconstantinou, P. S. Lomdahl, G. Dimonte, B. L. Holian, and B. J. Alder, “Nanohydrodynamics simulations: An atomistic view of the rayleigh-taylor instability,” Proceedings of the National Academy of Sciences of the United States of America 101, 5851–5855 (2004).
  * Hoshoudy (2011) G. A. Hoshoudy, “Quantum effects on the rayleigh-taylor instability of viscoelastic plasma model through a porous medium.” Journal of Modern Physics 2 (2011).
  * Srinivasan and Tang (2013) B. Srinivasan and X.-Z. Tang, “The mitigating effect of magnetic fields on rayleigh-taylor unstable inertial confinement fusion plasmas,” Physics of Plasmas 20, 056307 (2013).
  * Debnath (1994) L. Debnath,  _Nonlinear water waves_ (Academic Press, 1994).
  * Abarzhi (2010) S. I. Abarzhi, “Review of theoretical modelling approaches of rayleigh–taylor instabilities and turbulent mixing,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 368, 1809–1828 (2010).
  * Abarzhi (2008) S. I. Abarzhi, “Coherent structures and pattern formation in the rayleigh-taylor turbulent mixing,” Phys. Scr. 78 (2008).
  * Abarzhi and Rosner (2010) S. I. Abarzhi and R. Rosner, “Comparative study of approaches for modeling rayleigh-taylor turbulent mixing,” Phys. Scr. T142, 1 (2010).
  * Frenkel (1955) J. Frenkel,  _Kinetic Theory Of Liquids_ (Dover Publications, 1955).
  * Kaus and Becker (2007) B. J. Kaus and T. W. Becker, “Effects of elasticity on the rayleigh–taylor instability: implications for large-scale geodynamics,” Geophysical Journal International 168, 843–862 (2007).
  * Guido  _et al._ (2010) B. Guido, A. Mazzino, S. Musacchio, and L. Vozella, “Rayleigh–taylor instability in a viscoelastic binary fluid,” Journal of Fluid Mechanics 643, 127–136 (2010).
  * Kaw and Sen (1998) P. K. Kaw and A. Sen, “Low frequency modes in strongly coupled dusty plasmas,” Physics of Plasmas 5, 3552–3559 (1998).
  * Kaw (2001) P. K. Kaw, “Collective modes in a strongly coupled dusty plasma,” Physics of Plasmas 8, 1870–1878 (2001).
  * Tiwari  _et al._ (2012) S. K. Tiwari, A. Das, D. Angom, B. G. Patel, and P. Kaw, “Kelvin-helmholtz instability in a strongly coupled dusty plasma medium,” Physics of Plasmas 19, 073703 (2012).
  * S. Dharodi, K. Tiwari, and Das (2014) V. S. Dharodi, S. K. Tiwari, and A. Das, “Visco-elastic fluid simulations of coherent structures in strongly coupled dusty plasma medium,” Physics of Plasmas 21, 073705 (2014).
  * Tiwari  _et al._ (2014) S. K. Tiwari, V. S. Dharodi, A. Das, B. G. Patel, and P. Kaw, “Evolution of sheared flow structure in visco-elastic fluids,” AIP Conf. Proc 1582, 55 (2014).
  * Dharodi  _et al._ (2016) V. S. Dharodi, A. Das, B. G. Patel, and P. K. Kaw, “Sub-and super-luminar propagation of structures satisfying poynting-like theorem for incompressible generalized hydrodynamic fluid model depicting strongly coupled dusty plasma medium,” Physics of Plasmas 23, 013707 (2016).
  * Das and Kaw (2014) A. Das and P. K. Kaw, “Suppression of rayleigh-taylor instability in strongly coupled plasmas,” Physics of Plasmas 21 (2014).
  * Dharodi (2020) V. S. Dharodi, “Rotating vortices in two-dimensional inhomogeneous strongly coupled dusty plasmas: shear and spiral-density waves,” arXiv preprint arXiv:2006.04545 (2020).
  * Biot (1965) M. A. Biot,  _Mechanics of incremental deformations_ (1965).
  * Biot and Odé (1965) M. A. Biot and H. Odé, “Theory of gravity instability with variable overburden and compaction,” Geophysics 30, 213–227 (1965).
  * Odé (1966) H. Odé, “Gravitational instability of a multilayered system of high viscosity: Verhandelingen der koninklijke nederlandse akademie van wetenschappen, afdeeling natuurkunde: Reeks 1, wiskunde, natuurkunde, scheikunde, aardkunde, technische wetenschappen,” (1966).
  * Poliakov  _et al._ (1993) A. Poliakov, P. Cundall, Y. Podladchikov, and V. Lyakhovsky, “An explicit inertial method for the simulation of viscoelastic flow: an evaluation of elastic effects on diapiric flow in two-and three-layers models,” in  _Flow and Creep in the Solar System: observations, modeling and Theory_ (Springer, 1993) pp. 175–195.
  * Naimark and Ismail-Zadeh (1994) B. Naimark and A. Ismail-Zadeh, “Gravitational instability of maxwell upper mantle,” Comput. Seis. Geodyn 1, 36–42 (1994).
  * Avinash and Sen (2015) K. Avinash and A. Sen, “Rayleigh-taylor instability in dusty plasma experiment,” Physics of Plasmas 22, 083707 (2015).
  * Dharodi and Das (2014) V. S. Dharodi and A. Das, “Rayleigh-taylor instability in a visco-elastic medium using generalized hydrodynamic model,” (15-19 September 2014), 17th international congress on plasma physics at Instituto Superior Tecnico, Lisbon, Portugal.
  * Boris  _et al._ (1993) J. P. Boris, A. M. Landsberg, E. S. Oran, and J. H. Gardner,  _LCPFCT A flux-corrected transport algorithm for solving generalized continuity equations_ (Technical Report NRL Memorandum Report 93-7192, Naval Research Laboratory, 1993).
  * Swarztrauber, Sweet, and Adams (1999) P. Swarztrauber, R. Sweet, and J. C. Adams, “Fishpack: Efficient fortran subprograms for the solution of elliptic partial differential equations,” UCAR Publication, July (1999).
  * Hosseini  _et al._ (2017) B. S. Hosseini, S. Turek, M. Möller, and C. Palmes, “Isogeometric analysis of the navier–stokes–cahn–hilliard equations with application to incompressible two-phase flows,” Journal of Computational Physics 348, 171–194 (2017).
  * Talat  _et al._ (2018) N. Talat, B. Mavrič, V. Hatić, S. Bajt, and B. Šarler, “Phase field simulation of rayleigh–taylor instability with a meshless method,” Engineering Analysis with Boundary Elements 87, 78–89 (2018).
  * Pacha  _et al._ (2012) K. A. Pacha, J. R. Heinrich, S.-H. Kim, and R. L. Merlino, “Observation of the taylor instability in a dusty plasma,” Physics of Plasmas 19, 014501 (2012).
  * Horstmann  _et al._ (2014) J. T. Horstmann, P. Henningsson, A. L. Thomas, and R. J. Bomphrey, “Wake development behind paired wings with tip and root trailing vortices: consequences for animal flight force estimates,” PloS one 9 (2014).
  * Youngs (1991) D. L. Youngs, “Three-dimensional numerical simulation of turbulent mixing by rayleigh–taylor instability,” Physics of Fluids A: Fluid Dynamics 3, 1312–1320 (1991).
  * Lee and Kim (2013) H. G. Lee and J. Kim, “Numerical simulation of the three-dimensional rayleigh–taylor instability,” Computers & Mathematics with Applications 66, 1466–1474 (2013).


# A Stock Trading System for a Medium Volatile Asset using Multi Layer
Perceptron

Ivan Letteri
Department of Engineering, Computer Science and Maths
University of L’Aquila
L’Aquila, Italy, 67100
<EMAIL_ADDRESS>
Giuseppe Della Penna
Department of Engineering, Computer Science and Maths
University of L’Aquila
L’Aquila, Italy, 67100
<EMAIL_ADDRESS>
Giovanni De Gasperis
Department of Engineering, Computer Science and Maths
University of L’Aquila
L’Aquila, Italy, 67100
<EMAIL_ADDRESS>
Abeer Dyoub
Department of Engineering, Computer Science and Maths
University of L’Aquila
L’Aquila, Italy, 67100
<EMAIL_ADDRESS>
https://www.ivanletteri.it

###### Abstract

Stock market forecasting is a lucrative field of interest with promising
profits but not without its difficulties and for some people could be even
causes of failure. Financial markets by their nature are complex, non-linear
and chaotic, which implies that accurately predicting the prices of assets
that are part of it becomes very complicated. In this paper we propose a stock
trading system having as main core the feed-forward deep neural networks (DNN)
to predict the price for the next 30 days of open market, of the shares issued
by Abercrombie & Fitch Co. (ANF) in the stock market of the New York Stock
Exchange (NYSE).

The system we have elaborated calculates the most effective technical
indicator, applying it to the predictions computed by the DNNs, for generating
trades. The results showed an increase in values such as Expectancy Ratio of
2.112% of profitable trades with Sharpe, Sortino, and Calmar Ratios of 2.194,
3.340, and 12.403 respectively. As a verification, we adopted a backtracking
simulation module in our system, which maps trades to actual test data
consisting of the last 30 days of open market on the ANF asset. Overall, the
results were promising bringing a total profit factor of 3.2% in just one
month from a very modest budget of $100. This was possible because the system
reduced the number of trades by choosing the most effective and efficient
trades, saving on commissions and slippage costs.

_Keywords_ Deep Learning $\cdot$ Statistical Learning $\cdot$ Stock Market
Prediction $\cdot$ Trading System $\cdot$ Algorithmic Trading

##  1 Introduction

Stock market forecasting is considered a research field with promising returns
for investors. However, there are considerable challenges to predicting stock
market trends accurately and precisely enough due to their chaotic and non-
linear nature, and also complexity. At the state of the art, there are
numerous approaches for generating, processing, and optimizing a dataset, such
as (Letteri et al. (2020a)).

Many artificial intelligence methods have been employed to classify cyber
attacks (Letteri et al. (2018)), predict network traffic anomalies (Letteri et
al. (2019a)), course of a disease, and even the price trend in the stock
market. Artificial neural networks (ANNs) to date remain a fairly popular
choice for these kinds of tasks and are widely studied (Letteri et al.
(2019b)), having exhibited good performance (Yetis et al. (2014)).

In particular, although deep neural networks (DNNs) are predominantly used for
image recognition and natural language processing, showing surprising results
(Soniya et al. (2015)), they have also been applied to financial markets for
stock price prediction using textual news analysis (Day and Lee (2016)).

In this paper, we propose a trading system for the stock market that uses
forecasts by four Price Action (PA)-based DNNs to generate trading signals,
and subsequently evaluate their performance in the stock market. This
constitutes part of a work that aims to convert the DNNs forecasts into a
profitable investment and trading system that can be automated, using a robot
advisor.

##  2 Data Set

This section describes the methodology with which the dataset is collected,
the criteria adopted for the statistical analysis of the data collected, the
training process and the series of tests conducted on the DNNs with optimized
processes (Letteri et al. (2020b)). In the next section we will show the
trading rules applied on the forecast output with the related tests on
profitability verified by our backtracking system.

###  2.1 Technical Analysis

Technical analysis (TA) constitutes the type of investment analysis that uses
simple mathematical formulations or graphical representations of the time
series of financial assets to explore trading opportunities. In its
algorithmic form, TA uses the analysis of asset price history series (Wang et
al. (2021)), defined as OHLC, i.e., the opening, highers, lowest and closing
prices of an asset, typically represented with candlesticks charts (see, e.g.
in fig. 1). For each timeframe $t$, the OHLC of an asset is represented as a
4-dimensional vector
$X_{t}=(x^{(o)}_{t},x^{(h)}_{t},x^{(l)}_{t},x^{(c)}_{t})^{T}$, where
$x_{t}^{(l)}>0$, $x_{t}^{(l)}<x_{t}^{(h)}$ and
$x_{t}^{(o)},x_{t}^{(c)}\in[x_{t}^{(l)},x_{t}^{(h)}]$.

Figure 1: Example of candlestick chart.

Some stocks are more volatile than others. For example the shares of a large
blue-chip company may not have large price fluctuations and are therefore said
to have low volatility, whereas the shares of a tech stock may fluctuate often
and therefore have high volatility. There are also stocks with medium
volatility111https://www.fool.com/investing/how-to-invest/stocks/stock-market-
volatility/, and this is the case of the asset that we studied in this
article.

The dataset used in this work consists of historical OHLC prices data from the
New York Stock Exchange (NYSE), the world’s largest stock exchange by trading
volume and the second largest by number of listed companies. Its share volume
was surpassed by NASDAQ in the 1990s, but the total capitalization of the 2800
companies on the NYSE is five times that of the competing technology exchange.
The asset on which we conducted our study is the shares issued by Abercrombie
& Fitch Co. (ANF). From a summary of results for the third quarter ended
October 30, 2021 compared to the third quarter ended October 31, 2020, ANF
reported net sales of $905 million, up 10% year-over-year and 5% compared to
net sales in the third quarter pre-COVID 2019. Digital net sales of $413
million increased 8% compared to last year and 55% compared to pre-COVID 2019
third quarter digital net sales. Gross profit rate of 63.7%, down
approximately 30 basis points from last year and up approximately 360 basis
points from 2019 (source GlobeNewswire222www.globenewswire.com/news-
release/2021/11/23/2339734/0/en/Abercrombie-Fitch-Co-Reports-Third-Quarter-
Results.html).

Figure 2: ANF all trend.

ANF was born at the end of September 1996, and from April to October 2011 it
came close to the All Time High (ATH) at various times, meeting punctually a
resistance. This peculiarity, together with the others previously mentioned,
has led us to observe this stock with attention starting from the end of
October 2011 (see fig. 2) when the downtrend period below $70 begins, with the
further intent to determine a potential opportune moment for a market entry by
buying this asset.

It is certainly an asset with a sometimes controversial trend and consequently
well profitable if rightly analyzed, especially in the particular period of
2020/2021 due to the global pandemic (see fig. 3(a)) and because it lacks the
leverage effect (Black (1976)) that explains the negative correlation between
equity returns and return volatility (see fig. 3(b) and 3(c)) where most of
the measures of the volatility of the asset are negatively correlated with its
return.

Furthermore, in this long period analysed, ANF is not the classic always
profitable stock (e.g. Tesla, Apple, Microsoft, or Bitcoin) to which trivially
apply a passive Buy and Hold
strategy333https://www.investopedia.com/terms/b/buyandhold.asp for gains.

Figure 3: ANF price (a), return (b), and volatility (c) from October 30th 2011
to November 30th 2021.

The dataset consists of the time series of the price of the above mentioned
stocks, over the time period from October 16, 2011 to November 20, 2021, for a
total period of 2537 open market days. The time series of price observations
can be downloaded from Yahoo
Finance444https://it.finance.yahoo.com/quote/ANF?p=ANF&.tsrc and then split
into train and test sets, or directly from our github
repository555https://github.com/IvanLetteri/StockTradingSystem-MLP-Regressor.

Figure 4: OHLC scatter matrix of Validation set from 16th October to 20th
November 2021).

In particular, only the 30 days of the market open was used for the validation
of the values composed by the features characterized by the OHLC prices
aforementioned. In fig. 4, we can see the scatter matrix, which indicates that
there are few positive correlations between them, even looking at the
distribution of each individual variable along the diagonal represented with
the Kernel Density Estimation (KDE) (see Chen (2017)), which denotes a
smoother distribution surface across the plots. Intuitively, KDE has the
effect of smoothing out each data point into a smooth bump. The shape of curve
is determined by the kernel function $K(x)$ which sums overall the bumps to
obtain the density estimator. The regions with many observations, KDE yield a
large value because there are many bumps around. To the other side, for
regions with few observations, the density value from summing over the bumps
is low, because only few bumps contribute to the density estimate.

####  2.1.1 Stationarity Test

A time series is considered stationary when statistical properties such as
mean, variance, and covariance are constant over time. For making predictions,
the stationarity is a preferred characteristic, to the other side non-
stationarity cannot exploit valuable time-dependent models, and variance can
be misspecified by them.

To test whether our dataset is composed of a stationary time series, we used
the stationary unit root as a statistical test. Given the time series
$y_{t}=a*y_{t-1}+\epsilon_{t}$, where $y_{t}$ is the value at instant $t$ and
$\epsilon_{t}$ is the error term. We evaluate, for all observations
$y_{t}=a^{p}*y_{t-p}+\sum\epsilon_{t-i}*a^{i}$, if the value of $a$ is 1
(unity), then the predictions will equal $y_{t-p}$ and the sum of all errors
from $t-p$ to $t$, meaning that the variance increases over time.

To verify the unit root stationary tests, we used the Augmented Dickey Fuller
(ADF) test by analyzing the test statistic, p-value, and critical value at 1%,
5%, and 10% confidence intervals.

Table 1: ADF test stationarity with AIC optimization. | Test Statistic | p-value | Lags | Observations | Critical Value
---|---|---|---|---|---
|  |  |  |  | 1% | 5% | 10%
ANF | -2.302 | 0.171 | 5 | 2529 | -3.432 | -2.863 | -2.567

Tab. 1 shows the number of lags considered, automatically selected based on
the Akaike Information Criterion (AIC) (see Akaike (1974)) on ANF stock
prices. The p-value result above the threshold (such as 5% or 1%) suggests
rejecting the null hypothesis, so the time series turns out not to be
stationary.

##  3 Methodology

In this paper, ARIMA model as a classical method and ANN as a deep learning
model are chosen for ANF stocks price prediction. In this section, we analyze
the models, comparing their results and selecting the best model based on the
performance measures.

###  3.1 Classic Method

A statistical model is autoregressive (AR) if it predicts future values based
on past values. By the term lag window, we denote the set of previous
observations (lags) used to estimate the value at time $t$. For example, an
autoregressive model of order $p$ ($AR(p)$), uses previous observations $p$ to
predict the current value, so its lag window contains
$y_{t-1},\ldots,y_{t-p}$.

The Auto Correlation Function (ACF) shows how each observation is correlated
$\phi_{p}=corr(y_{t},y_{t-p})$ with its previous one, by measuring the linear
relationship between observations at time $t$ and observations at the previous
time $t-p$.

$y_{t}=\phi_{1}y_{t-1}+\dots+\phi_{p}y_{t-p}$

In our case, we can see from the ACF plot in Fig. 5, the plot of lag values on
the $x$-axis and how they increase, and at the same time the correlation
between price and lagged price on the $y$-axis deteriorates to $lag=750$. This
means that the values of the historical series are strongly correlated with
those of the lagged series for only an initial period, then the correlation
decreases faster and faster (see fig. 6(a)).

Figure 5: ANF Values Correlation.

Although ACF shows that there is no high correlation between $y_{t-p}$ and
$y_{t}$, it can happen that this correlation depends on the fact that both
delays can be correlated with some intermediate delay and their correlation
depends on that. To avoid this, we made further considerations by also
evaluating the PACF that computes the correlation between $y_{t-p}$ and
$y_{t}$ net of the intermediate variables (i.e., excluding correlations with
them by $y_{t-p}$ and $y_{t}$).

Figure 6: ACF and Partial ACF of ANF price.

Thus, if we transform the time series by filtering out the influence of
observations at times $t-1\ldots t-(p-1)$, leaving only the association
between $y_{t}$ and $y_{t-p}$, then the ACF becomes the Partial Auto
Correlation Function (PACF), which denotes by $\alpha(p)$ a very small
correlation (see fig. 6(b)).

Since there is no specific trend in these plots, it can be assumed that the
model was selected correctly and can be used for prediction.

Statistical models combined with machine learning have become ubiquitous in
modern life.

The well-known linear regression (LR) model has various forms, such as the
Autoregressive (AR) model, the moving average (MA) model, the autoregressive
moving average (ARMA) model, and its evolution, the autoregressive integrated
moving average (ARIMA) model (Marquez (1995)).

As a benchmark from a statistical point of view, we used the Auto-ARIMA model
by pmdarima 666https://alkaline-ml.com/pmdarima/setup.html library. In the
basic ARIMA model, we need to provide the $p$, $d$, and $q$ values using
statistical techniques by performing the difference to eliminate the non-
stationarity and plotting ACF and PACF graphs. In Auto-ARIMA, the model itself
automatically discovers the optimal order for an ARIMA model by conducting
differentiation tests (i.e., Kwiatkowski-Phillips-Schmidt-Shin, Augmented
Dickey-Fuller, or Phillips-Perron) to determine the order of differentiation,
$d$, and then fitting the models within defined ranges of beginning $p$,
maximum $p$, beginning $q$, maximum $q$. If the seasonal option is enabled,
Auto-ARIMA also attempts to identify the optimal hyper-parameters $P$ and $Q$
after conducting the Canova-Hansen (Kurozumi (2002)) to determine the optimal
order of seasonal differentiation, $D$.

Tab. 2 shows that the best model is ARIMA(0,1,0) by performing the stepwise
search to minimize AIC (Stoica and Selen (2004)). Auto-ARIMA apply only one
order of differenciation chosing the simple ARIMA called Random Walking (see
Danyliv et al. (2019)) with $p=0$, $d=1$, and $q=0$.

Table 2: ARIMA (p,d,q) models performing stepwise search to minimize Akaike Information Criterion. ARIMA | AIC |
---|---|---
(2,1,2) | 6831.648 |
(0,1,0) | 6830.378 |
(1,1,0) | 6829.872 |
(0,1,1) | 6829.776 |
(0,1,0) | 6828.926 |
(1,1,1) | 6830.470 |

In fig. 7 is shown the prediction made with the Auto ARIMA model to predict
the ANF trend of the $n=30$ days following the training date. The training
period falls in the time interval of 2507 days of market opening from October
30, 2011 to October 16, 2021.

Respect to the DNN used in the next section 4, ARIMA predictions are
completely out of measure, as shown Tab. 4, relatively to the closing price.
It not intercept any type of trend in the interval considered, as specified in
fig. 8, and with detail in Tab. 3).

For the sake of brevity we omit the graphs of Low, High, and Open prices which
are also far from the correct values. For this reason the next step was to
turn to more refined models of machine learning, such as Deep Neural Networks
(DNN).

Figure 7: ARIMA forecast of the last 30 days of ANF stock with p=0, d=1, q=0. Table 3: ARIMA forecast of the last 30 days of ANF stock price trend. ARIMA | mean | std | min | max
---|---|---|---|---
(p=0, d=1, q=0) | 38.94 | 0.123 | 38.73 | 39.14
Table 4: Error metrics of ARIMA on ANF stock price prediction. ARIMA | MSE | RMSE | MAE | MAPE | EVS
---|---|---|---|---|---
(p=0, d=1, q=0) | 25.43 | 5.04 | 3.90 | 0.09 | -0.03

Figure 8: Detail of ARIMA forecast of the last 30 days of ANF stock price.

###  3.2 Deep Learning Method

Deep neural networks are a special type of artificial neural network
characterized by an architecture consisting of a greater number of neurons and
hidden layers than a conventional neural network can be. For each added hidden
layer, in principle, the network will have a greater ability to extract high-
level features.

DNNs have proven to be very efficient at solving nonlinear problems, including
real-world problems. This is in contrast to many traditional techniques for
time series forecasting, such as ARIMA, which assume that the series are
generated by linear processes and consequently may be inappropriate for most
real-world problems that are overwhelmingly nonlinear. While artificial neural
networks (ANNs), as a soft computing technique, are the most accurate and
widely used as prediction models in many areas, including social, engineering,
economic, business, financial, foreign exchange and stock market problems.

An important work comparing the action prediction performance of ANN models
versus ARIMA models was conducted by Yao et al. (1999), they showed that ANN
model performed better than conventional ARIMA models. Equally decisive work
was carried out by Hansen et al. (1999) they compared the prediction
capabilities of ANNs and ARIMA focusing on time series prediction, they showed
also that ANNs outperformed ARIMA in predicting the direction of stock
movement as the latter was able to detect patterns hidden in the data used.

From empirical results obtained, it was shown that ANN model is superior to
ARIMA model and to these we add our work. Much other literature has
demonstrated the prevalent use of ANNs as an effective tool for stock price
prediction including Giordano et al. (2007). This makes ANNs a promising
technique or potential hybrid for predicting time series movement.

In this section, we maintain the same objective of forecasting OHLC values of
the $n=30$ days following the training date, as already done with the ARIMA
statistical model (see sect. 3.1). In this case, however, observations are
grouped in batches of 5 days (batch-size $bs$), corresponding to the working
days in an open market week.

When building a neural network suitable for financial applications, the main
consideration is to find a compromise between generalisation and convergence.
Of particular importance is to try not to have too many nodes in the hidden
layer, as this could lead the DNN to learn only by not acquiring the ability
to generalise.

The prediction method that determines the number of input neurons to DNNs is
based on the lag windowing criterion. This is a technique that consists in
windowing the autocorrelation coefficients before estimating the linear
prediction coefficients (LPC).

For the generation of the feature target, we therefore used a lag window
($t=5$ days) identical to the optimal one found with the Auto ARIMA
statistical model to maintain consistency for comparison. The size of the
number of input neurons will be 5 for each of the four DNNs used on OHLC
prices, and all with identical architecture. Actually, a very important part
of this work, has been the identification of the best network geometry that
allows the lowest error when applied on all 4 DNNs, which has been obtained by
making a python module that generates different network geometries in
combination with the sklearn GridSearchCV777https://scikit-
learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html,
for the setting of the hyper-parameters, propagated on all the DNNs.
Specifically, following the input layer, we added 2 hidden layers composed of
$n*t$ neurons in the first inner layer, reducing to $n*bs$ neurons in the
second inner layer as in Letteri et al. (2018) and in Letteri et al. (2019a).
Each neural network produces predictions of the price of the two considered
stocks over the next $n$ days (see fig. 9) via a single neuron in the output
layer.

Figure 9: Feed-forward DNN architecture.

In Figure 9 is shown, an example of the network architecture used to which we
applied a dropout of 0.2% on each of the two internal layers to reduce
overfitting (Hinton et al. (2012)). In addition, to introduce non-linearity
between layers, we used ReLU as the activation function for overall neurons,
which turns out to be better than a tanh or sigmoid (Krizhevsky et al.
(2012)), despite the fact that the depth of the network consists of only a few
internal layers.

We applied an extension of stochastic gradient descent (SGD) called Adaptive
moment (Adam) as an optimization algorithm to minimize the Loss function, in
combination with the Back Propagation (BP) automatic differentiation algorithm
to compute gradients of DNN weights, with respect to a training dataset.

The loss function applied is the Mean Absolute Error (MAE), also called L1
Loss, since it is considered more robust to outliers, and useful in our case,
since the price fluctuation of the asset ANF has some outliers due to its high
volatility (see fig. 3(c)).

Formally speaking, the MFA calculates the average of the sum of the absolute
differences between the actual and predicted values as follows:
$L1(y,\hat{y})=|y-\hat{y}|$, where $y$ represents the true value and $\hat{y}$
the predicted value. MAE checks the size of the errors in a set of predicted
values, without worrying about their positive or negative direction. If the
absolute values of the errors are not used, then the negative values may
cancel out the positive values.

All DNNs trained can be downloaded from our
repository888https://www.ivanletteri.it/dnnModels/.

###  3.3 Multi-layer Perceptron Regressor Training

The strategies of a trading system tell the investor when to buy or sell
shares in such a way that the sequence of these operations is profitable.

The trading system we propose makes forecasts of OHLC share prices for the $30
open market days following the training date ending on October 16, 2021. The
system’s trading rules follow the following logic: when the expected closing
price is higher than the current opening price, it is time to buy, and vice
versa, to sell all the shares (when they are in the portfolio) at the moment
when the expected closing price is lower. We have developed similar criteria
considering the forecasts of the last month with the trained DNNs by defining
a trading plan with a set of entry and exit rules.

We developed similar criteria by considering the forecasts of the last month
with the trained DNNs defining a trading plan with a set of rules entry and
exit.

##  4 DNN Results

The evaluation is conducted through two components. The first one forecasts
the OHLC prices of the asset, the second one calculates the performance of the
trading system.

###  4.1 Price Action Forecasting

The prediction procedure, of the 30 days of open market, is done by
forecasting one day at a time starting from the 5 days known previous days, as
shown in fig. 10. The single neuron makes a single step forward, the network
asks itself what the stock will do tomorrow, and this it does overall the next
30 days ahead using each time the data of the previous day.

Figure 10: DNN Input/Output.

Of continuation we bring back the graphs on the Price Action, in the 30 days
successive the training, from which it is possible to compare the effective
price regarding the previewed price from the DNNs.

Figure 11: DNNs prediction of OHLC prices on validation set.

####  4.1.1 Error Metrics

The evaluation metrics used for the predictive models are:

  * •

Mean Square Error (MSE): is calculated as the mean of the squared forecast
error values. The square root of the forecast error values dictates that they
are positive and induces the effect of giving greater weight to larger errors.
When the forecast errors are very large or outliers are squared, which has the
effect of dragging the mean of the squared forecast errors outwards, this
converges to a larger mean square error score. The MSE tends to give a worse
score to models that make large forecast errors. The MSE is expressed as
follows:

$MSE(y,\hat{y})=\frac{1}{N}\sum_{t=1}^{N}(y_{t}-\hat{y_{t}})^{2}$

  * •

Root Mean Square Error (RMSE): is the standard deviation of the residuals
(prediction errors). Residuals are a measure of how far from the regression
line data points are. RMSE measures of how spread out these residuals are as
follow:

$RMSE(y,\hat{y})=\sqrt{\sum_{t=1}^{N}(y_{t}-\hat{y_{t}})^{2}}.$

  * •

Mean Absolute Error (MAE): is a measure of the errors between pairs of
observations expressing the same phenomenon, calculated as:

$MAE(y,\hat{y})=\frac{1}{N}\sum_{t=1}^{N}|(y_{t}-\hat{y_{t}})|$

  * •

Mean Absolute Percentage Error (MAPE): measures this accuracy as a percentage
and can be calculated as the average absolute percentage error for each time
period minus the actual values divided by the actual values as follows:

$MAPE(y,\hat{y})=\frac{1}{N}\sum_{t=1}^{N}\big{|}\frac{y_{t}-\hat{y_{t}}}{y_{t}}\big{|}$

where: $n$ is the number of fitted points, $y$ is the actual value, $\hat{y}$
is the forecast value. $\sum$ is summation notation (the absolute value is
summed for every forecasted point in time). MAPE is the most common measure
used to forecast error, and works best if there are no extremes to the data
(and no zeros).

  * •

Explained Variance Score (EVS): explains the dispersion of errors of a given
dataset as follows:

$EVS(y,\hat{y})=1-\frac{Var(y-\hat{y})}{Var(y)}$

where $Var(y-\hat{y})$, and $Var(y)$ is the variance of prediction errors and
actual values respectively. Scores close to 1.0 are highly desired, indicating
better squares of standard deviations of errors.

Predictive accuracy of all four neural network model evaluated as the
coefficient determination i.e. MSE, RMSE, MAE, MAPE. All results reported are
for the test set. The predictive estimation results are summarized in Table 5.

Table 5: Error metrics of ANF stock price prediction of validation set. DNN | MSE | RMSE | MAE | MAPE | EVS
---|---|---|---|---|---
Open | 0.82 | 0.91 | 0.77 | 0.02 | 0.94
High | 1.14 | 0.08 | 0.90 | 0.02 | 0.93
Low | 0.93 | 0.97 | 0.83 | 0.02 | 0.93
Close | 1.75 | 1.32 | 1.07 | 0.02 | 0.91

##  5 Backtracking Trading System

Given a trading strategy that defines entry and exit points, backtesting is
the process of evaluating the performance of a trading strategy based on
historical data by analysing price charts to find the gross return, recording
all trades, etc…

###  5.1 Profit and Risk Metrics

In this section, we set out the metrics related to downside risk by estimating
the potential loss in value of the stock.

  * •

Maximum drawdown (MDD)999https://www.investopedia.com/terms/m/maximum-
drawdown-mdd.asp which measures the largest decline from the peak in the whole
trading period to show the worst case as follow:
$MDD=max_{\tau\in(0,t)}[max_{t\in(0,\tau)}\frac{n_{t}-n_{\tau}}{n_{t}}]$.

  * •

Sharpe ratio (SR)101010https://www.investopedia.com/terms/s/sharperatio.asp as
the risk-adjusted profit measure, which refers to the return per unit of
deviation as follow: $SR=\frac{\mathbb{E}[r]}{[r]}$.

  * •

Sortino ratio (SoR)111111https://www.investopedia.com/terms/s/sortinoratio.asp
is a variant of risk-adjusted profit measure, which applies DD as risk
measure: $SoR=\frac{\mathbb{E}[r]}{DD}$.

  * •

Calmar ratio (CR)121212https://www.investopedia.com/terms/c/calmarratio.asp is
another variant of risk-adjusted profit measure, which applies MDD as risk
measure: $CR=\frac{\mathbb{E}[r]}{MDD}$.

To check the goodness of trades, we mainly focused on the Total Returns
$R_{k}(t)$ for each ANF stocks $(k=1,...,p)$ in the time interval
$(t=1,...,n)$

$R_{k}(t)=\frac{Z_{k}(t+\Delta t)-Z_{k}(t)}{Z_{k}(t)}$ (1)

and furthermore analysing the standardized returns
$r_{k}=(R_{k}-\mu_{k})/\sigma_{k},$ with $(k=1,...,p)$, where $\sigma_{k}$ is
the standard deviation of $R_{k}$, e $\mu_{k}$ denote the average overtime for
the studied period.

###  5.2 Triple EMA strategy for Entry and Exit rules

For the selection of trading strategies we tested a set of 12 technical
indicators, three for each of the following areas: trend, momentum,
volatility, and volume indicators. The best of them all turned out to be the
Triple EMA. The Triple EMA strtegy is a trend following system, based of the
Exponential Moving Average (EMA)
$EMA_{n_{t}}=EMA_{t-1}+(\frac{2}{n+1})[x_{t}^{(ohlc)}-EMA_{t-1}]$.

For the sake of brevity we omit the results of the other 11 indicators,
obtained after the backtracking tests, reporting only the best (Triple EMA).
The triple EMA (Tsantekidis et al. (2017)) is generally used to make short-
term and medium-term intraday trading decisions to enter long or short trades,
depending on bullish or bearish signals. Below is the formula for the triple
exponential moving average:

$TripleEMA_{n_{t}}=(3*EMA_{n_{t}})-(3*EMA(EMA_{n_{t}}))+EMA(EMA(EMA_{n_{t}}))$

Determined the Triple EMA based on the expected OHLC with the DNNs, we define
the entry/exit rules as follows:

  * •

Entry rule
if $(x_{n_{t}}^{(l)}<TripleEMA_{n_{t}}^{(l)})$ or
$(x_{n_{t}}^{(h)}<TripleEMA_{n_{t}}^{(h)})$ and
$(x_{n_{t}}^{(c)}<TripleEMA_{n_{t}}^{(c)})$ or
$(x_{n_{t}}^{(o)}<TripleEMA_{n_{t}}^{(o)})$ is True.

  * •

Exit rule
if$(x_{n_{t}}^{(l)}>TripleEMA_{n_{t}}^{(l)})$ or
$(x_{n_{t}}^{(h)}>TripleEMA_{n_{t}}^{(h)})$ and
$(x_{n_{t}}^{(c)}>TripleEMA_{n_{t}}^{(c)})$ or
$(x_{n_{t}}^{(o)}>TripleEMA_{n_{t}}^{(o)})$ is True.

In this paper, we have used a backtracking system to estimate: the total
profit from a budget of $100 reinvested in compounded mode by calculating the
Expectancy Ratio (ER) as follows: $ER=WR*\mu_{win}-(1-WR)*|\mu_{loss}|$, where
$WR$ is the win rate that evaluates the proportion of the trading period with
positive profit among all trading intervals. The $\mu_{win}$ is the average of
profitable trades and $\mu_{loss}$ is the average of the break
even131313https://www.investopedia.com/terms/b/breakevenanalysis.asp. All of
these metrics, with details about Sharp, Sortino, and Calmar ratios, are shown
in the tab. 6.

Table 6: Performance of trading system using different entry rules. Validation set |  | Num. of
---
Trades
| Total
---
Return ($)
| Expectancy
---
Ratio
| Sharpe
---
Ratio
| Sortino
---
Ratio
| Calmar
---
Ratio
actual | 4 | 4.220 | 2.022 | 1.280 | 2.220 | 5.459
predicted | 3 | 6.126 | 2.112 | 2.194 | 3.340 | 12.403

From the tab. 6, we can see that the system, applying the Triple EMA to the
OHLC price with real values, makes 4 trades (see fig. 12(a)), while using the
Triple EMA on the OHLC prices predicted by the DNN the system chooses to
reduce the number of trades to three. The latter is more effective, resulting
in a higher total profit.

Figure 12: Performance results on real OHLC values on the left and predicted
OHLC values to the right.

In detail, we have that on 3 trades (see fig. 12(b)) a Win Rate of 66.6667%, a
Best Trade 12.611%, Worst Trade -6.5578 %, and a Avg. Trade 2.37017%. We also
verified that the strategy wins on the classic Buy & Hold quote, which has a
negative return of -0.904755%.

After determining the best trades, based on the Triple EMA, we applied them to
the forecast obtained with the DNNs. Subsequently, these same trades were
applied to the actual data and calculated the Cumulative Total Profit.

Figure 13: Trades predicted by DNN mapped on real OHLC price of the validation
set.

In fig. 13, we have mapped the 3 trades onto the actual values of the
validation set and reported in tab. 7 the dates of the trades to the market
price in the period of the trades. In particular, but the sum of the profits
and losses deriving from the spread of each trade allowed us to obtain a total
return of 3.2% in one month on an invested capital of $100.

Table 7: Details of trades. Trades |  | entry
---
[ date | price($)]
| exit
---
[ date | price($)]
Profit ($)
1 | 1st Nov. | 40.46 | 5th Nov. | 45.74 | 5.28
2 | 11th Nov. | 43.36 | 19th Nov. | 44.13 | 0.77
3 | 22th Nov. | 44.75 | 23th Nov. | 41.90 | -2.85

##  6 Conclusion

In this paper, we propose a stock market trading system that uses four deep
neural networks as part of its main components. After a forecasting activity
with the regressive Auto ARIMA model, which proved inconclusive, we used DNNs
with OHLC prices from historical data to forecast the trend of the ANF stock
in the 30 days of open market, following the training date. These four ML
models are incorporated into the trading system to make automated buy and sell
decisions using a robot advisor.

For the sake of simplicity, we have omitted experiments conducted by reducing
and increasing the number of days to predict, because they showed that
predictions further into the future produced less accurate results, whereas 30
days is a good compromise and on average gives the trading system sufficient
time to identify a trade.

The consequences of changing the trading rules were also studied by testing
multiple technical indicators and combinations of them. The best one turned
out to be the Triple EMA, which allowed us to make a profit of $3.2 in one
month by applying the trades obtained from crossing the 3-period Triple EMA
with the OHLC prices predicted by the DNNs, starting from a budget of only
$100.

In the near future we plan to explore additional technical indicators and more
exploratory data analysis following a proven process (see Letteri et al.
(2020a)). The new datasets will certainly require a balancing of buy, sell and
hold trades. This balancing will take place according to the algorithm called
Generative 1 Nearest Neighbours (Letteri et al. (2020c)) and targeted
oversampling (Letteri et al. (2021)).

Furthermore, we plan to study other types of deep learning algorithms by
combining RNN and CNN. We also plan to apply better trading rules that can use
more advanced methods of assigning stop limits and take profits. XAI criteria
inherent to feature selection as in Letteri et al. (2020c); Dyoub et al.
(2020) will be established.

Artificial Intelligence (AI) advisors in general, and AI buy & sell advisors
in particular, introduce ethical issues and concerns. There behaviours should
be monitored for ethical violation (Dyoub et al. (2021a), Dyoub et al.
(2021b)). As another future work, we plan to consider the ethical aspects
related to AI buy & sell advisors.

Given the enormous potential and thanks also to the extendability and
robustness of this trading system, we will further test its reliability on
other stock markets using other inputs to the dataset such as cryptocurrencies
and defi-tokens, also varying the timeframes for day trading and scalping
activities.

## References

  * Letteri et al. [2020a] Ivan Letteri, Giuseppe Della Penna, Luca Di Vita, and Maria Teresa Grifa.  Mta-kdd’19: A dataset for malware traffic detection.  2597:153–165, 2020a.  URL http://ceur-ws.org/Vol-2597/paper-14.pdf.
  * Letteri et al. [2018] Ivan Letteri, Giuseppe Della Penna, and Giovanni De Gasperis.  Botnet detection in software defined networks by deep learning techniques.  11161:49–62, 2018.  doi:10.1007/978-3-030-01689-0_4.  URL https://doi.org/10.1007/978-3-030-01689-0_4.
  * Letteri et al. [2019a] Ivan Letteri, Giuseppe Della Penna, and Giovanni De Gasperis.  Security in the internet of things: botnet detection in software-defined networks by deep learning techniques.  _Int. J. High Perform. Comput. Netw._ , 15(3/4):170–182, 2019a.  doi:10.1504/IJHPCN.2019.106095.  URL https://doi.org/10.1504/IJHPCN.2019.106095.
  * Letteri et al. [2019b] Ivan Letteri, Giuseppe Della Penna, and Pasquale Caianiello.  Feature selection strategies for HTTP botnet traffic detection.  pages 202–210, 2019b.  doi:10.1109/EuroSPW.2019.00029.  URL https://doi.org/10.1109/EuroSPW.2019.00029.
  * Yetis et al. [2014] Yunus Yetis, Halid Kaplan, and Mo M. Jamshidi.  Stock market prediction by using artificial neural network.  _2014 World Automation Congress (WAC)_ , pages 718–722, 2014.
  * Soniya et al. [2015] Soniya, Sandeep Paul, and Lotika Singh.  A review on advances in deep learning.  _2015 IEEE Workshop on Computational Intelligence: Theories, Applications and Future Directions (WCI)_ , pages 1–6, 2015.
  * Day and Lee [2016] Min-Yuh Day and Chia-Chou Lee.  Deep learning for financial sentiment analysis on finance news providers.  _2016 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM)_ , pages 1127–1134, 2016.
  * Letteri et al. [2020b] Ivan Letteri, Antonio Di Cecco, Abeer Dyoub, and Giuseppe Della Penna.  A novel resampling technique for imbalanced dataset optimization.  _CoRR_ , abs/2012.15231, 2020b.  URL https://arxiv.org/abs/2012.15231.
  * Wang et al. [2021] Huiwen Wang, Wenyang Huang, and Shanshan Wang.  Forecasting open-high-low-close data contained in candlestick chart, 2021\.
  * Black [1976] F. Black.  Studies of stock price volatility changes.  3002:177–181, 1976.
  * Chen [2017] Yen-Chi Chen.  A tutorial on kernel density estimation and recent advances, 2017.
  * Akaike [1974] H. Akaike.  A new look at the statistical model identification.  _IEEE Transactions on Automatic Control_ , 19(6):716–723, 1974.  doi:10.1109/TAC.1974.1100705.
  * Marquez [1995] Jamie Marquez.  Time series analysis : James d. hamilton, 1994, (princeton university press, princeton, nj).  _International Journal of Forecasting_ , 11(3):494–495, September 1995.
  * Kurozumi [2002] Eiji Kurozumi.  The limiting properties of the canova and hansen test under local alternatives.  _Econometric Theory_ , 18(5):1197–1220, 2002\.  ISSN 02664666, 14694360.  URL http://www.jstor.org/stable/3533370.
  * Stoica and Selen [2004] P. Stoica and Y. Selen.  Model-order selection: a review of information criterion rules.  _IEEE Signal Processing Magazine_ , 21(4):36–47, 2004.  doi:10.1109/MSP.2004.1311138.
  * Danyliv et al. [2019] Oleh Danyliv, Bruce Bland, and Alexandre Argenson.  Random walk model from the point of view of algorithmic trading, 2019\.
  * Yao et al. [1999] Jingtao Yao, Chew Lim Tan, and Hean-Lee Poh.  Neural networks for technical analysis: a study on klci.  _International journal of theoretical and applied finance_ , 2(02):221–241, 1999.
  * Hansen et al. [1999] James V Hansen, James B McDonald, and Ray D Nelson.  Time series prediction with genetic-algorithm designed neural networks: An empirical comparison with modern statistical models.  _Computational Intelligence_ , 15(3):171–184, 1999.
  * Giordano et al. [2007] Francesco Giordano, Michele La Rocca, and Cira Perna.  Forecasting nonlinear time series with neural network sieve bootstrap.  _Computational Statistics & Data Analysis_, 51(8):3871–3884, 2007.  ISSN 0167-9473.  doi:https://doi.org/10.1016/j.csda.2006.03.003.  URL https://www.sciencedirect.com/science/article/pii/S0167947306000740.
  * Hinton et al. [2012] Geoffrey E. Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  _CoRR_ , abs/1207.0580, 2012.  URL http://arxiv.org/abs/1207.0580.
  * Krizhevsky et al. [2012] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  Imagenet classification with deep convolutional neural networks.  In _Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1_ , NIPS’12, page 1097–1105, Red Hook, NY, USA, 2012. Curran Associates Inc.
  * Tsantekidis et al. [2017] Avraam Tsantekidis, Nikolaos Passalis, Anastasios Tefas, Juho Kanniainen, Moncef Gabbouj, and Alexandros Iosifidis.  Forecasting stock prices from the limit order book using convolutional neural networks.  In _2017 IEEE 19th Conference on Business Informatics (CBI)_ , volume 01, pages 7–12, 2017.  doi:10.1109/CBI.2017.23.
  * Letteri et al. [2020c] Ivan Letteri, Antonio Di Cecco, and Giuseppe Della Penna.  Dataset optimization strategies for malware traffic detection.  _CoRR_ , abs/2009.11347, 2020c.  URL https://arxiv.org/abs/2009.11347.
  * Letteri et al. [2021] Ivan Letteri, Antonio Di Cecco, Abeer Dyoub, and Giuseppe Della Penna.  Imbalanced dataset optimization with new resampling techniques.  In Kohei Arai, editor, _Intelligent Systems and Applications - Proceedings of the 2021 Intelligent Systems Conference, IntelliSys 2021, Amsterdam, The Netherlands, 2-3 September, 2021, Volume 2_ , volume 295 of _Lecture Notes in Networks and Systems_ , pages 199–215. Springer, 2021.  doi:10.1007/978-3-030-82196-8_15.  URL https://doi.org/10.1007/978-3-030-82196-8_15.
  * Dyoub et al. [2020] Abeer Dyoub, Stefania Costantini, Francesca Alessandra Lisi, and Ivan Letteri.  Logic-based machine learning for transparent ethical agents.  In Francesco Calimeri, Simona Perri, and Ester Zumpano, editors, _Proceedings of the 35th Italian Conference on Computational Logic - CILC 2020, Rende, Italy, October 13-15, 2020_ , volume 2710 of _CEUR Workshop Proceedings_ , pages 169–183. CEUR-WS.org, 2020.
  * Dyoub et al. [2021a] Abeer Dyoub, Stefania Costantini, Francesca A. Lisi, and Ivan Letteri.  Ethical monitoring and evaluation of dialogues with a mas.  3002:158–172, 2021a.
  * Dyoub et al. [2021b] Abeer Dyoub, Stefania Costantini, Ivan Letteri, and Francesca A. Lisi.  A logic-based multi-agent system for ethical monitoring and evaluation of dialogues.  In Andrea Formisano, Yanhong Annie Liu, Bart Bogaerts, Alex Brik, Verónica Dahl, Carmine Dodaro, Paul Fodor, Gian Luca Pozzato, Joost Vennekens, and Neng-Fa Zhou, editors, _Proceedings 37th International Conference on Logic Programming (Technical Communications), ICLP Technical Communications 2021, Porto (virtual event), 20-27th September 2021_ , volume 345 of _EPTCS_ , pages 182–188, 2021b.  doi:10.4204/EPTCS.345.32.  URL https://doi.org/10.4204/EPTCS.345.32.