# Clustering and chaotic motion of heavy inertial particles in an isolated
non-axisymmetric vortex

Anu V. S. Nath1  Anubhab Roy1<EMAIL_ADDRESS>1Department of Applied
Mechanics, Indian Institute of Technology Madras, Chennai 600036

###### Abstract

We investigate the dynamics of heavy inertial particles in a flow field due to
an isolated, non-axisymmetric vortex. For our study, we consider a canonical
elliptical vortex - the Kirchhoff vortex and its strained variant, the Kida
vortex. Contrary to the anticipated centrifugal dispersion of inertial
particles, which is typical in open vortical flows, we observe the clustering
of particles around co-rotating attractors near the Kirchhoff vortex due to
its non-axisymmetric nature. We analyze the inertia-modified stability
characteristics of the fixed points, highlighting how some of the fixed points
migrate in physical space, collide and then annihilate with increasing
particle inertia. The introduction of external straining, the Kida vortex
being an example, introduces chaotic tracer transport. Using a Melnikov
analysis, we show that particle inertia and external straining can compete,
where chaotic transport can be suppressed beyond a critical value of particle
inertia.

##  1 Introduction

Coherent vortical structures are ubiquitous in nature. In the planetary
context, the flows are predominantly two-dimensional due to the strong
influence of rotation. These turbulent flows naturally evolve into long-lived
isolated eddies/vortices (McWilliams, 1984), like in oceans and planetary
atmospheres — for example, the Great Red Spot on Jupiter, tropical cyclones,
and Gulf Stream rings. The inverse cascade of energy and planetary rotation
plays a prominent role in forming these coherent vortices (Vallis, 2017).The
two-dimensional (2D) turbulent flows in a rotating background will eventually
self-organize into concentrated vortical lumps as observed in various
simulations and observations (Fornberg, 1977; Basdevant et al., 1981; Babiano
et al., 1987; Benzi et al., 1988). In three-dimensional turbulent flows as
well, the formation of vortical filaments has been similarly identified
(Siggia, 1981; Vincent & Meneguzzi, 1991; Bartello et al., 1994). Thus, the
later evolution of turbulent flows is strongly influenced by the underlying
vorticity dynamics of these coherent eddies.

Turbulent flows in geophysical and astrophysical contexts are often dispersed
with particulate matter - water droplets and ice crystals in clouds (Shaw,
2003), pyroclastic flows (Dufek, 2016), wind-sand interactions in aeolian
processes (Kok et al., 2012), and dust in protoplanetary disks (Armitage,
2020). The coherent vortices embedded in these flows greatly influence
particle transport. Tracer particles can get trapped by the vortices for a
long time, much larger than the eddy turnover time, and get transported across
the distances over which the eddy travels (Elhmaïdi et al., 1993). The trapped
particles can be released only after the disruption of the vortex itself.
However, the particulate matter does not necessarily have negligible inertia;
the finite inertia aspect of particles can make the dynamics of the suspended
phase more complex with aspects of clustering (Bec, 2003, 2005; Sapsis &
Haller, 2010), and caustics (Crisanti et al., 1992; Falkovich et al., 2002;
Wilkinson & Mehlig, 2005). Heavy inertial particles are centrifuged away by
the vortex cores and get accumulated in the straining regions of the flow
(Maxey, 1987). However, in a rotating background, the heavy inertial particles
get pushed by the Coriolis force into the cores of anticyclonic vortices,
which is hypothesized to trigger the formation of planetesimals in the
astrophysical context (Tanga et al., 1996; Chavanis, 2000).

The transport of particles by various vortical structures has been extensively
modelled and studied in the past few decades. Batchelor & Nitsche (1994)
investigated the expulsion of heavy inertial particles from a rising bubble,
considering the combined effects of gravitational sedimentation and the
toroidal circulation of gas inside. Marcu et al. (1995) explored the transport
of inertial particles near a Burger’s vortex with and without the influence of
gravity. In the absence of gravity, particles with sufficiently small inertia
were captured by the vortex centre, while those with large inertia exhibited
stable limit cycle dynamics. The inclusion of gravity altered the dynamics,
generating additional fixed points that could capture the particles. Raju &
Meiburg (1997) have studied the transport of inertial particles with varying
density ratios in three model flows: a solid-body vortex, a point vortex, and
a stagnation point flow. Eames & Gilbertson (2004) considered the
sedimentation and dispersion of inertial particles past an isolated spherical
vortex and a random distribution of spherical vortices, revealing that the
interplay between particle inertia and stagnation points significantly
increased the vertical dispersivity of dense particles compared to tracers.
Hunt et al. (2007) investigated inertial particle transport near a vortex tube
and steadily propagating vortex rings, providing analytical treatment and
experimental comparisons. The clustering of heavy inertial particles in a pair
of co-rotating vortices was explored by Angilella (2010); Ravichandran et al.
(2014). Ravichandran & Govindarajan (2015) have studied the clustering of
inertial particles and the subsequent emergence of caustics in a point vortex
and a system of point vortices. In the context of airborne pathogen transport
through the atmosphere, recent studies by Dagan (2021) and Avni & Dagan (2022)
have modelled evaporating droplets as advected by a Lamb-Chaplygin vortex
dipole, revealing that the interaction with the vortex enhanced droplets’
settling time and transported them over large distances in the air. Motivated
by the dispersion of droplets in warm cumulus clouds during the condensational
phase, Nath et al. (2022) investigated the dispersion of condensing droplets
in a background flow modelled as an array of Taylor-Green vortices; they
demonstrated a significant enhancement in droplet dispersion as they acquire
more inertia by condensation.

Neighbouring coherent vortices in a turbulent flow can interact with each
other and induce shearing, which could disrupt them (Legras et al., 2001).
According to Reinaud et al. (2003), vortices that can withstand the highest
levels of strain are those most likely to be found in an actual turbulent
flow. An elliptic vortex patch of constant vorticity is an exact solution of
the incompressible 2D Euler equation (see Kirchhoff, 1876; Lamb, 1945;
Saffman, 1995); below aspect ratio of 3 the vortices are both linearly (Love,
1893) and nonlinearly (Wan, 1986; Tang, 1987) stable. Elliptic vortex patches
and their interactions have been extensively studied to understand better the
stability and evolution of vortices in ideal fluid (see Moore & Saffman, 1971;
Kida, 1981; Dritschel, 1990; Legras & Dritschel, 1991; Dritschel & JuáRez,
1996; Mitchell & Rossi, 2008), with motivations from geophysical turbulent
flows (Dritschel, 1995). Due to the non-axisymmetric vorticity distribution,
an elliptic patch of uniform vorticity ($\omega_{0}$) in an irrotational
background will rotate with constant angular velocity
$\Omega=\omega_{0}\,a\,b/(a+b)^{2}$, where $a$ and $b$ are semi-major and
semi-minor axes of the elliptic patch (see figure 1); the size and shape of
the elliptic patch is preserved during the rotation. The configuration, widely
known as the Kirchhoff vortex, is given by the stream function (as observed by
a stationary observer)

$\psi_{v}^{\prime}=\left\\{\begin{array}[]{ll}-\frac{\omega_{0}}{2\,(a+b)}\,(b\,x^{2}+a\,y^{2}),&\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}<1\vspace{0.35cm}\\\
-\frac{a\,b\,\omega_{0}}{4}\,\left(2\,\xi+e^{-2\,\xi}\,\cos(2\,\eta)\right),&\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}>1~{}.\\\
\end{array}\right.$ (1)

Figure 1: (a) Schematic showing the elliptic vortex patch with stationary
reference frame $X^{\prime}-Y^{\prime}$ and the co-rotating reference frame
$X-Y$. (b) The (steady) streamlines of the Kirchhoff vortex in the co-rotating
reference frame.

In the co-rotating frame with the vortex, the stream function is
$\psi=\psi^{\prime}+\frac{\Omega}{2}\,(x^{2}+y^{2})$, where $(x,y)$ and
$(\xi,\eta)$ are respectively the Cartesian and elliptic coordinates measured
in a co-rotating frame with the ellipse, which are inter-related as
$x=\sqrt{a^{2}-b^{2}}\,\cosh\xi\,\cos\eta$ and
$y=\sqrt{a^{2}-b^{2}}\,\sinh\xi\,\sin\eta$, with $\xi\geq 0$ and
$\eta\in[0,2\pi)$. We use primed (′) variables to denote quantities in the
stationary frame, whereas un-primed variables represent quantities in the co-
rotating frame. We follow the same convention throughout this paper unless
specified otherwise. The stationary reference frame ($X^{\prime}-Y^{\prime}$)
and the co-rotating reference frame ($X-Y$) are shown schematically in figure
1(a), making an instantaneous angle $\theta$, relates them as
$x^{\prime}=x\,\cos\theta-y\,\sin\theta$ and
$y^{\prime}=x\,\sin\theta+y\,\cos\theta$, where $\frac{d\theta}{dt}=\Omega$.
An advantage of choosing a co-rotating frame is that the velocity field is
steady in the co-rotating frame. The corresponding streamlines are shown in
figure 1(b). The tracer dynamics in the co-rotating frame is thus governed by
a 2D, time-independent dynamical system, which guarantees non-chaotic fluid
pathlines. Inside the ellipse, the flow field is a solid-body rotation, and
far away in the outer region, it resembles a decaying point vortex. The flow-
field is a tripole structure (see Viúdez, 2021; Xu & Krasny, 2023) with five
fixed points (where the flow velocity is zero) A, B, C, D and O as marked in
the figure. The origin O ($\overline{\xi}_{0}$, $\overline{\eta}_{0}$) is an
elliptic fixed point; the pair A and B ($\overline{\xi}_{1}^{\pm}$,
$\overline{\eta}_{1}^{\pm}$) located along the major axis line outside the
ellipse are hyperbolic type fixed points; the pair C and D
($\overline{\xi}_{2}^{\pm}$, $\overline{\eta}_{2}^{\pm}$) located along the
minor axis line outside the ellipse (inside the lobes) are elliptic type fixed
points (see Kawakami & Funakoshi, 1999). The hyperbolic fixed points are
interconnected by two pairs of heteroclinic orbits, denoted as $H_{1}^{\pm}$
and $H_{2}^{\pm}$. As we proceed further into our analysis, we will learn the
significance of these heteroclinic orbits and their possible perturbation,
which is crucial for the onset of chaos. The fixed points are ‘fixed points’
only for a co-rotating observer; however, for a stationary observer, they
resemble ‘Lagrange points’ of celestial mechanics (Fitzpatrick, 2012).

Figure 2: Schematic showing three major kinds of dynamics exhibited by a
sheared elliptical vortex: (a) initial vortex patch of anti-clockwise
vorticity content sheared by an external shear flow of opposite sense of
vorticity content, (b-d) rotation: full rotation and periodic straining due to
weak shear compared to vorticity (e-g) nutation: back and forth oscillation
and periodic straining due to medium shear rate and (h-j) elongation: the
dominant shear flow stretches the vortex patch indefinitely and destroys it.

To study vortex interactions Kida (1981) proposed a model of a vortex tube in
a uniform shear flow; the effects of the other vortices on a certain vortex
tube may be replaced, in the first approximation, by a linear flow. Moore &
Saffman (1971) had earlier studied steady elliptic vortex patches in uniform
shear; Kida (1981) generalized the solutions to include exact unsteady
elliptic vortices. However, we would like to mention that the first study of
an elliptic vortex in a specific linear flow, a simple shear flow, was carried
out by Chaplygin (see Chaplygin, 1899; Meleshko & Van Heijst, 1994). When an
external linear flow $\psi^{\prime}_{e}=\frac{s}{4}\,(-x^{\prime 2}+y^{\prime
2})-\frac{\gamma}{4}\,(x^{\prime 2}+y^{\prime 2})~{}$, where $s$ is the strain
rate and $\gamma$ is the vorticity, is superimposed on the Kirchhoff vortex
($\psi^{\prime}_{v}$) the configuration is known to have three kinds of impact
on the elliptic vortex patch: (i) rotation: full rotation of ellipse with
periodically changing aspect ratio and angular velocity, happens when the
external strain rate is weak, (ii) nutation: back and forth oscillatory
angular motion of elliptic patch, happens when the external vorticity
component is strong enough and in opposite sense of the vorticity of the
elliptic patch, and (iii) elongation: complete straining and disruption of
vortex patch due to strong external straining (see Kida, 1981; Dritschel,
1990); also, see the Schematic in figure 2. Though the area of the ellipse is
preserved, the unsteady rotation of the ellipse with changing aspect ratio
creates an unsteady flow-field around, even in the co-rotating frame (Kida,
1981), which is referred to as the ‘Kida vortex’ in this paper. The original
motivation for studying sheared vortical patches was to understand vortex
interactions better. However, when analyzed from the perspective of tracer
transport, sheared elliptic vortices exhibited chaotic Lagrangian trajectories
(Polvani & Wisdom, 1990; Dahleh, 1992). Even a minor imposed shear induces
periodic unsteadiness in the deformation of the Kida vortex, disrupting the
Hamiltonian integrability. The hyperbolic fixed points and heteroclinic
connections experience perturbations, making them susceptible to transverse
intersections and thus allow for the possibility of chaotic dynamics (Smale,
1967; Bertozzi, 1987). A comprehensive investigation was conducted into the
impact of unsteady perturbations on tracer transport in an otherwise
integrable system of a pair of oppositely signed point vortices by Rom-Kedar
et al. (1990). In the absence of perturbation, the vortex pair translates with
a constant velocity, resulting in a steady flow field in the co-moving frame
with the vortices. However, the introduction of an external periodic strain
field, even in the co-moving frame, makes the flow field unsteady, causing the
tangling heteroclinic orbits in the flow field and leading to the chaotic
transport of certain passive tracers. This phenomenon also results in fluid
entrainment by the vortex system, enhancing mixing. Notably, this study
represents one of the early applications of tools such as the Melnikov
analysis from dynamical systems to analyze and quantify chaos and mixing in a
fluid flow problem. In the current context of inertial particle transport in
the Kida vortex, we apply some of the techniques derived from their work.

The dispersed phase embedded in the coherent structures in the various
geophysical and astrophysical flows is rarely inertialess - dust, bubbles,
planktons. This raises the question of how particle inertia alters particle
dispersion in the neighbourhood of vortices. Here, we are interested in the
dynamics of small, heavy inertial particles, thus ignoring the additional
physics of added mass effects, the Basset history term, convective inertia,
and Faxen corrections. Studies have shown that particle inertia can suppress
the chaotic transport in vortical flows (Angilella, 2010; Angilella et al.,
2014). However, we have demonstrated recently (Nath et al., 2024) that
particle inertia can induce chaotic dynamics and lead to non-ergodic dynamics
due to the ‘scattering’ interaction of inertial particles with an ordered
array of stagnation points. Particle inertia modifies the fluid tracer fixed
points and their homoclinic/heteroclinic connections, which subsequently play
a prominent role in long-time particle transport. In this paper, we study the
transport of heavy inertial particles near a non-axisymmetric vortex patch -
first in the Kirchhoff vortex and then in its strained variant, the Kida
vortex. The configuration chosen is a simple scenario of an isolated vortex
where we can study the modification of the heteroclinic tangles by particle
inertia analytically and comment on the competing roles of background shear
and particle inertia to promote or suppress chaotic transport.

An earlier investigation on particle transport in a strained elliptical vortex
is documented in the work by Chavanis (2000). This study specifically focuses
on the trapping of dust by anticyclonic vortices in Keplerian proto-planetary
disks, proposing it as a mechanism for planet formation. In this context, the
vortices experience Keplerian shear, leading to the survival of only
anticyclonic vortices that achieve a steady elliptic configuration. The
strength of the vorticity is determined by Keplerian shear, employing a
solution provided by Moore & Saffman (1971). The analysis of inertial dust
particle transport is conducted in a frame co-rotating with the vortex. The
findings indicate that particles with small inertia approximately follow an
elliptic path, drifting inwards due to drag and the Coriolis force, ultimately
being captured by the vortex centre. On the other hand, particles with larger
inertia exhibit an epicyclic motion but eventually sink into the vortex.
Particles with substantial inertia may even escape the vortex. Despite some
apparent similarities with the second part of our work, which involves the
transport of inertial particles by a strained elliptic vortex, there are
notable differences. Our study considers an elliptic vortex model for a
coherent vortex in turbulence, allowing for any general value of the strain
rate it experiences. Consequently, the vortex is not in a steady state but
unsteady motion, leading to the intriguing particle dynamics discussed in this
paper.

The remainder of this paper is organized as follows. Section 2 provides an
overview of the dynamics of heavy inertial particles in the Kirchhoff vortex,
detailing their clustering in various fixed points and presenting a stability
analysis of these fixed points. In Section 3, we see the modification to these
dynamics in a strained Kirchhoff vortex due to an imposed weak planar
extensional flow, i.e., in a Kida vortex. We analyze the perturbative changes
from stable fixed points to stable limit cycles. Additionally, a Melnikov
analysis is employed on saddle points to demonstrate the existence of chaotic
dynamics for inertial particles with sufficiently small inertia. The
dispersion characteristics of heavy inertial particles in both Kirchhoff and
Kida vortex are discussed in Section 4. We conclude in Section 5.

