2.1. The derivation of classic LPM

For the purpose of clarity in thesubsequent discussion, the basic idea of the LPM is first introduced for the 2-D case. Specifically, two incompressible inviscid fluids are considered in a vertically infinite channel with the left and right boundaries satisfying the no-penetration condition. The initial material interface is characterised by a single-mode sinusoidal perturbation. The LPM for the RTI problem can be described by the following governing equations (Goncharov Reference Goncharov2002; Betti & Sanz Reference Betti and Sanz2006):

(2.1) \begin{align} {{\boldsymbol{\nabla }} \boldsymbol{\cdot }\boldsymbol{{v}} _{i} = 0}, \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} {\partial _{t}\eta + \boldsymbol{{v}}_{i}\boldsymbol{\cdot }{\boldsymbol{\nabla }}\eta -\boldsymbol{{v}} _{i}\boldsymbol{\cdot }\boldsymbol{{e}} _{z}=0\quad \text{at } z=\eta }, \\[-28pt] \nonumber \end{align}

(2.3) \begin{align} {\rho _i(\partial _{t}\boldsymbol{{v}}_{i}+\boldsymbol{{v}} _{i}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{{v}} _{i}+g\boldsymbol{{e}} _{z}) = -{\boldsymbol{\nabla }}p_{i}}, \\[-28pt] \nonumber \end{align}

(2.4) \begin{align} {p_{1}=p_{2}\quad\text{at } z=\eta }. \\[6pt] \nonumber \end{align}

Here, $x$ and $z$ are the coordinates perpendicular and parallel to the acceleration $g$ , and $\boldsymbol{{e}} _{x}$ and $\boldsymbol{{e}} _{z}$ are the unit vectors along the $x$ and $z$ directions, respectively. The subscripts $i=1,2$ denote heavy and light fluids, respectively. The velocity vector $\boldsymbol{{v}} _{i}$ is related to its potential function $\varphi _i$ by ${\boldsymbol{{v}} _{i}}\equiv \boldsymbol{\nabla }\varphi _i$ (Pimentel-Garcia Reference Pimentel-Garcia2023), and $p_{i}$ and $\rho _i$ represent the corresponding pressure and density, respectively. The interface between the two fluids is located at $z=\eta ( x,t )$ . According to the LPM, the bubble is approximated as a parabola of the form $\eta ( x,t )=z_0(t) + \xi (t)x^2$ , where $z_0(t)$ and $\xi (t)$ are the vertical position and curvature of the bubble, respectively. For the ideal incompressible fluids consider in this study, both the kinematic and dynamic conditions are imposed at the interface of the two fluids. The kinematic condition at the interface ensures that fluid particles cannot penetrate the interface, i.e. the normal velocities of the fluid particles on both sides of the interface are equal to each other. This can be represented by the kinematic equation of the interface (2.2), which is derived by establishing the material derivative $\mathrm{D}F(x,z,t)/\mathrm{D}t=0$ with $F(x,z,t)=z-\eta (x,t)=0$ since the two fluids share the same material interface. Applying (2.2) to both fluids can readily yield a zero jump in normal velocity across the interface (denoted by $[\boldsymbol{{v}} \boldsymbol{\cdot }\boldsymbol{{n}}]=0$ , with $\boldsymbol{{n}}$ being the unit normal vector at the interface). For the dynamic condition at the interface, the continuity of pressure is ensured across the interface (denoted by $[p]=0$ or (2.4)). In addition, the fluid is stationary and undisturbed at the far-field boundaries along the flow direction ( $z$ ), i.e. $\boldsymbol{{v}}_i{(z\to \pm \infty )}=\boldsymbol{\nabla }\varphi _i{(z\to \pm \infty )}=0$ . In the spanwise direction ( $x$ ), only the velocity component along the flow direction is allowed at the boundary since the evolution of a complete single-mode interface is considered in the present study. Therefore, the spanwise velocity boundary condition can be described as a no-penetration condition, i.e. $v_{i,x}{(x=\pm ({\lambda }/{2}))}=\partial _{x}\varphi _i{(x=\pm ({\lambda }/{2}))}=0$ . Thus, given a reasonable initial spatial distribution of the flow field, the bubble evolution can be obtained from the interface kinematic condition (2.2) and the momentum equation(2.3) under the pressure continuity condition (2.4) at the interface.

