Instability of a dusty vortex

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


Introduction
Vortical flows with heavy dispersed particles occur in diverse scenarios, both of natural and engineering origins (Balachandar & Eaton 2010;Guazzelli & Morris 2011). To cite a few examples, they include atmospheric funnel-type vortices (such as 'dust devil') (Bluestein et al. 2004), centrifugal separation devices, swirling biomass combustors, and aircraft trailing vortices with condensation droplets (Paoli & Shariff 2016;Paoli et al. 2008). The rudimentary 2 S. Shuai, D. J. Dhas, A. Roy, and M. H. Kasbaoui understanding of the process is that particles that are denser than the surrounding fluid drift away from intense vortical regions at a rate controlled by particle inertia. However, even in elementary vortical flows such as columnar vortices, the dispersion of inertial particles (micron-sized solid spheres, or liquid droplets, sufficiently small to remain spherical under the action of surface tension) is difficult to predict accurately. This uncertainty is due to models and simulations often omitting the particle feedback force. Yet, the latter may be significant even in dilute conditions (particle volume fractions is 10 −5 < < 10 −3 ) provided that the density ratio is high ( / = (1000)), e.g. when the carrier fluid is gas. In this study, we shall show that particle feedback arising from two-way momentum coupling activates a novel vortex instability, and this instability is controlled by particle inertia and vortex structure.
To discuss the stability properties in a vortex flow, we adopt the Rankine vortex as a typical vortex tube, i.e., one with uniform vorticity within the core and zero outside (Saffman 1995). Many measurements of atmospheric phenomena show that the Rankine vortex is suitable for describing velocity structure of naturally occuring vortical flows. Doppler radar measurements prove that tornados are well approximated by the Rankine vortex profile (Bluestein et al. 2003). A mesocyclone is another example of atmospheric flow that is well approximated by a Rankine vortex with an inner core of solid rotation as large as 5km (Brown et al. 2005).
The simple velocity profile of the Rankine vortex makes stability investigations amenable to analytical treatments, thus, often helping in better identifying the underlying physics. Lord Kelvin (Thomson 1880; Lamb 1993) analyzed the linear stability of an incompressible Rankine vortex. He showed that for 2D perturbations characterized by mode number , only one real eigenvalue exists, indicating that the perturbation will propagate with no growth or decay, corresponding to neutral stability. Besides these neutrally stable discrete modes, a Rankine vortex also supports a continuous spectrum of modes whose combination may lead to algebraic growth for short times (Roy & Subramanian 2014). The inclusion of additional physics can destabilize a vortex column. For instance, acoustic waves radiating to infinity can destabilize a Rankine vortex by extracting energy from the vortex core (Broadbent & Moore 1979). Analogous destabilization mechanisms have been shown to exist for shallow water flows (Ford 1994) and stratified flows (Schecter & Montgomery 2004;Le Dizès & Billant 2009), where outgoing waves draw energy from the mean flow. Recently a new instability mechanism has been found for a Rankine vortex in dilute polymer solutions (Roy et al. 2022). The instability occurs due to the resonant interaction of a pair of elastic shear waves aided by the differential convection of the irrotational shearing flow outside the vortex core. Compressibility effects may also destabilize a Rankine vortex as shown by Sozou & Lighthill (1987).
The configuration of a radially stratified vortex has special relevance, since, as we shall discuss in detail later, it bears analogy with a particle-laden vortex in the limit of zero particle inertia. If density increases monotonically with radius, Fung & Kurzweg (1975) showed that the flow is stable to both axisymmetric and non-axisymmetric modes. In contrast, a heavy-cored vortex may become unstable, as shown by Sipp et al. (2005a), through two mechanisms, one due to a centrifugal instability involving short-axial wavelength modes, and the second due to a Rayleigh-Taylor instability involving two-dimensional modes. Similarly, Joly et al. (2005) found that heavy-cored two-dimensional vortices are subject to a Rayleigh-Taylor instability. They attributed the basic mechanism to baroclinic vorticity generation caused by the misalignment between density gradient and centripetal acceleration. Using numerical simulations, Joly et al. (2005) showed that the unstable modes lead to the formation of spiraling arms in the density and vorticity fields that roll-up eventually as nonlinear effects emerge. Dixit & Govindarajan (2011) analyzed the case of a radially stratified vortex with a density jump. They showed that both heavy-cored and light-cored vortices could be unstable if the density jump is sufficiently steep. They proposed that the instability is triggered by the wave-interaction mechanism between the Kelvin wave located at the vortex core edge and counter-propagating internal waves caused by the centrifugal force (see Carpenter et al. (2011);Balmforth et al. (2012) for discussions on instabilities originating from interaction between waves riding on vorticity/density jumps). When the vorticity and density jumps are sufficiently separated, they act as "independent" neutral modes; the two interfacial displacements are in phase with each other and out of phase with the radial velocity. Thus, the radial velocity neither aid nor retard the interfacial displacement. When the jumps are located close to each other, the interfacial displacements are out of phase, and the radial velocity reinforces the interfacial displacement at the density jump resulting in a growing mode. Dixit & Govindarajan (2011) also found that a heavy-cored vortex is stabilized when the density jump is placed in a region immediately outside the vortex core. Our analysis in section §3.3 will reveal that, when we neglect particle inertia but account for non-zero mass loading, the system analyzed in the present study mimics a non-Bousinessq fluid. Thus, the observed instability has an interesting analogy with the physical mechanisms elucidated in the investigations mentioned above.
When inertial particles heavier than the carrier fluid disperse in the flow, they tend to migrate from regions of high rotation rate to regions of high strain rate. This mechanism known, as preferential concentration (Squires & Eaton 1991;Wang & Maxey 1993), results from a slip velocity between particles and flow, which increases with increasing particle inertia. Preferential concentrations leads to the clustering of inertial particles around vortical cores.
Although particle dispersion in vortical flows exhibits rich dynamics, the disperse phase does not alter the fundamental stability of the flow unless two-way coupling is considered. Recently there has been a strong interest in the modulation of flow instabilities due to the feedback force from small, heavy particles. Magnani et al. (2021) considered the effect of particle inertia in two-way coupled dusty Rayleigh-Taylor turbulence. They observed that the system behaves similarly to an equivalent denser fluid for low inertia particles. When particle inertia increases, turbulent mixing gets delayed. The non-monotonic role of the disperse phase was also observed by Sozza et al. (2022) in their stability study of two-way coupled dusty Kolmogorov flow. When particle inertia is weak, the instability is enhanced. However, for larger Stokes number, particles can both stabilize or destabilize the Kolmogorov flow, with a non-monotonic dependence on the mass fraction. Inertial particles can also trigger instabilities that would not otherwise exist in particle-free flows. Kasbaoui et al. (2015) showed that the interplay between preferential concentration and gravitation settling triggers a non-modal instability in simple shear flow. This instability may bootstrap additional Rayleigh-Taylor and particle-trajectory crossing instabilities (Kasbaoui et al. 2019a) and plays a role in the attenuation of turbulent fluctuations in particle-laden homogeneously sheared turbulence (Kasbaoui et al. 2019b;Kasbaoui 2019).
In the context of vortical flows, few studies considered the impact of two-way coupling on the stability of such flows. Marshall (2005) studied the effect of turbulence on dispersed particles in a vortex column under one-way coupling, i.e., neglecting the particle feedback force. He showed that inertial particles initially inside the vortex core are driven out, forming concentrated ring-like structures. Turbulence breaks these structures into smaller sections. Ravichandran & Govindarajan (2015) studied the formation of particle clusters around isolated and ensembles of two-dimensional vortices in one-way coupling. They showed that particles initially within a critical radius, that depends on the vortex circulation and particle response time, form caustics, thus, potentially leading to higher collision rates. When twoway coupling is considered, the dispersion of inertial particles may lead to a modulation of 4 S. Shuai, D. J. Dhas, A. Roy, and M. H. Kasbaoui the carrier flow. Druzhinin (1994) showed that the outward ejection of inertial particles from a particle-rich vortex core causes the attenuation of vorticity in the core. Considering twodimensional axisymmetric perturbations, he showed that the particles form a ring-shaped cluster that expands outwardly, akin to a concentration wave, at a rate controlled by Stokes number St Γ = / , where is particle response time, and is fluid characteristic timescale. These prior observations were corroborated recently by Shuai & Kasbaoui (2022) in two-way coupled Eulerian-Lagrangian simulations of particle-laden Lamb-Oseen vortex. The authors found that particle rings grow on a timescale = /St Γ , and proposed an expression for their expansion rate. Similar to (Druzhinin 1994), Shuai & Kasbaoui (2022) observed that particle-feedback on the fluid drives a faster decay of the vortex tube. Further, the observation of naturally emerging azimuthal perturbations in the vorticity field suggests that an instability activated by two-way coupling may be at play. These observations motivate the present study on the stability of a particle-laden vortex.
This paper is organized as follows. Section §2 presents evidence of instability in a two-way coupled particle-laden vortex. The governing equation and linear stability analysis for small inertia particles are presented in §3. In §4, we compare the analytical solutions with the results from Eulerian-Lagrangian simulations of the two-way coupled system. Concluding remarks are given in §5.

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

