## Abstract

A theoretical framework to study linear and nonlinear Richtmyer–Meshkov instability (RMI) is presented. This instability typically develops when an incident shock crosses a corrugated material interface separating two fluids with different thermodynamic properties. Because the contact surface is rippled, the transmitted and reflected wavefronts are also corrugated, and some circulation is generated at the material boundary. The velocity circulation is progressively modified by the sound wave field radiated by the wavefronts, and ripple growth at the contact surface reaches a constant asymptotic normal velocity when the shocks/rarefactions are distant enough. The instability growth is driven by two effects: an initial deposition of velocity circulation at the material interface by the corrugated shock fronts and its subsequent variation in time due to the sonic field of pressure perturbations radiated by the deformed shocks. First, an exact analytical model to determine the asymptotic linear growth rate is presented and its dependence on the governing parameters is briefly discussed. Instabilities referred to as RM-like, driven by localized non-uniform vorticity, also exist; they are either initially deposited or supplied by external sources. Ablative RMI and its stabilization mechanisms are discussed as an example. When the ripple amplitude increases and becomes comparable to the perturbation wavelength, the instability enters the nonlinear phase and the perturbation velocity starts to decrease. An analytical model to describe this second stage of instability evolution is presented within the limit of incompressible and irrotational fluids, based on the dynamics of the contact surface circulation. RMI in solids and liquids is also presented via molecular dynamics simulations for planar and cylindrical geometries, where we show the generation of vorticity even in viscid materials.

## 1. Introduction

The Richtmyer–Meshkov instability (RMI) was predicted by Richtmyer more than 50 years ago; since then it has attracted the attention of scientists worldwide because of its importance in different scenarios (Richtmyer 1960; Meshkov 1969; Fraley 1986; Holmes *et al.* 1999; Zabusky 1999; Velikovich *et al.* 2000, 2007). Richtmyer studied the problem numerically for a planar shock refracted normally at a contact surface separating two different ideal gases, assuming that the material interface had an initial sinusoidal corrugation. He developed the exact linear perturbation analysis of the perturbed problem between the reflected and transmitted shock fronts. He found that the initial ripple would grow in an oscillatory way until a constant asymptotic normal velocity is achieved when both wavefronts are at sufficient distance from the contact surface. Richtmyer (1960) considered only the case in which a shock is reflected. A decade later, Meyer & Blewett (1972) made a numerical analysis of the case in which a rarefaction is reflected after the incident shock interacts with the material surface.

In his pioneering work, Richtmyer provided a heuristic formula with which to estimate the asymptotic rate of growth at the surface ripple, known as the impulsive approximation. Namely, the effect of the incident shock upon the contact surface ripple might be equivalent to the action of an impulse acceleration of very short duration. Later, Meyer and Blewett generalized that formula to handle a reflected rarefaction, and derived an ad hoc prescription based on the results inferred from their numerical analysis. So far, the impulsive approximation is known to be valid only for weak incident shocks (Yang *et al.* 1994; Velikovich 1996; Wouchuk & Nishihara 1996, 1997; Velikovich *et al.* 2000; Wouchuk 2001*a*,*b*).

The hydrodynamic flows associated with this instability are important in all physical environments in which corrugated shocks/rarefaction fronts are the vehicles that transfer energy to the material. RMI-like flows appear, for example, in astrophysics (Andreopoulos *et al.* 2000) and in the shock remnants following a supernova explosion. It is also an important phenomenon in inertial confinement fusion (ICF), as it invariably develops during the early stages of target irradiation, when a sequence of shocks travels through the target. As the shock fronts are not uniform, their corrugations will induce a highly perturbed flow behind them, which will be the seed for the Rayleigh–Taylor instability (RTI) that occurs later during the acceleration phase of the target (Ishizaki & Nishihara 1997; Piriz *et al.* 1997; Nishihara *et al.* 1998; Goncharov 1999; Velikovich *et al.* 2000). The RMI has also been thoroughly studied in experiments designed in shock tubes (Jones & Jacobs 1997; Zabusky 1999; Brouillette 2002) and has also become a fascinating research topic in the field of high-energy density experiments in matter. Furthermore, this type of perturbation flow may be used to our advantage by creating those perturbed flows artificially in the laboratory to excite the related instabilities and using them as a tool to extract information from matter at extreme conditions of pressure, density and temperature. (e.g. the possibility of structural and phase transformations in solid materials (Mikhailov 2007)).

Because it is pervasive in many important subfields within the limits of the validity of the physics of fluids and plasmas, the need for an accurate understanding of the mechanisms that drive the flow dynamics in such situations is clear.

When dealing with RMI-like flows, the compressibility of materials at the high pressures and densities achieved in experiments becomes an important factor that cannot be ignored. Otherwise, the information that we could extract from experiments will be limited to narrow ranges of shock strength, precluding the exploration of interesting regimes achievable at higher compressions. Compressibility in this type of problem manifests itself as the essential ingredient for describing the interaction between the corrugated shocks and the boundaries that limit the domain of interest. This interaction is mediated by sound waves. Besides compressibility, other interesting physical phenomena, for example shocks propagating into viscous fluids (Miller & Ahrens 1991), the influence of elastic/plastic properties of solid materials (Zybin *et al.* 2006; Piriz *et al.* 2008) and the effect of magnetic fields on the dynamics of corrugated shocks in plasmas (Wheatley *et al.* 2005), are certainly worth studying. In fact, when dealing with more complex properties and states of matter, as well as the inclusion of magnetic fields for shocks propagating in plasmas, the interaction between the shock fronts and boundaries may be sustained not only by the longitudinal adiabatic sonic waves, but with another family of wave fronts that propagate at speeds other than the speed of sound in ordinary gases. These new problems would certainly provide an exciting playground to test the current understanding of RMI flows, widening even more the possibilities of exploring the behaviour of matter under extreme conditions, challenging our comprehension of basic physics and stimulating the development of additional knowledge. We are convinced that the study of the perturbation dynamics that develops behind corrugated shocks/rarefaction fronts is a fascinating topic in itself, and also that it will still challenge theoreticians and experimenters alike for years to come.

An impulsive model for quasi-incompressible flows was set up in Vandenboomgaerde *et al.* (1998). Viscosity effects have been dealt in Carles & Popinet (2001) and Griffond (2005). The initial phase of the RMI has also been investigated in Griffond (2006). Short distance shock–interface interaction has been recently investigated in Lombardini & Pullin (2009).

In the analytical models in this short review, however, we focus on a much humbler problem, limiting ourselves to the study of the linear and nonlinear evolution of RMI in inviscid fluids, within the realms of gas dynamics. We concentrate on ideal gases, as they represent the simplest model that can be used to minimize mathematical complications. However, it is necessary to emphasize that the important results shown here are in no way restricted by the use of an ideal gas equation of state (EOS). Note, however, that the results of molecular dynamic (MD) simulations will be briefly presented for solids and liquids showing such RMI-induced phase transitions.

The RMI is essentially driven by the vorticity deposited at a material interface in a shock–interface interaction. In linear theory, the initial circulation is modified by the pressure waves emitted by the corrugated shocks. Additional vorticity is generated in the bulk of the fluids if the shock is strong and the fluids are highly compressible (Wouchuk 2001*a*). The vorticity localized at the interface can be deposited in many other ways. This is how we generalize RM-like instabilities (Velikovich *et al.* 2000). They are defined as interfacial instabilities driven by the same physical mechanisms as RMI in shock–interface interaction, but differing in the source of the unstable configuration, and sometimes as well by some other factors, affecting perturbation evolution. It should be noted that, in RMI and RM-like instabilities, the vorticity is non-uniformly distributed along the interface, which is different from Kelvin–Helmholtz (KH) instability.

It is known that the vortices interact with each other (Saffman 1992). Namely, the interface dynamics has a non-local character, which plays an important role in RMI nonlinear evolution. The bulk vorticity behind the shocks may be important in the nonlinear phase when the interface circulation interacts with neighbouring vortices on both sides. Another important feature of RMI nonlinear evolution is creation and annihilation of the vorticity, due to the finite density jump at the interface. Therefore, the non-uniformly distributed vorticity results in complex dynamics of the interface in the nonlinear phase. Zabusky *et al.* visualized its complex dynamics with the use of hydrodynamic simulations, and coined the terms ‘vortex paradigm’ and ‘vortex projectile’ (Hawley & Zabusky 1989; Zabusky & Zeng 1998; Zabusky 1999; Zabusky & Zhang 2002). In tank experiments (Jacobs & Sheeley 1996), they also visualized complex dynamics of the interface, such as the spiral structure of a spike. It should, however, be mentioned that in this review we consider only sinusoidal perturbations as initial interface corrugations and exclude subsequent turbulent regimes often observed in shock tube experiments (Houas & Chemouni 1996; Poggi *et al.* 1998).

