In recent years, with the advent of high-power lasers, laser-driven shocks have become widely used to determine equations of state (EOS) for various materials^1–13 of fundamental interest to inertial confinement fusion (ICF).^14–19 According to the Rankine–Hugoniot relations, the determination of all the flow variables of a shock front in such experiments can be reduced to the determination of two kinematic parameters, namely, the propagation speed v_s of the shock wave and the average particle velocity u behind the shock front.²⁰ Although v_s can be measured with relative ease, determining u is more difficult experimentally. One commonly used indirect method is the so-called free-surface approximation,²⁰ which assumes that the velocity of the free surface is twice the average particle velocity, i.e., v_f = 2u, which is also known as the doubling rule. It has been shown both theoretically²⁰ and experimentally^21–25 that this approximation works well at relatively low Mach number. However, the way in which v_f departs from this rule under strong shocks remains an open question. Better understanding of this question is of great interest to laser-driven EOS experiments, in which ultra-strong shocks are often generated.

The most frequently used technique for measuring the velocities of free surfaces and interfaces is the velocity interferometer system for any reflector (VISAR), which has been improved^26–30 and applied in many experiments.^31–35 An important mechanism underlying the use of this system is the reflection of the detecting laser from a moving surface. However, since a perfect two-dimensional surface does not exist in reality, the positions where the light is reflected and what the device really measures are still a matter of debate.³⁶ A clear understanding of the detailed structure and dynamical evolution of particles in the vicinity of the shock front can help provide answers to these questions.

Our study focuses on the free-surface velocity v_f under ultra-strong shocks. With nonequilibrium molecular dynamics (NEMD) simulation methods, we find that the free surface no longer exists under ultra-strong shocks. It is replaced by a transient regime of expanding plasma, and v_f is actually the velocity of the critical surface of the reflected detecting laser light. The evolution of v_f exhibits a single-peaked structure as time passes by. v_f only approximately agrees with the doubling rule in a short transient period immediately after the shock breaks out. This is then followed by a short acceleration process to reach its maximum. After that, v_f slowly decays with the expansion of the material. During the decay process, v_f gradually converges to a value that depends on the wavelength of the detecting laser, but is independent of the atomic mass (or density) of the medium.

Simulation results further show that this trend of v_f is the result of two competing mechanisms. One is the hydrodynamic acceleration provided by the pressure gradient. The other is the backward shift of the reflecting surface caused by the decreasing density of particles during the expansion of the transient plasma surface. It is also found that the dynamical behavior of v_f can be described by a scaling law, which may be used to improve the measurement of u in laser-driven EOS experiments.

The rest of this article is organized as follows. Details of the NEMD simulations are presented in Sec. II. The main results and corresponding analyses are presented in Sec. III. A short summary is given in Sec. IV.

In this work, the NEMD method³⁷ is employed to study the structure and dynamical process near the surface of a medium surface under ultra-strong shocks. This method removes the limitation imposed by the assumption of local thermal equilibrium (LTE),³⁸ making it suitable for studies of ultra-strong shock compression in dense materials. In addition, molecular dynamics (MD) simulations can directly provide microscopic information about the structure of shock waves, which may not be accessed by hydrodynamical calculations or other methods.

The dimensions of the simulation box are 50 × 50 × 3000 ų. As displayed schematically in Fig. 1(a), 62 800 X atoms are assembled in the left 400 Å of this box. A moving reflective boundary is used as a piston to generate a shock wave propagating along the long side of the box, denoted as the z axis. The piston is placed at z = 0 at t = 0 and then moves at a constant velocity v_p = 75 km/s, 50 km/s, or 25 km/s. Along the two perpendicular axes, denoted as the x and y axes, periodic boundaries are used to prevent finite-size effects of boundary reflections. In addition, the simulation box is divided into slices of 2 Å thickness along the z axis to calculate the statistics of macroscopic observables and thermodynamic variables. Details of the simulation can also be found in our previous work.³⁸ To reveal the effect of relative atomic mass m_X (also the effect of density) on v_f, we carry out the simulation with several isotopes of X, such as ³X, ⁴X, ⁶X, ⁸X, and ¹⁰X, but keep the initial number density unchanged.