##  2  Dynamics of heavy inertial particles in an elliptic vortex

The dynamics of heavy inertial point particles in a background flow can be
studied using the Maxey-Riley equation (see Maxey & Riley, 1983). In a
rotating reference frame with the ellipse (of angular velocity $\Omega$), the
modified form of the Maxey-Riley equation by accounting for the pseudo forces
reads (Tanga et al., 1996), in nondimensional form

$\dot{\textbf{v}}=\frac{\textbf{u}(\textbf{x})-\textbf{v}}{St}+\textbf{x}\,\Omega^{2}-2\,\Omega\,\hat{\textbf{e}}_{z}\times\textbf{v}~{},$
(2)

where v is the particle velocity, $\textbf{u}(\textbf{x})$ is the velocity of
the Kirchhoff vortex (steady, in the co-rotating frame) evaluated at the
particle location x, $\hat{\textbf{e}}_{z}$ is the unit vector along the $z$
axis perpendicular to the plane. We use the length scale $\sqrt{a\,b}$ and the
rotation time scale $\omega_{0}^{-1}$ to nondimensionalize the system. Another
relevant time scale in the problem is $\tau_{p}=2/9\rho_{p}a^{2}/\rho_{g}\nu$
\- the relaxation time scale for a particle of characteristic size $a$ and
density $\rho_{p}$ navigating in a carrier phase of density $\rho_{g}$ and
kinematic viscosity $\nu$ respectively. The Stokes number
($St=\tau_{p}\,\omega_{0}$) quantifies the relative magnitude of the two time
scales and thus provides a nondimensional measure of particle inertia. We
denote the nondimensional time derivative using overdot ($\,\,\dot{}$ ). The
nondimensional quantities are represented with the same notation as
dimensional quantities, as we deal only with nondimensional quantities from
here onwards (unless specified explicitly). The elliptic vortex is
characterized by the nondimensional aspect ratio $r=b/a$, and without loss of
generality, we only focus here on the cases of $0\leq r<1$.

The first term on the right-hand side of 2 is the Stokes drag, the second is
the centrifugal force, and the third is the Coriolis force. The Euler force
will not appear here since the rotation rate is uniform (i.e.,
$\dot{\Omega}=0$). However, we will see later, in the context of the Kida
vortex, the Euler force needs to be accounted for due to its non-uniform
rotation rate. It is assumed that the particles are point size and thus in a
Stokes flow limit ($Re=\Omega a^{2}/\nu\ll 1$). Also, we are dealing with
heavy particles (i.e., much denser than the background fluid, as is the case
for dust particles or water droplets in the air). Thus, as was mentioned
earlier, the physics of added mass force and the Basset history effect is
negligible. However, suppose one were to study the role of fluid inertia with
the motivation of the current study and explore the long-time dispersion
dynamics of particles in the vicinity of coherent structures. In that case,
the inclusion of convective inertia ($Re\neq 0$) is expected to play a more
prominent role than the Basset history effect (see Lovalenti & Brady, 1993;
Dorgan & Loth, 2007). However, for the sake of simplicity, we restrict
ourselves to the $Re\ll 1$ regime and ignore the physics of both unsteady and
convective inertia.

In our study, we initialize a circular patch of inertial particles (randomly
distributed) with zero initial velocity in the Kirchhoff vortex around the
ellipse (see figure 3(a)). We let the particles evolve and track them using
the dynamic equation (2) along with the kinematic equation
$\dot{\textbf{x}}=\textbf{v}$. We integrate the system using the ODE113
routine in Matlab with a relative error of $10^{-12}$, absolute error of
$10^{-12}$ and a maximum time step of $1/10$ th of the Stokes number. The
typical evolution of $St=0.5$ particles in the Kirchhoff vortex of aspect
ratio $r=0.5$ observed from the co-rotating frame is shown as the snapshots in
figure 3.

Figure 3: Snapshots showing the evolution of $10^{5}$ particles of $St=0.5$ in
a Kirchhoff vortex of $r=0.5$ at (a) $t=0$, (b) $t=50$, (c) $t=100$ and (d)
$t=150$ nondimensional time units, obtained from numerical simulation. The
particles are initialized with zero velocity and randomly distributed inside a
circle of nondimensional radius $1.76$, enclosing the ellipse.

As expected, the snapshots show that the inertial particles are getting
centrifuged away from the central ellipse. However, unlike in the case of an
axisymmetric vortex (like point vortex or Rankine vortex), here we see that
some of the particles are getting attracted towards a pair of fixed points
outside the ellipse, within each lobe denoted as C and D. Note that, the
elliptic fixed points of fluid tracers in Kirchhoff vortex have already been
denoted as C and D in figure 1(b). The same choice of notation here for
attracting fixed points will make sense in the analytical exploration of the
system in the upcoming section.

###  2.1  Analytical evaluation of fixed points and their stability

The fixed points, as seen by an inertial particle in the Kirchhoff vortex, can
be evaluated by setting its velocity and acceleration to be zero, i.e.
$\textbf{v}=\dot{\textbf{v}}=\textbf{0}$. Substituting this in equation (2)
gives that the locations of the fixed points are the solution of the equation
$\textbf{u}(\textbf{x})+\textbf{x}\,St\,\Omega^{2}=\textbf{0}$, i.e., the
fixed points are formed by the balance between centrifugal force and Stokes
drag. Note that the fixed point equation is modified from that of fluid
tracers ($\textbf{u}(\textbf{x})=0$) due to the finite inertia effect as a
$St$ dependent term. For analytical treatment to be made accessible, we may
choose to rewrite the equation (2) in elliptic coordinates, in component form,
as

$\displaystyle\ddot{\xi}$ $\displaystyle=$
$\displaystyle\frac{h\,k^{-1}\,u_{\xi}-\dot{\xi}}{St}+2\,\Omega\,\dot{\eta}+h^{2}\,\left\\{\frac{\dot{\eta}^{2}-\dot{\xi}^{2}+\Omega^{2}}{2}\sinh
2\xi-\dot{\eta}\,\dot{\xi}\,\sin 2\eta\right\\}~{},$ (3a)
$\displaystyle\ddot{\eta}$ $\displaystyle=$
$\displaystyle\frac{h\,k^{-1}\,u_{\eta}-\dot{\eta}}{St}-2\,\Omega\,\dot{\xi}-h^{2}\,\left\\{\frac{\dot{\eta}^{2}-\dot{\xi}^{2}+\Omega^{2}}{2}\sin
2\eta+\dot{\eta}\,\dot{\xi}\,\sinh 2\xi\right\\}~{},$ (3b)

where $x=k\,\cosh\xi\,\cos\eta$ and $y=k\,\sinh\xi\,\sin\eta$, and the fluid
velocity components can be obtained from the corresponding stream function in
co-rotating frame ($\psi$) as
$u_{\xi}=\frac{h}{k}\,\frac{\partial\psi}{\partial\eta}$ and
$u_{\eta}=-\frac{h}{k}\,\frac{\partial\psi}{\partial\xi}$. For the Kirchhoff
vortex (of stream function
$\psi_{v}=\psi_{v}^{\prime}+\frac{\Omega}{2}\,(x^{2}+y^{2})$), these
components can be evaluated as

$\displaystyle u_{\xi}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{h}{2\,k}\,\left(e^{-2\xi}-k^{2}\,\Omega\right)\,\sin
2\eta,&\quad\quad\quad\quad\quad\tanh\xi>r\vspace{0.25cm}\\\
\frac{h}{4}\,k^{3}\,\Omega\,(\Lambda-\cosh 2\xi)\,\sin
2\eta,&\quad\quad\quad\quad\quad\tanh\xi<r\\\ \end{array}\right.$ (4c)
$\displaystyle u_{\eta}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{h}{2\,k}\,\left(1-e^{-2\xi}\,\cos
2\eta-k^{2}\,\Omega\,\sinh 2\xi\right),&\quad\,\tanh\xi>r\vspace{0.25cm}\\\
\frac{h}{4}\,k^{3}\,\Omega\,(\Lambda-\cos 2\eta)\,\sinh
2\xi,&\quad\,\tanh\xi<r~{},\\\ \end{array}\right.$ (4f)

where the scale factor $h=(\cosh^{2}\xi-\cos^{2}\eta)^{-1/2}$, the parameters
$k^{2}=(1/r-r)$, $\Lambda=(1+r^{2})/(1-r^{2})$ and the angular velocity
$\Omega=\dot{\theta}=r/(r+1)^{2}$. Here, $\tanh\xi>r$ indicates the region
outside the ellipse and $\tanh\xi<r$ indicates the region inside the ellipse
since $\tanh\xi=r$ defines the boundary of the ellipse itself. The same
equations in Cartesian coordinates are convenient in numerical simulations and
can be found in Appendix A.

The system of equations (3), along with the appropriate velocity field,
describes the trajectory of an inertial particle in an elliptic vortex. It is
a nonlinear coupled dynamical system in a four-dimensional phase-space on the
variables $\xi,\eta,\dot{\xi}$ and $\dot{\eta}$. The fixed points of the
system can be obtained by solving the equations $\dot{\xi}=0$, $\dot{\eta}=0$,
$\ddot{\xi}=0$ and $\ddot{\eta}=0$ simultaneously. From equations (3), this
gives the trivial criteria for all fixed points, i.e.
$\dot{\overline{\xi}}=\dot{\overline{\eta}}=0$; however, their locations
$(\overline{\xi},\overline{\eta})$ should be obtained by solving the
transcendental equations

$\displaystyle 2\,u_{\xi}+St\,h\,k\,\Omega^{2}\,\sinh 2\xi=0~{},$ (5a)
$\displaystyle 2\,u_{\eta}-St\,h\,k\,\Omega^{2}\,\sin 2\eta=0~{}.$ (5b)

These equations must be solved separately inside and outside the ellipse to
identify the fixed points since the velocity field is known in a piecewise
manner. The identification of fixed points and their stability analysis are
discussed in the following Sections 2.1.1 and 2.1.2. Note that setting $St=0$
should recover the dynamics of the passive fluid tracers. By doing so, one
could see that the fixed point equations (5) reduce to $u_{\xi}=0$ and
$u_{\eta}=0$ \- which will retrieve the five classical fixed points of the
Kirchhoff vortex mentioned in the introduction.

####  2.1.1  Fixed points outside the ellipse

Outside the ellipse ($\tanh\xi>r$), the fixed point equations (5) can be
written by substituting the appropriate velocity expressions from equations
(4) as

$\displaystyle(St\,\Omega^{2}\,\sinh 2\xi+(k^{-2}\,e^{-2\xi}-\Omega)\,\sin
2\eta)\,h^{2}=0~{},$ (6a) $\displaystyle(St\,\Omega^{2}\,\sin
2\eta+\Omega\,\sinh 2\xi+k^{-2}\,(e^{-2\xi}\,\cos 2\eta-1))\,h^{2}=0~{}.$ (6b)

These are modified equations for fixed points outside the elliptic vortex
accounting for the effect of finite $St$. The solutions to equations (6) gives
the fixed points ($\overline{\xi},\overline{\eta}$) with
$p=\tanh\overline{\xi}$ and $q=\tan\overline{\eta}$ are given by

$\displaystyle p^{\pm}$ $\displaystyle=$
$\displaystyle\frac{2-(\alpha\pm\beta)\,k^{2}\,\Omega^{2}}{2\,(1+k^{2}\,\Omega)}~{},$
(7a) $\displaystyle q^{\pm}$ $\displaystyle=$
$\displaystyle\frac{\alpha\mp\beta}{2\,St}~{},$ (7b)

where $\alpha=k^{2}\,(1+St^{2}\,\Omega^{2})$ and
$\beta=\sqrt{\alpha^{2}-4\,St^{2}}$. These solutions form a set of four fixed
points outside the ellipse, located at ($\overline{\xi}_{1}^{+}$,
$\overline{\eta}_{1}^{+}$) =$(\tanh^{-1}p^{+},\tan^{-1}q^{+})$,
($\overline{\xi}_{1}^{-}$, $\overline{\eta}_{1}^{-}$)
=$(\tanh^{-1}p^{+},-\pi+\tan^{-1}q^{+})$, ($\overline{\xi}_{2}^{+}$,
$\overline{\eta}_{2}^{+}$) = $(\tanh^{-1}p^{-},\tan^{-1}q^{-})$ and
($\overline{\xi}_{2}^{-}$, $\overline{\eta}_{2}^{-}$)
=$(\tanh^{-1}p^{-},-\pi+\tan^{-1}q^{-})$. In the limit of $St\rightarrow 0$,
we may deduce that these fixed points coincide with the Kirchhoff vortex’s
four classical fixed points A, B, C and D. For simplicity, we chose to call
these finite $St$ modified fixed points with the same name as that
corresponding to the fluid tracers (A, B, C, and D). The finite $St$
symmetrically displaces these fixed points, as shown in figure 4(a): A and B
shift counter-clockwise, while C and D shift clockwise as $St$ increases. As
$St$ increases, the fixed points A($\overline{\xi}_{1}^{+}$,
$\overline{\eta}_{1}^{+}$) and C($\overline{\xi}_{2}^{+}$,
$\overline{\eta}_{2}^{+}$) approaches each other. The same thing happens for
the counterpart fixed points B($\overline{\xi}_{1}^{-}$,
$\overline{\eta}_{1}^{-}$) and D($\overline{\xi}_{2}^{-}$,
$\overline{\eta}_{2}^{-}$) as well.

In the limit of small Stokes number ($St\ll 1$), the governing equation (2)
can be reduced to the slow-manifold form as

$\textbf{v}=\textbf{u}-St\,\left\\{\frac{\partial\textbf{u}}{\partial
t}+\textbf{u}\cdot\boldsymbol{\nabla}\textbf{u}-\textbf{x}\,\Omega^{2}+2\,\Omega\,\hat{\textbf{e}}_{z}\times\textbf{u}\right\\}+\textit{O}(St^{2})~{}.$
(8)

using which the particle trajectories are obtained and shown in In figure
4(b). By comparing these particle trajectories with the streamlines for fluid
tracers shown in figure 1(b), it is evident that the inertial particles
perceive different fixed points compared to fluid tracers. They match exactly
with the solution of equation (6) we obtained analytically (as marked in blue
and red). The trajectories also imply that the fixed points A and B remain
hyperbolic fixed points (saddles); however, the fixed points C and D behave as
stable spirals for inertial particles. The dissipative nature of the system
(3) for finite inertia particles is responsible for the elliptic fixed points
becoming spiral attractors. However, the hyperbolic fixed points remain
intact, though their position got displaced.

To analytically show this, we use the linear stability analysis and
systematically study the effect of finite $St$ on the stability of the fixed
points. The Jacobian matrix for the system of differential equations (3) is
evaluated about the fixed points is

$\mathsfbi{J}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\
\Gamma&-\Delta&-\frac{1}{St}&2\,\Omega\\\
\Delta-\frac{2\,k^{-2}\,(1-p)\,(p+q^{2})}{St\,(p^{2}+q^{2})}&\Gamma+\frac{2\,k^{-2}\,q\,(1-p)^{2}}{St\,(p^{2}+q^{2})}&-2\,\Omega&-\frac{1}{St}\end{pmatrix}~{},$
(9)

Figure 4: (a) Variation of the fixed points perceived by inertial particles in
a Kirchhoff vortex of $r=0.5$ as $St$ changes (for log-spaced distribution of
$St$). Elliptic and hyperbolic fixed points outside the ellipse merge at
$St\approx 0.772$. (b) The trajectories of $St=0.2$ particles in a Kirchhoff
of $r=0.5$ in the co-rotating reference frame, obtained using slow manifold
expression. Red indicates unstable fixed points, and blue indicates stable
fixed points. Figure 5: Variation of all eigenvalues with $St$ corresponding
to the typical fixed points of inertial particles in a Kirchhoff vortex of
$r=0.5$ : (a,d) Real and imaginary parts of the eigenvalues of fixed points A
or B. (b,e) Real and imaginary parts of the eigenvalues of fixed point C or D.
(c,f) Real and imaginary parts of the eigenvalues of fixed point O.

