2.3.1. Fixed points outside the ellipse

Outside the ellipse ($\tanh \xi > r$), (2.5) can be written by substituting the appropriate velocity expressions from (2.4) as

(2.6a)\begin{gather} (St \varOmega^2\sinh 2\xi+(k^{{-}2}\, {\rm e}^{{-}2\xi}-\varOmega)\sin 2\eta) h^2 = 0, \end{gather}

(2.6b)\begin{gather}(St \varOmega^2\sin 2\eta+\varOmega\sinh 2\xi+k^{{-}2}\, ({\rm e}^{{-}2\xi}\cos 2\eta-1)) h^2 = 0. \end{gather}

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

(2.7a)\begin{gather} p^{{\pm}} = \frac{2-(\alpha \pm \beta) k^2 \varOmega^2}{2(1+k^2 \varOmega)}, \end{gather}

(2.7b)\begin{gather}q^{{\pm}} = \frac{\alpha \mp \beta}{2\, St}, \end{gather}

where $\alpha = k^2 (1+St^2 \varOmega ^2)$ and $\beta = \sqrt {\alpha ^2-4 St^2}$. These solutions form a set of four fixed points outside the ellipse, located at $(\bar {\xi }_1^{+}, \bar {\eta }_1^{+}) =(\tanh ^{-1}p^+,\tan ^{-1}q^+)$, $(\bar {\xi }_1^{-}, \bar {\eta }_1^{-}) =(\tanh ^{-1}p^+,-{\rm \pi} +\tan ^{-1}q^+)$, $(\bar {\xi }_2^{+}, \bar {\eta }_2^{+}) = (\tanh ^{-1}p^-,\tan ^{-1}q^-)$ and $(\bar {\xi }_2^{-}, \bar {\eta }_2^{-}) =(\tanh ^{-1}$ $p^-,-{\rm \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 choose 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 3(a): A and B shift counterclockwise, while C and D shift clockwise as $St$ increases. As $St$ increases, the fixed points A($\bar {\xi }_1^{+}$, $\bar {\eta }_1^{+}$) and C($\bar {\xi }_2^{+}$, $\bar {\eta }_2^{+}$) approach each other. The same thing happens for the counterpart fixed points B($\bar {\xi }_1^{-}$, $\bar {\eta }_1^{-}$) and D($\bar {\xi }_2^{-}$, $\bar {\eta }_2^{-}$) as well.

Figure 3. (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.

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

(2.8)\begin{equation} \boldsymbol{v} = \boldsymbol{u}-St \left\{\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u} - \boldsymbol{x} \varOmega^2 + 2 \varOmega\hat{\boldsymbol{e}}_z \times \boldsymbol{u}\right\}+\textit{O}(St^2), \end{equation}

using which the particle trajectories are obtained and shown in the figure 3(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 than fluid tracers. They match exactly with the solution of (2.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 divergence of the system (2.2) in the four-dimensional phase space is given by $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}+\boldsymbol {\nabla }_{\boldsymbol {v}} \boldsymbol {\cdot } \dot {\boldsymbol {v}} = -2/St$, a negative quantity, indicating the dissipative nature of the dynamical system for any finite inertia particles. Here, $\boldsymbol {\nabla } \boldsymbol {\cdot } (\ )$ represents the divergence in physical space, while $\boldsymbol {\nabla }_{\boldsymbol {v}} \boldsymbol {\cdot } ()$ represents the divergence in velocity/momentum space. Even in the $St \ll 1$ limit, the divergence of the reduced system (2.8) in the 2-D phase space yields $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} =\textit {O}(St) \neq 0$, indicating a compressible nature of inertial particle flow in phase space unlike the fluid tracers where $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$. The dissipative nature of the system for heavy, finite inertia particles is responsible for the elliptic fixed points becoming spiral attractors. However, the hyperbolic fixed points remain intact, though their position gets 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 (2.3), evaluated at the fixed point ($\bar {\xi },\bar {\eta }$) is

(2.9)\begin{equation} \boldsymbol{\mathsf{J}} = \begin{pmatrix} 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ \varGamma & -\varDelta & -\dfrac{1}{St} & 2\, \varOmega\\ \varDelta-\dfrac{2 k^{{-}2}(1-p) (p+q^2) }{St (p^2+q^2)} & \varGamma+\dfrac{2 k^{{-}2} q (1-p)^2}{St (p^2+q^2)} & -2 \varOmega & -\dfrac{1}{St} \end{pmatrix}, \end{equation}

where $\varGamma =(p^2-1)\{\varOmega (1+q^2) (St \varOmega (p^2-q^2)-2pq)+2k^{-2}(p+q^2)q(1-p)\}/\{St (p^2+q^2)^2\}$ and $\varDelta = (1-p) (1+q^2) \{\varOmega (1+p)(p^2-q^2+2 p q\, St\, \varOmega )-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 $\boldsymbol{\mathsf{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 inertia particle. For instance, the variation of all the eigenvalues with Stokes number for each fixed point type is shown in figure 4. From figure 4(a,d), it can be seen that any finite $St$ particle will perceive fixed points A and B with purely real-valued eigenvalues with signature $(+,-,-,-)$, as saddles in the four-dimensional phase space. However, for the fixed point C or D, as shown in figure 4(b,e), the eigenvalue is 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 2). Note that, for the fixed point C or D, for $r=0.5$, when $St \gtrsim 0.771$ (shown in the insets of figure 4b,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 the quartic eigenvalue polynomial for $\boldsymbol{\mathsf{J}}$ (which can be obtained from (2.9) or (A4), by setting $\lvert \boldsymbol{\mathsf{J}}-\lambda \, \boldsymbol{\mathsf{I}} \rvert = 0$). By examining the discriminant of the solution for $\lambda$, one may find its behaviour changes when $St^2$ satisfies the cubic equation

(2.10)$$\begin{align} &St^6+\frac{3(1+r)^4(1+4r+r^2)}{2 r^3} St^4+\frac{(1+r)^8(1+6r+r^2)(5+22r+5 r^2)}{16 r^6}St^2\nonumber\\ &\qquad-\frac{(1+r)^{12}(1-r)^2(3+r)(1+3r)}{16 r^8}=0. \end{align}$$

The real-valued, non-negative solution of (2.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 figure 4(b,e).

Figure 4. 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.

As $St$ increases, saddle type and stable spiral type fixed points approach closer, as mentioned earlier. If we keep increasing $St$, we find that they merge (i.e. $p^+ = p^-$ and $q^+=q^-$) and vanish at a critical inertia value. From (2.7), one could deduce that this happens only if $\beta = 0$. Using the expression $\beta ^2 = k^4 (1+St^2 \varOmega ^2)^2-4 St^2$ and solving for Stokes number, one would find that 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 figure 4 for the Stokes number until $0.7721$.

2.3.2. Fixed points inside the ellipse

The origin O located inside the ellipse ($\tanh \xi < r$) continues to exist as a fixed point for heavy inertial particles (see figure 3b); however, O behaves as an unstable spiral. To show this, we follow the same procedure mentioned in the previous § 2.3.1, using the flow field inside the ellipse. By substituting appropriate expressions for $u_{\xi }$ and $u_{\eta }$ from (2.4) in (2.5), we obtain

(2.11a)$$\begin{gather} (2 St\, \varOmega\sinh 2\xi+k^{2}(\varLambda -\cosh 2\xi)\sin 2\eta) h^2 = 0, \end{gather}$$

(2.11b)$$\begin{gather}(2 St\, \varOmega\sin 2\eta-k^2\, (\varLambda -\cos 2\eta)\sinh 2\xi) h^2 = 0. \end{gather}$$

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$, (2.11) has a single real-valued solution $(\bar {\xi }_0,\bar {\eta }_0) = (0,{\rm \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

(2.12)\begin{equation} \boldsymbol{\mathsf{J}} = \frac{1}{St}\,\begin{pmatrix} 0 & 0 & St & 0\\ 0 & 0 & 0 & St\\ St\varOmega^2 & -r \varOmega & -1 & 2 St \varOmega\\ \varOmega/r & St \varOmega^2 & -2 St \varOmega & -1 \end{pmatrix}, \end{equation}

which has two pairs of complex conjugate eigenvalues. One of these pairs has a non-negative real part (see figure 4c,f), indicating unstable spiral behaviour for any non-zero $St$. Thus, the particles will be centrifuged away from the origin spirally. Some of them, starting within certain basins (see § 2.5) accumulate in the stable fixed points C and D outside the ellipse, as they are stable fixed points, as shown in the previous § 2.3.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.

The nature of fixed points and phase space topology is depicted schematically in figure 5 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 and D exhibit a ‘2-spiral sink’ behaviour 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 and B demonstrate a ‘$3:1$ saddle’ behaviour for $St < \mathcal {S}_2$. Additionally, fixed point A (B) annihilates with C (D) and ceases 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, Rieck & Sadlo (Reference Hofmann, Rieck and Sadlo2018).

Figure 5. 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.

3.1.1. Modification to the fixed points of the Kirchhoff vortex due to external straining

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

(3.4a)$$\begin{align} \ddot{\xi} &= \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 \nonumber\\ &\quad+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\}, \end{align}$$

(3.4b)$$\begin{align} \ddot{\eta} &= \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}}{2k\, St}\sin 2\eta\nonumber\\&\quad-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\}. \end{align}$$

Note that when $\dot {k} = \ddot {\theta } = 0$, these equations will reduce to the case of the Kirchhoff vortex as in (2.3); however, which is not the case for any finite strain $s \neq 0$. The flow velocity fields can be obtained from the total streamfunction ($\psi = \psi _v'+\psi _e'+({\dot {\theta }}/{2}) (x^2+y^2)$) as $u_\xi = ({h}/{k}) ({\partial \psi }/{\partial \eta })$ and $u_\eta = -({h}/{k}) ({\partial \psi }/{\partial \xi })$. One must remember that the ellipse's evolution (3.1) 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 (Reference IJzermans and Hagmeijer2006) as a periodic modification to the fixed points of the unperturbed system. Let ($\bar {\xi },\bar {\eta }$) denote the fixed points for any inertial particle in the Kirchhoff vortex, a time-dependent perturbation of the form $\xi = \bar {\xi }+s \xi '(t) + \textit {O}(s^2)$ and $\eta = \bar {\eta }+s \eta '(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 (streamfunction). For low strain rates, a regular perturbation approach yields the solution of (3.1) as $r = r_0+s r_1+\textit {O}(s^2)$ and $\theta = \theta _0+\varOmega _0 t+s \theta _1+\textit {O}(s^2)$, where

(3.5a)$$\begin{gather} r_1 ={-} \frac{r_0}{\varOmega_0}\sin(\varOmega_0 t)\sin(\varOmega_0 t+2 \theta_0) , \end{gather}$$

(3.5b)$$\begin{gather}\theta_1 = \frac{(1+r_0^3)}{ r_0 (1-r_0)}\sin(\varOmega_0 t)\cos(\varOmega_0 t+2 \theta_0)-\frac{t}{2}\, \frac{(1-r_0)}{(1+r_0)}\cos 2\theta_0 , \end{gather}$$

with $r_0 = r(t = 0)$, $\theta _0 = \theta (t = 0)$ and $\varOmega _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 \varLambda _0 r_1/(2 r_0)+\textit {O}(s^2)$, where $k_0^2 = (1/r_0-r_0)$ and $\varLambda _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 coordinates ($\xi = \bar {\xi }+s \xi '(t) + \textit {O}(s^2)$ and $\eta = \bar {\eta }+s \eta '(t) + \textit {O}(s^2)$), substituted in the governing equations (3.4) will reduce to the matrix equation

(3.6)\begin{equation} St\, \ddot{\chi}+\boldsymbol{\mathsf{L}}\,\dot{\chi}+\boldsymbol{\mathsf{K}}\,\chi ={-}\boldsymbol{\mathsf{M}}\, \lambda -N, \end{equation}

where $\chi = [\xi '(t),\eta '(t)]^\textrm {T}$, $\lambda = [\cos (2 \varOmega _0 t+2 \theta _0), \sin (2 \varOmega _0 t+2 \theta _0)]^\textrm {T}$, $N = [N_1, N_2]^\textrm {T}$ are vectors, and $\boldsymbol{\mathsf{K}}$, $\boldsymbol{\mathsf{L}}$ and $\boldsymbol{\mathsf{M}}$ are $2 \times 2$ matrices. The entries $N$, $\boldsymbol{\mathsf{K}}$, $\boldsymbol{\mathsf{L}}$ and $\boldsymbol{\mathsf{M}}$ are determined by the Kirchhoff vortex and its fixed points ($\bar {\xi },\bar {\eta }$) and are listed in Appendix B. Equation (3.6) 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 a periodic solution of the form $\chi = F_0 +\boldsymbol{\mathsf{F}}_1\, \lambda$ and the substitution yields two independent matrix equations (which are related to the particular solution of (3.6)),

(3.7a)$$\begin{gather} \boldsymbol{\mathsf{K}}\, F_0 ={-}N, \end{gather}$$

(3.7b)$$\begin{gather}\left(4\, \varOmega_0^2\, St\, \boldsymbol{\mathsf{I}}-2\, \varOmega_0\, \boldsymbol{\mathsf{L}}\, \boldsymbol{\varSigma}-\boldsymbol{\mathsf{K}} \right) \boldsymbol{\mathsf{F}}_1 = \boldsymbol{\mathsf{M}}, \end{gather}$$

where $\boldsymbol{\mathsf{I}} = \left (\begin{smallmatrix} 1 & 0\\ 0 & 1 \end{smallmatrix}\right )$ is a $2 \times 2$ identity matrix and $\boldsymbol {\varSigma } = \left (\begin{smallmatrix} 0 & -1\\ 1 & 0 \end{smallmatrix}\right )$. These equations are linear and a simple inversion yields the unknowns $F_0$ and $\boldsymbol{\mathsf{F}}_1$ and, thus, the limit cycles. However, the stability of the limit cycle is determined by the transient dynamics of the solutions of (3.6). It is determined by the eigenvalues of the matrix $\boldsymbol{\mathsf{P}}$ that satisfy the quadratic matrix equation $St\, \boldsymbol{\mathsf{P}}^2+\boldsymbol{\mathsf{L}}\boldsymbol {\cdot } \boldsymbol{\mathsf{P}}+\boldsymbol{\mathsf{K}} = \boldsymbol{\mathsf{0}}$ (related to the homogeneous solution of (3.6)). 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 are derived as a perturbation to the fixed points of the Kirchhoff vortex as ($[\xi,\eta ] = [\bar {\xi },\bar {\eta }]^T + s\, \chi +\textit {O}(s^2)$). They are then evaluated for the case of $r_0 = 0.5$ and $\theta _0 = 0$ and is shown in figure 12 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 11(d,h). While theoretically predicting that the hyperbolic fixed points (A and B) transition into unstable limit cycles under strain 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 § 3.1.2.

Figure 12. 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 (cf. figure 3a). 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$.

The perturbation to the fixed point at the origin O yields all the entries of $\boldsymbol{\mathsf{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 heavy inertial particles in the Kida vortex, which centrifuges them away from the centre.

In the inset of figure 11(d,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 11(h) that this method could not capture.

3.1.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 Reference Polvani and Wisdom1990). 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 Reference Smale1967; Ott Reference Ott2002), 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 Reference Bertozzi1987). 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 Reference Angilella2010; Angilella et al. Reference Angilella, Vilela and Motter2014). 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 (3.4) 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 the Kida vortex as

(3.8a)$$\begin{gather} \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), \end{gather}$$

(3.8b)$$\begin{gather}\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), \end{gather}$$

where the relevant functions $\hat {f}_1$ to $\hat {\phi }_2$ outside and inside the ellipse are listed in Appendix C. In (3.8) 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 that 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. (3.8) 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 Reference Bertozzi1987; Wiggins Reference Wiggins1990). For fluid tracers in the Kida vortex, Ngan et al. (Reference Ngan, Meacham and Morrison1996); Kawakami & Funakoshi (Reference Kawakami and Funakoshi1999) have evaluated the Melnikov function, which is shown to satisfy the requirements for tangle formation for any non-zero 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 the 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, will not 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 (3.8), the system can be rewritten up to $\textit {O}(s^2)$ as