2.2. Analysis of the classic LPM framework

Within the framework of the classical LPM, the sinusoidal potential function is employed to approximate the flow dynamics near the fingertip and in its immediate vicinity. While this approach approximately satisfies the dynamic condition of constant surface pressure, it assumes a time-invariant fingertip shape when predicting the growth rate. This assumption, however, conflicts with numerical simulations reported by Sohn (Reference Sohn2004) and Ramaprabhu & Dimonte (Reference Ramaprabhu and Dimonte2005), leading to unphysical predictions of the growth rate (Krechetnikov Reference Krechetnikov2009). Furthermore, the theoretical validity of potential functions without harmonics or singularities remains questionable as pointed out by Kull (Reference Kull1983). Additionally, the gravitational acceleration predicted by the sinusoidal potential-flow model exhibits a rapid increase below the fingertip, with the trend intensified for larger initial amplitudes (Kull Reference Kull1986). This result highlights a significant discrepancy considering the fact that the predicted acceleration is non-uniform across the flow field, whereas a conservative force such as gravitational acceleration should remain uniform. Consequently, the LPM struggles to capture the global flow characteristics such as the intense fluid motion localised near the interface and the negligible motion far from the interface.

Despite these limitations, the classical LPM performs reasonably well in describing the growth rates, which depend only weakly on acceleration uniformity. However, it fails to accurately predict the curvature of the fingertip, a more sensitive parameter that varies with the Atwood number $A$ . According to Chan et al. (Reference Chan, Jain, Hwang, Naveh and Abarzhi2023), the LPM predicts a less curved bubble tip but a more curved spike compared with numerical simulations. Although some experiments and numerical simulations appear to agree with LPM predictions, such agreement might be superficial as these investigations usually feature low spatial and temporal resolutions, limited observation time scales (approximately one decade) and a restricted set of diagnostic parameters.

This paper aims to provide an objective analysis of the RTI during the RA and SVS stages, emphasising the intrinsic physics of these processes. To this end, a dual-source model is proposed by introducing a dual-source velocity potential function in a uniform incoming flow. This approach accurately captures the streamline distribution and significantly improves the predictive performance of the potential-flow model, especially under low-density ratio conditions, as demonstrated by Liu et al. (Reference Liu, Wu-Wang, Zhang and Xiao2023a ). However, it should be mentioned that our approach has its limitations, which include the lack of rigorous justification for dual-source hypothesis, the approximate nature of the solution along a ‘line’ of vertices and the primary focus on local rather than global flow characteristics. Nevertheless, as a robust theoretical framework that has been widely employed by researchers over the past decades, the LPM provides a valuable tool for the present study.

2.3. The LPM for the later stages

For convenience in the following discussion, we show in figures 2(a) and 2(b) the flow fields at two extremely important and distinctive stages, i.e. the QS and RA stages, of RTI evolution with an Atwood number of $A=0.3$ , which are obtained in anumerical simulation using an in-house code CFD $^{2}$ (Zhang et al. Reference Zhang, He, Xie, Xiao and Tian2020). The Atwood number is defined as the relative density ratio in the form $A\equiv ( \rho _{1} - \rho _{2} )/ ( \rho _{1} + \rho _{2} )\in [ 0,1 ]$ , which is frequently used for a unified description of the RTI flows (Zhang & Guo Reference Zhang and Guo2016). Detailed introductions of the numerical methods and input parameters are provided in § 3.

Figure 2.

Vorticity distributions of the RTI flow field at the QS stage (a) and SVS stage (b) obtained from the present numerical simulation ( $A=0.3$ ), which are overlaid on contours of mass fraction. The finger is generated by the penetration of heavy fluid (red) into light fluid (blue). The black solid lines visualise the contours of vorticity. The black dashed line indicates the initial interface.