where
$\Gamma=(p^{2}-1)\\{\Omega\,(1+q^{2})\,(St\,\Omega\,(p^{2}-q^{2})-2\,p\,q)+2\,k^{-2}\,(p+q^{2})\,q\,(1-p)\\}/\\{St\,(p^{2}+q^{2})^{2}\\}$
and
$\Delta=(1-p)\,(1+q^{2})\,\\{\Omega\,(1+p)\,(p^{2}-q^{2}+2\,p\,q\,St\,\Omega)-k^{-2}\,(1-p)\,(p^{2}-q^{2})\\}/\\{St\,(p^{2}+q^{2})^{2}\\}$.
The same matrix evaluated in Cartesian coordinates can be found in Appendix A.
By evaluating the eigenvalues for the matrix $\mathsfbi{J}$, we can identify
that the fixed point pairs A and B are saddles, and the pairs C and D are
stable spirals for a finite inertial particle. For instance, the variation of
all the eigenvalues with Stokes number for each fixed point type is shown in
figure 5. From figures 5(a) and (d), it can be seen that any finite $St$
particle will perceive fixed points A & B with purely real-valued eigenvalues
with signature $(+,-,-,-)$ indicating saddles in the four-dimensional phase
space. However, for the fixed point C/D, as shown in figures 5(b) and (e), the
eigenvalue will be a pair of complex conjugates with a negative real part,
indicating a stable spiral in the four-dimensional phase space. Thus, we
conclude that the suspended heavy inertial particles will spiral and cluster
towards the fixed points C and D outside the elliptic vortex as time
progresses, which we have observed in numerical simulations (see figure 3).
Note that, for the fixed point C/D, for $r=0.5$, when $St\gtrsim 0.771$ (shown
in the insets of figures 5(b) and (e)), we may see that there occurs a
bifurcation and the eigenvalues are no more conjugate pairs, rather purely
real-valued with signature $(-,-,-,-)$, indicating stable node/sink in the
four-dimensional phase space. Consequently, the particles will radially move
towards the fixed points instead of executing a spiral motion. For an elliptic
vortex of arbitrary aspect ratio $r$, the critical Stokes number
$\mathcal{S}_{1}$ at which the spiral to node transition happens can be
identified from the discriminant of the solutions of quartic eigenvalue
polynomial for $\mathsfbi{J}$ (which can be obtained from equation (9) or
equation (25), by setting $\lvert\mathsfbi{J}-\lambda\,\mathsfbi{I}\rvert=0$).
By examining the discriminant, one may find its behaviour changes when
$St^{2}$ satisfies the cubic equation

$\begin{split}St^{6}+\frac{3\,(1+r)^{4}\,(1+4\,r+r^{2})}{2\,r^{3}}\,St^{4}+\frac{(1+r)^{8}\,(1+6\,r+r^{2})\,(5+22\,r+5\,r^{2})}{16\,r^{6}}\,St^{2}\\\
-\frac{(1+r)^{12}\,(1-r)^{2}\,(3+r)\,(1+3\,r)}{16\,r^{8}}=0~{}.\end{split}$
(10)

The real-valued, non-negative solution of the equation (10) gives the critical
Stokes number $St=\mathcal{S}_{1}$ for any $r$. For the case of $r=0.5$, one
could verify that this solution is $\mathcal{S}_{1}\approx 0.7713$, matching
with the prediction from figures 5(b) and (e).

As $St$ increases, saddle type and stable spiral type fixed points approach
closer, as mentioned earlier. If we keep increasing $St$, we will find that
they merge (i.e., $p^{+}=p^{-}$ and $q^{+}=q^{-}$) and vanish at a critical
inertia value. From equation (7), one could deduce that this happens only if
$\beta=0$. Using the expression for $\beta$ and solving for Stokes number, one
would find the critical value at which merging happens is
$St=\mathcal{S}_{2}=(1+r)^{2}\,(1-\sqrt{r})/(r\,(1+\sqrt{r}))$. i.e., the
fixed points outside the elliptic vortex (A, B, C and D) exist only for
particles with $0\leq St<\mathcal{S}_{2}$. For an elliptic vortex with
$r=0.5$, the critical value can be evaluated as $\mathcal{S}_{2}\approx
0.7721$, i.e. A merges with C and B merges with D at this critical value and
disappears. Thus, we have only shown the eigenvalues in the figure 5 for
Stokes number until $0.7721$.

####  2.1.2  Fixed points inside the ellipse

The origin O located inside the ellipse ($\tanh\xi<r$) continues to exist for
inertial particles (see figure 4(b)). However, for inertial particles, O
behaves as an unstable spiral. To show this, we follow the same procedure
mentioned in the previous Section 2.1.1, using the flow field inside the
ellipse. By substituting appropriate expressions for $u_{\xi}$ and $u_{\eta}$
from equations (4) in the fixed point equations (5), we obtain

$\displaystyle(2\,St\,\Omega\,\sinh 2\xi+k^{2}\,(\Lambda-\cosh 2\xi)\,\sin
2\eta)\,h^{2}=0~{},$ (11a) $\displaystyle(2\,St\,\Omega\,\sin
2\eta-k^{2}\,(\Lambda-\cos 2\eta)\,\sinh 2\xi)\,h^{2}=0~{}.$ (11b)

In the $St\rightarrow 0$ limit, these equations yield the elliptic fixed point
for fluid tracers at the origin. Moreover, irrespective of the value of $St$,
equations (11) has a single real-valued solution
$(\overline{\xi}_{0},\overline{\eta}_{0})=(0,\pi/2)$, indicating that the
origin O remains to be the fixed point for any inertial particle as well. The
Jacobian matrix for linear stability for the fixed point at origin O is,

$\mathsfbi{J}=\frac{1}{St}\,\begin{pmatrix}0&0&St&0\\\ 0&0&0&St\\\
St\,\Omega^{2}&-r\,\Omega&-1&2\,St\,\Omega\\\
\Omega/r&St\,\Omega^{2}&-2\,St\,\Omega&-1\end{pmatrix}~{},$ (12)

which has two pairs of complex conjugate eigenvalues. One of these pairs has a
non-negative real part (see figures 5(c) and (f)), indicating unstable spiral
behaviour for any nonzero $St$. Thus, the particles will be centrifuged away
from the origin spirally. Some of them, starting within certain basins (see
Section 2.3) accumulate in the stable fixed points C and D outside the
ellipse, as they are stable fixed points, as shown in the previous Section
2.1.1. Note that as we increase the particle inertia, even after the mutual
annihilation of fixed points A, B, C, and D, fixed point O exists and remains
an unstable spiral. Thus, beyond the $\mathcal{S}_{2}$, one can observe that
all particles merely centrifuge away without getting trapped in any Lagrange
points.

Figure 6: A schematic showing the projection of a four-dimensional phase space
topology into three dimensions, illustrating various fixed points and their
transitions at critical Stokes numbers.

The nature of fixed points and phase space topology is depicted schematically
in figure 6 as the Stokes number varies. It is important to note that the
phase space is four-dimensional; however, the schematic only presents a three-
dimensional projection. The fixed points C/D exhibit a ‘2-spiral sink’
behavior when $St<\mathcal{S}_{1}$ and transition to a ‘sink’ when
$\mathcal{S}_{1}<St<\mathcal{S}_{2}$. On the other hand, fixed points A/B
demonstrate a ‘3:1 saddle’ behaviour for $St<\mathcal{S}_{2}$. Additionally,
fixed points A (B) annihilate with C (D) and cease to exist when
$St=\mathcal{S}_{2}$. Meanwhile, fixed point O persists and behaves as a
‘2-spiral saddle’ for all Stokes numbers. For more on terminology, see Hofmann
et al. (2018).

Figure 7: The curves in the $r-St$ plane demarcate regions where the change in
the number of fixed points (blue colour) and the change in the nature of
stable fixed points (red colour) happens. The representative elliptic patch of
vorticity is also shown in greenish yellow color for three different aspect
ratios.

###  2.2  Effect of the aspect ratio on the distribution of fixed points

Until now, all the analyses have been performed explicitly for the elliptic
vortex with aspect ratio $r=0.5$. However, if the vortex has a different
aspect ratio, the distribution and stability nature of fixed points will
change. Solving equations (6) and 11), we obtain the fixed points inside and
outside the ellipse, respectively, for various $St$ and $r$ values. There is
no qualitative change in the distribution of fixed points, such that the fixed
point inside the ellipse will always be at the origin, and there will always
be four fixed points outside the ellipse below $\mathcal{S}_{2}$. In the
$r-St$ plane (see figure 7), it is shown that the critical curve (blue colour)
which demarcates the region where all the five fixed points (A, B, C, D and O)
coexists and the region where only the one at the origin only exists.

By analysing the stability of the fixed points, we have already seen that the
fixed point at the origin remains an unstable spiral for any $St>0$ in a
Kirchhoff vortex of $r=0.5$. One can verify that this will also be valid for
any $r\in(0,1)$. The pair A and B remain saddles for all aspect ratios
provided they exist (indicating the robustness of hyperbolic fixed points).
However, the stable spiral fixed points can change their behaviour if the
aspect ratio is above some critical value. The same fact has already been
discussed towards the end of Section 2.1.1 specifically for the aspect ratio
$r=0.5$. For any particular aspect ratio, some critical Stokes number
$\mathcal{S}_{1}$ exists above which the stable spirals become stable
nodes/sinks. As mentioned in Section 2.1.1, by solving for the non-negative
real-valued solution of equation (10), one can find the critical pairs of $St$
and $r$ at which this behaviour change happens, and is shown in the $r-St$
plane (see figure 7, red colour). From asymptotic analysis of equation (10),
one may show that
$\mathcal{S}_{1}=r^{-1}\,\sqrt{3/5}-298/(25\,\sqrt{15})+\textit{O}(r)$ for
$r\rightarrow 0$ and
$\mathcal{S}_{1}=(1-r)+\frac{1}{2}\,(1-r)^{2}+\textit{O}(1-r)^{3}$ for
$r\rightarrow 1$, and the asymptotes are also shown in the figure. On the
other hand the exact expression for
$\mathcal{S}_{2}=(1+r)^{2}\,(1-\sqrt{r})/(r\,(1+\sqrt{r}))$, mentioned at the
end of Section 2.1.1 has asymptotic forms,
$\mathcal{S}_{2}=r^{-1}-2/\sqrt{r}+4+\textit{O}(\sqrt{r})$ for $r\rightarrow
0$ and $\mathcal{S}_{2}=(1-r)+\frac{1}{2}\,(1-r)^{2}+\textit{O}((1-r)^{3})$
for $r\rightarrow 1$. We find that $\mathcal{S}_{1}<\mathcal{S}_{2}$ for any
$r\in(0,1)$. Note that the red curve is very close to the blue curve; however,
they never intersect, indicating that for any fixed $St>0$, there always
exists a narrow parameter regime bounded by both red and blue curves
($\mathcal{S}_{1}<St<\mathcal{S}_{2}$) inside which the fixed points C and D
will behave as nodes/sinks.

For any finite $r$ value, one may note that both $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ values are finite. However, as $r\rightarrow 0$, both
diverge as $\textit{O}(1/r)$, becoming a larger value. This divergence
suggests that particles must possess high inertia for critical dynamical
behavioural changes to occur. However, note that when $r\rightarrow 0$, the
vortex rotates slowly ($\Omega=r/(1+r)^{2}\rightarrow 0$). Thus, the relevant
timescale in the problem becomes the angular velocity $\Omega$ rather than the
vorticity of the patch. Consequently, the Stokes number based on the angular
velocity of the vortex, $St_{\Omega}=St\,r/(r+1)^{2}$, evaluated for critical
behaviours $\mathcal{S}_{\Omega,1}\sim\sqrt{3/5}$ and
$\mathcal{S}_{\Omega,2}\sim 1$, remains constant and bounded in the limit of
$r\rightarrow 0$, indicating that the divergence was merely an artefact of the
adopted scaling.

###  2.3  The basin of attraction of the fixed points

From numerical simulations, we have seen that all the particles repelled away
from the unstable fixed points (A, B and O) are not attracted to the stable
fixed points (C and D); instead, some spiral away to infinity. The initial
particle locations that would result in them getting attracted to any of the
two stable fixed points are shown in figure 8, representing the basin of
attraction of those fixed points. The basins are coloured to distinguish each
other. The region outside these coloured regions indicates the basin of
attraction of infinity; i.e., particles initialized in these regions will
eventually spiral away to infinity and never get attracted to any of the
stable fixed points. We have used $10^{5}$ particles randomly initialized in
the flow field to generate the figure. We tracked their evolution numerically
and identified those that would approach any fixed points over a long time. We
marked their initial locations using specific colours to distinguish between
the basins of each fixed point. As visible in figure 8, the basin of
attraction of the stable fixed points shrinks as $St$ increases and vanishes.
It can be verified that the basins will disappear beyond the critical value
$\mathcal{S}_{2}$, indicating the absence of stable fixed points. The figure
shows that the stable fixed points C and D are enclosed within the
corresponding basin of attraction, which is obvious. However, the unstable
fixed points A and B appear to fall at the edge of the basins. The regular
nature of the basins shows non-chaotic dynamics of inertial particles in the
elliptic vortex. Close to origin O, the basin of attraction of stable fixed
points (coloured regions), and that of infinity (empty region) forms an inter-
twisting pattern. As a result, a particle starting closer to the origin will
have dynamics sensitive to the initial condition. A small change in the
initial position could lead the particle to end up either in the fixed points
C or D or spiral away to infinity.

Here, we would like to highlight prior studies on the dynamics of inertial
particles in the neighbourhood of a pair of like-signed point vortices by
Angilella (2010); Ravichandran et al. (2014); Zhao et al. (2024). The system
of like-signed point vortices bears similarities with that of the Kirchhoff
vortex. The inertial particles exhibit qualitatively identical behaviour near
the point vortex pair. The types of fixed points, stability characteristics,
and the basin of attractions show qualitative similarities to those observed
in the case of the elliptic vortex. Though it seems like an elliptic vortex
can be simply replaced by a pair of like-signed vortices, it is worthwhile to
focus on the subtle differences. The flow field in both cases matches well in
the far field only. The Kirchhoff vortex and the pair vortices have different
features in the near-field. Notably, the origin behaves as a saddle in the
case of pair vortices, whereas in the elliptic vortex, it acts as a
center/spiral.

Figure 8: The basin of attraction for inertial particles of (a) $St=0.2$, (b)
$St=0.5$ and (c) $St=0.77$ in a Kirchhoff vortex of $r=0.5$. The corresponding
fixed points A, B, C and D are marked in each figure using the ‘$\times$’
symbol. The streamlines of the Kirchhoff vortex are shown in grey in the
background. Figure 9: Snapshots showing the evolution of $2\times 10^{5}$
inertial particles of $St=0.1$ in a Kida vortex of $r_{0}=0.5$, $\theta_{0}=0$
for (a-d) $s=0.01$ and (e-h) $s=0.035$, observed in a co-rotating frame,
obtained from numerical simulation. The time stamps are: (a,e) $t=50$, (b,f)
$t=100$, (c,g) $t=500$, (d,h) $t=2000$.The particles are initialised with zero
velocity inside a circular region of radius $2.5$, enclosing the ellipse. The
instantaneous streamlines of the corresponding Kida vortex are shown in grey
in the background. The insets of (d) and (h) are the enlarged views at
$t=2000$ showing the accumulation of particles in limit cycle trajectories,
obtained numerically, shown in red. The blue colour indicates the limit cycles
evaluated using a perturbative approach as in Section 3.1.

##  3  Dynamics of heavy inertial particles in a strained elliptic vortex

Here, we consider a special case of the Kida vortex, a Kirchhoff vortex
superimposed with weak pure straining/planar extensional flow ($\gamma=0$ and
$s\ll 1$) so that the ellipse will do a full rotation with changing aspect
ratio and angular velocity as governed by the equations

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle-s\,r\,\sin 2\theta~{},$
(13a) $\displaystyle\dot{\theta}$ $\displaystyle=$
$\displaystyle\Omega+\frac{s}{2}\,\Lambda\,\cos 2\theta~{},$ (13b)

where the external strain rate ($s$) and vorticity ($\gamma$) are
nondimensionalised with vorticity of the core $\omega_{0}$. Note that the form
of $\Omega=r/(r+1)^{2}$, $\Lambda=(1+r^{2})/(1-r^{2})$ and $k^{2}=(1-r^{2})/r$
remains the same as defined for Kirchhoff vortex; however, for Kida vortex,
they are defined using the instantaneous aspect ratio ($r=r(t)$), which is
evolving with time. The inertial particles in the co-rotating frame with the
ellipse are tracked using the modified form of the Maxey-Riley equation

$\dot{\textbf{v}}=\frac{\textbf{u}(\textbf{x},t)-\textbf{v}}{St}+\textbf{x}\,\Omega^{2}-2\,\Omega\,\hat{\textbf{e}}_{z}\times\textbf{v}-\dot{\Omega}\,\hat{\textbf{e}}_{z}\times\textbf{x}~{}.$
(14)

Note that the unsteady nature of the Kida vortex, even in the co-rotating
frame, is accounted for in the Stokes drag term, and the non-uniform rate of
rotation of the frame responsible for the Euler force – the last term on the
right-hand side. We simulate the dynamics of $St=0.1$ point particles in a
Kida vortex of initial aspect ratio $r(t=0)=r_{0}=0$ and orientation
$\theta(t=0)=\theta_{0}=0$. The simulation is performed for two different
strain rate values, and the particles’ evolution is shown in figure 9 as
snapshots, as observed in a co-rotating frame. For a smaller value of strain
rate ($s=0.01$), the overall dynamics resemble that in the case of the
Kirchhoff vortex, as can be seen from figure 9(a-d); most of the particles are
getting centrifuged away, and a fraction of them attracted towards a pair of
stable attractors. Unlike the Kirchhoff vortex, here the attractors are not
fixed points but limit cycles in extended phase space (including the time
axis), which can be seen as the red curve in the inset of figure 9(d),
obtained by tracking the particles for a long time.