This review is organized as follows: in §2, the dynamics of the interaction between the rippled shock and the contact surface is briefly discussed as a preliminary to understanding the physical mechanisms that drive RMI. The two main scenarios of shocks and rarefactions reflected at the material interface are succinctly presented. An exact analytical formula valid for any values of the governing parameters (shock strength, fluid compressibility and density ratio) is very briefly explained. Comparisons with recent experiments and simulations are shown. We also discuss some RM-like instabilities, including ablative RMI (Ishizaki & Nishihara 1997; Goncharov 1999), which show stabilization of the instability owing to mass flow through the surface and dynamical overpressure in a corona.

In §3, we investigate weakly and fully (arbitrary amplitude) nonlinear evolution of the interface in RMI with planar and cylindrical geometries within a limit of incompressible fluids. Weakly nonlinear stability analysis is performed to present the growth rate of a bubble and spike at the nonlinear stage. We also present fully (arbitrary amplitude) nonlinear evolution, such as the roll-up of the interface, by the vortex method (Krasny 1987). The Atwood number dependence on the nonlinear growth rate and various differences between RMI in planar and cylindrical geometries are discussed.

Over the last decade, a great advance in computing power has created feasible opportunities for modelling hydrodynamics problems via MD simulations. With a few billions of atoms simulated, MD can grasp at the atomic level the essentials of many non-equilibrium processes in material science, shock physics, detonation and hydrodynamic instabilities such as RM and RT (Germann *et al.* 2004; Kadau *et al.* 2004, 2008) and KH (Glosli *et al.* 2007). Some recent MD simulation results of RMI in planar and cylindrical geometries are discussed in §4.

## 2. Linear perturbation growth

### (a) Rippled shock dynamics and generation of tangential velocities

We consider the situation shown in figure 1*a*. Here, an incident planar shock hits a planar contact surface separating two fluids with different thermodynamical properties at the instant of time *t*=0. We suppose that the fluids are ideal gases with adiabatic exponents *γ*_{a} and *γ*_{b}. (However, the results of this section are not restricted to an ideal gas EOS.) The incident shock velocity is in the laboratory system of reference. This shock compresses the fluid on the right from the initial density *ρ*_{b0} to the value *ρ*_{b1} as shown in figure 1*a*. The fluid velocity behind the incident shock is . Let *p*_{0} be the pressure of both fluids before shock compression and *p*_{1} be the pressure that drives the incident shock. Immediately after the incident shock is totally refracted at *x*=0, the reflected shock forms and starts moving to the right. Similarly, a transmitted shock moves to the left inside the other fluid (figure 1*b*). The reflected shock compresses the initial fluid from *ρ*_{b1} to the final value *ρ*_{bf}. Conversely, the transmitted shock compresses the fluid to the left from *ρ*_{a0} to *ρ*_{af}. As a consequence of shock refraction, the material surface acquires a velocity . The reflected front moves with velocity and the transmitted front moves with velocity in the laboratory reference frame. The incident shock Mach number with respect to the uncompressed fluid is defined as *M*_{i}=*D*_{i}/*c*_{b0}>1, where *c*_{b0} is the initial sound speed of fluid *b*. The Mach numbers of the reflected and transmitted shocks with respect to the compressed fluids are: *β*_{r}=(*D*_{r}+*U*_{i})/*c*_{bf}<1 and *β*_{t}=(*D*_{t}−*U*_{i})/*c*_{af}<1, where *c*_{af} and *c*_{bf} are the final values of the compressed sound speeds of fluids *a* and *b*. The final pressure at both sides of the contact surface is *p*_{f}.

Let us now consider the situation in which the contact surface that separates the fluids is initially corrugated. We assume the initial ripple of the material interface to be of the form . Hence, just after shock refraction at *t*=0+, the reflected and transmitted wavefronts are deformed and their initial ripples resemble that at the contact surface, but with different amplitudes (figure 2). A simple kinematic calculation shows that the initial shock ripple amplitudes are (Richtmyer 1960)2.1for the reflected and transmitted fronts, respectively. These initial ripples are the ingredients that actually drive the perturbed flow in the space between the shocks and the material interface, providing the mechanism that drives later contact surface ripple growth. As is well known, the tangential velocity must be continuous across any of the deformed shock fronts (Landau & Lifshitz 1987). For the geometry of figure 1*b*, this can be satisfied, admitting the existence of lateral perturbed velocities in each fluid, only behind either the reflected or the transmitted shock front. We call those initial tangential velocities and . Using simple algebra, their values can be found to be equal to2.2

These two tangential velocities provide an initial circulation that is distributed along the contact surface ripple proportionally to , which is actually the initial cause of the growth of the contact surface ripple at *x*=0 (Velikovich 1996; Wouchuk & Nishihara 2004; Herrmann *et al.* 2008). As the shocks separate away from the interface, their corrugation amplitude will oscillate and decrease in time. This well-known fact was predicted theoretically and observed in experiments/simulations and documented over more than 50 years of research on corrugated shocks, either in shock tubes or, recently, with high-power lasers (Briscoe & Kovitz 1968; Fraley 1986; Velikovich & Phillips 1996; Wouchuk & Nishihara 1996). For an ideal gas EOS, it can be rigorously proven that the amplitude of those oscillations decreases asymptotically in time as *t*^{−3/2} for any shock with finite intensity (Zaidel 1960; Fraley 1986; Wouchuk & Cavada 2004). In figure 3, we show the behaviour of the reflected shock ripple as a function of time for a typical RMI problem. Stable constant-amplitude oscillations due to the D’yakov (1955) and Kontorovich (1959) instability have been seen in numerical simulations in Bates & Montgomery (2000), and studied analytically in Wouchuk & Cavada (2004).

As the shock ripple evolves in time, pressure perturbations are generated exactly behind the shock front and propagate through the compressed fluid with the local sound speed (Landau & Lifshitz 1987). For the RMI problem at hand, it can be seen that the pressure fluctuations are actually evanescent waves emitted by the corrugated shock fronts, which decay in space and time. Thus, the evolution of the perturbations at the material contact surface ripple, located at *x*=0, will be influenced by the time history of the pressure fluctuations since *t*=0+, during the lapse of time in which the shocks are not far from that surface. In other words, the early time history of the pressure perturbation field in the space between the corrugated shocks is important for getting an exact description of the subsequent linear growth of the material interface ripple.

### (b) Generation of vorticity by the corrugated shock fronts

As the corrugated shocks move into the undisturbed fluid and their ripples oscillate in time, they leave behind them a wake of tangential velocity arrows with alternating signs distributed inside the fluids. It can be proved that the velocity field associated with those perturbations is actually rotational. That is, a non-zero vorticity is frozen to the compressed fluid particles, which is spread at both sides of the material surface. Of course, as the shocks separate away from *x*=0, the amplitude of the oscillations of the shocks’ ripple decreases and the intensity of the shock-generated vorticity decreases with distance as well. This ‘*bulk*’ vorticity is therefore only appreciable near the contact surface. Its effect on the instability’s evolution is more important for strong shocks and very compressible gases. It must be emphasized that the characteristic length of the vorticity field generated by the deformed shocks is not determined only by the contact surface perturbation wavelength *λ*=2*π*/*k*, it is also strongly dependent on the shock intensity (Wouchuk 2001*a*,*b*). In fact, for very strong shocks and highly compressible fluids, it can be significantly smaller than *λ*. In a reference system co-moving with the compressed fluid, behind each of the corrugated fronts, vorticity is generated at the position (*x*,*y*) at the time *t*_{0}(*x*)=|*x*|/(*c*_{mf}*β*_{s}), where *m*=*a*,*b* indicates any one of the gases, and *s*=*r*,*t* refers to the reflected or transmitted shock front. This vorticity perturbation remains frozen to that fluid element in the absence of viscosity. The amount of vorticity generated at point *x*, which we call *δω*, can be seen to be proportional to the value of the instantaneous pressure perturbation behind the shock front, at the time *t*=*t*_{0}(*x*) (Wouchuk & Nishihara 1996, 1997; Kevlahan 1997; Wouchuk 2001*a*,*b*),2.3