In most of previous investigations, the focus is placed on the bubble evolution up to the QS stage. The physical image at the QS stage (see figure 2 a) shows that the vorticity generated at the interface remains localised in its neighbourhood, which implies that the vorticity-induced velocity at the bubble front is negligible (Layzer Reference Layzer1955). For this reason, the classical LPMs consider an irrotational flow system to further simplify the governing equations (2.1)–(2.4) by introducing the velocity potential function $\varphi$ . However, as shown in figure 2(b), the vorticity has already accumulated around the bubble when the RTI evolution enters the RA stage. At this point, the frontal vorticity-induced velocity cannot be ignored, which breaks the standard equilibrium (Betti & Sanz Reference Betti and Sanz2006) between buoyancy and drag forces, and the bubble velocity no longer follows the prediction of the classical LPM. As discussed above, there is no available analytical theory that can properly predict the bubble evolution beyond the QS stage to date. To this end, a novel theory is proposed in this paper by considering a bilaterally rotational flow system.

The velocity vector in the RTI can be decomposed into two parts according to the Helmholtz theorem:

(2.5) \begin{equation} \boldsymbol{v}_i = \boldsymbol{\nabla }\varphi _i + \boldsymbol{\nabla }\times \mathbf{\psi }_i, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\mathbf{\psi }_i = 0. \end{equation}

Here $\varphi _i$ is the scalar potential and $\mathbf{\psi }_i$ is the vector potential. The vorticity field is then given by

(2.6) \begin{equation} \boldsymbol{{\omega }_i} = \boldsymbol{\nabla }\times \boldsymbol{v}_i = \boldsymbol{\nabla }\times (\boldsymbol{\nabla }\times \mathbf{\psi }_i). \end{equation}

For incompressible flow, (2.6) is further simplified to

(2.7) \begin{equation} \boldsymbol{{\omega }_i} = -{\nabla} ^2 \mathbf{\psi }_i. \end{equation}

For a given vorticity distribution, one can solve for the vector potential from (2.7) with prescribed boundary conditions. Alternatively, one can obtain the vorticity distribution readily from (2.7) for a constructed vector potential. For the RTI, the vorticity effect can be ignored before the RA stage, and the velocity field can automatically degenerate into the solution of a classical LPM. After the QS stage, however, the modulation of vorticity-induced velocity becomes prominent with the increase in vorticity accumulation effect. Moreover, although the bubble undergoes a complex nonlinear evolution, the pressure at the interface always satisfies the continuity condition to avoid breakup of the interface, which essentially imposes an important constraint for the flow around the interface. Based on these considerations, the present modelling approach comprises two key steps. A coherent velocity field is first constructed by incorporating features of both potential-flow and rotational flow regimes on both sides of the interface, and then the relationship between the velocity fields on either side is determined through the pressure continuity condition at the interface.

A qualitative analysis of figure 2 can be instructive to properly construct the model of the flow field. As is clearly seen, the development of potential velocity (in the bulk) and vortex structures is indeed confined to the vicinity of the interface, and their intensities are localised around the bubble vertex. This is because the dynamics of the RTI in the early stage are inherently concentrated at the interface and its vicinity, with their effects rapidly diminishing at locations far way from the interface. In addition, the vorticities are symmetrically distributed around the vertex and cancel out each other at the bubble vertex. In this sense, the derivation of the present model should be based on these physical pictures.

For a rotational flow system, it is assumed that the velocity field is given by the superposition of a potential flow and a prescribed rotational flow. Formally, the potential flow can be described by the general velocity potential function that satisfies the Laplace equation $\boldsymbol{\nabla} ^2 \varphi _{i}=0$ with the corresponding boundary conditions as aforementioned . In the LPM, the 2-D and 3D potential functions $\varphi _{i}$ are written as (Liu et al. Reference Liu, Zhang and Xiao2023b )

(2.8) \begin{equation} \varphi _i = Q_i(t) \sum _{n=0}^\infty \frac {\cos [(2n+1)kx]}{(2n+1)} {\mathrm{e}}^{(-1)^i (2n+1)k(z - z_0 + a_i(t))} \quad \text{(2D)}, \end{equation}

and

