# Kelvin-Helmholtz instability in strongly coupled dusty plasma with
rotational shear flows and tracer transport

Vikram Dharodi<EMAIL_ADDRESS>Institute for Plasma Research, HBNI, Bhat,
Gandhinagar 382 428, India  Bhavesh Patel  Institute for Plasma Research,
HBNI, Bhat, Gandhinagar 382 428, India  Amita Das  Department of Physics,
Indian Institute of Technology, New Delhi, 110016, Delhi, India

###### Abstract

Kelvin-Helmholtz (KH) instability plays a significant role in transport and
mixing properties of any medium. In this paper, we numerically explore this
instability for a two-dimensional strongly coupled dusty plasma with
rotational shear flows. We study this medium using generalized hydrodynamic
fluid model which treats it as viscoelastic fluid. We consider the specific
cases of rotating vorticity with abrupt radial profiles of rotation. In
particular: single-circulation, and multi-circulation vorticity shell profiles
have been chosen. We observe the KH vortices at each circular interface
between two relative rotating flows along with a pair of ingoing and outgoing
wavefronts of transverse shear waves. Our studies show that due to the
interplay between KH vortices and shear waves in the strongly coupled medium,
the mixing and transport behaviour are much better than inviscid hydrodynamic
fluids. In interests of substantiating the mixing and transport behaviour, the
generalized hydrodynamic fluid model is extended to include the Lagrangian
tracer particles. The numerical dispersion of these tracer particles in a flow
provides an estimate of the diffusion in such a medium. We present the
preliminary observations of tracers distribution (cluster formation) and their
diffusion (mean square displacement) across the medium.

###### pacs:

##  I Introduction

The Kelvin-Helmholtz (KH) instability has been ubiquitously observed in
hydrodynamic fluids drazin1970kelvin ; chandrhd1981hydro , plasmas
horton1987kelvin , geophysical flows van2010deep , and astrophysical
situations foullon2011magnetic . This instability occurs in form of vortices
at the interface between two flows in the presence of a velocity shear. The
interactions between these KH vortices govern the transport processes like
mixing and diffusion smyth2012ocean . The main objectives of this study are to
explore the formation and evolution of KH vortices of a rotating dusty plasma
and quantifying the mixing they generate using the tracer particles
simulation. The formation and evolution of these vortices depend on the shear
and nature of the medium (here the viscoelasticity of the system). The KH
instability for dust flows has also been investigated theoretically
d1990kelvin ; rawat1993kelvin ; pandey2012shear ; banerjee2012shear ;
dolai2016effect  and numerically ashwin2010kelvin ; tiwari2012kelvina ;
tiwari2012kelvinb ; tiwari2014kelvin  as well experimentally luo2001kelvin
for planar sheared flows. The rotating vortex flows have been studied
considerably konopka2000rigid ; klindworth2000laser ; schablinski2014sheared ;
choudhary2020three , but to our knowledge, no prior studies have explicitly
examined the KH instability for such flows, except dharodi2014visco ;
dharodi2020rotating . Dharodi $et~{}al.$ dharodi2014visco  numerically studied
shearless smooth rotating flows to avoid KH destabilization and also
considered sharp rotating flows where KH arises in homogeneous medium. Dharodi
in dharodi2020rotating  has explored these rotating sheared flows in
heterogeneous medium.

A dusty plasma can exist in strong coupling state quite easily because of high
charged dust particles, which is called strongly coupled dusty plasma (SCDP).
The SCDP has been modelled under the formalism of generalized hydrodynamic
(GHD) fluid model. This model treats the SCDP as a VE fluid and characterizes
its VE effects through the coupling strength parameter which is often measured
as the ratio ${\eta/\tau_{m}}$, the coupling parameters $\eta$ and $\tau_{m}$
are the shear viscosity and the Maxwell relaxation time, respectively. Thus,
here, the effects of viscoelasticity on KH vortices are observed by varying
the ratio ${\eta/\tau_{m}}$. We consider the incompressible limit of GHD model
which in addition to the evolution of hydrodynamic KH instability, also
supports transverse shear (TS) waves that propagate at phase velocity
$\sqrt{\eta/\tau_{m}}$ kaw1998low ; dharodi2014visco ; dharodi2016sub . These
propagating shear waves have the same symmetry as that of their source
structure, until there is no boundary effect or no interaction with other
waves or obstacle like vortex. Since our interest is in KH instability of
rotating SCDPs, we consider the specific cases of sharp vorticity patches: (i)
single-circulation, and (ii) multi-circulation vorticity shell profiles. The
single-circulation case has already been somewhat discussed in
dharodi2014visco , in the present paper we explore it in more detail. We
observe the KH instability at each circular interface between two relative
rotating flows in form of small vortices along with a pair of ingoing and
outgoing wavefronts of TS waves. The interactions between interacting KH
vortices and TS waves help the VE fluid in better mixing than standard HD
fluids where the only interactions between the KH vortices take place. To
substantiate this observation, the passive tracers are dispersed throughout
the medium.

In the context of fluid mechanics the tracer transport has been studied
extensively for flow visualization douady1991direct ; fessler1994preferential
with the help of theoretical balkovsky2001intermittent ;
falkovich2004intermittent ; zaichik2003two ; bec2006acceleration ,
computational squires1991preferential ; boivin1998direct ; reade2000numerical
; zhou2001modelling ; yeung2001lagrangian ; yeung2002lagrangian ;
ishihara2002relative ; collins2004reynolds ; biferale2004multifractal ;
chun2005clustering ; biferale2005particle ; cencini2006dynamics  and
experimental ott2000experimental ; la2001fluid ; mordant2001measurement ;
voth2002measurement ; sawford2003conditional  approaches. This technique is
also used in complex fluids (polymers, colloids and biological materials)
mason1997particle ; waigh2005microrheology . An analysis of the separation of
the particle trajectories with the 2D hydrodynamic fluid is also being carried
out by Falkovich et al. falkovich2001particles . In dusty plasma, Schwabe et
al. schwabe2014collective  observed the vortex movements by adding some micro
particles around the void. Here, we consider two kinds of point-like tracer
particles, (i) non-inertial, and (ii) inertial tracer particles. The tracer
dynamics is simulated using a one-way coupled Lagrangian point-particle
approach riley1974diffusion ; mclaughlin1989aerosol  which means the tracers
are affected by the fluid flow, but not vice-versa. In case of multi-
circulation vortex profiles, at intermediate time range, a complete picture of
a turbulent flow is observed which has a collection of small KH vortices and
waves. When the system is left for a very long time, it ultimately settles
down to a single vortex faster than in HD fluid. It is observed that the
relaxing rate of such turbulent medium increases with the increasing coupling
strength.

This paper has been organized as follows. In section II, the GHD model
especially developed for the study of SCDPs is extended to include the
transport of passive Lagrangian tracers. This extended model is referred as
incompressible Generalized Hydrodynamic Tracer Transport (i-GHTT) model.
Section III presents the numerical procedure in order to solve the set of
equations of GHTT model. In Sec. IV, we numerically explore the evolution of
different types of sharp rotating vorticities in VE fluids and quantifying the
mixing they generate using the tracer particles simulations. Finally, Section
V contains the discussion and the conclusions.

##  II Generalized Hydrodynamic Tracer Transport (GHTT) Model

A dusty plasma can be prepared or found as a strongly coupled dusty plasma
(SCDP) rather easily because of high charge on the micron-sized dust
particles. Below the crystallization limit, a SCDP behaves like a viscoelastic
fluid which favors both the incompressible transverse shear modes and the
compressible longitudinal modes. To study such SCDP the generalized
hydrodynamic (GHD) fluid model is found to be quite suitable which takes into
account both types of modes. To scrutinize the effect of transverse modes and
to abate the possible pairing with the longitudinal mode, we consider the
incompressible limit of GHD (i-GHD) model. Thus, i-GHD model represents the
incompressible SCDPs which support transverse modes only. The momentum and
continuity equations for i-GHD of homogeneous strongly coupled dusty plasma
can be written as:

$\left[1+\tau_{m}\left(\frac{\partial}{\partial
t}+\vec{v}_{d}\cdot\nabla\right)\right]\left[{\frac{\partial\vec{v}_{d}}{\partial
t}+\vec{v}_{d}\cdot\nabla\vec{v}_{d}}+\frac{\nabla
P}{n_{d}}-{\nabla\phi_{d}}\right]={\eta}{\nabla^{2}}{\vec{v}_{d}}{,}$ (1)

and

$\nabla\cdot\vec{v}_{d}=0{,}$ (2)

respectively. $n_{d}$ is the number density which is normalized by its
respective equilibrium value $n_{d0}$. The scalar potential $\phi_{d}$ in the
dusty plasma system is normalized by ${{K_{B}}{T_{i}}}/{e}$. The parameters
$e$, $T_{i}$ and $K_{B}$ are the electronic charge, ion temperature and
Boltzmann constant, respectively. The charge fluctuation over each dust grain
has been ignored. The time and length are normalized by inverse of dust plasma
frequency
$\omega^{-1}_{pd}=\left({4\pi(Z_{d}e)^{2}n_{d0}}/{m_{d0}}\right)^{-1/2}$,
plasma Debye length
$\lambda_{d}=\left({K_{B}T_{i}}/{4{\pi}{Z_{d}}{n_{d0}}{e^{2}}}\right)^{1/2}$,
respectively. In incompressible limit the Poisson equation has been replaced
by the quasineutrality condition. The dust fluid velocity $\vec{v}_{d}$ is
normalized by ${\lambda_{d}}{\omega_{pd}}$. The term
${\tau_{m}}({\vec{v}_{d}}{\cdot}{\nabla})$ in the generalized momentum
equation is responsible for introducing the collective behavior in the medium.
When this term becomes zero ($\tau_{m}$=0), the momentum equation becomes
Navier-Stokes equation. In other words, the GHD model turns into a standard
hydrodynamic fluid model. Moreover, the presence of this term conserves the
Gallilean invariance tiwari2012longitudinal .

In our earlier manuscript dharodi2014visco , we proposed an idea of extending
the GHD model by including the passive tracer particles in interest to
estimate of the diffusion in a medium. In order to accomplish this, we
consider two kinds of point-like tracer particles, (i) non-inertial tracers,
and (ii) inertial tracers. The non-inertial tracers follow the flow exactly
while the velocity of inertial tracers differs from flow velocity due to
viscous drag force sapsis2010clustering . The tracer particles dynamics is
simulated using a one-way coupled Lagrangian point-particle approach. One way
coupled particle approach means the tracers do not affect the fluid motion, in
other words the tracers are passive. We also neglect any kind of interaction
between particles and gravity effect on their dynamics. We further assume that
the inertial tracer particles with density $\rho_{p}$ different than the
density $\rho_{d}$ of the fluid. Under these assumptions, the particles are
transported in the flow according to the equations:

$\displaystyle\frac{d{\vec{v}_{pi}}}{dt}$
$\displaystyle=\frac{1}{\tau_{s}}({\vec{v}_{d}}({\vec{r}_{pi}})-{\vec{v}_{pi}}){,}$
(3) $\displaystyle\frac{d{\vec{r}_{pi}}}{dt}$
$\displaystyle={\vec{v}_{pi}}{,}$ (4)