The proportionality constant *Ω*_{m} depends on the EOS of the gas (Wouchuk & Nishihara 1996, 1997). For an ideal gas EOS, it is for very weak shocks, but saturates at a finite non-negligible value for very strong shocks (Wouchuk 2001*a*,*b*).

### (c) Contact surface ripple: asymptotic linear rate of growth

The evolution of the velocity and pressure perturbations generated behind the corrugated shock fronts determine the initial linear growth of the contact surface ripple at *x*=0. As the reflected and transmitted shocks separate away, their ripples oscillate, generating additional vorticity in the interior of both gases and at the same time an evanescent pressure perturbation field, which extends behind each shock surface. Because the gases have an ideal gas EOS, the shocks’ ripple amplitudes tend asymptotically to zero and the shocks become virtually flat once they are several wavelengths away from the material interface. Thus, asymptotically in time, no more pressure perturbations are generated by the shock fronts in the space between the shocks and the contact interface. Furthermore, the pressure fluctuations that had been generated immediately after *t*=0+ would have decayed in the space near the surface ripple, and the velocity fields continue evolving by their own inertia, without significant pressure gradients, always within the limits of validity of the linear theory. That is, for , the pressure perturbation field is practically zero, and, owing to the conservation of mass, the velocity field becomes asymptotically incompressible, independently of the initial shock intensity. This last property is a very important consequence of the fact that the shocks become asymptotically flat.

Before proceeding to calculate the contact surface tangential and normal velocities, we must emphasize that this asymptotic incompressibility does not exclude the possibility of density perturbations remaining in the gases. In fact, while the shocks are corrugated, entropy perturbations are generated behind them; these entropic disturbances remain frozen to the compressed fluid particles just as vorticity perturbations do. The entropy fluctuations are associated with density perturbations that follow the fluid elements’ trajectories in the absence of dissipative processes. These entropic density perturbations clearly do not travel with the sound waves; they are frozen to the compressed fluid particles and their spatial dependence will be the same as that of the vorticity. In other words, in a system of reference co-moving with the compressed fluid(s), the particles preserve the entropic and vorticity perturbations generated by the deformed shock waves. It is only the perturbation field associated with the evanescent sound field that also decays asymptotically in time. Therefore, the evolution of velocity perturbations near the contact surface clearly depends on the dynamics of the shocks’ ripples through the intervention of both the sound pressure field and the vorticity generated at both sides, since *t*=0+. To examine the situation at the contact surface ripple, we need to integrate only the tangential fluid equation of motion at *x*=0. In fact, we write for fluid *m* (*m*=*a*,*b*)2.4where we have assumed that , and . This equation is valid at both sides of the contact surface at any time *t*>0+. If we use the pressure continuity at *x*=0 and integrate in the time interval , we obtain the jump relationship2.5where is the amplitude of the asymptotic tangential velocity at the surface ripple in fluid *m*. First, we realize that the final tangential velocity at the interface ripple is dependent on the integral of the pressure fluctuations at that point. The pressure fluctuations have been generated at the corrugated shock fronts, which propagate with the local speed of sound, arrive at the interface, modify the local circulation and get refracted there, bouncing between the corrugated fronts. This process of velocity variation at the interface owing to the sound wave field occurs at any shock intensity (weak or strong shocks). Additionally, as the shocks separate from the contact surface, their oscillations will induce the generation of vorticity in the bulk (Kevlahan 1997; Wouchuk & Nishihara 1997; Wouchuk 2001*a*). This bulk vorticity is negligible compared with the interface circulation for weak shocks and/or very incompressible fluids. In contrast, it cannot be neglected for situations in which the fluid’s compressibility is important, for example for very strong shocks. This circumstance has a definite effect on the normal velocity measured at the interface ripple, as will be discussed later. Our main objective in this section is to calculate the value of the asymptotic normal velocity perturbation at *x*=0, which we call . If we add and substract on each side of equation (2.5), we get2.6where and .

Equation (2.6) is the exact analytical expression for the linear asymptotic normal velocity at the contact surface ripple, valid for any values of the parameters that define the problem: fluid compressibilities, initial density ratio and incident shock Mach number. Although we have started our discussion considering ideal gases, the validity of equation (2.6) is general and holds for any EOS of the intervening fluids; i.e. to get from equation (2.6), we need to know the values of *F*_{a} and *F*_{b}. These parameters are also dependent on the pressure perturbations generated at the shock fronts. Incidentally, after some algebra, it can be shown that the quantities *F*_{m} can also be written as weighted spatial averages of the vorticity profiles generated by each corrugated shock in the corresponding fluid, i.e. it can be seen that (Wouchuk & Nishihara 1997; Wouchuk 2001*a*)2.7The velocity fluctuations emanating from *x*=0 would always decay to zero for |*x*|≫*λ*, and, at the same time, those velocity fluctuations must match the rotational velocity field in the bulk gases generated by the shocks. This constraint is guaranteed once we impose the condition in equation (2.7). The details of the arguments leading to equation (2.7) cannot be shown here because of lack of space. They can be found in Wouchuk & Nishihara (1997) and Wouchuk (2001*a*). An accurate and highly efficient method for calculating the quantities *F*_{m} is fully described in Wouchuk (2001*a*). We only add here that, for very weak shocks (*M*_{i}−1≪1), the quantities , and they can be safely neglected as discussed earlier. In fact, for weak shocks with *M*_{i}≤1.5, the first term on the right-hand side of equation (2.6) is enough to get an accurate estimate of the asymptotic normal velocity. Excellent agreement with linear numerical simulations was confirmed (Wouchuk 2001*a*). It was also seen that the asymptotic velocity can be negative for specific choices of the gas and shock parameters. This means that, in some regions of the parameter space, it is possible to have a strictly zero asymptotic normal velocity. This situation is commonly called *freeze-out* (Mikaelian 1994). The exact determination of the conditions that result in freeze-out is not a trivial problem. Those conditions are determined by requiring in equation (2.6). This requires that the first term in equation (2.6) cancels out the second term. That is, in freeze-out situations, the bulk vorticity generated by the corrugated shocks can be matched with the circulation at the interface ripple only if the normal velocity is zero.The details of the corresponding calculations can be found in Wouchuk & Nishihara (2004).

### (d) Rarefaction reflected case