(2.9) \begin{equation} \varphi _i = Q_i(t) \sum _{n=0}^\infty \frac {\cos \big[(2n+1)k_{3D}^x x \big]\cos \big[(2n+1)k_{3D}^y y \big]}{(2n+1)}{\mathrm{e}}^{(-1)^i (2n+1)k_{3D}(z - z_0 + a_i(t))} \quad \text{(3D)}, \end{equation}

respectively, where $k\equiv 2\pi /\lambda$ is the wavenumber. However, it poses big challenge in practical computations if one retains the complete form of the Fourier series. The series truncation in LPMs has been a long-standing challenge in modelling nonlinear RTI dynamics. It is noted that Abarzhi (Reference Abarzhi1995, Reference Abarzhi1996a , Reference Abarzhib ) is among the first to systematically analyse the truncation effects based on the group theory approach, demonstrating that the optimal truncation order depends on the interface symmetry and higher-order harmonics must satisfy conservation laws. It is also shown that the truncation error scales as $\mathcal{O}(n^{-3/2}$ ) for $n$ retained modes. These findings help fundamentally shape the subsequent developments in potential-flow models, including the present work. While different mathematical tools are employed in the present LPM-based approach, the use of a physical principle for mode selection can be of particular importance based on energy conservation within the group theory framework. Besides, the series truncation in LPMs has been further addressed via group-theoretic principles (Abarzhi Reference Abarzhi1998, Reference Abarzhi1999). We retain terms up to $n=2$ in the Fourier series for $\varphi _i$ and $\omega _i$ , ensuring a third-order spatial accuracy. For $n=2$ , the relative error in $v_{\textit{RT}}$ is reduced to less than $3\,\%$ , validating the convergence of the truncated system. Considering that the LPM method essentially provides an approximate solution along a ‘line’ of vertices, Goncharov (Reference Goncharov2002) demonstrates that higher-order series expansions can improve model accuracy without significantly altering the predictive trends. Specifically, retaining only the first two terms in each series can provide a solution with up to $95\,\%$ accuracy. This indicates that including only the first two orders of the Fourier series is sufficient to capture the main features of large-scale coherent motion in the RTI. Therefore, to simplify the calculation, we take the first term $n=0$ for $\varphi _{1}$ and the first two terms $n=0,1$ for $\varphi _{2}$ to describe the complex flow structures inside the bubble. As for the rotational flow part, the vorticity is related to velocity by ${\boldsymbol{\nabla }} \times \boldsymbol{{v}} _{i} = \boldsymbol{{\omega }} _{i}$ . Here ${\omega _{i}}$ is the vorticities for species $i$ (the corresponding flow patterns can be seen in figures 2 a and 2 b), which vanishes at the peaks and troughs of the interface and is strongest roughly at a middle position between them, but lightly skewed to the spike side.

Therefore, it is assumed that the interfacial vorticity distribution should be approximated as a sinusoidal function, with its extreme points located near the interfacial centres and its zero point located at the finger front of bubbles and spikes. According to the numerical results (Glimm, Li & Lin Reference Glimm, Li and Lin2002; Sohn Reference Sohn2011), the vorticity is stronger at the spike side than at the bubble side. In this study, the sinusoidal function also considers the variation of vorticity in the $z$ direction to make vorticity generation more reasonable near the spike side. In addition, it is noted in figure 2 that the evolution of vorticity is concentrated near the interface and decays rapidly to zero far away from the interface. This behaviour is consistent with that of the velocity potential function field. It indicates that the dynamic characteristics of the RTI are primarily localised at the interface, with intense fluid motion occurring in its vicinity and successively weak motions farther away from the interface. Based on this phenomenological analysis, it is conjectured that the vorticity takes the simple form $\boldsymbol{{\omega _{i}}} = - \omega _{i}(t) \sin (kx)\exp[(-1)^ik(z-z_0)] \boldsymbol{{e}} _{y}$ .

The vorticity ansatz $\omega _i = -\omega _i(t)\sin (kx)e^{(-1)^ik(z-z_0)}$ in our model ensures that vorticity is confined to the interface region and decays exponentially with distance from the interface. The characteristic decay length $L_d \sim k^{-1}$ is determined solely by the initial perturbation wavenumber $k$ , consistent with rigorous theoretical analyses of an inviscid RTI (Abarzhi Reference Abarzhi2008; Williams & Abarzhi Reference Williams and Abarzhi2022). This aligns with the group-theoretic framework, where the interfacial vorticity structures are localised within a distance proportional to the initial wavelength $\lambda = 2\pi /k$ , independent of viscous effects.