where ${\vec{r}_{pi}}$ and $\vec{v}_{pi}$ are the position and velocity of the
$i$th particle, respectively. ${\vec{v}_{d}}({\vec{r}_{pi}})$ is the dust
fluid velocity at the particle position, ${\vec{r}_{pi}}$ which is obtained by
solving the set of Eqs. (8) and (9). These equations are a simplified
approximation of the Maxey-Riley equations maxey1983equation . Although, the
neglected effects might have significant impacts in real flows, but these
could be incorporated into the future study because even after neglecting
them, the Eq. (3) describe an enough complex system and it is worth studying
to set the foundation for the future research. The particle time-scale
$\tau_{s}={2{r^{2}_{0}}{\rho_{p}}}/{(9{\eta}{\rho_{d}})}$ denotes the response
time of the particles is known as Stokes time guha2008transport . $r_{0}$ is
the radius of a particle.

Although the particles are assumed as point particles, but they do have finite
mass and therefore finite inertia. The ratio of a particle time-scale to a
fluid time-scale is known as Stokes number (St). The effect of particle
inertia is often given by using the $St$ or $\tau_{s}$. On the basis of St or
$\tau_{s}$ or inertia: For low value, the particles are predicted to follow
the fluid flow passively like fluid particles, while at very high value almost
the particles remain unaffected by the medium fluctuations. In-between these
two limits, when the particle and fluid time-scales are comparable, the
particles respond in fast and strong manner to the fluctuations.

The non-inertial tracers follow the fluid flow exactly and can be considered
as attached to fluid surface. They are characterised by their position
${\vec{r}_{pi}}$, and velocity ${\vec{v}_{pi}}={\vec{v}_{d}}({\vec{r}_{pi}})$
that is the dust fluid velocity at their position. Their equation of motion
corresponds to the limit $\tau_{s}\rightarrow$ 0 in the set of Eqs. (3) and
(4) becomes:

$\displaystyle\frac{d{\vec{r}_{pi}}}{dt}={\vec{v}_{d}}({\vec{r}_{pi}}){,}$ (5)

The set of Eqs. (1)-(2)-(3)-(4) and Eqs. (1)-(2)-(5) represent the
viscoelastic model for inertial and non-inertial tracer particles,
respectively, in which the tracers follow the evolution of fluid with time in
Lagrangian way. Both the set of equations would be referred as incompressible
Generalized Hydrodynamic Tracer Transport (i-GHT2 or i-GHTT) model henceforth
in the article. It should be noted that in i-GHTT model, in interest to
include both the compressible longitudinal and incompressible transverse modes
just replace the i-GHD model ((set of Eqs. (1)-(2)) with complete GHD model
(set of Eqs. (5)- (6)-(7) in kaw1998low ), say GHTT model .

###  II.1 Transport Properties

To quantify the average diffusion of tracer particles, the ensemble averaged
mean square displacement of tracers is measured which is associated with the
mixing of the fluid jeon2013anomalous . The MSD is defined as,

$MSD(t)={\frac{1}{N}}{\sum^{N}_{j=1}}{(r_{j}(t)-r_{j}(0))^{2}}{,}$

Here, $r_{j}(0)$ is initial position of $j$th particle at t=0 and $r_{j}(t)$
is position at time $t$. N is the total number of tracers in the ensemble. The
slope of the MSD versus time is proportional to the diffusion coefficient of
tracer particles which in turn is supposed to measure the mixing performance
of the carrier fluid. Thus, the mixing performance of the carrier fluid can be
quantified through the time-MSD slope, the larger slope means the carrier
fluid is a better mixture.

##  III Simulation methodology

For the numerical simulations, first we need to express the model Eq. (1) as
per requirements of simulation software, LCPFCT (Boris $etal.$ boris_book ).
To fulfill these requirements split the Eq. (1) in following two coupled
equations,

$\displaystyle{{\frac{\partial\vec{v}_{d}}{\partial
t}+\vec{v}_{d}\cdot\nabla\vec{v}_{d}}+\frac{\nabla
P}{n_{d}}-\nabla\phi_{d}}={\vec{\psi}}{,}$ (6)
$\displaystyle\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}}}{.}$
(7)

For two-dimensional (2D) studies all the variables are functions of $x$ and
$y$ only. The new introduced quantity ${\vec{\psi}}(x,y)$ in LHS of Eq. (6)
represents the strain induced in the elastic medium by the time-varying
velocity fields. Next, the gradient terms are eliminated by taking the curl of
Eq. (6) which yields an equation for the evolution of the vorticity field. So
the coupled set of Eqs. (6)-(7) has been recast in the following form:

$\frac{\partial{\xi}_{z}}{\partial
t}+\left(\vec{v}_{d}\cdot\vec{\nabla}\right){{\xi}_{z}}={\vec{\nabla}}{\times}{\vec{\psi}}{,}$
(8) $\displaystyle\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}}}{.}$
(9)

Here, ${\xi_{z}}={\nabla}{\times}{\vec{v}_{d}}$ is the vorticity.
${\vec{\xi}}$ is normalised with dust plasma frequency. The LCPFCT software
boris_book  is based on a finite difference method has been used to solve the
coupled set of Eqs. (8) and (9). Taking the curl of relation
${\xi_{z}}={\nabla}{\times}{\vec{v}_{d}}$ and using incompressible condition
given by Eq. (2) i.e. $\nabla\cdot\vec{v}_{d}=0$, we get

$\displaystyle{\nabla}^{2}{\vec{v}_{d}}=-{\nabla}{\times}{\xi_{z}}$ (10)

This velocity-vorticity relation is used to update the fluid velocity at each
time step using FISHPACK swarztrauber1999fishpack . The validation of our
numerical code has been done in our earlier papers dharodi2014visco ;
dharodi2020rotating .

Further details on simulation procedure in advancing the tracer particles with
flow of a VE fluid. We have the dust fluid velocity for each particle at their
respective position at each time step. This dust velocity is going to be used
in particle momentum Eq. (3). Equations (3) and (4) are numerically integrated
together to find the new position and velocity at the end of each time step.
This integration is based on the first order Runge-Kutta method. The particle
velocity $\vec{v}_{p}$ is calculated by interpolating the velocity defined on
nearby grid points, based on first-order Lagrangian interpolation scheme. The
particles are advanced with fluid time step.

##  IV Numerical results and discussion

The prime objectives of this section are to numerically explore the evolution
of different types of sharp rotating vorticity flows in VE fluids and
quantifying the mixing they generate using the tracer particles simulations.
For each type of flow, the simulations are performed for the varying coupling
strength of VE fluid which is usually measured as the ratio ${\eta/\tau_{m}}$.
All the simulations are performed with periodic boundary conditions in both
the x and y directions on simulation box.

###  IV.1 KH instability of rotating SCDPs

Since our interest is in KH instability of rotating SCDPs, we consider the
specific cases of sharp vorticity patches: (i) single-circulation, and (ii)
multi-circulation vorticity shell profiles. Before proceeding with the direct
assessment of the evolution of KH instability through the numerical results,
it is good to understand the process of formation of these vortices for sharp
rotating flows through the schematic picture illustrated in Fig. 1. The multi-
circulation vorticity profile can be depicted as core-shell fluid flows. In
Fig. 1, one of rotating flows forms the inner core (yellow color regime) and
the others make the outer shells (cyan and green color regimes). Each rotating
flow is divided by a sharp interface; the fluid on either side rotates in
opposite directions. The black circle with arrows in the flow regime indicate
its direction of rotation.

Figure 1: Formation of KH vortices at the circular sharp interfaces of core-
shell flow regions (without any corresponding scale). One of flows forms the
inner core and the others make the outer shells. The black circles with arrows
in the flow regions indicates the direction of rotation. At core-shell
interface, the KH votices rotate anti-clockwise (blue color vortices with blue
curved arrows) while at shell-shell interface rotate clockwise (red color
vortices with red curved arrows).

At each of the interfaces, the counter-rotating flows create a region of high
shear which immediately evolve into small co-rotating (like-sign) KH vortices.
The schematic diagram clearly illustrates how the direction of rotation of
these vortices depend on the relative motion between the alternative flows. At
inner (core-shell) interface, the KH votices rotate anti-clockwise (blue color
vortices with blue curved arrows) while at outer (shell-shell) interface
rotate clockwise (red color vortices with red curved arrows). Thus far, these
appraisals are particularly true for an inviscid HD fluid where no source term
exists. Whereas an incompressible VE fluid, besides KH instability, would also
support the TS waves that propagate through the medium at phase velocity
$v_{p}=\sqrt{\eta/\tau_{m}}$. A medium with higher coupling strength (ratio
${\eta/\tau_{m}}$) favors the faster and stronger TS waves. A stronger wave
shows a less fall in amplitude with time in comparison to the lower one (see
Figs. 2 and 4 in dharodi2014visco , and Fig. 4 in dharodi2020rotating ).
Since, the present simulations have been carried out in the x-y plane, which
is the plane of rotation of vorticity profile, the shear waves emitted from
each interface should be cylindrical in shape. Thus, a stronger cylindrical
wave tries to dominate over KH instability and attempts to keep the
cylindrical symmetry during the evolution.

In evolution of these KH vortices the transport processes like merging and
convection become important. When two co-rotating (like-sign) vortices are
brought sufficiently close to each other they start to rotate around one
another and eventually merge to form a single vortex while in convection the
counter-rotating (unlike-sign) propagate together as single structure (dipole)
to convect the fluid.

####  IV.1.1 Evolution of single-circulation sharp vorticity vortex

The velocity profile for the single-circulation sharp vorticity vortex is
given as follows

