Effects of Particle Migration on the Relaxation of Shock Wave Collisions

TL;DR

This paper explores how particle movement affects the relaxation process after shock wave collisions using molecular dynamics simulations.

Contribution

The study introduces a new perspective on shock relaxation by highlighting the role of particle migration in energy changes.

Findings

Particle migration significantly influences energy changes during shock relaxation.

The collision of two shock waves shows unique non-equilibrium features compared to single shock propagation.

Microscopic particle tracking reveals insights into the relaxation mechanisms of shock wave collisions.

Abstract

The non-equilibrium characteristics during the shock relaxation process hold a foundational position in various fields. In contrast to the propagation of a single shock wave, the collision process of two shock waves exhibits distinct non-equilibrium features. Employing non-equilibrium molecular dynamics, we simulated the collision of ultra-strong shock waves in a classical gas system, investigating the relationship between equilibrium relaxation time and shock intensity. Tracking the spatial migration of microscopic particles in the shock collision region during the relaxation process, we observed a significant contribution of particle migration to the average energy changes during relaxation. The discussion on particle migration provides a valuable new perspective for understanding the microscopic mechanisms of the relaxation process.

Figures

Full text

1. Introduction

Shock waves play a crucial role in diverse fields, including inertial confinement fusion (ICF) [1,2,3,4], astrophysics [5,6,7,8], hypersonic flight [9,10,11], and more [12,13,14,15]. Achieving ultra-high density and temperature for fusion fuels in ICF relies on generating superintense shock waves with powerful drivers like lasers or X-ray sources. Similarly, astrophysical investigations often require creating extreme states of matter through shock loading in labs. The non-equilibrium effects induced by shock waves significantly influence system evolution at the microscale. Understanding these dynamics is vital, especially in scenarios involving extreme states, where shock-induced non-equilibrium effects critically shape the system’s behavior.

The collision of two shock waves plays a crucial foundational role in issues such as the ICF shock ignition scheme [16,17,18] and the double cone collision ignition scheme [19]. The collision of two shock waves is a key process in elevating the density and temperature of fusion fuel. Non-equilibrium features in this shock wave collision may differ significantly from those observed in single shock wave scenarios. Current research on shock-induced non-equilibrium effects is primarily focused on individual shock waves, with limited investigation into non-equilibrium relaxation processes during the collision of two shock waves using microscopic methods.

The interaction between shock waves and matter will compress and heat the material in an extremely short period. For a stable, uniformly propagating shock wave, focusing our attention on a specific piece of material reveals that during the passage of the shock wave, the material undergoes a thermodynamic transition from one state, $[eqn]$ , to another state, $[eqn]$ . Consequently, there will be an equilibrium relaxation process from the $[eqn]$ state to the $[eqn]$ state. During this relaxation process, the material’s state deviates significantly from equilibrium, displaying strong non-equilibrium characteristics. Due to the steady and uniform propagation of the shock wave, when we establish the reference frame on the shock wave front and move with it, the equilibrium relaxation process experienced by the material over time will form a stable spatial distribution. This distribution can be referred to as the shock wave front structure.

The shock wave front structures operate on an extremely small spatial scale, comparable to the mean free path of molecules. The macroscopic continuum mechanics struggle to capture these structures. In continuum methods, shock waves are often depicted as infinitesimally thin discontinuities, with thermodynamic quantities interconnected through fluid conservation relations. In classical computational fluid dynamics, abrupt thermodynamic transitions across shock wave fronts cause problems like numerical oscillations, necessitating specialized treatments such as reducing computational scheme orders and introducing artificial viscosity [20]. A comprehensive understanding of shock waves necessitates a closer examination of their microscopic structure and the diverse non-equilibrium dynamic processes that occur within them.

A more accurate depiction of shock waves demands an exploration from the viewpoint of microscopic particles. Early investigations into the microscopic structure of shock wave fronts primarily employed kinetic theory methods, focusing on variations in the distribution function of microparticles under shock conditions. A pivotal theoretical contribution in this domain was the Bimodal approximation proposed by Mott-Smith [21] for the distribution function at shock wave fronts. This model subsequently laid the foundation for a series of subsequent studies [22,23,24]. Kinetic analyses are based on the Boltzmann equation, where the choice of collision models plays a crucial role in the system’s evolution. In scenarios where systems are sparsely populated and slightly deviate from equilibrium, collisional influences are minor, enabling kinetic analyses to effectively portray the evolution of distribution functions. However, for dense systems and higher degrees of non-equilibrium, the significance of collisional effects amplifies, rendering simplistic collision models inadequate in kinetic analyses.