For a larger value of the strain rate ($s=0.035$), the particle dynamics
become intricate, as evident in corresponding figures 9(e-h). We notice that
some particles are drawn towards an additional, larger-sized limit cycle,
depicted in red in the inset of figure 9(h). Contrary to a simple spiralling
dispersion, the remaining particles become trapped in stationary attractors
(for a co-rotating observer, which appears to be counter-rotating as in the
snapshots). Further details regarding this attractor and particle dispersion
can be found in Section 4. Additionally, the particle dynamics exhibit chaotic
behaviour, consistent with the effect of external strain on fluid tracers
(Polvani & Wisdom, 1990; Kawakami & Funakoshi, 1999). In this context, a
systematic analytical study is presented in the upcoming Sections 3.1 and 3.2
to illustrate that the interplay of inertia ($St$) and straining ($s$) leads
to complex dynamics, as depicted in figure 9.

###  3.1  Modification to the fixed points of the Kirchhoff vortex due to
external straining

For the analytical treatment to be accessible, the governing equation (14) is
written in component form using elliptic coordinates as

$\displaystyle\ddot{\xi}$ $\displaystyle=$
$\displaystyle\frac{h\,k^{-1}\,u_{\xi}-\dot{\xi}}{St}-\frac{2\,\dot{k}}{k}\,\dot{\xi}+2\,\dot{\theta}\,\dot{\eta}-\frac{h^{2}\,\dot{k}}{2\,k\,St}\,\sinh
2\xi$ (15a) $\displaystyle+$ $\displaystyle
h^{2}\,\left\\{-\left(\dot{\eta}\,\dot{\xi}+\frac{\dot{k}\,\dot{\theta}}{k}+\frac{\ddot{\theta}}{2}\right)\,\sin
2\eta+\left(\frac{\dot{\eta}^{2}-\dot{\xi}^{2}+\dot{\theta}^{2}}{2}-\frac{\ddot{k}}{2\,k}\right)\,\sinh
2\xi\right\\}~{},$ $\displaystyle\ddot{\eta}$ $\displaystyle=$
$\displaystyle\frac{h\,k^{-1}\,u_{\eta}-\dot{\eta}}{St}-\frac{2\,\dot{k}}{k}\,\dot{\eta}-2\,\dot{\theta}\,\dot{\xi}+\frac{h^{2}\,\dot{k}}{2\,k\,St}\,\sin
2\eta$ (15b) $\displaystyle-$ $\displaystyle
h^{2}\,\left\\{\left(\dot{\eta}\,\dot{\xi}+\frac{\dot{k}\,\dot{\theta}}{k}+\frac{\ddot{\theta}}{2}\right)\,\sinh
2\xi+\left(\frac{\dot{\eta}^{2}-\dot{\xi}^{2}+\dot{\theta}^{2}}{2}-\frac{\ddot{k}}{2\,k}\right)\,\sin
2\eta\right\\}~{}.$

Note that when $\dot{k}=\ddot{\theta}=0$, these equations will reduce to the
case of Kirchhoff vortex as in equations (3); however, which is not the case
for any finite strain $s\neq 0$. The flow velocity fields can be obtained from
the total stream function
($\psi=\psi_{v}^{\prime}+\psi_{e}^{\prime}+\frac{\dot{\theta}}{2}\,(x^{2}+y^{2})$)
as $u_{\xi}=\frac{h}{k}\,\frac{\partial\psi}{\partial\eta}$ and
$u_{\eta}=-\frac{h}{k}\,\frac{\partial\psi}{\partial\xi}$. One must remember
that the ellipse’s evolution equations (13) must also be integrated along to
obtain the particle trajectory. Since we know that the attractors of the
system are not fixed points but limit cycles from the simulations, we may
search for them using the approach by IJzermans & Hagmeijer (2006) as a
periodic modification to the fixed points of the unperturbed system. Let
($\overline{\xi},\overline{\eta}$) denote the fixed points for any inertial
particle in Kirchhoff vortex, a time-dependent perturbation of the form
$\xi=\overline{\xi}+s\,\xi^{\prime}(t)+\textit{O}(s^{2})$ and
$\eta=\overline{\eta}+s\,\eta^{\prime}(t)+\textit{O}(s^{2})$ may denote the
limit cycles for any small straining $s\ll 1$.

The system behaves as a perturbed Kirchhoff vortex in an external straining
flow of $s\ll 1$. A time-periodic, strain-dependent term will perturb the
Kirchhoff vortex’s Hamiltonian (stream function). For low strain rates, a
regular perturbation approach yields the solution of equations (13) as
$r=r_{0}+s\,r_{1}+\textit{O}(s^{2})$ and
$\theta=\theta_{0}+\Omega_{0}\,t+s\,\theta_{1}+\textit{O}(s^{2})$ where

$\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle-\frac{r_{0}}{\Omega_{0}}\,\sin(\Omega_{0}\,t)\,\sin(\Omega_{0}\,t+2\,\theta_{0})~{},$
(16a) $\displaystyle\theta_{1}$ $\displaystyle=$
$\displaystyle\frac{(1+r_{0}^{3})}{r_{0}\,(1-r_{0})}\,\sin(\Omega_{0}\,t)\,\cos(\Omega_{0}\,t+2\,\theta_{0})-\frac{t}{2}\,\frac{(1-r_{0})}{(1+r_{0})}\,\cos
2\theta_{0}~{},$ (16b)

Figure 10: The variation of limit cycles corresponds to the fixed points of
inertial particles of various $St$ in a Kida vortex of $r_{0}=0.5$,
$\theta_{0}=0$ and $s=0.035$, obtained using the perturbative approach (c.f.
figure 4(a)). The stable limit cycles are represented using continuous lines,
and the unstable limit cycles are represented using dashed lines. Only the
first quadrant is shown here. A similar pair of limit cycles also exists in
the third quadrant (not shown here). The perturbative approach does not
capture the larger limit cycles and thus is not shown here. The grey curves in
the background indicate the streamlines of the Kida vortex, and the black
curve indicates the ellipse at $t=0$.

with $r_{0}=r(t=0)$, $\theta_{0}=\theta(t=0)$ and
$\Omega_{0}=r_{0}/(r_{0}+1)^{2}$. Note that the leading order solutions (i.e.,
solutions when $s=0$) match the Kirchhoff vortex, indicating a constant aspect
ratio and uniform rotation rate. Since the parameter $k$ is dependent on
aspect ratio $r$, $k$ can also be expanded as a perturbation in $s$ as
$k=k_{0}-s\,k_{0}\,\Lambda_{0}\,r_{1}/(2\,r_{0})+\textit{O}(s^{2})$, where
$k_{0}^{2}=(1/r_{0}-r_{0})$ and $\Lambda_{0}=(1+r_{0}^{2})/(1-r_{0}^{2})$.
Similarly, the velocity field $(u_{\xi},u_{\eta})$ can also represent a strain
perturbation to that of the Kirchhoff vortex.

By using these perturbation expansions along with the perturbation ansatz for
the limit cycle co-ordinates
($\xi=\overline{\xi}+s\,\xi^{\prime}(t)+\textit{O}(s^{2})$ and
$\eta=\overline{\eta}+s\,\eta^{\prime}(t)+\textit{O}(s^{2})$), substituted in
the governing equations (15) will reduce to the matrix equation

$St\,\ddot{\chi}+\mathsfbi{L}\,\dot{\chi}+\mathsfbi{K}\,\chi=-\mathsfbi{M}\,\lambda-N$
(17)

where $\chi=[\xi^{\prime}(t),\eta^{\prime}(t)]^{\textrm{T}}$,
$\lambda=[\cos(2\,\Omega_{0}\,t+2\,\theta_{0}),\sin(2\,\Omega_{0}\,t+2\,\theta_{0})]^{\textrm{T}}$,
$N=[N_{1},N_{2}]^{\textrm{T}}$ are vectors, and $\mathsfbi{K}$, $\mathsfbi{L}$
and $\mathsfbi{M}$ are $2\times 2$ matrices. The entries $N$, $\mathsfbi{K}$,
$\mathsfbi{L}$ and $\mathsfbi{M}$ are determined by the Kirchhoff vortex and
its fixed points ($\overline{\xi},\overline{\eta}$) and are listed in the
Appendix B. The equation (17) represents a pair of coupled, forced, second-
order, linear ordinary differential equations. We are interested in the large-
time/steady-state solutions as they would tell us about the limit cycles. For
this, we assume the periodic solution of the form
$\chi=F_{0}+\mathsfbi{F}_{1}\,\lambda$ and a substitution yields two
independent matrix equations (which are related to the particular solution of
equation (17))

$\displaystyle\mathsfbi{K}\,F_{0}$ $\displaystyle=$ $\displaystyle-N~{},$
(18a)
$\displaystyle\left(4\,\Omega_{0}^{2}\,St\,\mathsfbi{I}-2\,\Omega_{0}\,\mathsfbi{L}\,\boldsymbol{\Sigma}-\mathsfbi{K}\right)\,\mathsfbi{F}_{1}$
$\displaystyle=$ $\displaystyle\mathsfbi{M}~{},$ (18b)

where $\mathsfbi{I}=\big{(}\begin{smallmatrix}1&0\\\
0&1\end{smallmatrix}\big{)}$ is $2\times 2$ identity matrix and
$\boldsymbol{\Sigma}=\big{(}\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\big{)}$. These equations are linear and a simple
inversion yields the unknowns $F_{0}$ and $\mathsfbi{F}_{1}$, thus the limit
cycles. However, the stability of the limit cycle is determined by the
transient dynamics of the solutions of equation (17). It is determined by the
eigenvalues of the matrix $\mathsfbi{P}$ which satisfy the quadratic matrix
equation
$St\,\mathsfbi{P}^{2}+\mathsfbi{L}\cdot\mathsfbi{P}+\mathsfbi{K}=\mathsfbi{0}$
(related to the homogeneous solution of equation (17)). Any eigenvalue with a
positive real part indicates an unstable limit cycle; if all the eigenvalues
have a negative real part, it is stable. The limit cycles as a perturbation to
the fixed points of Kirchhoff vortex
($[\xi,\eta]=[\overline{\xi},\overline{\eta}]^{T}+s\,\chi+\textit{O}(s^{2})$)
are then evaluated for the case of $r_{0}=0.5$ and $\theta_{0}=0$ and is shown
in figure 10 as $St$ varies. Note that the originally stable fixed points (C
and D) transform into stable limit cycles, as confirmed through numerical
verification by tracking particles over an extended period, as depicted in the
insets of Figure 9(d) and (h). While theoretically predicting that the
hyperbolic fixed points (A and B) transition into unstable limit cycles under
shear perturbation, it may be more accurate to state that the hyperbolic fixed
points retain their hyperbolic nature. However, the stable and unstable
manifolds connecting them are now susceptible to the oscillatory effects
induced by this unstable limit cycle, potentially leading to their
intersections and the formation of tangles. A detailed analysis of the
perturbations to the hyperbolic fixed points is presented in the upcoming
Section 3.2.

The perturbation to the fixed point at the origin O yields all the entries of
$\mathsfbi{M}$ and $N$ to be zero, indicating that no limit cycle perturbation
exists (either stable or unstable) to the origin. The origin remains an
unstable spiral for inertial particles in the Kida vortex, which centrifuges
them away from the centre.

In the inset of figure 9(d) and (h), the perturbatively obtained limit cycle
solutions are shown in blue and agree with the respective actual limit cycles
(red colour, obtained numerically). However, for the case of $s=0.035$, an
additional pair of larger limit cycles exists, as can be seen in the inset of
figure 9(h) that this method could not capture.

###  3.2  Non-integrable perturbation to saddles and criteria for chaos

The external straining ($s\ll 1$) introduces a time-periodic perturbation
(non-integrable) to the Kirchhoff vortex (integrable Hamiltonian system). This
perturbation leads to Lagrangian chaos of fluid tracers ($St=0$) in the Kida
vortex (Polvani & Wisdom, 1990). The perturbations affect the hyperbolic fixed
points of the Kirchhoff vortex, leading to the heteroclinic orbits to split
into stable and unstable manifolds, which transversely intersect and form
tangles, as shown by the zero crossings of the corresponding Melnikov
functions. The result is a topological similarity to a horseshoe map and a
dense set near the intersections (see Smale, 1967; Ott, 2002), which
guarantees the existence of a strange attractor and chaotic dynamics. Though
the original theorem is formulated for homoclinic orbits, a later extension
for heteroclinic orbits conveys a similar concept (see Bertozzi, 1987). Here,
we study the modification to the dynamics of suspended heavy inertial
particles in the Kida vortex due to their finite inertia using the same
approach. A competing effect may arise for inertial particles due to their
finite inertia and suppress chaos as seen from similar studies (Angilella,
2010; Angilella et al., 2014). To analytically study this, we consider
modifying the particle dynamics in the Kida vortex (of $s\ll 1$) due to weak
inertial particles ($St\ll 1$).

Using the slow manifold expansion of equations (15) for $St\ll 1$, we obtain
the modified kinematic equations for $\dot{\xi}$ and $\dot{\eta}$ as
perturbations in $St$. By combining the $s\ll 1$ expansion of instantaneous
aspect ratio $r$ and angle of rotation $\theta$ with this slow manifold
equation, we can express the modified kinematic equations describing the
trajectory of an inertial particle in Kida vortex as

$\displaystyle\dot{\xi}=\hat{f}_{1}(\xi,\eta)+s\,\hat{g}_{1}(\xi,\eta;t)+St\,\hat{\phi}_{1}(\xi,\eta)+\textit{O}(s^{2},St^{2},s\,St)~{},$
(19a)
$\displaystyle\dot{\eta}=\hat{f}_{2}(\xi,\eta)+s\,\hat{g}_{2}(\xi,\eta;t)+St\,\hat{\phi}_{2}(\xi,\eta)+\textit{O}(s^{2},St^{2},s\,St)~{},$
(19b)

where the relevant functions $\hat{f}_{1}$ to $\hat{\phi}_{2}$ outside and
inside the ellipse are listed in Appendix C. In equations (19), the functions
$\hat{f}_{1}$, and $\hat{f}_{2}$ describe the dynamics resulting from the
dominant Kirchhoff vortex (integrable); $\hat{g}_{1}$ and $\hat{g}_{2}$ are
time-periodic perturbations due to external straining flow; $\hat{\phi}_{1}$
and $\hat{\phi}_{2}$ are the time-independent perturbations due to particle
inertia which is novel in our work. Thus, the dynamics of weak inertial
particles in a Kida vortex will be perturbations to that of fluid tracers in a
Kirchhoff vortex, where the perturbations are together contributed from
external straining and particle inertia.

The dynamics of fluid tracers in a Kirchhoff vortex (i.e., equations (19) with
$s=0$ and $St=0$) is integrable and thus represents a Hamiltonian system. The
introduction of an external straining will not affect the integrability of the
system inside the elliptic vortex as the flow-field inside is pure solid body
rotation and devoid of any hyperbolic fixed points. However, the flow-field
outside the ellipse has two hyperbolic fixed points, A and B, connected by the
two pairs of heteroclinic orbits $\textrm{H}_{1}^{\pm}$ and
$\textrm{H}_{2}^{\pm}$ as shown in figure 1(b). The strain perturbation can
lead to the formation of heteroclinic tangles. To identify that, we use
Melnikov’s method. According to the method, one must evaluate the Melnikov
function for all hyperbolic fixed points, which measures the distance between
stable and unstable manifolds in the Poincaré section. Odd zeros of the
Melnikov function indicate a transverse crossing of the corresponding stable
and unstable manifolds (see Bertozzi, 1987; Wiggins, 1990). For fluid tracers
in the Kida vortex, Kawakami & Funakoshi (1999) have evaluated the Melnikov
function, which is shown to satisfy the requirements for tangle formation for
any nonzero straining. They also showed the space-filling nature of the
Poincaré section and thus concluded that the transport of fluid tracers in the
Kida vortex can be chaotic.

Here, we use the same analytical tools to identify the dynamical behaviour of
heavy inertial particles in the Kida vortex. For very small values of Stokes
number, the particles are expected to behave like fluid tracers and thus can
be chaotic. On the contrary, when the strain rate is small, the inertial
particles can simply be attracted to the fixed points of the dominant
Kirchhoff vortex and thus won’t be chaotic. It is when the Stokes number and
strain rate are of the same order that the competition between these effects
takes place, and complicated things can happen. Thus, for the case of $St\ll
1$ and $s\ll 1$, we consider $St=q\,s$, where $q$ is an $\textit{O}(1)$
quantity. By substituting the same in equations (19), the system can be re-
written up to $\textit{O}(s^{2})$ as

$\displaystyle\dot{\xi}=\hat{f}_{1}(\xi,\eta)+s\,\left\\{\hat{g}_{1}(\xi,\eta;t)+q\,\hat{\phi}_{1}(\xi,\eta)\right\\}~{},$
(20a)
$\displaystyle\dot{\eta}=\hat{f}_{2}(\xi,\eta)+s\,\left\\{\hat{g}_{2}(\xi,\eta;t)+q\,\hat{\phi}_{2}(\xi,\eta)\right\\}~{}.$
(20b)

The system is perturbed from the integrable Kirchhoff vortex in $s$ and thus
may not be integrable in general. Following Bertozzi (1988), the Melnikov
functions associated with the hyperbolic fixed points of the system can be
evaluated as