Eulerian-Lagrangian method
The Eulerian-Lagrangian simulations presented in this work are based on the volume-filtered formulation (Anderson & Jackson 1967;Jackson 2000;Capecelatro & Desjardins 2013). In this formulation, the fluid phase is treated in the Eulerian frame, whereas the particle phase is treated in the Lagrangian frame, i.e., each particle is individually tracked.
The mass and momentum conservation equations for the carrier phase in the semi-dilute regime, are given by the incompressible Navier-Stokes equations, where u is the fluid velocity, is the pressure, is the fluid density, and is the fluid viscosity. The term F represents momentum exchange between particles and fluid, which is expressed as where = − I + (∇ + ∇ )/2 is the total fluid stresses, is the particle volume fraction, u is the Eulerian particle velocity, and (·) | indicates fluid quantities at the particle locations. The first term in (2.3) is the stresses exerted by the undisturbed flow at the particle location. The second term is the stresses caused by the presence of particles which are represented using Stokes drag. When the density ratio is large ( / 1), as is the case presently, Stokes drag dominates the momentum exchange.
From a scaling analysis of the drag force in (2.3), it is clear that the mass loading = Focus on Fluids articles must not exceed this page length 5 / determines the strength of the coupling. Thus, if the mass loading is vanishingly small, the particle phase has little effect on the stability of the Rankine vortex. Conversely, if the mass loading is large, the interaction between two phases triggers a significant source or sink of energy that may enhance or attenuate perturbations in the flow (Kasbaoui et al. 2015).
The particle phase is described in the Lagrangian frame. The motion of the -th particle is given by (Maxey & Riley 1983) x where x , and u , are the position and velocity of the " "-th particle, respectively, = 2 /(18 ) is the particle response time, and the particle diameter. In this study, gravity is ignored in order to focus on inertial effects. Other interactions, including collisions, are negligible due to the large density ratio and low volume fraction. Note that in the equations above, the instantaneous particle volume fraction and Eulerian particle velocity field are obtained from Lagrangian quantities using where = 3 /6 is the particle volume, represents a Gaussian filter kernel of size = 3Δ , where Δ is the grid spacing (Capecelatro & Desjardins 2013). This computational framework was recently applied by the authors in a configuration similar to the one considered here. In Shuai & Kasbaoui (2022), the dynamics of a particleladen Lamb-Oseen vortex at moderate Stokes numbers are investigated using the Eulerian-Lagrangian methodology presented here. Readers interested in further details about the numerical approach are referred to this study.