For some choices of the incident shock Mach number (*M*_{i}), gas compressibilities (*γ*_{a}, *γ*_{b}) and initial fluid density ratio (*R*_{0}=*ρ*_{a0}/*ρ*_{b0}), the outcome of the *incident shock/contact surface* interaction may actually be a centred rarefaction wave moving to the right plus a transmitted shock moving to the left. This typically happens when *R*_{0}=*ρ*_{a0}/*ρ*_{b0}<1. The conditions that might give rise to either a reflected shock or a reflected rarefaction can be found in Yang *et al.* (1994) and Velikovich & Phillips (1996). In the general case, a shock impinging towards a lower density fluid will produce a reflected rarefaction travelling back into the denser fluid. The rarefaction that moves to the right has a more complicated spatial structure than a discontinuous shock-compressing front. Unlike the shock wave, a rarefaction wave is an isentropic expansion of the fluid elements; the expansion zone is limited by two fronts: a rarefaction head and a rarefaction tail, each moving with the local sound speed with respect to the fluid. This situation is shown in figure 4*a*. The different quantities inside the rarefaction region decrease continuously from their values at the rarefaction head up to the rarefaction tail. However, the fluid bounded by the rarefaction tail and the contact surface is uniform.

There is a striking difference between the evolution of the RMI in the reflected rarefaction and reflected shock cases. In the rarefaction situation, the growth rate is negative in general. This difference can be easily understood if we look at the way the incident shock gets refracted at the contact surface. In fact, as the transmitted shock velocity is generally larger than the incident shock speed (because generally *ρ*_{a0}<*ρ*_{b0}), the phase of the initial transmitted shock ripple is opposite to that of the contact surface ripple. Besides, the tangential velocity behind the rarefaction tail is also negative sign (Velikovich 1996; Velikovich & Phillips 1996; Wouchuk & Carretero 2003). The structure of the corrugated wavefronts at *t*=0+ is shown in figure 4*b*. The tangential velocities on both sides of the contact surface corrugation have a different orientation in this case from that of the shock reflected case. This process enhances ripple growth with an inverted phase. Of course, the formula given in equation (2.6) also applies exactly in this case, with one additional simplification: the term *F*_{b}=0 because there is no vorticity between the contact surface and the rarefaction tail. This last property is a consequence of the fact that the fluid in front of the rarefaction head is trivially irrotational, and no vorticity is generated inside the expanding fluid because the rarefaction wave is isentropic (Velikovich 1996; Velikovich & Phillips 1996; Wouchuk 2001*b*; Wouchuk & Carretero 2003). Additional details for this particular case of RMI cannot be discussed in any depth here owing to lack of space. The interested reader is referred to the published literature (Yang *et al.* 1994; Velikovich 1996; Wouchuk & Nishihara 1996, 1997; Wouchuk 2001*b*; Wouchuk & Carretero 2003). In figure 5, we show a comparison with more or less recent experiments done with the Nova laser, the details of which can be found in Dimonte *et al.* (1996). The dotted line corresponds to perfect agreement between the experimentally measured values and the theoretical prediction of equation (2.6). The calculations based on equation (2.6) for these experiments are fully described in Wouchuk (2001*b*).

### (e) RM-like instabilities and ablative RMI

We have shown in previous sections how the vorticity is deposited at an interface in a shock–interface interaction, which drives RMI. The vorticity localized at the interface can be deposited in many other ways. It can either be initially deposited or supplied by external sources, and the initially deposited vorticity can be affected by the rippled shock(s) through sound waves. A good example of vorticity deposited at the interface is the interfacial instability between immiscible fluids driven by a free-falling tank bouncing off a fixed spring for a short interval (Jacobs & Sheeley 1996). They give a sinusoidal interface perturbation, for example , between two fluids in the tank, using novel experimental techniques. During free-fall time intervals, no acceleration is felt by the fluids in the tank frame of reference. However, in the interval of time in which the tank makes contact with the spring, the system suffers an impulse that changes the velocity. Namely, vorticity is deposited at the interface during this short but finite contact time. The contact duration can be determined from the spring constant *K*, the mass of the container *M*, its initial velocity at contact *v*_{0} and the earth gravity *g*_{0}, as *t*_{f}=(*π*+2*δ*)/*ω*_{0}, where and . The motion of the interface perturbation is given as2.8and2.9according to Jacobs & Sheeley (1996), where *A* is the interface Atwood number, *A*=(*ρ*_{h}−*ρ*_{l})/(*ρ*_{h}+*ρ*_{l}), and *ρ*_{h} is the density of the heavier fluid and *ρ*_{l} that of the lighter fluid. By integrating the equation of motion, we calculate the normal velocity of the interface from the relation of *u*_{n0}=*dξ*/*dt* at time *t*_{f}. Since both fluids are incompressible, ∇⋅** u**=0, the fluid velocity at time

*t*

_{f}can be estimated as in the fluids. The interaction of the tank with the spring deposits a sinusoidal vorticity at the interface corresponding to the initial interface corrugation. Nonlinear dynamics of the unstable interface will be discussed in §3.

As another example of RM-like instability, we discuss the stability of a laser ablation surface. Laser ablation surfaces have characteristics similar to those of flames or a deflagration wave from a hydrodynamics point of view (Bobin 1971; Takabe *et al.* 1978), although electron thermal conduction plays an important role, as will be discussed later. A flame, acting on a gas ahead in a manner similar to that of a piston, can drive a shock wave (Landau & Lifshitz 1987). A laser ablation also acts as a piston and drives a shock wave propagating through the target, although some mass flow exists across the surface (called mass ablation). In laser–plasma interactions, a laser light can penetrate into a plasma up to the critical density where the electron plasma frequency equals the laser frequency. The critical density is much lower than the solid density. Electrons absorb laser energy and an electron heat wave propagates into the high-density region. The heat wave steepens due to electron nonlinear thermal conductivity , where *T*_{e} is the electron temperature. The electron energy is deposited to ions near the heat wavefront at high density, where the electron–ion energy exchange time, proportional to *T*_{e}^{3/2}/*n*, becomes very short. Here, *n* is the plasma density. This energy deposition in laser ablation corresponds to that due to chemical reactions in a flame. The heated ions are then accelerated to vacuum through the ablation layer. In the approximation of steady ablation, momentum flow is conserved through the ablation layer. Static pressure, called ablation pressure, dominates at the ablation surface, since the mass flow across the surface is subsonic.

When laser light is irradiated on a target with a corrugated surface, a rippled shock wave is launched corresponding to the corrugated surface. The interaction between the rippled shock front and the ablation surface through sound waves may induce RM-like instability. It is important to study such growth of RM-like instabilities because they determine the seed of the RT instability that develops during the acceleration phase in ICF (Emery *et al.* 1991; Desselberger *et al.* 1995; Endo *et al.* 1995; Azechi *et al.* 1997; Ishizaki & Nishihara 1997, 1998; Taylor *et al.* 1997; Nishihara *et al.* 1998; Goncharov 1999; Matsui *et al.* 1999; Goncharov *et al.* 2000; Velikovich *et al.* 2000).

A simple analytical model was developed to study propagation of a rippled shock associated with the initial surface roughness of a target in Ishizaki & Nishihara (1997, 1998). They sketch a schematic profile of a stationary shock wave driven by a steady laser ablation, shown in figure 6. The domain can be separated into four regions using the shock front, ablation front and sonic point. These regions are labelled 0, 1, A, and 2 from right to left. Region 0 is in a uniform state ahead of the shock; region 1 is the shock-compressed region; region A is an ablation layer between the ablation front and the sonic point; region 2 beyond the sonic point is an isothermal rarefaction region. The Rankine–Hugoniot and Chapman–Jouguet deflagration jump conditions were applied at the shock and ablation fronts, respectively, in Ishizaki & Nishihara (1997, 1998). They also assumed that the fluid velocity of the sonic point relative to the ablation surface is equal to the sound speed of the sonic point and that the density of the sonic point is the critical density (Manheimer & Colombant 1984). The zeroth-order variables in each region can then be uniquely determined using these conditions and the absorbed laser intensity.

First, we consider a rippled shock wave caused by the initial surface roughness of a target. The surface modulation of the target is assumed to be given as . The wave equation for the pressure perturbation in the shock-compressed region can be solved analytically with proper boundary conditions as described earlier (Richtmyer 1960; Zaidel 1960; Briscoe & Kovitz 1968; Ishizaki *et al.* 1996; Ishizaki & Nishihara 1997, 1998). The boundary conditions at the ablation surface are obtained by linearizing the CJ jump conditions, which give the relationships between the perturbations at the ablation front and those at the sonic point (Wouchuk & Piriz 1995; Ishizaki & Nishihara 1997, 1998). The first-order quantities at the sonic point are assumed to satisfy the condition that the local Mach number is equal to unity. In addition, the density perturbation at the sonic point is assumed to be zero because the density at the sonic point is taken to be the laser cut-off density (Manheimer & Colombant 1984). Because the flow in region 2 expands supersonically, neither sound nor entropy waves can cross the ablation front and affect the flow in region 1.

Figure 7 shows the shock front ripple (solid line), *a*_{s}/*a*_{0}, as a function of and the ablation surface deformation (dot-dashed line), *a*_{a}/*a*_{0}, as a function of , where *u*_{s} and *u*_{a} are the shock and ablation velocities, respectively, and *v*_{x1} and *c*_{1} are the fluid velocity and sound speed in region 1, respectively. The parameters used are the absorbed laser intensity, *I*=4×10^{13} *W* *cm*^{−2}, the laser wavelength, *λ*_{L}=0.53 μm, the density, *ρ*_{0}=1.06 *g* *cm*^{−3} (CH target), the pressure, *p*_{0}=0.703 *Mbar* (equivalent to the temperature, *T*_{0}=1 *eV*), the perturbation wavelength, *λ*=100 *μm*, and the isentropic exponents, *γ*_{0}=3, *γ*_{1}=3 and , where the subscripts 0, 1 and 2 denote the values in regions 0, 1 and 2, respectively. Once a rippled shock is launched, a pressure perturbation is induced by lateral fluid motion behind the shock, and the pressure perturbation causes the ripple of the shock front to be reversed and subsequently oscillate, similar to the shock–interface interaction discussed before. However, compared with the rippled shock driven by a rippled rigid piston (dotted line in figure 7; Briscoe & Kovitz 1968; Ishizaki *et al.* 1996), the amplitude of the shock surface ripple driven by laser ablation decays much faster than that driven by a rigid piston. The pressure perturbation increases the deformation of the ablation front monotonously. However, the deformation amplitude of the ablation surface (not the growth rate) reaches an asymptotic amplitude, which is quite different from the case of the shock–interface interaction. Namely, ablative RMI in this case is totally stabilized. In the weak shock limit, the amplitude of the shock front ripple *a*_{s} is given by2.10where *J* is the Bessel function and *M*_{s} is the shock Mach number. This approximate formula is shown by circles in figure 7, and agrees quite well with the exact solution even for relatively large *M*_{s}(=3.74). An asymptotic value of the ablation deformation for *a*_{a} is given by2.11where *v*_{a}/*c*_{1} represents the ablation Mach number. Note that, as shown in equation (2.11), the asymptotic amplitude decreases with increasing ablation velocity. The dependence of the amplitude on the incident laser intensity can be found in Ishizaki & Nishihara (1997, 1998). Those results indicate that ablative RMI exhibits stabilization owing to mass ablation because it carries away the vorticity and pressure perturbations from the ablation surface (Ishizaki & Nishihara 1997, 1998; Taylor *et al.* 1997; Nishihara *et al.* 1998; Velikovich *et al.* 1998).

In figure 8, the theoretical values are compared with the experimental results (Endo *et al.* 1995). The parameters used are the same as those in figure 7. Both the oscillation period and decay rate of the rippled shock front agree quite well with the experimental results, as shown in figure 8.

However, the deformation of the ablation surface given by equation (2.11) is much greater than experimental results with much longer pulse duration (Aglitskiy *et al.* 2001, 2002) and with higher laser intensity (Gotchev *et al.* 2006) and simulations (Taylor *et al.* 1997) because the above-mentioned model takes into account the stabilization mechanism only due to the mass ablation. Another stabilization mechanism in ablative RMI, called dynamic overpressure, is discussed in Goncharov (1999) and Goncharov *et al.* (2000), in which the sharp-boundary model (SBM) was used instead of the deflagration model. The SBM was originally developed to study ablative RTI (Piriz *et al.* 1997). It was shown that, using the isothermal approximation of the ablation front, the SBM reproduces the results of the self-consistent RTI theory (Sanz 1994, 1996). The dynamic overpressure stabilization can be explained as follows. When the ablation front moves to the corona (i.e. high-temperature region), its temperature does not increase because the corona region is isothermal; however, the temperature gradient increases. Then, the tops of the ripples, which come closer to the laser, receive more energy than the bottoms. The increased temperature gradient enhances local heat flux and then increases the rate of mass ablation, producing the dynamic pressure that pushes the front back. This results in an oscillation of the ablation surface ripple, and the peak value of the ripple amplitude is determined from the overpressure (Goncharov 1999; Goncharov *et al.* 2000).

Non-uniform laser irradiation also induces a deformation of the ablation surface and a rippled shock propagating ahead. This is called laser imprint, which also causes a RTI seed in ICF implosion (Emery *et al.* 1991; Ishizaki & Nishihara 1997, 1998; Taylor *et al.* 1997; Nishihara *et al.* 1998; Matsui *et al.* 1999; Goncharov *et al.* 2000; Velikovich *et al.* 2000). A rippled shock wave and ablation surface deformation caused by non-uniform laser irradiation on a smooth target were also evaluated with the deflagration model (Ishizaki & Nishihara 1997, 1998). In a weak shock limit, the approximate formulas for both the shock front ripple, *a*_{s}(*t*), and the asymptotic growth rate of the ablation surface, , are given by2.12and2.13where *β*_{s} is the shock Mach number with respect to the fluid behind the shock, and *K*_{2}=(1−*β*_{s})/(1+*β*_{s}). Figure 9 shows the normalized shock front ripple, *ka*_{s}/(*δI*/*I*), as a function of *r*_{s} and the normalized growth rate of the ablation surface, , as a function of *r*_{a}. The parameters used are the same as those of the previous problem. The shock front ripple increases with time at first and subsequently oscillates. The first maximum of the dimensionless shock ripple reaches approximately 0.65 for the parameters used. In this problem, since the non-uniformity is continuously supplied by laser, a finite asymptotic growth rate of the ablation front exists, in contrast to the previous case. However, since the growth rate is proportional to the inverse of the ablation velocity *v*_{a} as shown in equation (2.13), mass ablation also suppresses the forced ablative RMI. In a weak shock limit, the growth rate is proportional to the perturbation amplitude of the incident laser light as given in equation (2.13). However, the increase of the growth rate is relatively small for high-intensity laser light (Ishizaki & Nishihara 1998).