$\begin{split}M_{j}^{\pm}(\phi)&=\int_{-\infty}^{\infty}\bigg{[}\hat{f}_{1}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t))\,\left\\{\hat{g}_{2}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t);t+\phi)+q\,\hat{\phi}_{2}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t))\right\\}\\\
&-\hat{f}_{2}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t))\,\left\\{\hat{g}_{1}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t);t+\phi)+q\,\hat{\phi}_{1}(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t))\right\\}\bigg{]}\,e^{-\int_{0}^{t}\textrm{tr}(Df)\,dt^{\prime}}\,dt~{},\end{split}$
(21)

where the subscript $j(=1,2)$ indicates the heteroclinic orbits
$(\xi_{j}^{\pm}(t),\eta_{j}^{\pm}(t))$ of the unperturbed system (Kirchhoff
vortex) connected to the hyperbolic fixed points for which the Melnikov
function is evaluated. Here the term
$\textrm{tr}(Df)=\frac{\partial\hat{f}_{1}}{\partial\xi}+\frac{\partial\hat{f}_{2}}{\partial\eta}$
is the trace of the Jacobian matrix of the unperturbed system evaluated along
the same heteroclinic orbit, which is not zero in general. One may note that
the non-zero trace results from writing the dynamical system in non-standard
variables $(\xi,\eta)$. However, the divergence of the corresponding velocity
vector field in the elliptical coordinate system is different from this trace,
evaluated as
$\boldsymbol{\nabla}_{\xi,\eta}\cdot[u_{\xi},u_{\eta}]=h^{2}\,\left(\frac{\partial(\hat{f}_{1}/h^{2})}{\partial\xi}+\frac{\partial(\hat{f}_{2}/h^{2})}{\partial\eta}\right)$,
indeed is zero, in accordance with the Hamiltonian nature of tracer dynamics
in Kirchhoff vortex.

After splitting the integrals, it can be seen that the time-periodic
perturbations (due to straining) $\hat{g}_{1}$ and $\hat{g}_{2}$ will result
in a periodic term in $M$ of period $\pi/\Omega_{0}$, however the time-
independent perturbations $\hat{\phi}_{1}$ and $\hat{\phi}_{2}$ will
contribute as a constant term, i.e. $M(\phi)=m_{0}(\phi)+q\,m_{1}$ where
$m_{0}(\phi+\pi/\Omega_{0})=m_{0}(\phi)$ is the periodic function and $m_{1}$
is a constant, which are dependent on the aspect ratio of the ellipse. When
the particles are inertialess (i.e., $q=0$), the resulting Melnikov function
is periodic $m_{0}(\phi)$, known to have zeros for any finite $s\ll 1$.
However, for inertial particles ($q\neq 0$), the term $m_{1}$ effectively
shifts the periodic function $m_{0}(\phi)$ up or down when plotted against
$\phi$, depending on the sign of $m_{1}$. Thus, there can exist a critical
value of $q=\lvert\max(m_{0}(\phi))/m_{1}\rvert$ such that beyond which the
Melnikov function $M(\phi)$ ceases to have transverse zeros. This indicates
that there can be a critical value of $St$ for a given strain rate $s$ above
which the heteroclinic tangles will disappear and possibly the chaos too.
Since, in our system, there is more than one heteroclinic orbit, this critical
value will be the greatest of all the critical values corresponding to each
orbit so that the complete disappearance of heteroclinic tangles will be
ensured.

Similar dynamical behaviour can also be observed in systems where two
competing physics leads to a critical parameter value at which the
appearance/disappearance of chaos occurs, for example, (i) the settling of
inertial particles in a flow-field created by a pair of like-signed vortices
(see Angilella, 2010), and (ii) the transport of passive fluid tracers in the
same flow-field near a wall (see Angilella et al., 2014). In either case, the
zeros of the sinusoidal Melnikov function are affected by the vertical shift
to the Melnikov function due to the competing parameter in the problem.

Figure 11: The Melnikov functions (a) $M_{1}^{+}(\phi)$ and (b)
$M_{2}^{+}(\phi)$ for a Kida vortex of $r_{0}=0.5$ and $\theta_{0}=0$ are
plotted against $\phi$ for various $q$ values. The abscissa is shown as a
dashed-dotted line to identify the zeros of the Melnikov functions.  Figure
12: Schematic showing the unstable and stable manifolds of saddles A and B,
respectively in the colour red and blue for various scenarios: (a) for tracers
in Kirchhoff vortex ($s=0,St=0$), (b) for tracers in Kida vortex ($s\ll
1,St=0$), (c) for inertial particles in Kirchhoff vortex ($s=0,St\ll 1$),
(d-e) inertial particles in Kida vortex, where (d) $s\ll
1,St<q_{\textrm{cr},1}\,s$ and (e) $s\ll 1,St>q_{\textrm{cr},1}\,s$.

For the Kida vortex, we explain the scenario using an example. Consider a
Kirchhoff vortex with $r_{0}=0.5$ and $\theta_{0}=0$, perturbed with a planar
extensional flow to form the Kida vortex. The heteroclinic orbits
$\textrm{H}_{1}^{\pm}$ and $\textrm{H}_{2}^{\pm}$ are numerically obtained for
the case. By substituting them in equation (21), we evaluate the Melnikov
function for each orbit. For $M_{1}^{\pm}(\phi)$, the amplitude of the
periodic function $m_{0}(\phi)$ is evaluated to be approximately $0.0082$ and
$m_{1}\approx-0.4889$. Similarly, for $M_{2}^{\pm}(\phi)$, the amplitude of
corresponding $m_{0}(\phi)$ is evaluated to be approximately $0.8267$ and
$m_{1}\approx 0.1950$. The variation of the Melnikov functions with $\phi$
over a period $\pi/\Omega_{0}$ is shown in figure 11. We can evaluate the
critical value of $q$ at which the zeros of the Melnikov function disappear as
$q_{\textrm{cr},1}\approx 0.0168$ and $q_{\textrm{cr},2}\approx 4.2395$, which
is independent of whether it is $`+^{\prime}$ or $`-^{\prime}$ type
heteroclinic orbit (due to their symmetry). Thus, we may conclude that the
heteroclinic tangles of $H_{1}^{\pm}$ disappear if $St>q_{\textrm{cr},1}\,s$
and that of $H_{2}^{\pm}$ disappear if $St>q_{\textrm{cr},2}\,s$. When
$St>q_{\textrm{cr},2}\,s$, the complete disappearance of heteroclinic tangles
occurs. In other words, there will not be any chaos when the straining is less
than the critical value $s_{\textrm{cr},2}=St/q_{\textrm{cr},2}$.

Figure 12 illustrates the schematic representation of tangles associated with
the saddles A and B and the heteroclinic connection $H_{1}^{\pm}$ in various
contexts. In each scenario, the red curve denotes the unstable manifold of
saddle A, while the blue curve represents the stable manifold of saddle B. For
tracers in the Kirchhoff vortex, as depicted in figure 12(a), the stable and
unstable manifolds coincide without any transverse intersection. With a small
strain perturbation, as shown in figure 12(b), the stable and unstable
manifolds form folds and transversely intersect, leading to heteroclinic
tangles. This tangle is responsible for the onset of chaotic dynamics of
tracer particles in a Kida vortex. We now concentrate on the effect of
particle inertia on the dynamics in the absence of external straining -
inertial particles in the Kirchhoff vortex. Figure 12(c) shows the effect of
particle inertia, where inertial particles in the Kirchhoff vortex have wide-
open, stable and unstable manifolds due to centrifugal effects. Folds are
absent as there is no strain perturbation, preventing the formation of
tangles. Figures 12(d) and (e) examine the combined effect of strain
perturbation and particle inertia. In figure 12(d), particles with weak
inertia ($St<q_{\textrm{cr},1}\,s$) allow for weak centrifuging, enabling the
transverse intersection of the folds to form tangles. Conversely, in figure
12(e), dominant particle inertia ($St>q_{\textrm{cr},1}\,s$) causes manifolds
to spread apart and prevent the intersection of the folds.

For inertial particles with $St=0.1$, the critical strain rates
$s_{\textrm{cr}}=St/q_{\textrm{cr}}$ can be obtained as
$s_{\textrm{cr},1}\approx 5.9524$ and $s_{\textrm{cr},2}\approx 0.0236$, such
that when $s\lessapprox 0.0236$, there will be no chaotic dynamics. In figures
13(a) and (b), we have shown the basin of attraction (evaluated using the full
nonlinear system of equations (15)) of $St=0.1$ particles in Kida vortex with
$r_{0}=0.5$, $\theta_{0}=0$ of $s=0.02$ and $s=0.03$ respectively. The basin
of attraction has clear boundaries for the case of $s=0.02$, indicating
regular dynamics; however, for $s=0.03$, the boundaries are mixed/fractal,
indicating an underlying chaotic set near the separatrices. Note that the
existence of zeros of the Melnikov function is a necessary criterion for chaos
but not a sufficient one. Apart from a dense set, the heteroclinic tangles may
also result in (i) a countable infinity of periodic orbits with arbitrarily
high period and (ii) an uncountable infinity of periodic orbits (Smale, 1967).
Thus, for our system, whenever it is suspected, the existence of chaos needs
to be shown using other methods, like Lyapunov exponents, frequency spectra,
and Poincaré sections. However, from the investigations of chaotic dynamics of
tracers in Kida vortex (Polvani & Wisdom, 1990; Dahleh, 1992; Kawakami &
Funakoshi, 1999), we assume the same continues to exist for inertial particles
until they have the critical inertia value.

Figure 13: The basin of attraction of $St=0.1$ particles in a Kida vortex (of
$r_{0}=0.5$ and $\theta_{0}=0$) with strain rates (a,c) $s=0.02$ indicating
smooth basin boundaries and (b,d) $s=0.03$ indicating fractal basin
boundaries. Figures are generated using the full system of equations (15)
along with equations (13).  Figure 14: (a) The critical curves which demarcate
the regions where the heteroclinic tangles can and cannot exist for
$H_{1}^{\pm}$ and $H_{2}^{\pm}$ are shown in the parameter plot of $St/s$
versus $r_{0}$, for a Kida vortex, obtained from the Melnikov analysis. The
initial orientation of Kida vortex $\theta_{0}$ will not affect these curves.
(b) Curves in the $s-St$ plane demarcate the regions where heteroclinic
tangles can and cannot exist for $H_{p}^{\pm}$, for various $r_{0}$ values.

The critical ratio between $St$ and $s$ for the occurrence of heteroclinic
tangles of both $H_{1}^{\pm}$ and $H_{2}^{\pm}$ have been evaluated
numerically for various initial aspect ratios ($r_{0}$) of the elliptic vortex
and plotted in figure 14. The initial orientation ($\theta_{0}$) of the
elliptic vortex will only phase shift the Melnikov functions; thus, these
critical values won’t be affected. The regions in the parametric plot where
the heteroclinic tangles occur are identified and marked. The region above the
red curve is where no heteroclinic tangles form, indicating complete
suppression of chaos by inertia. However, below the red curve, chaotic
transport may occur. The degree of chaos may vary across the blue curve.

It is important to remember that these critical curves are derived from a
perturbative method, assuming both $St\ll 1$ and $s\ll 1$. However, note that
the red curve reaches a value of $St/s=\textit{O}(10^{2})$ when $r_{0}=0.1$.
For $s=\textit{O}(10^{-2})$, this results in the critical
$St_{\textrm{cr},2}\sim 1$. The caveat here is that the Kirchhoff vortex’s
saddles disappear for $St>\mathcal{S}_{2}$; on the other hand, the Kida vortex
itself is destroyed beyond critical straining. Thus, the validity of our
analysis is restricted to sufficiently smaller $St$ and $s$, although their
ratio can vary within the large range as shown in figure 14(a).

####  3.2.1  Chaotic dynamics far away from the ellipse

Unlike the Kirchhoff vortex, far away from the Kida vortex (distance of
$\textit{O}(s^{-1/2})$), the decaying vortex field and external straining flow
can balance each other and form a pair of hyperbolic fixed points/saddles. The
saddles remain stationary when observed from the lab reference frame and
exhibit counter-rotation when observed from a co-rotating perspective. The
snapshots in figure 15 show the evolution of streamlines in the far-field of
the Kida vortex over a period of revolution of the ellipse for $s=0.035$. The
saddle points $A_{p}$ and $B_{p}$ are connected by a pair of heteroclinic
orbits $H_{p}^{\pm}$ and are perturbed due to the unsteady oscillations of the
central ellipse, which is shown to bring chaotic dynamics for fluid tracers,
far away from the ellipse (Kawakami & Funakoshi, 1999). As done in the
previous Section 3.2, we modify the tracer analysis by incorporating the
effect of particle inertia. The detailed analysis of evaluating the Melnikov
function for $H_{p}^{\pm}$ can be found in Appendix D. However, our findings
are similar to the case for $H_{1}^{\pm}$ and $H_{2}^{\pm}$; i.e., finite
inertia competes with strain perturbation and suppresses the chaotic
transport.

Using the Melnikov analysis, the parametric regime in the $s-St$ plane, where
$H_{p}^{\pm}$ can or cannot form a tangle, is identified and illustrated in
figure 14(b). For $r_{0}=0.5$ and $St=0.1$, the critical value of the strain
rate is found to be $s_{\textrm{cr},p}\approx 0.0341$, below which particles
cannot exhibit chaotic dynamics near $H_{p}^{\pm}$. However, it is important
to note a caveat: the Melnikov analysis for $H_{p}^{\pm}$ assumes that $St\sim
1$ and $s\ll 1$ (see Appendix D); nonetheless, the critical curve in Figure
14(b) covers even smaller values of $St$.

Figure 15: The typical streamlines of the Kida vortex, as observed in a
stationary reference frame, show the hyperbolic fixed points far away for (a)
$t=0$, (b) $t=7$, (c) $t=11$ and (d) $t=27$, shown in grey. The strain rate
here is $s=0.035$, and the period of revolution of the ellipse is $\approx
29.034$. The rotating ellipse at corresponding instants is shown as the black
dashed curve.

##  4  Large time dispersion characteristics of particles

One of the motivations for studying the dynamics of particles embedded in a
sea of coherent structures is to learn about their long-time dispersion and
the rate of capture. As we have seen in the case of the Kirchhoff vortex,
heavy inertial particles can have two fates at large times: those starting
within the basin of attraction of a fixed point will get attracted towards the
corresponding fixed point, and those starting outside will spiral off to
infinity. We can quantify the dispersion of a particle using the quantity
squared distance (SD),
$\sigma^{2}(t)=\lVert\textbf{x}(t)-\textbf{x}(0)\rVert^{2}$, which is the
squared distance the particle travelled from its initial location. The
variation of SD with time can give information about the nature of dispersion
the particle follows. SD increasing quadratically with time indicates
ballistic dispersion of the particle, while a linear in-time behaviour of SD
indicates diffusive dispersion of the particle; if SD saturates with time, it
suggests the particle is attracted to some fixed point and will have limited
dispersion. The variation of SD with time for typical particles in a Kirchhoff
vortex is numerically obtained and shown in figure 16(a). The SD saturates at
a large time for a particle starting within a basin of attraction (blue
curve), indicating the attraction towards the respective fixed point. However,
SD increases with time, indicating dispersion to infinity for a particle
starting outside the basins (black curve). The oscillations in the SD curves
indicate the spiralling nature of trajectories. Despite that, the average
growth of the black curve can be shown to scale as
$\sigma^{2}(t)\sim\sqrt{t}$, indicating that the particle dispersion is slower
than a diffusive process. The same scaling has already been obtained by
(Ravichandran et al., 2014) for the dispersion of inertial particles in a pair
of like-signed vortices. In both cases, the particles perceive a point vortex
field far away. Appendix E provides a detailed derivation for the scaling.

Figure 16: The variation of SD (evaluated in a co-rotating frame) with time
for typical inertial particles of $St=0.5$ starting within and outside the
basin of attraction of stable fixed points/limit cycles in (a) Kirchhoff
vortex of $r=0.5$ and (b) Kida vortex of $r_{0}=0.5$ and $s=0.01$. The scaling
behaviour of SD is shown using dashed lines.

In the case of a Kida vortex, as the fixed points become limit cycles,
corresponding SD saturates but with sustained periodic oscillations, as shown
in figure 16(b) in blue. On the other hand, a particle centrifuged away will
be affected by the heteroclinic orbits $H_{p}^{\pm}$. These orbits limit the
centrifuging, trap the particles along its stable and unstable manifolds, and
limit their dispersion. One may remember the far-field attractor mentioned at
the beginning of Section 3. Since the fixed points $A_{p}$ and $B_{p}$ and the
associated heteroclinic connections $H_{p}^{\pm}$ are stationary in the lab
reference frame, for a co-rotating observer, this attractor seems to be
counter-rotating, as seen in figure 9(e-h). Below the associated critical
Stokes number $St_{\textrm{cr},p}$, the particles may get transported
chaotically and above it can be transported regularly near these orbits.
However, at large time, the inertial particles can leak near through the
saddles $A_{p}$ and $B_{p}$. This behaviour can be either due to the
oscillations induced to the particle trajectory by the time-periodic
perturbation by the central rotating ellipse or due to the under-damped
oscillations of inertial particle trajectories near the stagnation points
(Nath et al., 2022). The leaked inertial particles are tempted to travel along
the extensional axis of the straining flow exponentially with time and thus
will have an enhanced dispersion later. The situation can be seen in the
snapshot in figure 9(h), where the centrifuged inertial particles are aligned
along the orbits $H_{p}^{\pm}$; the particle trajectory oscillations near
saddles $A_{p}$ and $B_{p}$, and their leakage along the extensional axis can
be observed in a co-rotating frame. The black curve in Figure 16(b) shows the
typical SD for such a particle. Due to the heteroclinic orbit barrier, the
particle has a limited dispersion at intermediate times with large amplitude
oscillations. At large times, the dispersion becomes exponentially fast (but
anisotropic) as they get caught and confined along the unstable manifold of
the heteroclinic orbit. This behaviour can be quantified by the scaling
relation $\sigma^{2}(t)\sim\exp(t\,s)$ for small $St$.