Before the arrival of the shock front, the interfacial regime between the medium and vacuum can be viewed as a perfect two-dimensional surface. However, after the breakout of the shock, the free surface is replaced by a regime of expanding plasma, and the position depends on how one measures it. Experimental devices, such as VISAR, rely on the Doppler effect to obtain the velocity of the interface.²⁶ The position at which the reflection of the detecting laser takes place is essential to the measurements. In plasma environments, it is the critical number density of electrons n_ec that determines the position at which the detecting laser is reflected.⁴² The critical number density can be calculated as⁴³n_ec = (2πc/λ)²m_eϵ₀/e², where λ is the wavelength of the laser, c is the speed of light, e and m_e are the electron charge and mass, and ϵ₀ is the vacuum permittivity.

Here, we may assume the critical number density of atoms to be n_c = n_ec/2 for each X atom with two electrons, considering all electrons to be ionized under an ultra-strong shock as an illustrative simplification. Although the specific value of n_c may deviate from this estimate owing to different degrees of ionization and separation between free electrons and ions,⁴⁴ the behavior of v_f will be quite similar so long as the value of n_c is of the same order of magnitude, as we shall show in Sec. III. After n_c has been determined, the slices along the z axis with their number density n ∈ [0.75n_c, 1.25n_c] are considered as the “reflecting surface” with a finite thickness. A more accurate estimate of the value of n_c is beyond the scope of the present work. Such an estimate might be obtained with simulation techniques including a faithful description of excited electrons as well as the coupling effect between the radiation field and the hot dense plasma.

We first present a microscopic picture of what happens after an ultra-strong shock arrives at the surface. With the NEMD simulations, one can have atomic resolution for the fine structure of shock fronts. Figure 1(a) displays the entire process of an ultra-strong shock meeting with and being bounced back by the initial free surface. Before the shock front, indicated by the red vertical dashed lines in the figure, touches the initial free surface, the shock wave propagates inside the medium with its microscopic structure depicted well in previous studies.^38,44 The velocity of the shock wave is about v_s = 68.9 km/s for m_X = 4, corresponding to a piston speed v_p = 50.0 km/s, a compression ratio η = 3.64, and a thickness of the shock front δ of about 24 Å. Note that the thickness of the shock front varies depending on its definition. Here, we follow the definition used in our previous work.³⁸ After the shock front meets the initial free surface at z = 400 Å, a rarefaction wave is formed as the result of surface reflection. The formation of this rarefaction wave is illustrated in Fig. 1(b), where we plot the distributions of density and velocity along the z axis right after the shock breaks out. In the rarefaction wave regime (the rectangular area in the figure), the velocity is lower than that in the bulk medium, and the density is gradually increasing. These are typical features of a rarefaction wave.

The appearance of a rarefaction wave is expected as a direct consequence of continuum hydrodynamic theories. However, continuum theories do not provide any information on how the materials ahead of the initial free surface move, and such information is critical for the measurement of v_f. It can be seen from Fig. 1(a) that the transient regime between bulk medium and vacuum starts to move forward along the z axis after t = 0.6 ps, but the boundary is vague. From the distribution of velocity in Fig. 1(b), one can see that the velocity of the particles outside the free surface increases when the position is away from the surface. Meanwhile, the density decreases. At a distance far from the surface, the distribution of velocity consists of a number of spikes, which implies that a continuum description of the material is no longer valid.

Nevertheless, the distributions of velocity and density provide a physical picture of how material leaves the initial free surface. Shortly after the shock front meets the initial free surface, the particles outside can be considered as leaving the surface at the same moment at which the shock breaks out. The density distribution f(z) in Fig. 1(b) actually reflects the number distribution f(v_z) in velocity space at the very moment when the shock front meets the initial free surface. It follows an exponentially decay law similar to a Maxwellian distribution: $f(z)=p_1\exp [-p_2{(z-z_0)}^2]$, with fitting parameters z₀ = 383.6 Å, p₁ = 2.256 × 10²³ cm⁻³, and p₂ = 3.404 × 10⁻³ Å⁻². Thus, the process of material leaving the surface into vacuum is quite similar to the process of evaporation, as long as one observes it in the reference frame moving with the bulk material at velocity u.

