1. Introduction
Rayleigh–Taylor instability (RTI) evolves at the interface of two fluids with different densities when a heavier fluid is above a lighter one under an external acceleration normal to the interface, typically gravity. Similar to Rayleigh–Bénard convection – where buoyancy drives motion in a fluid layer heated from below (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Shishkina Reference Shishkina2021; Lohse & Shishkina Reference Lohse and Shishkina2023, 2024; Xia et al. Reference Xia, Huang, Xie and Zhang2023; Xia et al. Reference Xia, Chong, Ding and Zhang2025) – RTI is a fundamental buoyancy-driven instability. The occurrence of RTI can be seen in natural systems and engineering applications (Zhou Reference Zhou2017a , Reference Zhoub ), including the geophysical formation of volcanic islands (McGuire Reference McGuire1996; Zhang et al. Reference Zhang, Guo, Liu, Liu, Zhang, Lei, Zhao, Ma, Sepe and Ventura2018), atmospheric circulation (Wunsch & Ferrari Reference Wunsch and Ferrari2004; Boffetta & Mazzino Reference Boffetta and Mazzino2017), and the hydrodynamic behaviour of liquid films (Livescu et al. Reference Livescu, Roy and Schwartz2011).
The traditional studies on RTI typically rely on the continuum assumption. However, rarefaction effects become significant in many emerging problems, where the Knudsen number (
$Kn=l/L$
) quantifies the degree of rarefaction based on the molecular mean free path
$l$
and a characteristic length scale
$L$
. This is particularly relevant in astrophysics (Janka Reference Janka2012; Isobe et al. Reference Isobe, Miyagoshi, Shibata and Yokoyama2015; Ma et al. Reference Ma, Delamere and Otto2015), nuclear collisions (Bouras et al. Reference Bouras, El, Fochler, Niemi, Xu and Greiner2012), inertial confinement fusion (Betti & Hurricane Reference Betti and Hurricane2016; Larroche et al. Reference Larroche, Rinderknecht, Rosenberg, Hoffman, Atzeni, Petrasso, Amendt and Séguin2016), and the expansion of tin plasma under a low pressure environment in extreme-ultraviolet light lithography (Theofanous et al. Reference Theofanous, Li and Dinh2004; Bakshi Reference Bakshi2006; Rao & Basu Reference Rao and Basu2026; Chong et al. Reference Chong, Wang, Wang, Yan, Ding, Sun and Zhou2026). Overall, the continuum descriptions become invalid when
$Kn$
is sufficiently large, which can be achieved in high vacuum or microscale environments.
As early as the 1880s, Rayleigh (Reference Rayleigh1883) showed mathematically that when a heavier fluid overlies a lighter one, the interface between them becomes unstable. Subsequent studies demonstrated that the evolution of this interfacial instability depends on gravitational acceleration, density contrast and initial perturbations. Building on Rayleigh’s analysis, Taylor (Reference Taylor1950) experimentally quantified how the instability growth rate varies with perturbation wavelength, acceleration and density contrast, thereby laying the foundation for modern hydrodynamic instability theory. The evolution of RTI is commonly divided into three stages based on growth behaviour: the linear stage, the nonlinear growth stage, and the self-similar regime (Ristorcelli & Clark Reference Ristorcelli and Clark2004; Zhou & Jiang Reference Zhou and Jiang2016; Boffetta & Mazzino Reference Boffetta and Mazzino2017; Zhou et al. Reference Zhou2021). The relationship between the mixing layer width and time across these stages has been confirmed through both theoretical and experimental results (Layzer Reference Layzer1955; Zhou et al. Reference Zhou, Huang, Lu, Liu and Ni2016). Wei & Livescu (Reference Wei and Livescu2012) employed direct numerical simulations (DNS) to simulate two-dimensional RTI with a single-mode, fixed-wavelength perturbation at the initial interface, providing a detailed analysis of how the Reynolds number influences RTI growth rates across different stages.
Moreover, it was found that at low Reynolds numbers, both dynamic viscosity and diffusion exert a strong inhibitory effect on the growth of instabilities. Accordingly, Duff et al. (Reference Duff, Harlow and Hirt1962) derived a modified growth rate expression for gas-phase RTI that accounts for viscous and diffusive effects. To investigate the microscopic mechanisms underlying interfacial instability formation, Kadau et al. (Reference Kadau, Germann, Hadjiconstantinou, Lomdahl, Dimonte, Holian and Alder2004, Reference Kadau, Rosenblatt, Barber, Germann, Huang, Carlès and Alder2007, Reference Kadau, Barber, Germann and Alder2008, Reference Kadau, Barber, Germann, Holian and Alder2010) used molecular dynamics and direct simulation Monte Carlo (DSMC) to simulate RTI with an initially flat interface in continuous flow. Their results were consistent with theoretical predictions and experimental observations, while also showing that molecular scale fluctuations play an increasingly important role as system size decreases. Gallis et al. (Reference Gallis, Koehler, Torczynski and Plimpton2016) further applied DSMC to simulate RTI flows, and examined the evolution of flat and single-mode initial interfaces under varying Atwood numbers and gravitational accelerations at low
$Kn$
values. Their study demonstrated that across the linear, nonlinear and self-similar regimes, DSMC agree well with theoretical and experimental results, and that molecular level perturbations significantly influence the growth of interfacial perturbations.
Although extensive research has been conducted on RTI, the vast majority of studies rely on the continuum assumption. Consequently, studies on RTI in rarefied environments remain limited, despite their emerging relevance in modern applications. Kashkovsky et al. (Reference Kashkovsky, Kudryavtsev and Shershnev2019) employed the DSMC method to simulate an RTI system with an initially flat interface under low pressure conditions. At
$Kn=0.0036$
, the instability growth characteristics were found to closely match continuum flow predictions, a conclusion also supported by Shershnev et al. (Reference Shershnev, Kudryavtsev and Kashkovsky2023, Reference Shershnev, Kudryavtsev, Kashkovsky, Borisov and Poleshkin2024). Although their RTI study used a low-pressure gas, the Knudsen numbers were relatively small (
$Kn\lt 0.01$
). Consequently, these cases exhibit only relatively weak rarefaction effects. Sagert et al. (Reference Sagert, Bauer, Colbry, Howell, Staber and Strother2014, Reference Sagert, Howell, Staber, Strother, Colbry and Bauer2015) used an extended DSMC formulation to simulate RTI with fixed wavelength initial perturbations under different
$Kn$
. Their results indicated that at large
$Kn$
, enhanced molecular diffusion suppresses secondary growth of the interfacial instability. However, these studies did not account for the thermally driven most unstable wavelength, which increases with
$Kn$
(Duff et al. Reference Duff, Harlow and Hirt1962). When this most unstable wavelength exceeds the initially fixed perturbation wavelength, it dominates the RTI evolution, but this behaviour cannot be captured within a fixed computational domain. Only a limited number of studies have examined rarefaction effects in RTI. In particular, the coupling between molecular scale fluctuations in transitional flows and macroscopic hydrodynamic instabilities is poorly understood, making it challenging to assess how rarefaction alters RTI evolution.
According to the classical classification of rarefied gas flows (Bird Reference Bird1994; Gad-el-Hak Reference Gad-el-Hak1999), the continuum assumption is generally valid only for
$Kn \lt 0.001$
. In the slip flow regime (
$Kn \gt 0.001$
), DSMC results show better agreement with experimental measurements than the corresponding Navier–Stokes (NS) predictions (Stefanov et al. Reference Stefanov, Barber and Ota2005; Gallis et al.
Reference Gallis, Torczynski, Rader, Tij and Santos2006; Li et al. Reference Li, Peng, Wu, Ma, Tang and Liang2019). To systematically study the evolution of RTI in rarefied gases, we employ the DSMC method to simulate a two-dimensional Rayleigh–Taylor (RT) system with an initially flat interface over a range of Knudsen numbers. Moreover, considering that DSMC require a huge number of particles (in excess of hundreds of millions) to reduce statistical noise, a graphics processing unit (GPU) parallelised DSMC method is adopted to enhance computational efficiency (Gladkov et al. Reference Gladkov, Tapia, Alberts and D’Souza2012; Norman et al. Reference Norman, Valentini and Schwartzentruber2013). The objective of this study is to examine how the degree of rarefaction influences the evolution of RTI and the threshold mechanisms for instability onset in rarefied gases.
The remainder of this paper is organised as follows. Section 2 presents the governing equations and numerical set-up. Section 3 investigates the evolution of RTI under varying Knudsen numbers, including molecular transport behaviour, macroscopic flow characteristics, and the competition between diffusive damping and instability growth. The influence of the computational height is further examined by varying the Froude number
$Fr$
. Section 4 summarises the main conclusions and outlook.
2. Numerical methods and problem set-up
2.1. The DSMC
Rarefied gas flows have been investigated using a variety of numerical approaches in recent years (Xu & Huang Reference Xu and Huang2010; Gorji et al. Reference Gorji, Torrilhon and Jenny2011; Fei et al. Reference Fei, Zhang, Li and Liu2020; Su et al. Reference Su, Zhu, Wang, Zhang and Wu2020). Among these methods, the DSMC method developed by Bird (Reference Bird1963, Reference Bird1994) is particularly well suited for problems involving strong non-equilibrium effects, as it directly resolves molecular transport processes such as viscosity, thermal conductivity and mass diffusion. Owing to these characteristics, DSMC has been extensively validated against theoretical predictions over a wide range of Knudsen numbers, from near-continuum to highly rarefied regimes (Garcia Reference Garcia1986; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; McMullen et al. Reference McMullen, Krygier, Torczynski and Gallis2022; Ma et al. Reference Ma, Yang, Chen, Feng, Cui and Zhang2024). In the context of RTI at finite Knudsen numbers, where compressibility and non-equilibrium molecular effects may play a significant role, DSMC provides a physically consistent framework for capturing the underlying gas dynamics.
Based on the above considerations, DSMC is employed in this work to simulate a two-dimensional RTI with an initially flat interface at various Knudsen numbers. A schematic of the computational set-up and the corresponding initial and boundary conditions is shown in figure 1. The heavy fluid initially occupies the upper half of the computational domain, while the light fluid occupies the lower half. Gravity acts in the negative
$y$
-direction and normal to the interface. The boundaries at
$y = 0$
and
$y = H$
are treated as solid walls, whereas periodic boundary conditions are prescribed at
$x = 0$
and
$x = L$
.
Sketch of the set-up of the RTI system.
In DSMC, particle transport and collisions are treated separately, which requires the time step to be smaller than the mean collision time, and the cell size to be smaller than the molecular mean free path
$l$
(Bird Reference Bird1994; Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020; Ma et al. Reference Ma, Yang, Chen, Feng, Cui and Zhang2024). The initial particle positions and velocity distributions are prescribed according to the macroscopic properties of the gas. To ensure both spatial uniformity and randomness at the particle level, equal numbers of particles are placed in each computational cell, and their positions are randomly assigned within the cell. The particle velocities are sampled according to the Maxwell–Boltzmann distribution. The corresponding probability density function
$f(u)$
for the particle speed is given by (Cercignani Reference Cercignani1988)
\begin{align} f(u) = 4\pi \left (\frac {m_i}{2\pi k_B T}\right )^{3/2} \exp^{-\tfrac {m_i u^2}{2k_B T}} u^2, \end{align}
where
$m_i$
is the mass of a molecule of species
$i$
,
$k_B$
is the Boltzmann constant,
$T$
is the absolute temperature, and
$u$
denotes the magnitude of the molecular speed.
Since particle transport and collisions are decoupled in DSMC, particle motion within a time step
$\Delta t$
is treated as ballistic motion. The Maxwell diffuse reflection model can be used to describe particle–wall interactions (Bird Reference Bird1994; Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020; McMullen et al. Reference McMullen, Krygier, Torczynski and Gallis2023). After particle transport and boundary processing, particles are reassigned to their corresponding collision cells. To model particle–particle collisions, each collision cell is further subdivided into several subcells. Within each subcell, particle pairs are selected probabilistically according to the adopted collision model. In this study, the variable soft-sphere (VSS) model is employed to calculate the particle collision frequency (Koura & Matsumoto Reference Koura and Matsumoto1992; Stevens & Hrenya Reference Stevens and Hrenya2005; Gallis et al. Reference Gallis, Torczynski, Krygier, Bitter and Plimpton2021; Ma et al. Reference Ma, Yang, Chen, Feng, Cui and Zhang2024). For a multi-component gas mixture consisting of
$s$
species, the total number of collisions per unit volume and per unit time is given by
\begin{align} N_c=\frac {1}{2}\sum _{i=1}^s(n_i z_i), \end{align}
where
$n_i$
denotes the number density of molecules of component
$i$
, and
$z_i$
denotes the collision frequency of molecules of component
$i$
. The collision frequency is expressed as
\begin{align} z_i=\sum _{j=1}^s(n_j\, \overline {\pi d_{ij}^2 u_{r_{ij}}}), \end{align}
where
$d_{ij}$
is the effective collision diameter between molecular species
$i$
and
$j$
, and
$u_{r_{ij}}$
is the relative speed between pairs of collision particles. The overline denotes an ensemble average over the distribution of molecular velocities.
To ensure that molecular diffusion and viscosity are correctly represented in the DSMC framework, the VSS model combined with the first-order Chapman–Enskog expansion is adopted to relate collision parameters to macroscopic transport coefficients (Koura & Matsumoto Reference Koura and Matsumoto1992; Bird Reference Bird1994). The dynamic viscosity
$\mu$
is given by
\begin{align} \mu =\frac {5(\alpha +1)(\alpha +2)(\pi m k_B)^{1 / 2}(4 k_B / m)^\xi T^{1 / 2+\xi }}{16 \alpha \varGamma (4-\xi ) \sigma _{\mathit{T}, \mathit{ref}} u_{\mathit{r}, \mathit{ref}}^{2 \xi }},\end{align}
and the binary diffusion coefficient
$D_{12}$
is expressed as
\begin{align} D_{12}=\frac {3\left (\alpha _{12}+1\right ) \pi ^{1 / 2}\left (2 k_B T / m_{\mathit{r}}\right )^{\omega _{12}}}{16 \varGamma \left (7 / 2-\omega _{12}\right ) n\left (\sigma _{\mathit{T}, \mathit{ref}}\right )_{12} u_{\mathit{r}, \mathit{ref}}^{2(\omega _{12}-1 / 2)}} , \end{align}
where
$\alpha$
and
$\omega$
are the VSS scattering parameters, subscripts 1 and 2 denote the two molecular species, and
$\omega = 1/2 + \xi$
. Here, the reduced mass is
$m_r = m_1 m_2 / (m_1 + m_2)$
, the reference collision cross-section is
$\sigma _{\mathit{T,ref}} = \pi d^2$
, and
$u_{\mathit{r,ref}}$
is the reference relative velocity. Additional parameter definitions and the properties of the
$\varGamma$
function are provided in Bird (Reference Bird1994).
A known limitation of the DSMC method is the presence of statistical noise arising from the finite number of simulation particles. When each particle represents a real molecule, the statistical noise inherently accounts for the influence of molecular thermal fluctuations on the macroscopic flow field (Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Liao & Qin Reference Liao and Qin2022, Reference Liao and Qin2025). The statistical uncertainty decreases with increasing sample size, and can be reduced through ensemble averaging. In this study, all particles residing in a computational cell over a prescribed sampling interval are collected for ensemble averaging, from which macroscopic quantities such as density, velocity, pressure and temperature are evaluated. Since the RT system considered here involves two molecular species, their contributions are treated separately in both collision calculations and statistical evaluations of macroscopic variables.
2.2. Verification of the DSMC code
The evolution of the RTI is governed by gravitational acceleration and fluid properties. According to linear stability theory (Taylor Reference Taylor1950; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Boffetta & Mazzino Reference Boffetta and Mazzino2017), the perturbation growth rate is
\begin{align} \gamma _0 \approx \sqrt {A_t g k},\quad k = 2\pi /\lambda , \end{align}
where
$g$
is the gravitational acceleration,
$k$
is the wavenumber, and
$\lambda$
is the perturbation wavelength. The Atwood number
$A_t = (\rho _h - \rho _l)/(\rho _h + \rho _l)$
quantifies the density contrast between the heavy (
$\rho _h$
) and light (
$\rho _l$
) fluids. The corresponding characteristic time scale is
$\tau = 1/\gamma _0$
(Youngs Reference Youngs1984). Equation (2.6) indicates that shorter perturbation wavelengths lead to larger growth rates. However, viscosity and diffusion inhibit the growth of short-wavelength perturbations. For an initially flat interface, molecular thermal fluctuations are amplified and trigger the instability. These perturbations encompass a wide range of wavelengths, with the most unstable mode exhibiting the fastest growth rate. The most unstable wavelength
$\lambda _m$
of the RT system is given by (Sharp Reference Sharp1984; Youngs Reference Youngs1984)
\begin{align} \lambda _m \approx 4\pi \left (\frac {\bar {\nu }^2}{A_t g}\right )^{1/3}, \end{align}
where
$\bar {\nu }= (\mu _h + \mu _l)/(\rho _h + \rho _l)$
is the mean kinematic viscosity, and
$\mu _h$
and
$\mu _l$
denote the dynamic viscosities of the heavy and light fluids, respectively. During the early linear stage of the instability, the perturbation amplitude
$a$
evolves according to linear stability theory given by (Taylor Reference Taylor1950; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016)
\begin{align} a = a_0 \cosh (t \sqrt {A_t g k}), \end{align}
where
$a_0$
is the initial perturbation amplitude.
In a gaseous RT system, the kinematic viscosity, diffusion coefficient and density contrast strongly influence the growth of initial interfacial perturbations. Duff et al. (Reference Duff, Harlow and Hirt1962) corrected the growth rate by accounting for the influence of diffusion on interface instability. The growth rate
$\gamma _0$
can be expressed as
\begin{align} \gamma _0 = \sqrt {A_t g k + \bar {\nu }^2 k^4} - (\bar {\nu } + D_{12}) k^2, \end{align}
where
$\gamma _0$
has maximum value at a specific wavenumber
$k$
. This wavenumber corresponds to the most unstable wavelength
$\lambda _m$
for an initially flat interface. Accordingly, (2.9) is used to determine
$\lambda _m$
in this study.
In the nonlinear stage (
$ka \sim 1$
), following Layzer (Reference Layzer1955), Goncharov (Reference Goncharov2002) and Mikaelian (Reference Mikaelian2014), the bubble amplitude
$a_b$
and spike amplitude
$a_s$
are expressed as
\begin{align} a_b(A_t) = a_0 + \left (\frac {3 + A_t}{3(1 + A_t) k}\right ) \ln \left [ \cosh \left ( t \frac {\sqrt {6 g k A_t (1 + A_t)}}{3 + A_t} \right ) \right ], \end{align}
\begin{align} a_s(A_t) = a_b(A_t) \left ( 1 + \left [ 0.4 + 0.6 A_t^{10} \right ] \left [ \left ( \frac {a_s(1)}{a_b(A_t)} \right )^{A_t} - 1 \right ] \right ). \\[9pt] \nonumber \end{align}
In the self-similarity stage, experimental data indicate that the growth rates of bubbles and spikes can be described by a simple model (Cook et al. Reference Cook, Cabot and Miller2004; Ristorcelli & Clark Reference Ristorcelli and Clark2004; Boffetta & Mazzino Reference Boffetta and Mazzino2017)
\begin{align} a_{b,s} = \alpha _{b,s} A_t g t^2, \end{align}
where
$\alpha _b$
is approximately 0.005, and
$\alpha _s$
depends on
$A_t$
. The total interface amplitude is given by
$a = a_b + a_s$
.
To validate the code, a GPU parallel DSMC method was employed to simulate the RT system with initially fixed perturbation wavelengths
$\lambda$
. The relevant gas and RT system parameters are summarised in tables 1 and 2, respectively (Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016). The evolution of the interface height over time is shown in figure 2. The amplitude of interface perturbations as a function of time agrees well with both linear and nonlinear theories, and is consistent with the DSMC results (Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016). Further comparison with NS results at different Knudsen numbers shows that DSMC and NS predictions are in good agreement under continuum flow conditions, as shown in figure 2(c). However, in the slip flow regime (
$0.001 \lt Kn \lt 0.013$
), the NS results deviate significantly from the DSMC predictions, as shown in figure 2(d). This discrepancy arises because the NS equations cannot accurately capture the non-equilibrium molecular transport processes occurring in the slip flow regime, thereby demonstrating the advantage of DSMC in this range. To further verify mesh independence, the RT system with an initially flat interface was simulated using meshes of different sizes. The relevant RT system parameters are summarised in table 3 (case 3). The results show that smaller grid sizes do not affect the characteristics of the interface evolution, as shown in figure 2(b). Therefore, this study adopts the standard DSMC requirements for mesh size and time step. The cell size satisfies
$\Delta {x} \lt l/2$
, and the simulation time step satisfies
$\Delta {t}\lt 0.2\Delta {x}/V_{th}$
(Bird Reference Bird1994; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020), where the molecular mean free path is given by
$l=1/\sqrt {2}\,\pi \bar {d}^2n_{h,l}$
and
$\bar {d}=0.5(d_h+d_l)$
.
Properties of pure gas at
$P=101\,325$
Pa and
$T=273.15$
K. Here,
$\omega$
and
$\alpha$
are the parameters of the VSS molecular collision model (Bird Reference Bird1994).
Table 1. Long description
A table comparing properties of two gases, Helium and Argon, under specific conditions. The table has 8 columns and 2 rows. Column headers are Gas, m (kg), d (m), ρ (kg m^-3), μ (Pa s), v (m^2 s^-1), ω, and α. Row 1: Gas, He; m (kg), 6.65 x 10^-27; d (m), 2.30 x 10^-10; ρ (kg m^-3), 0.179; μ (Pa s), 1.865 x 10^-5; v (m^2 s^-1), 10.4 x 10^-5; ω, 1.26; α, 0.66. Row 2: Gas, Ar; m (kg), 66.35 x 10^-27; d (m), 4.11 x 10^-10; ρ (kg m^-3), 1.781; μ (Pa s), 2.117 x 10^-5; v (m^2 s^-1), 1.19 x 10^-5; ω, 1.4; α, 0.81.
Validation of RTI simulation at 1 mm
$\times$
4 mm domain. Here,
$P_{\mathit{ref}}$
is the initial interface reference pressure, and
$\lambda$
is the initial interface perturbation wavelength.
Table 2. Long description
The table has one row and seven columns. The columns are labeled Gases, A I, g (m s^-2), λ (mm), a 0 (μm), P ref (Pa), and T (K). The row contains the following values: Ar/He, 0.82, 1.0 x 10^8, 1.00, 10.0, 101325, 273.15.
Verification of the DSMC. (a) The normalised perturbation amplitude
$k(a-a_0)$
as a function of time
$t/\tau$
. (b) Evolution of the interface height in an RT system with an initially flat interface under different grid resolutions for
$Kn = 0.005$
. (c) Comparison of DSMC and NS predictions for the interface evolution height at
$Kn = 0.001$
. (d) Comparison of DSMC and NS predictions for the interface evolution height at
$Kn = 0.009$
.
2.3. Problem configuration
In this study, the kinematic viscosity and diffusion coefficient increase significantly as the gas becomes more rarefied. According to (2.9), the most unstable wavelength
$\lambda _m$
also increases. We first examined the sensitivity to the domain height ratio
$H/\lambda _m$
. According to Wei & Livescu (Reference Wei and Livescu2012), linear stability theory is valid only when the interface height satisfies
$a/\lambda _m \leqslant 0.4$
in the continuum regime. To investigate flow evolution beyond this linear regime, we fixed
$H/\lambda _m=0.5$
and simulated RTI for various
$Kn$
. The specific parameters for these simulations are summarised in cases 1–6 of table 3. To further examine the influence of the domain height, Betti et al. (Reference Betti, Goncharov, McCrory and Verdon1998) introduced the Froude number
$Fr$
as a control parameter for perturbation growth, defined as
$Fr=V_{th}^2/(gH)$
. Similarly, Zhou et al. (Reference Zhou, Huang, Lu, Liu and Ni2016) used the Rayleigh number to characterise the relationship between domain size and instability. Following this approach, simulations were performed for different
$Fr$
and
$Kn$
, as summarised in cases 7–48 of table 3, where
$V_{th} = \sqrt {2k_BT_{\mathit{ref}}/m}$
is the molecular thermal velocity, and
$L=1.5\lambda _m$
. In the RT system, the heavy and light fluids are argon (Ar) and helium (He), respectively, and their properties at standard conditions are listed in table 1. The reference interface temperature is
$T_{\mathit{ref}} = 273.15\ \mathrm{K}$
. This significantly enhanced gravitational acceleration is applied to the system in order to accelerate the growth of RTI within the confined computational domain (Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020). The initially flat interface is located at
$y = H/2$
. Then the pressure distribution is prescribed as
\begin{align} P_{h,l}(y) = P_{\mathit{ref}} + \rho _{h,l} g (0.5H -y). \end{align}
Using the ideal gas equation of state, the temperature distribution along
$y$
is given by
\begin{align} k_B\, T_{h,l}(y) = P_{h,l}(y)/n_{h,l}, \end{align}
where
$n_h = n_l$
at the reference state. The number densities are computed from the ideal gas relation
$n_{h,l} = P_{\mathit{ref}}/(k_B T_{\mathit{ref}})$
. Collision sampling at each time step is based on the ‘no time counter’ scheme (Bird Reference Bird1989; Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020; Ma et al. Reference Ma, Yang, Chen, Feng, Cui and Zhang2024). To reduce statistical uncertainty, macroscopic variables are obtained by averaging 20 time step samples, and performing eight independent ensemble implementations (Goshayeshi et al. Reference Goshayeshi, Di Staso, Toschi and Clercx2020; Ma et al. Reference Ma, Yang, Chen, Feng, Cui and Zhang2024). The simulation time exceeds the characteristic time at which interface growth ceases. A quasi-two-dimensional computational domain with thickness 100 nm is used to model the RT system, with each computational particle representing 1000 real gas molecules.
Detailed parameters of the numerical simulations for different
$Kn$
and
$Fr$
. Here,
$l$
is the molecular mean free path, and
$\tau$
is the characteristic time scale of the most unstable wavelength.
Table 3. Long description
A table with three columns and multiple rows. The columns are labeled with symbols representing different parameters. The first column lists values for a specific parameter, the second column lists corresponding values for another parameter, and the third column provides additional related values. The table details parameters for numerical simulations under varying conditions, focusing on molecular mean free path and characteristic time scale.
3. Results and discussion
In this section, we first investigate the influence of the Knudsen number on the RTI evolution under a fixed geometric condition
$H/\lambda _m=0.5$
(§§ 3.1–3.4). Based on these results, we further extend the analysis to examine the combined effects of the Knudsen number
$Kn$
and
$Fr$
on the instability dynamics (§ 3.5).
3.1. Evolution of the interface at different
$Kn$
In RTI, the concentration field better captures the morphology of the evolving mixing layer between the heavy and light gases; we hereafter refer to this as the ‘interface’. Figure 3 shows instantaneous snapshots of the normalised concentration field
$\alpha _v$
for different
$Kn$
. Here,
$\alpha _v=n_h/(n_h+n_l)$
, where
$n_h$
and
$n_l$
denote the number densities of the heavy and light fluids, respectively. Each snapshot displays a subregion of the computational domain with a 3 : 2 aspect ratio. As shown in figures 3(
$a{,}d{,}g{,}j$
), when
$Kn=0.001$
, the interface evolves from its initially flat state into an unstable one as the fluctuating amplitude grows over time. Figures 3(
$j{,}k{,}l$
) demonstrate the reduction in the fluctuating amplitude as
$Kn$
increases, indicating a weakening of the interfacial instability.
At
$Kn=0.011$
, the instability is no longer observable, and a reversal of the concentration field occurs at
$t=10\tau$
. Our results suggest that the instability at the interface weakens with increasing
$Kn$
, and is totally suppressed beyond a critical value. Additionally, the interface thickness also increases with
$Kn$
, as is evident from figures 3(d–f), reflecting enhanced diffusion between light and heavy molecules at large
$Kn$
. This behaviour is consistent with the findings of Sagert et al. (Reference Sagert, Howell, Staber, Strother, Colbry and Bauer2015).
Typical snapshots of the normalised concentration field
$\alpha _v$
at different
$Kn$
. The black solid line represents the interface position at
$\alpha _v=0.5$
, and
$\alpha _v$
represents the ratio of the number density
$n_h$
of the heavy fluid to the total number density
$n$
(
$\alpha _v=n_h/(n_h+n_l)$
). These correspond to
$t=4\tau$
in (
$a{-}c$
),
$t=6\tau$
in (
$d{-}f$
),
$t=8\tau$
in (
$g{-}i$
),
$t=10\tau$
in (
$j{-}l$
). Here, (
$a, d, g, j$
), (
$b, e, h, k$
) and (
$c, f, i, l$
) correspond to
$Kn=0.001$
, 0.005 and 0.011, respectively.
To quantitatively assess the influence of rarefaction on the instability, the fluctuating amplitude
$a$
as a function of time for various
$Kn$
is shown in figure 4. The fluctuating amplitude
$a$
is computed as the peak-to-trough distance of the low-pass filtered
$\alpha _v = 0.5$
contour (Harrington Reference Harrington1993), where the filter window is
$\lambda _m/20$
, and
$\lambda _m$
is the most unstable wavelength. Figure 4(a) demonstrates that the growth rate of the fluctuating amplitude decreases as
$Kn$
increases. In the near-continuum regime, the growth of the fluctuating amplitude transitions from an initial diffusive regime, where
$a/\lambda _m \sim t^{1/2}$
, to a subsequent mixing regime with
$a/\lambda _m \sim t^{2}$
(Taylor Reference Taylor1950; Wei & Livescu Reference Wei and Livescu2012; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Zhou et al. Reference Zhou2021).
For
$Kn \gtrsim 0.011$
, however, the initial amplitude growth is followed by decay, and the scaling of the mixing regime is absent. As shown in figures 3(
$i{,}l$
), this decay results from the interaction of the perturbation with the upper and lower domain boundaries. The maximum achieved amplitude
$a_{max }$
, plotted against
$Kn$
in figure 4(
$b$
), decreases with increasing
$Kn$
, and plateaus for
$Kn\gtrsim 0.011$
. This demonstrates that a larger Knudsen number suppresses perturbation growth, establishing
$Kn \approx 0.011$
as the critical value for the onset of the RTI in this system.
(
$a$
) Evolution of the perturbation amplitude
$a$
as a function of time for various
$Kn$
. The dashed lines represent the diffusive-growth scaling law
$a \sim t^{0.5}$
, and the self-similar growth scaling law
$a \sim t^{2}$
. (
$b$
) Maximum perturbation amplitude
$a_{max }$
at different
$Kn$
.
Figure 4. Long description
Panel A: A line graph shows the evolution of the perturbation amplitude as a function of time for various values. The x-axis represents time normalized by a factor, and the y-axis represents the perturbation amplitude normalized by a wavelength. Different symbols and colors represent different values, ranging from 0.001 to 0.013. The dashed lines represent the diffusive-growth scaling law and the self-similar growth scaling law. Panel B: A scatter plot shows the maximum perturbation amplitude at different values. The x-axis represents the values, and the y-axis represents the maximum perturbation amplitude normalized by a wavelength. Different symbols and colors represent different values, similar to Panel A.
3.2. Mean squared displacement of molecules at different
$Kn$
To further elucidate how rarefied gases suppress interfacial instability in the RT system, we examine the molecular translational motion using the mean squared displacement (MSD), a standard measure of particle or object motion in dynamical systems (Michalet Reference Michalet2010; Huang et al. Reference Huang, Chavez, Taute, Lukić, Jeney, Raizen and Florin2011; Peter Reference Peter2011; Schirripa & Luin Reference Schirripa and Luin2024). The molecular MSD is defined as
\begin{align} \langle \boldsymbol{r}^2(t)\rangle =\frac {1}{N} \sum _{i=1}^N \left [ \boldsymbol{x}_i(t+\Delta t)-\boldsymbol{x}_i(t) \right ]^2, \end{align}
where
$N$
is the number of particles,
$\Delta t$
is the lag time, and
$\boldsymbol{x}_i(t)$
denotes the position of particle
$i$
at time
$t$
. Since gravitational acceleration and wall confinement strongly affect motion in the
$y$
-direction, the MSD is evaluated using only the
$x$
-direction displacement, which is periodic. This isolates molecular motion from gravitational and boundary influences, providing a clear characterisation of rarefaction-induced transport.
Figures 5(
$a{,}b$
) show the molecular trajectories over the same characteristic time
$t=\tau$
for
$Kn=0.001$
and
$Kn=0.013$
, respectively. The red and cyan dots denote the final positions of heavy and light fluid molecules, superimposed on the
$\alpha _v$
field at
$t=\tau$
. Because the thermal velocity of the heavy molecules is lower than that of the light molecules, their displacement over the same time interval is correspondingly smaller. At
$Kn = 0.001$
, molecular trajectories display frequent directional changes, characteristic of collision-dominated molecular motion. Conversely, at
$Kn = 0.013$
, the trajectories are predominantly straight, indicating a transition towards ballistic transport over the same time scale.
Trajectories of molecules: (
$a$
)
$Kn=0.001$
and (
$b$
)
$Kn=0.013$
. The red and cyan dots represent the ending points of the trajectories of heavy and light fluid molecules, respectively. The background is the
$\alpha _v$
field at
$t=\tau$
. (
$c$
) The MSD of molecules with time
$t/\tau$
at different
$Kn$
. The blue and orange dashed lines indicate the diffusive scaling
$\langle \boldsymbol{r}^2\rangle \sim t$
and ballistic scaling
$\langle \boldsymbol{r}^2\rangle \sim t^2$
, respectively. (
$d$
) Normalised MSD of molecules with time
$t/\tau$
at different
$Kn$
.
Figure 5. Long description
Panel A: A scatter plot shows the trajectories of molecules with red and cyan dots representing the ending points of heavy and light fluid molecules, respectively. The background is the field at a specific time. Panel B: Another scatter plot shows the trajectories of molecules at a different Knudsen number. Panel C: A line graph displays the mean squared displacement (MSD) of molecules over time at different Knudsen numbers, with blue and orange dashed lines indicating diffusive and ballistic scaling. Panel D: A line graph shows the normalized MSD of molecules over time at different Knudsen numbers.
This qualitative contrast is quantified in figure 5(
$c$
), which compares the MSD across different Knudsen numbers. At early times, the MSD follows the ballistic scaling, after which it gradually transitions to diffusive behaviour once a sufficient number of collisions have occurred. The diffusion coefficient
$D$
can be obtained via the classical Einstein relation for diffusion (Einstein Reference Einstein1905):
\begin{align} \langle \boldsymbol{r}^2(t)\rangle =2Dt. \end{align}
The MSD at different
$Kn$
is normalised by the corresponding diffusion coefficient, as shown in figure 5(
$d$
). As
$Kn$
increases, the characteristic time required for molecules to transition from ballistic to diffusive motion becomes longer. This behaviour indicates that intermolecular collisions occur less frequently at large
$Kn$
, resulting in an extended duration of super-diffusive motion (
$\langle \boldsymbol{r}^2(t)\rangle \sim t^{\zeta },\ \zeta \gt 1$
).
3.3. Flow field characteristics at different
$Kn$
How do molecular level alterations in rarefied gases influence macroscopic flow field characteristics? The density distribution from the flow field provides a direct measure of these molecular level effects. Figure 6(
$a$
) shows the flow field density distribution at
$t=5\tau$
for different
$Kn$
. Here,
$\langle \rho \rangle _x$
denotes the density averaged in the
$x$
-direction, and
$\rho _h$
is the heavy fluid density at
$t=0$
. The inset presents the temporal evolution of the dimensionless density at
$Kn=0.005$
.
At the same characteristic time, the interfacial density gradient decreases as
$Kn$
increases, indicating that large
$Kn$
enhances diffusion within the RT system. This observation is consistent with the results shown in figure 3. Combined with the MSD results in figure 5, increasing the
$Kn$
promotes the super-diffusive behaviour of molecules, thereby enhancing interfacial mixing between light and heavy molecules. Macroscopically, this results in faster diffusion and a more rapid decay of the density gradient.
(
$a$
) Dimensionless density distribution with different
$Kn$
at
$t=5\tau$
, where
$\langle \rho \rangle _x$
is the average result in the
$x$
-direction. The inset shows the dimensionless density distribution with different times at
$Kn=0.005$
. (
$b$
) Temporal evolution of the effective Atwood number
$A_{t,\mathit{eff}}$
at different
$Kn$
.
The effect of density distribution on instability can be quantified by the effective Atwood number
$A_{t,\mathit{eff}}$
. According to (2.6),
$A_t$
also determines the growth rate of interface perturbations. The effective Atwood number is defined as
\begin{align} A_{t,\mathit{eff}}=\frac {\langle \rho \rangle _t-\langle \rho \rangle _b}{\langle \rho \rangle _t+\langle \rho \rangle _b}, \end{align}
where
$\langle \rho \rangle _t$
and
$\langle \rho \rangle _b$
denote the spatially averaged densities in
$y\gt 0.5H$
and
$y\lt 0.5H$
, respectively. The temporal evolution of
$A_{t,\mathit{eff}}$
for different
$Kn$
is shown in figure 6(
$b$
). The value of
$A_{t,\mathit{eff}}$
gradually decreases over time, and this decay becomes steeper as
$Kn$
increases. For
$Kn\geqslant 0.008$
,
$A_{t,\mathit{eff}}$
can even become negative. When
$A_{t,\mathit{eff}}\leqslant 0$
, the density distribution is completely reversed relative to the initial state, as confirmed by figure 3(
$l$
). From (2.6), it is evident that when
$A_{t,\mathit{eff}} \lt 0$
, the growth rate of the perturbation becomes imaginary. As demonstrated in § 3.2, the molecular MSD reveals that rarefaction significantly enhances the super-diffusive transport of molecules. Specifically, at higher
$Kn$
, the reduced collision frequency extends the duration of the ballistic regime, allowing molecules to travel longer distances without scattering. This enhanced molecular mobility facilitates rapid interfacial mixing between the heavy and light fluids. Consequently, this vigorous mixing at the molecular scale macroscopically manifests as an accelerated decay of the density gradient. The system transitions from a buoyancy-dominated regime to a diffusion-dominated regime so rapidly that the initial density stratification is destroyed before the RTI can develop, leading to the observed inversion of the concentration field.
Viscosity also plays an important role in the RT system. At a fixed reference temperature, the kinematic viscosity is inversely proportional to the density. Therefore, when the Knudsen number increases (due to the decrease in pressure), the density decreases, resulting in an increase in kinematic viscosity. This further suppresses the growth of interfacial instability perturbations for large Knudsen numbers. However, the condition whether
$A_t$
vanishes is not a sufficient criterion for stability, as instability may have already developed prior to the decay of the Atwood number. The key is the relative timing of instability onset and diffusive decay, a competition governed by the Knudsen number.
3.4. Stability analysis
In continuum RTI, interface growth transitions from an initial diffusive regime to a mixing regime (Wei & Livescu Reference Wei and Livescu2012; Gallis et al. Reference Gallis, Koehler, Torczynski and Plimpton2016; Zhou Reference Zhou2017a
; Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, Knaepen and De Wit2018). These regimes correspond physically to mass transport by diffusion and by buoyancy-driven convection. Consequently, mass transport in RTI evolves from diffusion-dominated to instability-dominated dynamics. However, at large Knudsen numbers, intermolecular diffusion is enhanced (Sagert et al. Reference Sagert, Howell, Staber, Strother, Colbry and Bauer2015; Keenan & Sauppe Reference Keenan and Sauppe2023; Guo et al. Reference Guo, Wu and Zhang2024), raising the question of whether instability can emerge against strong diffusive damping. To resolve this, we characterise the competition between diffusion and instability by defining two time scales: a diffusion time scale
$\tau _D$
, and an instability time scale
$\tau _I$
.
Through dimensional analysis of the diffusion equation, the diffusion width satisfies
$a^2(t) \sim 2Dt$
, where
$D$
is the diffusion coefficient. Since the present study considers diffusion primarily in the vertical direction, the corresponding diffusion time scale can be expressed as
$\tau _D \approx a^2(t)/(2D_{12})$
. Here,
$a(t)$
is the instantaneous perturbation amplitude of the interface. This quantity is directly extracted from the DSMC results as the peak to trough vertical distance of the low-pass filtered
$\alpha _v = 0.5$
concentration contour, as described in
$\S$
3.1. During the diffusive growth regime, the diffusion width
$\sqrt {2D_{12}t}$
and the perturbation amplitude
$a(t)$
remain of the same order of magnitude. The diffusion time scale
$\tau _D$
characterises the time required for molecular diffusion to produce an interface displacement comparable to the instantaneous perturbation amplitude.
The instability time scale
$\tau _I$
can be characterised by the growth rate
$\gamma$
(Youngs Reference Youngs1984). Taking into account the effect of viscosity and diffusion, the growth rate can be expressed as (Duff et al. Reference Duff, Harlow and Hirt1962)
\begin{align} \gamma =\sqrt {\frac {A_tgk}{\psi }+\bar {\nu }^2k^4}-(\bar {\nu }+D_{12})k^2, \end{align}
where
$\psi$
is the growth suppression factor, with
$\psi =1$
at
$t=0$
. Here,
$\psi$
represents the influence of the interface density gradient on stability, and it is obtained by solving the inviscid Rayleigh equation (Rayleigh Reference Rayleigh1883):
\begin{align} \frac {{\rm d}}{{\rm d}y}\left (\rho \frac {{\rm d}w}{{\rm d}y}\right )=wk^2 \left (\rho -\frac {g}{\gamma ^2_0}\frac {{\rm d}\rho }{{\rm d}y} \right ),\quad \gamma ^2=\frac {A_tgk}{\psi }, \end{align}
where
$w$
is the perturbation velocity. The density
$\rho$
and its gradient
${\rm d}\rho /{\rm d}y$
are obtained using a self-similar model of density diffusion:
\begin{align} \rho (y,t)=\frac {1}{2}(\rho _h+\rho _l) \left [1+2A_t\sqrt {\pi }\int _0^{{y}/{\epsilon }}\exp {(-\xi ^2)}\,{\rm d}\xi \right ], \end{align}
where the integral represents the error function with respect to
$y/ \epsilon$
, and the self-similar variable
$\epsilon =\sqrt {4Dt}$
denotes the diffusion width at time scale
$t$
. Equation (3.6) is substituted into (3.5), and
$\psi$
is obtained numerically, following the procedure given in Duff et al. (Reference Duff, Harlow and Hirt1962). Therefore, the growth rate
$\gamma$
is obtained, and the instability time scale is defined as
$\tau _I=1/\gamma$
.
To calculate the growth suppression factor
$\psi$
at different
$Kn$
, the density distribution of the RT system varies differently over time at different
$Kn$
, therefore the self-similarity model for density distribution cannot be applied. The density and density gradient in (3.5) are obtained from the numerical results of DSMC simulations. The growth suppression factor
$\psi$
is computed through numerical iteration applied to (3.5) (Duff et al. Reference Duff, Harlow and Hirt1962).
Figure 7(
$a$
) shows the growth-suppression factor
$\psi$
as a function of the dimensionless diffusion thickness
$k \epsilon$
for different
$Kn$
, where the theoretical values are obtained using the self-similar density diffusion model in (3.6). For comparison, the corresponding continuum DNS results are obtained by solving the NS equations under the same initial and boundary conditions. The growth inhibition factor
$\psi$
obtained from DSMC and DNS agrees well with the theoretical prediction during the early stage of instability evolution in the near continuum regime. From a molecular perspective,
$\psi$
physically represents the efficiency of momentum transfer across the interface. In the near continuum regime (low
$Kn$
), frequent intermolecular collisions maintain a sharp density gradient, resulting in a lower
$\psi$
, and allowing hydrodynamic instabilities to develop. However, as rarefaction increases, the enhanced molecular mean free paths (as shown in the MSD analysis in § 3.2) promote direct molecular transport across the interface. This microscopic mixing diffuses the density gradient, effectively increasing
$\psi$
and damping instability growth. To calculate the growth rate
$\gamma$
at different
$Kn$
, the kinematic viscosity
$\bar {\nu }$
and diffusion coefficient
$D_{12}$
in the rarefied gas are considered, which can be obtained from (2.4) and (2.5). The growth rate
$\gamma$
for different
$Kn$
was obtained by substituting the calculated
$\psi$
into (3.4).
Figure 7(
$b$
) shows the growth rate
$\gamma$
at different
$Kn$
. The growth rate decreases over time, and the larger the
$Kn$
value, the more rapidly the growth rate decays to zero. At late times, the growth rate becomes negative, indicating a decay of the perturbation amplitude. This behaviour is attributed to enhanced interfacial mixing under rarefied conditions, which weakens the density gradient and reduces the effective Atwood number, eventually causing diffusive transport to dominate instability growth. This implies that large
$Kn$
suppresses the growth of perturbations, which yields the same conclusion as that obtained in figure 6(
$b$
). Notably, the growth rate
$\gamma$
decreases as
$Kn$
increases at
$t=0$
. This is because a large
$Kn$
corresponds to a large most unstable wavelength, which has a small wavenumber and thus exhibits a small growth rate.
(
$a$
) Growth-suppression factor
$\psi$
as a function of
$k\epsilon$
for different
$Kn$
. (
$b$
) The growth rate
$\gamma$
with time
$t$
at different
$Kn$
. (
$c$
) The ratio
$\beta =\tau _D/\tau _I$
of the time scale of diffusion to the time scale of instability with time
$t$
at different
$Kn$
. (
$d$
) The maximum value
$a_{max}/\lambda _m$
of the dimensionless interface fluctuating amplitude versus
$\beta _{max}$
for different
$Kn$
.
Figure 7. Long description
Panel A: A line graph shows the growth-suppression factor as a function of kε for different Knudsen numbers (Kn). The x-axis is labeled kε and the y-axis is labeled ψ. Multiple lines represent different Kn values, with a black line indicating theoretical values by Duff et al., 1962. Panel B: A line graph displays the growth rate (γ) over time (t/τ) for different Kn values. The x-axis is labeled t/τ and the y-axis is labeled γ. Different colored lines represent various Kn values. Panel C: A scatter plot shows the ratio of the time scale of diffusion to the time scale of instability (β) over time (t/τ) for different Kn values. The x-axis is labeled t/τ and the y-axis is labeled β. Different symbols and colors represent various Kn values. Panel D: A scatter plot illustrates the maximum value of the dimensionless interface fluctuating amplitude (a_max/λm) versus the maximum value of β (β_max) for different Kn values. The x-axis is labeled β_max and the y-axis is labeled a_max/λm. Different symbols and colors represent various Kn values.
To determine whether instability overcomes diffusive damping, we examine the dimensionless time scale ratio
$\beta \equiv \tau _D / \tau _I$
, where
$\tau _D$
is the diffusive time scale and
$\tau _I$
is the instability time scale. When
$\beta \lt 1$
, the diffusion time scale is shorter than the instability time scale, indicating that mass transport is dominated by diffusion at the interface. Conversely,
$\beta \gt 1$
indicates that the instability dominates. The temporal evolution of
$\beta$
at different
$Kn$
is shown in figure 7(
$c$
). For
$Kn\lt 0.011$
, the value of
$\beta$
increases from less than 1 to greater than 1 over time, indicating that mass transport is dominated by a transition from diffusion to instability at the interface. For
$Kn \gtrsim 0.011$
,
$\beta$
remains below 1 throughout the evolution of the interface, indicating that mass transport is dominated by diffusion and does not transition to instability. This indicates that the critical Knudsen number for instability suppression in the present RT system is approximately
$Kn_{cr} = 0.011$
.
Figure 7(
$d$
) shows the dimensionless maximum value of the interface fluctuating amplitude (
$a_{max}/\lambda _m$
) as a function of the maximum time scale ratio
$\beta _{max}$
. Here,
$a_{max}/\lambda _m$
starts increasing when
$\beta _{max}$
becomes larger than 1. This again demonstrates that the proposed competition governs the emergence of instability.
3.5. Influence of Froude number on RTI
The previous subsections (§§ 3.1–3.4) examined the effect of the Knudsen number on the evolution of RTI at a fixed geometric ratio
$H/\lambda _m=0.5$
, showing that rarefaction suppresses the instability, with a critical threshold at
$Kn\approx 0.011$
. In practical applications, however, the domain height
$H$
significantly affects the strength of the driving force, with its effect often represented by
$Fr$
. We therefore extend the analysis to investigate the combined effects of
$Fr$
and
$Kn$
on the stability of the rarefied RT system.
Figure 8 presents the instability regime diagram in the
$(Kn, Fr)$
parameter space. The parameter space is partitioned into three distinct regimes: (i) no instability development, (ii) rarefaction-induced suppression, and (iii) unstable growth. When
$H/\lambda _m \lesssim 0.4$
, the RTI fails to develop beyond the linear stage in the continuum regime (Wei & Livescu Reference Wei and Livescu2012). In this regime, the constrained domain size prevents perturbations from transitioning to the nonlinear growth. For
$H/\lambda _m \gtrsim 0.4$
, perturbations can evolve into the nonlinear regime, and instability suppression is governed predominantly by the Knudsen number. As established in § 3.4, when
$Kn \gtrsim 0.011$
, the instability is completely suppressed, irrespective of the domain size or Froude number. This behaviour aligns with the time scale competition framework, wherein molecular diffusion remains faster than instability growth throughout the evolution (
$\beta = \tau _D/\tau _I \lt 1$
).
Instability regime diagram in the
$(Kn,Fr)$
parameter space for the rarefied RT system. The pink region denotes conditions where instability grows, while the blue region indicates suppression by rarefaction effects. The blue region represents cases where the domain height is insufficient to allow perturbations to evolve into the nonlinear stage before interacting with the boundaries. Solid squares indicate simulations where instability is observed, whereas hollow squares denote cases where no instability is detected. The green triangles highlight supplementary cases (cases 42–48) where the domain height
$H$
was increased and gravitational acceleration
$g$
was reduced proportionally to maintain a constant
$Fr$
, thereby verifying that the observed threshold is independent of absolute length scales.
To verify that this threshold is not an artefact of the specific domain configuration, additional simulations (cases 42–48) were performed and compared with cases 21–27 at fixed
$Fr$
. In these supplementary cases, indicated by the green triangles in figure 8, the gravitational acceleration was reduced by a factor of 5, while the domain height was increased proportionally. This approach maintained a constant
$Fr$
while varying the absolute length scales. The results confirm that the instability remains completely suppressed for
$Kn\gtrsim 0.011$
, consistent with the primary simulation set. This demonstrates that within the parameter range considered, the critical Knudsen number for instability suppression is independent of
$Fr$
, reinforcing that the stability threshold is determined by the competition between diffusive and instability time scales, rather than by geometric confinement effects.
4. Conclusions
In summary, we have systematically investigated the suppression of the Rayleigh–Taylor instability (RTI) in rarefied gases using GPU parallelised DSMC simulations across a broad range of Knudsen numbers. Our results establish
$Kn \approx 0.011$
as a critical threshold; for
$Kn \gtrsim 0.011$
, interfacial instability is completely suppressed, and the maximum perturbation amplitude remains constant while the concentration profile undergoes rapid diffusion-driven reversal.
The consequence of rarefaction operates across scales. At the molecular level, increasing
$Kn$
extends the duration of super-diffusive motion, as is evident from the molecular mean squared displacement, delaying the transition to normal diffusion. At the macroscopic scale, this enhances interfacial mixing, rapidly weakening the density gradient and reducing the effective Atwood number that drives the instability.
Furthermore, the underlying physics of instability emergence is quantitatively captured by a time scale competition framework. We introduce a diffusion time scale
$\tau _D$
and an instability time scale
$\tau _I$
. Instability emerges only when
$\tau _I \lt \tau _D$
(i.e.
$\beta = \tau _D/\tau _I \gt 1$
). For
$Kn \gtrsim 0.011$
,
$\tau _D$
remains shorter than
$\tau _I$
throughout the evolution (
$\beta \lt 1$
), confirming that diffusion perpetually dominates, and instability cannot develop – consistent with the observed saturation of the perturbation amplitude. Thus the transition to instability in rarefied RTI is governed by the competition between molecular diffusion and buoyant growth, leading to a critical
$Kn$
that dictates the stability of the system.
Finally, we have also identified two distinct regimes for RTI suppression: geometric confinement (
$H/\lambda _m\lesssim 0.4$
) prevents instability by restricting perturbation growth, while for larger domains (
$H/\lambda _m\gtrsim 0.4$
), suppression occurs when
$Kn\gtrsim 0.011$
regardless of
$Fr$
. Validation across different gravitational accelerations confirms that the critical Knudsen number is independent of
$Fr$
, reinforcing that instability onset is governed by the competition between diffusive and instability time scales.
This study establishes a framework for predicting flow stability in rarefied gases by linking molecular transport to macroscopic instability thresholds. The proposed mechanism can be extended to other instabilities, such as Richtmyer–Meshkov, Kelvin–Helmholtz and shear-flow transitions, offering a general approach to probe the microscopic origins of hydrodynamic instability. A long-term challenge is to construct a unified theory connecting discrete molecular dynamics to continuum instability phenomena across microscopic, mesoscopic and macroscopic scales.
Funding
This work is supported by the NSFC Excellence Research Group Programme for ‘Multiscale problems in nonlinear mechanics’ (no. 12588201). It was also supported by the Natural Science Foundation of China under grant nos 12372219, 52450304, 12432011, 12422208, 12421002 and 12372220.
Declaration of interests
The authors report no conflict of interest.
P=101325
T=273.15
ω
α
×
Pref
λ
k(a−a0)
t/τ
Kn=0.005
Kn=0.001
Kn=0.009
Kn
Fr
l
τ
αv
Kn
αv=0.5
αv
nh
n
αv=nh/(nh+nl)
t=4τ
a−c
t=6τ
d−f
t=8τ
g−i
t=10τ
j−l
a,d,g,j
b,e,h,k
c,f,i,l
Kn=0.001
a
a
Kn
a∼t0.5
a∼t2
b
amax
Kn
a
Kn=0.001
b
Kn=0.013
αv
t=τ
c
t/τ
Kn
⟨r2⟩∼t
⟨r2⟩∼t2
d
t/τ
Kn
a
Kn
t=5τ
⟨ρ⟩x
x
Kn=0.005
b
At,eff
Kn
a
ψ
kϵ
Kn
b
γ
t
Kn
c
β=τD/τI
t
Kn
d
amax/λm
βmax
Kn
(Kn,Fr)
H
g
Fr
