# Instability of a dusty vortex

Shuai Shuai 1 and Darish Jeswin Dhas 2 and Anubhab Roy2 and M. Houssem
Kasbaoui 1<EMAIL_ADDRESS>1School for Engineering of Matter, Transport and
Energy, Arizona State University, Tempe, AZ 85281, USA. 2Department of Applied
Mechanics, Indian Institute of Technology Madras, Chennai, 600036, India.

###### Abstract

We investigate the effect of inertial particles dispersed in a circular patch
of finite radius on the stability of a two-dimensional Rankine vortex in semi-
dilute dusty flows. Unlike the particle-free case where no unstable modes
exist, we show that the feedback force from the particles triggers a novel
instability. The mechanisms driving the instability are characterized using
linear stability analysis for weakly inertial particles and further validated
against Eulerian-Lagrangian simulations. We show that the particle-laden
vortex is always unstable if the mass loading $M>0$. Surprisingly, even non-
inertial particles destabilize the vortex by a mechanism analogous to the
centrifugal Rayleigh-Taylor instability in radially stratified vortex with
density jump. We identify a critical mass loading above which an eigenmode m
becomes unstable. This critical mass loading drops to zero as m increases.
When particles are inertial, modes that fall below the critical mass loading
become unstable, whereas, modes above it remain unstable but with lower growth
rates compared to the non-inertial case. Comparison with Eulerian-Lagrangian
simulations shows that growth rates computed from simulations match well the
theoretical predictions. Past the linear stage, we observe the emergence of
high-wavenumber modes that turn into spiraling arms of concentrated particles
emanating out of the core, while regions of particle-free flow are sucked
inward. The vorticity field displays similar pattern which lead to the
breakdown of the initial Rankine structure. This novel instability for a dusty
vortex highlights how the feedback force from the disperse phase can induce
the breakdown of an otherwise resilient vortical structure.

###### keywords:

keyword 1, keyword 2, keyword 3

##  1 Introduction

Vortical flows with heavy dispersed particles occur in diverse scenarios, both
of natural and engineering origins (Balachandar & Eaton, 2010; Guazzelli &
Morris, 2011). To cite a few examples, they include atmospheric funnel-type
vortices (such as ‘dust devil’)(Bluestein et al., 2004), centrifugal
separation devices, swirling biomass combustors, and aircraft trailing
vortices with condensation droplets (Paoli & Shariff, 2016; Paoli et al.,
2008). The rudimentary understanding of the process is that particles that are
denser than the surrounding fluid drift away from intense vortical regions at
a rate controlled by particle inertia. However, even in elementary vortical
flows such as columnar vortices, the dispersion of inertial particles (micron-
sized solid spheres, or liquid droplets, sufficiently small to remain
spherical under the action of surface tension) is difficult to predict
accurately. This uncertainty is due to models and simulations often omitting
the particle feedback force. Yet, the latter may be significant even in dilute
conditions (particle volume fractions is
$10^{-5}<\langle\phi_{p}\rangle<10^{-3}$) provided that the density ratio is
high ($\rho_{p}/\rho_{f}=O(1000)$), e.g. when the carrier fluid is gas. In
this study, we shall show that particle feedback arising from two-way momentum
coupling activates a novel vortex instability, and this instability is
controlled by particle inertia and vortex structure.

To discuss the stability properties in a vortex flow, we adopt the Rankine
vortex as a typical vortex tube, i.e., one with uniform vorticity within the
core and zero outside (Saffman, 1995). Many measurements of atmospheric
phenomena show that the Rankine vortex is suitable for describing velocity
structure of naturally occuring vortical flows. Doppler radar measurements
prove that tornados are well approximated by the Rankine vortex profile
(Bluestein et al., 2003). A mesocyclone is another example of atmospheric flow
that is well approximated by a Rankine vortex with an inner core of solid
rotation as large as 5km (Brown et al., 2005).

The simple velocity profile of the Rankine vortex makes stability
investigations amenable to analytical treatments, thus, often help in better
identifying the underlying physics. Lord Kelvin (Thomson, 1880; Lamb, 1993)
analyzed the linear stability of an incompressible Rankine vortex. He showed
that for 2D perturbations characterized by mode number $m$, only one real
eigenvalue exists, indicating that the perturbation will propagate with no
growth or decay, corresponding to neutral stability. Besides these neutrally
stable discrete modes, a Rankine vortex also supports continuous spectrum
modes whose combination may lead to algebraic growth for short times (Roy &
Subramanian, 2014). The inclusion of additional physics can destabilize a
vortex column. For instance, acoustic waves radiating to infinity can
destabilize a Rankine vortex by extracting energy from the vortex core
(Broadbent & Moore, 1979). Analogous destabilization mechanisms have been
shown to exist for shallow water flows (Ford, 1994) and stratified flows
(Schecter & Montgomery, 2004; Le Dizès & Billant, 2009), where outgoing waves
draw energy from the mean flow. Recently a new instability mechanism has been
found for a Rankine vortex in dilute polymer solutions (Roy et al., 2021). The
instability occurs due to the resonant interaction of a pair of elastic shear
waves aided by the differential convection of the irrotational shearing flow
outside the vortex core. Compressibility effects may also destabilize a
Rankine vortex as shown by Sozou & Lighthill (1987).

The configuration of a radially stratified vortex has special relevance,
since, as we shall discuss in detail later, it bears analogy with a particle-
laden vortex in the limit of zero particle inertia. If density increases
monotonically with radius, Fung & Kurzweg (1975) showed that the flow is
stable to both axisymmetric and non-axisymmetric modes. In contrast, a heavy-
cored vortex may become unstable, as shown by Sipp et al. (2005a), through two
mechanisms, one due to a centrifugal instability involving short-axial
wavelength modes, and the second due to a Rayleigh-Taylor instability
involving two-dimensional modes. Similarly, Joly et al. (2005) found that
heavy-cored two-dimensional vortices are subject to a Rayleigh-Taylor
instability. They attributed the basic mechanism to baroclinic vorticity
generation caused by the misalignment between density gradient and centripetal
acceleration. Using numerical simulations, Joly et al. (2005) showed that the
unstable modes lead to the formation of spiraling arms in the density and
vorticity fields that roll-up eventually as nonlinear effects emerge. Dixit &
Govindarajan (2011) analyzed the case of a radially stratified vortex with a
density jump. They showed that both heavy-cored and light-cored vortices could
be unstable if the density jump is sufficiently steep. They proposed that the
instability is triggered by the wave-interaction mechanism between the Kelvin
wave located at the vortex core edge and counter-propagating internal waves
caused by the centrifugal force (see Carpenter et al. (2011); Balmforth et al.
(2012) for discussions on instabilities originating from interaction between
waves riding on vorticity/density jumps). When the vorticity and density jumps
are sufficiently separated, they act as “independent” neutral modes; the two
interfacial displacements are in phase with each other and out of phase with
the radial velocity. Thus, the radial velocity neither aid nor retard the
interfacial displacement. When the jumps are located close to each other, the
interfacial displacements are out of phase, and the radial velocity reinforces
the interfacial displacement at the density jump resulting in a growing mode.
Dixit & Govindarajan (2011) also found that a heavy-cored vortex is stabilized
when the density jump is placed in a region immediately outside the vortex
core. Our analysis in section §3.3 will reveal that, when we neglect particle
inertia but account for non-zero mass loading, the system analyzed in the
present study mimics a non-Bousinessq fluid. Thus, the observed instability
has an interesting analogy with the physical mechanisms elucidated in the
investigations mentioned above.

When inertial particles heavier than the carrier fluid disperse in the flow,
they tend to migrate from regions of high rotation rate to regions of high
strain rate. This mechanism known, as preferential concentration (Squires &
Eaton, 1990, 2014; Wang & Maxey, 1993), results from a slip velocity between
particles and flow, which increases with increasing particle inertia.
Preferential concentrations leads to the clustering of inertial particles
around vortical cores.

Although particle dispersion in vortical flows exhibits rich dynamics, the
disperse phase does not alter the fundamentally stability of the flow unless
two-way coupling is considered. Recently there has been a strong interest in
the modulation of flow instabilities due to the feedback force from small,
heavy particles. Magnani et al. (2021) considered the effect of particle
inertia in two-way coupled dusty Rayleigh-Taylor turbulence. They observed
that the system behaves similarly to an equivalent denser fluid for low
inertia particles. When particle inertia increases, turbulent mixing gets
delayed. The non-monotonic role of the disperse phase was also observed by
Sozza et al. (2022) in their stability study of two-way coupled dusty
Kolmogorov flow. When particle inertia is weak, the instability is enhanced.
However, for larger Stokes number, particles can both stabilize or destabilize
the Kolmogorov flow, with a non-monotonic dependence on the mass fraction.
Inertial particles can also trigger instabilities that would not otherwise
exist in particle-free flows. Kasbaoui et al. (2015) showed that the interplay
between preferential concentration and gravitation settling triggers a non-
modal instability in simple shear flow. This instability may bootstrap
additional Rayleigh-Taylor and particle-trajectory crossing instabilities
(Kasbaoui et al., 2019a) and plays a role in the attenuation of turbulent
fluctuations in particle-laden homogeneously sheared turbulence (Kasbaoui et
al., 2019b; Kasbaoui, 2019).

In the context of vortical flows, few studies considered the impact of two-way
coupling on the stability of such flows. Marshall (2005) studied the effect of
turbulence on dispersed particles in a vortex column under one-way coupling,
i.e., neglecting the particle feedback force. He showed that inertial
particles initially inside the vortex core are driven out, forming
concentrated ring-like structures. Turbulence breaks these structures into
smaller sections. Ravichandran & Govindarajan (2015) studied the formation of
particle clusters around isolated and ensembles of two-dimensional vortices in
one-way coupling. They showed that particles initially within a critical
radius, that depends on the vortex circulation and particle response time,
form caustics, thus, potentially leading to higher collision rates. When two-
way coupling is considered, the dispersion of inertial particles may lead to a
modulation of the carrier flow. Druzhinin (1994) showed that the outward
ejection of inertial particles from a particle-rich vortex core causes the
attenuation of vorticity in the core. Considering two-dimensional axisymmetric
perturbations, he showed that the particles form a ring-shaped cluster that
expands outwardly, akin to a concentration wave, at a rate controlled by
Stokes number $\mathrm{St}_{\Gamma}=\tau_{p}/\tau_{f}$, where $\tau_{p}$ is
particle response time, and $\tau_{f}$ is fluid characteristic timescale.
These prior observations were corroborated recently by Shuai & Kasbaoui (2022)
in two-way coupled Eulerian-Lagrangian simulations of particle-laden Lamb-
Oseen vortex. The authors found that particle rings grow on a timescale
$\tau_{c}=\tau_{f}/\mathrm{St}_{\Gamma}$, and proposed an expression for their
expansion rate. Similar to (Druzhinin, 1994), Shuai & Kasbaoui (2022) observed
that particle-feedback on the fluid drives a faster decay of the vortex tube.
Further, the observation of naturally emerging azimuthal perturbations in the
vorticity field suggests that an instability activated by two-way coupling may
be at play. These observations motivate the present study on the stability of
a particle-laden vortex.