(3.9a)$$\begin{gather} \dot{\xi} = \hat{f}_1(\xi,\eta)+s \left\{\hat{g}_1(\xi,\eta;t)+q \hat{\phi}_1(\xi,\eta)\right\}, \end{gather}$$

(3.9b)$$\begin{gather}\dot{\eta} = \hat{f}_2(\xi,\eta)+s \left\{\hat{g}_2(\xi,\eta;t)+q \hat{\phi}_2(\xi,\eta) \right\}. \end{gather}$$

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

(3.10)$$\begin{align} M_j^{{\pm}}(\phi) &= \int_{-\infty}^{\infty}\left[\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\}\right.\nonumber\\ &\quad\left.-\,\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\}\right]\,\nonumber\\ &\quad\times{\rm e}^{-\int_0^t \textrm{tr}(Df)\,{\rm d}t'}\, {\rm d}t, \end{align}$$

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) = {\partial \hat {f}_1}/{\partial \xi }+{\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, evaluated as $\boldsymbol {\nabla }_{\xi,\eta } \boldsymbol {\cdot } [u_{\xi },u_{\eta }] = h^2\, ( {\partial (\,\hat {f}_1/h^2)}/{\partial \xi }+{\partial (\,\hat {f}_2/h^2)}/{\partial \eta } )$, is different from the trace mentioned above; $\boldsymbol {\nabla }_{\xi,\eta } \boldsymbol {\cdot } [u_{\xi },u_{\eta }]=0$ in accordance with the Hamiltonian nature of tracer dynamics in the 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 ${\rm \pi} /\varOmega _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 +{\rm \pi} /\varOmega _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 Reference Angilella2010), and (ii) the transport of passive fluid tracers in the same flow field near a wall (see Angilella et al. Reference Angilella, Vilela and Motter2014). 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.

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 pure-strain 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 (3.10), 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 ${\rm \pi} /\varOmega _0$ is shown in figure 13. We can evaluate the critical value of $q$ at which the zeros of the Melnikov function disappear as $q_{{cr},1} \approx 0.0168$ and $q_{{cr},2} \approx 4.2395$, which is independent of whether it is a $`+'$ or $`-'$ type heteroclinic orbit (due to their symmetry). Thus, we may conclude that the heteroclinic tangles of $H_1^{\pm }$ disappear if $St > q_{{cr},1}\, s$ and that of $H_2^{\pm }$ disappear if $St > q_{{cr},2}\, s$. When $St > q_{{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_{{cr},2} = St/q_{{cr},2}$.