Illustration of the instability
To illustrate the unstable dynamics activated by two-way coupling, we consider a Rankine vortex at the circulation Reynolds number Re Γ = Γ/(2 ) = 1000, laden with mono-disperse particles having diameter at the Stokes number St Γ = / = 0.025, average volume fraction = 1.2 10 −3 , density ratio / = 830, and mass loading = 1. Here, Γ corresponds to the vortex circulation, is the initial vortex core radius, = 2 2 /Γ is the characteristic fluid velocity, = 2 /(18 ) is the particle response time, and is the particle diameter. The particles are seeded randomly within the core region with velocities matching the fluid velocity at their respective locations. We use the Eulerian-Lagrangian framework described above on uniform Cartesian grid with resolution /Δ = 50. The simulations performed here are termed "pseudo-2D", meaning that the axial direction is resolved with one grid point and taken periodic with a thickness Δ = 3 . This allows the definition of volumetric quantities such as particle volume and volume fraction. Despite considering pseudo-2D simulations, a large number of particles must be tracker, here = 376612, due to the large scale ratio / = 1400. The remaining numerical parameters are as in Shuai & Kasbaoui (2022).  In addition to two-way coupled simulations, we perform one-way coupled simulations where the momentum exchange term (2.3) is purposely set to zero. In this way, comparing oneway and two-way simulations reveals the effect of particle feed-back on the flow dynamics. Figure 1 shows the evolution of the isocontours of axial vorticity normalized by the initial center vorticity = ( = 0, = 0). Successive snapshots are given between / = 0 and 48. When the particle feedback force is neglected, i.e., in one-way coupling, we observe that the vorticity field remain axisymmetric at all times. Although some viscous diffusion can be observed near the vorticity jump ∼ , the vorticity magnitude within the core remains mostly flat as viscous effects are too small to cause significant diffusion within the time frame 0 / 48. Contrary to the case of one-way coupling, we observe significant distortion of the flow field when two-way coupling is accounted for. The vorticity field quickly loses its cylindrical symmetry due to the emergence of azimuthal perturbations. The latter grow into vorticity filaments by / ∼ 12 that are gradually convected away from the vortex core. By / ∼ 48, we observe several well-established spiral arms emanating from the vortex core. In this region, vorticity assumes a diffuse profile in contrast to the nearly flat profile seen in one-way coupled simulations.
Particle dispersion is also significantly impacted by the interaction between the two phases. Figure 2 shows the iso-contours of particle volume fraction. When the particle feedback is neglected, the particles gradually accumulate into a ring-shaped cluster of diameter about 1.5× . These clustering dynamics were reported in the earlier work of Druzhinin (Druzhinin 1994(Druzhinin , 1995, Marshall (2005), and later by Shuai & Kasbaoui (2022). The outward particle ejection is due to preferential concentration (Squires & Eaton 1991) whereby inertial particles are driven out of strongly vortical regions. Like the carrier flow, the particle distribution is axisymmetric only in one-way coupled simulations. In contrast to these prior observations, two-way coupling leads to a loss of symmetry along with fast and wide dispersion of the particles. The two-way coupled simulations in figure 2 reveal the presence of azimuthal perturbations superimposed on the particle patch at / = 12. The perturbations extend into spiraling particle filaments reminiscent of density-driven radial Rayleigh-Taylor instability observed by Dixit &Govindarajan (2011) andJoly et al. (2005). These perturbations destroy the ring structure observed in one-way coupled simulations.
The simulations shown in figure 1 and 2 suggest that semi-dilute inertial particles trigger a Figure 2: Iso-contours of normalized particle volume fraction ( / ) with St Γ = 0.025 and Re Γ = 1000, for one-way coupling and two-way coupling at five non-dimensional times / . modal instability. The latter cannot be attributed to purely hydrodynamic nor purely granular effects since the Rankine vortex is neutrally stable and particle-particle interactions are neglected. Rather, it is an instability that arises from the two-way momentum exchange between the two phases. The instability observed here represents a distinct mechanism from the one studied by Kasbaoui et al. (2015). There, the authors show that the interplay between particle settling and preferential concentration activates a non-modal instability in two-way coupled shear flows. It results in the formation of sheets of concentrated particle clusters that rotate to progressively align with the direction of the shear flow. Since particle settling is required for perturbations to grow, the resulting growth rates depend on the gravitational acceleration. In the present work, a modal instability arises regardless of gravitational effects.

Linear stability analysis for weakly inertial particles
In this section, we establish a theoretical grounding to this novel two-phase instability through rigorous linear stability analysis (LSA). The analysis is intended to reveal the mechanisms that cause the 2D particle-laden vortex flow to be unstable, and their functional dependence on the non-dimensional numbers at hand. This study is restricted to the assumption of small, but non-zero particle inertia such that St Γ < 0.1.
In the following, we adopt an Eulerian-Eulerian description of the particle-laden flow based on the Two-Fluid model (Marble 1970;Druzhinin 1995;Jackson 2000;Kasbaoui et al. 2015Kasbaoui et al. , 2019b. In this framework, the conservation equations for the fluid phase and particle 8 S. Shuai, D. J. Dhas, A. Roy, and M. H. Kasbaoui phase in the two-fluid model read ∇ · = 0, (3.1) where is the fluid velocity, is the particle velocity, is the particle volume and = / is the particle number density. Equations (3.1) and (3.2) express mass and momentum conservation for the fluid phase, respectively. Equations (3.3) and (3.4) express mass and momentum conservation for the particle phase. In these equations, it is assumed that the particles experience a hydrodynamic force equal to Stokes drag. The two-phases are coupled through Stokes drag as seen in equations (3.2) and (3.4). If the particle feedback force is dropped from the fluid momentum equation (3.2), the evolution of the fluid becomes decoupled from that of the particles. The fluid evolves as in single-phase, and the flow may not become unstable as shown in §2 and Michalke & Timme (1967). The momentum exchange between the two phases is critical for the development of an instability.
Compared to the Eulerian-Lagrangian framework presented in §2, the form of the differential equations in the Two-Fluid model is more suitable for theoretical analysis using LSA. Still the two approaches yield qualitatively and quantitatively similar evolution of particle-laden flows in the semi-dilute regime provided that particle inertia is small. Further details on the comparison of the two formulations can be found in (Kasbaoui et al. 2019a). Thus, we restrict the analysis to cases where St Γ 1. Under the assumption of weakly inertial particles, it is possible to further simplify the governing equations by adopting the fast equilibrium approximation (see work by Balachandar and coworkers (Ferry & Balachandar 2002, 2001Rani & Balachandar 2003) and Maxey (Maxey 1987)). Provided that particle inertia is low enough, it is possible to express the particle velocity field as where it is seen that the particles deviate from the fluid streamlines by a small inertial correction. Owing to the preferential concentration mechanism, particles suspended in a vortex do not follow the closed-loop streamlines of the carrier flow. Instead, they tend to migrate outward radially at a rate controlled by their Stokes number. Next, we combine equation (3.5) with the conservation equations (3.1), (3.2) and (3.3). Using , Γ/(2 ), and the core number density 0 as the reference length scale, velocity scale, and number density, respectively, the following non-dimensionalized equations are obtained ∇ · * = 0, (3.6) (1 + * ) * * + * · ∇ * = −∇ * + where the stared variables denote non-dimensional quantities, Re Γ = Γ/(2 ) is the 9 (a) (b) Figure 3: Configuration studied using linear stability analysis: (a) a Rankine vortex with core radius is seeded with particles within the region ; (b) the resulting preferential concentration term is negative within the vortex core which causes inertial particles in this region to be ejected outward.
circulation Reynolds number, St Γ = 2 /Γ is the circulation Stokes number, and = 0 / is the mass loading. Equation (3.7) is derived by injecting the small Stokes expansion (3.5) in the momentum exchange term in (3.2). Likewise, equation (3.8) is obtained from (3.3) by replacing the particle velocity field with the expansion (3.5). For ease of notation, we drop the stars in the rest of the analysis.
In this form, the particle phase is coupled with the fluid through the preferential concentration term appearing in the right hand side of (3.8). Since Squires & Eaton (1991), preferential concentration term has been widely studied (Squires & Eaton 1991;Wang & Maxey 1993;Batchelor & Nitsche 1991;Druzhinin 1995). It is understood that this term relates to compressibility effects of the particle phase since ∇ · − ∇ : ∇ for low inertia particles. This term may be further expressed as ∇ : ∇ = 2 − 2 , where 2 and 2 represent the second invariants of the fluid strain rate tensor = (∇u + ∇u )/2 and rotation tensor = (∇u − ∇u )/2 (Squires & Eaton 1991). When the local strain exceeds the local rotation, this term is positive leading to particle accumulation. Conversely, when this term is negative, the particles are expelled. As a result, inertial particles are progressively depleted from rotational regions, such as vortex cores. The formation of a ring cluster shown in figure 2 is a direct effect of preferential concentration.
The number density term in (3.7) suggests that the particles may lead to non-Boussinesq dynamics similar to those observed in stratified fluids. In particular, Dixit & Govindarajan (2010) showed that density gradients could lead to a destabilization of 2D vortex. Sipp et al. (2005b) showed that density gradients produce two distinct kinds of instability, namely centrifugal instability (CTI) which mainly affects axisymmetric eigenmodes and Rayleigh-Taylor instability (RTI) which mainly affects non-axisymmetric eigenmodes. If particle inertia is negligible, one may ignore the preferential concentration term, and treat the suspension as a fluid with effective density eff = + . In such case, the suspension would be susceptible to the instabilities of variable density vortex flows aforementioned. These effects are revisited for the case of non-inertial particles in §3.3, and then extended to the case of weakly inertial particles in §3.4.

Base state
In the remainder of the manuscript, we perform a linear stability analysis for a base state initially comprising monodisperse particles seeded randomly within a disk of nondimensional size and a Rankine vortex. A schematic of the configuration is given in figure  3a. The non-dimensional base fluid azimuthal velocity field , and number density fields are given by (3.10) Note that the task of defining a base state requires careful consideration, when the suspended particles are inertial. Even when the base flow is steady, the particle distribution may be unsteady due to the outward ejection of particles caused by the preferential concentration term in equation (3.8). The latter is plotted in figure 3b for the Rankine vortex in equation (3.9). The preferential concentration term is negative within the core ( < 1) showing that it acts as a sink in this region. Conversely, the preferential concentration term is positive away from the core ( > 1), corresponding to regions where the particles will accumulate. Thus, the base number density in (3.10) is expected to vary as time progresses.
The time evolution of the base state may be obtained by solving equations (3.6), (3.7), and (3.8). Considering an inviscid and axisymmetric state ( ()/ = 0), the equations dictating the evolution of the base state are = 0, (3.11) Following Druzhinin (1994), it is possible to obtain an analytical solution for the unsteady number density by applying the method of characteristics. Denoting the initial particle number density in (3.10) by 0 = ( , 0), the unsteady number density is written as where cr = ( 4 − 1)/4St Γ and 0 ( , ) is obtained from the equation 4 log 0 + 4 0 = 4 − 4St Γ . Despite the explicit time-dependence of the base number density field , a quasi-steady assumption may be employed for weakly inertial particles. Under this assumption, the base state is assumed frozen at = 0. This assumption is justified provided that perturbations grow on a time scale that is much faster than the characteristic evolution time of the number density field. From solution (3.13), it can be seen that the timescale for the development of number density inhomogeneity is (2022) showed from Eulerian-Lagrangian simulations that the time required to eject most particles from the core of a vortex is ∼ 3 . This timescale may be significantly larger than the flow inertial timescale when St Γ 1. For such weakly inertial particles with low inertia, particle clustering is slow compared to inertial effects of the carrier vortex flow. Since we expect the instability to develop on a timescale comparable with the flow inertial timescale, the assumption of quasi-steady number density is justified.

Linearized equations
We assume that the base state discussed previously is subject to small perturbations and such that the total number density and fluid velocity are = + and = + . Taking the inviscid limit (Re Γ → ∞) and linearizing equations (3.6) and (3.7), the fluid perturbation is subject to the following equations 1 ( ) + 1 = 0, (3.14) (1 + ) where equation (3.14) represents mass conservation in cylindrical coordinates, and equations (3.15) and (3.16) represent the radial and azimuthal fluid momentum conservation respectively. Introducing the base state angular velocity Ω = / and axial vorticity = ( ( )/ )/ , equations (3.15) and (3.16) become The evolution of the perturbation axial vorticity, i.e., is obtained by combing equations (3.17) and (3.18), (3.20) The number density perturbation is solution to the following equation, .
Here, D denotes the derivative with respect to and = 1 + . Note that equations (3.23) and (3.24) retain the two-way coupling required for an instability.
Using the assumption of quasi-steady base state discussed above, we perform a modal stability analysis with the base state being frozen in time by considering perturbations that have an exponential dependence on time. Under this consideration, the perturbations can be written as The above pair of equations can be treated as an eigenvalue problem to obtain the stability parameters.

Limit of inertia-less particles (St Γ = 0) and non-Boussinesq effects
As previously discussed in this section, with the absence of particle inertia (St Γ = 0), the preferential concentration term drops out of the governing equations. Also, the base state in the case of St Γ = 0 remains constant over time as shown in equations (3.11) and (3.12). A simplified linear equations can be obtained as . (3.29) The above linearized equations are identical to the equations for a vortex with a radially stratified density written by Dixit & Govindarajan (2011). Without the particle inertia, the problem reduces to that of a vortex with a density stratification induced by particles. Following the approach of Dixit & Govindarajan (2011), it is possible solve equations (3.28) and (3.29) analytically and obtain a dispersion relation for a Rankine vortex. This is done by computing the perturbation velocity and number density fields separately in regions of < 1, 1 < < and > . The constants subsequently obtained are found using the continuity condition at the jump locations = 1 and = . The dispersion relation, thus obtained, is written as (3.33) where At = /(2 + ) is the Atwood number. The latter is a non-dimensional number that commonly arises in density stratified flows. Here, considering the particle-fluid mixture as a fluid with effective density eff = 1 + ( / ) = 1 + , we see that in the core eff ( ) = 1 + = max for 1, while eff = 1 = min away from the core ( 1). Thus, the Atwood number is At = ( max − min )/( max + min ) which measures the relative magnitude of the density jump between the inner and outer regions. Figures 4 (a) and (b) show contours of growth rates for wavenumbers = 2 and 4, respectively, in the At − plane. The white space labelled "S" represents the region of stability. Thus there exists a critical radius, cr , such that instability would only be observed when the density jump is located beyond cr . Using the property of the cubic equation (3.30) one can obtain a relation for cr as a function of and At, which on further expanding for At 1 yields, where Kelvin = √︁ /( − 1) is the critical radius corresponding to the discrete Kelvin mode (Roy & Subramanian 2014). As is evident from the complete numerical solution and the above small-At asymptotic expression, the region of stability shrinks as the wavenumber increases, suggesting that higher modes are more unstable. In addition, higher values of growth rates are achieved with increasing values of Atwood number and with the value of closer to one. This suggests that the configuration at = 1 with infinitely heavy core ( → ∞) is the most unstable configuration.
To further study the influence of the Atwood number on the critical wavenumber, we investigate the specific case of = 1, for which the dispersion relation reduces to (1 + ) 2 (2 + ) . (3.36) The above expression shows that increasing the Atwood number, or equivalently the mass loading, promotes the instability by making lower wavenumber modes unstable. In the limit case where → ∞, corresponding to At → 1, the critical wavenumber becomes cr → 1 showing that all modes are unstable. Conversely, when → 0, none of the modes are unstable since cr → ∞. Thus, we recover the known behavior for an unladen Rankine vortex ( → 0). From expression (3.35), the growth rate of unstable modes ( cr ) is given by From this relationship, we draw two conclusions. First, a naturally developing instability, i.e., unforced, will tend to emerge at high wavenumbers since ( ) increases with . Viscous effects are expected to curb ( ) beyond a certain mode number . Because viscous effects are not accounted for in the present formulation, it is not possible to estimate a priori which mode would emerge naturally. The second conclusion is that the growth rate increases with mass loading. The maximum growth rate obtained when

Effects of particle inertia
Next, we consider the effects of non-zero particle inertia (St Γ ≠ 0) on the stability of the particle-laden Rankine vortex. For this, we first perform a modal stability analysis by solving equations (3.26) and (3.27) with a frozen base state. For the case where the Rankine vortex and particle core have equal radii, i.e., = 1, a dispersion relation can be obtained analytically: Here, ( − 1) denotes the Dirac delta function. These functions correspond to perturbations that are concentrated at the number density and vorticity jumps = = 1. To further study the stabilizing/destabilizing effect of particle inertia, we expand the expression for growth rate (equation (3.39)) in a series of St Γ as ( 3.42) where Δ = (1 + At) 2 − 4 At. For St Γ = 0, Δ < 0 provides the condition for instability for a fixed , that is, a mode is stable if At < At cr = 2 − 1 − 2 √︁ ( − 1). For St Γ ≠ 0, particle inertia can play contrasting roles. For At < At cr , the dispersed phase destabilizes the vortex while for At > At cr , the dispersed phase stabilizes the vortex (3.44) Figure 5 (a) shows the variation of growth rate plotted over a range of Stokes numbers for wavenumber = 2. With the smaller value of Atwood number At = 0.1, the growth rates decrease with increasing Stokes numbers. However, the trends get reversed for higher Atwood numbers, with the growth rates increasing with increasing Stokes numbers, confirming the predictions from the asymptotic calculations (equations (3.43) and (3.44)). Figure 5 (b) plots the growth rates over a range of Atwood numbers for three values of Stokes numbers. It is found that the growth rate is always positive as long as the Atwood number At and Stokes number St Γ are non-zero. This indicates that destabilization of the vortex by the particles is guaranteed as long as the particle inertia is non-zero. For cases where > 1, we solve equations (3.26) and (3.27) numerically using Chebyshev spectral collocation method (Trefethen 2000). To handle the discontinuities, we smooth the base state using hyperbolic tangents of thickness Δ = 0.01. Figure 6 (a) and (b) show the contours of growth rate in At − plane for Stokes numbers 0.01 and 0.1 respectively. Unlike the case of St Γ = 0 (see figure 4), there is no region of stability indicating that mode = 2 is always unstable when the particles are inertial. Further, the particle-laden vortex is unstable even for values of larger than 1. Whereas the results for non-inertial particles are similar to those obtained by Dixit & Govindarajan (2011), inclusion of particle inertia breaks the equivalence with a radially startified Rankine vortex. For the case St Γ = 0.01, the optimal value of leading to the largest growth rate is slightly in excess of 1 for low values of At, and approaches = 1 as At increases. For the larger inertia case St Γ = 0.1, = 1 leads to the largest growth rates for all values of At.

Comparison of Linear Stability Analysis with Euler-Lagrange simulations
In order to validate the linear stability analysis, we compare the analytical results with those from Eulerian-Lagrangian simulations of the two-way coupled system. To make the comparison tractable and accessible, the case = 1 is used as the benchmark since the growth rates can be computed analytically from equation (3.38). Table 1 lists a summary of non-dimensional parameters for 8 cases presently considered. In all these simulations, the Stokes number, circulation Reynolds number Re Γ and ratio / are fixed at 0.002, 5000 and 1, respectively. The Stokes number is sufficiently small such that the initial state can be considered quasi-steady and comparisons with LSA are permissible. In dimensional terms, the carrier fluid has density = 1.2 kg/m 3 and viscosity = 1.8 × 10 −5 kg · m −1 · s −1 . The initial vortex core is = 0.385 m and circulation is Γ = 9.42 × 10 −2 m 2 · s −1 . Particles are seeded with density = 1000 kg/m 3 and average  volume fraction = 1.2 × 10 −3 . In all cases chosen, the ratio / is large (of the order 10 4 ), such that fluctuations due to the discrete particle forcing fall well below the viscous dissipation scale. In this way, the coupling between the fluid and particles occurs primarily through collective particle dynamics at scales comparable with , rather than discrete effects at scales comparable with the inter-particle distance.
The initial conditions of the numerical simulations correspond to a superposition of the base Rankine vortex and the perturbation eigenmodes (3.40) and (3.41). In order to capture the linear regime described by the linear stability analysis, the perturbations are initialized with a small amplitude = 0.03. The Dirac delta function that appears in the eigenfunctions (3.40) and (3.41) is discretized on the Cartesian grid according to the radial distance (i=1,2,3,4) between the vortex core and the four vertices of each cell. Cells that are intersected by the discontinuity (max( ) > 1 and min( ) < 1) have ( −1) = 1/(Δ Δ ), whereas ( −1) = 0 everywhere else. The Lagrangian particles are generated randomly within each cell based on the number density computed with the aforementioned discrete perturbation. The particle velocities are set to match the fluid velocity at their locations. Figures 7a and 7b show number density and axial vorticity isocontours of a representative perturbation with mode = 4 and disturbance amplitude = 0.03 .

18
S. Shuai,D. J. Dhas,A. Roy,and M. H. Kasbaoui Figure 8: Iso-contours of normalized particle volume fraction ( / ) for cases A, B, C, D, E (see table 1) at five instant time. As time progresses, the inertial particles are expelled out of the vortex cores, forming cluster arms of which the number equals to the mode number. The small perturbations continue growing and ultimately destroy the structure. Figure 8 shows the time evolution of the particle volume fraction between = 0 and = 6 for modes = 2−6. All time values in this section are non-dimensionalized with vortex response time . For 3.0, the perturbed modes grow linearly resulting in the deformation of the interface at the number density jump. The latter, which is initially circular, develops azimuthal perturbations depending on which mode has been triggered initially. Around 3.0, we notice the emergence of fast growing high wavenumber modes ∼ 24 − 26 that quickly overtake the initially seeded low wavenumber modes by = 4.5. These high wavenumber perturbations are triggered by the discrete Lagrangian particle forcing and errors associated with the discretization of the discontinuous eigenfunctions (3.40) and (3.41). As shown by the linear stability analysis, the perturbation growth rate increases with wavenumber , such that it may be anticipated that such modes will eventually dominate. Thus, whereas the earlystages are controlled by the seeded low-wavenumber modes, the later-stages are dominated by the high-wavenumber perturbations. By = 6, these modes evolve into long spiraling arms of inertial particles that expand radially outward, which indicates strong non-linearity.
Regions of unladen fluid in between the particle arms are sucked inward forming mushroomlike structures. The fact that the final stages seen in figure 8 show similar number density patterns irrespective of the seeded mode, suggests that the mode that would emerge naturally is a high wavenumber mode that results from the balance between two-way coupling and viscous dissipation. These dynamics may be related to the misalignment between the particle density interface and centripetal acceleration field (Joly et al. 2005), which is similar to the mechanism of wave interaction between density interface and vorticity interface (Dixit & Govindarajan 2011). The baroclinic torque caused by the misalignment concentrates vorticity on the spiral arms and destroys it in between, resulting in radial filaments.
In order to determine the total perturbation growth rates from the present simulations, we compute the perturbation kinetic energy, where is the base state velocity corresponding to a Rankine vortex. Figure 9a shows the evolution of total perturbation energy in logarithmic scale for simulations with perturbed modes = 2 -6. The perturbation energy is seen to increase with increasing mode number , reflecting the larger growth rate for higher mode numbers. The instability growth rate is obtained from Eulerian-Lagrangian simulations using the relation Figure 9b shows the time evolution of the growth rate given by (4.2). Unlike the linear stability analysis which predicts a constant value, the instability growth rate varies with time in Euler-Lagrange simulations. Figure 9b shows that the growth rate reaches peaks between ∼ 2 to ∼ 4, depending on the seeded mode. The early transient is due to imperfect initial conditions which causes the particles to be accelerated or decelerated until they reach their equilibrium velocity. The time window around the peak growth rate is associated with the exponential growth of kinetic energy seen in figure 9a, and is thus, closest to the dynamics predicted by the linear stability analysis. The growth rate in figure 9b eventually drops, as the perturbation kinetic energy in 9a saturates. The deviation from exponential growth correlates (a) (b) (c) Figure 12: Comparison of instability growth rates between linear stability analysis (LSA) and Eulerian-Lagrangian (EL) simulations at St Γ = 0.002 and Re Γ = 5000. (a) Growth rates for mode numbers =2 -6 at At = 1/3 ( = 1). (b) Growth rates for mode numbers =2 -6 at At = 1/5 ( = 1/2). (c) Growth rates for mode = 3 at various Atwood numbers. Symbols: u, LSA; 6, EL-peak growth rate of the total perturbation energy; @, EL-average growth rate of the total perturbation energy; A, EL-peak growth rate of the seeded mode.
with the emergence of non-linear effects in the flow and number density fields as seen in figure 8.
In order to quantify the contribution of a specific azimuthal mode to the overall dynamics, we compute the kinetic energy associated with this mode using energy is the sum of the energies of all Fourier modes, the mode initially seeded (i.e., = 2 for case A, = 3 for case B, and so on) accounts for most of the perturbation energy during the early evolution 4. Thus, the growth rate computed from the total perturbation energy, shown in figure 9b, matches the growth rate of the seeded mode early on. This is further evidenced in figure 11 showing the evolution of the kinetic energy and growth rate associated with the seeded mode. The latter peaks between = 2 and 3 depending on the wavenumber, at a level sensitively close to the peak growth rate of the total perturbation in figure 9b. At later times, the seeded mode is no longer dominant as the energy of other modes rises to a comparable level. Further, as time progresses, figure 10 shows that the energy of several lowwavenumber modes drops, which corresponds to negative growth rates, also seen in figure  11b. This indicates a non-linear energy transfer from low to higher wavenumber modes. Figure 12 shows comparisons between growth rates predicted by the linear stability analysis and those obtained from Euler-Lagrange simulations (equation (4.2)). From the latter, we report the total perturbation's peak growth rate and average value in a time window of size Δ = 1 centered on the peak (see figure 9b). We also report the peak growth rate of the seeded mode (see figure 11b). The difference between these values represents a confidence interval for the comparison of linear stability analysis and Euler-Lagrange simulations. Figure 12a and 12b show the results at various mode number =2 -6 for Atwood numbers At = 1/3 and 1/5, respectively. Both Euler-Lagrange simulations and linear stability analysis capture the growth rate increase with increasing mode number . The quantitative agreement is good given the different methodologies and underlying assumptions. The growth rates at various Atwood number for mode = 3 are shown in figure 12c. Here the Euler-Lagrange simulations and linear stability analysis reproduce the same trend, that is the instability growth rate increases with Atwood number At. Further, the peak growth rates obtained by the fully non-linear numerical simulation has excellent agreement with linear stability predictions.
While the growth rate obtained from numerical simulations and linear stability analysis are within close quantitative agreement, some of the differences can be related to the discretization of the eigenmodes and Lagrangian representation of the particles in the numerical simulations. Whereas eigenmodes derived from linear stability analysis require strict discontinuities in number density and vorticity, these effects are necessarily smoothed out when discretized on a mesh. Further results from the linear stability analysis are derived on the premise that viscous effects and discrete particle effects are negligible. However, these effects cannot be ruled out in the Euler-Lagrange simulations where viscosity still plays a role despite the large Re Γ = 5000 considered, and Lagrangian particles trigger high wavenumber modes. Despite the differences between the linear stability analysis and the Euler-Lagrange model, the agreement between the two approaches is satisfactory.

Conclusion
In this work, we have shown that introducing heavy particles in the core of a vortex causes the breakdown of the flow structure due to an instability activated by two-way coupling, i.e., momentum exchange between the two phases. If particle feedback is neglected (oneway coupling), no instability is observed: the vortex retains its coherent axisymmetric structure while inertial particles are slowly expelled outward forming a ring of clustered particles around the core. The inclusion of two-way coupling breaks the axisymmetry and gives rise to azimuthal perturbations in both number density and vorticity fields. As these perturbations develop, they turn into spiraling arms of concentrated particles emanating out of the core while regions of particle-free flow are sucked inward. The vorticity field displays similar pattern which cause the breakdown of the initial Rankine structure. Remarkably, this breakdown occurs even for inertia-less particles, i.e., zero-Stokes number, provided that the mass loading is not vanishingly small. Particle inertia plays a dual role, in that, it destabilizes low wavenumber modes, but reduces the growth rate of high wavenumber ones.
The mechanisms driving the instability are characterized using linear stability analysis of a base state consisting of a Rankine vortex with core radius and a circular patch of particles with radius . To describe the two-phase flow, we use an inviscid Two-Fluid model, and assume that the particles are weakly inertial, i.e., their Stokes number St Γ is small but not zero. This assumption has two merits. First, it allows us to approximate the particle velocity field in terms of the carrier fluid velocity field and an inertial correction, thus, reducing the number of equations to solve. Second, it justifies a quasi-steady approach in which the time-dependence of the base state is ignored enabling us to carry out a classic linear stability analysis. This step is justified by the fact that variations of the base state have a timescale = /St Γ which becomes very large in the limit St Γ 1 (Shuai & Kasbaoui 2022). Analysis of growth rates obtained from linear stability analysis reveals the existence of unstable modes regardless of the value of the Stokes number St Γ and as long as the mass loading ≠ 0. Generally, the most unstable configuration corresponds to matching particle patch and vortex core radii, i.e., = . For St Γ = 0 particles, the problem becomes analogous to a vortex with discontinuous radial density stratification described by Dixit & Govindarajan (2011). Modes with wavenumber above a critical wavenumber cr = (1 + ) 2 /((2 + ) ) are unstable and have a growth rate { } ∼ √ At when 1, and where the Atwood number At = /(2 + ) provides a measure of the relative magnitude of the density jump between inner (particle-rich) and outer (particle-free) regions (0 At 1). When the particles are inertial (St Γ ≠ 0) all modes become unstable regardless St Γ and ≠ 0, even those with < cr . However, the growth rate of modes cr is reduced compared to the case St Γ = 0. For large wavenumber modes 1, the growth rate approaches { } ∼ √ At(1 − 2St Γ √ At) + (St 2 Γ ). For both inertial and non-inertial particles, the tendency of growth rates to increase with wavenumber suggests that modes that naturally emerge in experiments and simulations shall have a high-wavenumber. This mode is likely determined by a balance between viscous effects and preferential concentration, although this cannot be ascertained as we have not considered viscous effects in our linear stability analysis.
In addition to the linear stability analysis, Eulerian-Lagrangian simulations are carried out at St Γ = 0.002, Re = 5000, and mass loading = 0.4 -1. The simulations are initialized with the base state superimposed with eigenmodes found from linear stability analysis. The simulations allow us to explore the transition from linear to nonlinear regimes. During the early evolution of the perturbations, growth rates computed from the perturbation kinetic energy in the Eulerian-Lagrangian simulations show excellent agreement with growth rates predicted by the linear stability analysis for modes = 2 -6 and mass loadings = 0.4 -1 (or At = 1/6 -1/3). When the flow enters a nonlinear stage, we observe the emergence of a high-wave number mode with ∼ 24 -26. This mode develops into spiraling arms of number density and vorticity that ultimately cause the breakdown of the initial base state.
One must note that the present study has a key difference from prior investigations of dusty Kolmogorov flow and Rayleigh-Taylor turbulence -a particle-free vortex is always stable. Our study emphasizes that the disperse phase can also be the source of a new instability besides modifying an existing one. The novel instability for a dusty vortex identified here highlights how the feedback force from a disperse phase can induce a breakdown of an otherwise resilient vortical structure.