Richtmyer–Meshkov instability when a shock is reflected for fluids with arbitrary equation of state

Abstract

First predicted by Richtmyer in 1960 and experimentally confirmed by Meshkov in 1969, the Richtmyer–Meshkov instability (RMI) is crucial in fields such as physics, astrophysics, inertial confinement fusion and high-energy-density physics. These disciplines often deal with strong shocks moving through condensed materials or high-pressure plasmas that exhibit non-ideal equations of state (EoS), thus requiring theoretical models with realistic fluid EoS for accurate RMI simulations. Approximate formulae for asymptotic growth rates, like those proposed by Richtmyer, are helpful but rely on heuristic prescriptions for compressible materials. These prescriptions can sometimes approximate the RMI growth rate well, but their accuracy remains uncertain without exact solutions, as the fully compressible RMI growth rate is influenced by both vorticity deposited during shock refraction and multiple sonic wave refractions. This study advances previous work by presenting an analytic, fully compressible theory of RMI for reflected shocks with arbitrary EoS. It compares theoretical predictions with heuristic prescriptions using ideal gas, van der Waals gas and three-term constitutive equations for simple metals, the latter being analysed with detailed and simplified ideal-gas-like EoS. We additionally offer an alternative explicit approximate formula for the asymptotic growth rate. The comprehensive model also incorporates the effects of constant-amplitude acoustic waves at the interface, associated with the D'yakov–Kontorovich instability in shocks.

1. Introduction

The Richtmyer–Meshkov instability (RMI) is a fundamental interfacial instability in hydrodynamics. It develops whenever a planar shock wave interacts with a material interface, separating two fluids with distinct thermodynamic properties. First predicted through numerical modelling by Richtmyer (1960) and subsequently experimentally observed by Meshkov (1969), this phenomenon holds paramount significance across various realms of physics and astrophysics (see Inogamov 1999; Brouillette 2002; Aglitskiy et al. 2010; Abarzhi, Gauthier & Sreenivasan 2013; Clavin & Searby 2016; Zhou 2017a,b; Zhou et al. 2021; Zhou 2024; Zhou, Sadler & Hurricane 2024) and references therein. In the domain of inertial confinement fusion and the associated area of high-energy-density physics (HEDP), the strong shocks propagate through condensed materials or high-pressure plasmas whose equations of state (EoS) under extreme conditions are typically non-ideal. Consequently, any theoretical framework aspiring to model HEDP conditions and replicate the RMI development must incorporate realistic fluid EoS.

Studies of interfacial instabilities begin with the small-amplitude, linear phase of their development, when the interfacial perturbations can be decomposed into Fourier series, and each Fourier component is dealt with separately. When the unperturbed state is time independent, which is the case for the Kelvin–Helmholtz (von Helmholtz 1868; Thomson 1871) and the Rayleigh–Taylor (Rayleigh 1882; Taylor 1950) interfacial instabilities, it is possible to use the normal-mode decomposition to determine the exponential growth rate. The RMI analysis is more challenging because the unperturbed flow is not steady (it involves a transmitted shock wave and a reflected shock or rarefaction wave propagating away from the shocked material interface), and, therefore, ‘one cannot separate the time dependence from the linearised perturbation equations’ (Brouillette 2002). This is why Richtmyer made his discovery by solving these equations numerically with the most powerful computer of the time, the MANIAC (mathematical analyser, numerical integrator and automatic computer) at Los Alamos. His linear analysis ‘places no restriction on the EoS, provided only that the values of several constants are known. But the only calculations performed with this code so far were for $\gamma$-law gases, initially cold, with the same value of $\gamma$’ (Richtmyer 1960). He was far ahead of his time: a code numerically solving the linear RMI perturbation equations was recreated only 40 years after Richtmyer's work by Yang, Zhang & Sharp (1994). It could handle arbitrary EoS, although its published output (see also Holmes et al. (1999)) still refers to ideal gases, generalising Richtmyer's work in several other ways (the case of reflected rarefaction wave, finite incident shock strength, different $\gamma$s.) The same applies to the analytical and semianalytical solutions of the linear RMI problem, starting from those published by Fraley (1986) and Velikovich (1996) for the cases of reflected shock and reflected rarefaction, respectively.

Given the relative complexity of finding and using the exact solutions of the linear RMI perturbation problem, it is helpful, particularly for theoretical analysis (see Mikaelian (1991, 1995), and references therein) to have an approximate but explicit formula for evaluating the constant asymptotic growth rate. Such a formula, also proposed by Richtmyer, is derived by treating the RMI as a limiting case of the Rayleigh–Taylor instability, with instantaneous, rather than constant, acceleration: $g=U_2\delta (t)$, where $U_2$ is the velocity instantly acquired by the shocked material interface, see figure 1, and $\delta (t)$ is the Dirac delta function. Then, the asymptotic growth rate is estimated as $\delta u_i^\infty =k \psi _i^0 U_2 \mathcal {A}_2$, where $\psi _i^0$ is the interfacial ripple amplitude, $\mathcal {A}_2$ is the Atwood number at the interface and $k$ is the interfacial ripple wavenumber. For an incompressible RMI flow triggered with a $\delta$-function acceleration pulse, this is indeed an exact result. Predictions of the linear compressible theory converge to it in the limits of low compressibility and/or weak shock, see the discussion by Velikovich, Herrmann & Abarzhi (2014). However, suppose the compressibility of the materials separated by the shocked interface is relevant, which is the case in most strong-shock RMI situations. In that case, the preshock values of the Atwood number $\mathcal {A}_0$ and the ripple amplitude, $\psi _i^{0-}$, differ from the respective postshock values, $\mathcal {A}_2$ and $\psi _i^{0}$. Then, the governing incompressible equations do not help determine which values to substitute into the above formula, and the choice has to be made heuristically, resulting in so-called prescriptions. Richtmyer, comparing the above formula with his numerical results, formulated the prescription involving the postshock values. Meyer & Blewett (1972) similarly developed a slightly different prescription, featuring the average of preshock and postshock interfacial ripple amplitudes and the postshock Atwood number. Vandenboomgaerde, Mügler & Gauthier (1998) suggested yet another prescription, averaging the preshock and postshock products of the Atwood number and the ripple amplitude. Of course, in the incompressible and weak-shock limits all prescriptions converge to the same, exact value of the RMI growth rate.

Figure 1. Sketch of the RMI when a shock is reflected. Panels (a,b) correspond to the preshock and postshock interaction conditions, respectively: (a) $t<0$; (b) $t>0$.

Sometimes, the heuristic prescriptions provide an excellent approximation to the asymptotic RMI growth rate; sometimes, the approximation is not good at all, and, unfortunately, it is impossible to predict which way it is without solving the perturbation problem exactly. The physical reason for this uncertainty is that the fully compressible RMI growth rate is not settled at the instant when the rippled interface is shocked. As demonstrated by Wouchuk & Nishihara (1997), the RMI growth rate is determined by two comparable contributions to the interfacial vorticity. One is the vorticity deposited at the instant of incident shock refraction. The other comes from the multiple refractions of sonic waves emitted by transmitted and reflected ripple shock fronts (or only transmitted shock front, in the case of reflected rarefaction). Both shock fronts are in causal interaction with the material interface, which plays the role of a rippled piston driving them both, and, through the piston, with each other. The sonic signals received and emitted by the interface become gradually weaker as the distance between the outgoing shock fronts and the interface increases, and this contribution to the RMI growth rates asymptotically approaches a constant value.

Due to the essential role of compressibility, this contribution cannot be captured by any of the prescriptions; one needs to solve the fully compressible perturbation problem. But when the EoS of shocked materials are not realistic, then the applicability of such analysis is still uncertain. Ideal-gas EoS, with $\gamma$s and the preshock pressure $p_0$ chosen as fitting parameters to match the one-dimensional (1-D) dynamics of observed shock flows, are often used for the compressible stability analysis of HEDP experiments with solid targets. For example, Ishizaki & Nishihara (1997) selected $\gamma =3$ for shock-compressed plastic and a fitting preshock pressure of $p_0=0.7$ Mbar to theoretically describe the dynamics of the observed rippled-shock oscillations in the experiments by Endo et al. (1995). Similarly, Holmes et al. (1999) selected the values of $\gamma =1.45$ and $1.8$ for low-density foam and beryllium, respectively, with the preshock pressure $p_0=0.1$ Mbar roughly corresponded to the x-ray preheat for their analysis of the RMI experimental data by Dimonte & Remington (1993) and Dimonte et al. (1996). Both Ishizaki & Nishihara (1997) and Holmes et al. (1999) successfully reproduced the results of their respective experiments during the small-amplitude stage of the perturbation development. However, since matching of the 1-D parameters does not necessarily translate into matching of the parameters determining the frequency and decay of the shock-fronts’ oscillations, which were observed by Endo et al. (1995) and affected the RMI growth measured by Dimonte & Remington (1993) and Dimonte et al. (1996), one cannot evaluate the accuracy of the compressible linear theory. This uncertainty is similar to that associated with the use of prescriptions. To eliminate it, the fully compressible linear RMI theory must be constructed for the realistic EoS, as envisioned by Richtmyer (1960).