One may note that at initial times, the dispersion of all particles starting
from rest gets centrifuged from the centre, scaling identically as
$\sigma^{2}(t)\sim t^{4}$ (super-ballistic), as illustrated in the figure.
This behaviour is identical for both the Kirchhoff and Kida vortex. It occurs
due to the constant initial acceleration experienced by the particles
initially, akin to the scenario elucidated in Nath et al. (2024). Since the
particles start with zero initial velocity ($\textbf{v}(t=0)=\mathbf{0}$),
equation (14) provides the initial acceleration
$\textbf{a}_{0}=\ddot{\textbf{x}}(t=0)=\textbf{u}(\textbf{x}_{0},t=0)/St+\textbf{x}_{0}\,\Omega^{2}-\dot{\Omega}\,\hat{\textbf{e}}_{z}\times\textbf{x}_{0}$,
where $\textbf{x}_{0}=\textbf{x}(t=0)$ denotes the initial particle location.
Depending on the particle location, it may undergo an approximately constant
initial acceleration contributed by the flow and the centrifugal force in a
Kirchhoff vortex and, in addition, the Coriolis force in a Kida vortex. For
$t\ll 1$, the approximate solution can be obtained as
$\textbf{x}(t)\approx\textbf{x}_{0}+\textbf{a}_{0}\,t^{2}/2$. This solution
yields $\sigma^{2}(t)\approx\lvert\textbf{a}_{0}\rvert^{2}\,t^{4}/4$, i.e.,
$\sigma\sim t^{4}$ when $t\ll 1$.

##  5  Conclusion

We have studied the dispersion of heavy inertial particles in (i) an elliptic
Kirchhoff vortex and (ii) a Kida vortex. This study has offered insights into
particle transport in two-dimensional flows relevant to geophysical and
astrophysical scenarios, which can be approximated as the resultant flow from
isolated vortex patches and their interactions. The Kirchhoff vortex is an
elliptic patch that self-rotates and creates an unsteady flow around it. The
stability characteristics of the fixed points govern the dynamics of the
inertial particles in this flow field. When observed from a co-rotating frame
with the vortex, some particles are attracted to stable fixed points in the
flow field, while the remaining ones are centrifuged to infinity. The location
of the fixed points depends on the inertia of the particles. Thus, in a
polydisperse system, each particle will cluster into its own fixed points,
which can result in the segregation of particles according to their inertia. A
critical Stokes number exists above which all the particles only get
centrifuged to infinity.

The introduction of a weak extensional flow to the system (owing to the
interaction from other vortices), known as the Kida vortex, disrupts the
particle trajectories in the Kirchhoff vortex. Inertial particles in the Kida
vortex can have limit cycle trajectories as well as chaotic trajectories in an
extended phase space over time. However, a Melnikov analysis revealed that
particles with sufficiently high inertia are less prone to chaotic transport.

The present study offers insights into the clustering and dispersion of
particles in an ambient vortical flow field. As mentioned earlier, this has
been of interest to understanding the trapping of dust, which has the
potential for forming planetesimals. However, we have exclusively considered
one-way coupling between the particle and fluid phases in the current study.
The feedback force from the dispersed phase to the carrier phase can be
significant, particularly during the clustering phase. Allowing for two-way
coupling introduces additional instabilities, such as streaming instabilities
in the context of proto-planetary disks (Youdin & Goodman, 2005; Youdin &
Johansen, 2007). In the axisymmetric vortex system context, two-way coupling
results in the expulsion of particle clumps in a spiralling arm manner, even
in the limit of vanishing particle inertia and leads to a dramatic
destabilization of an otherwise stable isolated vortex due to a particle-
induced baroclinic torque (Shuai et al., 2022). As mentioned earlier, a non-
axisymmetric vortex system formed by like-signed point vortices exhibits
topological and dynamical similarities to the system under consideration.
While it is known that like-signed vortex pairs are prone to merging below a
critical separation (Griffiths & Hopfinger, 1987; Dritschel, 2002; Cerretelli
& Williamson, 2003), the inclusion of two-way coupling introduces novel
vorticity dynamics before their merging (Shuai et al., 2024). At the leading
order, momentum coupling from fluid to particle induces clustering dynamics,
as demonstrated in this paper. However, once clustering is achieved, the
particle concentration at those locations becomes sufficiently high that the
feedback from the particle phase to the fluid phase cannot be neglected. In
two-dimensional dusty turbulence, this feedback has produced a non-universal
energy spectrum that depends on particle inertia and mass loading (Pandey et
al., 2019). Therefore, it is crucial to approach the results presented in this
paper with caution, acknowledging that their applicability may be limited as
long as the particle feedback to the fluid is negligible, particularly in a
dilute mixture, in the earlier stages of dynamics before particle clustering.
After the onset of clustering, the feedback force would alter the morphology
of the underlying coherent structures; subsequently, delineating a causal
effect between clustering and two-way coupling may become difficult.

[Acknowledgements]A.V.S.N. thanks the Prime Minister’s Research Fellows (PMRF)
scheme, Ministry of Education, Government of India. A.R. and A.V.S.N.
acknowledge the support of the Centre for Atmospheric and Climate Sciences
(CACS) at IIT Madras. A.R. and A.V.S.N. also acknowledge the support of the
Geophysical Flows Lab (GFL) at IIT Madras.

[Author ORCIDs]A. V. S. Nath, https://orcid.org/0000-0003-2144-2978; A. Roy,
https://orcid.org/0000-0002-0049-2653

##  Appendix A  Governing equations in the co-rotating Cartesian reference
frame

The dynamic equations (3) governing the transport of heavy inertial particles
in the Kirchhoff vortex are obtained from the Maxey-Riley equation equation
(2) by writing it in elliptic coordinates ($\xi,\eta$), which is helpful in
analytical exercises. However, for numerical simulations, it is convenient to
use these equations in the co-rotating Cartesian reference frame ($x,y$),
which are

$\displaystyle\dot{v}_{x}=\frac{u_{x}-v_{x}}{St}+x\,\Omega^{2}+2\,\Omega\,v_{y}~{},\quad\dot{v}_{x}=v_{x}~{},$
(22a)
$\displaystyle\dot{v}_{y}=\frac{u_{y}-v_{y}}{St}+y\,\Omega^{2}-2\,\Omega\,v_{x}~{},\quad\dot{v}_{y}=v_{y}~{}.$
(22b)

The flow velocity field in the stationary reference frame can be obtained from
the stream function as $u_{x}^{\prime}=\frac{\partial\psi^{\prime}}{\partial
y^{\prime}}$ and $u_{y}^{\prime}=-\frac{\partial\psi^{\prime}}{\partial
x^{\prime}}$, which is unsteady. The velocity field in the co-rotating frame
can be obtained from the following transformation

$\displaystyle
u_{x}=u_{x}^{\prime}\,\cos\theta+u_{y}^{\prime}\,\sin\theta+\Omega\,y~{},$
(23a) $\displaystyle
u_{y}=-u_{x}^{\prime}\,\sin\theta+u_{y}^{\prime}\,\cos\theta-\Omega\,x~{}.$
(23b)

For Kirchhoff vortex, we get

$\displaystyle u_{x}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}k^{-1}\,\sin\eta\,\\{-\cosh\xi+(1+k^{2}\,\Omega)\,\sinh\xi\\},&\quad\quad\tanh\xi>r\vspace{0.25cm}\\\
-k\,\Omega\,r^{-1}\,\sinh\xi\,\sin\eta,&\quad\quad\tanh\xi<r\\\
\end{array}\right.$ (24c) $\displaystyle u_{y}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}-k^{-1}\,\cos\eta\,\\{\sinh\xi+(-1+k^{2}\,\Omega)\,\cosh\xi\\},&\quad\,\,\tanh\xi>r\vspace{0.25cm}\\\
k\,\Omega\,r\,\cosh\xi\,\cos\eta,&\quad\,\,\,\tanh\xi<r~{},\\\
\end{array}\right.$ (24f)

where the variables in both the co-ordinate system are related as
$x=k\,\cosh\xi\,\cos\eta$ and $y=k\,\sinh\xi\,\sin\eta$. The dynamical system
remains four-dimensional; however, now in the phase space of variables
$x,y,v_{x}$ and $v_{y}$. The fixed points of the system can be obtained by
solving for $(\dot{x},\dot{y},\dot{v}_{x},\dot{v}_{y})=(0,0,0,0)$. Trivially
it implies that $\overline{v}_{x}=\overline{v}_{y}=0$ and the fixed point
positions $\overline{\textbf{x}}=(\overline{x},\overline{y})$ are solutions of
the vector equation $\textbf{u}(\textbf{x})+St\,\Omega^{2}\,\textbf{x}=0$.
Using numerical solvers like ‘fsolve’ in MATLAB or Newton-Raphson method, one
can obtain the fixed points and verify that they exactly match the analytical
expressions (7). For the Cartesian equations (22), the entries in the
stability matrix will be much simpler than equation (9) and are

$\displaystyle\mathsfbi{J}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\
\Omega^{2}+\frac{1}{St}\,\frac{\partial u_{x}}{\partial
x}&\frac{1}{St}\,\frac{\partial u_{x}}{\partial y}&-\frac{1}{St}&2\,\Omega\\\
\frac{1}{St}\,\frac{\partial u_{y}}{\partial
x}&\Omega^{2}+\frac{1}{St}\,\frac{\partial u_{y}}{\partial
y}&-2\,\Omega&-\frac{1}{St}\end{pmatrix}$ (25) $\displaystyle=$
$\displaystyle\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\
\Omega^{2}+\frac{h^{2}\,\sin
2\eta}{2\,k^{2}\,St}&\frac{\Omega}{St}+\frac{2-h^{2}\,\sinh
2\xi}{2\,k^{2}\,St}&-\frac{1}{St}&2\,\Omega\\\
-\frac{\Omega}{St}+\frac{2-h^{2}\,\sinh
2\xi}{2\,k^{2}\,St}&\Omega^{2}-\frac{h^{2}\,\sin
2\eta}{2\,k^{2}\,St}&-2\,\Omega&-\frac{1}{St}\end{pmatrix}~{}.$

Note that the eigenvalues of the Jacobian matrix in equation (25) and equation
(9) will be the same since they are invariants of the same dynamical system.
However, the components of the eigenvectors differ by the coordinate
transformation.

##  Appendix B  Perturbation analysis for the limit cycles for a Kida vortex -
entries of $\mathsfbi{K}$, $\mathsfbi{L}$, $\mathsfbi{M}$ and $N$

The strain perturbation to the system affects the fixed points
($\overline{\xi},\overline{\eta}$) of a Kirchhoff vortex. Using the method of
IJzermans & Hagmeijer (2006), we assume the fixed points are perturbed to a
limit cycle and solve for its coordinates
($\overline{\xi}+s\,\xi^{\prime}(t)+\mathcal{O}(s^{2}),\overline{\eta}+s\,\eta^{\prime}(t)+\mathcal{O}(s^{2})$)
using equation 17. The entries in the coefficient matrices depend on the fixed
points of the Kirchhof vortex and are listed below:

$\mathsfbi{L}=\left[\begin{array}[]{cc}1&-2\,St\,\Omega_{0}\\\
2\,St\,\Omega_{0}&1\\\
\end{array}\right],\quad\mathsfbi{K}=\left[\begin{array}[]{cc}K_{11}&K_{12}\\\
K_{21}&K_{22}\\\
\end{array}\right],\quad\textrm{and}\,\,\,\mathsfbi{M}=\left[\begin{array}[]{cc}M_{11}&M_{12}\\\
M_{21}&M_{22}\\\ \end{array}\right].$ (26)

Outside the ellipse ($\tanh\xi>r$), corresponding to the fixed points A or B
or C or D, denoted as $(\overline{\xi},\overline{\eta})$, the elements of
$\mathsfbi{K}$, $\mathsfbi{M}$ and $N$ are

$\displaystyle K_{11}$ $\displaystyle=$
$\displaystyle\overline{h}^{2}\,\Omega_{0}\left(\sin(2\,\overline{\eta})-e^{2\,\overline{\xi}}\,St\,\Omega_{0}\right)~{},$
(27a) $\displaystyle K_{12}$ $\displaystyle=$
$\displaystyle\,\overline{h}^{2}\,\cos(2\,\overline{\eta})\left(\,\Omega_{0}-k_{0}^{-2}\,e^{-2\,\overline{\xi}}\right)~{},$
(27b) $\displaystyle K_{21}$ $\displaystyle=$
$\displaystyle\overline{h}^{2}\left(\Omega_{0}\,\cosh(2\,\overline{\xi})-k_{0}^{-2}\,e^{-2\,\overline{\xi}}\,\cos(2\,\overline{\eta})\right)~{},$
(27c) $\displaystyle K_{22}$ $\displaystyle=$
$\displaystyle\overline{h}^{2}\,\Omega_{0}\left(St\,\Omega_{0}\,\cos(2\,\overline{\eta})-l_{1}^{-1}\,e^{-2\,\overline{\xi}}\sin(2\,\overline{\eta})\right)~{},$
(27d) $\displaystyle M_{11}$ $\displaystyle=$
$\displaystyle\frac{\overline{h}^{2}}{4}\,\left\\{\sin(2\,\overline{\eta})\,(l_{1}-\cosh(2\,\overline{\xi}))-\,l_{3}\,St\,\Omega_{0}\sinh(2\,\overline{\xi})\right\\}~{},$
(27e) $\displaystyle M_{12}$ $\displaystyle=$
$\displaystyle\frac{\overline{h}^{2}}{4}\,\left\\{\sinh(2\,\overline{\xi})(\Lambda_{0}-\cos(2\,\overline{\eta}))-2\,l_{1}\,St\,\Omega_{0}\sin(2\,\overline{\eta})\right\\}~{},$
(27f) $\displaystyle M_{21}$ $\displaystyle=$
$\displaystyle\frac{\overline{h}^{2}}{4}\,\left\\{\sinh(2\,\overline{\xi})(l_{1}-\cos(2\,\overline{\eta}))+l_{3}\,St\,\Omega_{0}\sin(2\,\overline{\eta})\right\\}~{},$
(27g) $\displaystyle M_{22}$ $\displaystyle=$
$\displaystyle-\frac{\overline{h}^{2}}{4}\,\left\\{\sin(2\,\overline{\eta})(\Lambda_{0}-\cosh(2\,\overline{\xi}))+2\,l_{1}\,St\,\Omega_{0}\sinh(2\,\overline{\xi})\right\\}~{},$
(27h) $\displaystyle N_{1}$ $\displaystyle=$
$\displaystyle\frac{\overline{h}^{2}}{4}\,\cos(2\,\theta_{0})\left(2\,k_{0}^{-2}\,\sin(2\,\overline{\eta})+l_{3}\,St\,\Omega_{0}\sinh(2\,\overline{\xi})\right)~{},$
(27i) $\displaystyle N_{2}$ $\displaystyle=$
$\displaystyle\frac{\overline{h}^{2}}{4}\,\cos(2\,\theta_{0})\left(2\,k_{0}^{-2}\,\sinh(2\,\overline{\xi})-l_{3}\,St\,\Omega_{0}\sin(2\,\overline{\eta})\right)~{},$
(27j)

where $\overline{h}=(\cosh^{2}\overline{\xi}-\cos^{2}\overline{\eta})^{-1/2}$,
$l_{1}=(1-r_{0})/(1+r_{0})$ and
$l_{3}=(1-4\,r_{0}+r_{0}^{2})/(1-r_{0}^{2})=2\,l_{1}-\Lambda_{0}$. Inside the
ellipse ($\tanh\xi<r$), corresponding to the fixed point at the origin O, the
corresponding elements are

$\displaystyle K_{11}$ $\displaystyle=$ $\displaystyle r_{0}\,\Omega_{0}~{},$
(28a) $\displaystyle K_{12}$ $\displaystyle=$ $\displaystyle
K_{21}=-St\,\Omega_{0}^{2}~{},$ (28b) $\displaystyle K_{22}$ $\displaystyle=$
$\displaystyle-\Omega_{0}/r_{0}~{},$ (28c) $\displaystyle M_{11}$
$\displaystyle=$ $\displaystyle M_{12}=M_{21}=M_{22}=N_{1}=N_{2}=0~{}.$ (28d)

##  Appendix C  Terms required for the evaluation of the Melnikov function for
weakly inertial particles in a Kida vortex

Outside the ellipse $\tanh\xi>r_{0}+s\,r_{1}+\textit{O}(s^{2})$,