This paper is organized as follows. Section §2 presents evidence of
instability in a two-way coupled particle-laden vortex. The governing equation
and linear stability analysis for small inertia particles are presented in §3.
In §4, we compare the analytical solutions with the results from Eulerian-
Lagrangian simulations of the two-way coupled system. Concluding remarks are
given in §5.

##  2 Evidence of instability in a two-way coupled particle-laden vortex

In this section, we show that the modulation of a prototypical vortex, the
Rankine vortex here, is driven by an instability activated by two-way
coupling. We illustrate this behavior in a sample flow using Eulerian-
Lagrangian simulations.

###  2.1 Eulerian-Lagrangian method

The Eulerian-Lagrangian simulations presented in this work are based on the
volume-filtered formulation (Anderson & Jackson, 1967; Jackson, 2000;
Capecelatro & Desjardins, 2013). In this formulation, the fluid phase is
treated in the Eulerian frame, whereas the particle phase is treated in the
Lagrangian frame, i.e., each particle is individually tracked.

The mass and momentum conservation equations for the carrier phase in the
semi-dilute regime, are given by the incompressible Navier-Stokes equations,

$\displaystyle\nabla\cdot\mathbf{u}_{f}$ $\displaystyle=$ $\displaystyle 0,$
(1) $\displaystyle\rho_{f}\left(\frac{\partial\mathbf{u}_{f}}{\partial
t}+\nabla\cdot\left(\mathbf{u}_{f}\mathbf{u}_{f}\right)\right)$
$\displaystyle=$ $\displaystyle-\nabla
p+\mu_{f}\nabla^{2}\mathbf{u}_{f}+\mathbf{F}_{p},$ (2)

where $\mathbf{u}_{f}$ is the fluid velocity, $p$ is the pressure, $\rho_{f}$
is the fluid density, and $\mu_{f}$ is the fluid viscosity. The term
$\mathbf{F}_{p}$ represents momentum exchange between particles and fluid,
which is expressed as

$\mathbf{F}_{p}=-\phi_{p}\nabla\cdot\tau|_{p}+\rho_{p}\phi_{p}\frac{(\mathbf{u}_{p}-\mathbf{u}_{f|p})}{\tau_{p}},$
(3)

where $\tau=-p\mathcal{I}+\mu(\nabla u_{f}+\nabla u_{f}^{T})/2$ is the total
fluid stresses, $\phi_{p}$ is the particle volume fraction, $\mathbf{u}_{p}$
is the Eulerian particle velocity, and $(\cdot)_{|p}$ indicates fluid
quantities at the particle locations. The first term in (3) is the stresses
exerted by the undisturbed flow at the particle location. The second term is
the stresses caused by the presence of particles which are represented using
Stokes drag. When the density ratio is large ($\rho_{p}/\rho_{f}\gg 1$), as is
the case presently, Stokes drag dominates the momentum exchange.

From a scaling analysis of the drag force in (3), it is clear that the mass
loading $M=\rho_{p}\langle\phi_{p}\rangle/\rho_{f}$ determines the strength of
the coupling. Thus, if the mass loading is vanishingly small, the particle
phase has little effect on the stability of the Rankine vortex. Conversely, if
the mass loading is large, the interaction between two phases triggers a
significant source or sink of energy that may enhance or attenuate
perturbations in the flow (Kasbaoui et al., 2015).

The particle phase is described in the Lagrangian frame. The motion of the
$i$-th particle is given by (Maxey & Riley, 1983)

$\displaystyle\frac{d\mathbf{x}_{p}^{i}}{dt}(t)$ $\displaystyle=$
$\displaystyle{\mathbf{u}_{p}^{i}}(t),$ (4)
$\displaystyle\frac{d\mathbf{u}_{p}^{i}}{dt}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{\rho_{p}}\nabla\cdot\tau(\mathbf{x}_{p}^{i},t)+\frac{\mathbf{u}_{f}(\mathbf{x}_{p}^{i},t)-\mathbf{u}^{i}_{p}(t)}{\tau_{p}},$
(5)

where $\mathbf{x}_{p}^{i}$, and $\mathbf{u}_{p}^{i}$, are the position and
velocity of the “$i$”-th particle, respectively,
$\tau_{p}=\rho_{p}d_{p}^{2}/(18\mu)$ is the particle response time, and
$d_{p}$ the particle diameter. In this study, gravity is ignored in order to
focus on inertial effects. Other interactions, including collisions, are
negligible due to the large density ratio and low volume fraction. Note that
in the equations above, the instantaneous particle volume fraction and
Eulerian particle velocity field are obtained from Lagrangian quantities using

$\displaystyle\phi_{p}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}V_{p}g\left(\left|\left|\mathbf{x}-\mathbf{x}_{p}^{i}\right|\right|\right),$
(6) $\displaystyle\phi_{p}\mathbf{u}_{p}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\mathbf{u}_{p}^{i}(t)V_{p}g\left(\left|\left|\mathbf{x}-\mathbf{x}_{p}^{i}\right|\right|\right),$
(7)

where $V_{p}=\pi d_{p}^{3}/6$ is the particle volume, $g$ represents a
Gaussian filter kernel of size $\delta_{f}=3\Delta x$, where $\Delta x$ is the
grid spacing (Capecelatro & Desjardins, 2013).

This computational framework was recently applied by the authors in a
configuration similar to the one considered here. In Shuai & Kasbaoui (2022),
the dynamics of a particle-laden Lamb-Oseen vortex at moderate Stokes numbers
are investigated using the Eulerian-Lagrangian methodology presented here.
Readers interested in further details about the numerical approach are
referred to this study.

###  2.2 Illustration of the instability

To illustrate the unstable dynamics activated by two-way coupling, we consider
a Rankine vortex at the circulation Reynolds number
$\mathrm{Re}_{\Gamma}=\Gamma/(2\pi\nu)=1000$, laden with mono-disperse
particles having diameter at the Stokes number
$\mathrm{St}_{\Gamma}=\tau_{p}/\tau_{f}=0.5$, average volume fraction
$\phi_{p}=1.2\,10^{-3}$, density ratio $\rho_{p}/\rho_{f}=830$, and mass
loading $M=1$. Here, $\Gamma$ corresponds to the vortex circulation, $r_{c}$
is the initial vortex core radius, $\tau_{f}=2\pi r_{c}^{2}/\Gamma$ is the
characteristic fluid velocity, $\tau_{p}=\rho_{p}d_{p}^{2}/(18\mu_{f})$ is the
particle response time, and $d_{p}$ is the particle diameter. The particles
are seeded randomly within the core region $r\leq r_{c}$ with velocities
matching the fluid velocity at their respective locations. We use the
Eulerian-Lagrangian framework described above on uniform Cartesian grid with
resolution $r_{c}/\Delta x=50$. The simulations performed here are termed
“pseudo-2D”, meaning that the axial direction $z$ is resolved with one grid
point and taken periodic with a thickness $\Delta z=3d_{p}$. This allows the
definition of volumetric quantities such as particle volume and volume
fraction. Despite considering pseudo-2D simulations, a large number of
particles must be tracker, here $N=376612$, due to the large scale ratio
$r_{c}/d_{p}=1400$. The remaining numerical parameters are as in Shuai &
Kasbaoui (2022).

In addition to two-way coupled simulations, we perform one-way coupled
simulations where the momentum exchange term (3) is purposely set to zero. In
this way, comparing one-way and two-way simulations reveals the effect of
particle feed-back on the flow dynamics.

Figure 1: Iso-contours of normalized vorticity ($\omega_{z}/\omega_{c}$) with
$\mathrm{St}_{\Gamma}=0.025$ and $\mathrm{Re}_{\Gamma}=1000$, for one-way
coupling and two-way coupling at five non-dimensional time $t/\tau_{f}$.

Figure 1 shows the evolution of the isocontours of axial vorticity
$\omega_{z}$ normalized by the initial center vorticity
$\omega_{c}=\omega_{z}(r=0,t=0)$. Successive snapshots are given between
$t/\tau_{f}=0$ and 48. When the particle feedback force is neglected, i.e., in
one-way coupling, we observe that the vorticity field remain axisymmetric at
all times. Although some viscous diffusion can be observed near the vorticity
jump $r\sim r_{c}$, the vorticity magnitude within the core remains mostly
flat as viscous effects are too small to cause significant diffusion within
the time frame $0\leq t/\tau_{f}\leq 48$. Contrary to the case of one-way
coupling, we observe significant distortion of the flow field when two-way
coupling is accounted for. The vorticity field quickly loses its cylindrical
symmetry due to the emergence of azimuthal perturbations. The latter grow into
vorticity filaments by $t/\tau_{f}\sim 12$ that are gradually convected away
from the vortex core. By $t/\tau_{f}\sim 48$, we observe several well-
established spiral arms emanating from the vortex core. In this region,
vorticity assumes a diffuse profile in contrast to the nearly flat profile
seen in one-way coupled simulations.

Figure 2: Iso-contours of normalized particle volume fraction with
$\mathrm{St}_{\Gamma}=0.025$ and $\mathrm{Re}_{\Gamma}=1000$, for one-way
coupling and two-way coupling at five non-dimensional time $t/\tau_{f}$.

Particle dispersion is also significantly impacted by the interaction between
the two phases. Figure 2 shows the iso-contours of particle volume fraction.
When the particle feedback is neglected, the particles gradually accumulate
into a ring-shaped cluster of diameter about $1.5\times r_{c}$. These
clustering dynamics were reported in the earlier work of Druzhinin (Druzhinin,
1994, 1995), Marshall (2005), and later by Shuai & Kasbaoui (2022). The
outward particle ejection is due to preferential concentration (Squires &
Eaton, 1990) whereby inertial particles are driven out of strongly vortical
regions. Like the carrier flow, the particle distribution is axisymmetric only
in one-way coupled simulations. In contrast to these prior observations, two-
way coupling leads to a loss of symmetry along with fast and wide dispersion
of the particles. The two-way coupled simulations in figure 2 reveal the
presence of azimuthal perturbations superimposed on the particle patch at
$t/\tau_{f}=12$. The perturbations extend into spiraling particle filaments
reminiscent of density-driven radial Rayleigh-Taylor instability observed by
Dixit & Govindarajan (2011) and Joly et al. (2005). These perturbations
destroy the ring structure observed in one-way coupled simulations.