Numerical analysis of this problem was published by Ward & Pullin (2011) for Mie–Grüneisen EoS; for details, see Ward (2011), and also by Tahir et al. (2011) for the case of aluminium and copper. Here, we report the analytic, fully compressible theory of the RMI for the reflected shock case with an arbitrary EoS, advancing the results of the earlier work by Wouchuk (2001), Campos & Wouchuk (2014) and Cobos-Campos & Wouchuk (2016). As examples of realistic EoS, we chose the venerable van der Waals (vdW) EoS for a non-ideal gas and the three-term, Grüneisen-type three-term EoS (Zel'dovich & Raizer 2002, Chapter XI, § 6), applicable to simple metals at pressures not exceeding $5$ Mbar. Of course, these results are readily modified for any EoS presented in analytical or tabular form. We demonstrate the specifics of the realistic-EoS, fully compressible RMI analysis, discuss various limiting cases and compare the theoretical predictions of the asymptotic growth rates with various prescriptions. The lack of Riemann problem similarity for different EoS (Quartapelle et al. 2003) is addressed by finding the best fit for both EoS zero-order profiles. With this approach, the incident shock strength and the effective $\gamma$ values are selected to yield closely matching values for the target parameters: shock density and pressure jumps, postshock Atwood number and postshock dimensionless contact interface amplitude.

In this work, in addition to the extension of the formulation to be applicable to arbitrary EoS, we use a different resolution method compared with previous theoretical studies (Wouchuk 2001; Campos & Wouchuk 2014; Cobos-Campos & Wouchuk 2016). Specifically, we address the initial-value problem by employing the direct inverse Laplace transform (ILT), which allows the final solution to be expressed as a closed integral. This approach overcomes issues related to the slow convergence of Fourier–Bessel coefficients and provides deeper insights into shock-to-shock acoustic coupling interactions, as demonstrated in Napieralski et al. (2024) for RMI developed in shock-flame interactions. Understanding the singularities in the complex plane is crucial for the ILT of the pressure field. For instance, the presence of non-removable singularities, as noted by Calvo-Rivera, Velikovich & Huete (2023), indicates the possibility of spontaneous acoustic emission (SAE) for one of the shocks. This allows the model to be extended to conditions where the interface is continuously excited by pressure disturbances.

The remainder of this paper is structured as follows. The problem formulation is presented in § 2 for the case of reflected and transmitted shocks. The following § 3 is devoted to the description of the pressure field comprised by the two postinteraction shocks. The analysis includes the possibility of SAE in the transmitted shock. Section 4 describes the velocity field split into acoustic and rotational perturbations. The asymptotic velocity field is also analytically derived. The RMI evolution of the interface is addressed in § 5 for different EoS and shock conditions. Results are compared with known prescribed formulae for the asymptotic growth rate and with $\gamma$-law EoS for metals. Lastly, conclusions are given in § 6. The resolution of the pressure field using hyperbolic coordinate transformation, followed by the Laplace transform and its inverse, is thoroughly detailed in Appendix A. The base-flow resolution particularised for vdW EoS, aluminium and copper is offered in the supplementary material available at https://doi.org/10.1017/jfm.2024.1010. Mathematica codes are also provided for both the base-flow properties and the perturbed RMI problem.

2. Problem formulation in Cartesian coordinates

2.1. Base-flow equilibrium condition for arbitrary EoS

The perturbation-free problem is readily given by the conservation equations across the shocks and the mechanical equilibrium condition at the interface. Because of there are five distinct flow regions, it is necessary to choose a clear notation criterion to facilitate their identification. In our case, as described in figure 1, numbers $0$, $1$ and $2$ are used to identify the preshock fluids, the flow behind the incident shock and the flow behind the transmitted and reflected shocks, respectively. An additional subindex is necessary to identify the position relative to the contact interface: $r$ and $t$ for the regions along where the reflected and transmitted shocks travel. The mechanical equilibrium applies at the interface. Before the interaction with the incident shock the equilibrium is dictated by the continuity of pressure $p_{0r}=p_{0t}$, provided that both fluids move at the same speed, i.e. null in the preshock materials’ reference frame.

In order to calculate the base flow equilibrium it is necessary to solve the conservation equations. In the absence of disturbances, the different regions separated by each shock wave and the contact interface are steady and uniform. Across each shock the flow variables changes are related to the Rankine–Hugoniot (RH) equations, namely

(2.1a––c)\begin{equation} \left[\rho u\right]=0,\quad \left[p+\rho u^{2}\right]=0,\quad \left[E+\dfrac{p}{\rho}+\dfrac{1}{2} u^{2}\right]=0, \end{equation}

for mass, momentum and energy, respectively; where the square brackets $[\ldots ]$ denote the difference between the indicated quantities’ values ahead and behind the shock front. In (2.1a–c) $\rho$, $u$, $p$ and $E$ mean, respectively, density, velocity, pressure and internal energy evaluated in front and behind each shock wave, with the material velocity measured in the corresponding shock reference frame, so that $U=u+{\rm const}$., in relation to the material velocity indicated in figure 1.

The zero-order problem formulation is completed after providing the equilibrium conditions that must be imposed at the planar interface: $p_{2t}=p_{2r}$ and $U_2(x<0)=U_2(x>0)$. Furthermore, along with the isobaric condition, an additional parameter is needed to relate the two initial thermodynamic states of the fluids at both sides of the interface. Typically, the parameter employed is the preinteraction Atwood number

(2.2)\begin{equation} \mathcal{A}_0=\dfrac{\rho_{0t}-\rho_{0r}}{\rho_{0t}+\rho_{0r}} \end{equation}

that is bounded $-1<\mathcal {A}_0<1$. The present study does not address the case when a rarefaction is reflected back after the incident shock refraction, which might manifest itself when $\mathcal {A}_0$ takes low (likely negative) values. However, the kind of the wave reflected (shock or rarefaction) cannot be predicted solely based on this parameter unless certain conditions are met, namely when the two materials share identical thermodynamic coefficients. So, only under these circumstances, a negative value for $\mathcal {A}_0$ ensures the occurrence of a reflected rarefaction. For an interface separating two ideal gases with similar adiabatic indices, the density jump $\rho _{0t}/\rho _{0r}$ is the only parameter that determines the problem configuration, rendering a shock-reflected problem whenever $\rho _{0t}/\rho _{0r}>1$. For the ideal-gas shock-interface interaction, the case of reflected rarefaction occurs if (5a)–(5b) in Yang et al. (1994) are satisfied. This requirement can be presented as a single inequality

(2.3)\begin{equation} \dfrac{1}{\gamma_r+\gamma_t+2}<\dfrac{\mathcal{A}_0}{\gamma_r-\gamma_t}<\dfrac{1}{\gamma_r+\gamma_t}, \end{equation}

as derived in Velikovich (1996), which is applicable to both $\gamma _r<\gamma _t$ and $\gamma _r>\gamma _t$ cases. In the present case, the boundary between the shock-reflected and rarefaction-reflected cases, typically referred as ‘total transmission’ since there are no waves reflected back, must be computed numerically for each specific EoS and flow conditions, as detailed in the supplementary material.

In dimensionless form, the pressure continuity at the contact interface, $p_{2t}=p_{2r}$, can be written as a function of the shock pressure jumps

(2.4)\begin{equation} \mathcal{P}_r\mathcal{P}_i=\mathcal{P}_t, \end{equation}

where $\mathcal {P}_r=p_{2r}/p_{1r}$, $\mathcal {P}_i=p_{1r}/p_{0r}$ and $\mathcal {P}_t=p_{2t}/p_{0t}$ correspond to the pressure jump functions across the reflected, incident and transmitted shocks, respectively. On the other hand, the continuity in the streamwise velocity $U_2(x<0)=U_2(x>0)$ renders

(2.5)\begin{equation} \mathcal{M}_{i}(1-\mathcal{R}_i^{{-}1})-\mathcal{M}_{r}(1-\mathcal{R}_r^{{-}1})\sqrt{\dfrac{\kappa_{1r}}{\kappa_r}\dfrac{\mathcal{P}_i}{\mathcal{R}_i}} =\mathcal{M}_{t}(1-\mathcal{R}_t^{{-}1})\sqrt{\dfrac{\kappa_t}{\kappa_r}\dfrac{1-\mathcal{A}_0}{1+\mathcal{A}_0}} \end{equation}

that, as it can be seen, involves more parameters. The functions $\mathcal {M}_r=(D_{R}+U_1)/c_{1r}$, $\mathcal {M}_i=D_{I}/c_{0r}$ and $\mathcal {M}_t=D_{T}/c_{0t}$ stand for the reflected, incident and transmitted shock Mach numbers, respectively. The $c$ functions identify the corresponding sound speeds, that must be computed upon determination of the associated EoS. Equation (2.5) also depends on the corresponding density jump functions, which read, $\mathcal {R}_r=\rho _{2r}/\rho _{1r}$, $\mathcal {R}_i=\rho _{1r}/\rho _{0r}$ and $\mathcal {R}_t=\rho _{2t}/\rho _{0t}$. Additionally, the associated sonic coefficients, that in the general case are state functions, are defined as $\kappa = \rho c^2 / p$, where, again, $c$ represents the speed of sound. Finally, the post shock Mach numbers are defined as $\mathcal {M}_{2r}=(D_{R}+U_2)/c_{2r}$, $\mathcal {M}_{2i}=(D_{I}-U_1)/c_{1r}$ and $\mathcal {M}_{2t}=(D_{T}-U_2)/c_{2t}$.

Therefore, (2.4) and (2.5), along with the corresponding RH relations (2.1a–c), the EoS and the preshock parameters, can be solved for a given incident shock intensity. This process has been conducted in the supplemental material for the vdW EoS (including the specific case of the ideal gas) and for simple metals modelled with the three-term EoS (applied to aluminium and copper). Our analysis presumes the existence of a unique solution to this problem, although this is not always guaranteed. In most cases involving ‘regular’ EoS, such as those discussed in the supplementary material, a unique solution either exists for the transmitted and reflected shocks. The case that renders the reflected rarefaction wave being reflected is left for a separate analysis. However, since we are considering an arbitrary EoS, we must acknowledge the possibility of anomalous situations. For instance, an exponentially unstable shock scenario might arise, where, instead of a single transmitted shock front, two fronts are produced. Such situations are discussed in Menikoff & Plohr (1989) and Kuznetsov (1989), and they are beyond the scope of this work.

2.2. Linearised perturbation equations

We consider a perturbed interface, separating two stagnant materials, that is weakly curved along the transverse coordinate $y$, before the incident shock first reaches it ($t\leqslant 0^{-}$), in the form of $\psi _i(t\leqslant 0^{-},y)=\psi _{i}^{0^-}\cos (ky)$, where the amplitude of the perturbation is $\psi _{i}^{0^-}$ and $k$ stands for the corresponding wavenumber. For linear theory to be applicable, the condition $\psi _{i}^0k\sim \varepsilon \ll 1$ must hold, with $\varepsilon$ being defined as the small parameter characterising the first-order perturbations magnitude. In what follows, all perturbation variables are of the same order $\varepsilon$. The dimensionless spatial coordinates are defined by $\bar {x}=kx$ and $\bar {y}=ky$, with the system of reference being fixed at the unperturbed interface, which moves with velocity $U_2 \,\hat {x}_0$ in the reference frame attached to the upstream stagnant gas.

When the incident shock wave has totally crossed the rippled interface, $t\geqslant 0^{+}$, two corrugated shock waves, the transmitted and the reflected one, are generated. The interfacial ripple amplitude, conveniently normalised with respect to its preshock value,

(2.6)\begin{equation} \frac{\psi_{i}(t,\bar{y})}{\psi_{i}^{0^-}}=\xi_{i}(t)\cos(\kern0.7pt\bar{y}), \end{equation}

is initially compressed by the shock passage, thereby modifying its initial amplitude to $\xi _{i}^{0}=\xi _{i}(t=0^+)=1-U_2/D_{I}$, provided in (S.2.11) of the supplementary material. Note that the postshock amplitude differs from the preshock amplitude due to compressible effects of the shock. Since $\rho U ={\rm const}$, across the shocks, changes in density lead to corresponding adjustments in velocity. The corresponding shock front amplitudes are also scaled with $\psi _{i}^{0^-}$, thus, the transmitted and reflected shock ripples are equally reduced:

(2.7a,b)\begin{equation} \frac{\psi_{t}(t,y)}{\psi_{i}^{0^-}}=\xi_{t}(t)\cos(\kern0.7pt\bar{y}) \quad\text{and}\quad \frac{\psi_{r}(t,y)}{\psi_{i}^{0^-}}=\xi_{r}(t)\cos(\kern0.7pt\bar{y}). \end{equation}

It is straightforward to see that, after incident shock refraction, the initial ripple amplitude of the transmitted and the reflected shocks can be determined by the relative distances between the most advanced parts of each front at the time the shock fully crosses the interface, thereby setting in their initial conditions through $\xi _{t}(t=0^+)=\xi _{t}^0=1-D_{T}/D_{I}$ and $\xi _{r}(t=0^+)=\xi _{r}^0=1+D_{R}/D_{I}$, where the relations $D_{T}/D_{I}$ and $D_{R}/D_{I}$ depend on the initial conditions and the corresponding EoS, as explicitly shown in (S.2.10) of the supplementary material.

All flow variables are disturbed behind the reflected and transmitted shocks, with the corresponding amplitude being of the same order as the perturbation parameter $\varepsilon$. This is used to define the dimensionless order-of-unity functions (bar symbols) for the pressure, density, longitudinal velocity and transverse velocity variables, namely

(2.8a)$$\begin{align} \frac{\delta p_{j}}{\rho_{j}c_{j}^2}&=\varepsilon{\bar{p}}_{j}({\bar{x},\tau_j})\cos(\kern0.7pt{\bar{y}}), \quad \frac{\delta u_{j}}{c_{j}}=\varepsilon{\bar{u}}_{j}({\bar{x},\tau_j})\cos(\kern0.7pt{\bar{y}}), \end{align}$$

(2.8b)$$\begin{align}\frac{\delta \rho_{j}}{\rho_{j}}&=\varepsilon{\bar{\rho}}_{j}({\bar{x},\tau_j})\cos(\kern0.7pt{\bar{y}}), \quad\frac{\delta v_{j}}{c_{j}}=\varepsilon{\bar{v}}_{j}({\bar{x},\tau_j})\sin(\kern0.7pt{\bar{y}}), \end{align}$$

where subscript ‘$j$’ indicates whether we are identifying the perturbations behind the transmitted shock (denoted as ‘$j=t$’) or the reflected shock (denoted as ‘$j=r$’). Note that definitions in (2.8) utilise the periodic symmetry of the flow within the linear regime. Note that the spatial coordinates and velocities are measured in the reference frame of the interface, which moves together with the compressed material for $\tau _j>0$, where $\tau _{j}$ is dimensionless time constructed with the aid of the local speed of sound: $\tau _{j}=k c_{2j} t$. This dimensionless transformation proves to be advantageous when expressing the linearised Euler equations since they become devoid of parameters. As usual, the linear Euler equations have been obtained by retaining the first-order perturbations $O(\varepsilon )$ and neglecting those of higher order on either side of the interface, see

(2.9a)$$\begin{gather} \frac{\partial {\bar{\rho}}_{j}}{\partial \tau_{j}}+\frac{\partial {\bar{u}}_{j}}{\partial {\bar{x}}}+{\bar{v}}_j=0, \end{gather}$$

(2.9b)$$\begin{gather}\frac{\partial {\bar{u}}_{j}}{\partial \tau_{j}}+\frac{\partial {\bar{p}}_{j}}{\partial {\bar{x}}}=0, \end{gather}$$

(2.9c)$$\begin{gather}{\frac{\partial {\bar{v}}_{j}}{\partial \tau_{j}}-{\bar{p}}_{j}}=0, \end{gather}$$

(2.9d)$$\begin{gather}\frac{\partial {\bar{p}}_{j}}{\partial \tau_{j}}-\frac{\partial {\bar \rho}_{j}}{\partial \tau_{j}}=0, \end{gather}$$

for the mass, streamwise momentum, transverse momentum and energy conservation equations, respectively, with the latter being reduced to the isentropic condition. The Euler equations can be rearranged to give the sound wave equation which, in this particular case, takes the form of the Klein–Gordon equation,

(2.10)\begin{equation} \frac{\partial^2{\bar{p}}_{j}}{\partial \tau^{2}_{j}}=\frac{\partial^2{\bar{p}}_{j}}{\partial {\bar{x}}^{2}_{j}}-{\bar{p}}_{j}. \end{equation}

The noticeable advantage of working with two free-parameter sound wave equations comes with the price of handling two temporal scales, which are interconnected by the factor

(2.11)\begin{equation} \beta=\frac{\tau_{r}}{\tau_{t}}=\frac{c_{2r}}{c_{2t}}. \end{equation}

This parameter plays a crucial role in the acoustic coupling between the reflected and transmitted shocks, as will be discussed in detail later.

The sound wave equations at both sides of the interface must be integrated, upon providing the corresponding initial and boundary conditions. The initial conditions are given by the initial shock ripple amplitudes, $\xi _t^0$ and $\xi _r^0$, provided in (2.7a,b) along with the isobaric condition ${\bar {p}}_{j}(\tau _j=0^+)=0$. The boundary conditions, on the other hand, are dictated by the perturbed RH equations and the equilibrium condition at the interface. For the reflected shock, at $\bar {x}_{rs}=-\mathcal {M}_{2r}\tau _r$, and the transmitted shock, at $\bar {x}_{ts}=\mathcal {M}_{2t}\tau _t$, respectively, they read as

(2.12a)$$\begin{gather} \dfrac{{\rm d} \xi_r}{{\rm d} \tau_r}={-}\frac{(1+h_{r})\mathcal{R}_{r}}{2\mathcal{M}_{2r}(\mathcal{R}_{r}-1)}{\bar{p}}_{rs}\quad\text{and} \quad\dfrac{{\rm d} \xi_t}{{\rm d} \tau_t}=\frac{(1+h_{t})\mathcal{R}_{t}}{2\mathcal{M}_{2t} (\mathcal{R}_{t}-1)}{\bar{p}}_{ts}, \end{gather}$$

(2.12b)$$\begin{gather}{\bar{u}}_{rs}={-}\frac{1-h_{r}}{2\mathcal{M}_{2r}}{\bar{p}}_{rs}\quad\text{and}\quad{\bar{u}}_{ts}= \frac{1-h_{t}}{2\mathcal{M}_{2t}}{\bar{p}}_{ts}, \end{gather}$$

(2.12c)$$\begin{gather}{\bar{\rho}}_{rs}={-}\frac{h_{r}}{\mathcal{M}_{2r}}{\bar{p}}_{rs}\quad\text{and}\quad{\bar{\rho}}_{ts}={-}\frac{h_{t}}{\mathcal{M}_{2t}}{\bar{p}}_{ts}, \end{gather}$$

(2.12d)$$\begin{gather}{\bar{v}}_{rs}=\mathcal{M}_{2r}(\mathcal{R}_{r}-1)\xi_{r}\quad\text{and}\quad{\bar{v}}_{ts}={-} \mathcal{M}_{2t}(\mathcal{R}_{t}-1)\xi_{t}, \end{gather}$$

where the subscript ‘$s$’ is used here to indicate that perturbations are evaluated at their respective reflected and transmitted shock position. The functions $\mathcal {M}_{2j}$, $\mathcal {R}_{j}$ and $h_{j}$, refer to the associated postshock Mach number, density jump and the D'yakov–Kontorovich (DK) parameter, named after the pioneering works of D'yakov (1954) and Kontorovich (1957). For the sake of generality, they are included here as independent parameters, yet they are ultimately determined by the initial conditions, the EoS of the fluids and incident shock intensity, as demonstrated in the supplementary material for vdW EoS and three-term EoS for condensed materials.

The RH conditions (2.12) can be further reduced by taking the total derivative of the streamwise velocities $\bar {u}_j$ along the corresponding shock trajectories, together with the linear Euler equations (2.9). They read as

(2.13a)$$\begin{gather} \left.\left(\frac{1-h_{t}}{2\mathcal{M}_{2t}}+\mathcal{M}_{2t}\right) \frac{\partial {\bar{p}}_{t}}{\partial\tau_{t}}\right|_{s}+\left.\left( \frac{1-h_{t}}{2}+1\right)\frac{\partial {\bar{p}}_{t}}{\partial{\bar{x}}}\right|_{s}={-}\mathcal{M}^{2}_{2t}(\mathcal{R}_{t}-1)\xi_{t}, \end{gather}$$

(2.13b)$$\begin{gather}\left.\left(\frac{1-h_{r}}{2\mathcal{M}_{2r}}+\mathcal{M}_{2r}\right) \frac{\partial {\bar{p}}_{r}}{\partial\tau_{r}}\right|_{s}-\left.\left( \frac{1-h_{r}}{2}+1\right)\frac{\partial {\bar{p}}_{r}}{\partial{\bar{x}}}\right|_{s}=\mathcal{M}^{2}_{2r}(\mathcal{R}_{r}-1)\xi_{r}, \end{gather}$$

for the transmitted and reflected shocks, respectively. Equations (2.13a) and (2.13b), complemented with (2.12a) for the shock ripples evolution, form a closed system of equations for the perturbed shocks.

The problem formulation is completed after providing the mechanical equilibrium condition at the contact interface: pressure and velocity are continuous across the interface. In terms of pressure, they read as

(2.14a)$$\begin{gather} {\bar{p}}_{t}(\tau_{t},0)=\vartheta\,{\bar{p}}_{r}(\tau_{r},0), \end{gather}$$

(2.14b)$$\begin{gather}\left.\frac{\partial {\bar{p}}_{t}}{\partial {\bar{x}}}\right|_{\bar{x}=0}=\left.\beta^{2}\frac{\partial {\bar{p}}_{r}}{\partial {\bar{x}}}\right|_{\bar{x}=0}, \end{gather}$$

where $\vartheta =\kappa _{2r}/\kappa _{2t}$ is the ratio of the sonic coefficients at both sides of the interface. For ideal gas EoS, this ratio reduces to $\gamma _r/\gamma _t={\rm const}.$, but this function depends on the corresponding thermodynamic state of the material, and therefore on the shock intensity, for more complex EoS, as provided in the supplementary material.

This problem can be readily numerically integrated using the method of characteristics, as done in Calvo-Rivera et al. (2023) for piston-driven shocks, yielding the spatial and temporal pressure distribution within the region bounded by the reflected and transmitted shocks. Nevertheless, the theoretical developments in the following sections, based on the formulation presented in Wouchuk (2001), will provide analytical expressions for key quantities, including the asymptotic growth rate of the RMI.

3. The pressure field

3.1. The resolution of the pressure field

This section examines the physical properties of the pressure field, drawing on the mathematical properties of the corresponding Laplace transform detailed in Appendix A. Therefore, although the discussion primarily addresses physical concepts, it is necessary to explicitly reference some of the properties of the associated Laplace transform functions. In particular, the evolution of the pressure field is derived by applying the inverse transformation to the Laplace transform functions described in Appendix A.2).