Advancements in computer technology and molecular dynamics methods have liberated researchers from the constraints of kinetic theory, enabling more direct and detailed descriptions of particle collisions at the microscale [25,26,27]. Consequently, a series of computational studies employing the molecular dynamics method have emerged, focusing on the microscopic structure of shock wave fronts and the non-equilibrium processes induced by shock effects. Research endeavors span shock wave front thickness [28,29,30,31,32,33], distribution functions [24,34], and extend to more intricate microscopic dynamic phenomena including shock-induced phase transitions [35,36,37,38], molecular dissociation, ionization, and other complex phenomena [39,40].

In this work, we systematically simulated the collision process of strong shock waves in dense gas systems using the Non-Equilibrium Molecular Dynamics (NEMD) method. To form a clear and analyzable physical picture and eliminate the possible complex effects such as ionization, dissociation, and radiation under strong shocks, we restricted our study to a shock system that only includes a gas system with classical two-body collisions. In this paper, we conducted simulation studies using a classical helium gas system as a representative. The focus of this work is on the impact of shock wave intensity on non-equilibrium characteristics; thus, scenarios with different shock wave velocities were set. The selected range of shock wave velocities covers typical speeds observed in ICF implosion conditions.

This work will explore the relationship between the equilibrium relaxation time of various thermodynamic quantities and the shock velocity. By tracking and statistically analyzing the collision and motion of microscopic particles in the collision region of shock waves, the impact of particle collisions and migration mechanisms on the relaxation process will be examined and discussed. The statistical results indicate that as the shock velocity increases, the relaxation time monotonically decreases. Furthermore, the analysis of particle migration behavior suggests that particle migration plays a crucial role in accelerating the relaxation process.

The subsequent arrangement of this paper is outlined as follows. In Section 2, we will provide an exposition of the specific configurations of the molecular dynamics simulations and the statistical methodologies employed. Section 3 will be dedicated to the discussion of the obtained simulation results. Finally, Section 4 will encompass a comprehensive summary of the entire study.

2. Computational Details

2.1. NEMD Configurations

The non-equilibrium molecular dynamics (NEMD) calculations in this study were conducted using the LAMMPS [41] code (lammps-7Aug19). The simulated shock medium was helium, corresponding to an initial temperature of $[eqn]$ K (above the gas–liquid phase transition temperature of helium) and an initial density of $[eqn]$ g/cm^3^ (referencing parameters from Nellis’ shock experiments [42]). The system was equilibrated to the predefined initial state using the Nose–Hoover thermostat [43,44,45] and later set to the NVE ensemble in LAMMPS.

In the calculations, the X-axis was chosen as the direction of shock wave propagation. To drive a stable, propagating pair of shock waves within the system, we established a pair of reflecting boundaries as shock pistons. These pistons move uniformly in opposite directions along the x-axis at a predetermined piston velocity $[eqn]$ . To ensure sufficient evolution time for both shock waves to enter the stable propagation stage before colliding, the computational box was configured as an elongated region with dimensions of 4065 Å(x) × 203 Å(y) × 203 Å(z), containing over 2.5 million atoms, as depicted in Figure 1.

In this study, the classical Lennard-Jones (LJ) potential function is employed to describe the pairwise interactions between helium atoms, which could be written as

$$[eqn]$$

where r is the distance between particles, $[eqn]$ is the strength of potential energy, and $[eqn]$ is the distance scale between particles. According to LJ potential [46], the values of $[eqn]$ and $[eqn]$ are set as $[eqn]$ kcal/mol and $[eqn]$ Å. As mentioned in the introduction, the subject of this study is a classical gas dominated by two-body collisions. Therefore, the application of the Lennard-Jones (LJ) potential can reveal the main physical picture. For the primary calculations, we also performed repeated simulations using the Beck [47] and Morse [48] pair potentials for helium. The comparison shows that the computational results are consistent across different potential functions.