The simulations shown in figure 1 and 2 suggest that semi-dilute inertial
particles trigger a modal instability. The latter cannot be attributed to
purely hydrodynamic nor purely granular effects since the Rankine vortex is
neutrally stable and particle-particle interactions are neglected. Rather, it
is an instability that arises from the two-way momentum exchange between the
two phases.

The instability observed here represents a distinct mechanism from the one
studied by Kasbaoui et al. (2015). There, the authors show that the interplay
between particle settling and preferential concentration activates a non-modal
instability in two-way coupled shear flows. It results in the formation of
sheets of concentrated particle clusters that rotate to progressively align
with the direction of the shear flow. Since particle settling is required for
perturbations to grow, the resulting growth rates depend on the gravitational
acceleration. In the present work, a modal instability arises regardless of
gravitational effects.

##  3 Linear stability analysis for weakly inertial particles

In this section, we establish a theoretical grounding to this novel two-phase
instability through rigorous linear stability analysis (LSA). The analysis is
intended to reveal the mechanisms that cause particle-laden vortex tube to be
unstable, and their functional dependence on the non-dimensional numbers at
hand. This study is restricted to the assumption of small, but non-zero
particle inertia such that $\mathrm{St}_{\Gamma}<0.1$.

In the following, we adopt an Eulerian-Eulerian description of the particle-
laden flow based on the Two-Fluid model (Marble, 1970; Druzhinin, 1995;
Jackson, 2000; Kasbaoui et al., 2015, 2019b). In this framework, the
conservation equations for the fluid phase and particle phase in the two-fluid
model read

$\displaystyle\nabla\cdot\bm{u}_{f}$ $\displaystyle=$ $\displaystyle 0,$ (8)
$\displaystyle\rho_{f}\left(\frac{\partial\bm{u}_{f}}{\partial
t}+\bm{u}_{f}\cdot\nabla\bm{u}_{f}\right)$ $\displaystyle=$
$\displaystyle-\nabla
p+\mu\nabla^{2}\bm{u}_{f}+\frac{\rho_{p}V_{p}n}{\tau_{p}}(\bm{u}_{p}-\bm{u}_{f}),$
(9) $\displaystyle\frac{\partial\rho_{p}n}{\partial
t}+\nabla\cdot\left(\rho_{p}n\bm{u}_{p}\right)$ $\displaystyle=$
$\displaystyle 0,$ (10)
$\displaystyle\frac{\partial\rho_{p}n\bm{u}_{p}}{\partial
t}+\nabla\cdot\left(\rho_{p}n\bm{u}_{p}\bm{u}_{p}\right)$ $\displaystyle=$
$\displaystyle-\frac{\rho_{p}V_{p}n}{\tau_{p}}(\bm{u}_{p}-\bm{u}_{f}),$ (11)

where $\bm{u}_{f}$ is the fluid velocity, $\bm{u}_{p}$ is the particle
velocity, $V_{p}$ is the particle volume and $n=\phi_{p}/(\rho_{p}V_{p})$ is
the particle number density. Equations (8) and (9) express mass and momentum
conservation for the fluid phase, respectively. Equations (10) and (11)
express mass and momentum conservation for the particle phase. In these
equations, it is assumed that the particles experience a hydrodynamic force
equal to Stokes drag. The two-phases are coupled through Stokes drag as seen
in equations (9) and (11). If the particle feedback force is dropped from the
fluid momentum equation (9), the evolution of the fluid becomes decoupled from
that of the particles. The fluid evolves as in single-phase, and the flow may
not become unstable as shown in §2 and Michalke & Timme (1967). The momentum
exchange between the two phases is critical for the development of an
instability.

Compared to the Eulerian-Lagrangian framework presented in §2, the form of the
differential equations in the Two-Fluid model is more suitable for theoretical
analysis using LSA. Still the two approaches yield qualitatively and
quantitatively similar evolution of particle-laden flows in the semi-dilute
regime provided that particle inertia is small. Further details on the
comparison of the two formulations can be found in (Kasbaoui et al., 2019a).
Thus, we restrict the analysis to cases where $\mathrm{St}_{\Gamma}\ll 1$.

Under the assumption of weakly inertial particles, it is possible to further
simplify the governing equations by adopting the fast equilibrium
approximation (see work by Balachandar and coworkers (Ferry & Balachandar,
2002, 2001; Rani & Balachandar, 2003)). Provided that particle inertia is low
enough, it is possible to express the particle velocity field as

$\displaystyle\bm{u}_{p}=\bm{u}_{f}-\tau_{p}\left(\frac{\partial\bm{u}_{f}}{\partial
t}+\bm{u}_{f}\cdot\nabla\bm{u}_{f}\right)+O(\tau_{p}^{2}),$ (12)

where it is seen that the particles deviate from the fluid streamlines by a
small inertial correction. Owing to the preferential concentration mechanism,
particles suspended in a vortex tube do not follow the closed-loop streamlines
of the carrier flow. Instead, they tend to migrate outward radially at a rate
controlled by their Stokes number.

Next, we combine equation (12) with the conservation equations (8), (9), (10),
and (11). Using $r_{c}$, $\Gamma/(2\pi r_{c})$ as the reference length and
velocity scales, respectively, the following non-dimensionalized equations are
obtained

$\displaystyle\nabla\cdot\bm{u}^{*}$ $\displaystyle=$ $\displaystyle 0,$ (13)
$\displaystyle(1+Mn^{*})\left(\frac{\partial\bm{u}^{*}}{\partial
t^{*}}+\bm{u}^{*}\cdot\nabla\bm{u}^{*}\right)$ $\displaystyle=$
$\displaystyle-\nabla
p^{*}+\frac{1}{\mathrm{Re}_{\Gamma}}\nabla^{2}\bm{u}^{*},$ (14)
$\displaystyle\frac{\partial n^{*}}{\partial
t^{*}}+\left(\bm{u}^{*}-\mathrm{St}_{\Gamma}\left(\frac{\partial\bm{u}^{*}}{\partial
t^{*}}+\bm{u}^{*}\cdot\nabla\bm{u}^{*}\right)\right)\cdot\nabla n^{*}$
$\displaystyle=$
$\displaystyle\mathrm{St}_{\Gamma}n^{*}\nabla\bm{u}^{*}:\nabla\bm{u}^{*},$
(15)

where the stared variables denote non-dimensional quantities,
$\mathrm{Re}_{\Gamma}=\Gamma/(2\pi\nu)$ is the circulation Reynolds number,
$\mathrm{St}_{\Gamma}=\tau_{p}r_{c}^{2}/\Gamma$ is the circulation Stokes
number, and $M=\rho_{p}V_{p}n_{0}/\rho_{f}$ is the mass loading. For ease of
notation, we drop the stars in the rest of the analysis.

In this form, the particle phase is coupled with the fluid through the
preferential concentration term appearing in the right hand side of (15).
Since Squires & Eaton (1990), preferential concentration term has been widely
studied (Squires & Eaton, 2014, 1991; Wang & Maxey, 1993; Batchelor & Nitsche,
1991; Druzhinin, 1995). It is understood that this term relates to
compressibility effects of the particle phase since
$\nabla\cdot\bm{u}_{p}\simeq-\tau_{p}\nabla\bm{u}_{f}:\nabla\bm{u}_{f}$ for
low inertia particles. This term may be further expressed as
$\nabla\bm{u}_{f}:\nabla\bm{u}_{f}=S^{2}-R^{2}$, where $S^{2}$ and $R^{2}$
represent the second invariants of the fluid strain rate tensor
$S=(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})/2$ and rotation tensor
$R=(\nabla\mathbf{u}-\nabla\mathbf{u}^{T})/2$ (Squires & Eaton, 2014). When
the local strain exceeds the local rotation, this term is positive leading to
particle accumulation. Conversely, when this term is negative, the particles
are expelled. As a result, inertial particles are progressively depleted from
rotational regions, such as vortex cores. The formation of a ring cluster
shown in figure 2 is a direct effect of preferential concentration.

The number density term in (14) suggests that the particles may lead to non-
Boussinesq dynamics similar to those observed in stratified fluids. In
particular, Dixit & Govindarajan (2010) showed that density gradients could
lead to a destabilization of 2D vortex. Sipp et al. (2005b) showed that
density gradients produce two distinct kinds of instability, namely
centrifugal instability (CTI) which mainly affects axisymmetric eigenmodes and
Rayleigh–Taylor instability (RTI) which mainly affects non-axisymmetric
eigenmodes. If particle inertia is negligible, one may ignore the preferential
concentration term, and treat the suspension as a fluid with effective density
$\rho_{\mathrm{eff}}=\rho_{f}+\rho_{p}V_{p}n$. In such case, the suspension
would be susceptible to the instabilities of variable density vortex flows
aforementioned. These effects are revisited for the case of non-inertial
particles in §3.3, and then extended to the case of weakly inertial particles
in §3.4.

###  3.1 Base state

(a)

(b)

Figure 3: Figure showing the configuration of rankine vortex flow and particle
distribution

In the remainder of the manuscript, we perform a linear stability analysis for
a base state initially comprising monodisperse particles seeded randomly
within a disk of non-dimensional size $r_{p}$ and a Rankine vortex. A
schematic of the configuration is given in figure 3(a). The non-dimensional
base fluid azimuthal velocity field $u^{b}_{f,\theta}$ and number density
fields $n^{b}$ are given by

$\displaystyle u^{b}_{f,\theta}(r,t=0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lr}r,&r\leq 1\\\
1/r,&r>1\end{array}\right.$ (18) $\displaystyle n^{b}(r,t=0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l l}1,&\mbox{if }r\leq r_{p}\\\
0,&\mbox{otherwise.}\end{array}\right.$ (21)

Note that the task of defining a base state requires careful consideration,
when the suspended particles are inertial. Even when the base flow is steady,
the particle distribution may be unsteady due to the outward ejection of
particles caused by the preferential concentration term in equation (15). The
latter is plotted in figure 3(b) for the Rankine vortex in equation (18). The
preferential concentration term is negative within the core ($r<1$) showing
that it acts as a sink in this region. Conversely, the preferential
concentration term is positive away from the core ($r>1$), corresponding to
regions where the particles will accumulate. Thus, the base number density
$n^{b}$ in (21) is expected to vary as time progresses.