In the transmitted-shock and reflected shock regions, it can be computed with the aid of

(3.1a)\begin{equation} \bar{p}_{j}(\eta_j,\chi_j)=\frac{1}{2{\rm \pi} \mathrm{i}}\int_{r-\mathrm{i}\infty}^{r+{\rm i}\infty}\mathbb{P}_{j}(s_j,\chi_j){\rm e}^{s_j \eta_j}\,{\rm d}s_j, \end{equation}

where $\mathbb {P}_{j}$ is the pressure Laplace transform with respect to the coordinate $\eta _j$, $s_j$ is the Laplace variable, $\mathrm {i}=\sqrt {-1}$ is the imaginary unit and $r$ is a real number to the right of the singularities of $\mathbb {P}_{j}(s_j)$. The independent variables $\eta _j$ and $\chi _j$ refer to the hyperbolic coordinates

(3.2a,b)\begin{equation} {\bar{x}}=\eta_{j}\sinh\chi_{j}\quad\text{and}\quad \tau_{j}=\eta_{j}\cosh\chi_{j} \end{equation}

defined conveniently to transform the triangular integration domain into a rectangular one (Zaidel’ 1960; Briscoe & Kovitz 1968).

Keeping in mind the analysis of the singularities within the complex plane detailed in Appendix A.4, there are two contributions to the Bromwich integral, the first one is the integral around the integration contour in which we introduce the branch cut to avoid the multivalued nature of the complex functions (A10a) and (A10b), in particular from the line integrals along the imaginary axis, specifically between the branch points. The second contribution comes from the existence of imaginary poles inside the contour, which are calculated by the application of the Cauchy residue theorem. To clarify the analysis, we split the contributions in the form

(3.3)\begin{equation} {\bar{p}}_{j}\left(\eta_{j},\chi_{j}\right)=\lim_{R\rightarrow\infty}\frac{1}{2{\rm \pi} \mathrm{i}}\left\{\oint_{\varGamma}{\mathbb{P}}_{j}(s_{j},\chi_{j}){\rm e}^{\eta_{j} s_j}\,{\rm d}s_j-\int_{C_j}{\mathbb{P}}_{j}(s_{j},\chi_{j}){\rm e}^{\eta_{j} s_j}\,{\rm d}s_j\right\}, \end{equation}

where the first term of the right-hand side, obtained by the residue theorem, is associated with the non-decaying oscillations of a shock front producing SAE. The second term, which requires the evaluation of the integral along each path of the contour excluding the vertical line that pass through $r$, as sketched in Appendix A.4, is associated with the damped oscillations.

Beginning with the latter, it should be noted that the value of the ratio of the speeds of sound, $\beta$, determines the branch cut on each side. When $\beta >1$, the influence of the reflected shock oscillations involves a frequency that is higher than that of the transmitted shock oscillations. Therefore, the branch cut must be extended to range $\{-\mathrm {i}\beta,\mathrm {i}\beta \}$ for the transmitted side but not for the reflected one, see Appendix A.4. When $\beta <1$, the opposite happens, the oscillations of the transmitted shock exhibit a higher frequency than the reflected shock, signifying into an extension of the branch cut of the reflected side with the range $\{-\mathrm {i}\beta ^{-1},\mathrm {i}\beta ^{-1}\}$ due to the change in the dominant evanescent acoustic waves. It results in the closed integral expression for the pressure field which read as

(3.4)\begin{equation} {\bar{p}}_{j}^{{dec}}(\eta_{j},\chi_{j})=\frac{2}{\rm \pi}\int_{0}^{\zeta_{j}}{\rm Im}\{\mathbb{P}_j\left(-\mathrm{i} z ,\chi_j\right)\}\sin(z \eta_j)\,{\rm d}z,\end{equation}

where $\mathbb {P}_{j}$ is given by the solution of the functional system corresponding to each zone. The domain of validity of (3.4) goes $0\leqslant \chi _t \leqslant \tanh ^{-1}\mathcal {M}_{2t}$, for the transmitted side, and $-\tanh ^{-1}\mathcal {M}_{2r} \leqslant \chi _r \leqslant 0$ for the shock-reflected region within the domain $0\leqslant \chi _t \leqslant \tanh ^{-1}\mathcal {M}_{2t}$. Also $\zeta _j$ depends on the value of $\beta$, being

(3.5a,b)\begin{equation} \zeta_t=\begin{cases} 1 & \text{if } \beta \leqslant 1 \\ \beta & \text{if } \beta > 1 \end{cases}, \quad \zeta_r= \begin{cases} \beta^{{-}1} & \text{if } \beta < 1 \\ 1 & \text{if } \beta \geqslant 1 \end{cases}, \end{equation}

for the transmitted and reflected cases, respectively. The superscript ${dec}$ is used in (3.4) because the branch cut contribution always decays, similar to Bessel functions. These terms represent the typical transient growth contribution for the RMI and are consistently reported in the literature.

To obtain the contribution from the residues, it is necessary to consider the imaginary poles that arise in the complex plane, which are the isolated singularities of (A21a) and (A21b) for the transmitted and reflected cases, respectively. Such singularities are the complex roots of $D(s_{j})=0$ in (A22). The amplitude of the non-decaying oscillations is readily given by the residue theorem

(3.6)\begin{equation} \bar{p}_{js}^{osc}\left(\tau_{j}\right)=\frac{1}{2{\rm \pi} \mathrm{i}}\oint_{\varGamma}{\mathbb{P}}_{js}(s_{j},\chi_{j}){\rm e}^{\tau_{j} s_j}\,{\rm d}s_{j}=\sum_{k=1}^{n}{\rm Res}\{{\rm e}^{\tau_{j} s_j}{\mathbb{P}}_{js}\left(s_{j}\right),s_{I,jk}\}, \end{equation}

where $s_{I,jk}$ is the $k$th imaginary pole in (A24).

We consider the case of the shock natural frequency, i.e. there are no secondary modes associated with reverberations with the interface and/or the reflected shock. In this case,

(3.7) \begin{align} &\lim_{s\rightarrow {\rm i}s_{I,j}}(s-{\rm i}s_{I,j}){\mathbb{P}}_{js}(s_{j})\exp({s_{j}\tau_{j}\cosh{\chi_{js}}})+ \lim_{s\rightarrow -{\rm i}s_{I,j}}(s+{\rm i}s_{{I},j}){\mathbb{P}}_{js}(s_{j})\exp({-s_{j}\tau_{j}\cosh{\chi_{js}}})\nonumber\\ &\quad = {\rm \pi}_{j}\sin{(\omega_j\tau_{j})}, \end{align}

with use made of the symmetry property

(3.8)\begin{equation} {\rm Im}\left\{{\mathbb{P}}_{js}\left[q_{j}\left({\rm i}s_{I,j}\right)\right]\right\}={-}{\rm Im}\left\{{\mathbb{P}}_{js}\left[q_{j}\left(-{\rm i}s_{I,j}\right)\right]\right\}. \end{equation}

Using the same decomposition of contributions as in (A21), the shock pressure amplitude associated with non-decaying oscillations can be expressed as

(3.9)\begin{equation} {\rm \pi}_{j}=\dfrac{\sigma_{cj}-\sigma_{bj}s_{I,j}^{2}+s_{I}\sqrt{s_{I,j}^{2}-1}}{\left(s_{I,j}^{2}+s_{R,j}^{2}\right)\left({\sigma_{bj}^{2}-1}\right)}\left[ \frac{\sigma_{dj}}{s_{I,j}}+ {\rm Im}\left\{F_{j}^{{\pm}}\left[q\left({\rm i}s_{I,j}\right)\pm\chi_{js}\right]\right\}\right]. \end{equation}

3.1.1. Pressure perturbations at the shock fronts

The evolution of the pressure field at the shocks is readily determined by the superposition of the decaying and non-decaying oscillating contributions, namely

(3.10)\begin{equation} {{\bar{p}}_{js}(\tau_{j})= \bar{p}}_{j}^{{dec}}(\eta_{j}=\sqrt{1-\mathcal{M}_{2j}^2}\,\tau_j,\chi_{j}=\chi_{js})+{\rm \pi}_{j}\sin(\omega_{js}\tau_{j}) \end{equation}

when that SAE condition is met. Otherwise, only the decaying term contributes.

To illustrate the different possibility that may arise depending on the shock properties and the EoS parameters, figure 2 illustrates the pressure perturbation history at the shocks for four distinct cases. The transmitted shock coordinate corresponds to the positions $\bar {x}_{ts}=\mathcal {M}_{2t} \tau _t >0$ (right half-space) and the reflected shock to $\bar {x}_{rs}=-\mathcal {M}_{2r} \tau _r <0$ (left half-space). In figure 2(a), the pressure field is computed for the same case as in Yang et al. (1994), representing a shock wave propagating from air to ${\rm SF}_{6}$ with an incident Mach number $\mathcal {M}_{i}=1.24$. Figure 2(b) depicts a scenario with similar gas conditions but a stronger incident shock, $\mathcal {M}_{i}=3$. Here, the pressure field of the transmitted shock changes its initial slope sign and exhibits a sudden rise that subsides quickly. Figure 2(c) shows computations where SAE occurs at the transmitted shock side. The asymptotic non-decaying contribution is indicated with a dashed line. This case involves a highly compressible vdW gas with significant covolume and Coulomb contributions on the heavy gas side, similar to those described in Bates & Montgomery (1999), Huete et al. (2021) and Calvo-Rivera et al. (2023). The conditions are those marked with a triangle in figure S2 of the supplementary material. Figure 2(d) presents a computation where the parameter $\beta$ is less than unity, necessitating a different contour path in its resolution via the ILT. For this case, it is noticeable the change in the amplitude of the oscillations of the reflected shock. Finally, figure 2(e,f) illustrate the results for a shock wave propagating through aluminium and impacting copper, leading to the reflection of a backward-moving shock wave. The incident Mach numbers are $\mathcal {M}_{i}=1.3$ and $\mathcal {M}_{i}=2.5$, corresponding approximately to the weak and strong test cases analysed in Tahir et al. (2011), with incident shock pressures of approximately $p_1= 0.28$ Mbar and $2.9$ Mbar, respectively. Base-flow conditions are described in figure S3 of the supplementary material.

Figure 2. Pressure perturbation at the respective shocks. The transmitted shock corresponds to the positions $\bar {x}_{ts}=\mathcal {M}_{2t} \tau _t >0$ and reflected shock to $\bar {x}_{rs}=-\mathcal {M}_{2r} \tau _r <0$. Computations are done for different shock strengths and for ideal-gas (a,b), for vdW gas (c,d) and simple metals (e,f).

The six panels in figure 2 demonstrate the robustness of the ILT method compared with the previously used separation of variables technique, which describes the pressure field in terms of Fourier–Bessel series. Unlike the separation of variables method, the use of which is limited by the requirement of convergence of its coefficients, the ILT method consistently provides accurate solutions even under varying conditions and long-time regimes, this is due to the boundary conditions being time dependent which makes it necessary to obtain a number of coefficients that tends to infinity for larger times, as seen in the first chapter of Duffy (2004).

3.1.2. Influence of the pressure disturbances on the interface

Once we understand the evolution of the pressure field at the shock fronts and the potential for SAE, it is natural to consider the impact on the interface. Ultimately, this affects the development of the RMI, which is quantified in § 5. Starting from the decaying contribution, a direct variable change transformation in (3.4) yields the pressure field in Cartesian coordinates,

(3.11)\begin{equation} {\bar{p}}_{j}^{{dec}}(\bar{x}_{j},\tau_{j})={\bar{p}}_{j}^{{dec}}\left[\eta_{j}= \sqrt{\tau_j^2-\bar{x}^2},\chi_{j}=\tanh^{{-}1}\left({\bar{x}}/{\tau_j}\right)\right],\end{equation}

which is applicable within the domain enclosed by the corresponding shock and the interface.