The exponential decay profile ( $\sim{\mathrm{e}}^{-k|z|}$ ) of the proposed vorticity aligns with the asymptotic analytical solutions from Zhou et al. (Reference Zhou, Clark, Clark, Gail Glendinning, Aaron Skinner, Huntington, Hurricane, Dimits and Remington2019), who derived interfacial vorticity confinement in an inviscid RTI. This also confirms that the vorticity decay length scale is governed solely by the initial conditions, not viscosity.

However, due to the complexity in solving the vorticity-induced velocity and the fact that the proposed model focuses on the flow in the vicinity of the bubble vertex, the expressions for the velocity field are only given at the vertex $(z=\eta )$ . By combing both flow regimes, the velocity at the bubble front $(z=\eta )$ can be written as

(2.10) \begin{align} \boldsymbol{{v}} _{i} = {\boldsymbol{\nabla }} \varphi _{i} + \frac {3\omega _{i}(t) \cos (kx)}{2k}\boldsymbol{{e}} _{z}. \end{align}

Given that the planar RTI evolves only along the flow direction, the 3-D model can be extended from the 2-D model by employing a centrally symmetric approach (Liu et al. Reference Liu, Wu-Wang, Zhang and Xiao2023a , Reference Liu, Zhang and Xiaob ). Therefore, for 3-D cases such as the RTI flow in a square box (Yan et al. Reference Yan, Betti, Sanz, Aluie, Liu and Frank2016; Guo & Zhang Reference Guo and Zhang2020), the symmetry of this specific configuration leads to a potential function $\varphi _{i}$ in a series form $\varphi _{i}=Q_{i}(t)\sum _{n=0}^{\infty }\{\cos [(2n+1)k_{3D}^{x}]\cos [(2n+1)k_{3D}^{y}]\exp$ $[(-1)^i(2n+1)k(z-z_0+a_{i}(t)) ]/(2n+1) \}$ , where, $k_{3D}^{x}=$ $k_{3D}^{y}=k/\sqrt {2}$ denotes the wave numbers in the two horizontal directions. Meanwhile, the corresponding vorticity can be written as $\boldsymbol{{\omega _{i}}} =- \omega _{i}(t)exp[(-1)^ik/\sqrt {2}(z-z_0)]$ $\sin (k_{3D}^{x})\sin (k_{3D}^{y}) (\boldsymbol{{e}} _{x} + \boldsymbol{{e}} _{y})$ . Other settings are the same as those in the 2-D case. Consequently, the velocities at the bubble front in the 3-D case read

(2.11) \begin{equation} \begin{aligned} \boldsymbol{{v}} _{i} = {\boldsymbol{\nabla }} \varphi _{i} + \frac {9\omega _{i}(t) \cos (k_{3D})}{2k_{3D}}\boldsymbol{{e}} _{z}. \end{aligned} \end{equation}

Figure 3.

Evolutions of (a) the bubble vorticities on both the light- ( $\omega _2$ ) and heavy-fluid ( $\omega _1$ ) sides and (b) the vorticity ratio around the bubble tip obtained in present numerical simulations of a 2-D RTI at $A=0.3$ using CFD $^{2}$ .

Considering the continuity of interfacial pressure, one can derive the constraint relation governing the velocity fields on both sides of the interface. As mentioned previously, the velocity field can be decomposed into two distinct components, i.e. the velocity induced by potential forces and that induced by vorticity. While the velocity component induced by constant gravitational acceleration $g$ has been extensively studied within the classical LPM, our focus is placed on analysis of the constraint imposed on the vorticity-induced velocity component. The vorticities are calculated around the bubble tip on both the light- and heavy-fluid sides, which are shown in figure 3(a) with respect to time. It is observed that the vorticity ( $\omega _2$ ) for the light fluid is much larger in the late stage than in the early stage. For the heavy fluid, however, the vorticity ( $\omega _1$ ) remains relatively small throughout the entire range of time. In fact, it is seen that $\omega _1$ and $\omega _2$ are comparable with each other during the early stage of the RTI evolution. In other words, the influence of vorticity around the bubble tip can be negligible in the early stage of the RTI development. The ratio $r(A) = \omega _1/\omega _2$ is derived up to third order by solving the expanded pressure continuity condition. The revised expression