Periodic boundary conditions were applied along the y-axis and z-axis directions to minimize the influence of lateral dimensions’ boundary effects and size effects on the propagation of shock waves [49]. The coordinate origin was set at the geometric center of the box. With the left and right pistons moving at the same velocity, the system possessed mirror symmetry about the x=0 plane. This allowed the precise determination that the collision point of the shock waves occurred within the central region of the system’s geometry.

In this work, we selected a range of piston velocities $[eqn]$ from 10 km/s to 100 km/s. As time evolves, the piston will drive two shock waves propagating in opposite directions within the system. The propagation speed of the first shock wave is denoted as $[eqn]$ , which causes the material behind the shock to move at the piston velocity, and the thermodynamic state variables of density, temperature, and pressure change to $[eqn]$ , $[eqn]$ , and $[eqn]$ , respectively. When the two shock waves collide in the central region, they generate a secondary shock wave propagating outward with a speed of $[eqn]$ , leading to changes in the thermodynamic state variables to $[eqn]$ , $[eqn]$ , and $[eqn]$ . Table 1 presents the thermodynamic state variables behind the shock front and shock wave speeds for both primary and secondary shocks in all cases considered in this work.

Since the chosen range of shock piston velocities for simulations spans a wide range, in order to maintain the continuity of atomic trajectories and use larger time steps $[eqn]$ whenever possible to conserve computational resources, the following equation is employed to estimate the optimized time step for different shock cases,

$$[eqn]$$

where $[eqn]$ and T represent the number density of particles and the temperature in the downstream region of the shock wave, respectively. $[eqn]$ is the Boltzmann constant, and $[eqn]$ denotes the mass of helium atoms. Although the number density $[eqn]$ and temperature T in the downstream region require accurate values from shock simulations, based on approximate relationships summarized in previous work [49], we can make a rough estimate of their likely ranges. In essence, the derived optimized time step $[eqn]$ from the equation above is approximately 1/80 of the time taken for helium atoms to traverse the average distance between adjacent particles at their average thermal velocity. Across all simulation cases in this study, the time step ranges from $[eqn]$ fs to $[eqn]$ fs.

2.2. Thermodynamic Statistical Methodologies

The results of MD simulations provide the temporal evolution of kinematic quantities for a multitude of microscopic particles. To acquire corresponding macroscopic thermodynamic quantities and delve into the evolution of thermodynamic states during non-equilibrium relaxation processes, statistical analysis of these microscopic particles is essential. Considering the propagation of the shock wave along the x-axis, while the distribution of thermodynamic quantities is uniform along the y and z directions, the computational domain is divided into numerous thin slices along the x-axis for thermodynamic averaging, as illustrated in Figure 1. The thickness of each slice is denoted as $[eqn]$ . The selection of slice thickness $[eqn]$ necessitates the consideration of two interrelated factors. On one hand, the slice should be sufficiently thick to ensure an ample number of microscopic particles within this region, meeting the fundamental requirements of thermodynamic statistics. On the other hand, given the small spatial scale of the shock process, if the slice is excessively thick, the thermodynamic quantities obtained from statistics might lack locality, rendering them incapable of discerning the details of equilibrium relaxation processes. In this study, a slice thickness of 1 Å was ultimately chosen.

For a specific slice, the statistical formula for various thermodynamic quantities is as follows [32,50]:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , u, T, and P represent the density, velocity, dynamic temperature, and pressure corresponding to the slice, respectively. N is the total number of particles in the slice, and $[eqn]$ is the volume of the slice. $[eqn]$ , $[eqn]$ , and $[eqn]$ represent the components of the velocity vector of atom i along the x, y, and z directions, respectively. $[eqn]$ represents the force generated between atom i and atom j due to the interaction potential, $[eqn]$ and $[eqn]$ are the position vector of atom i and atom j, respectively.

In the process of non-equilibrium relaxation, significant non-equilibrium effects will be observed, and the kinetic energy of microparticles’ thermal motion exhibits noticeable anisotropy. To quantitatively describe this non-equilibrium characteristic, based on the fundamental principles defined by dynamic temperature, we further decompose the dynamic temperature into temperature components, denoted as $[eqn]$ ( $[eqn]$ ) and $[eqn]$ (the average of $[eqn]$ and $[eqn]$ ), parallel and perpendicular to the direction of the shock, respectively. The corresponding formula for $[eqn]$ and $[eqn]$ can be written as [49]