To better understand the non-decaying oscillating contribution, we can isolate this effect by analysing a scenario where a shock oscillates with constant frequency and amplitude. First, we refer again to the dispersion relation (A22). From the analysis of the singularities in the complex plane performed before, purely imaginary poles arise when $\sigma _{cj}>\sigma _{bj}$, which ultimately translates into oscillations with frequency $\smash {\omega _{js}>(1-\mathcal {M}_{2j}^{2})^{1/2}}$. Using a normal mode decomposition for the acoustic frequency and wavenumber, referred to from now on as the acoustic eigenvalues, $\omega _{j}$ and ${k}_{j}$, the shock oscillation frequency is related to the acoustic eigenvalue through the sound wave equation $\smash {\omega _{j}^{2}={k}_{j}^{2}+1}$ and the compatibility condition $\smash {\omega _{js}=\omega _{j}-\mathcal {M}_{2j}{k}_{j}}$, to give

(3.12a,b)\begin{equation} {k}_{j}=\dfrac{\omega_{js} \mathcal{M}_{2j}-\sqrt{\omega_{js}^{2}-(1-\mathcal{M}_{2j}^{2})}}{1-\mathcal{M}_{2j}^{2}},\quad \omega_{j}=\dfrac{\omega_{js}- \mathcal{M}_{2j}\sqrt{\omega_{js}^{2}-(1-\mathcal{M}_{2j}^{2})}}{1-\mathcal{M}_{2j}^{2}}. \end{equation}

As the wrinkled shock oscillates, it generates pressure disturbances in the downstream fluid in the form of sound waves (see § 90 of Chapter IX in Landau & Lifshitz (1987)). A simple inspection of ${k}_{j}$ in (3.12a,b) reveals that, depending on the shock oscillation frequency, these sound waves may or may not directly impact the contact interface. Therefore, three distinct regimes can be identified, as follows.

1. (i) When $\smash {\omega _{ts}^2 < 1-\mathcal {M}_{2j}^{2}}$, pressure disturbances at the shock decay according to the power-law $\smash {t^{-3/2}}$ and diminish exponentially with distance from the shock, resulting in evanescent waves in the longitudinal direction. This scenario corresponds to the typical case studied in RMI for both transmitted and reflected shock fronts moving in ideal gases. Therefore, the acoustic influence only affects the interface at early times, as the exponential decay quickly reduces its relative contribution.

2. (ii) When $\smash {1-\mathcal {M}_{2j}^{2}\leqslant \omega _{ts}^2<1}$, the stability behaviour of the shock front switches to a neutrally stable regime, and non-decaying acoustic waves are radiated by the shock. This regime, referred as SAE, is a manifestation of the DK instability (Calvo-Rivera et al. 2023). In such case, the acoustic eigenvalues become real, with the wavevector pointing to the radiating shock direction. In the theory of dynamical systems, this corresponds to a transition through a Hopf-like bifurcation, see Clavin & Williams (2009, 2012), to a limit cycle which will exhibit permanent oscillations, as explained in pp. 252–257 of Strogatz (2015). However, regarding the interface behaviour, the contact interface will not be affected by the acoustic waves, as they cannot reach the hydrodynamic interface. The permanently oscillating waves travel to the right in the compressed gas (and interface) reference frame but not reaching the shock front or the interface, as sketched in figure 3 for the transmitted shock case.

3. (iii) When $\omega _{js} \geqslant 1$, the behaviour of the corresponding shock front dynamics remains the same, but the acoustic influence changes as the longitudinal projection of the wavenumber in the compressed material reference frame points backwards. Consequently, sound waves propagate to the left, reaching the interface, as sketched in figure 4 for the transmitted shock.

Further scenarios can also be explored as $\omega _{js}$ increases beyond unity: the radiated sound waves can reach the reflected shock, causing it to oscillate with constant amplitude and frequency, and/or they can be reflected back towards the transmitted shock, imposing a second mode of oscillation, similar to what occurs with piston-driven shocks (Calvo-Rivera et al. 2023). However, these cases are rare unless the shock oscillation frequency is externally imposed, as described in Velikovich et al. (2007), where the shock travels across a non-uniform medium.

Figure 3. Sketch of the acoustic influence of permanent oscillations for case (ii). Solid grey lines represent the acoustic characteristics in the contact surface reference frame, corresponding to the grey dashed lines. The trajectory of the constant-amplitude acoustic field is indicated by the orange dashed lines.

Figure 4. Sketch of the acoustic influence of permanent oscillations for case (iii). Solid grey lines represent the acoustic characteristics in the contact surface reference frame, corresponding to the grey dashed lines. The trajectory of the constant-amplitude acoustic field is indicated by the orange dashed lines (radiated by the transmitted shock), while the blue dashed lines and green lines denote the acoustic paths transmitted and reflected across the interface, respectively.

To facilitate the discussion, we consider the specific case where SAE occurs only in the transmitted shock, i.e. $\sigma _{bt}<\sigma _{ct}$, as depicted in figures 3 and 4 for cases (ii) and (iii), respectively. As per the former, the pressure field associated with the sound radiation by the transmitted shock, which also corresponds to the pressure field in the long-time regime in this case, can be expressed as a piecewise function of $\bar {x}/\tau _{t}$, namely

(3.13)\begin{equation} \bar{p}_{t}^{osc}(\tau_{t},{\bar{x}}) = \begin{cases} {\rm \pi}_t\sin\left(\omega_{t}\tau_{t}-{k}_{t} \bar{x}\right) & \text{for}\ \dfrac{k_{t}}{\omega_{t}}<\dfrac{\bar{x}}{\tau_{t}}<\mathcal{M}_{2t}\\[10pt] 0 & \text{for}\ -\beta\mathcal{M}_{2r}<\dfrac{\bar{x}}{\tau_{t}}<\dfrac{k_{t}}{\omega_{t}} \end{cases} \end{equation}

where the non-zero pressure field corresponds to the orange-shadowed region in figure 3. Notice the decaying contribution given in (3.11) is present within the whole domain. The two limiting cases that bound case (ii) are $\omega _{ts}\rightarrow (1-\mathcal {M}_{2j}^{2})^{1/2}$ and $\omega _{ts}\rightarrow 1$. Per the former, the value of $k_{t}/\omega _{t}\rightarrow \mathcal {M}_{2t}$, which corresponds to the trajectory of the transmitted shock wave, thereby eliminating the orange-shadowed region. Regarding the latter, the value of $k_{t}/\omega _{t}\rightarrow 0$, which corresponds to the trajectory of the entropy-vorticity waves (vertical lines in figure 3). Such a case renders a standing acoustic wave and the domain of influence for the orange-shadowed region is the space between the transmitted shock and the interface.

For the backwards running waves to reach the contact interface (case (iii)), the acoustic wave number ${k}_t$ must be negative, implying that $\omega _{ts} > 1$, which is a more stringent condition to satisfy if the shock oscillation frequency is not externally imposed. Nonetheless, it is worthwhile to outline the consequences if this condition is met, which can be readily anticipated by simple inspection of figure 4. It can be seen that four distinguished regions are identified, since the train of acoustic waves radiated by the transmitted shock are partially reflected from and transmitted through the contact interface, as detailed below. The following derivation is essentially the same as the one found in § 66 of Chapter 8 in Landau & Lifshitz (1987).

As in (3.13), the hyperbolic character of the problem allows as to write the pressure field as a piecewise function of $\bar {x}/\tau _{t}$, namely

(3.14) \begin{align} \bar{p}_{t}^{osc}(\tau_{t},{\bar{x}}) = \begin{cases} {\rm \pi}_t\sin\left(\omega_{t}\tau_{t}-{k}_{t} \bar{x}\right), & \text{for}\ -\dfrac{k_{t}}{\omega_{t}}<\dfrac{\bar{x}}{\tau_{t}}<\mathcal{M}_{2t}, \\[10pt] {\rm \pi}_t\left[\sin\left(\omega_{t}\tau_{t}-{k}_{t} \bar{x}\right)+\mathfrak{R} \sin\left(\omega_{t}\tau_{t}+{ k}_{t} \bar{x}\right) \right], & \text{for}\ 0<\dfrac{\bar{x}}{\tau_{t}}<{-}\dfrac{k_{t}}{\omega_{t}}, \\[10pt] {\rm \pi}_t\mathfrak{T}\sin\left(\omega_{t}\tau_{t}-{k}_{r} \bar{x}\right), & \text{for}\ \dfrac{{k}_{r}}{\omega_{t}}<\dfrac{\bar{x}}{\tau_{t}}<0, \\[10pt] 0, & \text{for}\ -\beta\mathcal{M}_{2r}<\dfrac{\bar{x}}{\tau_{t}}<\dfrac{{k}_{r}}{\omega_{t}}, \end{cases} \end{align}

where the first region corresponds to the orange-shadowed zone in figure 4 and that is similar to that in figure 3, with the coefficient ${\rm \pi} _t$ being obtained from the residue theorem described before. The following green-shadowed region is composed of the coexisting shock-radiated waves moving backwards and the interface-reflected waves moving forwards, with the amplitude of the latter being readily given by the reflection coefficient for the dimensionless pressure amplitude, $\mathfrak {R}$. And the blue-shadowed region corresponds to the pressure field that crosses the interface with an amplitude that is determined by the transmission coefficient, $\mathfrak {T}$. Such are obtained by applying the boundary conditions at the material surface, (2.14), taking the form of

(3.15a,b)\begin{equation} \mathfrak{R}=\dfrac{{k}_{t}\vartheta-{k}_{r}\beta^{2}}{{k}_{t}\vartheta+{k}_{r}\beta^{2}}\quad \text{and} \quad \mathfrak{T}=\dfrac{2 {k}_{t}}{{k}_{t}\vartheta+{k}_{r}\beta^{2}}, \end{equation}

that involve the wavenumber at the transmitted and reflected sides, respectively, which are attached to the compatibility condition, i.e. both of them must fulfil the sound wave equation being subjected to the same frequency, $\omega _{t}$ giving ${k}_{r}^2=\omega _{t}^{2}/\beta ^{2}-1$ and ${k}_{t}^{2}=\omega _{t}^{2}-1$, and it can be checked that they have the same form of those given in problem 2 of § 84, Chapter IX of Landau & Lifshitz (1987). Notice that, as ${k}_{t}={k}_{t}(\omega _{t})$ and ${k}_{r}={k}_{r}(\omega _{t})$, the acoustic impedance of the contact interface depends on the acoustic frequency because it changes the acoustic wavenumber, as expected. For instance, in the limiting case where $\omega _{ts}\rightarrow \infty$, the value of ${k}_{t}/\omega _{t}\rightarrow -1$ corresponds to the wavenumber of 1-D sound waves (represented by the grey lines). These scenarios extend beyond the current analysis because, with the modified wavenumber, the constant-amplitude sonic waves that are reflected and transmitted will reach both the transmitted and reflected shocks, respectively. This interaction induces at least secondary modes in the transmitted shock, as described by Calvo-Rivera et al. (2023) for piston-driven shocks. The present analysis remains consistent because, as $\beta \rightarrow 1$ and $\vartheta \rightarrow 1$, the coefficients $\mathfrak {R}\rightarrow 0$ and $\mathfrak {T}\rightarrow 1$, corresponding to the case of a weak discontinuity at the interface, addressing the problems already solved in Wouchuk, Huete Ruiz de Lira & Velikovich (2009) and Huete Ruiz de Lira, Velikovich & Wouchuk (2011).

The opposite limit, involving the total reflection of a sonic wave from a perfectly solid wall, cannot be directly addressed by simply taking limits of the parameters introduced here, i.e. $\mathfrak {T}\neq 0$. This is because the nature of the refracted sonic wave fundamentally changes. As explained by Haberman (2013), see Chapter 4, § 6, the refraction–reflection problem is characterised by the angle of incidence of the sonic wave defined by the projection of the wavevector in the transverse direction. In our formulation, this angle is just $\sin ^{-1}(1/\omega _{t})$. When this angle is greater than the critical value, i.e. when $\sin ^{-1}(1/\omega _{t})>\sin ^{-1}(1/\beta )$, the refracted (transmitted acoustic wave in the reflected-shock region) sonic wave is no longer a running wave but it becomes evanescent. Instead, an evanescent wave is transmitted, which decays exponentially with the distance from the interface, as readily seen by the compatibility condition ${k}_{r}^2=\omega _{t}^{2}/\beta ^{2}-1$, since ${k}_{r}$ turns complex when $\omega _t<\beta$. By introducing this consideration and applying the same boundary conditions for the pressure and its gradient, the wavenumber for the refracted wave becomes complex. Consequently, one finds that $\mathfrak {R}=-1$, similar to the piston-driven shocks discussed in Calvo-Rivera et al. (2023).

For the sake of completeness, we briefly comment on the possibility of having permanent acoustic waves at the interface from the Laplace transform equations. With this purpose, it has proved convenient to write the pressure field in the transmitted shock side as functions of the coordinate $\chi _t$ and variable $q_t$,

(3.16)\begin{align} {\mathbb{P}}_{t}(q_{t},\chi_{t})&=\dfrac{F_{t}^{+}\left(q_{t}+\chi_{t}\right)-F_{t}^{+}\left(q_{t}+2\chi_{t}-\chi_{ts}\right)}{\cosh{q_{t}}}+\dfrac{\cosh{(q_{t}+\chi_{ts}-\chi_{t})}}{\cosh{q_{t}}}\nonumber\\ &\quad\times\left[F_{t}^{+}\left(q_{t}+2\chi_{t}-\chi_{ts}\right)\alpha_{t}\left(q_{t}+\chi_{ts}-\chi_{t}\right)-{\mathbb{P}}^{iso}_{ts}\left(q_{t}+\chi_{ts}-\chi_{t}\right)\right], \end{align}

where the equation particularised at the shock (A21a) is recovered if $\chi _{t}=\chi _{ts}$. This equation is particularly useful because it can easily adapt to different Riemann-type configurations. For instance, if we assume that no disturbance can propagate back to the transmitted shock from behind, making, $F_{t}^{+}={\rm const}.$, the expression simplifies to match (48) of Wouchuk et al. (2009) for isolated shocks. Despite the evident differences, a similar analysis can be applied. The first term in (3.16) introduces poles corresponding to the secondary DK modes. Conversely, the second term in (3.16) contains a common denominator, whose zeros are the same as those previously derived from the dispersion relation, but shifted relative to the shock reference frame,

(3.17)\begin{equation} D_{t}\left(\tilde{ s}_{t}\right)=\tilde{s}_{t}\sqrt{\tilde{s}_{t}^{2}+1}+\sigma_{bt}\tilde{s}_{t}^{2}+\sigma_{ct}, \end{equation}