The time evolution of the base state may be obtained by solving equations
(13), (14), and (15). Considering an inviscid and axisymmetric state
($\partial()/\partial\theta=0$), the equations dictating the evolution of the
base state are

$\displaystyle\frac{\partial u_{\theta}^{b}}{\partial t}=0,$ (22)
$\displaystyle\frac{\partial n^{b}}{\partial
t}+\frac{\mathrm{St}_{\Gamma}}{r}\frac{\partial}{\partial
r}\left((u_{\theta}^{b})^{2}n^{b}\right)=0.$ (23)

Following Druzhinin (1994), it is possible to obtain an analytical solution
for the unsteady number density by applying the method of characteristics.
Denoting the initial particle number density in (21) by
$n_{0}^{b}=n^{b}(r,0)$, the unsteady number density is written as

$\displaystyle n^{b}(r,t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l
l}e^{-2\>\mathrm{St}_{\Gamma}\>t}\;n_{0}^{b}\left(re^{-\mathrm{St}_{\Gamma}\;t}\right),&\mbox{if
}r\leq 1\\\ \frac{\displaystyle
r^{2}}{\displaystyle\sqrt{r^{4}-4\mathrm{St}_{\Gamma}\>t}}\;n_{0}^{b}\left((r^{4}-4\mathrm{St}_{\Gamma}\>t)^{1/4}\right),&\mbox{if
}r>1\,\,\mbox{and}\,\,t\leq t_{\mbox{cr}}\\\
\left(rr_{0}\right)^{2}\;n_{0}^{b}\left(r_{0}\right),&\mbox{if
}r>1\,\,\mbox{and}\,\,t>t_{\mbox{cr}}\end{array}\right.$ (27)

where $t_{\mbox{cr}}=(r^{4}-1)/4\mathrm{St}_{\Gamma}$ and $r_{0}(r,t)$ is
obtained from the equation $4\log
r_{0}+r_{0}^{4}=r^{4}-4\mathrm{St}_{\Gamma}t$. Despite the explicit time-
dependence of the base number density field $n^{b}$, a quasi-steady assumption
may be employed for weakly inertial particles. Under this assumption, the base
state is assumed frozen at $t=0$. This assumption is justified provided that
perturbations grow on a time scale that is much faster than the characteristic
evolution time of the number density field. From solution (27), it can be seen
that the timescale for the development of number density inhomogeneity is
$\tau_{c}=\tau_{f}/\mathrm{St}_{\Gamma}$. Shuai & Kasbaoui (2022) showed from
Eulerian-Lagrangian simulations that the time required to eject most particles
from the core of a vortex is $\sim 3\tau_{c}$. This timescale may be
significantly larger than the flow inertial timescale $\tau_{f}$ when
$\mathrm{St}_{\Gamma}\ll 1$. For such weakly inertial particles with low
inertia, particle clustering is slow compared to inertial effects of the
carrier vortex flow. Since we expect the instability to develop on a timescale
comparable with the flow inertial timescale, the assumption of quasi-steady
number density is justified.

###  3.2 Linearized equations

We assume that the base state discussed previously is subject to small
perturbations $n^{\prime}$ and $\bm{u}^{\prime}$ such that the total number
density and fluid velocity are $n=n^{b}+n^{\prime}$ and
$\bm{u}=\bm{u}^{b}+\bm{u}^{\prime}$. Taking the inviscid limit
($\mathrm{Re}_{\Gamma}\rightarrow\infty$) and linearizing equations (13) and
(14), the fluid perturbation is subject to the following equations

$\displaystyle\frac{1}{r}\frac{\partial(ru_{r}^{\prime})}{\partial
r}+\frac{1}{r}\frac{\partial u_{\theta}^{\prime}}{\partial\theta}$
$\displaystyle=$ $\displaystyle 0,$ (28)
$\displaystyle(1+Mn^{b})\left(\frac{\partial u_{r}^{\prime}}{\partial
t}+\frac{u^{b}_{\theta}}{r}\frac{\partial
u_{r}^{\prime}}{\partial\theta}-2\frac{u^{b}_{\theta}u^{\prime}_{\theta}}{r}\right)-Mn^{\prime}\frac{(u^{b}_{\theta})^{2}}{r}$
$\displaystyle=$ $\displaystyle-\frac{\partial p^{\prime}}{\partial r},$ (29)
$\displaystyle(1+Mn^{b})\left(\frac{\partial u_{\theta}^{\prime}}{\partial
t}+\frac{u^{b}_{\theta}}{r}\frac{\partial
u_{\theta}^{\prime}}{\partial\theta}-u^{\prime}_{r}\frac{\partial
u^{b}_{\theta}}{\partial r}+\frac{u^{\prime}_{r}u^{b}_{\theta}}{r}\right)$
$\displaystyle=$ $\displaystyle-\frac{1}{r}\frac{\partial
p^{\prime}}{\partial\theta},$ (30)

where equation (28) represents mass conservation in cylindrical coordinates,
and equations (29) and (30) represent the radial and azimuthal fluid momentum
conservation respectively. Introducing the base state angular velocity
$\Omega^{b}=u^{b}_{\theta}/r$ and axial vorticity
$\omega^{b}_{z}=(\partial(ru^{b}_{\theta})/\partial r)/r$, equations (29) and
(30) become

$\displaystyle(1+Mn^{b})\left(\left(\frac{\partial}{\partial
t}+\Omega^{b}\frac{\partial}{\partial\theta}\right)u^{\prime}_{r}-2\Omega^{b}u^{\prime}_{\theta}\right)-Mn^{\prime}r(\Omega^{b})^{2}$
$\displaystyle=$ $\displaystyle-\frac{\partial p^{\prime}}{\partial r},$ (31)
$\displaystyle(1+Mn^{b})\left(\left(\frac{\partial}{\partial
t}+\Omega^{b}\frac{\partial}{\partial\theta}\right)u^{\prime}_{\theta}+\omega_{z}^{b}u^{\prime}_{r}\right)$
$\displaystyle=$ $\displaystyle-\frac{1}{r}\frac{\partial
p^{\prime}}{\partial\theta}.$ (32)

The evolution of the perturbation axial vorticity, i.e.,

$\displaystyle\omega^{\prime}_{z}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}(ru^{\prime}_{\theta})-\frac{1}{r}\frac{\partial
u^{\prime}_{r}}{\partial\theta},$ (33)

is obtained by combing equations (31) and (32),

$\displaystyle\left(\frac{\partial}{\partial
t}+\Omega^{b}\frac{\partial}{\partial\theta}\right)\omega^{\prime}_{z}+\frac{\partial\omega^{b}_{z}}{\partial
r}u^{\prime}_{r}$ $\displaystyle=$ (34)
$\displaystyle-\frac{M}{1+Mn^{b}}\left((\Omega^{b})^{2}\frac{\partial
n^{\prime}}{\partial\theta}+\frac{\partial n^{b}}{\partial
r}\left(\left(\frac{\partial}{\partial
t}+\Omega^{b}\frac{\partial}{\partial\theta}\right)u^{\prime}_{\theta}+\omega^{b}_{z}u^{\prime}_{r}\right)\right).$

The number density perturbation is solution to the following equation,

$\displaystyle\frac{\partial n^{\prime}}{\partial t}+\Omega^{b}\frac{\partial
n^{\prime}}{\partial\theta}+\frac{\partial n^{b}}{\partial r}u^{\prime}_{r}$
$\displaystyle=$ (35) $\displaystyle\mathrm{St}_{\Gamma}\left(\frac{\partial
n_{b}}{\partial r}\left(\frac{\partial u^{\prime}_{r}}{\partial
t}+\Omega^{b}\frac{\partial
u^{\prime}_{r}}{\partial\theta}\right)-r(\Omega^{b})^{2}\frac{\partial
n^{\prime}}{\partial r}-2\Omega^{b}\frac{\partial n^{b}}{\partial
r}u^{\prime}_{\theta}\right)$
$\displaystyle+\mathrm{St}_{\Gamma}n^{b}\frac{2}{r}\left(\frac{\partial(r\Omega)}{\partial
r}\left(\frac{\partial
u^{\prime}_{r}}{\partial\theta}-u^{\prime}_{\theta}\right)-r\Omega^{b}\frac{\partial
u^{\prime}_{\theta}}{\partial
r}-\mathrm{St}_{\Gamma}\frac{1}{r}\frac{\partial}{\partial
r}((r\Omega^{b})^{2})n^{\prime}\right).$

Next, the perturbations $n^{\prime}$, $u_{\theta}^{\prime}$, and
$u_{r}^{\prime}$ are decomposed using Fourier series as

$\left[u_{r}^{\prime}(r,\theta,t),u_{\theta}^{\prime}(r,\theta,t),n^{\prime}(r,\theta,t)\right]=\left[\hat{u}_{r}(r,t),\hat{u}_{\theta}(r,t),\hat{n}(r,t)\right]e^{im\theta},$
(36)

where $m$ represents the mode number. Combining this form of perturbations
with equations (34) and (35), the set of linear equations become

$\displaystyle\mathcal{D}\left(\bar{n}r\mathcal{D}\left(r\frac{\partial\hat{u}_{r}}{\partial
t}\right)\right)-\bar{n}m^{2}\frac{\partial\hat{u}_{r}}{\partial
t}+im\Omega^{b}\left(\mathcal{D}\left(\bar{n}r\mathcal{D}(r\hat{u}_{r})\right)-\bar{n}m^{2}\hat{u}_{r}\right)$
$\displaystyle-
imr\mathcal{D}(\bar{n}\omega_{z}^{b})\hat{u}_{r}+rm^{2}(\Omega^{b})^{2}M\hat{n}=0,$
(37) $\displaystyle\left(\frac{\partial n}{\partial
t}-\mathrm{St}_{\Gamma}\mathcal{D}(n^{b})\frac{\partial\hat{u}_{r}}{\partial
t}\right)+im\Omega^{b}\left(n-\mathrm{St}_{\Gamma}\mathcal{D}(n^{b})\hat{u}_{r}\right)+\mathcal{D}(n^{b})\hat{u}_{r}$
$\displaystyle-i\;\mathrm{St}_{\Gamma}\left\\{\frac{i}{r}\mathcal{D}(r^{2}(\Omega^{b})^{2}\hat{n})-\frac{2}{m}\Omega^{b}\mathcal{D}(n^{b})\mathcal{D}(r\hat{u}_{r})\right.$
$\displaystyle\left.{}-\frac{2n^{b}}{mr}\left[\mathcal{D}\left(\Omega^{b}r\mathcal{D}(r\hat{u}_{r})\right)-m^{2}\mathcal{D}(r\Omega^{b})\hat{u}_{r}\right]\right\\}=0.$
(38)

