1. Introduction
Richtmyer–Meshkov instability (RMI) arises when a perturbed interface separating fluids with disparate physical properties is accelerated by a shock wave (Richtmyer Reference Richtmyer1960; Meshkov Reference Meshkov1969). It is generally considered as an impulsive version of Rayleigh–Taylor instability (RTI), which occurs only when a heavy fluid is accelerated by a light fluid (Rayleigh Reference Rayleigh1882; Taylor Reference Taylor1950); RMI is an important phenomenon occurring in many applications (Prestridge Reference Prestridge2018; Zhou et al. Reference Zhou, Clark, Clark, Gail Glendinning, Aaron Skinner, Huntington, Hurricane, Dimits and Remington2019). For example, in inertial confinement fusion reactions, the mixing induced by RMI results in fuel contamination and limits the fusion energy gain (Chu et al. Reference Chu, Wang, Qi, Xu and Li2022; Hurricane et al. Reference Hurricane, Patel, Betti, Froula, Regan, Slutz, Gomez and Sweeney2023). It is therefore essential to comprehend the evolution of RMI on the target interface (Zhou Reference Zhou2017; Zhai et al. Reference Zhai, Zou, Wu and Luo2018; Liang & Luo Reference Liang and Luo2023).
The RMI of a single-mode interface has been widely investigated theoretically. In the linear stage, the interface profile is dominated by the fundamental mode, which exists on the initial interface. As the interface evolves into nonlinear stages, higher-order harmonics with wavenumbers being integer multiples of the fundamental mode wavenumber are generated and grow. This process is referred to as the self-coupling of the fundamental mode. Richtmyer (Reference Richtmyer1960) first derived the linearized governing equation for the single-mode RMI growth. Subsequently, various models for linear growth (Meyer & Blewett Reference Meyer and Blewett1972; Yang, Zhang & Sharp Reference Yang, Zhang and Sharp1994; Wouchuk & Nishihara Reference Wouchuk and Nishihara1997; Wouchuk Reference Wouchuk2001) and nonlinear growth (Dimonte & Ramaprabhu Reference Dimonte and Ramaprabhu2010; Zhang & Guo Reference Zhang and Guo2016) have been established through diverse physical hypotheses. Recently, a unified theoretical model for spatiotemporal development of Richtmyer–Meshkov fingers was proposed by combining the classical potential flow theory with a dual-source model (Liu, Zhang & Xiao Reference Liu, Zhang and Xiao2023). In addition, by incorporating different physical mechanisms into the analytical models, the three-dimensional (3-D) effect (Chapman & Jacobs Reference Chapman and Jacobs2006; Luo et al. Reference Luo, Guan, Si, Zhai and Wang2016), the geometric effect (Luo et al. Reference Luo, Li, Ding, Zhai and Si2019; Ge et al. Reference Ge, Li, Zhang and Tian2022), the high-amplitude effect (Wang et al. Reference Wang, Wang, Zhai and Luo2023b; Dimonte et al. Reference Dimonte, Nagori, Ramaprabhu and Boureima2024), the Atwood number (defined as 
$A=(\rho _2-\rho _1)/(\rho _2+\rho _1)$, with 
$\rho _2$ and 
$\rho _1$ being the fluid densities at each side of the interface) effect (Chen et al. Reference Chen, Li, Zhang and Tian2019, Reference Chen, Xing, Wang, Zhai and Luo2023), the interface coupling effect (Liang & Luo Reference Liang and Luo2023; Schalles Reference Schalles2023) and the Mach number effect (Motl et al. Reference Motl, Oakley, Ranjan, Weber, Anderson and Bonazza2009) have been evaluated.
In practical situations, random multi-mode perturbations are commonly present on initial interfaces, which means that more than one fundamental mode exists on the initial interface. Therefore, besides self-coupling of each individual fundamental mode, the different fundamental modes will also couple together in the multi-mode RMI, which is referred to as the mutual coupling between two fundamental modes. To characterize the perturbation growth of the multi-mode RMI, several theoretical approaches were adapted and improved from single-mode situations (Vandenboomgaerde, Gauthier & Mügler Reference Vandenboomgaerde, Gauthier and Mügler2002). A single bubble potential model was developed by Layzer (Reference Layzer1955), and various models have been derived since then through the potential flow method with different potential functions and mathematical treatments. Starting from the potential flow method, the perturbation expansion technique was employed to model the modal evolution of the single-mode RMI by considering the interface perturbation as a superposition of various Fourier modes (Zhang & Sohn Reference Zhang and Sohn1997). Besides, Haan (Reference Haan1991) proposed a second-order solution for multi-mode classical RTI growth (the Haan model). The Haan model and its variations have been verified extensively through simulations and high energy density experiments (Remington et al. Reference Remington, Weber, Marinak, Haan, Kilkenny, Wallace and Dimonte1995; Ofer et al. Reference Ofer, Alon, Shvarts, McCrory and Verdon1996; Di Stefano et al. Reference Di Stefano, Rasmus, Doss, Flippo, Hager, Kline and Bradley2017; Elbaz & Shvarts Reference Elbaz and Shvarts2018).
In the models mentioned above, the linear growth of the amplitude was satisfied after the shock impact, while the starting point of theoretical prediction was chosen artificially. Note that an implicit initial condition for shock–interface interaction is that the perturbation growth rate is zero, not a finite value, immediately after the shock impact (Fraley Reference Fraley1986; Yang et al. Reference Yang, Zhang and Sharp1994). This means that the amplitude growth rate will experience an acceleration from zero to an asymptotic value after the shock impact, which is referred to as the startup process. In some linear compressible analytical models (Richtmyer Reference Richtmyer1960; Yang et al. Reference Yang, Zhang and Sharp1994; Wouchuk & Nishihara Reference Wouchuk and Nishihara1997; Wouchuk Reference Wouchuk2001; Cobos Campos & Wouchuk Reference Cobos Campos and Wouchuk2014, Reference Cobos Campos and Wouchuk2016), the authors have already discussed the startup process, and have identified two main contributions to the asymptotic growth rate: the vorticity deposited at the interface by shock refraction, and the acoustic interaction of the perturbed interface with the shock front. This interaction is responsible for the oscillatory evolution of the growth rate. However, in some linear compressible analytical models (Richtmyer Reference Richtmyer1960; Yang et al. Reference Yang, Zhang and Sharp1994), the startup process was solved numerically, and there are no analytical solutions. Although there are infinite series solutions to the startup process in the linear model proposed by Wouchuk & Nishihara (Reference Wouchuk and Nishihara1997), the form of the solutions is complicated. In addition, the startup process is rarely considered in previous nonlinear analytical models, many of which were derived based on the initial conditions of linear growth. Vandenboomgaerde et al. (Reference Vandenboomgaerde, Gauthier and Mügler2002) matched their nonlinear theory to the early-time compressible phase without considering the startup process, which may lead to an obvious shift between theoretical predictions and numerical data, as given in figure 2(a) in their work. As a result, the nonlinear evolution of the perturbation immediately after the shock impact cannot be fully described by the existing nonlinear models due to the ignorance of the startup process.
Recently, Zhang, Deng & Guo (Reference Zhang, Deng and Guo2018) proposed a quantitative theory for the single-mode RMI that covers the entire time domain from early to late times based on the two-point Padé approximation method. This is the first nonlinear model that considers the startup process, to our best knowledge. However, it should be noted that the matching technique employed by Zhang et al. (Reference Zhang, Deng and Guo2018) introduces the nonlinear trend from the beginning of the interface evolution, not after the startup process, which may be inappropriate for the potential flow method. The startup process was investigated in detail, respectively, for a light–heavy interface (Lombardini & Pullin Reference Lombardini and Pullin2009) and a heavy–light interface (Li et al. Reference Li, Chen, Zhai and Luo2024), and it was concluded that the startup time is in proportion to the perturbation wavelength. Therefore, for an initially multi-mode interface, the fundamental modes have diverse startup times since the wavelengths for the fundamental modes are different. In particular, for a heavy–light interface, there is a unique phenomenon, called phase inversion, that may occur during or after the startup process (Li et al. Reference Li, Chen, Zhai and Luo2024). The prediction of the phase inversion has never been achieved in previous modal models. Considering the presence of complex processes such as phase inversion in the early evolution of RMI, the startup process must be considered in the model. Also, the lower accuracy of the modal model will cause divergence for predicting the amplitude growths of the high-order harmonics. Considering the startup process will not only make the model more approximate to the physical assumptions of initial conditions of linear growth, but also improve the accuracy of the theoretical predictions. Nevertheless, it is challenging to establish a model that can characterize the diverse startup times of different fundamental modes to give a complete mode evolution after the shock impact, and demonstrate the mode-coupling effects on different harmonics. These concerns motivate the current work.
In this study, an asymptotic matching modal model is established for RMI. The startup process is incorporated into the model to provide a complete description of the mode evolution after the shock impact. Through considering the feedback from high-order harmonic to the third-order harmonic, the model accuracy is improved and the model divergence is prevented. In addition, the model can evaluate the mutual-coupling effect on the amplitude variations of high-order harmonics besides the ‘beat modes’. In the following, first the derivation of the model is provided. Then experiments on single- and dual-mode interfaces subject to a shock wave are conducted to verify the model. In experiments, both the light–heavy and heavy–light interfaces are involved. The amplitudes of different modes are extracted from the experiments and compared with the predictions from the model.
2. Theoretical derivation
2.1. Classical modal models
Modal models are a class of potential flow models that are capable of describing the Fourier mode evolution. The amplitude of 
$j$th-order interface mode is defined as
\begin{equation} a_j(t)=\lambda^{{-}1}\int \exp({-2{\rm i} j{\rm \pi} x/\lambda})\,a(x,t)\,\mathrm{d}\kern 0.06em x, \end{equation}
where 
$\lambda$ is the perturbation wavelength for one period, 
$k_j=2j{\rm \pi} /\lambda$ is the 
$j$th-order perturbation wavenumber, 
${\rm i}$ is the imaginary unit, and 
$a(x,t)$ represents the interface profile. Starting from the potential flow method, the perturbation expansion technique was employed to model the modal evolution of the single-mode RMI by considering the interface perturbation as a superposition of various Fourier modes (Zhang & Sohn Reference Zhang and Sohn1997). The first three orders of this model, hereafter referred to as the ZS model, can be expressed as
$$\begin{gather} v_1^{ZS}(t)=v_1-\tfrac{1}{8}k_1^2v_1^2[(4A^2+1)v_1t^2+2a_{10}t], \end{gather}$$
$$\begin{gather}v_{2}^{ZS}(t)={-}Ak_1v_1^2t+\tfrac{1}{6}k_1^3v_1^2(8A^3v_1^2t^3+3Aa_{10}^2t), \end{gather}$$
$$\begin{gather}v_{3}^{ZS}(t)=\tfrac{3}{8}k_1^2v_1^2[(4A^2-1)v_1t^2-2a_{10}t], \end{gather}$$
where 
$v_1$ and 
$a_{10}$ are the initial velocity and amplitude for the linear RMI growth, respectively. Note that the direction of the incident shock movement is considered to be the positive direction of coordinate 
$z$. Consequently, the symbols for the even-order terms are different from those used by Zhang & Sohn (Reference Zhang and Sohn1997). However, except for the variable 
$v_3^{ZS}(t)$, the variables 
$v_1^{ZS}(t)$ and 
$v_2^{ZS}(t)$ can receive the feedback from higher-order variables under the accuracy considered by Zhang & Sohn (Reference Zhang and Sohn1997). This may result in the divergence for predicting 
$v_3^{ZS}(t)$ under certain circumstances.
To characterize the perturbation growth of the multi-mode RMI, several theoretical approaches were adapted and improved from single-mode situations. For example, the ZS model was extended by Vandenboomgaerde et al. (Reference Vandenboomgaerde, Gauthier and Mügler2002) to predict the early nonlinear amplitude growth of multi-mode interfaces. Besides, Haan (Reference Haan1991) proposed a second-order solution for multi-mode classical RTI growth, which is referred to as the Haan model. The Haan model in the ordinary differential equation form for RTI growth rate can be expressed as
\begin{align} \ddot{a}_j(t) &= \gamma^2k_j\,a_j(t)+Ak_j\sum_{l}\left[\ddot{a}_{l}(t)\, a_{j-l}(t)\,(1-\hat{\boldsymbol{k}}_l\boldsymbol{\cdot}\hat{\boldsymbol{k}}_j) \vphantom{\left(\frac{1}{2}-\hat{\boldsymbol{k}}_l\boldsymbol{\cdot} \hat{\boldsymbol{k}}_j-\frac{1}{2}\hat{\boldsymbol{k}}_l\boldsymbol{\cdot} \hat{\boldsymbol{k}'}\right)}\right.\nonumber\\ &\quad \left.{}+\dot{a}_{l}(t)\,\dot{a}_{j-l}(t)\left(\frac{1}{2}-\hat{\boldsymbol{k}}_l\boldsymbol{\cdot} \hat{\boldsymbol{k}}_j-\frac{1}{2}\hat{\boldsymbol{k}}_l\boldsymbol{\cdot} \hat{\boldsymbol{k}'}\right)\right], \end{align}where
\begin{equation} \boldsymbol{k}'=\boldsymbol{k}_j-\boldsymbol{k}_l, \end{equation}
and 
$\gamma (k)=\sqrt {gkA}$ is the linear growth rate for the RTI perturbation, with 
$g$ the acceleration. Here, 
$\boldsymbol {k}$ is the wave vector for mode 
$k$, and 
$\hat {\boldsymbol {k}}=\boldsymbol {k}/k$ is the unit vector. The subscripts 
$j$ and 
$l$ denote the corresponding modes. The first term on the right-hand side represents the contribution from the linear growth of the mode 
$k_j$ itself, and the other terms represent the contribution from the second-order mode coupling between modes 
$k_l$ and 
$k_{j-l}$. Ofer et al. (Reference Ofer, Alon, Shvarts, McCrory and Verdon1996) solved (2.5) under the assumption that 
$g$ is a constant, and a solution for 
$a_j(t)$ up to second-order accuracy was given as
\begin{equation} a_j(t)=a_j^{lin}(t)+\frac{1}{2}\,Ak_j\left(\sum_{l}a_{l}^{lin}(t)\,a_{j+l}^{lin}(t)- \frac{1}{2}\sum_{l< j}a_{l}^{lin}(t)\,a_{j-l}^{lin}(t)\right), \end{equation}
where 
$a^{lin}_j(t)$ denotes the first-order amplitude of mode 
$k_j$, which is zero when mode 
$k_j$ is not a fundamental mode. The second and third terms represent the generation of mode 
$k_j$ from shorter-wavelength modes (coupling between modes 
$k_l$ and 
$k_{j+l}$) and longer-wavelength modes (coupling between modes 
$k_l$ and 
$k_{j-l}$), respectively. Note that (2.5)–(2.7) have second-order accuracy for both self-coupling and mutual coupling. If 
$j=2l$, then the third term in (2.7) reduces to the first term in (2.3). Although RMI and RTI share the same linearized governing equations, whether the acceleration 
$g$ equals zero or not leads to significant differences in the properties of solutions (Zhou Reference Zhou2017). The linearized governing equation for RMI gives a linear solution (Richtmyer Reference Richtmyer1960). Specifically, 
$a^{lin}_j(t)=a_{j0}\cosh (\gamma (k_j)\,t)$ establishes when 
$g\neq 0$ in the RTI situation, while 
$a^{lin}_j(t)=a_{j0}+v_j t$ stands when 
$g=0$ in the RMI situation. The initial conditions 
$a_{j0}$ and 
$v_j$ are the amplitude and velocity of the mode 
$k_j$ at the time when the perturbation growth enters the linear regime. Note that the derivation from (2.5) to (2.7) is not restricted to exponential growth, and (2.7) holds for RMI too. Therefore, the solutions for the mutual-coupling part of RMI can be simplified into (Remington et al. Reference Remington, Weber, Marinak, Haan, Kilkenny, Wallace and Dimonte1995)
\begin{equation} a_{j\pm l}(t)\approx{\mp} \tfrac{1}{2}(k_j\pm k_l) v_{j}v_{l}t^2,\quad A\approx1, \end{equation}
with 
$v_j$ and 
$v_l$ being the linear growth rates of modes 
$k_j$ and 
$k_l$, respectively. The modes 
$k_j \pm k_l$ are referred to as the ‘beat modes’ generated from the coupling between two modes 
$k_j$ and 
$k_l$. Equation (2.8) equals 0 when mode 
$k_j$ or 
$k_l$ is not a fundamental mode. The Haan model has concise expressions and can predict the generations of beat modes.
For the dual-mode RMI (Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020), the Haan model is reformulated through some mathematical treatments to obtain the weakly nonlinear solutions when the fundamental modes include 
$k_1$ and 
$k_2$:
$$\begin{gather} v^{wn}_1(t)=v_1+k_1A\left(\sqrt{2}v_1a_{20}+\tfrac{3}{2}v_1v_2t\right), \end{gather}$$
$$\begin{gather}v^{wn}_2(t)=v_2-k_1A(v_1a_{10}+v_1^2t), \end{gather}$$
$$\begin{gather}v^{wn}_3(t)={-}\tfrac{3}{8}k_2A[2v_1v_2t+v_1a_{20}+(1+\sqrt{2})v_2a_{10}], \end{gather}$$
$$\begin{gather}v^{wn}_4(t)={-}k_2A(v_2^2t+2v_2a_{20}). \end{gather}$$However, these weakly nonlinear solutions were derived on the premise that the exponential growth is satisfied. This may lead to an overestimation of the weakly nonlinear solutions. To suppress this overestimation, a nonlinear model proposed by Zhang & Guo (Reference Zhang and Guo2016) for single-mode RMI (the ZG model), was introduced. The ZG model was used in previous work (Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020) to suppress the nonlinear terms, because previous work (Liu et al. Reference Liu, Liang, Ding, Liu and Luo2018) has verified that the ZG model can give reasonable predictions under similar physical parameters, and can suppress the rapid growth of the mode amplitudes from the weakly nonlinear stages. However, the ZG model is not a solution to the complete perturbation problem, because it stems from a plausible closure of the solution expanded near the tip of the bubble, which may or may not provide an accurate approximation to the actual solution. In some cases, the ZG model does not work and needs to be modified. The modified ZG model (mZG model) can be expressed as
\begin{equation} v^{mZG}_i(t)=\frac{v^{wn}_i(t)}{1+\hat{a}k_i\,v^{wn}_i(t)\,t}, \end{equation}
where 
$i=1,2,3,4$, and
\begin{equation} \hat{a}=\frac{3}{4}\,\frac{(1+A)(3+A)}{3+A+\sqrt{2}(1+A)^{1/2}}\, \frac{4(3+A)+\sqrt{2}(9+A)(1+A)^{1/2}}{(3+A)^2+2\sqrt{2}(3-A)(1+A)^{1/2}}. \end{equation}The Haan and mZG models are both capable of giving reasonable predictions for the evolution of RMI at a light–heavy interface. Under second-order accuracy, the divergence of these two models in late stages may be inevitable. Although the ZG model was introduced into the mZG model to suppress the overestimation, when the initial velocity is negative, the denominator of (2.13) may approach 0, which will also cause divergence.
2.2. Asymptotic matching modal model
We will establish a modal model based on the ideal gas equation of state for predicting the RMI development using an asymptotic matching method. After the shock impact, the amplitude growth will experience a startup process, which can be treated as a boundary layer on the time axis caused by the contradiction between the initial conditions and the linear solutions. This kind of problem can be solved by the singular perturbation method (Vasil'Eva, Butuzov & Kalachev Reference Vasil'Eva, Butuzov and Kalachev1995). Based on this property, the general solution 
$v^{gs}(t)$ can be decomposed into three parts: the internal solution 
$v^{in}(t)$, the external solution 
$v^{ex}(t)$, and the uniform approximation 
$v^{uni}(t)$ (Lin & Segel Reference Lin and Segel1988), i.e.
\begin{equation} v^{gs}(t)=v^{in}(t)+v^{ex}(t)-v^{uni}(t). \end{equation}
It is clear that the uniform approximation part is the linear growth rate. For a light–heavy interface accelerated by a weak shock wave, the impulsive model proposed by Richtmyer (Reference Richtmyer1960) has been proven to provide a good prediction to the linear growth (Chen et al. Reference Chen, Xing, Wang, Zhai and Luo2023). For a heavy–light interface accelerated by a weak shock wave, the irrotational model proposed by Wouchuk and Nishihara (WN model) can give an excellent prediction for the linear growth rate (Wouchuk Reference Wouchuk2001; Li et al. Reference Li, Chen, Zhai and Luo2024). As a result, the impulsive model and the WN model are used here to predict the linear growth rates 
$v^{uni}(t)$ of the light–heavy and heavy–light interfaces, respectively. The impulsive model can be written as
\begin{equation} v^{imp}=ka_0^+A\,\Delta V, \end{equation}and the WN model can be expressed as
\begin{equation} v^{WN}=k a_0^-\left[\frac{1-A}{2}\left(1+\frac{V_{RW}}{V_{IS}}\right)(U_1-\Delta V)+ \frac{1+A}{2} \left(1-\frac{V_{TS}}{V_{IS}}\right)\Delta V\right], \end{equation}
where 
$a_0^-$ is the initial amplitude, 
$a_0^+=(1-{\Delta V}/{V_{IS}})a_0^-$ is the post-shock amplitude, 
$V_{IS}$, 
$V_{TS}$, 
$V_{RW}$ and 
$\Delta V$ are the velocities of the incident shock, transmitted shock, reflected wave and shocked interface under laboratory coordinates, respectively, and 
$U_1$ is the flow velocity behind the incident shock. Given the incident shock Mach number and the physical properties of the fluids, these quantities can be solved by one-dimensional (1-D) gas dynamics theory (Anderson Reference Anderson1990).
The internal solution should capture the rapid acceleration of the amplitude growth immediately after the shock impact. The linear theory (Yang et al. Reference Yang, Zhang and Sharp1994) or the startup process model (SP model) (Li et al. Reference Li, Chen, Zhai and Luo2024) can be used to calculate the internal solution. The linear theory can capture the global rapid acceleration and the local vibration caused by pressure perturbations, but it has no explicit expression and must be solved numerically. The SP model has an explicit expression and can characterize the overall growth trend during the startup process. When the incident shock is not strong, and the local vibration caused by pressure perturbations is less obvious, the SP model is therefore chosen for clarity (Li et al. Reference Li, Chen, Zhai and Luo2024). The SP model can be expressed as
\begin{equation} v^{in}(t)=\frac{2 v^{uni}(t)}{(1-A)\coth [k\times L_{RW}(t)]+(1+A)\coth [k\times L_{TS}(t)]}, \end{equation}
where 
$v^{uni}(t)$ results from either (2.16) or (2.17), depending upon the interface type, and 
$L_{RW}(t)$ and 
$L_{TS}(t)$ are the distances of the reflected wave and transmitted shock from the shocked interface, which can be calculated as
\begin{equation} \left.\begin{matrix} L_{RW}(t)= (V_{RW}+\Delta V)t,\\[2pt] L_{TS}(t)=(V_{TS}-\Delta V)t. \end{matrix}\right\} \end{equation}The external solution should describe the nonlinear evolution following the startup process, namely, after the amplitude reaches linear growth. At this stage, the shock is sufficiently apart from the interface and compressibility can be ignored (Zhang et al. Reference Zhang, Deng and Guo2018). Previously, the approach proposed by Haan (Reference Haan1991) was employed to investigate the nonlinear evolution of multi-mode RMI. However, this approach has only second-order accuracy (Haan Reference Haan1991; Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020). In the present work, a perturbation expansion method is used to obtain a series solution for the external solution. One advantage of this method is that the accuracy can be controlled to any desired degree. Since the external solution should primarily capture the characteristics of incompressible nonlinear evolution in the late stage, the potential flow approach can be applied. This is done by starting from the governing equations for incompressible, inviscid, and irrotational fluids (Layzer Reference Layzer1955; Zhang & Sohn Reference Zhang and Sohn1997):
$$\begin{gather} \nabla^2 \phi_{q}(t,x,z)=0, \end{gather}$$
$$\begin{gather}\frac{\partial a}{\partial t}+\frac{\partial \phi_{q}}{\partial x}\,\frac{\partial a}{\partial x}-\frac{\partial \phi_{q}}{\partial z}=0\quad \textrm{at}\ z=a, \end{gather}$$
$$\begin{gather}\sum_{q=1}^{2} ({-}1)^q\rho_q \left(\frac{\partial \phi_q}{\partial t}+\frac{1}{2}\left[\left(\frac{\partial \phi_q}{\partial x}\right)^2+\left(\frac{\partial \phi_q}{\partial z}\right)^2\right]\right)=f(t)\quad \textrm{at}\ z=a. \end{gather}$$
Here, 
$\phi _q$ represents the velocity potential in fluid 
$q$, and the velocity field in fluid 
$q$ is given by 
$\boldsymbol {v}_q=\boldsymbol {\nabla }\phi _q$. The subscripts 
$q=1,2$ denote the fluids separated by the interface, and 
$z=a(x,t)$ represents the interface profile at a given time 
$t$, with 
$x$ and 
$z$ referring to the spanwise and normal directions of the interface, respectively. The selection of the velocity potential functions is significant for the potential flow method. By using the perturbation expansion technique, the interface and velocity potentials can be expanded into (Zhang & Sohn Reference Zhang and Sohn1997)
$$\begin{gather} a(x,t)\doteq \sum_{n=1}^{N} a_n(t)\cos (nkx)=\sum_{n=1}^{N}\varepsilon^n \sum_{m=0}^{\lfloor {(n+1)}/{2}\rfloor }a_{n,n-2m}(t)\cos(n-2m)kx, \end{gather}$$
$$\begin{gather}\phi_q(x,z,t)\doteq\sum_{n=1}^{N}\varepsilon^n \sum_{m=0}^{\lfloor {(n+1)}/{2} \rfloor }\phi_{i,n,n-2m}(t)\exp({({-}1)^{q+1}(n-2m)kz})\cos(n-2m)kx, \end{gather}$$
where 
$\varepsilon$ is the auxiliary perturbation quantity used to construct the 
$n$th-order perturbation expansion equations with dimension 1. This 
$\varepsilon$ will be eliminated in the procedure and will not appear in the final expressions. In addition, The amplitude 
$a(x,t)$ in (2.23) is defined only at the interface, thus the far-field conditions do not exist. The potential function 
$\phi _q(x,z,t)$ in (2.24) is defined in the fluids on both sides of the interface. The separation of variables form of (2.24) is set to be 
$\exp ({(-1)^{q+1}(n-2m)kz})$ for the normal direction coordinate 
$z$, which ensures that the perturbation potential diminishes far from the interface. Here, 
$N$ is the expansion order of the perturbation expansion method, and the ZS model is established based on 
$N=4$. The initial profile of a multi-mode interface can be expressed as
\begin{equation} a^-(x)=a(x,0^-)=\sum_{j} a_{j}^-\cos{(jkx)}. \end{equation}
We specify the mode 
$k_j$ as the fundamental mode for clarity hereafter because in the subsequent derivation, if mode 
$k_j$ is not the fundamental mode, then the evolution results will degenerate to 0. The initial conditions at the nonlinear starting point for the interface can be expressed as
$$\begin{gather} a_{n,n-2m}(t=0)=\sum_{j} a_{j0}\delta_{n,j} \delta_{n-2m,j}, \end{gather}$$
$$\begin{gather}v_{n,n-2m}(t=0)=\sum_{j} v_{j} \delta_{n,j} \delta_{n-2m,j}, \end{gather}$$
where 
$\delta$ is the Kronecker symbol. The value of startup time 
$t=\tau$ when the perturbation evolution enters the linear growth regime can be calculated as (Lombardini & Pullin Reference Lombardini and Pullin2009; Li et al. Reference Li, Chen, Zhai and Luo2024)
\begin{equation} \tau=\frac{1}{k}\left(\frac{1-A}{V_{RW}+\Delta V}+\frac{1+A}{V_{TS}-\Delta V}\right). \end{equation}
Because 
$\tau$ relates to 
$k$, the startup times are different for fundamental modes with different wavenumbers. The treatment for different startup times will be discussed later. Here, 
$v_j$ is equal to the linear growth rate, and 
$a_{j0}$ in previous work has different criteria such as 
$(a_j^++a_j^-)/2$ (Zhang et al. Reference Zhang, Deng and Guo2018) or 0 (Niederhaus & Jacobs Reference Niederhaus and Jacobs2003). To approximate the amplitude at 
$t=\tau$ as much as possible, 
$a_{j0}=a_j^++v_j\tau /2$ is used in the present work since the perturbation growth at the early stage of the startup process is approximately quadratic (Lombardini & Pullin Reference Lombardini and Pullin2009). The diversity among various criteria is small when the initial amplitude is small, but may be significant when the initial amplitude is high.
Subsequently, each order of equations with a balanced order of 
$\varepsilon$ can be provided, and the solutions can be derived directly. However, the expressions for the direct solutions are very complicated. In addition, the fundamental modes have varying startup times due to their disparate wavenumbers, and the solutions for the development of the fundamental modes have disparate orders of accuracy. To overcome these difficulties, the external solution can be decomposed into two components, i.e. the self-coupling of each fundamental mode and the mutual coupling between different fundamental modes. Here we define the self-coupling of the fundamental mode 
$k_j$ as the nonlinear evolution of the mode 
$k_j$ itself and the generation of its integer-multiple modes of mode 
$k_j$, namely, modes 
$2k_j$, 
$3k_j$, and so on under the single mode case. Mutual coupling between two fundamental modes is the additional effect when the initial interface has more than one fundamental mode.
By considering 
$N=4$, Zhang & Sohn (Reference Zhang and Sohn1997) proposed the ZS model for single-mode RMI development. The expressions for the first three modes of self-coupling of mode 
$k_j$ (
$\kern 1.5pt j=1$) are given by (2.2)–(2.4). As described earlier, because the feedback from high-order harmonics to the third-order harmonic is ignored, the ZS model may diverge in predicting the third-order harmonic growth. To prevent this divergence, the equation is expanded to 
$N=5$ to include higher-order terms, enabling the feedback from the fifth-order harmonic to the third-order harmonic. In addition, because 
$ka_i \ll O(1)$, the high-order 
$ka_i$ terms can be reduced for each order of 
$t$, thereby the simplified solutions can be provided. The self-coupling parts of our new model are
$$\begin{gather} v_j^{{(s)}}(t)=v_j-\tfrac{1}{8}k_j^2v_j^2[(4A^2+1)v_jt^2+2a_{j0}t], \end{gather}$$
$$\begin{gather}v_{2j}^{{(s)}}(t)={-}Ak_jv_j^2t+\tfrac{4}{3}A^3k_j^3v_j^4t^3, \end{gather}$$
\begin{gather}\begin{aligned} v_{3j}^{{(s)}}(t) &={-}\tfrac{3}{4}k_j^2a_{j0}v_j^2t+\tfrac{3}{8}k_j^2v_j^3(4A^2-1)t^2 \nonumber\\ &\quad +\tfrac{1}{32}(15+32A^2)k_j^4a_{j0}v_j^4t^3+\tfrac{3}{128}(5+40A^2-144A^4)k_j^4v_j^5t^4.\end{aligned} \end{gather}
The superscript ‘
${{(s)}}$’ represents the self-coupling of a single fundamental mode. Note that our model is improved based on the ZS model. In the ZS model, only the leading term feedbacks of high-order harmonics to the first- and second-order harmonics were considered. In our model, similarly, we consider only the leading term feedbacks of high-order harmonics in the case of self-coupling, i.e. only the feedback from the third-order harmonic to the first-order harmonic is considered, and only the feedback from fifth-order harmonic to the third-order harmonic is considered. As a result, the feedback from the fifth-order harmonic is present in (2.31) but is absent in (2.29). Actually, the divergence of the nonlinear perturbation series for self-coupling is a serious issue, and a partial sum of this diverging series is not an accurate approximation to the solution. However, it is difficult to determine the number of high-order terms to achieve the best prediction. Taking 
$k_j$ (
$\kern 1.5pt j=1, 2$) as an example of fundamental modes combination to illustrate mutual coupling and eliminating the small quantities by scale comparison, the generations from mutual coupling are
\begin{align} v_{112}^{{(m)}}(t) &= k_1 Av_1v_2 t+\tfrac{1}{4} k_1^2 ({-}2 A^2 -3 ) v_1 v_2^2 t^2 \nonumber\\ &\quad -\tfrac{1}{3} k_1^3 (A^3 +A) v_1^3 v_2 t^3, \end{align}
\begin{align} v_{212}^{{(m)}}(t) &={-}4A^2v_1^2 v_2k_1^2 t^2+\tfrac{1}{3} (4 A^3 +5 A)k_1^3 v_1^2 v_2^2 t^3 \nonumber\\ &\quad +\tfrac{1}{24} (72 A^4 +33 A^2 -4) k_1^4 v_1^4 v_2 t^4, \end{align}
\begin{align} v_{312}^{{(m)}}(t) &={-}3 k_1 A v_1 v_2 t+\tfrac{3}{8} k_1^2[(4 A^2- 1) v_1^2+ ({-}12 A^2+ 1) v_2^2]v_1t^2 \nonumber\\ &\quad +\tfrac{1}{4}k_1^3 [(54 A^3-11 A)v_1^2+(48 A^3+9 A)v_2^2]v_1 v_2 t^3, \end{align}
\begin{align} v_{412}^{{(m)}}(t) &= k_1^2 (8 A^2-2) v_1^2v_2 t^2+\tfrac{1}{3}k_1^3 (96 A^3 -28 A) v_1^2 v_2^2 t^3 \nonumber\\ &\quad +\tfrac{1}{12} k_1^4 ({-}512 A^4+223 A^2+3)v_1^4 v_2 t^4. \end{align}
The superscript ‘
${{(m)}}$’ represents the mutual coupling between two different fundamental modes. The subscript ‘
$ijl$’ denotes that the velocity of the mode 
$k_i$ is generated from mutual coupling between the fundamental modes 
$k_j$ and 
$k_l$. Three terms for each mode considering the mutual coupling are reserved. Series solutions with more than one term can provide convenience for further mathematical treatments such as Padé approximations. The previous studies have shown that the beat mode with wavenumber 
$k_j \pm k_l$ appears when more than one mode is initially present (Haan Reference Haan1991; Remington et al. Reference Remington, Weber, Marinak, Haan, Kilkenny, Wallace and Dimonte1995; Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020). At lowest order of accuracy of time 
$O(t)$ and initial amplitude 
$O[(ka_i)^0]$, this new model can reduce to (2.7). It can be found from the new model that the other modes will also be affected by mutual coupling at higher order. This phenomenon will be verified by experiments hereinafter. By deriving a higher-accuracy model, the divergence of solutions is prevented until the profile of interface is multi-valued.
The complete form for our new model to predict the mode 
$k_i$ evolution can be written as
\begin{align}
v^{gm}_i=\sum_j\left\{\left[v^{in}_i(t)-v_i^{uni}\right]\delta_{i,
j}+v_{i}^{{(s)}}\left[(t-\tau_j)H(t-\tau_j)\right] + \sum_{l\neq
j}v_{ijl}^{{(m)}}\left[(t-\tau_{jl}^{(m)})H(t-\tau_{jl}^{(m)})\right]\right\},
\end{align}
where 
$H(t)$ is the Heaviside function, 
$\tau _i$ is the time scaling for the startup process for each fundamental mode 
$k_i$, and 
$\tau ^{{(m)}}_{jl} = (\tau _j+\tau _l)/2$ is the startup time of the mode generated from the mutual coupling between the fundamental modes 
$k_j$ and 
$k_l$. The first term in the sum is used to match the startup process with the nonlinear growth; the second term represents the self-coupling of fundamental mode 
$k_j$, and it exists only when 
$k_i$ is an integer multiple of 
$k_j$; the last term represents the mutual coupling between two different fundamental modes 
$k_j$ and 
$k_l$. Note that this asymptotic matching method is not restricted to this nonlinear model, but is also applicable to other nonlinear models on RMI.
Note that the new model is established assuming that the different fundamental modes have small initial perturbation amplitudes, i.e. 
$ka_i\ll O(1)$, because there are existing theoretical models for predicting the linear and nonlinear growth rates of the amplitude. If the shock hits the interface with high initial perturbation amplitudes, then complicated phenomena, such as unsteady shock refraction and Mach stems, arise, and there are no analytical solutions for the linear and nonlinear growths, to the best of the authors’ knowledge. In the present work, the relative phase difference of two different fundamental modes should be integer multiples of 
${\rm \pi}$, such that the dual-mode interface can be expressed in the form 
$z=a_{01}\cos (k_1x)\pm a_{02}\cos (k_2x)$. For high-frequency multi-mode initial conditions (with fundamental modes greater than 2), provided that the relative phase differences of the different fundamental modes are integer multiples of 
${\rm \pi}$, the different startup processes of the different fundamental modes can still be captured by the asymptotic matching method proposed, and (2.36) still holds. Besides, the perturbation expansion method for nonlinear evolution can still be applied although the specific expressions for the nonlinear part are different.
In addition, the new model should be effective for a wide range of Atwood numbers as it is derived directly from potential flow models without empirical parameters and has a series form of Atwood number, but it may lose validity at extreme Atwood numbers such as 
$|A|<0.1$ or 
$|A|>0.9$. For light–heavy interfaces, the model may be applicable only for low shock Mach numbers because the shock proximity effect will significantly flatten the bubble front when the incident shock is strong (Motl et al. Reference Motl, Oakley, Ranjan, Weber, Anderson and Bonazza2009). The change in the bubble morphology will alter the modal evolution, which affects the validity of the model. For heavy–light interfaces, since the transmitted shock wave will quickly move away from the interface due to the high acoustic velocity of the fluid on the transmitted side, the shock proximity effect is very weak, and the model is expected to be effective for relatively high Mach numbers.
To verify the model proposed, experiments are performed. Note that very limited data are available on the startup process in the literature. Moreover, RMI of a dual-mode heavy–light interface has never been considered previously. Consequently, both single-mode and dual-mode interfaces are considered, and both light–heavy and heavy–light interfaces are involved in our present experiments. In the following, experimental methods and results are first described and discussed, then comparison of experimental results with theoretical results is made.
3. Experimental methods
In experiments, five different interface profiles are designed, including three kinds of single-mode interface with different wavenumbers, and two kinds of dual-mode interface with different phase combinations. The profiles of the initial interfaces can be expressed as
\begin{equation} \left.\begin{array}{lll} \mathrm{SM\_lh} & z=a_{01}\cos(k_1 x), & x \in [0,120]\,\textrm{mm},\\ \mathrm{SM\_hl} & z={-}a_{01}\cos(k_1 x), & x \in [{-}30,90]\,\textrm{mm},\\ \mathrm{AP\_lh} & z=a_{01}\cos(k_1 x)-a_{02}\cos(k_2 x), & x \in [{-}90,30]\,\textrm{mm},\\ \mathrm{IP\_lh} & z=a_{01}\cos(k_1 x)+a_{02}\cos(k_2 x), & x \in [{-}60,60]\,\textrm{mm},\\ \mathrm{AP\_hl} & z={-}a_{01}\cos(k_1 x)+a_{02}\cos(k_2 x), & x \in [{-}30,90]\,\textrm{mm},\\ \mathrm{IP\_hl} & z={-}a_{01}\cos(k_1 x)-a_{02}\cos(k_2 x), & x \in [{-}60,60]\,\textrm{mm},\\ \end{array}\right\} \end{equation}
where 
$k_1$ and 
$k_2$ are wavenumbers of two fundamental modes, and 
$a_{01}$ and 
$a_{02}$ are initial amplitudes of two fundamental modes at the boundary plane. Note that the initial negative amplitude of the mode 
$k_1$ in dual-mode heavy–light cases is designed to maintain a positive velocity. The width of the flow field is 140 mm, consisting of perturbation with two to four periods and two straight segments on both sides to eliminate the boundary effects (Luo et al. Reference Luo, Liang, Si and Zhai2018). For each interface configuration, both air–SF
$_6$ and SF
$_6$–air interfaces are involved.
The soap-film technique (Liu et al. Reference Liu, Liang, Ding, Liu and Luo2018; Gao et al. Reference Gao, Guo, Zhai and Luo2024) is used to form the discontinuous interfaces in this work. The super-hydrophobic-oleophobic surface technique developed by Li et al. (Reference Li, Cao, He, Zhai and Luo2023) instead of filaments is used to constrain the soap-film interface. Due to the surface tension, the soap-film interface is a 3-D surface with a minimum-surface feature (Luo, Wang & Si Reference Luo, Wang and Si2013). The multi-mode 3-D surface can be expressed as
\begin{equation} z=\sum_{j}f_j(y)\cos(k_jx). \end{equation}
When only one fundamental mode exists, the interface reduces to a single-mode one, and the governing function for 
$f_j(y)$, according to the work of Luo et al. (Reference Luo, Wang and Si2013), is
\begin{equation} f_j''=k^2f_j(1+f_j'^2),\quad f_j'|_{y=0}=0,\quad f_j|_{y=0}=a^{sp}_{j}, \end{equation}
where 
$a^{sp}_{j}$ is the perturbation amplitude of the fundamental mode 
$k_j$ at the symmetry plane 
$y=0$. Note that (3.3a–c) is a homogeneous ordinary differential equation, and each 
$f_j(y)$ should satisfy (3.3a–c) for a multi-mode interface. For different fundamental modes in a multi-mode soap-film interface, the 3-D effect is different, therefore the 3-D correction of the perturbation amplitude is obligatory. The average amplitude of the fundamental mode 
$k_j$ is calculated by
\begin{equation} a_{{v}j}=\frac{\displaystyle\int_{{-}h/2}^{h/2}f_j(y)\,\textrm{d}y}{h}. \end{equation}The initial interface parameters, including the average amplitudes of the fundamental modes, are provided in table 1.
Initial interface parameters, where ‘_lh’ or ‘_hl’ indicate that the interfaces are light–heavy or heavy–light, respectively. For single-mode cases, the parameters for light–heavy and heavy–light interfaces are the same except for 
$\theta _{b}$ and 
$\theta _{s}$. For a light–heavy interface, 
$\theta _{b}=0$ and 
$\theta _{s}={\rm \pi}$. Also, ‘AP’ and ‘IP’ denote anti-phase and in-phase of two basic modes, 
$\lambda$ is the perturbation wavelength, 
$a_{0}$ is the perturbation amplitude on the constraint boundary, and 
$a_{v}$ is the integral average amplitude of the 3-D surface. The subscripts 1 and 2 denote the first and second basic modes on the initial interface, 
$a_{b}$, 
$a_{s}$, 
$a$ are the initial heights of the bubble, spike and amplitude, respectively, and 
$\theta _{b}$ and 
$\theta _{s}$ are the angular locations of bubble and spike. The unit of length is mm.
The experiments are conducted in a horizontal shock tube (Guo et al. Reference Guo, Cong, Si and Luo2022). For air–SF
$_6$ interfaces, the Mach number of the incident shock moving in air is 
$1.26 \pm 0.01$. For SF
$_6$–air interfaces, the Mach number of the incident shock moving in SF
$_6$ is 
$1.31 \pm 0.01$. The physical parameters of the background flow of ten cases are listed in table 2. The experimental measurements agree well with the theoretical calculations from 1-D theory. The post-shock flow is recorded by high-speed schlieren photography. The frame rate of the high-speed video camera (FASTCAM SA-Z, Photron Limited) is 50 000 frames per second, and the exposure time is 
$1/2\,880\,000$ s. The spatial resolution of the schlieren image is 
$0.26 \pm 0.01$ mm pixel
$^{-1}$. The ambient pressure and temperature are 
$101.3 \pm 0.1$ kPa and 
$297.5 \pm 0.7$ K, respectively.
Physical parameters of background flow: 
$Ma$ is the incident shock Mach number, 
$\psi$ is the volume fraction of SF
$_6$ on the heavy gas side of the interface, and 
$V_{IS}$, 
$\Delta V$ and 
$V_{TS}$ are the velocities of the incident shock, shocked interface and transmitted shock, respectively. The superscripts 
$e$ and 
$1\textrm {-}D$ denote the data obtained from experimental measurements and 1-D gas dynamics theory, respectively. The unit of velocity is m s
$^{-1}$.
4. Results and discussion
4.1. Qualitative results
Developments of wave patterns and interface morphologies of the shock-accelerated air–SF
$_6$ interfaces and SF
$_6$–air interfaces obtained from experiments are shown in figures 1 and 2, respectively. The time origin 
$t=0$ is defined as the moment when the transmitted shock wave leaves the interface. When the interface and waves are away from the initial position, the shadows of the initial interface are erased from experimental images for clarity.
Sequences of schlieren images of air–SF
$_6$ interface evolution and wave propagation for five cases: (a–c) single-mode interfaces, (d,e) dual-mode interfaces. Here, IS means incident shock, RS means reflected shock, TS means transmitted shock, and SI means shocked interface.
Sequences of schlieren images of SF
$_6$–air interface evolution and wave propagation for five cases: (a–c) single-mode interfaces, (d,e) dual-mode interfaces.
The developments of a single-mode air–SF
$_6$ interface and SF
$_6$–air interface accelerated by a planar shock have been widely investigated, and only a brief description is provided here. Taking case SM40 as an example, after the incident shock passes through the air–SF
$_6$ interface, as shown in figure 1, the perturbation amplitude begins to grow, and the spikes (heavier fluid penetrating lighter fluid) and bubbles (lighter fluid penetrating heavier fluid) arise quickly (
$t = 1127\,\mathrm {\mu }$s). At the late stages, vortices are generated on both sides of the spikes (
$t = 1667\,\mathrm {\mu }$s). For a shock-accelerated SF
$_6$–air interface, as shown in figure 2, after the shock impact, reflected rarefaction waves (RW) are clearly observed, and become wider in width (
$t= 86\unicode{x2013} 206\,\mathrm {\mu }$s). No additional reflected waves are generated benefiting from the new interface formation method. The shocked interface experiences a phase inversion, then develops continuously. Later, vortices also arise at the spike head (
$t = 1446\,\mathrm {\mu }$s).
For the dual-mode interfaces, we take case AP as an example to detail the process. The air–SF
$_6$ interface in figure 1 has a large spike at the centre and small ones located on both sides. After the shock wave accelerates the interface, the large spike rolls up with a pair of vortices formed at its neck, and bubbles develop with an obvious inclination to the large spike (
$t = 855\,\mathrm {\mu }$s). At late times, the difference between the sizes of the different spikes becomes more significant, but the interface profile remains relatively thin and distinct (
$t = 1535\, \mathrm {\mu }$s). The SF
$_6$–air interface in figure 2 has a large bubble at the centre and small ones located on both sides. Because different fundamental modes have different phase-inversion times, a moment when the amplitude reduces to zero does not exist. The mode with a short wavelength finishes phase inversion earlier, and the interface appears to be a single-mode one with 
$\lambda = 60$ mm (
$t = 86\,\mathrm {\mu }$s). Then the mode with a long wavelength finishes phase inversion, causing the dual-mode interface to resemble a single-mode interface with 
$\lambda = 30$ mm (
$t = 166\, \mathrm {\mu }$s). Later, both large and small bubbles appear on the interface (
$t = 526\,\mathrm {\mu }$s), and the small spikes develop with an apparent inclination to the large bubble (
$t = 946\,\mathrm {\mu }$s).
To perform modal analysis, in the previous work, Fourier analysis on the interface profile is made. However, for each 
$x$ pixel, we will choose a certain single 
$z$ pixel during the data processing since Fourier analysis is executed on a single-valued function, and only one period of the interface is involved. In the present work, the cosine series function is used to fit the experimental profiles using the nonlinear least squares method. This method allows us to retain all the dark pixels near the interface during the binarization process rather than choosing an individual pixel for each 
$x$ coordinate. In addition, this procedure can be conducted over the whole observation area rather than for a single period. Therefore, more information from the experimental images is utilized, and the human subjectivity is reduced. A remarkable advantage of the modal model is that the direct profile of the interface can be predicted, whereas the other models based on a selected linear start point cannot give this prediction. Figure 3 shows comparisons between the experimental schlieren images and the interface profiles extracted from mode decomposition and predicted from the theoretical model for different cases. From comparison, the modal model can predict the profile very well, although some small localized features cannot be predicted due to the ignorance of higher-order harmonics in the model.
Comparisons between experimental schlieren images and interface profiles extracted from mode decomposition (solid lines) and theoretical prediction (dotted lines).
4.2. Quantitative comparisons
To verify our model, comparisons of the amplitude variations between experimental measurements and theoretical predictions from the new model for single-mode cases are given in dimensionless form in figure 4. Note that as the small-amplitude asymptotic RMI growth rate is itself proportional to the initial ripple amplitude, by the time the startup process ends, we need to ensure that the linear theory is still applicable, i.e. the small-amplitude conditions should be satisfied. For light–heavy interfaces, perturbation amplitudes will continue to grow after the shock impact. The amplitude (
$a_i$) at the time of the startup process ending (
$t=\tau _i$) will be larger than the value at the initial time. For heavy–light interfaces, however, due to the presence of the phase-inversion process, the amplitude will first reduce and then grow in the opposite phase, which leads to 
$|a_i(t=\tau _i)|<|a_i(t=0)|$. In our present single-mode and dual-mode experiments, the dimensionless amplitude of each fundamental mode after the startup process 
$|ka_i(t=\tau _i )|<0.3$ holds. Under this condition, the high-amplitude effect is negligible for the linear and nonlinear growths of both light–heavy and heavy–light interfaces (Dimonte & Ramaprabhu Reference Dimonte and Ramaprabhu2010; Wang et al. Reference Wang, Wang, Zhai and Luo2023a), and linear theory is still valid to provide a linear growth rate. Therefore, the small-amplitude conditions hold for the nonlinear modal evolution in our experiments. The amplitudes of the first three harmonics are provided since the amplitudes of other higher-order harmonics are too small. The time and 
$n$th-order harmonic amplitude are normalized as 
$kv_{lin}t$ and 
$ka_n$, respectively, where 
$v_{lin}$ is the linear growth rate of the total perturbation amplitude, which can be calculated by
\begin{equation} v_{lin}=\sum_j [v_j\cos(j\theta_{s})- v_j\cos(j\theta_{b})]/2, \end{equation}
where 
$\theta _{b}$ and 
$\theta _{s}$ are the phase positions of the bubble and spike, respectively, as listed in table 1. Benefiting from the asymptotic matching model, the development of harmonic amplitude can be predicted from the moment when the shock impact is finished, namely 
$t=0$, rather than from an artificially selected startup time. Note that the ZS model cannot predict the startup process, and the start point of the prediction by the ZS model is given by our new model. As mentioned before, relative to the ZS model, the correction is made only to the third harmonic in the new model. Therefore, for the amplitude growths of the first and second harmonics, the new model provides predictions similar to those of the ZS model. One can observe some obvious differences regarding normalized 
$a_1$ between experiments and theoretical predictions in light–heavy cases. Previous works have verified that the ZS model slightly overestimates the evolution of 
$a_1$ in single-mode situations with Atwood number smaller than 0.75 (Liu et al. Reference Liu, Liang, Ding, Liu and Luo2018; Chen et al. Reference Chen, Xing, Wang, Zhai and Luo2023), generally because the amplitude growth rate predicted by the ZS model satisfies a 
$1/t^2$ law when the Atwood number is smaller than 0.75, but the late-time 
$1/t$ law of the amplitude growth rate is expected from the potential flow model, and has been validated in previous numerical and experimental works (Dimonte & Ramaprabhu Reference Dimonte and Ramaprabhu2010; Mansoor et al. Reference Mansoor, Dalton, Martinez, Desjardins, Charonko and Prestridge2020). The early overestimation results in the divergence of the ZS model in predicting 
$a_1$. In figure 4, the predictions to the amplitude growth of the third harmonic from the ZS model and the new model are provided for comparison. For an air–SF
$_6$ interface, the predictions from the ZS model and the new model almost coincide except in the late stages. For SF
$_6$–air interfaces, the phase-inversion process can be well predicted by the model. For the third-order harmonic, the prediction from the ZS model coincides with the prediction from the new model before 
$kv_{lin}t=0.8$. Afterwards, the prediction from the ZS model diverges, whereas the new model can still give a reasonable prediction.
Comparisons of the normalized amplitude developments between experimental measurements and theoretical predictions for single-mode cases. The pink zone indicates the startup process of the fundamental mode.
Note that in the derived model, the divergence is well inhibited in the late stages by considering the fifth-order feedback to the third-order harmonic evolution. It is also possible for considering higher-order feedback to the higher-order (larger than 3) harmonic evolution, but it may not be meaningful. First and most importantly, the present modal model can already reasonably describe the profile evolution of the interface before the interface enters a non-single-valued period, and it is not necessary to consider more high-order harmonics. Second, the mathematical treatment becomes quite complicated if higher-order feedback to the higher-order (larger than 3) harmonic evolution is considered without obviously improving the accuracy. For further analysis of multi-mode perturbation fingers, it may be more effective to apply Padé approximations than higher-order expansion.
Comparisons of the amplitudes of the first four harmonics between experimental measurements and theoretical predictions for dual-mode cases are given in figure 5. The dimensionless time is shorter than that in single-mode cases because the interface enters a multi-valued period earlier. Here, the predictions from the Haan model (Haan Reference Haan1991) and from the mZG model (Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020) are provided as references, and they start from the end of the startup process, i.e. from 
$t=\tau$, whereas the predictions from the new model start from 
$t=0$. Note that the same asymptotic method is applied for our model as for the previous models. For air–SF
$_6$ interfaces, the predictions for the first harmonic from the new model are less accurate than those from the mZG model, which is probably ascribed to the overestimation of the ZS model to the first harmonic. Nevertheless, the new model generally gives good predictions to the amplitude developments of first four harmonics, especially for the high-order harmonics. For SF
$_6$–air interfaces, the new model provides better predictions to the amplitude developments of the first four harmonics than the other two models. Specifically, the different phase-inversion times for different fundamental modes are also well captured by the new model, which is never considered in the previous models, to the best of our knowledge. Note that the fourth-order harmonic is not the beat mode, therefore the Haan model cannot capture the mutual-coupling effect on the amplitude evolution of the fourth-order harmonic. In the Haan model, the fourth-order harmonic comes from the self-coupling of the fundamental second mode. Consequently, the Haan model gives the same results for the amplitude growths of the fourth-order harmonic in AP and IP cases. However, the amplitude growths of the fourth-order harmonic in AP and IP cases are obviously different in experiments. This fact shows that the mutual-coupling effect cannot be ignored when predicting the amplitude growth of the fourth-order harmonic, and the new model is able to quantify this effect.
Comparisons of the amplitude developments between experimental measurements and theoretical predictions for dual-mode interfaces. The coloured zones indicate the startup processes of the fundamental modes.
As for the mZG model (Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020), the amplitude growth rates of the third- and fourth-order harmonics can be calculated using (2.13) and (2.14). Here, when the fundamental modes are 
$k_j$ (
$\kern 1.5pt j=1,2$), the third- and fourth-order expressions for the weakly nonlinear solutions to initiate the mZG model are given by (2.11) and (2.12). Note that the mZG model provides different initial velocities for the amplitude growths of the third- and fourth-order harmonics compared to the early evolution trends of experimental results. For example, when the two fundamental modes of the dual-mode heavy–light interface are in phase, the initial velocity calculated by (2.11) is negative, while the initial velocities predicted by the Haan model and our new model are both zero. This is because the weakly nonlinear solution for the mZG model is based on the RTI frame, which leads to an artificial initial velocity for higher-order harmonics generated from mutual coupling. However, the linearized governing equation for the RMI problem is a homogeneous ordinary differential equation, which means that the RMI problem within the linear regime satisfies the linear superposition principle. Namely, the higher-order harmonics generated from mutual coupling do not own their initial velocities. This is also consistent with the experimental results and theoretical predictions from the Haan model and the new model.
5. Conclusions
For a shock-accelerated single-mode interface, previous works have shown that there is a startup process before the amplitude growth enters the linear and nonlinear regimes. In particular, there is a phase-inversion process for a shock-accelerated heavy–light interface. Previous models were generally established based on an artificially selected initial conditions to predict the linear and nonlinear growth rates. The startup process and the phase-inversion process cannot be described by the models. If multiple fundamental modes are present on the initial interface, then the startup times and phase-inversion times are distinct for different fundamental modes, which bring about more challenges for the model to predict each mode evolution.
In this work, an asymptotic matching modal model was established for predicting mode evolution in single- and dual-mode interfaces accelerated by a shock wave. The key to establishing the new model is the treatment of the startup process. The startup process can be treated as a boundary layer on the time axis caused by the contradiction between the initial conditions and the linear solutions. This kind of problem can be solved by the singular perturbation method. By decomposing the general solution into three parts – i.e. the internal solution, external solution and the uniform approximation – a modal model that can characterize the complete evolution of the perturbed interface after the shock impact is established. In the new model, the diverse startup times of the fundamental modes and higher-order harmonics are considered by asymptotically matching the nonlinear evolution of each mode separately. Relative to the previous modal model proposed by Haan (Reference Haan1991), in which only one term is reserved, three terms for each mode generated from mode coupling are reserved in the new model, which provides a convenience for further mathematical treatments. To the best of the authors’ knowledge, the new model can characterize the complete evolution of various-order modes for single- and multi-mode perturbations, whereas the previous models cannot. In addition, the asymptotic matching method established in this work is not restricted to the new model but is also applicable to other nonlinear models.
To validate the model, experiments on the developments of single- and dual-mode interfaces accelerated by a shock wave are conducted, and both light–heavy and heavy–light interfaces are involved. To obtain the interface profiles, the cosine series function is used to fit the experimental profile using the nonlinear least squares method. This method allows for utilizing more information from the experimental images, and reduces the human subjectivity. The interface profiles extracted from mode decomposition and predicted by the new model show high consistency with the experimental schlieren images.
For the shocked single-mode interface evolution, the ZS model (Zhang & Sohn Reference Zhang and Sohn1997) predicts the amplitude growth from the end of the startup process. The ZS model contains higher-order feedback only to the first- and second-order harmonics, which leads to the divergence in predicting the third-order harmonic evolution in the late stages. This divergence is well inhibited by considering the fifth-order feedback to the third-order harmonic evolution in the new model. For the shocked dual-mode interface evolution, the new model provides thorough solutions and characterizes the diverse startup times and phase inversion times for different fundamental modes for the first time. Relative to the modal model proposed by Haan (Reference Haan1991), the new model has a higher-order accuracy, and can evaluate the effect of mutual coupling between the fundamental modes on the amplitude variations of higher-order harmonics besides the ‘beat modes’. Relative to the mZG model (Luo et al. Reference Luo, Liu, Liang, Ding and Wen2020), the new model provides better predictions to the amplitude growths of the first four harmonics for the heavy–light interfaces, because the mZG model provides inappropriate velocities to initiate the developments of the higher-order harmonics, and gives completely different growth trends due to the presence of the phase-inversion process.
However, the present experiments were conducted only with weakly incident Mach number 1.2–1.3 and specific gas combinations of air and SF
$_6$ to validate the new model. In future work, more experiments, including different shock Mach numbers, different gas combinations, and more fundamental modes (greater than two) existing at the initial interface, will be conducted to further verify the new model.
Funding
This work was supported by the National Natural Science Foundation of China (nos 12372281 and 12388101), the Strategic Priority Research Program of Chinese Academy of Sciences (no. XDB0620201) and the Youth Innovation Promotion Association CAS.
Declaration of interests
The authors report no conflict of interest.