where $\tilde {s}_{t}=\sinh {(q_{t}+\chi _{ts}-\chi _{t})}$. The poles given for each $\chi _{t}$ value, associated with different relative speeds, can be expressed in the form of

(3.18)\begin{equation} s_{I}\left(\chi_{t}\right)=\cosh{\left[\cosh^{{-}1}(s_{I})-\chi_{ts}+\chi_{t}\right]}, \end{equation}

which, through manipulation and substitution of the hyperbolic change of variables, allows us to obtain the pressure field in the spatiotemporal domain,

(3.19)\begin{equation} \bar{p}_{t}={\rm \pi}_t\sin{\left[\cosh{\left(\cosh^{{-}1}(s_{I})-\chi_{ts}\right)}\tau_{t}-\sinh\left(\chi_{ts}-\cosh^{{-}1}(s_{I})\right)\bar{x}\right]}, \end{equation}

where

(3.20a,b)\begin{equation} \omega_{t}=\cosh{\left[\cosh^{{-}1}(s_{I})-\chi_{ts}\right]}\quad\text{and}\quad {k}_{t}=\sinh\left[\chi_{ts}-\cosh^{{-}1}(s_{I})\right], \end{equation}

as a proof of consistency between the Laplace transform and the normal mode analyses. An alternative scenario arises when the reflected shock satisfies $\omega _{rs} > 1$, with the acoustic radiation hitting on the interface and travels to the right from the reflected shock. Since the physical discussion is similar, this case will not be specifically considered in this work.

3.2. The infinitely dense limit for the heavy material: the sequential model

As commented on in Velikovich et al. (2020), when the Atwood number is sufficiently high $1-\mathcal {A}_0\ll 1$, the RMI problem can be reduced, to some degree, to a shock-wall reflection. However, it is convenient to take this limit cautiously because what really distinguishes the interaction is the postshock Atwood number, in our formulation it is $\mathcal {A}_2= (\beta ^2-\vartheta )/ (\beta ^2+\vartheta )$. The condition $1-\mathcal {A}_2\ll 1$ may not be easily anticipated by preshock conditions for complex EoS, so the infinitely dense limit is conveniently rewritten in terms of $\beta ^2/\vartheta \gg 1$. In such a case, the mathematical description of the pressure field can be significantly simplified. In terms of the mathematical formulation, this distinguished limit allows us to obtain analytic expressions for the pressure fields at both sides. When $\beta ^2/\vartheta \gg 1$, expression (A4a–c) for the pressure gradient becomes $\bar {l}_{ri}\simeq 0$, that is equivalent to that with the problem of shock wave reflected by a weakly perturbed wall, as studied in Briscoe & Kovitz (1968) and extended to arbitrary EoS in Calvo-Rivera et al. (2023).

The physics behind the decoupling of the perturbation evolution of reflected and transmitted shocks can be understood as follows if dimensional variables are considered. Due to the significantly higher speed of sound on the reflected-shock side, the corrugated contact surface evolves more slowly compared with the dynamics of the reflected shock, whose oscillations decay before the interface shape undergoes significant changes. As a result, on the time scale of the interface shape variation, the reflected-shock side maintains a constant pressure from the low-density left-hand side of the interface. For the high-density right-hand side, the situation corresponds to a rippled surface of a compressible fluid is instantly loaded with a constant, uniform pressure, resulting in an RMI-like instability, as first analysed by Nikolaev (1965) for the opposite limit $\beta ^2/\vartheta \ll 1$. In our dimensionless formulation, the reflected shock operates independently of any knowledge about the transmitted shock, while the latter is completely influenced by the former. Thus, the resolution occurs in a sequential manner: initially addressing the reflected shock, such as in the piston-reflection problem, and subsequently utilising the pressure field at the contact interface to resolve the transmitted shock. We refer to this as the sequential model. The mathematical details of the function decomposition in the Laplace transform are provided in Appendix A.5. By applying the same procedure used in the previously described ILT, the pressure field in the spatiotemporal domain can be computed.

For the sake of simplicity, the sequential model has been presented and computed in the absence of SAE, since the main purpose is the evaluation of the accuracy in regular conditions. For example, in figure 5(a), the pressure field is computed for the same case as in figure 2(a) corresponding to Yang et al. (1994). Figure 5(b) is a similar case but with higher Atwood number and shock strength. The orange line corresponds to the fully consistent model and the black-dashed line with the sequential model. It is observed that, even when the parameter $\beta ^{2}/\vartheta$ is not very large, the agreement in the shock pressure perturbations is fairly good. However, it should be noted that a relatively good agreement in the pressure evolution of the shocks does not necessarily guarantee a similar match in the evolution of the contact interface, since the simplification stems from a singular perturbation in neglecting the pressure gradient at the interface. Therefore, the validity of this method should be checked by comparing the growth rate given from the full formulation against the present limit, as done in § 5.

Figure 5. Pressure perturbation at the respective transmitted and reflected shocks for different incident Mach and Atwood conditions. Orange-solid line corresponds to the complete model and black-dashed line with the sequential model $\beta ^2/\vartheta \gg 1$.

4. The acoustic, rotational and asymptotic velocity fields

This section is dedicated to the analysis of the velocity field. As previously established in Landau & Lifshitz (1987) and Clavin & Williams (2009, 2012), the velocity field can be decomposed into two distinct types of waves: an acoustic wave, which is generated by the pressure sound waves emitted by each front, and an entropy-vorticity wave, which, under linear stability analysis, represents an isobaric and incompressible flow convected by the unperturbed flow velocity. The general equation for such velocity field can be given by the Helmholtz decomposition in the form of

(4.1)\begin{equation} \dfrac{\partial^{2}{\boldsymbol{\bar{v}}}_{j}}{\partial\tau_{j}^{2}}=\nabla^{2}{\boldsymbol{\bar{v}}}_{j}+ \boldsymbol{\nabla}\times\left(\boldsymbol{\nabla}\times{\boldsymbol{\bar{v}}}_{j}\right), \end{equation}

where $\boldsymbol {\nabla }\times (\boldsymbol {\nabla }\times {\boldsymbol {\bar {v}}}_{j})= \boldsymbol {\nabla }\times {\boldsymbol {\bar {w}}}_{j}$ is the curl of the vorticity field $\boldsymbol {\bar {w}}_j$. It is convenient to use Helmholtz decomposition of the velocity field by the superposition of its acoustic or potential part and its rotational or entropic contribution, namely ${\boldsymbol {\bar {v}}}_j(\boldsymbol {{\bar {x}}},\tau _j)={\boldsymbol {\bar {v}}}_j^{ac}(\boldsymbol {{\bar {x}}},\tau _j)+{\boldsymbol {\bar {v}}}_j^{rot}(\boldsymbol {{\bar {x}}})$. Note that the rotational contribution remains constant for each fluid particle, as derived from the linear Euler equations. Therefore, any changes in the rotational field (or vorticity) must be externally induced, as exemplified here by the corrugated oscillating shocks.

Similar to the pressure field that has already been calculated in the previous section, the velocity field can be expressed as a combination of the decaying and non-decaying components of that pressure field. First, we introduce the acoustic velocity field with contributions from the decaying part, followed by the permanent oscillatory component. This results in a piecewise function for the velocity field, depending on the domain of SAE influence.

4.1. The acoustic wave

The acoustic component of the velocity field can be calculated directly from Euler momentum equations (2.9). With use made of the hyperbolic transformation (3.2a,b), and the corresponding Laplace transforms, the streamwise component of the acoustic velocity field in the transmitted-shock and reflected-shock sides reads as

(4.2a)$$\begin{gather} \mathbb{U}^{ac}_{j}(s_{j},\chi_{j})=\sinh{\chi_{j}}\cosh{\chi_{j}}\mathbb{P}_{j}(s_{j},\chi_{j})- \cosh^{2}{\chi_{j}}\dfrac{\mathbb{L}_{j}(s_{j},\chi_{j})}{s_{j}}, \end{gather}$$

(4.2b)$$\begin{gather}\mathbb{V}^{ac}_{j}(s_{j},\chi_{j})=\cosh{\chi_{j}}\dfrac{\mathbb{L}_{j}(s_{j},\chi_{j})}{s_{j}}+ \dfrac{{\bar{v}}_{ji}^0}{s_{j}}, \end{gather}$$

where ${\bar {v}}_{ji}^0$ is suited for both sides, namely, ${\bar {v}}_{ti}^0=-\mathcal {M}_{2t}(\mathcal {R}_{t}-1)\xi _t^0$ and ${\bar {v}}_{ri}^0=\mathcal {M}_{2r}(\mathcal {R}_{r}-1) \xi _r^0$ are the initial value for the transverse velocity at the interface, which are proportional to the initial shock corrugation amplitudes, defined after (2.7a,b) and explicitly given in (S.2.10) of the supplementary material. The spatiotemporal evolution can be obtained with the ILT, similarly as done in (3.3) for the pressure field. For the streamwise components, we obtain

(4.3)\begin{align} {\bar{u}}_{j}^{ac}(\eta_j,\chi_{j})=\sinh{\chi_j}\cosh{\chi_{j}}\bar{p}^{dec}_{j}(\eta_{j},\chi_{j})- \dfrac{2\cosh^{2}{\chi_{j}}}{\rm \pi}\int_{0}^{\zeta_{j}}{\rm Im} \left\{\dfrac{\mathbb{L}_{j}(\mathrm{-i}z,\chi_{j})}{-\mathrm{i}z}\right\}\sin\left(\eta_{j}z\right){\rm d}z \end{align}

for the transmitted and reflected sides, respectively. On the other hand, for the transverse velocity contribution, we find

(4.4)\begin{equation} {\bar{v}}_{j}^{ac}(\eta_j,\chi_{j})=\dfrac{2\cosh{\chi_{j}}}{\rm \pi}\int_{0}^{\zeta_{j}}{\rm Im}\left\{\dfrac{\mathbb{P}_{j}(-\mathrm{i}z,\chi_{j})}{\mathrm{i}z}\right\}\sin\left(\eta_{j}z\right){\rm d}z+{\bar{v}}_{ji}^0. \end{equation}

In RMI, an important particular case involves the evolution of velocity profiles at the interface. This simplifies the aforementioned equations, as for $\bar {x}=\eta _{j}\sinh {0}=0$ and for the temporal variable $\tau _{j}=\eta _{j}\cosh {0}=\eta _{j}$, yielding

(4.5a)$$\begin{gather} {\bar{u}}_{ti}^{ac}(\eta_t)=\dfrac{2}{\rm \pi}\int_{0}^{\zeta_{t}}{\rm Im}\left\{\dfrac{\mathbb{L}_{t}(\mathrm{i}z,0)} {\mathrm{i}z}\right\}\sin\left(\eta_{t}z\right){\rm d}z , \end{gather}$$

(4.5b)$$\begin{gather}{\bar{v}}_{ti}^{ac}(\eta_t)={\bar{v}}_{t0}-\dfrac{2}{\rm \pi}\int_{0}^{\zeta_{t}} {\rm Im}\left\{\dfrac{\mathbb{P}_{t}(\mathrm{i}z,0)}{\mathrm{i}z}\right\}\sin\left(\eta_{t}z\right){\rm d}z, \end{gather}$$

for the streamwise and transverse velocity contributions, respectively.

As previously mentioned, the expressions presented above are only valid when the shock front oscillates with a frequency lower than $(1-\mathcal {M}_{2j}^{2})^{1/2}$, which was referred as case (i). When the frequency of the shock exceeds this limit, permanent sound pressure waves are emitted in the compressed fluid, as identified previously as cases (ii) and (iii).

For clarity, we will restrict our consideration to $\omega _{st}>1$, associated with cases where the transmitted shock front can enter the SAE regime and these waves reach the interface, i.e. case (iii). In this scenario, the longitudinal velocity field can be expressed as $\bar {u}_{j}^{ac}(\tau _{j},{\bar {x}}) =\bar {u}_{j}^{dec}(\tau _{j}, {\bar {x}}) + \bar {u}_{j}^{osc}(\tau _{j},{\bar {x}})$, where the decaying component is given in (4.3) with the change of variables to Cartesian coordinates. The non-decaying component is straightforward to calculate, as it only involves differentiating and integrating corresponding harmonic functions in (3.14). Consequently, the transmitted- and reflected-side normal velocity components can be expressed as

(4.6a)$$\begin{align} \bar{u}_{t}^{osc}(\tau_{t},{\bar{x}}) = \begin{cases} Q_{t}\sin\left(\omega_{t}\tau_{t}-{k}_{t} \bar{x}\right), & \text{for}\ -\dfrac{k_{t}}{\omega_{t}}<\dfrac{\bar{x}}{\tau_{t}}<\mathcal{M}_{2t}, \\[10pt] Q_{t}\left[\sin\left(\omega_{t}\tau_{t}-{k}_{t} \bar{x}\right)+\mathfrak{R} \sin\left(\omega_{t}\tau_{t}+{k}_{t} \bar{x}\right) \right], & \text{for}\ 0<\dfrac{\bar{x}}{\tau_{t}}<{-}\dfrac{k_{t}}{\omega_{t}}, \end{cases} \end{align}$$

(4.6b)$$\begin{align}\bar{u}_{r}^{osc}(\tau_{r},{\bar{x}}) = \begin{cases} \beta\dfrac{{k}_{r}}{{k}_{t}}Q_t\mathfrak{T}\sin\left(\omega_{t}\tau_{r}\beta^{{-}1}-{k}_{r} \bar{x}\right), & \text{for}\ \dfrac{{k}_{r}}{\omega_{t}}<\dfrac{\bar{x}}{\tau_{t}}<0, \\[10pt] 0, & \text{for}\ -\mathcal{M}_{2r}<\dfrac{\bar{x}}{\tau_{r}}<\dfrac{{k}_{r}}{\omega_{t}}, \end{cases} \end{align}$$

where $Q_{t}={\rm \pi} _{t}{k}_{t}/\omega _{t}$ is the amplitude of the non-decaying contribution. The transverse part is equally simple to obtain and will be omitted here for the sake of conciseness. Upon close examination of these expressions, the influence of the SAE at the interface becomes clear. When the transmitted front oscillates at a frequency higher than one, the wavenumber turns negative. This causes the non-decaying component of the longitudinal interface velocity to oscillate with a mean value around zero, while the decaying part increases from zero to a final constant value. As a result, the oscillatory and decaying components compete over time without damping the oscillations.

4.2. The vortex-entropic wave

Due to the steady state of the rotational field, the governing equation (4.1) can be expressed as $\nabla ^{2}\boldsymbol {{\bar {v}}}_{j}^{rot}=-\boldsymbol {\nabla }\times \boldsymbol {\bar {w}}_{j}$. This vector equation can be decomposed into its streamwise and transverse components,