$$[eqn]$$ $$[eqn]$$

It is worth noting that the concept of temperature requires a basis in thermodynamic equilibrium. In the non-equilibrium region of a shock wave, the ‘temperature’ and its components obtained from Equations (5), (7), and (8) cannot be directly interpreted as the temperature concept in equilibrium states. Instead, they manifest as parameters describing the characteristics of particle distribution, reflecting the deviation from Maxwell distribution under non-equilibrium conditions.

3. Results and Discussion

3.1. Overshoot Structure and Velocity Distribution

The shock wave front is notably characterized by its non-equilibrium nature, reflected on a microscopic level through deviations in the distribution function from the Maxwellian distribution. This non-equilibrium behavior extends to macroscopic thermodynamic quantities, such as the overshoot phenomenon in dynamic temperature distribution. Previous studies have extensively explored these aspects at individual shock wave fronts [24,31,34,49].

In this study, we further investigate non-equilibrium features in shock wave collisions by examining the evolution of distribution functions and the occurrence of the overshoot phenomenon during the collision of two shock waves. Building upon prior research, we focus on the most pronounced overshoot feature observed in the dynamic temperature component along the shock direction, denoted as $[eqn]$ . To illustrate, we analyze a case with a shock piston velocity of 10 km/s, presenting distribution plots of $[eqn]$ before and after the shock wave collision (see Figure 2).

The overshoot of $[eqn]$ is manifested as a peak structure in its spatial distribution profile, as illustrated in Figure 2. Notably, the peak structure caused by a single shock wave before collision is more pronounced than the one induced by the secondary shock wave after the collision. Following the methodology of prior work [49], we quantify the overshoot using the dimensionless parameter $[eqn]$ , representing the ratio of the peak value $[eqn]$ to the post-wave average value $[eqn]$ . Figure 3 presents the values of $[eqn]$ at the shock wave front for both single ( $[eqn]$ ) and secondary ( $[eqn]$ ) shock waves at different shock piston velocities, taking the situation illustrated in Figure 2 as an example, where $[eqn]$ and $[eqn]$ . In the case of a single shock wave, the results align with previous research [49], with $[eqn]$ fluctuating within a small range around 1.5. However, in the case of the secondary shock wave after collision, $[eqn]$ varies within a small range near 1.33, significantly lower than the single shock wave scenario.

The primary reason for this difference lies in the substantial thermodynamic differences between the shock front and the post-wave state during a single shock wave, with a density compression ratio close to 4 and a temperature increase by a factor of several thousand. In contrast, for the secondary shock wave after collision, the thermodynamic differences between the post-wave and pre-wave states decrease, with a density ratio of approximately 2 and a temperature ratio of about 4.

The overshoot phenomenon arises from the system undergoing a relaxation process, mirroring the evolution of the microparticle distribution function. We delve into this process by examining the particle velocity distribution function under shock waves, comparing functions at the overshoot peak and in the post-wave equilibrium region to understand the origins of the overshoot.

In Figure 4, panel (a) depicts velocity distribution function curves for $[eqn]$ at various positions under a single shock wave. The vertical lines correspond to average velocity values for the colored distribution functions. At the overshoot peak, the distribution function displays a high, narrow, low-temperature peak and a broad, high-energy tail, with the average velocity notably deviating from the low-temperature peak. The calculation formula for dynamic temperature Equation (5) indicates that this “low-temperature” particle portion significantly contributes to the overall temperature after accounting for translational velocity. Contrastingly, in panel (b), the distribution function at the overshoot peak on the shock wave front after the shock wave collision is closer to the Maxwell velocity distribution curve and the average velocity falls within the central peak region of the velocity distribution function. This results in a weaker overshoot effect compared to the single shock wave scenario.