$\displaystyle\hat{f}_{1}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{2}\,\left(-\Omega_{0}+k_{0}^{-2}\,e^{-2\xi}\right)\,\sin
2\eta~{},$ (29a) $\displaystyle\hat{f}_{2}$ $\displaystyle=$
$\displaystyle-\frac{h^{2}}{2}\,\left(\Omega_{0}\,\sinh
2\xi+k_{0}^{-2}\,\left[-1+e^{-2\xi}\,\cos 2\eta\right]\right)~{},$ (29b)
$\displaystyle\hat{g}_{1}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{4}\,\left\\{\mathcal{F}(\xi,\eta;r_{0},\Omega_{0}\,t+\theta_{0})-2\,\frac{\Omega_{0}}{r_{0}}\,r_{1}\,\left(l_{1}-l_{2}\,e^{-2\xi}\right)\,\sin
2\eta\right\\}~{},$ (29c) $\displaystyle\hat{g}_{2}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{4}\,\left\\{\mathcal{G}(\xi,\eta;r_{0},\Omega_{0}\,t+\theta_{0})-2\,\frac{\Omega_{0}}{r_{0}}\,r_{1}\,\left(l_{1}\,\sinh
2\xi+l_{2}\,[-1+e^{-2\xi}\,\cos 2\eta]\right)\right\\}~{},$ (29d)
$\displaystyle\hat{\phi}_{1}$ $\displaystyle=$ $\displaystyle
h^{2}\,k_{0}^{-2}\,e^{-2\xi}\,\left(k_{0}^{-2}-\Omega_{0}\,\cos 2\eta\right)$
(29e) $\displaystyle\hat{\phi}_{2}$ $\displaystyle=$
$\displaystyle-h^{2}\,k_{0}^{-2}\,e^{-2\xi}\,\Omega_{0}\,\sin 2\eta~{}.$ (29f)

Inside the elliptic region $\tanh\xi<r_{0}+s\,r_{1}+\textit{O}(s^{2})$,

$\displaystyle\hat{f}_{1}$ $\displaystyle=$
$\displaystyle\frac{h^{2}}{2}\,\left(-\Omega_{0}+\frac{1}{2}\,\left\\{1-l_{1}\,\cosh
2\xi\right\\}\right)\,\sin 2\eta~{},$ (30a) $\displaystyle\hat{f}_{2}$
$\displaystyle=$
$\displaystyle\frac{h^{2}}{2}\,\left(-\Omega_{0}+\frac{1}{2}\,\left\\{1-l_{1}\,\cos
2\eta\right\\}\right)\,\sinh 2\xi~{},$ (30b) $\displaystyle\hat{g}_{1}$
$\displaystyle=$
$\displaystyle\frac{h^{2}}{4}\,\left\\{\mathcal{F}(\xi,\eta;r_{0},\Omega_{0}\,t+\theta_{0})-2\,\frac{\Omega_{0}}{r_{0}}\,r_{1}\,\left(l_{1}-\cosh
2\xi\right)\,\sin 2\eta\right\\}~{},$ (30c) $\displaystyle\hat{g}_{2}$
$\displaystyle=$
$\displaystyle\frac{h^{2}}{4}\,\left\\{\mathcal{G}(\xi,\eta;r_{0},\Omega_{0}\,t+\theta_{0})-2\,\frac{\Omega_{0}}{r_{0}}\,r_{1}\,\left(l_{1}-\cos
2\eta\right)\,\sinh 2\xi\right\\}~{},$ (30d) $\displaystyle\hat{\phi}_{1}$
$\displaystyle=$ $\displaystyle\frac{h^{2}\,l_{1}}{2}\,\sinh
2\xi\,\left(k_{0}^{-2}-\Omega_{0}\,\cos 2\eta\right)~{},$ (30e)
$\displaystyle\hat{\phi}_{2}$ $\displaystyle=$
$\displaystyle-\frac{h^{2}\,l_{1}}{2}\,\sin
2\eta\,\left(k_{0}^{-2}-\Omega_{0}\,\cosh 2\xi\right)~{},$ (30f)

where $l_{2}=(1+r_{0}^{2})/(1-r_{0})^{2}$. The functions $\mathcal{F}$ and
$\mathcal{G}$ are defined as

$\displaystyle\mathcal{F}(\xi,\eta;r,\theta)$ $\displaystyle=$
$\displaystyle\cosh 2\xi\,\sin 2\eta\,\cos 2\theta+\sinh 2\xi\,\cos
2\eta\,\sin 2\theta$ (31a) $\displaystyle-$ $\displaystyle\Lambda\,\left(\sin
2\eta\,\cos 2\theta+\sinh 2\xi\,\sin 2\theta\right)~{},$
$\displaystyle\mathcal{G}(\xi,\eta;r,\theta)$ $\displaystyle=$
$\displaystyle\sinh 2\xi\,\cos 2\eta\,\cos 2\theta-\cosh 2\xi\,\sin
2\eta\,\sin 2\theta$ (31b) $\displaystyle+$ $\displaystyle\Lambda\,\left(\sin
2\eta\,\sin 2\theta-\sinh 2\xi\,\cos 2\theta\right)~{}.$

Expressions for $\hat{f}_{1}$ to $\hat{g}_{2}$ for a Kida vortex of
$\theta_{0}=\pi/4$ can also be found in Kawakami & Funakoshi (1999), but with
a sign mistake in $l_{1}$ and $l_{2}$.

##  Appendix D  Evaluation of the Melnikov function in the far-field of a Kida
vortex

The streamfunction for a Kida vortex in the stationary reference frame
($\psi^{\prime}=\psi^{\prime}_{v}+\psi^{\prime}_{e}$) can be written in far-
field (far from the central elliptic vortex) using polar coordinates as

$\psi^{\prime}=-\frac{1}{2}\,\log R-\frac{s\,R^{2}}{4}\,\cos
2\Theta-\frac{k^{2}}{16\,R^{2}}\,\cos(2\,\Theta-2\,\theta)+\textit{O}(R^{-4})~{},$
(32)

where the Cartesian coordinates in stationary frame $(x^{\prime},y^{\prime})$
are related to the polar coordinates in stationary frame $(R,\Theta)$ as
$x^{\prime}=R\,\cos\Theta$ and $y^{\prime}=R\,\sin\Theta$. The dominant point
vortex field balances with external straining to create an integrable system
far-field. However, the central rotating ellipse induces a time-periodic
perturbation to this system. We re-scale radial and time coordinates with
strain rate such that the balance between the point vortex and extensional
flow is apparent and the perturbation due to the elliptic vortex is clearly
visible. Following Kawakami & Funakoshi (1999), the appropriate scaling is
$\tilde{R}=R\,\sqrt{s}$, $\tilde{t}=t\,s$, which will reduce equation (32) to
the stream function in scaled variables as

$\psi^{\prime}=-\frac{1}{2}\,\log\tilde{R}-\frac{\tilde{R}^{2}}{4}\,\cos
2\Theta-\frac{s\,k^{2}}{16\,\tilde{R}^{2}}\,\cos(2\,\Theta-2\,\theta)+\textit{O}(s^{2})~{}.$
(33)

The first two terms on the right-hand side indicate the integrable system
formed by the point vortex and the extensional flow. The streamlines of this
Hamiltonian system are shown in figure 17, along with the heteroclinic orbits
$H_{p}^{\pm}$ and far-field fixed points $A_{p}$ and $B_{p}$ as marked.
Observe the similarity of the far-field streamlines to the snapshots in figure
15. The third term on the right-hand side is the perturbation due to the
central rotating elliptic vortex patch. Using the slow manifold equations for
inertial particles in the stationary reference frame

$\displaystyle\dot{x}^{\prime}=u_{x}^{\prime}-St\,\left[\frac{\partial
u_{x}^{\prime}}{\partial t}+u_{x}^{\prime}\,\frac{\partial
u_{x}^{\prime}}{\partial x^{\prime}}+u_{y}^{\prime}\,\frac{\partial
u_{x}^{\prime}}{\partial y^{\prime}}\right]+\textit{O}(St^{2})~{},$ (34a)
$\displaystyle\dot{y}^{\prime}=u_{y}^{\prime}-St\,\left[\frac{\partial
u_{y}^{\prime}}{\partial t}+u_{x}^{\prime}\,\frac{\partial
u_{y}^{\prime}}{\partial x^{\prime}}+u_{y}^{\prime}\,\frac{\partial
u_{y}^{\prime}}{\partial y^{\prime}}\right]+\textit{O}(St^{2})~{},$ (34b)

Figure 17: The streamlines of the flow-field result from the balance between
the irrotational (point) vortex and external extensional flow. The flow
pattern mimics the far-field flow around a Kida vortex (c.f figure 15). The
hyperbolic fixed points and connecting heteroclinic orbits are marked.

we may write the modified kinematic equations for a particle in polar
coordinates as

$\displaystyle\dot{R}$ $\displaystyle=$
$\displaystyle\frac{1}{R}\,\frac{\partial\psi^{\prime}}{\partial\Theta}+\frac{St}{R}\,\left\\{-\frac{\partial^{2}\psi^{\prime}}{\partial\Theta\,\partial
t}+\left(\frac{\partial\psi^{\prime}}{\partial
R}\right)^{2}+\frac{1}{R}\,\left[\frac{\partial^{2}\psi^{\prime}}{\partial\Theta^{2}}\,\frac{\partial\psi^{\prime}}{\partial
R}-\frac{\partial^{2}\psi^{\prime}}{\partial
R\,\partial\Theta}\,\frac{\partial\psi^{\prime}}{\partial\Theta}\right]+\frac{1}{R^{2}}\,\left(\frac{\partial\psi^{\prime}}{\partial\Theta}\right)^{2}\right\\}$
(35a) $\displaystyle+$ $\displaystyle\textit{O}(St^{2})~{},$
$\displaystyle\dot{\Theta}$ $\displaystyle=$
$\displaystyle-\frac{1}{R}\,\frac{\partial\psi^{\prime}}{\partial
R}+\frac{St}{R}\,\left\\{\frac{\partial^{2}\psi^{\prime}}{\partial R\,\partial
t}+\frac{1}{R}\,\left[\frac{\partial^{2}\psi^{\prime}}{\partial
R^{2}}\,\frac{\partial\psi^{\prime}}{\partial\Theta}-\frac{\partial^{2}\psi^{\prime}}{\partial
R\,\partial\Theta}\,\frac{\partial\psi^{\prime}}{\partial
R}\right]+\frac{1}{R^{2}}\,\frac{\partial\psi^{\prime}}{\partial\Theta}\,\frac{\partial\psi^{\prime}}{\partial
R}\right\\}$ (35b) $\displaystyle+$ $\displaystyle\textit{O}(St^{2})~{}.$

After substituting the small $s$ expansion of $k$ and $\theta$ (see Section
3.1) in equation (35), evolution equations for the scaled coordinates can be
obtained as

$\displaystyle\frac{d\tilde{R}}{d\tilde{t}}$ $\displaystyle=$
$\displaystyle\tilde{f}_{1}+s\,\left\\{\tilde{g}_{1}+St\,\tilde{\phi}_{1}\right\\}+\textit{O}(s^{2},St^{2}\,s,s^{2}\,St)~{},$
(36a) $\displaystyle\frac{d\Theta}{d\tilde{t}}$ $\displaystyle=$
$\displaystyle\tilde{f}_{2}+s\,\left\\{\tilde{g}_{2}+St\,\tilde{\phi}_{2}\right\\}+\textit{O}(s^{2},St^{2}\,s,s^{2}\,St)~{},$
(36b)

where

$\displaystyle\tilde{f}_{1}$ $\displaystyle=$
$\displaystyle\frac{\tilde{R}}{2}\,\sin 2\Theta~{},$ (37a)
$\displaystyle\tilde{f}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2\,\tilde{R}^{2}}+\frac{1}{2}\,\cos 2\Theta~{},$ (37b)
$\displaystyle\tilde{g}_{1}$ $\displaystyle=$
$\displaystyle\,\frac{k_{0}^{2}}{8\,\tilde{R}^{3}}\,\sin(2\,\Theta-2\,\theta_{0}-2\,\Omega_{0}\,\tilde{t}/s)~{},$
(37c) $\displaystyle\tilde{g}_{2}$ $\displaystyle=$
$\displaystyle-\frac{k_{0}^{2}}{8\,\tilde{R}^{4}}\,\cos(2\,\Theta-2\,\theta_{0}-2\,\Omega_{0}\,\tilde{t}/s)~{},$
(37d) $\displaystyle\tilde{\phi}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{4\,\tilde{R}^{3}}\,\left\\{1-\tilde{R}^{4}+k_{0}^{2}\,\Omega_{0}\,\cos(2\,\Theta-2\,\theta_{0}-2\,\Omega_{0}\,\tilde{t}/s)\right\\}~{},$
(37e) $\displaystyle\tilde{\phi}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2\,\tilde{R}^{2}}\,\left\\{\sin
2\Theta+\frac{k_{0}^{2}\,\Omega_{0}}{2\,\tilde{R}^{2}}\,\sin(2\,\Theta-2\,\theta_{0}-2\,\Omega_{0}\,\tilde{t}/s)\right\\}~{}.$
(37f)

Note that the above system of equations reduces to the form given in Kawakami
& Funakoshi (1999) for the case of $St=0$ and $\theta_{0}=0$ (with a sign
mistake in $\tilde{g}_{1}$ and $\tilde{g}_{2}$). For $s\ll 1$ and $St\ll 1$,
the Melnikov function associated with the far-field hyperbolic fixed points
$A_{p}$ and $B_{p}$ can be evaluated as

$\begin{split}M_{p}^{\pm}(\phi,s)=\int_{-\infty}^{\infty}\bigg{[}\tilde{f}_{1}(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t}))\,\left\\{\tilde{g}_{2}\left(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t});\frac{\tilde{t}}{s}+\phi\right)+St\,\tilde{\phi}_{2}\left(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t});\frac{\tilde{t}}{s}+\phi\right)\right\\}\\\
-\tilde{f}_{2}(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t}))\,\left\\{\tilde{g}_{1}\left(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t});\frac{\tilde{t}}{s}+\phi\right)+St\,\tilde{\phi}_{1}\left(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t});\frac{\tilde{t}}{s}+\phi\right)\right\\}\bigg{]}\,e^{-\int_{0}^{\tilde{t}}\textrm{tr}(D\tilde{f})\,dt^{\prime}}\,d\tilde{t}~{},\end{split}$
(38)

where the integration is performed along the heteroclinic orbits $H_{p}^{\pm}$
of the unperturbed integrable system parametrised by
$(\tilde{R}_{0}^{\pm}(\tilde{t}),\Theta_{0}^{\pm}(\tilde{t}))$. Here

$\textrm{tr}(D\tilde{f})=\frac{\partial\tilde{f}_{1}}{\partial\tilde{R}}+\frac{\partial\tilde{f}_{2}}{\partial\Theta}=-\frac{1}{2}\,\sin
2\Theta$ (39)

is the trace of the Jacobian matrix of the unperturbed system, which is not
zero in general. From equation (38), one can infer that $M_{p}^{\pm}$ is
periodic in $\phi$ with period $\pi/\Omega_{0}$. Unlike $M_{1}^{\pm}$ and
$M_{2}^{\pm}$, $M_{p}^{\pm}$ has explicit dependancy on the strain rate $s$.
However, it does not affect the validity of the Melnikov analysis and related
theorems as stated by Kawakami & Funakoshi (1999). The variation of the
Melnikov function $M_{p}^{\pm}$ with $\phi$ for various $s$ values for
$St=0.1$ particles in a Kida vortex is shown in figure 18. Since $M_{p}^{\pm}$
depends on $s$ implicitly, it will result in a nonlinear relation between the
critical values of $s$ and $St$, as can be seen in the figure 14(b).

Figure 18: The Melnikov function $M_{p}^{+}$ corresponding to a Kida vortex of
$r_{0}=0.5$ and $\theta_{0}=0$ for particles of $St=0.1$ is plotted against
$\phi$ for different values of strain rate $s$. The abscissa is shown as a
dashed-dotted line to identify the zeros of the Melnikov function.

##  Appendix E  Scaling law for large-time dispersion with time

The inertial particles centrifuged away to infinity by the Kirchhoff vortex
spiral outward and form a ring-like structure, as we see in the simulation
result (figure 3). At large times, these particles reach far away from the
central ellipse. Thus, the background flow they perceive can be approximated
to that created by an irrotational vortex (point vortex). Using the polar
coordinates $(R,\Theta)$ about the origin, the flow-field for the point vortex
can be represented in stationary frame as

$U_{R}=0~{},\quad U_{\Theta}=\frac{1}{2R}~{}.$ (40)

Then, the dynamics of a heavy inertial particle in the point vortex is given
by the Maxey-Riley equations in polar coordinates as (see Ravichandran &
Govindarajan, 2015)

$\displaystyle\dot{v}_{R}+\frac{v_{R}}{St}=\frac{v_{\Theta}^{2}}{R}~{},\quad\dot{R}=v_{R},$
(41a)
$\displaystyle\dot{(R\,v_{\Theta})}+\frac{R\,v_{\Theta}}{St}=\frac{R\,U_{\Theta}}{St}~{},\quad\dot{\Theta}=\frac{v_{\Theta}}{R}~{}.$
(41b)

Solving the first equation among equations (41b) along with equation (40)
yields the azimuthal velocity of particle as