Here, $\mathcal{D}$ denotes the derivative with respect to $r$ and
$\overline{n}=1+Mn^{b}$. Note that equations (37) and (38) retain the two-way
coupling required for an instability.

Using the assumption of quasi-steady base state discussed above, we perform a
modal stability analysis with the base state being frozen in time by
considering perturbations that have an exponential dependence on time. Under
this consideration, the perturbations can be written as

$\left[\hat{u}_{r}(r,t),\hat{u}_{\theta}(r,t),\hat{n}(r,t)\right]=\left[u_{r}(r),u_{\theta}(r),n(r)\right]e^{\sigma
t}$ (39)

where $\sigma$ is the complex growth rate. Perturbations whose real part of
$\sigma$, is positive grow without bound. Thus, the condition for an unstable
base state is $\Re(\sigma)>0$. Injecting the forms (39) in equations (37) and
(38), we obtain

$\displaystyle(\sigma+im\Omega^{b})\left[\mathcal{D}\left(\bar{n}r\mathcal{D}(ru_{r})\right)-\bar{n}m^{2}u_{r}\right]-imr\mathcal{D}(\bar{n}\omega_{z}^{b})u_{r}+rm^{2}(\Omega^{b})^{2}Mn=0,$
(40)
$\displaystyle(\sigma+im\Omega^{b})\left(n-\mathrm{St}_{\Gamma}\mathcal{D}(n^{b})u_{r}\right)+\mathcal{D}(n^{b})u_{r}+\mathrm{St}_{\Gamma}\left\\{\frac{1}{r}\mathcal{D}(r^{2}(\Omega^{b})^{2}n)+\frac{2i}{m}\Omega^{b}\mathcal{D}(n^{b})\mathcal{D}(ru_{r})\right.$
$\displaystyle\left.{}+\frac{2in^{b}}{mr}\left[\mathcal{D}\left(\Omega^{b}r\mathcal{D}(ru_{r})\right)-m^{2}\mathcal{D}(r\Omega^{b})u_{r}\right]\right\\}=0.$
(41)

The above pair of equations can be treated as an eigenvalue problem to obtain
the stability parameters.

###  3.3 Limit of inertia-less particles $(\mathrm{St}_{\Gamma}=0)$ and non-
Boussinesq effects

As previously discussed in this section, with the absence of particle inertia
($\mathrm{St}_{\Gamma}=0$), the preferential concentration term drops out of
the governing equations. Also, the base state in the case of
$\mathrm{St}_{\Gamma}=0$ remains constant over time as shown in equations (22)
and (23). A simplified linear equations can be obtained as

$\displaystyle(\sigma+im\Omega^{b})\left[\mathcal{D}\left(\bar{n}r\mathcal{D}(ru_{r})\right)-\bar{n}m^{2}u_{r}\right]-imr\mathcal{D}(\bar{n}\omega_{z}^{b})u_{r}+rm^{2}(\Omega^{b})^{2}Mn=0,$
(42) $\displaystyle n=-\frac{\mathcal{D}(n^{b})}{(\sigma+im\Omega^{b})}u_{r}.$
(43)

The above linearized equations are identical to the equations for a vortex
with a radially stratified density written by Dixit & Govindarajan (2011).
Without the particle inertia, the problem reduces to that of a vortex with a
density stratification induced by particles. Following the approach of Dixit &
Govindarajan (2011), it is possible solve equations (42) and (43) analytically
and obtain a dispersion relation for a Rankine vortex. This is done by
computing the perturbation velocity and number density fields separately in
regions of $r<1$, $1<r<r_{p}$ and $r>r_{p}$. The constants subsequently
obtained are found using the continuity condition at the jump locations $r=1$
and $r=r_{p}$. The dispersion relation, thus obtained, is written as

$\sigma^{3}+a_{2}\sigma^{2}+a_{1}\sigma+a_{0}=0,$ (44)

where

$\displaystyle a_{2}=i\left(m+2mr_{p}^{-2}-1-\mathrm{At}\;r_{p}^{-2m}\right),$
(45) $\displaystyle
a_{1}=mr_{p}^{-2}\left[2-m\left(2+r_{p}^{-2}\right)-\mathrm{At}\left(r_{p}^{-2}-2r_{p}^{-2m}\right)\right],$
(46) $\displaystyle
a_{0}=im(1-m)r_{p}^{-4}\left[m+\mathrm{At}\left(1-r_{p}^{-2m}\right)\right].$
(47)

where $\mathrm{At}=M/(2+M)$ is the Atwood number. The latter is a non-
dimensional number that commonly arises in density stratified flows. Here,
considering the particle-fluid mixture as a fluid with effective density
$\rho_{\mathrm{eff}}=1+(\rho_{p}/\rho_{f})V_{p}n=1+Mn$, we see that in the
core $\rho_{\mathrm{eff}}(r)=1+M=\rho_{\mathrm{max}}$ for $r\leq 1$, while
$\rho_{\mathrm{eff}}=1=\rho_{\mathrm{min}}$ away from the core ($r\geq 1$).
Thus, the Atwood number is
$\mathrm{At}=(\rho_{\mathrm{max}}-\rho_{\mathrm{min}})/(\rho_{\mathrm{max}}+\rho_{\mathrm{min}})$
which measures the relative magnitude of the density jump between the inner
and outer regions.

Figure 4 (a) and (b) show contours of growth rates for wavenumbers $m=2$ and
$4$, respectively, in the $\mathrm{At}-r_{p}$ plane. The white space labelled
“S” represents the region of stability. Thus there exists a critical radius,
$r_{p}^{\textrm{cr}}$, such that instability would only be observed when the
density jump is located beyond $r_{p}^{\textrm{cr}}$. Using the property of
the cubic equation (44) one can obtain a relation for $r_{p}^{\textrm{cr}}$ as
a function of $m$ and $\mathrm{At}$, which on further expanding for
$\mathrm{At}\ll 1$ yields,

$\displaystyle
r_{p}^{\textrm{cr}}=r^{\textrm{Kelvin}}\left[1-\frac{3}{2}\left\\{\frac{\mathrm{At}}{4m^{2}}\left(\frac{m-1}{m}\right)^{m-1}\right\\}^{1/3}\right]+\mathcal{O}(\mathrm{At}^{2/3})$
(48)

where $r^{\textrm{Kelvin}}=\sqrt{m/(m-1)}$ is the critical radius
corresponding to the discrete Kelvin mode (Roy & Subramanian, 2014). As is
evident from the complete numerical solution and the above small $\mathrm{At}$
asymptotic expression, the region of stability shrinks as the wavenumber
increases, suggesting that higher modes are more unstable. In addition, higher
values of growth rates are achieved with increasing values of Atwood number
and with the value of $r_{p}$ closer to one. This suggests that the
configuration at $r_{p}=1$ with infinitely heavy core ($M\rightarrow\infty$)
is the most unstable configuration.

To further study the influence of the Atwood number on the critical
wavenumber, we investigate the specific case of $r_{p}=1$, for which the
dispersion relation reduces to

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\frac{i}{2}(1+\mathrm{At})-im\pm\frac{i}{2}\sqrt{(1+\mathrm{At})^{2}-4m\mathrm{At}}.$
(49)

For $r_{p}=1$, the third root of equation (44) corresponds to a mode rotating
with core angular frequency ($\sigma=-im$.) The expression (49) indicates the
critical wavenumber ($m_{\mathrm{cr}}$) for unstable system is

$\displaystyle
m_{\mathrm{cr}}=\frac{(1+\mathrm{At})^{2}}{4\mathrm{At}}=\frac{(1+M)^{2}}{(2+M)M}.$
(50)

The above expression shows that increasing the Atwood number, or equivalently
the mass loading, promotes the instability by making lower wavenumber modes
unstable. In the limit case where $M\rightarrow\infty$, corresponding to
$\mathrm{At}\rightarrow 1$, the critical wavenumber becomes
$m_{\mathrm{cr}}\rightarrow 1$ showing that all modes are unstable.
Conversely, when $M\rightarrow 0$, none of the modes are unstable since
$m_{\mathrm{cr}}\rightarrow\infty$. Thus, we recover the known behavior for an
unladen Rankine vortex ($M\rightarrow 0$).

From expression (49), the growth rate of unstable modes ($m\geq
m_{\mathrm{cr}}$) is given by

$\Re(\sigma)=\frac{1}{2}\sqrt{4m\mathrm{At}-(1+\mathrm{At})^{2}}\simeq\sqrt{m\mathrm{At}}\quad\text{for
}m\gg 1.$ (51)

From this relationship, we draw two conclusions. First, a naturally developing
instability, i.e., unforced, will tend to emerge at high wavenumbers since
$\Re(\sigma)$ increases with $m$. Viscous effects are expected to curb
$\Re(\sigma)$ beyond a certain mode number $m$. Because viscous effects are
not accounted for in the present formulation, it is not possible to estimate a
priori which mode would emerge naturally. The second conclusion is that the
growth rate increases with mass loading. The maximum growth rate obtained when
$M\rightarrow\infty$ ($\mathrm{At}\rightarrow 1$) is
$\Re(\sigma)_{\mathrm{max}}=\sqrt{m}$.

###  3.4 Effects of particle inertia

(a) $m=2$

(b) $m=4$

Figure 4: Contours of the analytically obtained growth rate (equation (44)) in
the $\mathrm{At}-r_{p}$ plane for a Rankine vortex with
$\mathrm{St}_{\Gamma}=0$.

Next, we consider the effects of non-zero particle inertia
($\mathrm{St}_{\Gamma}\neq 0$) on the stability of the particle-laden Rankine
vortex. For this, we first perform a modal stability analysis by solving
equations (40) and (41) with a frozen base state.

For the case where the Rankine vortex and particle core have equal radii,
i.e., $r_{p}=1$, a dispersion relation can be obtained analytically:

$\sigma^{2}-i\left(1+\mathrm{At}-2m+im\mathrm{At}\mathrm{St}_{\Gamma}\right)\sigma-m(m-1-i(m-2)\mathrm{At}\mathrm{St}_{\Gamma})=0,$
(52)

which yields