The discussed overshoot phenomenon in temperature is not only observable in spatial distribution at a given moment but also in the temporal evolution within a specified region. Figure 5 illustrates the evolution of $[eqn]$ in the collision slice (−0.5 Å $[eqn]$ Å) over time for piston speeds $[eqn]$ ranging from 20 to 40 km/s. The curves are normalized using the downstream equilibrium temperature for ease of comparison, revealing that the overshoot proportion remains relatively consistent across different piston speeds, while the equilibrium relaxation time decreases with higher shock intensity.

In Figure 6, the blue dashed line depicts the particle velocity distribution function at the overshoot moment $[eqn]$ of $[eqn]$ within the collision region. The green dotted line represents the particle velocity distribution function at equilibrium after the collision at time $[eqn]$ , with (a) and (b) corresponding to the velocity distribution functions of $[eqn]$ and $[eqn]$ , respectively.

Firstly, in comparison to the equilibrium distribution function after the collision, the velocity distribution functions of $[eqn]$ and $[eqn]$ at the moment $[eqn]$ exhibit distinct low-temperature peak structures at v = 0 compared to the Maxwell distribution function at time $[eqn]$ . This suggests the presence of numerous low-temperature particles in the region at $[eqn]$ . Meanwhile, in Figure 6a, the blue curve on both sides in the high-speed range is higher than the green dotted line. This phenomenon is attributed to high-temperature particles behind shock waves carrying not only substantial thermal kinetic energy but also significant translational kinetic energy along the shock direction. At $[eqn]$ , this translational kinetic energy is primarily retained in the x-direction, resulting in a peak in the statistical analysis of the $[eqn]$ component. As the relaxation process unfolds, this translational kinetic energy gradually disperses uniformly in all three directions through particle collisions and scattering, ultimately reaching isotropic thermodynamic temperature. Consequently, $[eqn]$ decreases from its peak value.

In contrast, in Figure 6b, the transverse velocity distribution function in the high-speed range on both sides shows that the blue dashed line is lower than the green dotted line. There is no overshoot structure in the statistically obtained $[eqn]$ since there is no translational kinetic energy induced by the shock in the transverse direction.

3.2. Relaxation Time under Shock Collision

This work primarily investigates relaxation processes in shock wave collisions. Notably, the shock system undergoes distinct thermodynamic state changes: the first shock ( $[eqn]$ to $[eqn]$ ) and the second shock ( $[eqn]$ to $[eqn]$ ). A variable period of equilibrium at $[eqn]$ occurs between these shocks in most regions. However, at the central collision region, the transition is direct from $[eqn]$ to $[eqn]$ , causing notable differences in the equilibrium relaxation process compared to a single shock wave case. Our analysis will focus on the temporal evolution of thermodynamic quantities during the shock collision, exploring their interrelationships and connections to shock intensity.

Figure 7 illustrates the time evolution of density ( $[eqn]$ ), temperature (T), and parallel temperature components ( $[eqn]$ ) in the shock collision region, with a piston velocity ( $[eqn]$ ) of 20 km/s. Similar to the single shock scenario, the density curve lags behind the temperature increase, indicating that shock heating precedes shock compression. Due to distinct relaxation processes under shock, different thermodynamic quantities exhibit varying relaxation times.

Here, we obtain the relaxation time $[eqn]$ using the method illustrated in Figure 7. For example, the relaxation time ( $[eqn]$ ) for density $[eqn]$ is determined by considering the average density values before and after collision as reference points. Horizontal lines are drawn at these equilibrium states, and the deviation points on the density curve mark the start ( $[eqn]$ ) and end ( $[eqn]$ ) of the relaxation process, with the time difference representing $[eqn]$ . This study conducts a statistical analysis of relaxation processes for density, temperature, and its components at different shock velocities, presenting the obtained relaxation times in Table 2.

To visually depict the variation in relaxation times with shock piston velocity, we generated plots for the data. It is evident that with increasing shock intensity, the relaxation times of various thermodynamic parameters consistently decrease, as shown in Figure 8. Analyzing the decreasing trend in the $[eqn]$ - $[eqn]$ curve, it adheres to the characteristics of an inverse proportionality function. If the hypothesis of this inverse proportional function holds true, then the relationship between $[eqn]$ and $[eqn]$ can be explained as follows: the thickness of the shock wave collision region remains constant, and the propagation speed of the shock wave is approximately linearly related to $[eqn]$ . In this case, the relaxation time $[eqn]$ can be approximated as the time it takes for the shock wave to traverse this fixed thickness or is linearly related to this time.