$v_{\Theta}=\frac{1}{2\,R}\left(1+C_{1}\,e^{-t/St}\right)~{},$ (42)

where $C_{1}$ is an integration constant. The form of the solution says that
the azimuthal velocity of the particle approaches that of fluid tracers at a
large time, i.e. $v_{\Theta}=U_{\Theta}$ as $t\rightarrow\infty$. Substituting
this result in equations (41a) and solve for $R$ asymptotically
($\dot{v}_{R}\ll v_{R}$) yields, at large time $R(t)\sim(t\,St)^{1/4}$. I.e.,
in terms of the dispersion parameter SD,
$\sigma^{2}(t)=R(t)^{2}\sim\sqrt{t\,St}$ at large time.

The same scaling law can also be obtained from an analysis using the method of
characteristics for the evolution of the number density field of particles.
The single fluid model is an Eulerian model for particle-laden flows, where
the suspended particles are treated as a field, characterised by their number
density. The model is valid if the particles are only weakly inertial and
their concentration is just right. The model yields the evolution equation for
the axisymmetric number density field $n(t,R)$ as

$\frac{\partial n}{\partial t}+\frac{St}{R}\,\frac{\partial}{\partial
R}(U_{\Theta}^{2}\,n)=0~{}.$ (43)

The solution for this partial differential equation in a point vortex field
equations (40) can be obtained using the method of characteristics as,

$n(t,R)=\frac{R^{2}}{\sqrt{R^{4}-t\,St}}\,n_{0}\left(\left\\{R^{4}-t\,St\right\\}^{1/4}\right)~{},$
(44)

where $n(0,R)=n_{0}(R)$ is the initial number density distribution. Not that
the scaling $R^{4}\sim t\,St$ automatically appears in this solution.

## References

  * Angilella (2010) Angilella, Jean-Régis 2010 Dust trapping in vortex pairs. Physica D: Nonlinear Phenomena 239 (18), 1789–1797.
  * Angilella et al. (2014) Angilella, Jean-Régis, Vilela, Rafael D & Motter, Adilson E 2014 Inertial particle trapping in an open vortical flow. Journal of fluid mechanics 744, 183–216.
  * Armitage (2020) Armitage, Philip J 2020 Astrophysics of planet formation. Cambridge University Press.
  * Avni & Dagan (2022) Avni, Orr & Dagan, Yuval 2022 Dynamics of evaporating respiratory droplets in the vicinity of vortex dipoles. International Journal of Multiphase Flow 148, 103901.
  * Babiano et al. (1987) Babiano, Armando, Basdevant, Claude, Legras, Bernard & Sadourny, Robert 1987 Vorticity and passive-scalar dynamics in two-dimensional turbulence. Journal of Fluid Mechanics 183, 379–397.
  * Bartello et al. (1994) Bartello, Peter, Métais, Olivier & Lesieur, Marcel 1994 Coherent structures in rotating three-dimensional turbulence. Journal of Fluid Mechanics 273, 1–29.
  * Basdevant et al. (1981) Basdevant, Claude, Legras, Bernard, Sadourny, Robert & Béland, Me 1981 A study of barotropic model flows: intermittency, waves and predictability. Journal of Atmospheric Sciences 38 (11), 2305–2326.
  * Batchelor & Nitsche (1994) Batchelor, GK & Nitsche, JM 1994 Expulsion of particles from a buoyant blob in a fluidized bed. Journal of Fluid Mechanics 278, 63–81.
  * Bec (2003) Bec, Jeremie 2003 Fractal clustering of inertial particles in random flows. Physics of fluids 15 (11), L81–L84.
  * Bec (2005) Bec, Jérémie 2005 Multifractal concentrations of inertial particles in smooth random flows. Journal of Fluid Mechanics 528, 255–277.
  * Benzi et al. (1988) Benzi, R, Patarnello, S & Santangelo, P 1988 Self-similar coherent structures in two-dimensional decaying turbulence. Journal of Physics A: Mathematical and General 21 (5), 1221\.
  * Bertozzi (1987) Bertozzi, Andrea Louise 1987 An extension of the Smale-Birhoff homoclinic theorem, Melnikov’s method, and chaotic dynamics in incompressible fluids. PhD thesis, Princeton University.
  * Bertozzi (1988) Bertozzi, Andrea Louise 1988 Heteroclinic orbits and chaotic dynamics in planar fluid flows. SIAM Journal on Mathematical Analysis 19 (6), 1271–1294.
  * Cerretelli & Williamson (2003) Cerretelli, C & Williamson, CHK 2003 The physical mechanism for vortex merging. Journal of Fluid Mechanics 475, 41–77.
  * Chaplygin (1899) Chaplygin, SA 1899 On a pulsating cylindrical vortex. Trans. Phys. Sect. Imperial Moscow Soc. Friends Natural Sciences 10 (1), 13–22.
  * Chavanis (2000) Chavanis, Pierre-Henri 2000 Trapping of dust by coherent vortices in the solar nebula. Astronomy and Astrophysics 356, 1089–1111.
  * Crisanti et al. (1992) Crisanti, Andrea, Falcioni, Massimo, Provenzale, Antonello, Tanga, Paolo & Vulpiani, Angelo 1992 Dynamics of passively advected impurities in simple two-dimensional flow models. Physics of Fluids A: Fluid Dynamics 4 (8), 1805–1820.
  * Dagan (2021) Dagan, Yuval 2021 Settling of particles in the vicinity of vortex flows. Atomization and Sprays 31 (11).
  * Dahleh (1992) Dahleh, Marie D 1992 Exterior flow of the kida ellipse. Physics of Fluids A: Fluid Dynamics 4 (9), 1979–1985.
  * Dorgan & Loth (2007) Dorgan, AJ & Loth, E 2007 Efficient calculation of the history force at finite reynolds numbers. International journal of multiphase flow 33 (8), 833–848.
  * Dritschel (1990) Dritschel, David G 1990 The stability of elliptical vortices in an external straining flow. Journal of Fluid Mechanics 210, 223–261.
  * Dritschel (1995) Dritschel, David G 1995 A general theory for two-dimensional vortex interactions. Journal of Fluid Mechanics 293, 269–303.
  * Dritschel (2002) Dritschel, David G 2002 Vortex merger in rotating stratified flows. Journal of Fluid Mechanics 455, 83–101.
  * Dritschel & JuáRez (1996) Dritschel, David G & JuáRez, Manuel De La Torre 1996 The instability and breakdown of tall columnar vortices in a quasi-geostrophic fluid. Journal of Fluid Mechanics 328, 129–160.
  * Dufek (2016) Dufek, Josef 2016 The fluid mechanics of pyroclastic density currents. Annual Review of Fluid Mechanics 48, 459–485.
  * Eames & Gilbertson (2004) Eames, I & Gilbertson, MA 2004 The settling and dispersion of small dense particles by spherical vortices. Journal of Fluid Mechanics 498, 183–203.
  * Elhmaïdi et al. (1993) Elhmaïdi, Dalila, Provenzale, Antonello & Babiano, Armando 1993 Elementary topology of two-dimensional turbulence from a Lagrangian viewpoint and single-particle dispersion. Journal of Fluid Mechanics 257, 533–558.
  * Falkovich et al. (2002) Falkovich, Gregory, Fouxon, A & Stepanov, MG 2002 Acceleration of rain initiation by cloud turbulence. Nature 419 (6903), 151–154.
  * Fitzpatrick (2012) Fitzpatrick, Richard 2012 An introduction to celestial mechanics. Cambridge University Press.
  * Fornberg (1977) Fornberg, Bengt 1977 A numerical study of 2-D turbulence. Journal of Computational Physics 25 (1), 1–31.
  * Griffiths & Hopfinger (1987) Griffiths, RW & Hopfinger, EJ 1987 Coalescing of geostrophic vortices. Journal of Fluid Mechanics 178, 73–97.
  * Hofmann et al. (2018) Hofmann, Lutz, Rieck, Bastian & Sadlo, Filip 2018 Visualization of 4D vector field topology. In Computer Graphics Forum, , vol. 37, pp. 301–313. Wiley Online Library.
  * Hunt et al. (2007) Hunt, JCR, Delfos, René, Eames, Ian & Perkins, Richard J 2007 Vortices, complex flows and inertial particles. Flow, turbulence and combustion 79, 207–234.
  * IJzermans & Hagmeijer (2006) IJzermans, Rutger HA & Hagmeijer, Rob 2006 Accumulation of heavy particles in N-vortex flow on a disk. Physics of fluids 18 (6).
  * Kawakami & Funakoshi (1999) Kawakami, Akihiko & Funakoshi, Mitsuaki 1999 Chaotic motion of fluid particles around a rotating elliptic vortex in a linear shear flow. Fluid dynamics research 25 (4), 167.
  * Kida (1981) Kida, Shigeo 1981 Motion of an elliptic vortex in a uniform shear flow. Journal of the Physical Society of Japan 50 (10), 3517–3520.
  * Kirchhoff (1876) Kirchhoff, Gustav 1876 Mechanik. Vol. 1 of his Vorlesungen über mathematische Physik. Teubner.
  * Kok et al. (2012) Kok, Jasper F, Parteli, Eric JR, Michaels, Timothy I & Karam, Diana Bou 2012 The physics of wind-blown sand and dust. Reports on progress in Physics 75 (10), 106901.
  * Lamb (1945) Lamb, Horace 1945 Hydrodynamics. Dover, New York.
  * Legras & Dritschel (1991) Legras, Bernard & Dritschel, David G 1991 The elliptical model of two-dimensional vortex dynamics. I: The basic state. Physics of Fluids A: Fluid Dynamics 3 (5), 845–854.
  * Legras et al. (2001) Legras, Bernard, Dritschel, David G & Caillol, Philippe 2001 The erosion of a distributed two-dimensional vortex in a background straining flow. Journal of Fluid Mechanics 441, 369–398.
  * Lovalenti & Brady (1993) Lovalenti, Phillip M & Brady, John F 1993 The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motion at small reynolds number. Journal of Fluid Mechanics 256, 561–605.
  * Love (1893) Love, AEH 1893 On the stability of certain vortex motions. Proceedings of the London Mathematical Society 1 (1), 18–43.
  * Marcu et al. (1995) Marcu, B, Meiburg, E & Newton, PK 1995 Dynamics of heavy particles in a Burgers vortex. Physics of Fluids 7 (2), 400–410.
  * Maxey (1987) Maxey, Martin R 1987 The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. Journal of fluid mechanics 174, 441–465.
  * Maxey & Riley (1983) Maxey, Martin R & Riley, James J 1983 Equation of motion for a small rigid sphere in a nonuniform flow. The Physics of Fluids 26 (4), 883–889.
  * McWilliams (1984) McWilliams, James C 1984 The emergence of isolated coherent vortices in turbulent flow. Journal of Fluid Mechanics 146, 21–43.
  * Meleshko & Van Heijst (1994) Meleshko, VV & Van Heijst, GJF 1994 On Chaplygin’s investigations of two-dimensional vortex structures in an inviscid fluid. Journal of Fluid Mechanics 272, 157–182.
  * Mitchell & Rossi (2008) Mitchell, TB & Rossi, LF 2008 The evolution of Kirchhoff elliptic vortices. Physics of Fluids 20 (5).
  * Moore & Saffman (1971) Moore, DW & Saffman, P Gr 1971 Structure of a line vortex in an imposed strain. In Aircraft Wake Turbulence and Its Detection: Proceedings of a Symposium on Aircraft Wake Turbulence held in Seattle, Washington, September 1–3, 1970. Sponsored jointly by the Flight Sciences Laboratory, Boeing Scientific Research Laboratories and the Air Force Office of Scientific Research, pp. 339–354. Springer.
  * Nath et al. (2022) Nath, Anu VS, Roy, Anubhab, Govindarajan, Rama & Ravichandran, Sivaramakrishnan 2022 Transport of condensing droplets in Taylor-Green vortex flow in the presence of thermal noise. Physical Review E 105 (3), 035101.
  * Nath et al. (2024) Nath, Anu VS, Roy, Anubhab, Ravichandran, Sandhanakrishnan & Govindarajan, Rama 2024 Irregular dependence on Stokes number, and nonergodic transport, of heavy inertial particles in steady laminar flows. Physical Review Fluids 9 (1), 014302.
  * Ott (2002) Ott, Edward 2002 Chaos in dynamical systems. Cambridge university press.
  * Pandey et al. (2019) Pandey, Vikash, Perlekar, Prasad & Mitra, Dhrubaditya 2019 Clustering and energy spectra in two-dimensional dusty gas turbulence. Physical Review E 100 (1), 013114.
  * Polvani & Wisdom (1990) Polvani, Lorenzo M & Wisdom, J 1990 Chaotic Lagrangian trajectories around an elliptical vortex patch embedded in a constant and uniform background shear flow. Physics of Fluids A: Fluid Dynamics 2 (2), 123–126.
  * Raju & Meiburg (1997) Raju, N & Meiburg, E 1997 Dynamics of small, spherical particles in vortical and stagnation point flow fields. Physics of Fluids 9 (2), 299–314.
  * Ravichandran & Govindarajan (2015) Ravichandran, S & Govindarajan, Rama 2015 Caustics and clustering in the vicinity of a vortex. Physics of Fluids 27 (3), 033305.
  * Ravichandran et al. (2014) Ravichandran, S, Perlekar, Prasad & Govindarajan, Rama 2014 Attracting fixed points for heavy particles in the vicinity of a vortex pair. Physics of Fluids 26 (1), 013303.
  * Reinaud et al. (2003) Reinaud, Jean Noel, Dritschel, David Gerard & Koudella, CR 2003 The shape of vortices in quasi-geostrophic turbulence. Journal of Fluid Mechanics 474, 175–192.
  * Rom-Kedar et al. (1990) Rom-Kedar, V, Leonard, A & Wiggins, S 1990 An analytical study of transport, mixing and chaos in an unsteady vortical flow. Journal of Fluid Mechanics 214, 347–394.
  * Saffman (1995) Saffman, Philip G 1995 Vortex dynamics. Cambridge university press.
  * Sapsis & Haller (2010) Sapsis, Themistoklis & Haller, George 2010 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.
  * Shaw (2003) Shaw, Raymond A 2003 Particle-turbulence interactions in atmospheric clouds. Annual Review of Fluid Mechanics 35 (1), 183–227.
  * Shuai et al. (2022) Shuai, Shuai, Dhas, Darish Jeswin, Roy, Anubhab & Kasbaoui, M Houssem 2022 Instability of a dusty vortex. journal of fluid mechanics 948, A56.
  * Shuai et al. (2024) Shuai, Shuai, Roy, Anubhab & Kasbaoui, M Houssem 2024 The merger of co-rotating vortices in dusty flows. Journal of Fluid Mechanics 981, A27.
  * Siggia (1981) Siggia, Eric D 1981 Numerical study of small-scale intermittency in three-dimensional turbulence. Journal of Fluid Mechanics 107, 375–406.
  * Smale (1967) Smale, Stephen 1967 Differentiable dynamical systems. Bulletin of the American mathematical Society 73 (6), 747–817.
  * Tang (1987) Tang, Yun 1987 Nonlinear stability of vortex patches. Transactions of the American Mathematical Society 304 (2), 617–638.
  * Tanga et al. (1996) Tanga, P, Babiano, A, Dubrulle, B & Provenzale, A 1996 Forming planetesimals in vortices. Icarus 121 (1), 158–170.
  * Vallis (2017) Vallis, Geoffrey K 2017 Atmospheric and oceanic fluid dynamics. Cambridge University Press.
  * Vincent & Meneguzzi (1991) Vincent, Albert & Meneguzzi, Maria 1991 The spatial structure and statistical properties of homogeneous turbulence. Journal of Fluid Mechanics 225, 1–20.
  * Viúdez (2021) Viúdez, A 2021 Robust and unstable axisymmetric vortices, including neutral vortices, of a new two-dimensional vortex family. Physics of Fluids 33 (5).
  * Wan (1986) Wan, Yieh-Hei 1986 The stability of rotating vortex patches. Communications in mathematical physics 107 (1), 1–20.
  * Wiggins (1990) Wiggins, S 1990 Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer-Verlag.
  * Wilkinson & Mehlig (2005) Wilkinson, M & Mehlig, Bernhard 2005 Caustics in turbulent aerosols. Europhysics Letters 71 (2), 186.
  * Xu & Krasny (2023) Xu, Ling & Krasny, Robert 2023 Dynamics of elliptical vortices with continuous profiles. Physical Review Fluids 8 (2), 024702.
  * Youdin & Johansen (2007) Youdin, Andrew & Johansen, Anders 2007 Protoplanetary disk turbulence driven by the streaming instability: linear evolution and numerical methods. The Astrophysical Journal 662 (1), 613.
  * Youdin & Goodman (2005) Youdin, Andrew N & Goodman, Jeremy 2005 Streaming instabilities in protoplanetary disks. The Astrophysical Journal 620 (1), 459.
  * Zhao et al. (2024) Zhao, Zilong, Guo, Zhiwei, Zuo, Zhigang & Qian, Zhongdong 2024 Trapping of inertial particles in a two-dimensional unequal-strength counterrotating vortex pair flow. Physical Review Fluids 9 (2), 024307.