(4.7a)$$\begin{align} \dfrac{{\rm d}^{2}\bar{u}^{rot}_{j}}{{\rm d}\kern0.7pt {\bar{x}}^{2}}-\bar{u}^{rot}_{j}={-}{\bar{w}}_{j}, \end{align}$$

(4.7b)$$\begin{align}\dfrac{{\rm d}^{2}\bar{v}^{rot}_{j}}{{\rm d}\kern0.7pt {\bar{x}}^{2}}-\bar{v}^{rot}_{j}=\dfrac{{\rm d}{\bar{w}}_{j}}{{\rm d}\kern0.7pt {\bar{x}}}. \end{align}$$

Recognising the incompressibility of the linear rotational velocity field, the latter equation can be replaced with $({\rm d}{\bar {u}^{rot}_j})/({\rm d}\kern0.7pt {\bar {x}})=-{\bar {v}}_{j}^{rot}$, simplifying the differential to be solved to only (4.7a) for each side of the contact interface. This fact allows us to calculate the value of ${\bar {w}}_{j}$ using RH conditions and Euler equations, namely

(4.8)\begin{equation} {\bar{w}}_j ={\pm}\varOmega_j \, \bar{p}_{js}\left(\bar{x}\right)={\pm}\varOmega_j \, \bar{p}_j\left(\bar{x},\tau_j=\dfrac{\bar{x}}{\sinh{\chi_{js}}}\right), \end{equation}

where the positive sign must be chosen for the transmitted shock $j=t$ and the negative sign for the reflected shock $j=r$. Notice that $\tau _t=\bar {x}/\sinh {\chi _{ts}}$ and $\tau _r=\bar {x}/\sinh {\chi _{rs}}$ represent the time when the transmitted and reflected shock fronts, respectively, have just arrived at the evaluated position $\bar {x}$, where $\sinh {\chi _{js}}=\mathcal {M}_{2j}/(1-\mathcal {M}_{2j}^2)^{1/2}$. The amplitude of the vorticity disturbances relative to the shock pressure perturbations is

(4.9)\begin{equation} \varOmega_{j}=\dfrac{\left({\mathcal{R}_{j}}-1\right)\left(1+h_{j}\right)}{2\mathcal{M}_{2j}}. \end{equation}

Equation (4.7a) is a second-order linear differential equation, the solution of which involves a combination of the homogeneous and particular solutions. The particular solution can be determined using the method of undetermined coefficients, employing the integral equation as a proposed solution, as outlined in Wouchuk & Cobos-Campos (2018).

The homogeneous part can be easily solved by using an integrating factor of the form $A {\rm e}^{-\bar {x}}+B {\rm e}^{ \bar {x}}$ and neglecting the non-regular contributions of the general solution, in the case of the transmitted zone $B=0$ while for the reflected one $A=0$. In determining the integration constant, a boundary condition arises from the observation that the only contribution to the vorticity field stems from the oscillations of the shock. Hence, it becomes evident that at $\bar {x}=0$, the rotational contribution to the velocity field is zero. Thus, the ultimate expression for the normal component of the rotational velocity corresponds to the particular solution of (4.7a), which takes the form of

(4.10)\begin{equation} {\bar{u}}^{rot}_{j}(\bar{x})=\dfrac{2\varOmega_{j}\sinh^{2}{\chi_{js}}}{\rm \pi}\int_{0}^{\zeta_{j}}\dfrac{{\rm Im}\{\mathbb{P}_{js}(-{\rm i}z)\}}{z^{2}+\sinh^{2}{\chi_{js}}}\sin\left(z\dfrac{{\bar{x}}}{\sinh{\chi_{js}}}\right){\rm d}z. \end{equation}

Given the incompressibility of the linear rotational field, the transverse contribution can be determined as the derivative with respect to ${\bar {x}}$ of (4.10), that renders

(4.11)\begin{equation} {\bar{v}}_{j}^{rot}(\bar{x})={-}\dfrac{2\varOmega_{j}\sinh{\chi_{js}}}{\rm \pi}\int_{0}^{\zeta_{j}}\dfrac{{\rm Im}\{\mathbb{P}_{js}(-{\rm i}z)\}}{z^{2}+\sinh^{2}{\chi_{js}}}\cos\left(z\dfrac{{\bar{x}}}{\sinh{\chi_{js}}}\right) z\,{\rm d}z. \end{equation}

Notice that, as before the expressions are only valid when the frequency of the shock is lower than $\omega _{ts}<(1-\mathcal {M}^{2}_{2t})^{1/2}$. If one of the shock front exceeds that limit, the shock front enters in the SAE regime. In order to see what modification implicates such a change we consider the case presented in figure 2(c), i.e. the transmitted shock wave exhibits permanent oscillations. In terms of mathematical treatment, this only implies that the transmitted rotational velocity field should contains another particular solution proportional to $\sin {({k}^{rot}_{t} \bar {x})}$, where the rotational wavenumber is given by ${k}^{rot}_{t}=s_{I}(1-\mathcal {M}_{2t}^{2})/\mathcal {M}_{2t}$. The corresponding constant is obtained from the method of undetermined coefficients. Thus, it is possible to include the permanent oscillatory part of the rotational velocity field as

(4.12a)$$\begin{gather} {\bar{u}}^{rot}_{t}(\bar{x})=\dfrac{2\varOmega_{t}\sinh^{2}{\chi_{ts}}}{\rm \pi}\int_{0}^{\zeta_{t}} \dfrac{{\rm Im}\{\mathbb{P}_{ts}(-{\rm i}z)\}}{z^{2}+\sinh^{2}{\chi_{ts}}}\sin \left(z\dfrac{{\bar{x}}}{\sinh{\chi_{ts}}}\right){\rm d}z+\dfrac{{\rm \pi}_{t}\varOmega_{t}}{1+{k}^{rot}_{t}} \sin{({k}^{rot}_{t}\bar{x})}, \end{gather}$$

(4.12b)$$\begin{gather}{\bar{v}}_{t}^{rot}(\bar{x})={-}\dfrac{2\varOmega_{t}\sinh{\chi_{ts}}}{\rm \pi}\int_{0}^{\zeta_{t}}\dfrac{{\rm Im}\{\mathbb{P}_{ts}(-{\rm i}z)\}}{z^{2}+\sinh^{2}{\chi_{ts}}}\cos{\left(z\dfrac{{\bar{x}}}{\sinh{\chi_{ts}}}\right)} z \,{\rm d}z-\dfrac{{\rm \pi}_{t}\varOmega_{t}{k}^{rot}_{t}}{1+{k}^{rot}_{t}}\cos{({k}^{rot}_{t} \bar{x})}, \end{gather}$$

for the streamwise and lateral velocities associated with the transmitted-shock side, respectively. It is readily seen that SAE-induced rotational disturbances do not contribute to the velocity field at the interface $\bar {x}=0$, and therefore, they have no effect on the growth rate.

4.3. Asymptotic velocity field

The asymptotic velocity field depends on the shock oscillation frequency that ultimately determines the cases (i), (ii) and (iii) commented before. As per the former, as the shocks gradually attain planarity over extended periods, the acoustic contribution diminishes, resulting in an incompressible velocity field for $\tau _{j}\gg 1$. Consequently, (4.1) simplifies to

(4.13)\begin{equation} \nabla^{2}\boldsymbol{{\bar{v}}}_{j}^{\infty}={-}\boldsymbol{\nabla}\times\boldsymbol{{\bar{w}}}_{j}, \end{equation}

which is similar to that for the rotational field, albeit with the additional consideration of accounting for the history of the acoustic perturbation influence, which manifests as the homogeneous component of the solution to the differential equations. Consequently, the system of differential equations that must be solved for each side is written as follows:

(4.14a)$$\begin{gather} \dfrac{{\rm d}^{2}\bar{u}^{\infty}_{j}}{{\rm d}\kern0.7pt {\bar{x}}^{2}}-\bar{u}^{\infty}_{j}={-}{\bar{w}}_{j}, \end{gather}$$

(4.14b)$$\begin{gather}\dfrac{{\rm d}\bar{u}^{\infty}_{j}}{{\rm d}\kern0.7pt {\bar{x}}}={-}v_{j}^{\infty}. \end{gather}$$

The general solution for the streamwise and transverse asymptotic velocity field in the transmitted side can be written as

(4.15a)$$\begin{gather} {\bar{u}}_{t}({\bar{x}})={\bar{u}}_{ti}^{\infty}\mathrm{e}^{-{\bar{x}}}+{\bar{u}}_{t}^{rot}({\bar{x}}), \end{gather}$$

(4.15b)$$\begin{gather}{\bar{v}}_{t}({\bar{x}})={-}{\bar{u}}_{ti}^{\infty}\mathrm{e}^{-{\bar{x}}}+{\bar{v}}_{t}^{rot}({\bar{x}}), \end{gather}$$

where the rotational contribution corresponds to the particular solution. The homogeneous part, involving the exponentially decaying term, can be calculated via the final value theorem,

(4.16)\begin{equation} {\bar{u}}_{ti}^{\infty}=\lim_{\tau_{t}\to\infty}{\bar{u}}^{ac}_{t}(\tau_{t})=\lim_{s_{t}\to0}s \mathbb{U}^{ac}_{t}(s,0). \end{equation}

Notice that this term, commonly referred to as the asymptotic growth rate in the RMI context, is pivotal and will be thoroughly examined in § 5.

Similarly, the reflected field can be articulated as

(4.17a)$$\begin{gather} {\bar{u}}_{r}({\bar{x}})={\bar{u}}_{ri}^{\infty}\mathrm{e}^{{\bar{x}}}+{\bar{u}}_{r}^{rot}({\bar{x}}), \end{gather}$$

(4.17b)$$\begin{gather}{\bar{v}}_{r}({\bar{x}})={\bar{u}}_{ri}^{\infty}\mathrm{e}^{{\bar{x}}}-{\bar{v}}_{r}^{rot}({\bar{x}}), \end{gather}$$

where the homogeneous term involves ${\bar {u}}_{ri}^{\infty }={\bar {u}}_{ti}^{\infty }/\beta$.

As discussed in Wouchuk (2001) and Cobos-Campos & Wouchuk (2016), an alternative method exists for determining the asymptotic velocity, particularly the value of ${\bar {u}}_{ji}^{\infty }$. This involves substituting the independent variable ${\bar {x}}$ in (4.14a) with the shock position, denoted as ${\bar {x}}=\eta _{j}\sinh {\chi _{js}}$. Upon transformation in the Laplace variable, an expression in the Laplace variable for the asymptotic normal velocity can be found, namely

(4.18)\begin{equation} {\mathbb{U}}_{j}^{\infty}=\dfrac{s_{j} {\bar{u}}_{ti}^{\infty}-\sinh^{2}{\chi_{js}{\bar{v}}_{ti}^{\infty}}}{s_{j}^{2}-\sinh^{2}{\chi_{js}}}-\dfrac{{\bar{w}}_{j}\sinh^{2}{ \chi_{js}}}{s_{j}^{2}-\sinh^{2}{\chi_{js}}}, \end{equation}

for both the transmitted and the reflected sides. By focusing on the second term in the aforementioned expression, it becomes evident that it corresponds to the Laplace transform of the second term in (4.15a) and (4.17a). However, attempting to ILT the first term in (4.18) yields non-physical results due to the poles at $s_{j}=\pm \sinh {\chi _{js}}$. Thus, to obtain a regular solution, this term must be nullified at those points. This requirement for regularity provides us with a set of equations in the form

(4.19a)$$\begin{gather} {\bar{u}}_{ti}^{\infty}-{\bar{v}}_{ti}^{\infty}={\bar{w}}_{t}\sinh{\chi_{ts}}{\mathbb{P}}_{ts}\left(s_{t}=\sinh{\chi_{ts}}\right), \end{gather}$$

(4.19b)$$\begin{gather}{\bar{u}}_{ri}^{\infty}+{\bar{v}}_{ri}^{\infty}={\bar{w}}_{r}\sinh{\chi_{rs}}{\mathbb{P}}_{rs}\left(s_{r}={-} \sinh{\chi_{rs}}\right), \end{gather}$$

where ${\bar {u}}_{ti}^{\infty }$, ${\bar {u}}_{ri}^{\infty }$, ${\bar {u}}_{ti}^{\infty }$ and ${\bar {v}}_{ri}^{\infty }$ are unknown quantities referring to the asymptotic normal and lateral velocities at both sides of the contact interface. The system of equations completes with use made of the continuity of normal velocity, ${\bar {u}}_{ti}^{\infty }=\beta \,{\bar {u}}_{ri}^{\infty }$, along with the conservation of tangential momentum at large times, $\beta ^2( \bar {v}_{ti}^{\infty }-\bar {v}_{ti}^{0})=\vartheta (\bar {v}_{ri}^{\infty }-\bar {v}_{ri}^{0})$. Simple manipulation allows us to express the asymptotic normal velocity at the interface with

(4.20)\begin{align} {\bar{u}}_{ti}^{\infty}=\frac{\left(\beta {\bar{v}}_{ti}^{0}-\vartheta{\bar{v}}_{ri}^{0}\right)\beta}{\beta^{2}+\vartheta}+\frac{\left[\beta \varOmega_{t} \sinh{\left(\chi_{ts}\right){\mathbb{P}}_{ts}\left(\chi_{ts}\right)+\vartheta\varOmega_{r}\sinh{\left(\chi_{rs}\right)}} {\mathbb{P}}_{rs}\left(-\chi_{rs}\right)\right]\beta}{\beta^{2}+\vartheta},\end{align}

that renders the same value as (4.16), the latter being similar to (40) of Wouchuk & Nishihara (1997), which was further analysed in Wouchuk (2001).

If we examine the two terms in (4.20), the first term represents the initial tangential velocities at the interface, which induce torque and subsequently enhance the interfacial normal velocity, thereby promoting the growth of ripple amplitude. The second term refers to the pressure wave induced by the oscillations of the shock fronts. This term can either augment ${\bar {u}}_{ti}^{\infty }$ or subtract from it, sometimes resulting in a freeze-out at the interface and halting the growth of interfacial instability (Fraley 1986; Mikaelian 1994; Wouchuk & Nishihara 2004; Wouchuk & Sano 2015). One might be tempted to assert that the calculation of the growth rate is only valid in the absence of SAE at the shock fronts, suggesting that the presence of acoustic disturbances in the long-term regime, as seen in cases (ii) and (iii), can affect the asymptotic growth rate. For case (ii), since non-decaying acoustic contributions do not reach the interface, they have no impact on the contact interface dynamics. In case (iii), where constant-amplitude acoustic disturbances do reach and cross the contact interface, the result remains the same. The mean value of this contribution is strictly zero, meaning there is no imbalance in the history of the acoustic influence for this non-decaying mode. Therefore, the asymptotic growth rate is also $\bar {u}_{ti}^{\infty }$, though the asymptotic dynamic response is $\bar {u}_{ti}^{\infty }+Q_{t}\sin (\omega _{t}\tau _{t})$. Thus, even when the interface oscillates without damping, oscillations have a mean value, which is $\bar {u}_{ti}^{\infty }$.