Figure 13. 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 14 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 14(a), the stable and unstable manifolds coincide without any transverse intersection. With a small strain perturbation, as shown in figure 14(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 14(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. Figure 14(d,e) examines the combined effect of strain perturbation and particle inertia. In figure 14(d), particles with weak inertia ($St < q_{{cr},1} s$) allow for weak centrifuging, enabling the transverse intersection of the folds to form tangles. Conversely, in figure 14(e), dominant particle inertia ($St > q_{{cr},1} s$) causes manifolds to spread apart and prevent the intersection of the folds.

Figure 14. 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 the Kirchhoff vortex ($s =0, St=0$), (b) for tracers in the Kida vortex ($s \ll 1, St=0$), (c) for inertial particles in the Kirchhoff vortex ($s =0, St \ll 1$), (d–e) inertial particles in the Kida vortex, where (d) $s \ll 1, St < q_{{cr},1} s$ and (e) $s \ll 1, St > q_{{cr},1} s$.

For inertial particles with $St = 0.1$, the critical strain rates $s_{{cr}} = St/q_{{cr}}$ can be obtained as $s_{{cr},1} \approx 5.9524$ and $s_{{cr},2} \approx 0.0236$, such that when $s \lessapprox 0.0236$, there will be no chaotic dynamics. In figure 15(a,b) we have shown the basin of attraction (evaluated using the full nonlinear system of (3.4)) of $St = 0.1$ particles in the 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. The dimension of the basins is evaluated using a method similar to that used by Angilella et al. (Reference Angilella, Vilela and Motter2014). We identify a set of points on a line intersecting the basin boundary. We chose the line segment $x = 1.68$ and $-0.4 \leq y \leq 0.4$, close to the hyperbolic fixed point A. Pairs of random points on the line segment are checked to see if they belong to different basins. If they do not, they are discarded. If they do, the bisection method is used to refine these points closer until their inter-distance reaches a threshold of $10^{-8}$, ensuring they approach the basin boundary. Once a few thousand such points are identified, lying on both the line and the basin boundary, we use a box-counting algorithm to identify the scaling law $N(\epsilon ) \sim \epsilon ^{-D^{(1)}}$ to obtain the dimension $D^{(1)}$ in the limit of $\epsilon \rightarrow 0$. Here, $N(\epsilon )$ is the number of boxes of side length $\epsilon$ required to cover the sampled points in the basin boundary. The dimension obtained ($D^{(1)}$) is the dimension of the intersection set between the basin boundary and the initial line segment. To obtain the basin boundary dimension in 2-D space, one needs to evaluate $D^{(2)} = 1 + D^{(1)}$. Figure 16(a) shows the variation of $\epsilon$ vs $N(\epsilon )$ for the cases of $s=0.02$ and $s=0.03$, corresponding to the fractal basins in figure 15. A curve fitting reveals that the respective slopes are zero and $-0.7679$, as mentioned in the figure, which indicates that the respective fractal dimensions are $D^{(2)} = 1$ for the basin in figure 15(a) and $D^{(2)} \approx 1.7679$ for the basin in figure 15(b). Figure 16(b) shows the variation of $D^{(2)}$ with $s$ for the basin boundaries corresponding to a fixed $St = 0.1$. The critical value of $s$ for the onset of a fractal basin boundary is $s_{{crit}}\approx 0.0288$. The prediction of $s_{{crit}}$ from the Melnikov analysis of the overdamped system (3.9) is $s_{{cr,2}} \approx 0.0236$. We would expect the predictions of the two methods to agree better for smaller values of $St$. The overdamped system is strictly valid in the asymptotic limit of $St\rightarrow 0$; $St=0.1$ is significantly large from the perspective of its validity (see Nath et al. (Reference Nath, Roy, Govindarajan and Ravichandran2022) for discussion on flows with stagnation points). However, the overdamped system and the subsequent Melnikov analysis offer significant computational advantages over the full system, which tends to get very stiff for small $St$. It can be used to efficiently diagnose the chaotic dynamics when we are interested in first-order effects of particle inertia.

Figure 15. 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) $s = 0.02$ indicating smooth basin boundaries and (b) $s = 0.03$ indicating fractal basin boundaries. Figures are generated using the full system of (3.4) along with (3.1).

Figure 16. (a) Plot showing the variation of $\log (N(\epsilon ))$ with $\log (\epsilon )$, where $N(\epsilon )$ is the number of boxes of side length $\epsilon$ required to cover the set of points in the basin boundary that intersects with a line segment at $x = 1.68$ and $-0.4 \leq y \leq 0.4$. The cases shown correspond to the two typical strain rates in figure 15. The black straight lines represent the fitted curves with the specific slopes mentioned. (b) The variation of the dimension $D^{(2)}$ with the strain rate $s$. All cases here correspond to $St = 0.1$ particles in a Kida vortex with $r_0 = 0.5$ and $\theta _0 = 0$.

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 an arbitrarily high period and (ii) an uncountable infinity of periodic orbits (Smale Reference Smale1967). 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 a Kida vortex (Polvani & Wisdom Reference Polvani and Wisdom1990; Dahleh Reference Dahleh1992; Ngan et al. Reference Ngan, Meacham and Morrison1996; Kawakami & Funakoshi Reference Kawakami and Funakoshi1999), one can assume the same continues to exist for inertial particles until they have the critical inertia value. For further validation, we have evaluated the largest Lyapunov exponent (LLE) of the reduced dynamical system (3.9) to show the existence of chaos for weakly inertial particles and its disappearance for heavy inertial particles. We used the classical algorithm by Benettin et al. (Reference Benettin, Galgani, Giorgilli and Strelcyn1980) for this purpose. For a fixed strain rate of $s = 10^{-2}$, we evaluated the time variation of the LLE for various $St$ values for a point close to the hyperbolic fixed point A of the Kida vortex with $r_0 = 0.5$ and $\theta _0 = 0$. The results are shown in figure 17. Note that fixed point A periodically changes location as the system changes periodically with time. Thus, we chose the location of point A at $t = 0$ to evaluate the Lyapunov exponent. The figure shows that the system has a non-zero positive asymptotic value for the LLE for small values of $q$ (or equivalently $St$, as $St = q \, s$). As $q$ increases beyond the first critical value $q_{{cr},1} \approx 0.0168$, the variation of the LLE with time does not seem to saturate for large times, which indicates a weakening of chaotic behaviour. Similarly, as $q$ further increases to the next critical limit $q_{{cr},2} \approx 4.2395$, the LLE asymptotically approaches zero at large times, indicating a complete disappearance of chaos.

Figure 17. The variation of the largest Lyapunov exponent (LLE) with time for a Kida vortex with $r_0 = 0.5$, $\theta _0 = 0$ and $s = 10^{-2}$, corresponding to inertial particle trajectories of different $q$ (or $St$). Note that the case $q = 0$ indicates the LLE associated with tracers in the Kida vortex. The initial location of the particle is chosen to be the location of hyperbolic fixed point A at $t = 0$.

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 18. The initial orientation ($\theta _0$) of the elliptic vortex will only phase shift the Melnikov functions; thus, these critical values will not 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.

Figure 18. (a) The critical curves that 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$ vs $r_0$, for a Kida vortex, obtained from the Melnikov analysis. The initial orientation of the 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.

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_{{cr},2} \sim 1$. The caveat here is that the Kirchhoff vortex's saddles disappear for $St>\mathcal {S}_2$; on the other hand, as mentioned earlier, the Kida vortex itself is elongated and destroyed beyond critical strain rate $0.3$. 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 18(a).

3.1.3. 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 19 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 $\textrm {A}_p$ and $\textrm {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 (Ngan et al. Reference Ngan, Meacham and Morrison1996; Kawakami & Funakoshi Reference Kawakami and Funakoshi1999). As done in the previous § 3.1.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.

Figure 19. 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.

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 18(b). For $r_0=0.5$ and $St = 0.1$, the critical value of the strain rate is found to be $s_{{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 18(b) covers even smaller values of $St$.