To further clarify the relationship between relaxation times ( $[eqn]$ ) and shock piston velocity ( $[eqn]$ ), we initially attempted a direct fit of the data from Table 2 using the inverse proportionality function $[eqn]$ . However, the results deviated significantly from the MD simulation data points. This deviation indicates that the determination of the relaxation time for shock wave collisions cannot be explained by the simplified physical picture mentioned above alone. It necessitates consideration of the influences of more complex effects.

To improve fitting accuracy, we introduced a corrective term to the inverse proportionality function, resulting in the following fitting equation:

$$[eqn]$$

where a, b, and c represent the fitting parameters. The final fitting outcomes closely align with the MD simulation results, as depicted in Figure 8. The values of the ultimate fitting parameters, a, b, and c, for each thermodynamic parameter are provided in Table 3.

3.3. Particle Migration During Relaxation

In the following discussion, we focus primarily on two aspects. First, we explore the spatial migration patterns of particles during the relaxation process, with the relevant results shown in Figure 9 and Figure 10. Next, we examine the influence of particle migration on the energy changes of particles within the shock collision slice, with the corresponding results presented in Figure 11, Figure 12 and Figure 13.

To simplify the discussion, the following symbols will be used in the subsequent discourse. Before delving into the discussion, below is an explanation of the specific meanings associated with the symbols:

• $[eqn]$ : the starting point of the relaxation process.
• $[eqn]$ : the ending point of the relaxation process.
• $[eqn]$ : the moment when $[eqn]$ reaches the peak of the overshoot.
• G_0_ or G_1_: the group of particles located within the collision slice at time $[eqn]$ or $[eqn]$ .
• $[eqn]$ or $[eqn]$ : the number of G_0_ or G_1_ particles in a slice.
• N: the total number of all particles in a slice.

During the relaxation process, changes in macroscopic thermodynamic quantities within the collision slice correspond to variations in microscopic particle kinematics, a transformation typically accomplished through particle collisions. However, it is crucial to note that particles in G_0_ do not remain confined to the collision slice throughout the entire relaxation process; they disperse into a broader space over time. Concurrently, particles outside the collision slice at time $[eqn]$ may also enter the collision slice during relaxation, contributing significantly to the thermodynamic state change of this slice. In addition to the contribution from particle collisions, the entry and exit of these particles play a notable role in the transformation of the thermodynamic state of the collision slice.

In the context of thermodynamic equilibrium, the movement of particles in and out of a region is primarily governed by random thermal motion, leading to particle diffusion. However, in shock wave collision systems, particles near the collision region undergo not only random thermal motion but also translational acceleration along the shock direction due to the shock effect. In subsequent discussions, we term this spatial particle change during the relaxation process as the “migration” process.

Firstly, by tracking the spatial positions of particles in G_0_, we can obtain the temporal evolution of the spatial distribution function of these particles, as shown in Figure 9. It can be observed that, for a case with $[eqn]$ 10 km/s, within a relaxation time of 226 fs, the spatial distribution of particles in G_0_ evolves from a uniform distribution within a thin slice of thickness 1 Å to a Gaussian distribution with a half-width of approximately 7.9 Å. By comparing the results of cases with different $[eqn]$ , it can be observed that with an increase in shock intensity, the spatial distribution of particles at time $[eqn]$ will broaden to a larger range.

Similar to the above analysis, we incorporate particles within the slice at $[eqn]$ into the group G_1_ and attempt to explore the contributions of particle collisions and migration to the relaxation process by tracking the spatial positions and kinematic changes of these particles in the reverse direction of time. Figure 10a displays the spatial distribution function curves of G_1_ particles at different times during the relaxation process. Unlike the situation with G_0_, the spatial distribution function at time $[eqn]$ for G_1_ does not exhibit a Gaussian distribution but instead shows three distinct peak structures, as illustrated by the red curve in Figure 10a. The peak at $[eqn]$ is due to the same reasons that lead to the peak structure of the Gaussian distribution in G_0_. Under random walk motion, the probability of a particle remaining in place at the ending time is always the highest.