Laser smoothing techniques, such as induced spatial coherence (ISI) (Lehmberg *et al.* 1987) and smoothing by spectral dispersion (SSD; Skupsky *et al.* 1989), are implemented to reduce laser imprinting in ICF. Ishizaki & Nishihara (1998) also show a reduction in the growth rate by applying high-frequency modulation of the incident laser light. However, the perturbation growth of the ablation front has been observed to slow down instead of grow linearly even in cases of non-uniform lasers in experiments (Desselberger *et al.* 1995; Endo *et al.* 1995; Azechi *et al.* 1997). Moreover, simulations demonstrated that the mass perturbations exhibit decaying oscillations after reaching some maximum value (Nishihara *et al.* 1998). It should also be mentioned that the stabilization mechanism due to the mass ablation is considered only in equation (2.13), and that the overpressure effect leads to decaying oscillations of the ablation surface ripple rather than growth (Goncharov *et al.* 2000). The oscillations were also observed in experiments with a relatively long laser pulse (Aglitskiy *et al.* 2001, 2002).

Finite length effects of the thermal conduction region on the ablative RMI have been recently considered in detail (Abéguilé *et al.* 2006; Goncharov *et al.* 2006; Gotchev *et al.* 2006; Clarisse *et al.* 2008). Goncharov had pointed out that long-wavelength modes with *kD*_{c}<1 grow because of the combined effect of the RM-like and Landau–Darrieus instabilities (Landau & Lifshitz 1987), although the modes with *kD*_{c}>1 are stable and experience oscillatory behaviour by the dynamic overpressure, where *D*_{c} is the length of the thermal conduction region. Abéguilé *et al.* had obtained a self-similar model with continuous unsteady compressible ablation flows (Abéguilé *et al.* 2006; Clarisse *et al.* 2008). That model reveals the existence of three different regimes of *kD*_{c} for the evolution of the ablation surface deformation. The evolution is algebraic growth ∼*t*^{4/3} for *kD*_{c}≪1, modulated amplitude oscillation for 1<*kD*_{c}<10^{3} and damped oscillation for *kD*_{c}>10^{3}. Other examples of RM-like instabilities and many references can be found in Velikovich *et al.* (2000).

## 3. Nonlinear theory