(2.12) \begin{equation} r(A) = 0.2(1 - A) + 0.05(1 - A)^2 + \mathcal{O}((1 - A)^3) \end{equation}

improves agreement with numerical data (figure 3 b, dashed curve). Stability analysis confirms that $r(A)$ remains invariant under perturbations to the initial vorticity distribution. As one can see, the vorticities on both sides of the bubble tip are close to each other in the early stage (on the left-hand side of the vertical dashed line), especially in the QS stage where there exists a local plateau with a value of about 0.8. After the QS stage, however, the curve drops drastically and tends to converge to a constant value of 0.16 during the SVS stage. Although the vorticity ratio is time dependent from the early to late stages, a constant assumption during the SVS stage is reasonable for a given Atwood number. In light of the above numerical analysis, it is conjectured that the ratio of bubble-tip vorticites on the light- and heavy-fluid sides can be written as

(2.13) \begin{equation} r(A)=\frac {\omega _1(t)}{\omega _2(t)}, \end{equation}

which varies with $A$ . Here, we assume that the ratio $r(A)$ is a constant value throughout the entire stage, which hardly influences the prediction accuracy of the model in the early stage RTI evolution because the vorticity effect around the bubble tip is negligible during these stages. According to relation (2.13), the vorticity on the light-fluid side can be expressed in the form $\omega _2=r(A)\omega _{1}(t)$ .

In summary, there are three unknowns to be determined in the present model for prediction of the RTI interface evolution, i.e. $z_0$ , $\xi$ and $\omega _{1}$ , provided that the vorticity ratio is given a priori or empirically. Following the LPM, the present model only solves the RTI problem in the neighbourhood of the fingertip instead of the whole interface. Thus, the series expansions in the $x$ direction (as well as the potential functions and the interfacial vorticity functions) are valid in the neighbourhood of the bubble vertex. As a result, the auxiliary parameters $Q_{1}$ , $Q_{2}$ , $a_{1}$ and $a_{2}$ can be solved by substituting expansions of $\varphi _{1}$ and $\varphi _{2}$ into (2.2) and expanding the resultant equations to the order of $x^2$ . The amplitude of potential functions and interfacial vorticity functions can be obtained from (2.3) of the bubble vertex. Therefore, in order to close the current model, one needs to further integrate (2.3) and expand through orders of $x^2$ and $x^4$ in (2.4), respectively, to yield the following ODEs:

(2.14) \begin{align} & \dot \xi = F(k, r, \xi , v, \omega _{1}), \\[-12pt]\nonumber \end{align}

(2.15) \begin{align} & {S_1}g + {S_2}v^2 + {S_3}\dot v + {S_4}v\omega _{1} + {S_5}\omega _{1}^2 + {S_6} \dot \omega _{1} = 0, \\[-12pt]\nonumber \end{align}

(2.16) \begin{align} & {T_2}v^2 + {T_3}\dot v + {T_4}v\omega _{1} + {T_5}\omega _{1}^2 + {T_6} \dot \omega _{1} = 0. \\[7pt]\nonumber \end{align}

Here, $v={\mathrm{d}}z_0/{\mathrm{d}}t$ is the bubble velocity; $F(k, r, \xi , v, \omega _{1})$ can be obtained from (2.2) in both 2-D and 3-D cases; $S_{p}(A, k, r, \xi , v, \omega _{1})$ and $T_{p}(A, k, r, \xi , v, \omega _{1})$ $(p=1,\boldsymbol{\cdot}s ,6)$ are coefficients of the $x^2$ and $x^4$ terms in (2.3) and (2.4). The detailed expressions of these coefficients are obtained with the aid of Mathematica 11.2 (Wolfram Research, Inc., 2017) and provided in Appendix A for reference. It should be pointed out that the LPM focuses on the time evolution of the bubble and physical quantities in its neighbourhood. Hence, the interfacial vorticity in the LPM corresponds to the neighbouring vorticity around the bubble, which reflects the accumulation of vorticity therein.