With this physical picture in mind, one can now consider how to determine v_f. Since both density and velocity vary sharply with increasing distance from the initial free surface, as displayed in Fig. 1(b), v_f depends strongly on the position of the reflecting surface. Here, we take n_c = 2.348 × 10²¹ cm⁻³, which corresponds to λ = 5000 Å. Note that this n_c is the number density of X ions. The density of electrons is twice this quantity, assuming that both of the electrons are ionized. The evolution of v_f with respect to time is displayed in Figs. 2(a)–2(c) for three typical cases with v_p = 75 km/s, 50 km/s, and 25 km/s, respectively. It shows that the measured surface speed v_f agrees with the doubling rule only in a short transient period immediately after the shock front meets the initial free surface. For example, in Fig. 2(b), the average velocity along the z axis is 100 km/s at t = 0.6 ps, which is twice the velocity of particles behind the shock wave, namely, 50 km/s. However, the doubling rule remains valid only at times less than 200 fs, and therefore it cannot be picked up by experimental devices with a temporal resolution of several picoseconds.²⁹ After that, v_f increases further to reach a peak value and then slowly decays.

To find out how competition between kinetic and potential energies affects the evolution of v_f, we also vary the atomic mass of the X ions from 3 to 10 in the simulations. In this way, we can keep v_s roughly the same at a given piston speed v_p,^20,38 while changing the ratio between kinetic and potential energies. For example, for the case in which the piston is pushed forward at v_p = 50 km/s, the measured η = 3.64, and no apparent variation is observed. The measured v_s varies from 66.7 km/s (in the case of ¹⁰X) to 69.4 km/s (in the case of ³X), which fluctuate around the value v_s = 68.9 km/s calculated under the continuum condition. Figure 2 reveals that the potential energy affects the evolution of v_f only in the acceleration phase, which is a short period less than 2 ps. After that, the v_f values for different atomic masses converge together, displaying some features of scaling invariance.

Another important factor that may affect the measurement of v_f is the wavelength λ of the detecting laser. Since n_c is proportional to 1/λ², detecting lasers of different wavelengths will be reflected on the surface at points corresponding to their particular value of n_c. Thus, v_f of that surface depends on λ. The velocities of the reflecting surfaces for selected values of n_c are plotted in Fig. 3. It can be seen that v_f for each n_c has a single-peaked structure similar to those in Fig. 2. In addition, v_f will reach a converged value independent of the atomic mass of X (but depending on the detecting wavelength) after a period of acceleration, regardless of the wavelength of the detecting laser.

The trend that v_f converges to a value independent of the atomic mass of X can be attributed to two competing mechanisms. The acceleration mechanism comes mainly from the pressure gradient. The distributions of pressure along the z axis are shown in Fig. 4 for the case v_p = 50 km/s and a detecting wavelength of 5000 Å. For the purpose of comparison, the pressure contributions for ⁴X ions are plotted as solid curves, while those for ⁸X ions are plotted as dashed curves. Figure 4(a) shows the distribution of pressure contributions at t = 0.6 ps, at which v_f is in the acceleration phase, and Fig. 4(b) shows the distribution at a time (t = 2.0 ps) when v_f is not affected by the mass of X. The pressure contribution originating from kinetic energy is called the thermal pressure and denoted by p_th. It is proportional to the average kinetic energy E_k of particles as $p_{\text{th}}=\frac{2}{3}n{E}_k$.³⁸ The contribution from interaction is called the virial pressure and denoted by p_virial. The region between each pair of dashed vertical lines of the same color is the thin layer of the “reflecting surface” in which the number density of X is between 0.75n_c and 1.25n_c, as described in Sec. II. Note that this region is not the same for different atomic masses: the regions for ⁴X and ⁸X ions are distinguished by different colors. It can be seen from Fig. 4 that there is a pressure gradient across the thin layer, which pushes the layer forward with a certain acceleration. From the insets in Figs. 4(a) and 4(b), it can also be seen that this acceleration decreases with time, as the result of a decreasing pressure gradient.