When the shocks have travelled a distance greater than a wavelength, the system can be regarded as incompressible and irrotational (for weak shocks) except for the interface (contact surface) on which non-uniform vorticity is induced by shocks. Then we can treat the interface between two fluids as a non-uniform vortex sheet (Matsuoka *et al.* 2003; Matsuoka & Nishihara 2006*a*). In this section, we examine the nonlinear evolution of the interface analytically and numerically. We here restrict ourselves to nonlinear models for an initial sinusoidal perturbation without broad-band spectrum. It should be noted in the nonlinear analysis that the vorticity at the interface is assumed to be initially given in accordance with the asymptotic linear growth rates introduced in the previous section, , independently with the initial corrugation amplitude of the interface. The relation between the corrugation amplitude and the vorticity at the interface can be calculated using the linear theory in §2 for application to real systems. We also use time and space scales normalized by the linear growth rate as and by the initial wavelength for a planar case, *kx*, and the initial interface radius for a cylindrical case as *r*/*r*_{0}, which makes the results applicable to real systems in general. In §3*a*,*b*, we perform weakly nonlinear stability analysis and investigate the deviation from linear theory for both planar and cylindrical geometries. The fully (arbitrary amplitude) nonlinear evolution of the interface as a vortex sheet is presented in §3*c*,*d* by the vortex method (Krasny 1987).

### (a) Weakly nonlinear theory in planar geometry

We regard the interface as a curve in the (*x*,*y*) plane and parameterize it using a Lagrangian parameter *β*. The governing equations for the nonlinear analysis are the Bernoulli equation and kinematic boundary conditions. The Bernoulli equation, i.e. the pressure continuous condition at the interface, is given by3.1where *ρ*_{i} (*i*=1,2) is the density of fluid *i*, and the velocity potential *ϕ*_{i} is related to the fluid velocity *u*_{i} as *u*_{i}=∇*ϕ*_{i} in each region *i*. The interface *y*=*η*(*x*,*t*) satisfies the kinematic boundary condition, i.e. the normal velocity continuous condition at the interface,3.2where all quantities in equations (3.1) and (3.2) are estimated as deviations from *x*=0. Because the system is assumed to be incompressible, the velocity potential *ϕ*_{i} satisfies the Laplace equation △*ϕ*_{i}=0 in each fluid region *i* (*i*=1,2).

We normalize *β*→*kβ*, *x*→*kx*, *y*→*ky* and with the wavenumber *k* and the asymptotic linear growth rate in equation (2.6) in §2*c*. The velocity potential *ϕ*_{i} is also normalized by the asymptotic linear growth rate as .

With these dimensionless variables, we expand *ϕ*_{i} and *η* with a formal expansion parameter *ϵ* (*ϵ*≪1) asandwhere *ϕ*_{i} are expanded taking into account that they are solutions of the Laplace equation,3.3where Gauss’s symbol [(*m*+1)/2] denotes the maximum integer that does not exceed (*m*+1)/2, and is the amplitude of the (*m*−2*l*)th Fourier mode in the *m*th order of *ϵ* and the sign ∓ corresponds to the region *x*>0 or *x*<0. The circulation and vortex strength at the interface are given by *Γ*=*ϕ*_{1}−*ϕ*_{2} at *x*=0 and by *κ*=∇*Γ*, respectively, where ∇ should be taken along the interface. Therefore, the first order of the circulation at the interface, at *t*=0, corresponds to the initial circulation normalized by the wavenumber and the linear growth rate as . Note that the higher modes arise owing to the nonlinear mode coupling. The nonlinear evolution of RMI in incompressible fluids is thus determined from the initial circulation deposited at the interface. Expanding *ϕ*_{i} at *x*=0 asand substituting expansion (3.3) into equations (3.1) and (3.2), we obtain up to the third order of *O*(*ϵ*)3.4where *A*=(*ρ*_{2}−*ρ*_{1})/(*ρ*_{2}+*ρ*_{1}) is the Atwood number, where the fluid (*i*=2) is heavier than the fluid (*i*=1), and *a*_{0} is the initial amplitude of the interface normalized by the perturbation wavenumber. The first term of *η*^{(1)} is proportional to *t*, corresponding to linear growth. This result coincides with that obtained by Zhang & Sohn (1997), and a similar result is also obtained by Vandenboomgaerde *et al.* (2002) with different initial conditions. As mentioned before, the nonlinear theory can be developed for an arbitrary corrugation amplitude; however, it should be related to the initial circulation for a real system. The relation for real systems can be obtained from linear theory, namely according to an incident shock intensity and corrugation amplitude, the Atwood number and the material properties, as discussed in the previous section.

When we take expansion and contraction in the *x*-direction into account, we can obtain a perturbation solution in better agreement with the hydrodynamic simulation result (Matsuoka *et al.* 2003) using IMPACT2D (Sakagami & Nishihara 1990*b*) for a longer time than solution (3.4). In this case, we set *y*=*β*+*ξ*(*β*,*t*) and *x*=*η*(*β*,*t*), and substitute them into equations (3.1) and (3.2). Then we have a perturbation solution up to the third order of *O*(*ϵ*) (Matsuoka *et al.* 2003),3.5

Figure 10 compares the time-dependent weakly nonlinear growth rates with a hydrodynamic simulation. The time is normalized by the asymptotic linear growth rate . The parameters used are given in Matsuoka *et al.* (2003). Positive and negative values correspond to the spike and bubble growth rates, respectively, normalized by the linear growth rate. The analytical growth rates of the spike and bubble can be obtained from at *β*=0 and *β*=*π*/2, respectively. Solid lines with open and closed circles indicate the growth rate for the shocked interface in the hydrodynamic simulation and that for the vortex sheet, respectively. They agree fairly well except near the beginning, owing to the finite transition time required for the incident shock to pass through the interface with a finite corrugation amplitude and also that required for the instability to reach the asymptotic linear growth rate. The oscillation of the growth rate for the shocked case occurs because the sound waves propagate back and forth between the reflected or transmitted shocks and the interface. In the linear regime, the growth rates of the bubble and spike are the same. However, as shown in figure 10, they are not the same in the nonlinear regime, where the growth rate of the spike is larger than that of the bubble. The nonlinear growth of the spike increases from its initial value and reaches its maximum value around a time of in the case of the vortex sheet and in the case of the shocked interface, and then decays slowly compared with that of the bubble. As we see from figure 10, the nonlinear growth rate of the bubble decays faster than that of the spike in both cases.

The above-mentioned agreement between simulation and analytical solution indicates that RMI grows owing to the non-uniform velocity shear that the transmitted and reflected shocks leave at the corrugated interface. We would like to mention that the use of the linear asymptotic growth rate for the initial velocity shear gives not only qualitative but also quantitative agreement between numerical simulations and analytical solutions.

The difference between the simulation results and the analytical results can be improved by using the Padé approximation (Baker & Graves-Morris 1981; Zhang & Sohn 1997; Matsuoka *et al.* 2003). The result using the Padé approximants (Matsuoka *et al.* 2003) is also shown in figure 10. We can see that the analytical results of the nonlinear growth rates agree quite well with the simulation results for a relatively long time. All of the previous works that used the impulsive model had to introduce a free parameter to adjust the nonlinear growth rate. However, our model can uniquely determine the nonlinear evolution of RMI from the initial corrugation amplitude, incident shock intensity and Atwood number before the incident shock hits the interface.

We now compare our nonlinear prediction with experiments by Dimonte *et al.* (1996), in which a rarefaction wave was reflected and phase inversion of an interface occurred. Taking this into account, we evaluate the post-shock variables such as the Atwood number *A*=0.58426, the corrugation amplitude *a*_{0}=−2.0482/100 and the linear growth rate cm s^{−1}. In figure 11, we plot the bubble and spike amplitude as a function of time and compare them with the experimental points. Note that, in the experiments, they did not resolve the bubble and spike separately but showed only Fourier mode amplitudes. Therefore, we also plot the amplitudes taking an average between the bubble and the spike. Despite the fact that the incident shock Mach number in the experiments is very high, *M*=15.3, the analytical solutions agree fairly well with the experiment.

### (b) Weakly nonlinear theory in cylindrical geometry

In this section, we consider the nonlinear stability of an RMI interface for incompressible fluids with cylindrical geometry. The assumption of incompressible fluids may not be realistic for cylindrical geometry because it cannot provide shock bouncing at the centre. However, it still represents important features of cylindrical geometry, such as the existence of two independent spatial scales, radius and wavelength, and also ingoing and outgoing growth in addition to bubbles and spikes. Throughout this section and §3*d*, we determine the subscripts *i*=1 and 2 correspond to the inner and outer fluids; therefore, *A*>0 and *A*<0 correspond to a case of light inner fluid and heavy outer fluid, and vice versa.

The governing equations for stability analysis are equations (3.1) and (3.2). However, they must be rewritten in the circular coordinate system, i.e. (*x*,*y*)→(*r*,*θ*), and a cylindrical interface is given by *r*=*η*(*θ*,*t*). In a weakly nonlinear theory in cylindrical geometry, all quantities should be estimated at the initial interface *r*=*r*_{0}, where *r*_{0} is the initial radius of the interface.

We replace *η*, *ϕ*_{i}, *r* and time *t* with dimensionless variables *η*→*η*/*r*_{0}, *ϕ*_{i}→*ϕ*_{i}/(*r*_{0}*v*_{0}), *r*→*r*/*r*_{0} and *t*→*tv*_{0}/*r*_{0}, where *v*_{0} is the initial growth rate, which corresponds to the linear growth rate of a planar RMI for the wavelength *λ*=2*πr*_{0}/*n* with mode number *n* (Matsuoka & Nishihara 2006*b*). With these dimensionless variables, we expand *ϕ*_{i} (*i*=1,2) and *η* with a formal expansion parameter *ϵ* (*ϵ*≪1) as3.6where the sign, ±, corresponds to the region (*r*<*r*_{0}) or (*r*>*r*_{0}).

Expanding *ϕ*_{i} at the normalized initial radius *r*=1, we obtain a perturbation solution up to the third order of *O*(*ϵ*) (Matsuoka & Nishihara 2006*c*)3.7where *a*_{0} is the normalized initial amplitude, *a*_{0}/*r*_{0}, and we impose here the initial perturbation of the interface and its normalized normal perturbation velocity, . The corresponding initial circulation is given as at *t*=0; it corresponds to the initial circulation normalized by the initial interface radius and the linear growth rate as *Γ*/(*r*_{0}*v*_{0}).