The peaks on the left and right sides at time $[eqn]$ are mainly due to the non-uniformity of the particle number density in space at time $[eqn]$ . At time $[eqn]$ , the two shock waves have not yet fully collided; behind the shock wave front corresponds to the compressed high-density region, while the central slice, because the shock waves have not yet arrived, is in a low-density state. Therefore, the spatial profile of the particle number density will form a trough at the central slice at time $[eqn]$ , as shown in Figure 10b. Despite the fact that, according to the random walk model, particles farther from the collision slice at $[eqn]$ have a lower probability of migrating to the collision slice at $[eqn]$ , the higher particle density in the compressed regions on both sides ultimately leads to an increased number of particles migrating into the collision slice.

We calculate the probability of individual particles at different x at $[eqn]$ migrating to the central collision slice at $[eqn]$ by dividing the number of G_1_ particles within the slice at position x at $[eqn]$ by the total number of particles in the slice. This probability distribution with respect to spatial location is shown in Figure 10c. The probability distribution closely resembles a Gaussian curve, indicating that, from the perspective of time reversal, the particle migration process in the central collision slice during relaxation exhibits characteristics consistent with random walk motion [51], despite substantial changes in the macroscopic thermodynamic state of these particles throughout the relaxation process.

3.4. Changes in Particle Kinetic Energy during the Relaxation Process

The following provides explanations for the symbols that will be used in the upcoming discussion:

• $[eqn]$ or $[eqn]$ : the average kinetic energy of particles in G_0_ or G_1_.
• $[eqn]$ or $[eqn]$ : the component of average kinetic energy along axis x (the parallel direction) for particles in G_0_ or G_1_.
• $[eqn]$ or $[eqn]$ : the component of average kinetic energy along axis y (the perpendicular direction) for particles in G_0_ or G_1_.
• $[eqn]$ : the average kinetic energy of particles located within the collision slice.

Next, we aim to analyze the impact of particle collisions and migration on the relaxation process by examining changes in particle kinetic energy. In Figure 11a, the time evolution curves of $[eqn]$ (red curve) and its directional components $[eqn]$ (blue curve) and $[eqn]$ (green curve) from $[eqn]$ to $[eqn]$ are presented, alongside the time evolution curve of $[eqn]$ (black curve).

The trend of $[eqn]$ reveals a continuous increase in the kinetic energy of G_1_ particles throughout the relaxation process, achieved through collisions with other particles in the spatial region. The blue curve representing $[eqn]$ initially increases and then decreases during relaxation. This behavior is primarily due to the presence of some G_1_ particles in the post-shock wave region at $[eqn]$ , carrying both thermal and translational kinetic energy in the shock direction. As the relaxation progresses, collisions between particles moving in opposite directions result in a transfer of momentum from the parallel direction (x) to the perpendicular direction (y and z). This transition leads to a change in kinetic energy from parallel to perpendicular until reaching equilibrium with an isotropic distribution of kinetic energy components.

The above analysis of $[eqn]$ is supported by the statistical velocity distribution function of G_1_ particles in Figure 13. The velocity distribution function of G_1_ particles at time $[eqn]$ in the collision slice corresponds to the blue curve in Figure 13. Since at this time, G_1_ particles are distributed within a wide range of −30 Å to 30 Å in space, including both the pre-shock wave and post-shock wave regions, their velocity distribution function in the shock direction can be approximately regarded as a superposition of the cold Maxwell distribution function in the central region and two high-temperature Maxwell velocity distribution functions on both sides with translational velocities. As the relaxation process progresses, the peak value of the cold Maxwell peak in the central region gradually decreases, corresponding to the heating of the central region by the collision of the shock waves.

Due to particle migration, the energy evolution of G_1_ particles ( $[eqn]$ ) will show a significant difference from that of particles in the central collision slice ( $[eqn]$ ), and we can discuss this difference to understand the specific impact of migration on the relaxation process. In Figure 11a, the black dashed line shows the change in the average kinetic energy of central slice particles $[eqn]$ over time. Since the central slice has not been affected by the shock waves in the early stages of relaxation and has a lower temperature, its average kinetic energy value at time $[eqn]$ is much smaller than that of G_1_. As the relaxation process progresses, the black dashed curve gradually rises and coincides with the red curve at time $[eqn]$ .