To quantitatively compare the influence of the pressure gradient on the acceleration, all pressure contributions presented in Fig. 4 are divided by the relative atomic mass of the ions. It is found that the acceleration provided by thermal pressure is nearly the same, since the blue dashed curves are close to the solid black curves in Fig. 4. This is easy to understand, because p_th is proportional to atomic mass by definition. However, the effect of the virial pressure differs greatly among the different cases. In Fig. 4(a), where v_f is in the acceleration phase, p_virial is of the same order of magnitude as p_th in the region of bulk material far behind the reflecting surface. In the vicinity of the reflecting surface, p_virial of ⁴X is even larger than p_th, while p_virial of ⁸X is less than 1% of p_th. On the other hand, p_virial is less than 10% of p_th in the region containing bulk material in Fig. 4(b), where v_f is independent of the atomic mass of the ions. In the inset in Fig. 4(b), the contribution from interaction, p_virial, is almost unobservable near the reflecting surface. Therefore, one can conclude that the pressure gradient resulting from interaction, i.e., from p_virial, is the origin of the different accelerations of v_f in the acceleration phase shown in Fig. 2.

The average particle flow in the vicinity of the reflecting surface provides a deceleration mechanism. The physical picture is as follows. When the reflecting surface, which is actually a thin layer microscopically, moves forward, some particles leave the layer and enter the region ahead. The majority of these particles have a speed higher than v_f. Meanwhile, more particles join the layer from behind, most of which have speeds lower than v_f. The decelerating effect on v_f comes mainly from those particles joining the layer with lower speed and those leaving with higher speed. The overall effect is that the particles in the layer are gradually replaced by the particles behind the reflecting surface, and this eventually leads to a decrease in v_f.

The distributions of number density and velocity are shown in Fig. 5, where we again take ⁴X and ⁸X as examples. At t = 0.8 ps, both the local average velocity v_z along the z axis and the number density n differ between the two atomic masses, as can be seen from Figs. 5(a) and 5(b). However, at t = 2.0 ps, the distributions of v_z and n are almost the same for the two atomic masses. This suggests that the system evolves to a state in which the dynamics along the z axis is insensitive to atomic mass. This results from expansion of high-temperature plasma in the transient regime, which makes the potential energy much smaller than the kinetic energy. The dynamical behavior of the system is thus dominated by the kinetic energy, which accounts for the insensitivity of v_z and n to m_X in Figs. 5(c) and 5(d).

Figure 6 shows ${\tilde{v}}_f/\alpha$ obtained in all of our simulations as a function of $\tilde{t}$, with atomic mass varying from 3 to 10 and piston speed from 25 km/s to 75 km/s. The dimensionless function ${\tilde{v}}_f/\alpha$ with the best-fitting parameters is also displayed as a dashed curve. It can be seen that under the scaling transformations v_f → v_f/(αc_s) and t − t₀ → (t − t₀)c_s/l₀, v_f can be described well by a unified function of scaled time. The figure also reveals that v_f may be 30%–40% higher than the estimate from the doubling rule. This value may not be applied directly to improve experimental measurements, since it is obtained from simulations that do not take the motion of electrons into consideration. However, it does show that overestimation of particle velocity u in experimental measurements via the doubling rule should be considered, which leads further to overestimation of η. To use the scaling law to correct the overestimation in practice, one may first determine the forms of α(M) and $f(\tilde{t})$ in Eqs. (6) and (7) by comparing the experimental results and theoretical calculations for selected values of η and then obtain corrections for all experimental points.

In summary, we have used NEMD simulations to provide an atomic-scale picture of the dynamics of particles near the surface of a medium under ultra-strong shocks. This has revealed that v_f under ultra-strong shocks has a single-peaked structure. After a short period of acceleration, v_f exhibits a long slowly decaying tail, which is not sensitive to the atomic mass of the medium. A scaling law has also been revealed to describe the dynamics of v_f, and this may be used to improve the measurement of particle velocity u in future experiments.