Figure 6 shows the asymptotic velocity fields, with figure 6(a) displaying the longitudinal components and figure 6(b) showing the transverse components. These are further divided into potential (dashed orange line) and rotational (dashed green line) contributions. Figure 6(c) presents a vector plot of the asymptotic velocity field overlaid on the two-dimensional vorticity field. The parameters are chosen to meet Yang et al. (1994) results. Dimensionless velocity disturbances on the reflected side are scaled by $\beta$ to match the normalisation of those on the transmitted shock side. This adjustment ensures that the continuity of the longitudinal contribution across the interface is easily identifiable in figure 6(a). As expected, the velocity field near the interface, which is crucial for the development of the RMI, is governed by the potential contributions, while the far-field domain is dominated by the rotational contributions. This is also noticed in figure 6(c), where the misalignment between the eddies and the vorticity peaks close to the interface appears due to the potential contribution to the velocity field.

Figure 6. Longitudinal (a) and transverse (b) asymptotic velocity fields split into the potential (dashed-orange line) and rotational (dashed-green line) contributions. Panel (c) shows a vector plot of the asymptotic velocity disturbances over the vorticity field.

In the linear regime, the evolution of the rippled interface is symmetric. However, as the distortion increases, weakly nonlinear effects become significant, leading to the emergence of distinct bubbles and spikes. The former refers to the lighter fluid penetrating into the heavier fluid, corresponding the vertical positions where the velocity peaks and $\bar {u}_{ti}>0$, namely $\bar {y}=0, \bar {y}= \pm 2{\rm \pi},\ldots$, in our reference frame. The term spike refers to the heavier fluid penetrating into the lighter fluid, corresponding the vertical positions where the condition $\bar {u}_{ti}<0$ peaks: $\bar {y}=\pm {\rm \pi}, \bar {y}=\pm 3 {\rm \pi},\ldots$. To evaluate whether the distinguished acoustic fields presented in figure 2 actually affect the evolution of the nonlinear interface, an estimation of the characteristic times must be done. From one side, linear theory can estimate the time at which nonlinear effects enter into play, $\xi _i(\tau _t)\sim \varepsilon ^{-1}$, which is found to occur for $\tau _t\sim 1/ (\varepsilon \bar {u}_{ti}^{\infty })$. On the other hand, the time when the acoustic influence vanishes is of the order of unity, see Velikovich et al. (2014) and Probyn et al. (2021). Therefore, if $\varepsilon$ is sufficiently small, nonlinearity – manifested as the distinction between spikes and bubbles – becomes significant only after the asymptotic RMI growth rate is reached. Regardless of the scenario, the combined rotational and acoustic contributions to the velocity field will affect the evolution of bubbles and spikes differently depending on the side of the interface, with SAE amplifying this difference.

5. The RMI at the contact interface

5.1. Prescriptions and approximated models

For many decades, various models have been developed to estimate the linear asymptotic growth rate. These models are referred to as prescriptions. They provide approximate but explicit formulae for evaluating the constant asymptotic growth rate. Such a formula, also proposed by Richtmyer (1960), is derived by treating the RMI as a limiting case of the Rayleigh–Taylor instability, with instantaneous rather than constant acceleration, assuming that perturbation fields are incompressible. According to his model, Richtmyer's formula should be applied using preshock parameters, namely

(5.1)\begin{equation} \bar{u}_{i}^{R}=\xi_i^{0-}\mathcal{M}_{2t}{(\mathcal{R}_{t}-1)}\mathcal{A}_0, \end{equation}

where $\xi _i^{0-}=1$ with the normalisation chosen. However, to improve agreement with numerical calculations, Richtmyer himself suggested using postshock values instead of preshock ones, resulting in the widely known formula

(5.2)\begin{equation} \bar{u}_{i}^{R_2}=\xi_i^0 \mathcal{M}_{2t}{(\mathcal{R}_{t}-1)}\mathcal{A}_2 , \end{equation}

where $\mathcal {A}_2$ is the postshock Atwood number and $\xi _i^0$ is the postshock interface ripple amplitude explicitly given in (S.2.11) of the supplementary material. For ideal-gas configurations, the impulsive models work relatively well for weak shocks and/or high adiabatic indexes, i.e. when compressible effects are negligible, see figure 1 in Velikovich et al. (2014). For heavy-to-light accelerations, Meyer & Blewett (1972) empirically determined that $\xi _i^0$ should be replaced by the average of the initial unshocked and shocked amplitudes to align with their numerical results, $\bar {u}_{i}^{MB}= (\xi _i^{0-}+\xi _i^0)\mathcal {M}_{2t}\mathcal {A}_2/2$, where $\xi _i^{0-}=1$ with our variables normalisation. More accurate formulations have been proposed, as the prescription by Vandenboomgaerde et al. (1998) that combines both Richtmyer and Meyer–Blewett expressions,

(5.3)\begin{equation} \bar{u}_{i}^{VMG}= \mathcal{M}_{2t}\left[\tfrac{1}{2}\left(\mathcal{A}_2 \xi_i^{0}+\mathcal{A}_0 \xi_i^{0-}\right)-\tfrac{1}{6}\left(\mathcal{A}_2-\mathcal{A}_0\right)\left(\xi_i^{0}-\xi_i^{0-}\right) \right]. \end{equation}

For most cases tested in the literature, the second term of the right-hand side of (5.3) appears to be very small.

In addition to Richtmyer's impulsive model, another alternative impulsive model considers two different accelerations, as developed by Wouchuk & Nishihara (1996). In this case, the growth rate corresponds to the first term of (4.20), which accounts only for the initial velocity shear deposited at the interface after the incident shock refraction,

(5.4)\begin{equation} {\bar{u}}_{i}^{WN}= \frac{\beta^2 \mathcal{M}_{2t}(\mathcal{R}_{t}-1) \xi_t^0+\vartheta \beta \mathcal{M}_{2r}(\mathcal{R}_{r}-1) \xi_r^0\beta}{\beta^{2}+\vartheta},\end{equation}

where the contribution of transverse velocities at the interface is determined by the initial shock ripple amplitudes. This two-acceleration formula is an approximation that neglects the history of the transient acoustic disturbances, but it offers improved predictions for both reflected shocks and rarefaction cases compared with earlier incompressible models (Cobos-Campos & Wouchuk 2016; Probyn et al. 2021; Li et al. 2024). The irrotational Wouchuk–Nishihara model was later extended by the same authors to include bulk vorticity, as seen in (34) of Wouchuk & Nishihara (1997), which is formally similar to (4.20) and (4.16) presented in this work. However, while these equations are analytical and exact, they cannot be expressed in a fully explicit form due to the impossibility to get an explicit expression for the corresponding Laplace transform: it requires the iterative evaluation of the functions $\mathbb {P}_{ts}$. Previously, Yang et al. (1994) developed a compressible linear model derived from the Euler equations, which allows for the calculation of the exact linear growth rate. However, unlike the posterior works by Wouchuk & Nishihara (1997) and Wouchuk (2001), the model developed by Yang et al. (1994) is not analytical and requires a numerical evaluation, except for the weak-shock limit that allows for explicitness. There is still potential for (4.20) to be further developed to achieve more accurate predictions in an explicit form. For example, Cobos-Campos & Wouchuk (2016) proposed an approximation of the second term on the right-hand side of (4.20) based on the first guess of the iterative process. Optimising this choice for the term that accounts for the compressible history of the interface can lead to very accurate predictions for the asymptotic growth rate. The functions involved in this second term can be expressed as in (A21), where $\alpha _{j}$ and $\mathbb {P}^{iso}_{js}$ are explicit, accounting for the pressure field generated by the reflected and transmitted shocks as if they were isolated. Therefore, the only non-explicit expressions are the acoustic-coupling functions $F_{t}^{+}$ and $F_{r}^{-}$, which require an iterative scheme to be fully calculated, see (A19). For weak shocks, only the first order of iteration is needed, as the acoustic coupling is negligible. However, since the effect of pressure disturbances generated by the oscillating shocks is only important very close to the shock fronts, except for conditions where SAE occurs, this approximate solution offers very good results even for finite strength shocks, as shown below.

The explicit approximated solution where the pressure field generated by the shocks is taken as in isolated-shock cases reads as

(5.5)\begin{equation} {\bar{u}}_{ti}^{ISO}= {\bar{u}}_{i}^{WN}+\frac{\vartheta\beta\varOmega_{r}\sinh{\left(\chi_{rs}\right) \mathbb{P}_{rs}^{iso}\left(-\chi_{rs}\right)-\beta^2 \varOmega_{t} \sinh{\left(\chi_{ts}\right) \mathbb{P}_{ts}^{iso}\left(\chi_{ts}\right)}}}{\beta^{2}+\vartheta}, \end{equation}

where $\mathbb {P}_{js}^{iso}$ are given in (A21), this expression gives a quick and more accurate evaluation of the asymptotic growth rate in most of the cases. Another explicit expression can be constructed that renders the most accurate results compared with all other approximations for the asymptotic growth rate. It is readily obtained by considering the first guess $[k=0]$ to the acoustic functions, so that ${\bar {u}}_{ti}^{[0]}= {\bar {u}}_{i}^{ISO} + \delta \bar {u}_i^{[0]}$, where

(5.6)\begin{align} \delta \bar{u}_i^{[0]}=\frac{\beta^2 \varOmega_{t} \sinh{\left(\chi_{ts}\right)\alpha_{t}\left(\chi_{ts}\right)F_{t}^{+,[0]} \left(2\chi_{ts}\right)+\vartheta\beta\varOmega_{r}\sinh{\left(\chi_{rs}\right)}} \alpha_{r}\left(\chi_{rs}\right)F_{r}^{-,[0]}\left({-}2\chi_{rs}\right)}{\beta^{2}+\vartheta} , \end{align}

where $F_{t}^{+,[0]}$ and $F_{r}^{-,[0]}$ are given in (A20). The price to pay for this more complicated expression is compensated by the accuracy gained, as observed in figures 7 and 8, respectively. Nevertheless it is easy to see that the solution converges to the exact value of $\bar {u}_{ti}^\infty$ as the iteration index $k$ increases: ${\bar {u}}_{ti}^{\smash {[k]}}\rightarrow \bar {u}_{ti}^\infty$, for $k\gg 1$, so the expression converges to (4.20), or the equivalent (34) from Wouchuk & Nishihara (1997), as the iteration step advances. The advantage of ${\bar {u}}_{ti}^{\smash {[k]}}$ is that it starts from a physically coherent scenario (isolated shocks), and the effects of acoustic coupling are incorporated iteratively, with the first correction being fully explicit.

Figure 7. Comparison of the different prescription and approximated models with the exact solution (solid-black line).

Figure 8. Transient evolution of the growth rate (solid-black line) alongside different prescriptions and approximated models. Each panel corresponds to the respective case in figure 2. In panel (a), circles represent the data of the solution of Yang et al. (1994).

Another approximate solution is the one given by the sequential model developed in Appendix A.5, which formally applies for conditions where $\beta ^2/\vartheta \gg 1$. For that case, the value of the asymptotic growth rate is that given by (4.20), where the pressure functions in the second term of the right-hand side are given by (A27) for ${\mathbb {P}}_{rs}$, and by the corresponding combination involving the acoustic function (A28) for ${\mathbb {P}}_{ts}$. This sequential-model solution will be denoted as ${\bar {u}}_{ti}^{SM}$.

Before examining the temporal evolution of the unstable contact interface in the linear regime, we conduct a direct comparative analysis of various prescriptions and approximated models for the long-time asymptotic growth rate, including the exact solution provided by (4.16) (or (4.20), equivalently). This analysis is illustrated in figure 7, where we present a shock-strength dependence analysis for certain preshock thermodynamic conditions of the two materials or materials. Figure 7(a) corresponds to the gas-to-air SF$_6$ case, as shown in figures 2(a) and 2(b). Results are aligned with Li et al. (2024), where they find, for the heavy-to-light case, that ‘the incompressible linear model loses validity when the incident shock is strong’. Figures 7(b) and 7(c) illustrate the shock-strength variation for the cases computed in figures 2(c) and 2(d), which are associated with the vdW EoS. Finally, figure 7(d) presents a shock-strength dependence study for Al–Cu, as depicted in figures 2(e) and 2(f), where the values for $\mathcal {M}_{i}=1.3$ and $\mathcal {M}_{i}=2.5$ correspond to the linear growth rate of the nonlinear ‘light-to-heavy’ configuration addressed numerically by Tahir et al. (2011); see § 5.3 for a more elaborated discussion.

Through direct inspection, we can draw both general and specific conclusions. First, as expected, all solutions converge to the linear shock-strength dependence, $\bar {u}_{ti}^\infty \sim \mathcal {M}_i$, in the weak-shock limit $\mathcal {M}_i-1\ll 1$. Second, the approximated solution ${\bar {u}}_{ti}^{\smash {[0]}}$ although complex in form, is fully explicit and is the most accurate across the entire range of parameters explored with various EoS. In some cases, the discrepancies are nearly imperceptible within the order-of-unity scale range. Third, the simplest impulsive-model prescription, based on preshock properties, consistently renders the least accurate predictions. Further conclusions are more specific. The approximation ${\bar {u}}_{ti}^{WN}$ is generally the second-best, except for ideal gases in relatively strong-shock conditions, where ${\bar {u}}_{ti}^{SM}$ and ${\bar {u}}_{ti}^{R}$ provide slightly better approximations. The prescription ${\bar {u}}_{ti}^{VMG}$, commonly applied to the reflected rarefaction scenario, is demonstrated here only under conditions where $\beta < 1$ (see figure 7b,c). The results indicate that the accuracy of prescriptions varies depending on the specific conditions, with no definitive trend favouring one model over others. This variability suggests that caution should be taken when applying them, particularly in the context of non-ideal EoS.

Upon close inspection of the inset in figure 7(d), which corresponds to the three-term EoS for Al–Cu, it can be observed that the slope of the weak-shock linear dependence varies depending on the model being assessed. This variation arises because the three-term EoS lacks sufficient accuracy in the weak-shock limit, where the state functions fail to converge to the imposed cold conditions, such as the reference sound speeds, $c_{0r}$ and $c_{0t}$. The sensitivity of the cold-condition approach depends on specific parameters, such as $\mathcal {R}$ or $\mathcal {M}_2$, resulting in each model being affected differently in this weak-shock regime. However, this is not a major concern, as the weak-shock limit holds minimal relevance in solid RMI.