Figure 12 shows analytical growth rates for various Atwood numbers for modes *n*=1,2,3 (Matsuoka & Nishihara 2006*c*). Here, we set *a*_{0}=0, i.e. the initial interface is a circle with radius *r*=1, but the initial circulation is deposited along a uniform interface. The weakly nonlinear growth rates in the cylindrical case are quite different from those in the planar case. As shown in figure 12, the weakly nonlinear growth rates strongly depend on the mode number corresponding to two independent spatial parameters, radius and wavelength, and also on ingoing or outgoing growth. In cylindrical geometry, the nonlinear growth rate of the ingoing spikes increases with time, especially for cases of light inner fluid, *A*>0. The increase is more apparent with increasing density difference and for lower modes. Conversely, the nonlinear growth rate of the outgoing bubbles decreases with time; however, it is hardly affected by the Atwood number for *A*>0. For the cases of heavier inner fluid, *A*<0, the weakly nonlinear growth of spikes depends weakly on the density difference compared with that for *A*>0. The weakly nonlinear growth rates of the ingoing bubbles become either larger or smaller than the linear growth rates, depending on the density difference and mode number. When *A*=−1.0 and *n*=1 (figure 12*d*), the bubble and spike growth rates remain approximately at the value of the linear growth rate for a long time. These results clearly show that weakly nonlinear growth in cylindrical geometry is strongly affected by the direction of growth, i.e. ingoing or outgoing, in addition to the bubble and spike. The spikes penetrate into light fluid much faster than the linear growth rate especially for low mode even in RMI. This acceleration of spike penetration is also observed in two-dimensional and three-dimensional simulations of RTI (Sakagami & Nishihara 1990*a*,*b*).

### (c) Arbitrary amplitude theory in planar geometry

In this section, we numerically investigate the fully nonlinear behaviour of a planar interface in RMI with the use of the Birkhoff–Rott equation (Rott 1956; Birkhoff 1962; Saffman 1992). This equation and the vortex method (Krasny 1987; Kotelnikov *et al.* 2000; Sohn 2004; Matsuoka & Nishihara 2006*a*) enable us to calculate the roll-up of the interface. The velocity of the interface (*x*,*y*)=(*X*(*β*,*t*),*Y* (*β*,*t*)) is given as (Baker *et al.* 1982; Matsuoka & Nishihara 2006*a*)3.8where *α* is an artificial parameter such that *α*≠0 when the Atwood number *A*≠0 (Baker *et al.* 1982; Pullin 1982), the vortex sheet strength *κ* is defined by *κ*=∂*Γ*/∂*s*=*Γ*_{β}/*s*_{β}, *Γ*=*ϕ*_{1}−*ϕ*_{2}, and *s* is the arc length of the sheet. The subscript denotes differentiation with respect to the variable. The vortex-induced velocity *U*=*U*(*β*,*t*) and *V* =*V* (*β*,*t*) are given by the Birkhoff–Rott equation,3.9where we regularize the Cauchy integral using Krasny’s *δ* (Krasny 1987). The integration represents the non-local character of the interface dynamics. Equation (3.8) and the evolution equation for *κ*3.10derived from the Bernoulli equation (Matsuoka & Nishihara 2006*a*) are the governing equations of our numerical calculations. Details of the numerical methods for solving equations (3.8) and (3.10) are presented in Matsuoka & Nishihara (2006*a*). The initial condition for numerical calculations is given aswhere all variables are normalized by the wave number *k* and the linear growth rate in a manner similar to that in §3*a* and the above-mentioned equations represent the same initial conditions as given in that section.

Figure 13 shows the temporal evolution of an interface for Atwood number *A*=0.155, and initial amplitude *a*_{0}=0.2 over the normalized time 0≤*t*≤9.4. This corrugation amplitude, the Atwood number, the time and the regularized parameter *δ* (*δ*=0.15 here) are chosen so that the numerical computations are directly comparable with the experiment by Jacobs & Sheeley (1996). *a*_{0}=0.2 corresponds to the interface perturbation amplitude when a tank was detached from a coil spring in the experiment 0.3 cm in a real scale. We can also estimate *k*=0.81 cm^{−1} and s^{−1} from the linear theory described in §2*e*. With these values, the normalized time *t*=1 corresponds to a real time of 63.3 ms. Profiles (*a*),(*b*),(*c*),(*d*),(*e*) and (*f*) in figure 13 correspond to video images (*c*),(*e*),(*g*),(*h*),(*i*) and (*j*), respectively, in fig. 4 in Jacobs & Sheeley (1996). When we set the regularized parameter *δ*=0 in equation (3.9), interface roll-up does not occur and we need highly accurate numerical schemes. For details of this case, see Matsuoka & Nishihara (2006*a*).

An example of a multi-mode profile in the RMI with initial conditionsis shown in figure 14, where the calculation breaks down after a few time steps of *t*=1.0. Similar profiles have been obtained in experiments by Jacobs & Sheeley (1996). Various multi-mode profiles including the RTI are presented in Matsuoka & Nishihara (2008). The fully nonlinear theory can track real experiments.

Nonlinear coherent structure of bubbles and spikes in RMI has been investigated in detail numerically and analytically (Abarzhi 2002; Abarzhi *et al.* 2003; Herrmann *et al.* 2008), which shows the flattening of the bubble fronts as a consequence of multi-scale character of the coherent dynamics.

### (d) Arbitrary amplitude theory in cylindrical geometry

In this section we investigate the fully nonlinear motion of a circular interface using the same method as in §3*c*. The interface velocity in (*r*,*θ*) coordinates is given as (Matsuoka & Nishihara 2006*b*)3.11where *α*=*α*(*A*) is an artificial parameter, *β* is a Lagrangian parameter as in the previous sections, and the (regularized) vortex-induced velocity *q*^{r} and *q*^{θ} are given by the Biot–Savart Law3.12

The evolution equation for the sheet strength *κ* in cylindrical coordinates becomes3.13

By solving equations (3.11) and (3.13) numerically, we can determine the nonlinear motion of a vortex sheet in RMI in cylindrical geometry.

The normalized initial conditions are given asfor the mode number *n* (for details, see Matsuoka & Nishihara (2006*b*)).

Figure 15 shows the temporal evolution of the interfacial profiles for various modes *n*=1,2,3 and 8, and Atwood numbers *A*=±0.2. The number of grid points *N* in the computations is taken as *N*=256*n* for the mode number *n*. Here we take the regularized parameter *δ*=0 and *δ*=0.1 in the calculations for mode *n*=1 and other modes, respectively. For *A*>0 (*A*<0), the outgoing parts *β*=(2 *mπ*)/*n*, and the ingoing parts *β*=(1+2*m*)*π*/*n*,(*m*=0,…,*n*−1) correspond to bubbles (spikes) and spikes (bubbles). The interface does not cross the origin for any of the modes in an incompressible cylindrical fluid, since the velocity potential is given as (*r*<*r*_{0}), (*r*>*r*_{0}), and then the potentials diverge when the interface approaches the origin.

In a planar case, the profiles of bubbles and spikes are very different from each other. However, as shown in figure 15, the profile of the outgoing bubbles (ingoing spikes) for *A*>0 is not so different from those of the outgoing spikes (ingoing bubbles) for *A*<0. The interface profile in cylindrical geometry is, therefore, determined mainly from either outgoing or ingoing growth, regardless of bubbles and spikes, for the Atwood number chosen, *A*=±0.2. For a large absolute Atwood number, the nonlinear growth may also depend on either spikes or bubbles, as discussed in §3*b*. The interface profiles at *t*=2 in figure 15 are also quite similar for *n*≥3. These results indicate that the fully nonlinear evolution of growth in cylindrical geometry strongly depends on the mode number and is scaled by the initial growth rate *v*_{0} and radius, instead of by the perturbation wavelength at least for small mode numbers.

## 4. Molecular dynamics simulations of RMI