By solving the ODEs (2.14)–(2.16), it is found that there exist two fixed points with $\dot v=0$ , $\dot \xi =0$ and $\dot \omega _{1}=0$ . Accordingly, the bubble velocity, curvature and the neighbouring vorticity converge to constant values. For the first fixed point, the neighbouring vorticity around the bubble approaches zero (i.e. $\omega _{1} \sim 0$ ) and can be ignored. This corresponds to the classical asymptotic behaviour of the RTI in the QS stage. As for the second fixed point, the neighbouring vorticity converges to a non-zero constant, which is significantly larger than that for the first fixed point. The second plateau corresponds to a new asymptotic behaviour that is rarely reported in literature.

Linear stability analysis of the ODEs (2.14)–(2.16) reveals that the QS and SVS fixed points are asymptotically stable. The Jacobian matrix $\mathcal{J}$ at these points has eigenvalues with negative real parts (i.e. $\mathrm{Re}(\lambda _1, \lambda _2, \lambda _3) = (-0.12, -0.08, -0.05)$ ) for $A=0.3$ . This ensures that small perturbations decay exponentially, corroborating the numerical observations and confirming an asymptotic stability.

To explore the newly identified characteristic behaviour, comparisons are made between existing numerical simulations, which successfully capture the RA stage of a single-mode RTI (Glimm et al. Reference Glimm, Li and Lin2002; Ramaprabhu et al. Reference Ramaprabhu, Dimonte, Woodward, Fryer, Rockefeller, Muthuraman, Lin and Jayaraj2012; Wei & Livescu Reference Wei and Livescu2012; Bian et al. Reference Bian, Aluie, Zhao, Zhang and Livescu2020). The collective findings show that, at high Reynolds numbers, the bubble growth rate momentarily reaches a stationary state during the transition from the end of the RA stage to the early CD stage, but is characterised by fluctuations around a constant value (as depicted in figure 2b of Wei & Livescu Reference Wei and Livescu2012). After that, the growth rate continues to increase linearly with time during the CD stage. However, as the Reynolds number decreases, various stages of the single-mode RTI are observed to gradually disappear due to the viscous effects except in the nonlinear stage. Given the good consistency between the properties of a single-mode RTI by solving the Euler equation system and Navier–Stokes equation system at high Reynolds numbers, it is inferred that this new velocity saturation phenomenon corresponds to a transient stage from the end of the RA stage to the early CD stage, as illustrated in figure 1. During this stage, vortex pairs reach the mixing layers at the tip of the bubble, thereby limiting the increase in induced velocity. Meanwhile, the continuous augmentation of vortical motions near the centreline, where the shear velocity is the largest, quickly leads to complex vortex interactions, influencing the velocity induced by vorticity and causing significant fluctuations. Therefore, this phenomenon is termed the ‘secondary velocity saturation’ . To the best of the authors’ knowledge, this marks the first theoretical realisation of this particular stage.

It should be mentioned that the present model does not consider the effects of viscosity and surface tension, rendering it applicable only to situations at high Reynolds numbers. If the method introduced by Sohn (Reference Sohn2009) is employed to take into account the effects of viscosity and surface tension, the model can be further utilised to analyse the mechanism of the successive disappearance of various stages at low-Reynolds-number conditions.

Based on the physical implications of these two fixed points, it is further inferred that bubbles achieve the velocity saturation due to the standard equilibrium between buoyancy and drag forces at the QS stage. However, with the increase of neighbouring vorticity, the vorticity-induced velocity breaks this standard equilibrium and triggers the re-acceleration of bubbles. Furthermore, the combined effects of gravity and vorticity establish a new equilibrium at the interface, resulting in a novel SVS phenomenon in the early chaotic stage. Therefore, by simultaneously considering the potential force-induced component from potential flow and vorticity-induced component from rotational flow, the coupled system can accurately describe the bubble evolution up to the early chaotic stage of a single-mode RTI.