${\vec{v}_{0}}=\left\\{\begin{array}[]{ll}v_{x0}=-{\phi_{0}}{\frac{(y-y_{c})}{b}};\>v_{y0}={\phi_{0}}{\frac{(x-x_{c})}{a}}&\quad|r|\leq
6\\\ 0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (11)

The vorticity corresponding to the above velocity profile is given below

${\xi_{z0}}=\left\\{\begin{array}[]{ll}{\phi_{0}}{\left(\frac{1}{a}+\frac{1}{b}\right)}&\quad|r|\leq
6.0\\\ 0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (12)

Here
$|r|={\sqrt{{\left({(x-x_{c})/a}\right)}^{2}+{\left({(y-y_{c})/b}\right)}^{2}}}$,
$a$ and $b$ are the major and minor axes, respectively. $x_{c}$ and $y_{c}$
are the $x$ and $y$ coordinates of the center of the vorticity profile. The
vorticity will have a clockwise rotation if amplitude ${\phi_{0}}{>}0$, or
have a anti-clockwise rotation if ${\phi_{0}}{<}0$ and have no rotation if
${\phi_{0}}{=}0$.

We consider the clockwise rotating vorticity profile has circular symmetry
with parameters $a=b=1$, amplitude ${\phi_{0}}=1$, sharp cutoff at radial
distance $|r|$=6 units away from the centre of the vortex $(0,0)$. This
vorticity profile (yellow color) is shown in the following Fig. 2(a) which has
circular interface with surrounding stagnant fluid at $t=0$. The simulation
region is a square box of length 12$\pi$ units.

Figure 2: Initial rotating sharp vorticity profiles at t = 0 to study the KH
instability. The color scale corresponding to a particular vorticity profile
has been given in the following respective figure which shows its evolution.
The simulation region is a square box of length 12$\pi$ units for all these
systems.

The panels of Fig. 3 shows time evolution of vorticity profile given in Fig.
2(a) in the inviscid HD fluid. The sharpness of the clockwise rotating
vorticity profile generates a strong rotational sheared flow. This sheared
flow results in creation of small co-rotating (anti–clockwise) KH vortices
(dark blue color vortices) at the vorticity interface $|r|=6$. These like-sign
vortices start merging as rotation progresses that leads the fluid to evolve
into an anisotropic isolated structure.

Figure 3: An inviscid HD fluid. The time evolution of a sharp clockwise
rotating vorticity which has circular interface with surrounding stagnant
fluid. The sharpness of the vorticity profile generates small KH vortices at
the interface that results in a anisotropic isolated structure.

In Fig. 4 where ${\eta}=5$ and ${\tau_{m}}=20$, once the vortex begins to
rotate, we observe a pair of ingoing and outgoing cylindrical shear waves from
the interface along with these small like-sign KH vortices. During the
evolution, it is observe that both the waves carry the like-sign vortices
which interact with themselves in order to merge and have interplay with these
waves as well. The fluid within the inner region ($|r|\leq 6$) favors mixing
due to the ingoing waves, while the stationary fluid in the outer region
($|r|\geq 6$) gets mixing due to the outgoing waves. Thus, the TS waves assist
the process of fluid mixing by convecting it inside and outside the vortex
structure. Unlike the inviscid case, here, the KH vortices are confined to the
radial emitted waves. Since the wavefronts are cylindrical in shapes, the
interplay between waves and vortices leads the evolution of the VE fluid
toward an isotropic structure.

Figure 4: Viscoelastic fluid with coupling parameters ${\eta}=5$, and
${\tau_{m}}=20$. The time evolution of a sharp rotating vorticity which has
circular interface with surrounding stagnant fluid. The sharpness of the
vorticity profile generates small KH vortices at the interface along with a
pair of ingoing and outgoing wavefronts of TS waves that assist in fluid
mixing by convecting it inside and outside the vortex structure.

In Fig. 5, we have simulated another case of VE fluid which has less coupling
strength ($\eta/\tau_{m}$=0.125) with ${\eta}=2.5$ and ${\tau_{m}}$=20. Unlike
Fig. 4, since the TS waves are not strong enough that they can dominant over
the KH instability, less confinement of KH vortices to the waves. This is
evident from Fig. 5 where the evolution of medium towards an isotropic
structure and mixing process are much slower in comparison to Fig. 4.

Figure 5: The time evolution of a sharp rotating vorticity in VE fluid with
coupling parameters ${\eta}=2.5$, and ${\tau_{m}}=20$. Due to the less
coupling strength, the evolution of medium towards an isotropic structure and
mixing process are much slower in comparison to Fig. 4.

In conclusion, as a result of greater the coupling strength or stronger TS
wave, the medium evolution attempts to possess radial symmetry and shows
better mixing. The mixing is found minimal in inviscid fluid.

The single-circulation sharp profile carries single types of anti-clockwise
rotating KH vortices across its only interface. While a multi-circulation
profile with two or more than two interfaces produces clockwise as well as
anti-clockwise KH vortices. Thus, in multi-circulation case apart from merging
between like-sign vortices, the propagation of unlike-sign vortices as a
dipolar structure also becomes important which can assist to increasing the
spatial domain of mixing fluids.

####  IV.1.2 Multi-circulation vorticity shell profile

Here, we first consider the simplest case of multiple shells of vorticity
having two flows: Inner core flow and a outer shell flow, both flows have
reversal circulation. The velocity profile for this configuration is given
below.

${\vec{v}_{0}}=\left\\{\begin{array}[]{ll}v_{x0}=-{\phi_{0}}{(y-y_{c})};\>v_{y0}={\phi_{0}}{(x-x_{c})}&\quad|r|\leq
5\\\ v_{x0}={\phi_{0}}{(y-y_{c})};\>v_{y0}=-{\phi_{0}}{(x-x_{c})}&\quad
5<|r|\leq 10\\\ 0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (13)

The vorticity corresponding to the above velocity profile is given below,

${\xi_{z0}}=\left\\{\begin{array}[]{ll}2{\phi_{0}}&\quad|r|\leq 5\\\
-2{\phi_{0}}&\quad 5<|r|\leq 10\\\
0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (14)

With the parameter $\phi_{0}$=1, we consider both flows rotate with equal
rotation rates and in opposite directions (Fig. 2(b)). Figure 6 shows the
evolution of this vorticity profile for an inviscid fluid. At inner interface,
the KH vortices rotate anti-clockwise (blue color vortices) while at outer
interface rotate clockwise (red color vortices). As time goes on, the merging
between like-sign vortices take place at both the interfaces, and
simultaneously growing closeness between both the interfaces due to the radial
gradient in vorticity results in interactions between the counter-rotating
(red-blue) vortices as well. These counter-rotating vortices results in the
formation of propagating dipolar structures which help to convect the fluid
across the wider domain than the single interface (Fig. 3) which only favors
the merging process.

Figure 6: An inviscid HD fluid. Evolution of of vorticity having two circular
sharp flows, inner core and an outer shell flows have reversal circulation, in
time generates small KH vortices at each interface that results in a
anisotropic isolated structure.

In Fig. 7, we observe that a pair of ingoing and outgoing wavefronts emanates
from each of the two sharp interfaces of the vortex structure along with KH
vortices (see the first panel). It is observe that both the wavefronts at
inner interface carry the blue color (rotating anti-clockwise) vortices while
at outer interface carry red color (rotating clockwise) vortices. During the
emission of these shear waves like-sign vortices interact with themselves in
order to get merge.

Figure 7: Viscoelastic fluid with coupling parameters ${\eta}=5$, and
${\tau_{m}}=20$. The time evolution of vorticity having two circular sharp
interfaces, inner core and outer shell have reversal circulation, generates
small KH vortices at each interface along with a pair of ingoing and outgoing
TS waves fronts that results in well fluid mixing by convecting it inside and
outside the vortex structure.

The stagnant fluid in the outermost region ($|r|\geq 10$ ) undergoes mixing
due to the outgoing wave from the outermost interface at 10, while the
innermost vortex region undergoes mixing due to the ingoing wave emanating
from the sharp interface located at 5. Interestingly, the vortex region
confined within the two sharp interfaces ($5<|r|\leq 10$) undergoes mixing due
to the ingoing wave from the outermost interface and the outgoing wave from
the innermost interface. This region is probably convection dominating region
due to the higher possibility of formation of propagating dipolar structures.
As the results of multiple interaction processes the mixing becomes fast and
efficient than the cases discussed so far. Next, we consider the another VE
fluid in Fig. 8 having lower coupling strength i.e.
${\eta}=2.5,{\tau_{m}}=20,v_{p}$=0.35, it is observed that at each time step
in the spatial confinement of this medium is lower than Fig. 7. Moreover, the
comparison manifests that the mixing and evolution symmetry of a medium are
proportional to the coupling strength of that medium as observed for earlier
cases.

Figure 8: Viscoelastic fluid with coupling parameters ${\eta}=2.5$, and
${\tau_{m}}=20$. The time evolution of vorticity having two circular sharp
interfaces, both have reversal circulation. Due to the less coupling strength
the spatial confinement of this medium is lower than Fig. 7.

Next, in order to make a better comparative numerical analysis about the
mixing rate, we have considered a more complex scenario of multiple
circulations (each consecutive one having a reversal in its circulation)
having the following velocity flow profile:

${\vec{v}_{0}}=\left\\{\begin{array}[]{ll}v_{x0}=-{\phi_{0}}{(y-y_{c})};\>v_{y0}={\phi_{0}}{(x-x_{c})}&\quad|r|\leq
2.5\\\ v_{x0}={\phi_{0}}{(y-y_{c})};\>v_{y0}=-{\phi_{0}}{(x-x_{c})}&\quad
2.5<|r|\leq 5\\\
v_{x0}=-{\phi_{0}}{(y-y_{c})};\>v_{y0}={\phi_{0}}{(x-x_{c})}&\quad 5<|r|\leq
7.5\\\ v_{x0}={\phi_{0}}{(y-y_{c})};\>v_{y0}=-{\phi_{0}}{(x-x_{c})}&\quad
7.5<|r|\leq 10\\\ 0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (15)

The vorticity corresponding to the above velocity profile is given below,

${\xi_{z0}}=\left\\{\begin{array}[]{ll}2{\phi_{0}}&\quad|r|\leq 2.5\\\
-2{\phi_{0}}&\quad 2.5<|r|\leq 5\\\ 2{\phi_{0}}&\quad 5<|r|\leq 7.5\\\
-2{\phi_{0}}&\quad 5<|r|\leq 10\\\
0&\quad\mathrm{otherwise}{.}\end{array}\right.$ (16)

For the parameter $\phi_{0}$=1, the initial vorticity profile is shown in Fig.
2(c). The complexity of this motion of multi-circulation structure is evident
from the subplots of Fig. 9 for inviscid fluid. In initial time period, the
vortices of KH instability develop across the interface of each shell. At
intermediate time range, this evolution provides a complete picture of a
turbulent flow throughout the entire system which is collection of several
small symmetric and non-symmetric vortices. The transport phenomena like
merging between two co-rotating vortices, convection due to the evolution of
dipolar (two counter-rotating ) vortices, collision between these dipolar
vortices, and also the formation and evolution of triplor structures becomes
more frequent than above discussed cases.

Figure 9: An inviscid HD fluid. The time evolution of multi-circulation
vorticity profile, each consecutive one having a reversal in its circulation,
generates small KH vortices at each interface. This evolution provides a
complete picture of a turbulent flow which is collection of various kind of
vortices, and exhibits transport properties like convection and merging.

Figure 10, represents the evolution of same initial profile (Fig. 9) of
vorticity for VE fluid with coupling parameters ${\eta}=5$, ${\tau_{m}}$=20.
From the comparative observations between Fig. 9 and Fig. 10, it is
interesting to notice that the presence of TS waves leads to the relaxing of
the turbulent medium to a single vortex faster than in inviscid fluid.

Figure 10: The time evolution of multi-circulation vorticity profile in VE
fluid with ${\eta}=5$, and ${\tau_{m}}=20$. The strong interaction between KH
vortices and TS waves results in the relaxing the medium into a single vortex
faster than inviscid fluid.

In Fig. 11 (${\eta}=2.5,{\tau_{m}}=20,v_{p}$=0.35), since the TS waves are
weaker than Fig. 10, the less spatial confinement of the turbulent medium that
results in slower merging process, due to which the relaxation time of the
medium becomes longer.

Figure 11:  The time evolution of multi-circulation vorticity profile in VE
fluid with ${\eta}=2.5$, and ${\tau_{m}}=20$. Since the TS waves are weaker
than Fig. 10, the relaxation time of the medium becomes longer.

The comparison of Fig. 10 and Fig. 11 clearly displays these observations.

###  IV.2 Transport of tracer particles

In Section II we have stated that tracer particles with very low inertia
follow the flow passively, while particles with very high inertia will remain
almost unaffected by the medium fluctuations. In between these two limits
particles show the strongest response to the medium fluctuations. The
simulations are performed for all these three inertial particles: very low
($\tau_{s}$=0.05), intermediate ($\tau_{s}$=1), and very high inertia
($\tau_{s}$=50). To observe the exclusive effect of inertia, we transport
these particles through a similar smooth rotating vorticity vortex in an
inviscid fluid. An inviscid fluid has no source term which favors the emission
of TS waves or dissipative term like viscosity. And also, we choose the smooth
rotating vorticity vortex which does not satisfy the KH destabilization
condition anywhere in the vorticity patch. The equation of such smooth
rotating vorticity is given as

${\xi_{0}}={\Omega_{0}}exp\left(-\left({x^{2}+y^{2}}\right)/{a^{2}_{c}}\right){.}$
(17)

The evolution of same structure given by Eq. 17 for ${a_{c}}$=1.0,
${\Omega_{0}}=5$ in an inviscid fluid is shown in Figs. (12), (13), and (14).
From these figures it is clear that the rotating vortex keeps rotating without
any change.

Figure 12: An inviscid HD fluid carries tracer particles are shown as red
dots. The low inertial tracers ($\tau_{s}$=0.05) follow the dynamics along the
smooth rotating vortex.

Now, in order to see the response of tracer particles, initially $(t=0)$, we
have distributed 900 inertial particles (shown by red dots) homogeneously
throughout the domain.

Figure 13: An inviscid HD fluid carries tracer particles are shown as red
dots. The high inertial tracers ($\tau_{s}$=50) show negligible response to
the rotating vortex.

From Fig. 12, it is clear that low inertial particles ($\tau_{s}$ =0.05)
follow the dynamics along the rotating vortex, and the particles with higher
inertia $\tau_{s}$= 50 (Fig. 13) show negligible response to the vorticity
gradient.

Figure 14: An inviscid HD fluid carries tracer particles are shown by red
dots. Smooth rotating vorticity vortex with tracers having intermediate
inertia i.e $\tau_{s}$=1 show high response to the vorticity gradient, the
particles are pushed away where the flow is strong enough and get accumulate
in strain-dominated region.

In comparison to previous cases (Fig. 12 and Fig. 13), Fig. 14 shows that the
particles with intermediate value of $\tau_{s}$=1 counter a significant
outward push. It is because the particle and fluid time- scales are comparable
which results the particles experience a notable centrifugal force due to
vorticity gradient. Since, the inertial particles are pushed away from regions
where the flow is strong enough, these particles get accumulate in strain-
dominated regions. Thus, particles tend to leave regions of high vorticity and
cluster into regions of high strain ravichandran2017clustering . Note, we have
simulated a range of intermediate inertial particles, and observed the same
effect with varying outward push.

With the identification of intermediate inertial particles and understanding
of their evolution, next we compare the evolution of intermediate inertial
(typical, $\tau_{s}$=1) with non-inertial particles for the inviscid HD and VE
($\eta=5$, $\tau_{m}=20$) fluids in the following Fig. 15 and Fig. 16,
respectively. For this, we choose a sharp rotating vorticity profile which is
given by Eq. (12). The middle row in Fig. 15 and Fig. 16 represent the time
evolution of this profile for inviscid and VE fluids, respectively. As the
sharp vortex starts to rotate, produces larger strain (deformation) in medium
along the interface that results in formation of KH instability. The evolution
of this rotor has already been discussed in the earlier Section IV.1 for
inviscid fluid in Fig. 3 and for VE fluid in Fig. 4. Initially $(t=0)$, here,
we distribute 3600 tracer particles homogeneously throughout these fluids.

Figure 15: An inviscid HD fluid advect tracer particles. First and third row
show the temporal and spatial distribution of non-inertial particles (shown by
magenta dots) and inertial particles ($\tau_{s}$=1, shown by red dots),
respectively, corresponding to the sharp vorticity profile evolution shown in
second row.

In both Figs. 15 and 16, the first and the third row visualize the
distribution of non-inertial particles (shown by magenta dots) and inertial
particles (shown by red dots) respectively, advect by their respective fluid
flow shown in the middle row. The accumulation of non-inertial particles is
observed in rotation-dominated regions over the vortex structures. While the
inertial particles accumulate in strain-dominated regions along the
interfaces.

Figure 16: Viscoelastic fluid with coupling parameters ${\eta}=5$, and
${\tau_{m}}=20$ advect tracer particles. First and third row show the
spatiotemporal distribution of non-inertial particles (shown by magenta dots)
and inertial particles ($\tau_{s}$=1, shown by red dots), respectively,
corresponding to the sharp vorticity profile evolution given in second row.

This accumulation process leads to the spatial inhomogeneous distribution of
particles. This inhomogeneous distribution of the particles is known as
clustering or preferential concentration. Clustering is well-studied in the
case of inertial particles maxey1987gravitational ; squires1991preferential ;
fessler1994preferential ; balkovsky2001intermittent ;
falkovich2002acceleration ; goto2006self ; petersen2019experimental ;
oka2021generalized , and several studies are available for non-inertial
particles drotos2019inhomogeneities .

In real flows the tracer particles always have some inertial value, so we
compute the ensemble averaged MSD of intermediate inertial particles to
analyze the diffusion of particles in the carrier VE fluid. This diffusion of
particles is associated with the mixing of the fluid. For this, we advect the
inertial particles having $\tau_{s}$=1 using the sharp rotating flows
discussed above, inviscid fluid in Fig. 3/Fig. 15, VE fluid with $\eta=5$,
$\tau_{m}=20$ in Fig. 4/Fig. 16, and another VE fluid with $\eta=2.5$,
$\tau_{m}=20$ in Fig. 5. Initially t=0, we disperse 3600 tracer particles
homogeneously over these fluids.

Figure 17: The MSD as function of time for tracing particles. The sharp
rotating vorticity advect the particles with (a) $\tau_{s}$=1, and (b)
$\tau_{s}$=0.5. At later time, the MSD is proportional to the coupling
strength.

Figure 17(a) compares the ensemble averaged MSD values of inertial particles
of $\tau_{s}$=1 for all three types of fluids. The evolution of MSD occurs in
three stages. Initially (here, $0{\geq}t{\approx}1$), inertial particles do
not respond until their time-scale (here, $\tau_{s}=1$) is less than fluid
time-scale or Stokes number (St) is less than unity. During the second stage
the particles start to respond the flow ($1{\geq}t{\leq}3.8$), the slopes grow
at the almost same rate for all fluids. In the final stage ($t{\geq}3.8$), the
slope shows the larger value for the higher coupling strength and it is
minimal for inviscid fluid. The slope of the MSD versus time is proportional
to the diffusion coefficient of the tracer particles. Thus, the larger time-
MSD slope means the carrier fluid shows a better mixing performance. Thus,
this result is found consistent with our earlier observations made from the
comparative analysis of pictorial evolution of vorticities in Figs. 3, 4, and
5. In order to substantiate these observations we further compute the MSD of
tracers with $\tau_{s}$=0.5 (Figure 17(b)) for the all three flows under the
same flow conditions. From the plot, again the diffusion of particles
increasing with coupling strength. It should be noted that the MSD has been
calculated up to time before hitting the flow to the boundaries.

In the present work, although we have undertaken the simplest tracer model and
preliminary investigation of tracers distribution (cluster formation) and
their transport property (MSD), but it is worth attempting as a foundation for
the future study.

##  V Summary and conclusions

In this paper, we numerically explore the KH instability for a two-dimensional
rotating SCDP under the formalism of GHD fluid model. This model treats the
SCDP as a viscoelastic fluid. Here, we consider the specific cases of
vorticity with abrupt radial changes, in particular: single-circulation, and
multi-circulation vorticity shell profiles. We observe the KH instability in
form of small vortices at each interface along with a pair of ingoing and
outgoing wavefronts of TS waves. The interactions between KH vortices, and
with TS waves govern the mixing which has been quantified using the passive
tracer particles simulation. The advection of tracers with the flow has been
noticed. Some main observations are as follows.

  1. 1.

The interplay between the TS waves and interacting KH vortices results in
better mixing of VE fluids than standard hydrodynamic fluids where the only
interaction between KH vortices happens.

  2. 2.

By tweaking the coupling strength parameter (ratio $\eta/\tau_{m}$) which
usually represents the strength of viscoelasticity of the medium, one can
control the evolution of KH instability that results in control over the
transport properties like mixing and diffusion.

  3. 3.

The GHD model especially developed for the study of SCDPs is extended to
include the transport of passive Lagrangian inertial and non-inertial tracer
particles. This extended model is referred as Generalized Hydrodynamic Tracer
Transport (GHTT) model.

  4. 4.

We observe that the diffusion of intermediate inertial tracer particles in VE
fluids is proportional to the coupling strength. It is least for an inviscid
fluid.

  5. 5.

The multi-circulation vortex profiles, we found that the relaxing rate of a
turbulent medium increases with the increasing coupling strength.

The particle tracking model appears a suitable diagnostic to understand the
associated mixing in a fluid through the diffusion process, and nonlinear dust
fluid dynamics through the clustering of a tracers. In present paper, we just
discuss the preliminary simulations for tracers which are based on first order
schemes like to advance the tracers we use first order RK scheme, and for the
velocity the first order interpolation scheme. Further improvisations of this
model would be to include higher order schemes, finite size tracers, other
forces e.g. gravity, inter-particle interactions, feedback from tracers to the
carrier fluid, etc. These improvisations would be extremely fruitful which are
left to future publications.

Although, the rotational dust flows are extensively studied, but to our
knowledge, no prior studies have explicitly examined the KH instability for
the rotating flows, except dharodi2014visco ; dharodi2020rotating . An
experimental research effort is required in order to valid this numerical
work.

## References

  * [1] PG Drazin.  Kelvin-helmholtz instability of finite amplitude.  Journal of Fluid Mechanics, 42(2):321–335, 1970.
  * [2] S Chandrasekhar.  Hydrodynamic and hydromagnetic stability, dover, new york, 1981.  1981\.
  * [3] W Horton, T Tajima, and T Kamimura.  Kelvin–helmholtz instability and vortices in magnetized plasma.  The Physics of fluids, 30(11):3485–3495, 1987.
  * [4] Hans van Haren and Louis Gostiaux.  A deep-ocean kelvin-helmholtz billow train.  Geophysical Research Letters, 37(3), 2010.
  * [5] Claire Foullon, Erwin Verwichte, Valery M Nakariakov, Katariina Nykyri, and Charles J Farrugia.  Magnetic kelvin-helmholtz instability at the sun.  The Astrophysical Journal Letters, 729(1):L8, 2011.
  * [6] William D Smyth and James N Moum.  Ocean mixing by kelvin-helmholtz instability.  Oceanography, 25(2):140–149, 2012.
  * [7] N d’Angelo and Bin Song.  The kelvin-helmholtz instability in dusty plasmas.  Planetary and space science, 38(12):1577–1579, 1990.
  * [8] SPS Rawat and NN Rao.  Kelvin-helmholtz instability driven by sheared dust flow.  Planetary and space science, 41(2):137–140, 1993.
  * [9] BP Pandey, SV Vladimirov, and A Samarian.  Shear driven instabilities in dusty plasmas.  In Europhysics Conference on the Atomic and Molecular Physics of Ionized Gases (21st: 2012), pages 1–2. European Physical Society, 2012.
  * [10] D Banerjee, MS Janaki, and N Chakrabarti.  Shear flow instability in a strongly coupled dusty plasma.  Physical Review E, 85(6):066408, 2012.
  * [11] Bivash Dolai, RP Prajapati, and RK Chhajlani.  Effect of different dust flow velocities on combined kelvin-helmholtz and rayleigh-taylor instabilities in magnetized incompressible dusty fluids.  Physics of Plasmas, 23(11):113704, 2016.
  * [12] J Ashwin and R Ganesh.  Kelvin helmholtz instability in strongly coupled yukawa liquids.  Physical review letters, 104(21):215003, 2010.
  * [13] Sanat Kumar Tiwari, Amita Das, Predhiman Kaw, and Abhijit Sen.  Kelvin–helmholtz instability in a weakly coupled dust fluid.  Physics of Plasmas, 19(2):023703, 2012.
  * [14] Sanat Kumar Tiwari, Amita Das, Dilip Angom, Bhavesh G Patel, and Predhiman Kaw.  Kelvin-helmholtz instability in a strongly coupled dusty plasma medium.  Physics of Plasmas, 19(7):073703, 2012.
  * [15] Sanat Tiwari, Vikram Dharodi, Amita Das, Predhiman Kaw, Abhijit Sen, et al.  Kelvin-helmholtz instability in dusty plasma medium: Fluid and particle approach.  J. Plasma Phys, 80(6):817–823, 2014.
  * [16] QZ Luo, N D’Angelo, and RL Merlino.  The kelvin–helmholtz instability in a plasma with negatively charged dust.  Physics of Plasmas, 8(1):31–35, 2001.
  * [17] U Konopka, D Samsonov, AV Ivlev, J Goree, V Steinberg, and GE Morfill.  Rigid and differential plasma crystal rotation induced by magnetic fields.  Physical Review E, 61(2):1890, 2000.
  * [18] M Klindworth, A Melzer, A Piel, and VA Schweigert.  Laser-excited intershell rotation of finite coulomb clusters in a dusty plasma.  Physical Review B, 61(12):8404, 2000.
  * [19] Jan Schablinski, Dietmar Block, Jan Carstensen, Franko Greiner, and Alexander Piel.  Sheared and unsheared rotation of driven dust clusters.  Physics of Plasmas, 21(7):073701, 2014.
  * [20] Mangilal Choudhary, Roman Bergert, Slobodan Mitic, and Markus H Thoma.  Three-dimensional dusty plasma in a strong magnetic field: Observation of rotating dust tori.  Physics of Plasmas, 27(6):063701, 2020.
  * [21] Vikram Singh Dharodi, Sanat Kumar Tiwari, and Amita Das.  Visco-elastic fluid simulations of coherent structures in strongly coupled dusty plasma medium.  Physics of Plasmas, 21(7):073705, 2014.
  * [22] Vikram S Dharodi.  Rotating vortices in two-dimensional inhomogeneous strongly coupled dusty plasmas: Shear and spiral density waves.  Physical Review E, 102(4):043216, 2020.
  * [23] PK Kaw and A Sen.  Low frequency modes in strongly coupled dusty plasmas.  Physics of Plasmas, 5(10):3552–3559, 1998.
  * [24] Vikram Dharodi, Amita Das, Bhavesh Patel, and Predhiman 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(1):013707, 2016.
  * [25] S Douady, Y Couder, and ME Brachet.  Direct observation of the intermittency of intense vorticity filaments in turbulence.  Physical review letters, 67(8):983, 1991.
  * [26] John R Fessler, Jonathan D Kulick, and John K Eaton.  Preferential concentration of heavy particles in a turbulent channel flow.  Physics of Fluids, 6(11):3742–3749, 1994.
  * [27] E Balkovsky, Gregory Falkovich, and A Fouxon.  Intermittent distribution of inertial particles in turbulent flows.  Physical Review Letters, 86(13):2790, 2001.
  * [28] Gregory Falkovich and Alain Pumir.  Intermittent distribution of heavy particles in a turbulent flow.  Physics of Fluids, 16(7):L47–L50, 2004.
  * [29] Leonid I Zaichik, Olivier Simonin, and Vladimir M Alipchenkov.  Two statistical models for predicting collision rates of inertial particles in homogeneous isotropic turbulence.  Physics of Fluids, 15(10):2995–3005, 2003.
  * [30] Jeremy Bec, Luca Biferale, Guido Boffetta, Antonio Celani, Massimo Cencini, Alessandra Lanotte, S Musacchio, and Federico Toschi.  Acceleration statistics of heavy particles in turbulence.  Journal of Fluid Mechanics, 550:349–358, 2006.
  * [31] Kyle D Squires and John K Eaton.  Preferential concentration of particles by turbulence.  Physics of Fluids A: Fluid Dynamics, 3(5):1169–1178, 1991.
  * [32] Marc Boivin, Olivier Simonin, and Kyle D Squires.  Direct numerical simulation of turbulence modulation by particles in isotropic turbulence.  Journal of Fluid Mechanics, 375:235–263, 1998.
  * [33] Walter C Reade and Lance R Collins.  A numerical study of the particle size distribution of an aerosol undergoing turbulent coagulation.  Journal of Fluid Mechanics, 415:45–64, 2000.
  * [34] Yong Zhou, Anthony S Wexler, and Lian-Ping Wang.  Modelling turbulent collision of bidisperse inertial particles.  Journal of Fluid Mechanics, 433:77, 2001.
  * [35] PK Yeung.  Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations.  Journal of fluid mechanics, 427:241, 2001.
  * [36] PK Yeung.  Lagrangian investigations of turbulence.  Annual review of fluid mechanics, 34(1):115–142, 2002.
  * [37] Takashi Ishihara and Yukio Kaneda.  Relative diffusion of a pair of fluid particles in the inertial subrange of turbulence.  Physics of Fluids, 14(11):L69–L72, 2002.
  * [38] Lance R Collins and Arun Keswani.  Reynolds number scaling of particle clustering in turbulent aerosols.  New Journal of Physics, 6(1):119, 2004.
  * [39] L Biferale, Guido Boffetta, Antonio Celani, BJ Devenish, Alessandra Lanotte, and Federico Toschi.  Multifractal statistics of lagrangian velocity and acceleration in turbulence.  Physical review letters, 93(6):064502, 2004.
  * [40] Jaehun Chun, Donald L Koch, Sarma L Rani, Aruj Ahluwalia, and Lance R Collins.  Clustering of aerosol particles in isotropic turbulence.  Journal of Fluid Mechanics, 536:219, 2005.
  * [41] Luca Biferale, Guido Boffetta, Antonio Celani, Alessandra Lanotte, and Federico Toschi.  Particle trapping in three-dimensional fully developed turbulence.  Physics of Fluids, 17(2):021701, 2005.
  * [42] M Cencini, J Bec, L Biferale, Guido Boffetta, A Celani, AS Lanotte, S Musacchio, and F Toschi.  Dynamics and statistics of heavy particles in turbulent flows.  Journal of Turbulence, (7):N36, 2006.
  * [43] Søren Ott and Jakob Mann.  An experimental investigation of the relative diffusion of particle pairs in three-dimensional turbulent flow.  Journal of Fluid Mechanics, 422:207–223, 2000.
  * [44] Arthur La Porta, Greg A Voth, Alice M Crawford, Jim Alexander, and Eberhard Bodenschatz.  Fluid particle accelerations in fully developed turbulence.  Nature, 409(6823):1017–1019, 2001.
  * [45] Nicolas Mordant, Pascal Metz, Olivier Michel, and J-F Pinton.  Measurement of lagrangian velocity in fully developed turbulence.  Physical Review Letters, 87(21):214501, 2001.
  * [46] Greg A Voth, Arthur La Porta, Alice M Crawford, Jim Alexander, and Eberhard Bodenschatz.  Measurement of particle accelerations in fully developed turbulence.  Journal of Fluid Mechanics, 469:121–160, 2002.
  * [47] Brian Lewis Sawford, PK Yeung, Michael S Borgas, Prakash Vedula, Arthur La Porta, Alice M Crawford, and Eberhard Bodenschatz.  Conditional and unconditional acceleration statistics in turbulence.  Physics of Fluids, 15(11):3478–3489, 2003.
  * [48] Thomas G Mason, Kavita Ganesan, J Harry van Zanten, Denis Wirtz, and Scot C Kuo.  Particle tracking microrheology of complex fluids.  Physical review letters, 79(17):3282, 1997.
  * [49] Thomas A Waigh.  Microrheology of complex fluids.  Reports on progress in physics, 68(3):685, 2005.
  * [50] Gregory Falkovich, K Gawedzki, and Massimo Vergassola.  Particles and fields in fluid turbulence.  Reviews of modern Physics, 73(4):913, 2001.
  * [51] Mierk Schwabe, Sergey Zhdanov, Christoph Räth, David B Graves, Hubertus M Thomas, and Gregor E Morfill.  Collective effects in vortex movements in complex plasmas.  Physical review letters, 112(11):115002, 2014.
  * [52] James J Riley and GS Patterson Jr.  Diffusion experiments with numerically integrated isotropic turbulence.  The Physics of Fluids, 17(2):292–297, 1974.
  * [53] John B McLaughlin.  Aerosol particle deposition in numerically simulated channel flow.  Physics of Fluids A: Fluid Dynamics, 1(7):1211–1224, 1989.
  * [54] Sanat Kumar Tiwari, Amita Das, Predhiman Kaw, and Abhijit Sen.  Longitudinal singular response of dusty plasma medium in weak and strong coupling limits.  Physics of Plasmas, 19(1):013706, 2012.
  * [55] Themistoklis Sapsis and George Haller.  Clustering criterion for inertial particles in two-dimensional time-periodic and three-dimensional steady flows.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 20(1):017515, 2010.
  * [56] Martin R Maxey and James J Riley.  Equation of motion for a small rigid sphere in a nonuniform flow.  The Physics of Fluids, 26(4):883–889, 1983.
  * [57] Abhijit Guha.  Transport and deposition of particles in turbulent and laminar flow.  Annu. Rev. Fluid Mech., 40:311–341, 2008.
  * [58] Jae-Hyung Jeon, Natascha Leijnse, Lene B Oddershede, and Ralf Metzler.  Anomalous diffusion and power-law relaxation of the time averaged mean squared displacement in worm-like micellar solutions.  New Journal of Physics, 15(4):045011, 2013.
  * [59] 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.
  * [60] P Swarztrauber, R Sweet, and John C Adams.  Fishpack: Efficient fortran subprograms for the solution of elliptic partial differential equations.  UCAR Publication, July, 1999.
  * [61] S Ravichandran, P Deepu, and Rama Govindarajan.  Clustering of heavy particles in vortical flows: a selective review.  Sādhanā, 42(4):597–605, 2017.
  * [62] Martin R Maxey.  The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields.  Journal of fluid mechanics, 174:441–465, 1987.
  * [63] G Falkovich, A Fouxon, and MG Stepanov.  Acceleration of rain initiation by cloud turbulence.  Nature, 419(6903):151–154, 2002.
  * [64] Susumu Goto and JC Vassilicos.  Self-similar clustering of inertial particles and zero-acceleration points in fully developed two-dimensional turbulence.  Physics of Fluids, 18(11):115103, 2006.
  * [65] Alec J Petersen, Lucia Baker, and Filippo Coletti.  Experimental study of inertial particles clustering and settling in homogeneous turbulence.  Journal of Fluid Mechanics, 864:925–970, 2019.
  * [66] Sunao Oka and Susumu Goto.  Generalized sweep-stick mechanism of inertial-particle clustering in turbulence.  Physical Review Fluids, 6(4):044605, 2021.
  * [67] Gábor Drótos, Pedro Monroy, Emilio Hernández-García, and Cristóbal López.  Inhomogeneities and caustics in the sedimentation of noninertial particles in incompressible flows.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(1):013115, 2019.


# A State Feedback Controller for Mitigation of Continuous-Time Networked SIS
Epidemics

Yuan Wang  Sebin Gracy  César A. Uribe  Hideaki Ishii  and Karl Henrik
Johansson  Division of Decision and Control Systems, School of Electrical
Engineering and Computer Science, KTH Royal Institute of Technology, and
Digital Futures, Stockholm, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS>Department
of Electrical and Computer Engineering, Rice University, Houston, TX, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Computer Science, Tokyo
Institute of Technology, Yokohama, Japan<EMAIL_ADDRESS>
###### Abstract

The paper considers continuous-time networked susceptible-infected-susceptible
(SIS) diseases spreading over a population. Each agent represents a sub-
population and has its own healing rate and infection rate; the state of the
agent at a time instant denotes what fraction of the said sub-population is
infected with the disease at the said time instant. By taking account of the
changes in behaviors of the agents in response to the infection rates in real-
time, our goal is to devise a feedback strategy such that the infection level
for each agent strictly stays below a pre-specified value. Furthermore, we are
also interested in ensuring that the closed-loop system converges either to
the disease-free equilibrium or, when it exists, to the endemic equilibrium.
The upshot of devising such a strategy is that it allows health administration
officials to ensure that there is sufficient capacity in the healthcare system
to treat the most severe cases. We demonstrate the effectiveness of our
controller via numerical examples.

###### keywords:

Epidemic processes, SIS epidemics, Control of epidemics, Characterization of
endemic equilibrium, Suppressing endemic equilibrium

††thanks: This work was supported in part by the Knut and Alice Wallenberg
Foundation, Wallenberg Scholar under Grant 66469; by a Distinguished Professor
Grant from the Swedish Research Council (Org: JRL, project no: 3058); and by
JSPS under Grant-in-Aid for Scientific Research Grant No. 22H01508; and by the
National Science Foundation under Grants No. 2211815 and No. 2213568

##  1 Introduction

Mathematical epidemiology has a rich history, tracing its roots to the seminal
work by Daniel Bernoulli in (Bernoulli, 1760). Subsequent decades and
centuries have witnessed tremendous advancements in this area, with
contributions spanning across various fields such as physics (Van Mieghem et
al., 2008), computer science (Wang et al., 2003), mathematics (Hethcote,
2000), and so on. Some of the central questions pertain to the ascertainment
of conditions leading to a disease becoming extinct or persisting in the
population. To this end, several models abound in the literature. This paper
is focused on the susceptible-infected-susceptible (SIS) model.

The conventional SIS model in (Kermack and McKendrick, 1927) analyzes the
spread of epidemics from a macroscopic perspective. Recently, to study virus
spread across population subgroups, networked SIS models have been proposed
for both continuous-time systems (Khanafer et al., 2016) and discrete-time
systems (Paré et al., 2020b). In this paper, the focus is on _continuous-time
networked SIS_ models. In this model, an agent is either in the healthy or in
the infected state. A healthy agent could get infected depending on its
infection rate $\beta$, scaled by the interactions it has with its neighbors;
an infected agent recovers based on its healing rate $\gamma$. The terms
$\beta$ and $\gamma$ denote the (reciprocal of) time that it takes for an
agent to get infected, and to get recovered, respectively. We say that the
networked SIS model is in the healthy state if all the agents are healthy, or
equivalently, in the disease-free equilibrium (DFE). If the epidemic remains
persistent, we say that the networked SIS model is in the endemic state.

Stability analysis of continuous-time networked SIS models has been studied in
the literature (Liu et al., 2019; Ye et al., 2022; Fall et al., 2007). In the
context, the results in, among others, (Liu et al., 2019; Ye et al., 2022;
Fall et al., 2007) provide a foundational understanding on what causes an
epidemic to become extinct or to become persistent in the population. In other
words, the results in these works cater to the question: Assuming that there
are no interventions (pharmaceutical or otherwise), under what conditions
would one expect the number of infected individuals in a population to grow
exponentially? However, one of the main challenges that public health
officials confront is ensuring that the available healthcare facilities do not
get overwhelmed. In this context, this means that one of the questions of
interest is the following: What kind of interventions are needed to guarantee
that the proportion of infected in a district or city stays within an
acceptable limit? One approach towards addressing this question would be to
make the individuals aware of the proportion of infected in the population, so
that they can regulate their interactions with others and increase their means
of protection in real-time. From a systems theory viewpoint, the interest,
then, is in designing a (possibly closed-loop) control strategy so that the
state trajectories (which means the infected proportion dynamics of each
population) never exceed a certain pre-specified bound. In this paper, each
agent is equipped with a local pre-specified infection level which can be
different from agent to agent. From the practical viewpoint, such pre-
specified bound is considered as the overwhelmed limit for the healthcare
facilities in each agent.

Control strategies for epidemic models have been presented in (Pasqualetti et
al., 2014; Preciado et al., 2014; Liu et al., 2019; Ye et al., 2022; Morris et
al., 2021; Wang et al., 2021). More specifically, the control schemes proposed
in (Pasqualetti et al., 2014) are for a class of linearized SIS models. The
paper (Preciado et al., 2014) studies the problem of cost-optimal distribution
of resources (such as vaccines and antidotes) in networked SIS models,
whereas, by regarding the healing rate as the local control input, distributed
feedback controllers are proposed in (Liu et al., 2019; Ye et al., 2022). Note
that the goals in (Preciado et al., 2014; Ye et al., 2022; Liu et al., 2019)
do not necessarily ensure the infection levels always staying below a pre-
specified level. In the absence of drugs and vaccines, the control scheme in
(Morris et al., 2021) regards the effective infection rate of each agent as
the control input, and proposed an optimal strategy to reduce the peak in the
macroscopic susceptible-infected-recovered (SIR) model. A similar control
scheme was adopted in (Wang et al., 2022) in the context of opinion dynamics.

For discrete-time networked SIS models, a state feedback controller is
proposed in (Wang et al., 2021), to guarantee that the fraction of infected in
each sub-population stays below half for all time instants. Note that the
discrete-time SIS model is an approximation of the continuous-time SIS model
(Paré et al., 2020b), which in turn is a mean field approximation of a $2^{n}$
state Markov chain model. In order to ensure that the approximation error is
low, several additional assumptions are needed on the sampling period for the
model in (Wang et al., 2021), and furthermore the pre-specified level is the
same for all agents. It is natural to ask how well would the controller behave
in a continuous-time setting. For this reason, we consider a continuous-time
counterpart of the work (Wang et al., 2021) in this paper. Furthermore, we do
not insist on the pre-specified level being the same for all agents, and we
are interested in guaranteeing that even for the closed-loop SIS system, the
local proportion of infected individuals in each node converges either to zero
or to a scalar that is strictly less than said pre-specified level.

_Contributions:_ The main contribution of this paper is to devise a control
scheme for guaranteeing that the fraction of infected individuals in a sub-
population stays within a pre-specified level for all time instants. Our
approach is as follows: First, we modify the continuous-time SIS model in (Liu
et al., 2019) by introducing a parameter that scale the strength of
interconnections between agents, based on their current infection levels. We
then show that for this modified SIS model the following properties hold:

  1. (i)

The proportion of infected individuals in a sub-population does not exceed a
pre-specified level; see Proposition 2.

  2. (ii)

If the spectral abscissa of the state matrix linearized around the healthy
state is not greater than zero, then the DFE is asymptotically stable; see
Proposition 3

  3. (iii)

If the spectral abscissa of the state matrix linearized around the healthy
state is greater than zero, then there exists an endemic equilibrium to which
the dynamics converge for all non-zero initial conditions; see Proposition 4.

This paper is structured as follows. In Section 2, the continuous-time
networked SIS model and feedback control problem of interest are formally
introduced. Our main results and related discussions are given in Section 3. A
numerical example is presented to illustrate our theoretical findings in in
Section 4. Finally, we conclude this paper, and provide future directions in
Section 5.

_Notation:_ Let $\mathbb{R}_{+}$ denote the sets of non-negative real numbers.
We use $[n]$ to denote the set $\\{1,2,...,n\\}$ for any positive integer $n$.
For any two vectors $\mathbf{a},\mathbf{b}\in\mathbb{R}^{n}$, we write
$\mathbf{a}>\mathbf{b}$ if $a_{i}>b_{i}$ for every $i\in[n]$. For a matrix
$A$, let $s(A)$ denote the largest real part among the eigenvalues of $A$. A
diagonal matrix is denoted as $\operatorname{diag}(\cdot)$.

##  2 Problem Formulation

###  2.1 Open-loop networked SIS model

Consider a disease spreading over a network of $n$ agents. The interconnection
among the various agents is represented by a directed graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\\{1,2,\ldots,n\\}$ is the set of agents. Let
$A=[a_{ij}]_{n\times n}$ be the weighted adjacency matrix of $\mathcal{G}$.
The edge set $\mathcal{E}$ is defined as follows: $\mathcal{E}=\\{(j,i)\mid
a_{ij}>0\\}$. More specifically, if a pair of vertices $(j,i)\in\mathcal{E}$,
there is a directed edge from node $j$ to node $i$. The set of in-neighbors
for agent $i$ is denoted as $\mathcal{N}_{i}$. Throughout the rest of this
paper, an agent can be interpreted as a _sub-population_. The weight between
agents can be considered as the frequency of interaction between the sub-
populations. Each sub-population is comprised of individuals and the number of
individuals in a sub-population is fixed. Each individual is assumed to be
either Infected (I) with the disease, or not infected but Susceptible (S) to
the disease. The disease can spread both due to the interaction between
individuals in a sub-population and also across sub-populations.

Let $x_{i}(t)\in[0,1]$ denote the proportion of (sub)population $i\in[n]$ that
is infected with the disease at time $t\in\mathbb{R}_{+}$ (thus $1-x_{i}(t)$
is the proportion of (sub)population $i$ who are susceptible). Let $\beta_{i}$
and $\gamma_{i}$ denote the infection rate and healing rate, respectively, of
agent $i$. The evolution of the infection level of agent $i$ can, then, be
captured by the following differential equation (Fall et al., 2007):

$\dot{x}_{i}(t)=\beta_{i}(1-x_{i}(t))\textstyle\sum_{j=1}^{n}a_{ij}x_{j}(t)-\gamma_{i}x_{i}(t),$
(1)

where $\beta_{i}>0$ and $\gamma_{i}\geq 0$.

Define $\mathbf{x}(t):=[x_{1}(t),x_{2}(t),\ldots,x_{n}(t)]^{T}$,
$B:=\mathrm{\operatorname{diag}}(\beta_{i})$,
$\Gamma:=\mathrm{\operatorname{diag}}(\gamma_{i})$, and
$X(t):=\mathrm{\operatorname{diag}}(\mathbf{x}(t))$ being the diagonal
matrices. Then in vector form the open-loop system (1) can be written as

$\mathbf{\dot{x}}(t)=\big{[}(I-X(t))BA-\Gamma\big{]}\mathbf{x}(t).$ (2)

It is clear that $\mathbf{0}_{n}$ is an equilibrium of the open-loop system
(2). This equilibrium point is referred to as the healthy state, or the
disease-free equilibrium (DFE). Any equilibrium other than $\mathbf{0}_{n}$ is
referred to as the endemic equilibrium, under which the disease persists in at
least one agent.

We need the following assumptions to ensure that system (2) is well-defined.

###### Assumption 1

For each $i\in[n]$, $\beta_{i}>0$ and $\gamma_{i}\geq 0$.

###### Assumption 2

The matrix $A$ is irreducible.

Assumption 1 stipulates that all infection and healing rates must be positive
and non-negtaive, respectively. Assumption 2 states that the adjacency matrix
$A$ is irreducible, which is fulfilled if and only if the underlying graph
$\mathcal{G}$ is strongly connected. A graph $\mathcal{G}$ is strongly
connected if and only if each node $i\in\mathcal{V}$ has a path to every other
node $j\in\mathcal{V}$.

###  2.2 Modified SIS models with local control inputs

Our objective is to ensure that, for all $i\in[n]$ and $t>0$, $x_{i}(t)$
remains upper bounded by a constant ${1}/{c_{i}}$, for some $c_{i}>1$, and
that the dynamics converge to one of the equilibria of the closed-loop system.
That is, we are interested in ensuring that the fraction of infected
population in each node stays below a pre-specified value for all times.
Furthermore, we would like to ensure that
$\lim_{t\rightarrow\infty}x_{i}(t)=\bar{x}_{i}$, where $\bar{x}_{i}$ denotes
some equilibrium of the closed-loop system. To this end, each agent $i\in[n]$
is equipped with a local control input $u_{i}(t)$. As a consequence, System
(1) can be modified as follows:

$\displaystyle\dot{x}_{i}(t)=\beta_{i}(1-x_{i}(t))\textstyle\sum_{j=1}^{n}a_{ij}x_{j}(t)-\gamma_{i}x_{i}(t)+u_{i}(t),$
(3)

where $u_{i}(t)$ is a state feedback control of the form:

$\displaystyle
u_{i}(t)=f_{i}\left(x_{i}(t),\\{x_{j}(t)\\}_{j\in\mathcal{N}_{i}}\right).$ (4)

The state feedback controller $u_{i}(t)$ can be considered as the expected
reduction for growth rate of infection at time $t$ because of pharmaceutical
or non-pharmaceutical interventions. From a practical standpoint, the non-
pharmaceutical interventions may include hand-washing, mask-wearing, social
distancing, etc. The typical pharmaceutical intervention is the distribution
of immunity-boosting or therapeutic drugs. A rough approach to estimating such
$u_{i}(t)$ is to compare the growth rate of infection before and after the
interventions in a short time period. However, we note that these consequences
of the mentioned pharmaceutical or non-pharmaceutical interventions are very
hard to measure or estimate accurately. There are some practical limitations
on how the mentioned interventions can affect the spread. The effective
disease prevention measures, according to the type of disease, are quite
different. As a simple example, mask-wearing may slow down the spreading of
airborne diseases. However, it may not mitigate the infection of waterborne
diseases. The modeling for effects of detailed interventions are out of scope
in this paper. See Remarks 1 and 2 for a slightly more detailed explanation
for the proposed feedback controller.

We formally state the problem being investigated in this paper.

###### Problem 1

Consider the closed-loop System (3). Based on knowledge of $\beta_{i}$,
$x_{i}(t)$, $c_{i}$, $a_{ij}$ and $x_{j}(t)$ where $j\in\mathcal{N}_{i}$,
design a local state feedback controller in accordance with (4) such that

  1. i)

for all $i\in[n]$ and $t>0$, $x_{i}(t)\leq{1}/{c_{i}}$; and

  2. ii)

$\lim_{t\rightarrow\infty}x_{i}(t)=x_{i}^{*}$, where either $x_{i}^{*}=0$ or
$x_{i}^{*}>0$.

A discrete-time version of Problem 1 has been partially solved in (Wang et
al., 2021) with the local state feedback controller having a specific
structure, namely, i) $c_{i}=2$ for all $i\in[n]$, and ii) the strength of
interconnections modified as a function of the current infection levels of
each node and all of its neighbors.

Our approach towards addressing Problem 1 involves two steps: First, given
some $c_{i}>1$, we identify $u_{i}(t)$ that obeys (4), such that
$x_{i}(t)\leq{1}/{c_{i}}$ for all $t$. Second, we show that, for the proposed
feedback controller $u_{i}(t)$, the closed-loop system (1) always converges to
some equilibria. To elucidate more on the second part of our solution, we now
recall a result which fully characterizes the class of equilibria for the
open-loop system (1), and admits $\gamma_{i}=0$.

###### Proposition 1

(Liu et al., 2019, Propositions 2 and 3) Consider the open-loop System (2)
under Assumptions 1 and 2. The following statements hold:

  1. 1.

If $s(BA-\Gamma)\leq 0$, then, the disease-free equilibrium is the unique
equilibrium, and it is globally asymptotically stable.

  2. 2.

If $s(BA-\Gamma)>0$, then, other than the disease-free equilibrium, there
exists a unique endemic equilibrium, which is globally asymptotically stable
on $[0,1]^{n}\setminus\\{\mathbf{0}_{n}\\}$.

In the rest of this paper, we will show that the result in Proposition 1 can
be, with a suitable choice of $u_{i}(t)$ as in (4), extended for the closed-
loop system in (3).

##  3 Main results

In this section, we devise a local state feedback controller to solve Problem
1. Observe that for the closed-loop System (3), since the term
$\beta_{i}(1-x_{i}(t))\textstyle\sum_{j=1}^{n}a_{ij}x_{j}(t)$ is positive, the
terms that contribute to the increase in the state value are exactly the
aforementioned ones. Hence, in order to ensure that the state trajectories
remain below $1/c_{i}$ for all times, it is natural to consider the controller
that partially offsets the increase caused by the aforementioned terms. For
this reason, we consider the following control input

$\displaystyle
u_{i}(t)=-\beta_{i}c_{i}x_{i}(t)(1-x_{i}(t))\textstyle\sum_{j=1}^{n}a_{ij}x_{j}(t).$
(5)

Consequently, by plugging (5) into (3), we obtain:

$\displaystyle\dot{x}_{i}(t)=\beta_{i}(1-c_{i}x_{i}(t))(1-x_{i}(t))\textstyle\sum_{j=1}^{n}a_{ij}x_{j}(t)-\gamma_{i}x_{i}(k).$
(6)

Let $D:=\mathrm{\operatorname{diag}}(c_{i})$. Then, the closed-loop System (6)
can be written as

$\displaystyle\mathbf{\dot{x}}(t)=\big{[}(I-DX(t))(I-X(t))BA-\Gamma\big{]}\mathbf{x}(t).$
(7)

We refer to System (7) as the _controlled system_.

###### Remark 1

Define $b_{i}(t):=1-c_{i}x_{i}(t)$, and observe that $b_{i}(t)\in[0,1]$. That
is, the term $b_{i}(t)$ scales the infection rate for agent $i$ at time
instant $t$. More concretely, and especially from a practical standpoint, the
parameter $b_{i}(t)$ may be thought of as a term/quantity, possibly
communicated to all the agents by a central entity such as health
administration officials, that informs each agent as to by how much should it
modify the strength of its interactions with its neighbors. As a consequence,
the infection rate of agent $i$ is reduced from $\beta_{i}$ to
$b_{i}(t)\beta_{i}$.

###### Remark 2

The proposed feedback controller is based on the ideal environment, and there
are certain assumptions made for practical implementation as follows: Agents
must be able to get the exact states, including their own and their neighbors,
and the proposed control strategy must be faithfully implemented by agents. We
must implement that, in practice, there are asymptomatic members (who are in
the latent period), and it is hard to handle the exact proportion of infected
people. The design of a feedback controller which is based on noisy state
estimation, as well as the epidemic control policy within the antagonistic
environment (where some agents may not follow the control strategy) are left
for our future studies.

We have the following assumption.

###### Assumption 3

For each $i\in[n]$, $x_{i}(0)\in[0,1/c_{i}]$.

Assumption 3 ensures that all initial states are within the desired bound,
since otherwise the control objectives cannot be achieved.

With $u_{i}(t)$ as in (4), the infection levels in the closed-loop System (3)
do not exceed the pre-specified level, ${1}/{c_{i}}$. We have the following
propositions. We skip the proof of the propositions due to space reasons.

###### Proposition 2

Consider the controlled System (7) under Assumptions 1 and 3. Then
$x_{i}(t)\in[0,1/c_{i}]$ for all $i\in[n]$ and $t\geq 0$.

Proposition 2 establishes that the set $[0,1/c_{i}]$ is positively invariant
for each agent $i$ with respect to the closed-loop System (6). Positive
invariance of the set $[0,1]^{n}$ with respect to the uncontrolled System (1)
has been established in (Liu et al., 2019, Lemma 8). Under Assumptions 1 and
3, the proposed controller shrinks the positive invariant set for each agent
$i$ from $[0,1]$ to $[0,1/c_{i}]$. In a practical sense, the local infection
upper bound $1/c_{i}$ can be regarded as the capacity of healthcare facilities
in a local population $i$. As one of the main advantages of the controlled
system, if the initial infection level do not get overwhelmed, then the
epidemic remains manageable.

While Proposition 2 guarantees that the infection level does not exceed a
certain pre-specified level, it does not give any guarantees regarding the
long-term behavior of the system. That is, with the aforementioned controller
in place, does the disease die out, or does it become endemic? We address the
same in the rest of this subsection.

It turns out that convergence guarantees, similar to that for the open-loop
case discussed in Proposition 1, can be given even for the closed-loop System
(7). The next two propositions handle the case when $s(BA-\Gamma)\leq 0$ and
$s(BA-\Gamma)>0$, respectively.

###### Proposition 3

Consider the controlled System (7) under Assumptions 1–3. If $s(BA-\Gamma)\leq
0$, then the disease-free equilibrium is asymptotically stable with the domain
of attraction $[0,1/c_{i}]$ for each $i\in[n]$.

Proposition 3 states that, irrespective of whether an agent is initially
infected or healthy, as long as $s(BA-\Gamma)\leq 0$ the infection levels of
the closed-loop system converge to the disease-free equilibrium. Note that
Proposition 2 depends only on the healing rate, infection rate and the network
structure. If the healing rate of each agent dominates the infection rate
(which is scaled by the interconnection weights between agents), the epidemic
dies out naturally.

###### Proposition 4

Consider the controlled System (7) under Assumptions 1–3. If $s(BA-\Gamma)>0$,
then there exists a unique endemic equilibrium $\overline{\mathbf{x}}$ such
that $0<[\overline{\mathbf{x}}]_{i}\leq{1}/{c_{i}}$ for each $i\in[n]$, and
$c_{i}>1$. Furthermore, if for each $i\in[n]$, $c_{i}\geq 2$, then the endemic
equilibrium is asymptotically stable with the domain of attraction
$(0,1/c_{i}]$.

Proposition 4 states that when $s(BA-\Gamma)>0$ System (7) has two equilibria,
namely, the disease-free equilibrium and the endemic equilibrium
$\mathbf{\overline{x}}$. Furthermore, as long as the initial state is non-
trivial and $c_{i}\geq 2$, the state of the closed-loop system asymptotically
converges to $\mathbf{\overline{x}}$, which means that the disease-free
equilibrium is unstable. In other words, our proposed controller cannot
eradicate the disease; rather, it ensures that each population $i$ has a
stable non-trivial infection level strictly smaller than ${1}/{c_{i}}$.

We are now ready to state the main result of the paper. It summarizes the
results from Propositions 2–4. Recall that the overall closed-loop system in
(7) consists of the agent dynamics in (3) and the local control inputs in (5).

###### Theorem 1

Consider the closed-loop System (3) under Assumptions 1–3. The local nonlinear
state feedback control law (5), guarantees $x_{i}(t)\leq{1}/{c_{i}}$ for all
$i\in[n]$ and $t\geq 0$. Furthermore, if $s(BA-\Gamma)\leq 0$, then
$\lim_{t\rightarrow\infty}x_{i}(t)=0$; otherwise if $s(BA-\Gamma)>0$, and
$c_{i}\geq 2$ for all $i\in[n]$, then
$\lim_{t\rightarrow\infty}x_{i}(t)=x_{i}^{*}$, where $x_{i}^{*}>0$.

As a practical consequence, Theorem 1 indicates how the faithfully implemented
epidemic prevention interventions influence the infection spreading. Note that
the convergence result for the case when $s(BA-\Gamma)>0$, and $1<c_{i}<2$ is
not provided. From a technical viewpoint, the main hindrance is that the
current Lyapunov candidate is not strictly decreasing for the aforementioned
case. Simulations, however, indicate that even for the case when
$s(BA-\Gamma)>0$, and $1<c_{i}<2$ convergence to the endemic equilibrium can
be achieved. A rigorous analysis of the aforementioned case is left for future
work.

##  4 Numerical Example

In this section, we first illustrate the results from Section 3, and then show
some interesting behavior via simulations.

The simulations are presented in a multi-agent network with $n=100$ sub-
populations. The sub-populations are uniformly and randomly located in an area
of $100\times 100$. All agents within a distance of $r=25$ are allowed to
communicate with each other. Following a geometric graph model, for any two
sub-population $i$ and $j$ with the distance between them being less than $r$,
there exists an undirected edge between the said nodes, with the weight
$0<a_{ij}<1$. The nodes are indexed from $1$ to $100$. Each simulation follows
the same initial conditions, such that the nodes $1$ to $10$ have $10$ percent
of their population infected, and the other nodes, i.e., numbered $11$ to
$100$ are susceptible. The first topology of the multi-agent network is as
depicted in Fig. 1, where the nodes that are healthy (or susceptible) are
marked blue, while the nodes that are infected are marked red. It can be
easily verified that this network is strongly connected. The edge weights
$a_{ii}$ and $a_{ij}$ are set to $a_{ii}=0.3,a_{ij}=0.003$, for all connected
nodes $i$ and $j$ such that $a_{ij}\neq 0$.

Figure 1: Network topology 1 with communication radius $r=25$.

The pre-specified infection level $1/c_{i}$ is chosen as follows: For $1\leq
i\leq 20$, choose $1/c_{i}=0.5$; for $21\leq i\leq 40$, choose $1/c_{i}=0.45$;
for $41\leq i\leq 60$, choose $1/c_{i}=0.3$; for $61\leq i\leq 80$, choose
$1/c_{i}=0.25$; and for $81\leq i\leq 100$, choose $1/c_{i}=0.2$. Hence, for
all $i\in[n]$, the initial states, which are at most $0.1$, are strictly less
than $1/c_{i}$.

We illustrate the case $s(BA-\Gamma)\leq 0$ in the first simulation. Set
$\beta_{i}=0.3$, $\gamma_{i}=0.5$ for all $i$. Then, we have
$s(BA-\Gamma)=-0.313$ for the network in Fig. 1. The trajectories of the
states are shown in Fig. 2. Both the states in uncontrolled and controlled
systems stays strictly below the minimum pre-specific infection level
$1/c_{i}=0.2$. Furthermore, the dynamics of the closed-loop system converges
to the disease-free equilibrium, which is in line with the result in Theorem
1.

Figure 2: Time responses for both uncontrolled and controlled systems when
$s(BA-\Gamma)=-0.313$.

Then, we study the case where $s(BA-\Gamma)>0$. Set $\beta_{i}=0.8$,
$\gamma_{i}=0.3$ for all $i$, and consequently, $s(BA-\Gamma)=0.198$. The
disease persists in both the uncontrolled and controlled systems; see Fig. 3.
The uncontrolled system converges to the endemic equilibrium around $0.48$.
For the controlled system, the states asymptotically converge to the endemic
equilibrium, and they are all strictly below the minimum pre-specific
infection level $1/c_{i}=0.2$, which is consistent with the result in Theorem
1. In the controlled system, even though the disease persists, the infection
levels are significantly suppressed in all sub-populations. Comparing the
states with different pre-specific infection levels, the sub-populations that
choose a smaller $1/c_{i}$ generally have a lower infection level.

Figure 3: Time responses for both uncontrolled and controlled systems when
$s(BA-\Gamma)=0.198$, and $c_{i}\geq 2$.

We check the convergence result for the case when $s(BA-\Gamma)>0$, and
$1<c_{i}<2$. Set $\beta_{i}=0.8$, $\gamma_{i}=0.3$ for all $i$, and
consequently, $s(BA-\Gamma)=0.198$. For $1\leq i\leq 20$, we reset the pre-
specified infection level by $1/c_{i}=0.9$, and thus $1<c_{i}<2$ for these
agents. The agent dynamics are shown in Fig. 4. It can be seen that the
infection level for $1\leq i\leq 20$ (blue lines) reach the endemic
equilibrium.

Figure 4: Time responses for both uncontrolled and controlled systems when
$s(BA-\Gamma)=0.198$ and $c_{i}<2$.

We next study how the higher proximity between agents affects the spreading
dynamics. Keep the location of all agent as that in Fig. 1, and improve the
communication radius by $r=50$. Set $\beta_{i}=0.8$, $\gamma_{i}=0.3$ for all
$i$, and consequently, $s(BA-\Gamma)=0.276$. The time responses are shown in
Fig. 5. The disease spreading is strictly suppressed below the local pre-
specific infection levels in all sub-populations. However, comparing with the
time responses in Fig. 3 (where communication radius is $r=25$), the disease
persists with higher infection proportions in both the uncontrolled and
controlled systems.

Figure 5: Time responses for both uncontrolled and controlled systems when
$s(BA-\Gamma)=0.276$ and $r=50$.

Finally, we check the influence of higher concentration of agents. We test the
disease spreading within a network topology as shown in Fig 6, where the
communication radius is set by $r=25$. Comparing with that in Fig 1, there is
a cluster of sub-populations on the left bottom side. Set $a_{ii}=0.3$,
$a_{ij}=0.003$, and $\beta_{i}=0.8$, $\gamma_{i}=0.3$ for all $i$.
Consequently, $s(BA-\Gamma)=0.214$. The infection dynamics are shown in Fig.
7. Comparing with that in Fig. 3, it can be seen that with higher
concentration of agents, the disease persists with higher infection
proportions. As a practical consequence of the simulation results, avoiding
crowd aggregation is beneficial to the suppression of infectious diseases.

Figure 6: Network topology 2 with communication radius $r=25$. Figure 7: Time
responses within network topology 2 for both uncontrolled and controlled
systems when $s(BA-\Gamma)=0.214$ and $c_{i}\geq 2$.

##  5 Conclusion

We considered an SIS epidemic spreading over a networked population. We
designed a state feedback controller where the objective was to ensure that
the infection level of each agent never exceeded a pre-specified value. With
the controller in place, we identified a condition under which the disease
gets eradicated (resp. remains persistent in the population). One of the
future directions would be to extend both the feedback controller and analysis
approaches for the other disease-spreading models, such as SIR, SEIR, etc.
Another interesting direction would be to design state feedback controllers
that maintain the infection below a certain level for time-varying graphs, as
well as the epidemic control strategy with noisy data in the antagonistic
environment.

## References

  * Bernoulli (1760) Bernoulli, D. (1760).  Essai d’une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l’inoculation pour la prévenir.  _Histoire de l’Acad., Roy. Sci.(Paris) avec Mem_ , 1–45.
  * Fall et al. (2007) Fall, A., Iggidr, A., Sallet, G., and Tewa, J.J. (2007).  Epidemiological models and Lyapunov functions.  _Mathematical Modelling of Natural Phenomena_ , 2(1), 62–83.
  * Hethcote (2000) Hethcote, H.W. (2000).  The mathematics of infectious diseases.  _SIAM Review_ , 42(4), 599–653.
  * Kermack and McKendrick (1927) Kermack, W.O. and McKendrick, A.G. (1927).  A contribution to the mathematical theory of epidemics.  _Proceedings of the Royal Society of London. Series A_ , 115(772), 700–721.
  * Khanafer et al. (2016) Khanafer, A., Başar, T., and Gharesifard, B. (2016).  Stability of epidemic models over directed graphs: A positive systems approach.  _Automatica_ , 74, 126–134.
  * Liu et al. (2019) Liu, J., Paré, P.E., Nedić, A., Tang, C.Y., Beck, C.L., and Başar, T. (2019).  Analysis and control of a continuous-time bi-virus model.  _IEEE Transactions on Automatic Control_ , 64(12), 4891–4906.
  * Morris et al. (2021) Morris, D.H., Rossine, F.W., Plotkin, J.B., and Levin, S.A. (2021).  Optimal, near-optimal, and robust epidemic control.  _Communications Physics_ , 4(1), 1–8.
  * Paré et al. (2020b) Paré, P.E., Liu, J., Beck, C.L., Kirwan, B.E., and Başar, T. (2020b).  Analysis, estimation, and validation of discrete-time epidemic processes.  _IEEE Transactions on Control Systems Technology_ , 28(1), 79–93.
  * Pasqualetti et al. (2014) Pasqualetti, F., Zampieri, S., and Bullo, F. (2014).  Controllability metrics, limitations and algorithms for complex networks.  _IEEE Transactions on Control of Network Systems_ , 1(1), 40–52.
  * Preciado et al. (2014) Preciado, V.M., Zargham, M., Enyioha, C., Jadbabaie, A., and Pappas, G.J. (2014).  Optimal resource allocation for network protection against spreading processes.  _IEEE Transactions on Control of Network Systems_ , 1(1), 99–108.
  * Van Mieghem et al. (2008) Van Mieghem, P., Omic, J., and Kooij, R. (2008).  Virus spread in networks.  _IEEE/ACM Transactions On Networking_ , 17(1), 1–14.
  * Wang et al. (2003) Wang, Y., Chakrabarti, D., Wang, C., and Faloutsos, C. (2003).  Epidemic spreading in real networks: An eigenvalue viewpoint.  _in Proc. 22nd International Symposium on Reliable Distributed Systems_ , 25–34.
  * Wang et al. (2021) Wang, Y., Gracy, S., Ishii, H., and Johansson, K.H. (2021).  Suppressing the endemic equilibrium in SIS epidemics: A state dependent approach.  _in Proc. 11th IFAC Symposium on Biological and Medical Systems_ , 163–168.
  * Wang et al. (2022) Wang, Y., Ishii, H., Bonnet, F., and Défago, X. (2022).  Resilient consensus for multi-agent systems under adversarial spreading processes.  _IEEE Transactions on Network Science and Engineering_ , 9(5), 3316–3331.
  * Ye et al. (2022) Ye, M., Liu, J., Anderson, B.D.O., and Cao, M. (2022).  Applications of the Poincare-Hopf Theorem: Epidemic models and Lotka-Volterra systems.  _IEEE Transactions on Automatic Control_ , 67(4), 1609–1624.