In this process, on the one hand, particles within the collision slice increase their energy through collisions with high-energy particles from the central region. On the other hand, particle migration also makes a significant contribution to the increase in average kinetic energy. Here, the role of migration in increasing average kinetic energy $[eqn]$ includes both the contribution of external high-energy particles entering the slice and the contribution of low-temperature particles in the central slice migrating out of the slice. It should be noted that the specific particle migration process exhibits strong randomness. During the relaxation process, particles entering through migration may also be low-energy particles, and those leaving may be high-energy particles, thereby contributing to a decrease in the average particle kinetic energy. However, from a statistical perspective, the average kinetic energy of particles within the collision slice continues to increase during the relaxation process. The contribution of particle migration to the average particle kinetic energy, as obtained through statistical analysis, is also positive.

The change in $[eqn]$ from $[eqn]$ to $[eqn]$ is approximated as the contribution of particle collisions to the average particle energy change, denoted as $[eqn]$ . Meanwhile, the change in $[eqn]$ is considered the total change in average particle energy during the relaxation process, denoted as $[eqn]$ . We subtract from the two yields the change caused by migration, denoted as $[eqn]$ . This relationship can be expressed by the following formulas as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

By dividing the contributions of collision and migration by the total change, we can quantitatively assess the proportion of the contributions of particle migration and collision to the particle kinetic energy during the relaxation process. We systematically calculated the proportions for different $[eqn]$ cases, and the results are summarized in Table 4. It can be observed that the proportion of particle migration effects remains relatively constant with increasing $[eqn]$ , consistently exceeding 50%. This indicates that the impact of particle migration on the relaxation process is quite significant.

To quantify the impact of migration and collisions on the relaxation process, we calculated the time derivative of the average kinetic energy and its components, as shown in Figure 11b. The black dashed line represents the rate of energy change $[eqn]$ from both collisions and migrations. Initially, the energy change rate is low but peaks around 1160 fs when the shock waves begin to collide.

The red solid line shows the rate of change in energy due to particle collisions ( $[eqn]$ ), which remains stable throughout the relaxation period, only dropping to zero in the late phase.

Since the energy from collisions is distributed as either thermal motion kinetic energy or translational kinetic energy, we analyzed how the added energy is allocated between these two forms. Assuming isotropy in the thermal motion components, the thermal kinetic energy of G_1_ particles is approximated as $[eqn]$ and the translational kinetic energy as $[eqn]$ . By calculating their rates of change over time, shown in Figure 12, we observe that before 1200 fs, the rates for both energy types are similar, indicating that the additional energy from collisions is proportionally distributed between translational and thermal motion kinetic energy. This behavior aligns with that of a single shock wave in ideal gas [49]. After 1200 fs, the translational kinetic energy decreases as the shock waves begin to collide directly, converting the lost translational energy into thermal energy for G_1_ particles.

4. Conclusions

This study systematically investigates the collision process of strong shock waves propagating through a classical helium medium using the Non-Equilibrium Molecular Dynamics (NEMD) method. The simulation results reveal a decrease in relaxation times for various thermodynamic quantities as the shock piston velocity increases, and this relationship is accurately modeled by an inverse proportional function. We extensively discuss the influence of particle migration and collisions on the change in particle kinetic energy during the relaxation process. The findings suggest that particle migration significantly accelerates the relaxation process. Unlike typical relaxation process studies that primarily focus on microscopic particle collision mechanisms, the results of this work indicate that in the relaxation process during strong shock wave collisions, particle migration contributes significantly alongside particle collision mechanisms. This study provides valuable insights into the fundamental physics of shock wave collisions, particularly the discussion on particle migration effects, offering a novel perspective on understanding the evolution of relaxation processes in shock wave collisions.

Considering the substantial simplifications made in this work regarding the strong shock wave collision process, a more accurate quantitative assessment of the impact of particle migration on relaxation processes during strong shock wave collisions in real-world scenarios would require consideration of factors such as electron ionization and the potential function under extreme conditions. Building upon the results of this research, incorporating a deeper understanding of additional complex factors will allow for a more detailed exploration of the contributions of different factors to the relaxation process. This could serve as a basis for future research endeavors.