The MD approach is based on atom motion tracking and can hence provide all information regarding system behaviour, including phase transitions in matter, diffusion, thermoconductivity and other non-equilibrium processes far from the thermodynamic steady state. Because it directly computes atom–atom interactions, the MD method has fundamental advantages over numerical hydrodynamic methods, which assume shocks as a structureless discontinuity and requires an EOS. Different hydrodynamic codes sometimes show significant deviations in the simulation flow, and a higher mesh resolution cannot always improve the situation (Kamm 1999). Later we discuss recent MD simulations of RMI induced in shocked material in planar and cylindrical geometries.

In several works (Zhakhovskii *et al.* 2006; Zybin *et al.* 2006), MD simulations of RMI in solids were performed in planar geometry. For the solid–solid interface, the system was composed of inert gas crystals with different densities, light *ρ*_{1} and heavy *ρ*_{2}, with a density ratio of *ρ*_{2}/*ρ*_{1}=4. Two cases of shock propagation were considered: from the light to the heavy solid and from the heavy to the light solid. In addition, the simulation of a planar shock refracted through the perturbed interface between a copper crystal and a vacuum was done. Such a simulation corresponds to experiments in which a metal with corrugations at the free rear-side surface is shocked from the front side and supersonic jet formation is observed on the rear side.

For the solid–solid case, the atomic interaction was modelled via a modified Lennard-Jones (LJ) pair potential (Zhakhovskii *et al.* 2002), and for the solid-vacuum case, the Mishin’s EAM potential for copper was used (Mishin *et al.* 2003). The LJ system was built from two LJ solids that differ only in atomic mass: light *m*_{1}=*m*_{a} and heavy *m*_{2}=4*m*_{a}. About 20 million atoms were used to build the solid–solid system with dimensions 108×15×487 nm. Hereinafter we use MD units (mdu) based on the common dimensionless form of the LJ potential, where *σ* is a unit of length, *ϵ* is an energy unit and the unit of time is The copper crystal had about 6 million atoms and a size of 74×5.1×218 nm. Periodic boundaries were imposed in both transverse directions.

Zybin *et al.* (2006) showed that the growth rate of RMI strongly depends on the initial shock strength in a light-to-heavy solid–solid system. For instance, almost no growth occurs at the piston velocity *u*_{p}=0.7 mdu; however, at a higher piston velocity *u*_{p}=1 mdu, growth starts as soon as the shock passes through the interface. Such threshold behaviour in solids is essentially different from that for RMI in liquids and indicates strong dependence of the growth conditions on the material properties (e.g. shear strength), generation of dislocations, defects and plastic flow. Local melting may occur near the interface even though a chosen piston velocity does not reach some critical value at which the light solid could start to melt.

The evolution of the interface for a shock wave moving from a heavy to a light solid was simulated as well in Zybin *et al.* (2006). In this case, instead of a reflected shock, a non-planar rarefaction wave moves backwards into the first (heavy) material as shown in figure 16. The transmitted shock wave *u*_{s1} also becomes non-planar with front oscillations, but the oscillations soon decay. After the shock passage, the perturbation changes its phase, and the spike starts to grow. In contrast to the light-to-heavy system, in this case the growth rate of the spike amplitude is slowing down to zero when the rarefaction wave from the free right surface gets back to the piston. Multiple intersections of the transmitted non-planar shock waves (e.g. as seen in the second snapshot in figure 16) create local areas with increased stress and temperature around the growing spike (‘hot spots’) where the material might become either molten or amorphous.

In modelling the shock-induced instability of the perturbed copper rear surface (Zybin *et al.* 2006), a shock wave was generated along the [110] direction. A characteristic perturbation has a sinusoidal-like shape with wavelength *λ*=74 nm and the amplitude *ψ*_{0}=*λ*/4. Figure 17 shows two different examples of instability growth when the shock wave reaches the free surface and the rarefaction wave moves backwards into the bulk. For a large corrugation amplitude of *ψ*_{0}=*λ* (see figure 17*a*), the spike starts to grow for a piston velocity *u*_{p}>1 km s^{−1}. Remarkably, the cup of the spike displays a high density of defects compared with that of dislocations in the bulk crystal, indicating possible amorphization of the material inside of and beneath the spike, which can facilitate plastic flow around it and sustain its growth. In the case of a smaller perturbation *ψ*_{0}=*λ*/4, spike growth was observed for piston velocities above 2.0 km s^{−1}. The growth is associated with the appearance of local spots that might be nuclei of molten material, as shown in figure 17*b*. In addition, void formation and nucleation with incipient spall have been detected underneath the hills of the initial perturbation during its phase inversion and spike growth at the unloading surface. Void formation is not observed in the simulations of a shock reflected from a plane unperturbed surface; the authors related it to the intersections of non-planar rarefaction waves moving down from the perturbed interface (Zybin *et al.* 2006).

An MD study of RMI in a cylindrical target having an interface between two LJ liquids was presented in Zhakhovskii *et al.* (2002, 2006). A converging shock was induced by a cylindrical piston surrounding a target. At *t*=0, the piston starts to move with constant velocity *u*_{p}, and, at *t*=16 mdu, it is removed from the system. A shock Mach number of *M*=3.3 was measured for piston velocity *u*_{p}=1 mdu. Figure 18 shows snapshots of the density and vorticity fields for a shock propagating from heavy to light liquids with mass ratio *m*_{h}/*m*_{l}=9 and *u*_{p}=0.5 mdu. Maps of the flow velocity field for a simulation with *u*_{p}=1 mdu are presented in figure 19. When the shock front hits the interface, a transmitted shock and rarefaction wave are produced (figures 18 and 19*a*). The vorticities are clearly allocated around the material interface, causing a spike to form (figure 19*b*). The shock reflected from the centre of the cylinder hits the interface repeatedly and causes significant growth of the instability, which is quite different from the planar case. Figure 18*a*(iv),*b*(iv) shows the RMI late in the simulation, when the total width of the mixing zone goes into the asymptotic stage of growth.

The nonlinear evolution of the growth rates for different mode numbers was estimated from differentiation of the asymptotic value of the thickness of the mixing zone *Δ*∼*Δ*(*t*) obtained by fitting MD-simulated data. In the case of mode 8, the perturbations grow asymptotically as *Δ*_{8}∼*t*^{0.45}, and for modes 5 and 3 they grow as *Δ*_{5}∼*t*^{0.3} and , respectively. Hence the growth rates decay as *dΔ*_{8}/*dt*∼*t*^{−0.55}, *dΔ*_{5}/*dt*∼*t*^{−0.7} and *dΔ*_{3}/*dt*∼*t*^{−1}, respectively.

## 5. Conclusion

In this article, we re-examined RMI and RM-like instabilities in terms of vorticity. We may conclude that RMI and RM-like instabilities extend well beyond the strict limits of applicability in Richtmyer’s original analogy of the perturbed interface with a pendulum driven by an impulsive acceleration. The perturbation growth at an interface is driven by vorticity, which is either deposited initially or supplied by external sources. The vorticity is localized near the interface and distributed non-uniformly along the interface.

In the linear theory of shock–interface RMI, we described the deposition of vorticity through the interaction between shock fronts and an interface by sound waves, and gave an exact formula for the asymptotic growth rate. We also discussed RM-like instabilities including ablative RMI, which was stabilized owing to mass ablation and dynamical overpressure.

In the limit of incompressible approximation, we showed the weakly and fully nonlinear dynamics of an interface, along which an initial sinusoidal vorticity was deposited, for planar and cylindrical geometries. Existence of two independent spatial scales in the cylindrical geometry results in quite different nonlinear evolution for the planar case.

MD simulation reveals nonlinear RMI dynamics in solids and liquids with RMI-induced phase transition, although further study is required.

Although we limited our discussion in two-dimensional interfaces, there are theoretical works for three-dimensional interfaces using group theory (Abarzhi 2002, 2008). It should also be mentioned that it is possible to extend the theory in §3*c*,*d* for a case with three-dimensional interface and two-dimensional sheet strength. Although transition to a turbulent regime is not well understood yet, there has been an attempt to solve shock interaction with an isotropic turbulent vorticity field (Wouchuk *et al.* 2009).

## Acknowledgements

We express our sincere thanks to Dr S. Abarzhi for her encouragement to write this review. Part of this work is supported by a Grant-in-Aid for Scientific Research (C and B) from the Japan Society for the Promotion of Science.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond’.

- © 2010 The Royal Society