5.2. The exact transient model

We have observed that prescriptions and approximated models perform relatively well under weak-shock conditions. However, significant deviations from the exact asymptotic solution occur at finite shock strengths. Since the RMI pertains to a transient growth where nonlinear effects can become significant as the contact interface amplitude increases, it is crucial to include the exact temporal evolution in the study.

The exact evolution of the ripple amplitude, (2.6), can be derived from the integration of the velocity field evaluated at the interface, since we are interested in the point where the amplitude is greater the growth rate of the interface ripple is given by

(5.7)\begin{equation} \dfrac{{\rm d}\xi_{i}}{{\rm d}\tau_{t}}={\bar{u}}_{ti}(\tau_{t}), \end{equation}

taking the Laplace transform in time of the above equation the solution take the form of

(5.8)\begin{equation} \mathbb{S}_{i}=\dfrac{\xi_{i}^{0}}{s_t}+\frac{\mathbb{U}^{ac}_{t}(s_{t},0)}{{s_t}}, \end{equation}

where $\xi _{i}^{0}$ is the dimensionless postshock amplitude at $\tau _t=0^+$, provided in (S.2.11) of the supplementary material. With use made of the ILT, along with the convolution theorem for the second term, we obtain

(5.9)\begin{equation} \xi_{i}(\tau_{t})=\xi_{i}^{0}+\frac{2}{\rm \pi}\int_{0}^{\zeta_t}\mathrm{Im}\left\{\frac{\mathbb{L}_{t}(\mathrm{i}z,0)} {\mathrm{i}z}\right\}\frac{1-\cos(\tau_{t}z)}{z}\,{\rm d}z, \end{equation}

allowing us to write the temporal evolution of the ripple amplitude with an integral expression, since $\mathbb {L}_{ti}={\mathbb {L}}_{t}(s_{t},0)$ is already determined. However, for sufficiently large times $\tau _t\gg 1$, the asymptotic linear growth rate can be used to write

(5.10)\begin{equation} \xi_{i}(\tau_{t}\gg1)=\xi_{i0}^\infty + \bar{u}_t^\infty \tau_{t}, \end{equation}

where the offset

(5.11)\begin{equation} \xi_{i0}^{\infty}=\xi_{i}^{0^+}-\left.\dfrac{\partial\mathbb{L}_{ti}}{\partial s_{t}}\right|_{s_{t}=0}, \end{equation}

is given through the Taylor series expansion of (5.8) in the variable $s_{t}$. The value of the second term is obtained by taking the derivative of the functional expression $\mathbb {L}_{ti}$ evaluated at $s_{t}=0$, as done in Cobos-Campos & Wouchuk (2016).

The above expressions are valid for cases (i) and (ii). In case (iii), where SAE disturbances reach and cross the contact interface, the result for the ripple evolution is similar but includes a new oscillatory term. Specifically, it can be readily found through straightforward construction that

(5.12)\begin{equation} \xi_{i}(\tau_{t})=\xi_{i}^{0}+\frac{2}{\rm \pi}\int_{0}^{\zeta_t}\mathrm{Im}\left\{ \frac{\mathbb{L}_{t}(\mathrm{i}z,0)}{\mathrm{i}z}\right\}\frac{1-\cos(\tau_{t}z)}{z}{\rm d}z + \frac{Q_t}{\omega_t}\left[1-\cos(\omega_t \tau_t )\right], \end{equation}

i.e. an oscillatory term must be included. Likewise, we obtain

(5.13)\begin{equation} \xi_{i}(\tau_{t}\gg 1)=\xi_{i0}^\infty + \bar{u}_t^\infty \tau_{t} + \frac{Q_t}{\omega_t}\left[1-\cos(\omega_t \tau_t )\right] \end{equation}

for the asymptotic long-time regime.

To analyse the transient evolution of the contact interface, we perform calculations of the interface dynamics under the same conditions as in figure 2. The results, presented in figure 8, indicate that in most cases, the known expressions are not suitable for predicting the RMI, in agreement with the discussion following figure 7. Additionally, it is worth noting the multifrequency character of the solution associated with the acoustic coupling with the oscillating reflected and transmitted shocks. For example, the solid line in figure 8(a) shows the interface evolution in the linear regime, along with the circle symbols corresponding to the case computed in Yang et al. (1994). Given that the application of the ILT in this specific case is novel, we believe the Yang results, based on the same physical assumptions as our model but integrated with a different numerical method, are well-suited for benchmarking our model. Note that in our figure, the vertical axis represents the dimensionless interface ripple growth rate and incorporates a factor of $\mathcal {M}_i c_{0r}/c_{2t}=3.075$ compared with that in figure 7 of Yang et al. (1994). Meanwhile, the horizontal axis, associated with the temporal scale, is the inverse. Figure 8(b), depicting a similar air-SF$_6$ configuration but with an stronger incident shock, demonstrates similar qualitative dynamics. In contrast, figure 8(c,d), which correspond to cases using the vdW EoS, exhibit noticeable differences in the frequency and amplitude of oscillations. These differences arise from the distinct characteristics of the two configurations: figure 8(c) depicts a neutrally stable shock with an oscillating regime that does not decay over time, while figure 8(d) illustrates a postshock situation where $\beta <1$. For this case, is the reflected shock the one producing stronger oscillations, as seen in figure 2(d). Please note that we have not included the line corresponding to Richtmyer's postshock formula, as it yields a negative value ($\bar {u}_{i}^{R_2} = -0.0297$). The curve $\bar {u}_{i}^{SM}$ has also been omitted for the similar reasons. Finally, figure 8(e,f) pertain to Al–Cu configurations under different shock strengths. The stronger-shock case results in a much smoother curve with significantly lower oscillations, consistent with the pressure field results shown in figures 2(e) and 2(f). A more detailed explanation of the effect of the EoS modelling for Al–Cu configurations is provided below.

5.3. The impact of the EoS choice

In the previous subsection we have analysed the accuracy of the approximated models or prescriptions when compared with the exact solution, regardless of the EoS employed in the analysis. The present work, which is accommodated to employ any EoS, can be used to analysed whether $\gamma$-like assumptions are sufficiently accurate, see Holmes et al. (1999), Ishizaki & Nishihara (1997) and Endo et al. (1995), where a fitted or effective value of $\gamma$ and initial pressure $p_0$ in an ideal gas EoS are selected to match the post shock values of pressure $p_2$ and velocity $U_2$ observed in real experiments, or Ward & Pullin (2011), where the authors carry out a dimensional analysis in order to find complete similarity between Mie–Grüneisen and perfect gas EoS modelled flows.

To conduct this analysis, we examine the case of RMI in an Al–Cu system, similar to the numerical results presented by Tahir et al. (2011) (case II), which also correspond to the cases shown in figure 2(e,f). The incident shock Mach numbers associated with incident shock pressures of approximately $p_1= 0.28$ Mbar and $2.9$ Mbar are $\mathcal {M}_{i}\sim 1.3$ and $\mathcal {M}_{i}\sim 2.5$, respectively, according to the three-term EoS model. Although the base-flow conditions are detailed in figure S3 of the supplementary material, here we analyse the ‘equivalent’ ideal gas problem. Due to the impossibility of finding an initial preshock combination of parameters for $\gamma$-like EoS that exactly matches the base-flow results with more complex EoS (as the three-term description used in this work), i.e. due to the lack of Riemann problem similitude for different EoS (Quartapelle et al. 2003), we adopt an alternative approach. In particular, the procedure for setting up effective initial parameters for perfect-gas-like problems involves finding the best fit for both EoS zero-order profiles. The method is as follows: first, the preshock pressure $p_0$ and the adiabatic index $\gamma _r$ of the fluid through which the incident shock travels are determined by the density and pressure jumps across the shock. Then, the effective adiabatic indices $\gamma _t$, along with $\mathcal {M}_r$ and $\mathcal {M}_t$, are calculated to match the conservation of material velocity at the interface, given a fixed preshock Atwood number and the density jumps across the reflected and transmitted shocks. This procedure has been found to be very effective, at least for this Al–Cu combination considered in this work, as demonstrated in table 1, where all the base-flow density jumps across the shocks, and the postshock Atwood number and dimensionless interface corrugation are fairly similar between three-terms EoS and $\gamma$-like models. Specifically, for the Al–Cu case with $\mathcal {M}_i=1.3$, the equivalent postshock pressure conditions approximate those of an ideal gas with effective values of $\gamma _r=\gamma _{\rm Al}=3.64$, $\gamma _t=\gamma _{\rm Cu}=5.44$, with a preshock pressure $p_0=0.1366$ Mbar and the same initial density jump at the interface, characterised by $\mathcal {A}_0=0.53436$. Similarly, for the Al–Cu case with $\mathcal {M}_i=2.5$, the effective values are $\gamma _{\rm Al}=2.52$, $\gamma _{\rm Cu}=3.14$ and $p_0=0.3382$ Mbar.

Table 1. Matched initial preshock parameters between the three-terms EoS for Al–Cu cases for $\mathcal {M}_i=1.3$ and $2.5$ and their equivalent ideal gas ones. Preshock Atwood number is $\mathcal {A}_0=0.53436$. Dimensional pressures quantities are in measured in Mbar.

The exact temporal evolution of the interface ripple velocity and its corresponding asymptotic value are shown in figure 9 using the three-term EoS (black) and the ideal-gas simplification (red) for the cases represented in table 1. Precisely, in the cases presented in figure 9(a,b), the ideal gas model underestimates the growth rate by 8 % and overestimates it by 28 %, respectively. Additionally, for figure 9(a), the ideal gas model predicts a slightly higher frequency and amplitude in the ripple velocity time evolution compared with the three-terms model. In figure 9(b), the frequency difference is minimal, but the amplitude increases because the temporal evolution of the three-terms model becomes quasiplanar after the velocity reaches its asymptotic value. Moreover, if a prescription for the asymptotic growth rate is used in addition to the simplification of the EoS, the discrepancies are expected to increase further, as demonstrated in figure 8.

Figure 9. Temporal evolution of the interface ripple velocity and its asymptotic value with the three-terms EoS for Al–Cu cases – $\mathcal {M}_i=1.282$ in (a) and $2.532$ in (b) – and their equivalent ideal gas results using the effective preshock parameters listed in table 1.

Note that, although Tahir et al. (2011) examines a similar scenario of a shock incident on an Al–Cu interface, a direct comparison between the asymptotic growth rate and the numerical results there shown is not feasible due to differences in the regimes analysed. In particular, the numerical set-up in their study corresponds to a fully nonlinear regime where, in addition, elastic–plastic effects are considered. It is there shown that, for weak shocks, the interface grows from $\psi _i^{0^-}\sim 25\ \mathrm {\mu }{\rm m}$ to $\psi _i(t^*)\sim 77\ \mathrm {\mu }{\rm m}$, in a time scale that ranges $t^*\sim 800$ ns. In our dimensionless variables, this corresponds to $\xi (\tau ^*)\sim 3$, where the dimensionless time domain between frames corresponds to $\tau ^* \sim 270$, as $\lambda \sim 100\ \mathrm {\mu }{\rm m}$ and $c_{2t}\sim 5.5\ {\rm km}\ {\rm s}^{-1}$ for the weak shock case ($\mathcal {M}_i\sim 1.3$). Therefore, even if the initial amplitude were significantly smaller – which is not the case, as $\varepsilon \sim 1$ – the long-time behaviour would likely correspond to a nonlinear regime in practical conditions. The strong shock case considered in Tahir et al. (2011) pertains to the nonlinear regime too. Motivated by the lack small-amplitude studies, this work aims to address such shortage of exact verification solutions in RMI in solids. To facilitate this, we provide, not only an analytic description detailed here and in the supplementary material, but also Mathematica notebooks for both base-flow and perturbation problems. Note that, although the linear RMI problem presented here cannot be compared with the results in Tahir et al. (2011), the base-flow solutions do allow for some direct comparisons, as discussed in the supplementary material.

6. Conclusions

This work advances the well-studied problem of the RMI in the linear regime by presenting an analytical, fully compressible theory for RMI involving reflected shocks with arbitrary EoS. The key advancements are that on one hand, the initial value problem perturbation model is solved using the ILT. In contrast to the previously used separation of variables technique, which describes the pressure field in terms of Fourier–Bessel series, the method proposed in this work offers accurate solutions even under varying conditions and over long-time regimes. On the other hand, this study advances previous work by not being restricted to any particular EoS. Specific theoretical computations are conducted for an ideal gas, a vdW gas and three-term constitutive equations for simple metals. The problem formulation also includes the possibility of SAE for any of the reflected or transmitted shocks and the study of its influence on the compressed material and, eventually, on the contact interface. Additionally, an explicit expression for the asymptotic growth rate is presented, demonstrating excellent agreement across a wide variety of conditions, including variations in incident shock intensity and types of EoS. Two major conclusions are drawn, as follows.

1. (i) When comparing heuristic prescriptions with the exact model for the linear evolution of the interface, it is found that the former have a very narrow limit of validity when the shock is not sufficiently weak, regardless of the choice of EoS.

2. (ii) When comparing realistic EoS with $\gamma$-like EoS, artificially constructed to mimic zeroth-order conditions, it is found that such a EoS simplification offers low accuracy concerning the asymptotic growth rate and the characteristic dynamical properties of the transient evolution. Differences increase with the shock strength. Such comparison has been done with use made of the exact linear model.

3. (iii) The combined use of heuristic simplifications and $\gamma$-like EoS models rarely yields accurate results for the linear evolution of the interface. To reliably predict the contact interface dynamics and growth rates for the RMI, it is crucial to develop precise analytical solutions and realistic EoS models. As shown in this work, achieving this requires extensive algebraic calculations.

This study opens several opportunities for extensions and potential improvements that the authors consider highly significant. First, the heavy-to-light configuration should be addressed in similar terms. Second, for both the light-to-heavy and heavy-to-light cases, the effect of a non-ideal EoS in the nonlinear regime can be investigated using the incompressible model (see Velikovich & Dimonte (1996), Zhang & Sohn (1999), Matsuoka, Nishihara & Fukuda (2003b), Matsuoka, Nishihara & Fukuda (2003a), Velikovich et al. (2014), Matsuoka & Nishihara (2020) and references therein). Third, the use of more complex materials as those undergoing elastic–plastic regimes (Piriz et al. 2008). Finally, multimode and non-periodically symmetric configurations should be explored to better simulate realistic scenarios (see Groom & Thornber (2019), Thornber et al. (2019), Probyn et al. (2021) and Dimonte et al. (2024)).