$\sigma=\frac{i}{2}(1+\mathrm{At}+i\mathrm{At}\mathrm{St}_{\Gamma}m)-im\pm\frac{i}{2}\sqrt{(1+\mathrm{At}+i\mathrm{At}\mathrm{St}_{\Gamma}m)^{2}-4m(\mathrm{At}+2i\mathrm{At}\mathrm{St}_{\Gamma})}$
(53)

Note that expression (53) collapses on (49) when
$\mathrm{St}_{\Gamma}\rightarrow 0$. The corresponding expressions for the
eigenfunctions of a mode $m$ are

$\displaystyle u_{r}$ $\displaystyle=$ $\displaystyle\left\\{\begin{array}[]{l
l}r^{(m-1)},&\mbox{if }r\leq 1\\\
r^{-(m+1)},&\mbox{otherwise}\end{array}\right.$ (56) $\displaystyle n$
$\displaystyle=$
$\displaystyle\left(-\mathrm{St}_{\Gamma}+\frac{1+2i\;\mathrm{St}_{\Gamma}}{\sigma+im}\right)\delta(r-1).$
(57)

Here, $\delta(r-1)$ denotes the Dirac delta function. These functions
correspond to perturbations that are concentrated at the number density and
vorticity jumps $r=r_{p}=1$. To further study the stabilizing/destabilizing
effect of particle inertia, we expand the expression for growth rate (equation
(53)) in a series of $\mathrm{St}_{\Gamma}$ as

$\displaystyle\sigma=\sigma\big{|}_{\mathrm{St}_{\Gamma}=0}-\frac{m\mathrm{St}_{\Gamma}\mathrm{At}}{2}\pm\left\\{\frac{m\mathrm{At}(3-\mathrm{At})}{2\sqrt{\Delta}}\mathrm{St}_{\Gamma}+\frac{im^{2}\mathrm{At}^{2}((m-2)\mathrm{At}+2)}{\Delta^{3/2}}\mathrm{St}_{\Gamma}^{2}\right\\}+\mathcal{O}(\mathrm{St}_{\Gamma}^{3}).$

where $\Delta=(1+\mathrm{At})^{2}-4m\mathrm{At}$. For
$\mathrm{St}_{\Gamma}=0$, $\Delta<0$ provides the condition for instability
for a fixed $m$, that is, a mode $m$ is stable if
$\mathrm{At}<\mathrm{At}_{\textrm{cr}}=2m-1-2\sqrt{m(m-1)}$. For
$\mathrm{St}_{\Gamma}\neq 0$, particle inertia can play contrasting roles. For
$\mathrm{At}<\mathrm{At}_{\textrm{cr}}$, the dispersed phase destabilizes the
vortex

$\displaystyle\Re(\sigma)=\frac{m\mathrm{At}\mathrm{St}_{\Gamma}}{2}\left\\{\frac{3-\mathrm{At}}{\sqrt{\Delta}}-1\right\\}+\mathcal{O}(\mathrm{St}_{\Gamma}^{3}),$
(59)

while for $\mathrm{At}>\mathrm{At}_{\textrm{cr}}$, the dispersed phase
stabilizes the vortex

$\displaystyle\Re(\sigma)=\frac{1}{2}\sqrt{|\Delta|}-\frac{m\mathrm{St}_{\Gamma}\mathrm{At}}{2}+\mathcal{O}(\mathrm{St}_{\Gamma}^{2}).$
(60)

Figure 5 (a) shows the variation of growth rate plotted over a range of Stokes
numbers for wavenumber $m=2$. With the smaller value of Atwood number
$\mathrm{At}=0.1$, the growth rates decrease with increasing Stokes numbers.
However, the trends get reversed for higher Atwood numbers, with the growth
rates increasing with increasing Stokes numbers, confirming the predictions
from the asymptotic calculations (equations (59) and (60)). Figure 5 (b) plots
the growth rates over a range of Atwood numbers for three values of Stokes
numbers. It is found that the growth rate is always positive as long as the
Atwood number $\mathrm{At}$ and Stokes number $\mathrm{St}_{\Gamma}$ are non-
zero. This indicates that destabilization of the vortex by the particles is
guaranteed as long as the particle inertia is non-zero.

For cases where $r_{p}>1$, we solve equations (40) and (41) numerically using
Chebyshev spectral collocation method (Trefethen, 2000). To handle the
discontinuities, we smooth the base state using hyperbolic tangents of
thickness $\Delta=0.01$. Figure 6 (a) and (b) show the contours of growth rate
in $\mathrm{At}-r_{p}$ plane for Stokes numbers $0.01$ and $0.1$ respectively.
Unlike the case of $\mathrm{St}_{\Gamma}=0$ (see figure 4), there is no region
of stability indicating that mode $m=2$ is always unstable when the particles
are inertial. Further, the particle-laden vortex is unstable even for values
of $r_{p}$ larger than 1. Whereas the results for non-inertial particles are
similar to those obtained by Dixit & Govindarajan (2011), inclusion of
particle inertia breaks the equivalence with a radially startified Rankine
vortex. For the case $\mathrm{St}_{\Gamma}=0.01$, the optimal value of $r_{p}$
leading to the largest growth rate is slightly in excess of 1 for low values
of $\mathrm{At}$, and approaches $r_{p}=1$ as $\mathrm{At}$ increases. For the
larger inertia case $\mathrm{St}_{\Gamma}=0.1$, $r_{p}=1$ leads to the largest
growth rates for all values of $\mathrm{At}$.

(a)

(b)

Figure 5: Growth rates $\Re(\sigma)$ for a Rankine vortex with $r_{p}=1$ and
$m=2$ obtained analytically (equation (52)), (a)
$\Re(\sigma)-\Re(\sigma)|_{\mathrm{St}_{\Gamma}=0}$ plotted over a range of
Stokes numbers for $\mathrm{At}=0.1$ (—), $\mathrm{At}=0.3$ (- - -),
$\mathrm{At}=0.9$ ($\cdots$); (b) $\Re(\sigma)$ plotted over a range of Atwood
numbers for $\mathrm{St}_{\Gamma}=0$ (red), $\mathrm{St}_{\Gamma}=0.01$
(black) and $\mathrm{St}_{\Gamma}=0.1$ (blue).

(a) $\mathrm{St}_{\Gamma}=0.01$

(b) $\mathrm{St}_{\Gamma}=0.1$

Figure 6: Plot of the growth rate $\Re(\sigma)$ for a vortex using the smooth
profiles with $\Delta=0.01$ for varying $\mathrm{At}$ and $r_{p}$ for $m=2$.

##  4 Comparison of Linear Stability Analysis with Euler-Lagrange simulations

In order to validate the linear stability analysis, we compare the analytical
results with those from Eulerian-Lagrangian simulations of the two-way coupled
system. To make the comparison tractable and accessible, the case $r_{p}=1$ is
used as the benchmark since the growth rates can be computed analytically from
equation (52).

Table 1 lists a summary of non-dimensional parameters for 8 cases presently
considered. In all these simulations, the Stokes number, circulation Reynolds
number $\mathrm{Re}_{\Gamma}$ and ratio $r_{p}/r_{c}$ are fixed at $0.002$,
$5000$ and $1$, respectively. The Stokes number is sufficiently small such
that the initial state can be considered quasi-steady and comparisons with LSA
are permissible. In dimensional terms, the carrier fluid has density
$\rho_{f}=1.2\;\mathrm{kg/m^{3}}$ and viscosity $\mu_{f}=1.8\times
10^{-5}\;\mathrm{kg\cdot m^{-1}\cdot s^{-1}}$. The initial vortex core is
$r_{c}=0.385\ \mathrm{m}$ and circulation is $\Gamma=9.42\times
10^{-2}\;\mathrm{m^{2}\cdot s^{-1}}$. Particles are seeded with density
$\rho_{p}=1000\;\mathrm{kg/m^{3}}$ and average volume fraction
$\langle\phi_{p}\rangle=1.2\times 10^{-3}$. In all cases chosen, the ratio
$r_{c}/d_{p}$ is large (of the order $10^{4}$), such that fluctuations due to
the discrete particle forcing fall well below the viscous dissipation scale.
In this way, the coupling between the fluid and particles occurs primarily
through collective particle dynamics at scales comparable with $r_{c}$, rather
than discrete effects at scales comparable with the inter-particle distance.

The initial conditions of the numerical simulations correspond to a
superposition of the base Rankine vortex and the perturbation eigenmodes (56)
and (57). In order to capture the linear regime described by the linear
stability analysis, the perturbations are initialized with a small amplitude
$\epsilon=0.01$. The Dirac delta function that appears in the eigenfunctions
(56) and (57) is discretized on the Cartesian grid according to the radial
distance $l_{i}$ (i=1,2,3,4) between the vortex core and the four vertices of
each cell. Cells that are intersected by the discontinuity ($\max(l_{i})>1$
and $\min(l_{i})<1$) have $\delta(r-1)=1/(\Delta x\Delta y)$, whereas
$\delta(r-1)=0$ everywhere else. The Lagrangian particles are generated
randomly within each cell based on the number density computed with the
aforementioned discrete perturbation. The particle velocities are set to match
the fluid velocity at their locations. Figures 7(a) and 7(b) show number
density and axial vorticity isocontours of a representative perturbation with
mode $m=4$ and disturbance amplitude $\epsilon=0.01$ .

Case |  $m$ |  $\mathrm{Re}_{\Gamma}$ |  $\mathrm{At}$ |  $M$ |  $\langle\phi_{p}\rangle$ |  $r_{c}/d_{p}$
---|---|---|---|---|---|---
1 |  2 |  $5000$ |  1/3 |  1 |  $1.2\;10^{-3}$ |  10760
2 |  3 |  $5000$ |  1/3 |  1 |  $1.2\;10^{-3}$ |  10760
3 |  4 |  $5000$ |  1/3 |  1 |  $1.2\;10^{-3}$ |  10760
4 |  5 |  $5000$ |  1/3 |  1 |  $1.2\;10^{-3}$ |  10760
5 |  6 |  $5000$ |  1/3 |  1 |  $1.2\;10^{-3}$ |  10760
6 |  3 |  $5000$ |  1/6 |  2/5 |  $1.2\;10^{-3}$ |  10760
7 |  3 |  $5000$ |  1/5 |  1/2 |  $1.2\;10^{-3}$ |  10760
8 |  3 |  $5000$ |  1/4 |  2/3 |  $1.2\;10^{-3}$ |  10760
Table 1: Non-dimensional parameters considered in simulations.

(a)

(b)

Figure 7: Example of initial perturbation in the Eulerian-Lagrangian
simulations. Normalized number density (a) and axial vorticity (b)
perturbations based on eigenmodes (56) and (57) with mode number $m=4$.
Figure 8: Iso-contours of particle volume fraction normalized by the average
volume fraction for cases 1-5 at five instant time. As time progresses, the
inertial particles are expelled out of the vortex cores, forming cluster arms
of which the number equals to the mode number. The small perturbations
continue growing and ultimately destruct the structure. Figure 9: Time
evolution of the normalized perturbation kinetic energy in Eulerian-Lagragian
simulations of a Rankine vortex superimposed with eigenmodes $m=$ 2 – 6.

(a)

(b)

Figure 10: Comparison of instability growth rates between linear stability
analysis (LSA) and Eulerian-Lagrangian (E-L) simulations at
$\mathrm{St}_{\Gamma}=0.002$, $\mathrm{Re}_{\Gamma}=5000$, and $M=1$. (a)
Growth rates for mode numbers $m=$2 – 6 at $\mathrm{At}=1/3$ (LSA:●, EL:○) and
$\mathrm{At}=1/5$ (LSA:▴, EL:▵). (b) Growth rates for mode $m=3$ at various
Atwood numbers (LSA:●, EL:○).

Figure 8 shows the time evolution of the particle volume fraction between
$t=0$ and $t=60$ for modes $m=2-6$. During $0\leq t\leq 15$, the perturbed
modes grow linearly resulting in the deformation of the interface at the
number density jump. The latter, which is initially circular, develops
azimuthal perturbations depending on which mode has been triggered initially.
Notably, we observe the gradual appearance of high wavenumber modes with
$m\sim 24-26$ in addition to the low wavenumbers perturbed manually. High
wavenumber perturbations are expected because the discretization of the
eigenfunctions and discrete Lagrangian particles introduce high wavenumber
fluctations. As shown by the linear stability analysis, the perturbation
growth rate increases with wavenumber $m$, such that it may be anticipated
that such modes will eventually dominate. Indeed, the high wavenumbers
continue growing faster until they enter a non-linear stage by $t=30$. These
modes develop into long spiraling arms of inertial particles that expand
radially outward. Regions of unladen fluid in between the particle arms are
sucked inward forming mushroom-like structures seen at $t=45$ and $t=60$.

These dynamics may be related to the misalignment between the particle density
interface and centripetal acceleration field (Joly et al., 2005), which is
similar to the mechanism of wave interaction between density interface and
vorticity interface Dixit & Govindarajan (2011). The baroclinic torque caused
by the misalignment concentrates vorticity on the spiral arms and destroys it
in between, resulting in radial filaments.

In order to determine the perturbation growth rates from the present
simulations, we compute the perturbation kinetic energy

$E=\iint\frac{1}{2}(\bm{u}-\bm{u}_{f}^{b})\cdot(\bm{u}-\bm{u}_{f}^{b})dS$ (61)

where $\bm{u}_{f}^{b}$ is the base state velocity corresponding to a Rankine
vortex. Figure 9 shows the evolution of perturbation energy in logarithm scale
for simulations with perturbed modes $m=2$ –6. After a short early transient
adjustment of the flow, we observe an exponential growth of the perturbation
energy from $t\sim 1.5$ to $t\sim 3.0$. The perturbation energy is seen to
increase with increasing mode number $m$, reflecting the larger growth rate
for higher mode numbers. The exponential growth eventually tapers off and the
perturbation kinetic energy saturates. The deviation from exponential growth
correlates with the emergence of non-linear effects in the flow and number
density fields as seen in figure 8.

The instability growth rate is obtained from Eulerian-Lagragian simulations
using the relation

$\Re(\sigma)=\frac{1}{2E}\frac{dE}{dt},$ (62)

Perturbation growth rates are obtained from least squares regression between
$t=1.5$ and $t=3$. Figure 10(a) shows comparisons between growths rates
obtained from Euler-Lagrange simulations (equation (62)) and those predicted
by the linear stability analysis. Figure 10(a) shows the results for Atwood
numbers $\mathrm{At}=1/3$ and 1/5 and mode numbers $m=$2 – 6. Both Euler-
Lagrange simulations and linear stability analysis capture the growth rate
increase with increasing mode number $m$. The growth rates for mode $m=3$ and
varying Atwood number are shown in figure 10(b). Both Euler-Lagrange and
linear stability analysis show that the instability growth rate increases with
Atwood number $\mathrm{At}$. The growth rates predicted by the full numerical
simulation has excellent agreement with linear stability predictions.

While the growth rate obtained from numerical simulations and linear stability
analysis are within close quantitative agreement, some of the differences can
be related to the discretization of the eigenmodes and Lagrangian
representation of the particles in the numerical simulations. Whereas
eigenmodes derived from linear stability analysis require strict
discontinuities in number density and vorticity, these effects are necessarily
smoothed out when discretized on a mesh. Further results from the linear
stability analysis are derived on the premise that viscous effects and
discrete particle effects are negligible. However, these effects cannot be
ruled out in the Euler-Lagrange simulations where viscosity still plays a role
despite the large $\Rey_{\Gamma}=5000$ considered, and Lagrangian particles
trigger high wavenumber modes. Despite the differences between the linear
stability analysis and the Euler-Lagrange model, the agreement between the two
approaches is satisfactory.

##  5 Conclusion

In this work, we have shown that introducing heavy particles in the core of a
vortex causes the breakdown of the flow structure due to an instability
activated by two-way coupling, i.e., momentum exchange between the two phases.
If particle feedback is neglected (one-way coupling), no instability is
observed: the vortex retains its coherent axisymmetric structure while
inertial particles are slowly expelled outward forming a ring of clustered
particles around the core. The inclusion of two-way coupling breaks the
axisymmetry and gives rise to azimuthal perturbations in both number density
and vorticity fields. As these perturbations develop, they turn into spiraling
arms of concentrated particles emanating out of the core while regions of
particle-free flow are sucked inward. The vorticity field displays similar
pattern which cause the breakdown of the initial Rankine structure.
Remarkably, this breakdown occurs even for inertia-less particles, i.e., zero-
Stokes number, provided that the mass loading is not vanishingly small.
Particle inertia plays a dual role, in that, it destabilizes low wavenumber
modes, but reduces the growth rate of high wavenumber ones.

The mechanisms driving the instability are characterized using linear
stability analysis of a base state consisting of a Rankine vortex with core
radius $r_{c}$ and a circular patch of particles with radius $r_{p}$. To
describe the two-phase flow, we use an inviscid Two-Fluid model, and assume
that the particles are weakly inertial, i.e., their Stokes number
$\mathrm{St}_{\Gamma}$ is small but not zero. This assumption has two merits.
First, it allows us to approximate the particle velocity field in terms of the
carrier fluid velocity field and an inertial correction, thus, reducing the
number of equations to solve. Second, it justifies a quasi-steady approach in
which the time-dependence of the base state is ignored enabling us to carry
out a classic linear stability analysis. This step is justified by the fact
that variations of the base state have a timescale
$\tau_{c}=\tau_{f}/\mathrm{St}_{\Gamma}$ which becomes very large in the limit
$\mathrm{St}_{\Gamma}\ll 1$ (Shuai & Kasbaoui, 2022).

Analysis of growth rates obtained from linear stability analysis reveals the
existence of unstable modes regardless of the value of the Stokes number
$\mathrm{St}_{\Gamma}$ and as long as the mass loading $M\neq 0$. Generally,
the most unstable configuration corresponds to matching particle patch and
vortex core radii, i.e., $r_{c}=r_{p}$. For $\mathrm{St}_{\Gamma}=0$
particles, the problem becomes analogous to a vortex with discontinuous radial
density stratification described by Dixit & Govindarajan (2011). Modes with
wavenumber $m$ above a critical wavenumber
$m_{\mathrm{cr}}=(1+M)^{2}/((2+M)M)$ are unstable and have a growth rate
$\Re\\{\sigma\\}\sim\sqrt{m\mathrm{At}}$ when $m\gg 1$, and where the Atwood
number $\mathrm{At}=M/(2+M)$ provides a measure of the relative magnitude of
the density jump between inner (particle-rich) and outer (particle-free)
regions ($0\leq\mathrm{At}\leq 1$). When the particles are inertial
($\mathrm{St}_{\Gamma}\neq 0$) all modes become unstable regardless
$\mathrm{St}_{\Gamma}$ and $M\neq 0$, even those with $m<m_{\mathrm{cr}}$.
However, the growth rate of modes $m\geq m_{\mathrm{cr}}$ is reduced compared
to the case $\mathrm{St}_{\Gamma}=0$. For large wavenumber modes $m\gg 1$, the
growth rate approaches
$\Re\\{\sigma\\}\sim\sqrt{m\mathrm{At}}(1-2\mathrm{St}_{\Gamma}\sqrt{m\mathrm{At}})+O(\mathrm{St}_{\Gamma}^{2})$.
For both inertial and non-inertial particles, the tendency of growth rates to
increase with wavenumber $m$ suggests that modes that naturally emerge in
experiments and simulations shall have a high-wavenumber. This mode is likely
determined by a balance between viscous effects and preferential
concentration, although this cannot be ascertained as we have not considered
viscous effects in our linear stability analysis.

In addition to the linear stability analysis, Eulerian-Lagrangian simulations
are carried out at $\mathrm{St}_{\Gamma}=0.002$, $\Rey=5000$, and mass loading
$M=$ 0.4 – 1. The simulations are initialized with the base state superimposed
with eigenmodes found from linear stability analysis. The simulations allow us
to explore the transition from linear to nonlinear regimes. During the early
evolution of the perturbations, growth rates computed from the perturbation
kinetic energy in the Eulerian-Lagrangian simulations show excellent agreement
with growth rates predicted by the linear stability analysis for modes $m=$ 2
– 6 and mass loadings $M=$ 0.4 – 1 (or $\mathrm{At}=$ 1/6 – 1/3). When the
flow enters a nonlinear stage, we observe the emergence of a high-wave number
mode with $m\sim 24$ – 26. This mode develops into spiraling arms of number
density and vorticity that ultimately cause the breakdown of the initial base
state.

One must note that the present study has a key difference from prior
investigations of dusty Kolmogorov flow and Rayleigh-Taylor turbulence – a
particle-free vortex is always stable. Our study emphasizes that the disperse
phase can also be the source of a new instability besides modifying an
existing one. The novel instability for a dusty vortex identified here
highlights how the feedback force from a disperse phase can induce a breakdown
of an otherwise resilient vortical structure.

Acknowledgements. The authors acknowledge support from IIT Madras for the
“Geophysical Flows Lab” research initiative under the Institute of Eminence
framework.

Declaration of Interests. The authors report no conflict of interest.

## References

  * Anderson & Jackson (1967) Anderson, T. B. & Jackson, R. 1967 Fluid Mechanical Description of Fluidized Beds. Equations of Motion. Industrial & Engineering Chemistry Fundamentals 6 (4), 527–539.
  * Balachandar & Eaton (2010) Balachandar, S. & Eaton, J. K. 2010 Turbulent dispersed multiphase flow. Annual review of fluid mechanics 42, 111–133.
  * Balmforth et al. (2012) Balmforth, N.J., Roy, A & Caulfield, C.P. 2012 Dynamics of vorticity defects in stratified shear flow. Journal of Fluid Mechanics 694, 292–331.
  * Batchelor & Nitsche (1991) Batchelor, G. K. & Nitsche, J. M. 1991 Instability of stationary unbounded stratified fluid. Journal of Fluid Mechanics 227, 357–391.
  * Bluestein et al. (2003) Bluestein, H. B., W., Lee, Bell, M., Weiss, C. C. & Pazmany, A. L. 2003 Mobile Doppler Radar Observations of a Tornado in a Supercell near Bassett, Nebraska, on 5 June 1999. Part II: Tornado-Vortex Structure. Monthly Weather Review 131 (12), 2968–2984.
  * Bluestein et al. (2004) Bluestein, H. B., Weiss, C. C. & Pazmany, A. L. 2004 Doppler Radar Observations of Dust Devils in Texas. Monthly Weather Review 132 (1), 209–224.
  * Broadbent & Moore (1979) Broadbent, E. G. & Moore, D. W. 1979 Acoustic destablilization of vortices. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 290 (1372), 353–371.
  * Brown et al. (2005) Brown, R. A., Flickinger, B. A., Forren, E., Schultz, D. M. & al, et 2005 Improved Detection of Severe Storms Using Experimental Fine-Resolution WSR-88D Measurements. Weather and Forecasting 20 (1), 3–10,12–14.
  * Capecelatro & Desjardins (2013) Capecelatro, J. & Desjardins, O. 2013 An Euler–Lagrange strategy for simulating particle-laden flows. Journal of Computational Physics 238, 1–31.
  * Carpenter et al. (2011) Carpenter, Jeffrey R, Tedford, Edmund W, Heifetz, Eyal & Lawrence, Gregory A 2011 Instability in stratified shear flow: Review of a physical interpretation based on interacting waves. Applied Mechanics Reviews 64 (6).
  * Dixit & Govindarajan (2010) Dixit, H. N & Govindarajan, R. 2010 Vortex-induced instabilities and accelerated collapse due to inertial effects of density stratification. Journal of Fluid Mechanics 646, 415–439.
  * Dixit & Govindarajan (2011) Dixit, H. N. & Govindarajan, R. 2011 Stability of a vortex in radial density stratification: Role of wave interactions. Journal of Fluid Mechanics 679, 582–615.
  * Druzhinin (1994) Druzhinin, O. A. 1994 Concentration waves and flow modification in a particle-laden circular vortex. Physics of Fluids 6 (10), 3276–3284.
  * Druzhinin (1995) Druzhinin, O. A. 1995 On the two-way interaction in two-dimensional particle-laden flows: The accumulation of particles and flow modification. Journal of Fluid Mechanics 297, 49–76.
  * Ferry & Balachandar (2001) Ferry, J. & Balachandar, S. 2001 A fast Eulerian method for disperse two-phase flow. International Journal of Multiphase Flow 27 (7), 1199–1226.
  * Ferry & Balachandar (2002) Ferry, J. & Balachandar, S. 2002 Equilibrium expansion for the Eulerian velocity of small particles. Powder Technology 125 (2–3), 131–139.
  * Ford (1994) Ford, R. 1994 The instability of an axisymmetric vortex with monotonic potential vorticity in rotating shallow water. Journal of Fluid Mechanics 280, 303–334.
  * Fung & Kurzweg (1975) Fung, Y. T. & Kurzweg, U. H. 1975 Stability of swirling flows with radius-dependent density. Journal of Fluid Mechanics 72 (02), 243.
  * Guazzelli & Morris (2011) Guazzelli, E. & Morris, J. F 2011 A physical introduction to suspension dynamics, , vol. 45. Cambridge University Press.
  * Jackson (2000) Jackson, Roy 2000 The Dynamics of Fluidized Particles. Cambridge University Press.
  * Joly et al. (2005) Joly, L., Fontane, J. & Chassaing, P. 2005 The Rayleigh–Taylor instability of two-dimensional high-density vortices. Journal of Fluid Mechanics 537 (-1), 415.
  * Kasbaoui (2019) Kasbaoui, M. H. 2019 Turbulence modulation by settling inertial aerosols in Eulerian-Eulerian and Eulerian-Lagrangian simulations of homogeneously sheared turbulence. Physical Review Fluids 4 (12), 124308.
  * Kasbaoui et al. (2019a) Kasbaoui, M. H., Koch, D. L. & Desjardins, O. 2019a Clustering in Euler–Euler and Euler–Lagrange simulations of unbounded homogeneous particle-laden shear. Journal of Fluid Mechanics 859, 174–203.
  * Kasbaoui et al. (2019b) Kasbaoui, M. H., Koch, D. L. & Desjardins, O. 2019b The rapid distortion of two-way coupled particle-laden turbulence. Journal of Fluid Mechanics 877, 82–104.
  * Kasbaoui et al. (2015) Kasbaoui, M. H., Koch, D. L. & Subramanian, G.and Desjardins, O. 2015 Preferential concentration driven instability of sheared gas–solid suspensions. Journal of Fluid Mechanics 770, 85–123.
  * Lamb (1993) Lamb, H. 1993 Hydrodynamics. Cambridge; New York: Cambridge University Press.
  * Le Dizès & Billant (2009) Le Dizès, S. & Billant, P. 2009 Radiative instability in stratified vortices. Physics of Fluids 21 (9), 096602\.
  * Magnani et al. (2021) Magnani, M., Musacchio, S. & Boffetta, G. 2021 Inertial effects in dusty rayleigh–taylor turbulence. Journal of Fluid Mechanics 926.
  * Marble (1970) Marble, F E 1970 Dynamics of Dusty Gases. Annual Review of Fluid Mechanics 2 (1), 397–446.
  * Marshall (2005) Marshall, J. S. 2005 Particle dispersion in a turbulent vortex core. Physics of Fluids 17 (2), 025104.
  * Maxey & Riley (1983) Maxey, M. R. & Riley, J. J. 1983 Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids 26 (4), 883–889.
  * Michalke & Timme (1967) Michalke, Alfons & Timme, Adalbert 1967 On the inviscid instability of certain two-dimensional vortex-type flows. Journal of Fluid Mechanics 29 (4), 647–666.
  * Paoli & Shariff (2016) Paoli, R. & Shariff, K. 2016 Contrail Modeling and Simulation. Annual Review of Fluid Mechanics 48 (1), 393–427.
  * Paoli et al. (2008) Paoli, R., Vancassel, X., Garnier, F. & Mirabel, P. 2008 Large-eddy simulation of a turbulent jet and a vortex sheet interaction: Particle formation and evolution in the near field of an aircraft wake. Meteorologische Zeitschrift 17 (2), 131–144.
  * Rani & Balachandar (2003) Rani, S. L. & Balachandar, S. 2003 Evaluation of the equilibrium Eulerian approach for the evolution of particle concentration in isotropic turbulence. International Journal of Multiphase Flow 29 (12), 1793–1816.
  * Ravichandran & Govindarajan (2015) Ravichandran, S. & Govindarajan, R. 2015 Caustics and clustering in the vicinity of a vortex. Physics of Fluids 27 (3), 033305.
  * Roy et al. (2021) Roy, A., Garg, P., Reddy, J. S. & Subramanian, G. 2021 Inertio-elastic instability of a vortex column. arXiv preprint arXiv:2101.00805 .
  * Roy & Subramanian (2014) Roy, A. & Subramanian, G. 2014 Linearized oscillations of a vortex column: the singular eigenfunctions. Journal of Fluid Mechanics 741, 404–460.
  * Saffman (1995) Saffman, P. G. 1995 Vortex dynamics. Cambridge university press.
  * Schecter & Montgomery (2004) Schecter, D. A. & Montgomery, M. T. 2004 Damping and pumping of a vortex rossby wave in a monotonic cyclone: critical layer stirring versus inertia–buoyancy wave emission. Physics of Fluids 16 (5), 1334–1348.
  * Shuai & Kasbaoui (2022) Shuai, S. & Kasbaoui, M. Houssem 2022 Accelerated decay of a Lamb-Oseen vortex tube laden with inertial particles in Eulerian-Lagrangian simulations. Journal of Fluid Mechanics 936.
  * Sipp et al. (2005a) Sipp, D., Fabre, D., Michelin, S. & Jacquin, L. 2005a Stability of a vortex with a heavy core. Journal of Fluid Mechanics 526, 67–76.
  * Sipp et al. (2005b) Sipp, D., Fabre, D., Michelin, S. & Jacquin, L. 2005b Stability of a vortex with a heavy core. Journal of Fluid Mechanics 526, 67–76.
  * Sozou & Lighthill (1987) Sozou, C. & Lighthill, M. J. 1987 Adiabatic perturbations in an unbounded Rankine vortex. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 411 (1840), 207–224.
  * Sozza et al. (2022) Sozza, A., Cencini, M., Musacchio, S. & Boffetta, G. 2022 Instability of a dusty kolmogorov flow. Journal of Fluid Mechanics 931.
  * Squires & Eaton (1990) Squires, K. D. & Eaton, J. K. 1990 Particle response and turbulence modification in isotropic turbulence. Physics of Fluids A: Fluid Dynamics 2 (7), 1191–1203.
  * Squires & Eaton (1991) Squires, K. D. & Eaton, J. K. 1991 Preferential concentration of particles by turbulence. Physics of Fluids A: Fluid Dynamics 3 (5), 1169–1178.
  * Squires & Eaton (2014) Squires, K. D & Eaton, J. K 2014 Preferential concentration of particles by turbulence. Phys. Fluids A 3 (5), 11.
  * Thomson (1880) Thomson, W. 1880 XXIV. Vibrations of a columnar vortex. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 10 (61), 155–168.
  * Trefethen (2000) Trefethen, L. N. 2000 Spectral methods in MATLAB. SIAM.
  * Wang & Maxey (1993) Wang, L. & Maxey, M. R. 1993 Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence. Journal of Fluid Mechanics 256, 27–68.