Theory and simulation of shock waves: Entropyproduction and energy conversion

We have considered a shock wave as a surface of discontinuity and computed the entropy production using non-equilibrium thermodynamics for surfaces. The results from this method, which we call the"Gibbs excess method"(GEM), were compared with results from three alternative methods, all based on the entropy balance in the shock front region, but with different assumptions about local equilibrium. Non-equilibrium molecular dynamics (NEMD) simulations were used to simulate a thermal blast in a one-component gas consisting of particles interacting with the Lennard-Jones/spline potential. This provided data for the theoretical analysis. Two cases were studied, a weak shock with Mach number $M \approx 2$ and a strong shock with $M \approx 6$ and with a Prandtl number of the gas $Pr \approx 1.4$ in both cases. The four theoretical methods gave consistent results for the time-dependent surface excess entropy production for both Mach numbers. The internal energy was found to deviate only slightly from equilibrium values in the shock front. The pressure profile was found to be consistent with the Navier-Stokes equations. The entropy production in the weak and strong shocks were approximately proportional to the square of the Mach number and decayed with time at approximately the same relative rate. In both cases, some 97 \% of the total entropy production in the gas occurred in the shock wave. The GEM showed that most of the shock's kinetic energy was converted reversibly into enthalpy and entropy, and a small amount was dissipated as produced entropy. The shock waves traveled at almost constant speed and we found that the overpressure determined from NEMD simulations agreed well with the Rankine-Hugoniot conditions for steady-state shocks.


Introduction
The amount of energy carried by a shock wave is considerable, and the wave travels at supersonic speed. Shock waves produced from explosions, rapid phase transitions, sudden release of pressurized gas, or other blasts, are highly irreversible phenomena. Shock waves are therefore both interesting and challenging phenomena to understand and quantify. Several laboratory and large-scale field experiments have been carried out to determine the impact of blast waves as function of explosion type and strength, distance from the blast, and topology of the surroundings [1].
The basic equations describing the conditions for shock-waves in one dimension were developed in the late 19th century by Rankine and Hugoniot [2,3]. In these early studies, the shock wave was considered as a surface of discontinuity with conservation of mass, momentum, and energy. The shock's properties were described in terms of the differences between the upstream and downstream properties of the bulk fluids (see e.g. Hirschfelder et al. [4], Hoover and Hoover [5], and Uribe [6]). The Rankine-Hugoniot relations give a macroscopic description of the state variables in front of, and behind, the wave front, but not the details of e.g. energy dissipation at the front. Application of hydrodynamic theories in the early 20th century gave more details of the shock-wave front, such as its thickness [7,8]. The developments of kinetic theories at about the same time supported and examined the limitations of the Navier-Stokes equations as applied to shock waves [8]. It was found that the thickness of the shock-wave front given by the Navier-Stokes equations was too small compared with experiments and improved theories [9]. Questions were also raised about the consistency between the entropy profile showing a peak at the front position and the second law of thermodynamics [10]. The interest in blast-wave theory and experiments was high during and after the second world war, which led to significant progress in the understanding of shock waves [11,12,13]. The more recent progress in hydrodynamic theory [14], kinetic theory [15,16], and extended thermodynamics [17,18] has given substantial new insight into many properties of shock waves (see e.g. Section 9 in García-Colín et. al [19] for a good review).
The development of Direct Simulation Monte Carlo (DSMC) and molecular dynamics (MD) simulations provided the necessary link between experiments and theories. Klimenko and Dremin [20], Hoover [21], Holian [22], Salomons [23], and their coworkers were pioneers with the aims to clarify some of the puzzles in shock-wave theory and demonstrate the applicability of computer simulations. For instance, the exact thickness of a shock wave had been estimated [15], but was not exactly known until simulations produced accurate results [21,23]. With the increase in compute power, very detailed analyses can now be made [23,24,25,26]. The status is that several theories work for shocks with Mach numbers up to approximately 2 (which may be defined as "weak shocks") whereas theories for "strong shocks" (with higher Mach numbers) are still being developed.
Despite the fact that shock-wave propagation is an irreversible process, few papers have been concerned with the energy dissipation and entropy production in shock waves. Brinkley and Kirkwood presented a theory of non-steady shock waves in 1947 which included the concept of energy dissipation and wave speed retardation [27]. At about the same time, Tolman and Fine published a comprehensive paper on entropy production in irreversible processes, including shock waves [28]. It has been shown that kinetic energy is not equipartitioned in the shock-wave front [29], which is a good reason to question the local-equilibrium assumption made in non-equilibrium thermodynamics (NET). Nevertheless, this assumption was adopted by Velasco and Uribe, who used the Gibbs equation in the normal way for bulk fluids to obtain the entropy production in the shock front [30]. By introducing empirical temperature dependencies of the viscosity in combination with the Navier-Stokes-Fourier relations, they got good agreements with results from DSMC for Mach numbers between 1.55 and 9. There are, however, remaining questions, such as exactly how the kinetic and compression energy carried by a shock wave is dissipated or converted to other forms, in particular when the wave hits an obstacle or a body. Such questions are important for studying impact of detonations [31], in material science [32], formation and collapse of bubbles [33], and traumatic brain injuries from improvised explosive devices (IEDs) [1], to mention a few.
All approaches to shock-wave analyses use conservation of mass, momentum, and energy (see e.g. Landau and Lifshitz, ref), shown here for a plane shock wave moving with constant velocity in x-direction in a one-component, single-phase fluid: where ρ is the mass density, v and v s are the streaming velocity and the shock-wave velocity, respectively, in x-direction in the stationary frame of reference, P xx the pressure tensor component normal to the shock-wave front, u the internal energy per unit mass, and J q,x is the measurable heat flux in x-direction. For simplicity, we have not included gravity in these equations. Under these conditions, the c i 's are constants. In the classical treatment of shock waves in Newtonian fluids, Navier-Stokes and Fourier constitutive relations are introduced into the conservation laws (in addition to an equation of state) [23], viz. (4) P y y = P zz =p + where P y y = P zz is the pressure parallel with the shock front, p is the equilibrium pressure as given by the equation of state at the local conditions, η S and η B are the shear and bulk viscosities, respectively, λ the thermal conductivity, and T is the temperature. All the quantities in Eqs. (4) -(6) are in general functions of position x and time t . In the framework of Navier-Stokes, it follows from Eqs. (4) and (5) that if the viscous terms are small, the diagonal components of the pressure tensor are approximately equal to the equilibrium pressure. A shock wave may be characterized by a sharp front with significant changes in density and pressure over such a short distance that, at a macroscopic scale, it can be considered to be a discontinuity in the system's characteristic properties [4]. This is not unlike the case of a regular surface, e. g. a liquid/vapor surface. NET for surfaces has been developed by Kjelstrup and Bedeaux [34] and we recently reported that NET can successfully be applied to a weak shock wave (Mach number 2.1) [35]. In the present work, we develop the method further, analyze its basis in more detail, and use the results to describe the excess entropy production and energy conversions on the surface. Furthermore, we find a flux-force relation for mass transfer across the shock front in the surface description. The work presented here is a new approach to assess the properties of shock waves. This work includes data for a significantly stronger shock wave than the one we considered in our previous analysis [35].
We have used four different methods to determine the surface excess entropy production. The methods are based on different assumptions and the consistency between the methods are used as a criterion for the validity of their underlying assumptions. The concept of local equilibrium must be defined in different ways depending on the context of the methods we apply. All four methods are based on the balance equation for entropy across the shock-wave front. In the "bulk balance method" (BBM) and the linear irreversible thermodynamics (LIT) method we integrate the local entropy production over the shock wave thickness. The BBM uses the entropy balance directly whereas the LIT involve the Gibbs equation as well. The "surface balance method" (SBM) models the wave front as a surface of discontinuity and considers the entropy balance on the surface. In the "Gibbs excess method" (GEM), we use the Gibbs equation for surface variables [34,36] and derive a new, more detailed expression for entropy production and a new tool to analyze the energy conversions in a shock wave. The basic assumption in GEM is that this Gibbs equation for surface excess variables remains valid also when the system at large is out of equilibrium, as suggested by Bedeaux, Albano and Mazur in their construction of NET for heterogeneous systems [37,38]. Energy dissipation and conversion can be determined from this analysis, which leads to new information about energy conversion at the shock front. We show that a particularly simple constitutive relation can be used for shocks at, or close to, steady state. This relation has not been used in earlier work. At present, we restrict the discussion to one-component fluids.
The analysis starts with the balance equation for entropy in Section 2. The BBM and LIT are based on this balance. In Section 3 we establish the framework for the surface description of shock waves. This includes the definition of the Gibbs equimolar surface and reformulation of the balance equations and conservation laws in the surface description. The SBM is based on this surface description. The new relations based on NET are derived in Section 4 where we specify the conditions and assumptions made. The key result in this section is an expression for the surface excess entropy production. In order to quantify the theoretical results, we have done non-equilibrium molecular dynamics (NEMD) simulations of two shock waves, one weak and one strong (Mach number M ≈ 2 and M ≈ 6, respectively). Section 5 includes a description Figure 1: Illustration of the entropy density, ρ s , (solid line) and the entropy production, σ s , (dashed line) around a shockwave moving from left to right with velocity v s . The dash-dot-line is a macroscopic illustration of the shockwave with discontinuous fluid properties. of the model system used in the NEMD simulations and how the simulations were carried out. We show in Section 6 that NEMD is well suited to analyze the theories and provide unbiased data. The combination of NET theory and NEMD data lead to a new insight into the various contributions to the entropy production in a shock wave. Section 6 also includes descriptions of how the NEMD data were used to examine the assumptions mentioned above (Sections 6.2 and 6.1). Finally, we summarize the conclusions from this work in Section 7.

Basic entropy balance and the meanings of "local"
In addition to the conservation laws for mass, momentum, and energy, Eqs. (1)-(3), we consider the entropy over an infinitesimal control volume in one dimension: where ρ s , J s , and σ s are the density, flux, and production of entropy, respectively. Profiles of the entropy density and entropy production, centered around a shock wave front, are sketched in Figure 1. The entropy flux is a combination of heat conduction and entropy transported with the fluid flow, where v is the local fluid velocity. Whereas the heat conduction is independent of the frame of reference, the transported entropy is not, and is here given in the stationary frame of reference. In Eqs. (7) and (8) we have used the term "local" in the meaning "infinitesimally small domains in space and time". In the following, we shall use "local" in two different meanings, depending on whether we consider the microscopic picture and the wave front as a spatial domain with continuously varying properties, or the macroscopic picture and the wave front as a surface. "Local" in the former context means a small control volume, which in the numerical work is determined by the thickness of each layer in the NEMD simulations, typically about 30 molecular diameters or slightly more than the mean free path in the equilibrium gas ahead of the shock. "Local" in the latter sense means "on the surface", i.e. as given by the surface excess properties. The surface has no extension in x-direction. In both contexts, "local" in time means a time interval determined by the sampling time in the numerical simulations.
We shall now consider the four routes to the surface excess entropy production in the shock wave based on Eq. (7). The basic route is a direct integration of Eq. (7), which we call the "bulk balance method" (BBM). The results from BBM are discussed in Section 6.3. In two other routes we consider the shock wave front as a surface and apply NET for surfaces where we make use of the equilibrium bulk properties at both sides of the wave. The fourth method is the classic linear irreversible thermodynamics (LIT) method based on the assumption of local equilibrium and the Gibbs equation for bulk systems [36]. When applied to each local control volume in the system, this leads in the present context to where Π xx is the x-component of the viscous pressure tensor. The LIT method was recently used by Velasco and Uribe in an analysis of shock waves [30].

The shockwave as a surface
A typical surface has a thickness that is small compared to the thickness of the adjacent homogeneous phases, and the surface appears to be two-dimensional. In reality, a shock wave is several molecular diameters thick, depending on the system's thermodynamic state and speed of the wave [15], in this study about 100 molecular diameters. At the macroscopic scale, transport of heat and matter across a surface will give rise to discontinuities in intensive variables like temperature and chemical potential. The fluxes and forces may also become discontinuous. Our derivation of NET for shock waves builds on Gibbs' definition of a surface [39], and the assumption first made by Bedeaux, Albano and Mazur [37,38], viz. that thermodynamic relations between surface variables remain valid locally, also when the system overall is out of equilibrium. This assumption means that we define the interface as a separate and autonomous thermodynamic system [35]. The surface is assumed to possess a temperature, chemical potential, and other thermodynamic variables of its own. The assumption may seem drastic because the shock-wave front is in a non-equilibrium state without a corresponding equilibrium state (in contrast to the case for e.g. a liquid-vapor surface). We shall therefore examine this assumption in detail in Section 6.3.3. Thermodynamic properties of surfaces are well defined using Gibbs' surface excess densities of mass, entropy and energy [39]. Following the systematic procedure given by Albano et al. [38], we first derive the entropy production on the surface.

The Gibbs surface
Gibbs defined the equimolar surface as "a geometrical plane, going through points in the interfacial region, similarly situated with respect to conditions of adjacent matter" [39]. Many different positions can be chosen for a plane of this type. If the density of a quantity "a", ρ a , varies in x-direction according to ρ a (x, y, z), the excess density ρ s a is ρ s a (y, z) = Here, x 1 and x 2 are positions in the bulk phases and is the position of the equimolar surface. The superscripts "d" and "u" indicate a function used to extrapolate ρ a (x, y, z) from the bulk values on the left (downstream) and right (upstream) side, respectively, of the wave as illustrated by the dash-dot-line in Figure 1. The figure also illustrates the positions x 1 and x 2 . The Heaviside step function, Θ, is by definition unity when the argument is positive and zero when the argument is negative. Note that whereas the bulk density is per unit volume, the excess density is per unit area. The excess density is in general a function of the position (y, z) along the surface. We shall, however, only consider the case of constant properties in the y, z plane, and these coordinates will be omitted from here on. Moreover, the cross-sectional area perpendicular to the x-direction is independent of x. Examples of ρ a considered below are the mass, momentum, energy, and entropy densities.
All surface excess properties can be given by integrals like Eq. (10). We shall first consider the mass density ρ. Requiring the excess molar density to be zero defines the equimolar surface, which we shall use to define the position of the surface. 1 It is chosen such that the surplus of matter on one side of the surface is equal to the deficiency on the other side. The shockwave position is a function of time t , (t ). The velocity of the surface is given by in the stationary frame of reference. If v s (t ) is independent of t , the surface moves at steady state.
Like the mass density, other excess variables are given per unit area of the surface. They describe the surface and how it differs from the adjacent homogeneous phases. In particular, the surface excess mass density ρ s of the equimolar surface is zero. With the surface location so defined, other surface excess variables will in general be non-zero. Within reasonable limits, one may shift the positions x 1 and x 2 without changing the extrapolated values of interest. In this sense, the precise locations of x 1 and x 2 are not important for the value of the excess property as long as they are in the bulk phases near the surface as illustrated in Figure 2

Balance equations in the surface frame of reference
The density, flux, and production of a surface property "a", ρ a , J a , and σ a , respectively, obey dynamic balance equations. At the macroscopic scale, we may write ρ a in the form [34] The limit of the Gaussian-type distribution function shown as the dashed line in Figure 1 is the Dirac delta function in the surface description. Fluxes J a and source terms σ a are given by similar expressions. By substituting Eq. (12) into Eq. (7) and using Eq. (11) we obtain [37] , and ρ d a ( , t ) in the third line are the extrapolated values (to the surface position) of the respective quantities. In order for Eq. (13) to be correct, the sum of all terms inside each of the square brackets have to be zero. The first two brackets give equations for the bulk phases in the macroscopic description. For the surface we obtain from the third bracket where we have used the notation for the difference across the surface. Equation (14) shows that the accumulation of the property "a" on the surface is due to the difference in the flux (in the surface frame of reference) in and out of the surface plus the excess production. In particular, if we consider the mass density ρ, we find is the mass flux in the surface frame of reference. By construction, ρ s (t ) = 0. This gives the mass conservation, Eq. (1), in the surface description: The balance equation for the surface excess entropy production follows from Eq. (14): Here ρ s s is the surface excess entropy density (per unit surface area). Eq. (19) is given in terms of entropy fluxes into and out of the surface (in the surface frame of reference), and the excess entropy production, σ s s . Although Eq. (7) allows us to determine the entropy production in the entire system, we focus here on the surface for which x and t are related through the temporal position of the surface. Eq. (19) is basis for the SBM, one of the two surface methods we will use to determine σ s s quantitatively in Section 6. The fourth term in Eq. (13) gives which implies that the excess flux of the quantity "a" in the direction of the shock-wave propagation is equal to zero in the frame of reference that moves with the shock.

Conservation laws
We can now apply the general considerations in Section 3.2 to the conservation of mass, momentum, and energy. From the conservation of mass it follows that in the bulk phases. Eq. (18) describes conservation of mass for the surface. From conservation of momentum, it follows that in the bulk phases and d (ρv) s d t for the surface. In Eqs. (22) and (23), P xx = p +Π xx where p is the thermodynamic pressure and Π xx is the xx-component of the viscous pressure tensor. Making the Navier-Stokes assump- where the energy density ρ e is the sum of internal and kinetic energy density: ρ e = ρ u + ρ k , where ρ k = ρv 2 /2, and J q is the total heat flux in the barycentric frame of reference, all in the bulk phases. In the one-component system that we consider, J q = J q where J q is the measurable heat flux, which is independent of the frame of reference. Furthermore, for the surface. By analogy to the bulk energy density, ρ s e = ρ s u + ρ s k where ρ s k = (ρv 2 ) s /2 for the surface. The sum of properties in the bracket, ρ e (v −v s )+P xx v + J q , extrapolated to the surface, is the energy flux in the surface frame of reference. 2 For the excess internal energy density it follows from Eqs. (23) and (25) where the specific enthalpy is h = u + p.

The entropy production
So far, we have three routes to the surface excess entropy production, the BBM and the LIT (the integrals of Eq. (7) and Eq. (9), respectively, over the surface thickness), and the SBM, Eq. (19). We now proceed to find a third route using the Gibbs equation for the surface. We shall see that σ s s can to a very good approximation be written as the product of a mass flux and the entropy difference across the surface. Furthermore, the GEM provides detailed information about the energy conversions in the shock wave.
The integrated form of the Gibbs equation for a surface is [34], where T s is the surface temperature, 4 γ is defined by γ = (∂U s /∂Ω) {S s ,N s } where U s , S s and N s are the surface excess internal energy, entropy, and number of particles, respectively, 5 and µ s is the specific Gibbs energy of the surface. When Eq. (27) is combined with the Gibbs-Duhem equation, we find The second equality is due to the fact that the surface excess density ρ s is zero by construction. Eq. (29) is the statement of local equilibrium in the surface description. The statement implies that the surface excess properties are related in a way that can be used to assess the entropy production. The time derivative of the excess entropy density is: where ρ s s and ρ s u are determined from Eq. (10) using from the equimolar surface. Eq. (30) then gives the surface temperature T s .
By introducing Eq. (26) into Eq. (30), and comparing the result with the entropy balance, Eq. (19), we obtain the following expression for the excess entropy production, using the same bracket notation as in Eq. (15): where and where s is the specific entropy. Eqs. (32) and (33) contain quantities that are available from the equation of state plus the thermal conductivity and the viscosity. The results from this method will be compared with results from the other three methods in Section 6. The excess entropy production is independent of the frame of reference, but as a property of the surface, it will in general depend on how the surface is defined. It is in other words invariant under a coordinate transformation. We may therefore convert all fluxes and conjugate forces from any frame of reference, to the surface frame of reference and back, without changing the entropy production in the different phases, σ d s , σ s s , and σ u s .

Stationary shock front
If the shock front moves with a constant velocity, all shock-front variables, except for the position of the shock front, are independent of time. The conservation equations (23) and (26) then reduce to the Rankine-Hugoniot conditions, which in our notation are given by Eq. (35) can be used to eliminate T s in Eq. (33). Since the upstream system is at equilibrium, the upstream heat flux equals zero. The downstream heat flux is close to zero because the temperature gradient just behind the wave front is small (cf Figure 6). Under these conditions we can therefore neglect the contribution to the entropy production from the heat flux (Eq.

Constitutive equations
Eq. (36) gives a particularly simple flux-force relation in the surface description with just one flux ( j ), one force ([s] − ), and one transport coefficient (L): For the sake of completeness, we also include here the corresponding equations from the entropy production, Eq. (9). Note that the two terms on the right hand side of Eq. (9) are of different tensorial order and therefore do not couple, so that the constitutive equations are in this description: where the last terms on the right-hand side of both equations are the Fourier-Navier-Stokes constitutive equations. Eqs. (38) and (39) are valid for both steady and transient states.

Non-equilibrium molecular dynamics simulations of a blast wave
NEMD simulations were carried out with a Lennard-Jones/spline (LJ/s) model using an inhouse Fortran code. The model is defined by the pair potential where σ and ε are the usual Lennard-Jones potential parameters and a and b are coefficients in the spline function that truncates the potential smoothly between the potential's inflection point at r s and zero value at r c . The parameters a, b, and r c are determined such that the potential and its derivative are continuous at r s and r c . The LJ/s model has essentially the same features as the LJ model, but since the potenial is of shorter range, the thermodynamic properties are different. The shorter range of the LJ/s also leads to significantly shorter simulation times. Further details on the spline model and its thermodynamic properties can be found in refs. [40] and [41].
The system layout is shown in Figure 3. The simulations were made with a single component with N = 524, 288 particles in an elongated MD cell. Periodic boundary conditions were used in all three directions. The aspect ratio was set to L x /L y = L x /L z = 512 in order to give the shock wave enough distance in x-direction to separate the wave front from the heat diffusion from the blast. The number of layers was chosen so as to satisfy three criteria: (1) each layer should contain of the order 1,000 particles to ensure good signal-to-noise ratio for the properties computed in each layer, (2) the layer thickness, ∆x, should be at least of the order one Figure 3: MD cell layout. The bright orange regions at the ends is where the energy was added during the heat pulse. The aspect ratio was different from that shown in the illustration. molecular mean free path, and (3) the resolution in x-direction should be good enough to see details of the density-, pressure-, and temperature profiles in the shock front, i.e. ∆x should be at least 3-4 times smaller than the thickness of the shock front. The overall number density was set to n * = N σ 3 /V = 0.01, where V is the volume of the MD cell. 6 This low density allows us to use the virial expansion [41] as an accurate equation of state in the analysis of the shock-wave data. Condition (2) and (3) are counteractive in the sense that (2) favors a large ∆x whereas (3) favors a small ∆x. The thickness of the shock wave depends on its speed (the Mach number). An estimate is 5−10 times the molecular mean free path for Mach numbers ≈ 2 and smaller for higher Mach numbers [15]. An estimate of the mean free path based on elementary kinetic theory is λ ≈ V 2N πσ 2 , which at the actual density amounts to approximately 20 in Lennard-Jones units. The system was accordingly divided into 512 layers of equal thickness normal to the x− direction, so that each layer contained on average of the order 1024 particles. The layers were used as control volume for computing local properties of the system. With the density used in the simulations, this gives ∆x * = 29.5, i.e. approximately 50% larger than the mean free path in the equilibrium gas ahead of the wave.
The blast was generated by thermostating one layer at each end of the MD cell to a temperature T H by simple velocity rescaling [42]. The other 510 layers were not thermostated.
The simulations included 20 parallel runs. Each run was started from a configuration that was randomized with a Monte Carlo sequence of m steps, m = [1, 2, ...19, 20] × 10 5 followed by equilibrium simulations at T * = k B T /ε = 1.0. This temperature is slightly above the critical temperature for this model (T * c = 0.885, [41]) and the gas has a Prandtl number P r ≈ 1.4. The number density and the mass density are numerically identical in reduced LJ units. Each timestep was δt * = 0.002, with t * = t σ ε m 1/2 . The density and temperature used in this work correspond to argon at approximately 120 K and 4 bar (assuming the usual Lennard-Jones parameters for ε and σ, i.e. ε/k B = 124 K and σ = 3.418 Å).
Starting from the equilibrium state for each of the 20 equilibrated systems, energy was added as a pulse by setting the thermostats in the regions marked "H" in Fig. 3 for 2,000 time steps. Two cases were studied with T * H = 130 and 2080. This thermal blast and sudden increase in the local temperature and pressure at the ends of the MD cell generated pressure waves traveling in x-direction from the ends of the MD cell towards its center. After the initial 2,000 time steps, the simulation was continued as a NV E -simulation with the same time step δt * . Since the 20 equilibrium configurations were slightly different from each other, the energy inputs were also different, leading to slightly different propagations of the individual waves. This was particularly visible for the strong shock. The symmetry of the system was used to pool data from the two halves of the cell. The outcome of the 20 runs in each series was used for postprocessing of average properties and uncertainties.
A typical density profile for the left half-cell is shown in Figure 2 in reduced LJ units. The location of of the shock-wave front at a given time was determined from Eq. (10) using the excess mass density ρ s as illustrated in the insert in Figure 2. A linear function was fitted to the density profile ρ(x) for x < . For x > , the equilibrium density ρ * = 0.01 was used in the extrapolation. The condition ρ s ( ) = 0, which determines the location of the shock front, was solved by the "solve" function in MS Excel and Simpson-integration of the NEMD data.
The position of the shock front was recorded as function of time and the wave speed was computed from Eq. (11). Figure 4 shows the position as function of time for the strong shock generated in this work in comparison with the weak shock discussed in ref. [35]. Figure 5 shows that the speed decays with time, slowly for the weak shock and faster for the strong shock. This indicates that energy is dissipated faster in the strong shock than in the weak shock.
The speed of sound in the gas ahead of the wave was determined from by independent MD simulations of C p , C V , and ∂p/∂ρ T and found to be 1.298, which is essentially the ideal-gas value, 1.291, at T * = 1.0.
The NET-analysis of the entropy production requires information about the enthalpy, entropy, density, and kinetic energy in front of, and behind the shock wave. In addition, we also need the transport properties mass flux, measurable heat flux, and the x-component of the viscous pressure tensor, which includes the shear and bulk viscosities of the gas. These properties were computed as time-and spacial averages of NEMD and equilibrium MD results using the expressions shown below. The entropy was computed from the equation of state as explained in the supplementary material. We have used the kinetic temperature as a measure of the temperature in our analysis, i .e.
where k B is Boltzmann's constant, v i is the 3-dimensional velocity of particle i (all the particles have the same mass, m, in this one-component case), and v is the streaming velocity (the velocity of the local center of mass). The summation is done over all the N CV particles in the local control volume (CV), i .e. each layer in the MD cell. The local streaming velocity was determined as v = 1 where M CV = mN CV is the total mass in CV. Because the transport is in x-direction only, the y− and z− components of v are zero on average and the x−component is the local streaming velocity v, cf. Eqs. (1) -(3). The kinetic temperature in a shock wave front has different values in x-, y-, and z-directions, and is therefore a tensorial quantity [23]. The temperature was first computed in the MD frame of reference, The conversion to the kinetic temperature was done in the postprocessing using T y y = T MD y y All these quantities are local in the CV. The reported data for these quantities are space-and time averages of these quantities. The shock wave creates a sharp density gradient in the fluid. The pressure was therefore calculated using the coarse-grained version of the virial equation [43,44] with the Irving-Kirkwood contour C i j , the straight line between i and j [45]. We summarize the method here for a plane surface normal to the x-direction. Consider a pair of particles i j . One of them or both may be either inside or outside CV. The configurational contribution to the q qcomponent of the pressure in CV from that pair is where R is some point in space and l is a point on the contour C i j . In the present context, Eq. (49) reduces to where H (x i , x j ) is a book-keeping function that defines how much of the contour C i j that is inside the control volume. Further details of the algorithm were described by Ikeshoji et al. [44]. The kinetic contribution to the pressure was computed as In general, the fluxes depend on how we choose the frame of reference. In this context, there are three obvious choices, the MD cell coordinate system, the barycentric coordinate system, and the shockwave co-moving coordinate system. In MD simulations, fluxes are most conveniently computed in the MD cell (stationary) frame of reference. The mass flux j defined by Eq. (17) refers to the co-moving coordinate system and the total heat flux in Eq. (24) refers to the barycentric frame of reference. Conversion between different frames of reference was done in postprocessing as described in the following.
The local streaming velocity is given by Eq. (43). The total heat flux in x-direction in the barycentric frame of reference was given by Evans and Morriss [46]: where φ i is the potential energy of particle i in the field of all the other particles within range (including those outside the CV), f i j is the force acting on i due to j , and x i j = x j − x i is the distance from i to j in x-direction. The total heat flux in the barycentric frame of reference is equal to the measurable heat flux in the one-component system considered here. The corresponding energy flux in the MD cell frame of reference is found by setting v = 0 and v = {0, 0, 0} in Eq. (52): Eq. (53) introduced into Eq. (52) allows a separation of the heat flux into J MD q,x and the rest: where

Results and discussion
In this Section, we first discuss our findings for the kinetic properties, viz. the kinetic temperature and the velocity distributions. We show that the temperature is non-isotropic, in agreement with previous results [29]. In Section 6.2, we include the potential energy and the configurational contribution to the pressure and show that the non-equilibrium properties deviate from the equilibrium values in the microscopic description of the shock front. Section 6.3 is devoted to the entropy production computed by BBM and SBM. The GEM is a major contribution in this work and will be discussed in detail in Section 6.3.3, including an analysis of the energy conversions in the shock front. Finally, the four methods are compared in Section 6.4, where we discuss the validity of our calculation of the entropy production.

Temperature and velocity profiles
Shock waves are non-equilibrium and non-isotropic structures. For instance, the kinetic temperature in the shock front is non-isotropic as shown in previous simulations [5,24,29,40,47,48]. The insert of Fig. 6 agrees with these earlier simulations; the kinetic temperature is highly non-isotropic in the front of the strong shock, which indicates lack of local 7 equilibrium in the system. A peak in T xx is known to occur for strong shocks [24,29]. In this work, we have used We showed in a recent paper that the speed distribution in a weak shock front (Mach number 2.1) was a perfect Maxwell-Boltzmann distribution and concluded that this was consistent with a state of local equilibrium [35]. A comparison of the distribution functions for the weak and strong shocks is illustrated in Fig. 7, based on the speed of N CV ∼ 30, 000 particles (total from 20 runs) that were in a control volume of thickness ∆x * , centered at positions x * = 3434 and 6204 for the weak and strong shock, respectively, and at the end of each simulation run. 8 The mean free path is λ * ≈ 2∆x * /3 ahead of the wave and λ * ≈ ∆x * /3 behind the wave. The Figure 7: Particle speed distributions for the weak and strong shocks. The data for the weak shock were recorded at x * = 3434 and t * = 1, 000, and for the strong shock at x * = 6204 and t * = 600. The weak shock shows a perfect Maxwell-Boltzmann distribution, whereas the strong shock does not. The uncertainties are three standard errors.
fitted Maxwell-Boltzmann distribution gave a temperature T * = 1.92±0.01 for the weak shock, in fair agreement with the local kinetic temperature T * = 1.79 ±0.01. The corresponding numbers for the strong shock are T * = 3.4 ± 0.2 from the fitted distribution, in poor agreement with the local kinetic temperature T * = 5.1 ± 0.7 (uncertainties given as three standard errors of the mean).
To compute the surface excess entropy production with the GEM boils down to using Eqs. (31)- (33). It is worth noting that apart from the variable T s in these equations, the values of all the other properties are extrapolated values from the regions ahead of, and behind the shock front. The system ahead of the shock is in local (and global) equilibrium. The system immediately behind the shock is also in local equilibrium as shown by the temperature profiles. For the purpose of this work, we conclude that, despite the fact that the system in each control volume in the shock front is not in local equilibrium for the strong shock, the thermodynamic properties used in the GEM (adjacent to the shock front) are in local equilibrium.

Energy and pressure profiles
In this section, we consider the configurational contributions to local properties, in particular the internal energy and pressure. The system in question in this work is a moderately dense gas, so the configurational contributions are likely to be small. We will nevertheless assess to what extent the non-equilibrium configurational properties deviate from the equilibrium values with focus on the mechanical properties internal energy and pressure, which are easily obtained in NEMD. Irrespective of the main theme of this work, viz. the entropy production, this assessment provides interesting insight in the shock wave by itself.
The specific internal energy, u, for the strong shock is shown in Figure 8. The internal energy is completely dominated by the kinetic part, the potential (configurational) energy accounts for at most ≈ 0.3% of the total internal energy. 9 The potential energy is less negative at the downstream side of the surface because the particles are on average closer together there. 9 The potential energy contribution will clearly be larger at higher densities, such as in a liquid. The dots show the total internal energy determined by NEMD, u * (total), and the squares show the configurational (potential) part of it, u * p . The black and white squares show data from the non-equilibrium and equilibrium simulations, respectively. Note that the configurational contributions (referring to the right axis) are so small that they do not visibly separate the kinetic contributions to the internal energy from the total in u * . The error bars are three standard errors, the errors for the equilibrium results are smaller than the symbol size. The vertical dashed line shows the position of the Gibbs equimolar surface.
The potential energy does show a difference between the non-equilibrium energy and the energy determined by equilibrium simulations in the range 6, 200 ≤ x * ≤ 6, 300. The difference is ≈ 40% of the potential energy or ≈ 0.1% of the total internal energy. The equilibrium data were generated at the local non-equilibrium density and temperature T = 1 3 Tr(T) with the temperature components given by Eqs. (45) - (47). It is interesting to note that this difference occurs at the upstream, low-density side of the equimolar surface. Based on this result, we conclude that the internal energy density is very accurately given by the equilibrium values. Figure 9a shows the total pressure 1 3 Tr(P) as function of x in the shock front region for the strong shock at t * = 600. The non-equilibrium pressure is, assuming the Navier-Stokes relations, 1 3 Tr(P) + η B ∂v x ∂x ≈ 1 3 Tr(P). The bulk viscosity η B is small for a monatomic dilute gas. We have estimated, based on data from Hoheisel et al. for the Lennard-Jones fluid [49], that η * B < 10 −3 for our Lennard-Jones spline system in the actual states, which makes the contribution from the bulk viscosity to the pressure at least four orders of magnitude smaller than 1 3 Tr(P). We have therefore assumed that η B = 0 in this work. The total pressure may be split into a kinetic (ideal-gas) contribution, 1 3 Tr(P kin ) and a configurational term, 1 3 Tr(P conf ). For the kinetic term, we have used (in Lennard-Jones units) 1 3 Tr(P * kin ) = ρ * 1 3 Tr(T * ) = ρ * T * = p * ig . The configurational term was computed according to Eqs. (48)- (50). The pressure is almost zero ahead of the shock and increases monotonically through the front. The ideal-gas pressure accounts for 98 % and 96 % of the total pressure immediately upstream and downstream, respectively, of the shock. Figure 9a also shows a comparison between the non-equilibrium and equilibrium configurational pressures. The configurational pressure is slightly negative ahead of the shock wave and positive behind the wave where the gas is more compressed. Unlike the configurational energy, there is virtually no difference between the equilibrium and non-equilibrium pressures. This is consistent with the Navier-Stokes relations, Eqs. (4) and (5), which imply that 1 3  6154. (a) The black dots and triangles represent the total non-equilibrium pressure, 1 3 Tr(P), and its kinetic contributions, p ig , respectively. The squares are the corresponding equilibrium pressures p and p conf . Note that the kinetic contributions are equal in the two cases. (b) The x-component of the pressure tensor (P xx ) and the viscous pressure computed from the NEMD data (Π MD xx ) and Navier-Stokes (Π NS xx ).
Eqs. (4) and (5) may be combined to Π MD xx ≈ p for the viscous pressure, where superscript "MD" means "as determined from the NEMD simulations". The viscous pressure may also be determined as ∂x where we have used superscript "NS" to distinguish it from the NEMD results. The shear viscosity η S was determined by independent nonequilibrium MD simulations with LAMMPS [50] and ∂v x ∂x was taken from the velocity profile in the shockwave. Fig. 9b shows a comparison between Π MD xx and Π NS xx . The main observation is that Π xx contributes only in the shock front where ∂v x ∂x is significant. Here, Π xx accounts for up to 40% of P xx . The dominant viscous contribution is on the low-density side of the equimolar surface. The agreement between Π MD xx and Π NS xx is within the uncertainties, and again consis-tent with the Navier-Stokes equations. This indicates that the Navier-Stokes equation gives a correct description of the pressure profiles.

The entropy production
Based on the analyses in sections 6.1 and 6.2, we found that the assumption of local equilibrium is good as measured by the internal energy and pressure. Lacking values for the nonequilibrium entropy, we shall in the following assume that also the entropy can be estimated by the equilibrium values. We shall now use these results to determine the surface excess entropy production with the four methods; the BBM based on Eq. (7), the LIT based on Eq. (9), the SBM based on Eq. (19), and the GEM based on Eqs. (31) -(33).

The bulk balance method (BBM) and linear irreversible thermodynamics (LIT)
Integrating Eq. (7) over the thickness of the wave, we get the total entropy production in the wave, which is also the surface excess entropy production: In this method, we consider the properties of the wave as continuously changing over the wave front like in the bulk. The terms ∂ ∂t ρ s (x, t ) and ∂ ∂x J s (x, t ) in Eq. (7) were determined by fivepoint numerical differentiation with the results shown in Figure 10a. The two contributions are opposite in sign with a relatively small sum. The integrand σ s (x, t ) and the integration limits x 1 and x 2 are shown in Figure 10b. A Gaussian function was fitted to σ s (x, t ) to smooth the NEMD data and the fit was integrated analytically. The graph shows that the values of the integration limits were not critical, the entropy production occurs only in the shock front. The negative dip in σ s (x, t ) at the left side of the peak in Figure 10b is due to a slight mismatch in the peaks of ∂ ∂t ρ s (x, t ) and ∂ ∂x J s (x, t ). We believe this is not significant and an artifact of the five-point numerical differentiation methods we have used. The results from the BBM show that the entropy production is negligible a few mean free paths away from the shock front. The surface represents the dominant entropy production. This procedure was repeated for a series of times between t * = 200 and t * = 600. The surface excess entropy production is shown as function of time in Figure 15 and compared with data from the other three methods used.
We showed in Sections 6.1 and 6.2 that although the system is not in local equilibrium in the shock front region, it is close to being so. This is the basis for using the Gibbs equation in the normal way for bulk fluids [36], like Velasco and Uribe did [30]. The key result for the entropy production in this method is Eq. (9), which integrated over the shock thickness (Eq. (56)) gives σ s s . The σ s (x) determined in this way is shown in Figure 10b marked "LIT". The agreement with the BBM is very good, indicating that the local equilibrium assumption is good even for the strong shock. Whereas σ s (x) determined from Eq. (7) is a small difference between large numbers (cf. Figure 10a), when determined from Eq. (9), it is a sum of small numbers and therefore less noisy, especially downstream of the shock front. All the quantities on the RHS of Eq. (9) were determined directly from the NEMD results without any assumptions for the transport coefficients. Integrating Eq. (9) over the entire system showed that some 97 % of the total entropy production occurred in the shock front.

The surface balance method (SBM)
In the derivation of Eq. (19), we considered the wave front as a surface, but without employing the Gibbs equation. The term ρ s s (t ) was determined from Eq. (10) with the result shown as function of time in Figure 11. As input to Eq. (10), we used the entropy density given by the equation of state, using the local density and temperature as input. The equation of state we used was based on the virial expansion and is given in the supplementary information. The The integrand σ s (x, t ) and the fitted Gaussian. We also show results from the LIT method used by Velasco and Uribe [30], i.e. the Gibbs equation used locally for each control volume in the simulation (see text). time derivative was determined from a linear fit to ρ s s (t ). The surface excess entropy density varies little with time and contributes less than 2 % to σ s s . The term [J s ] − depends on the extrapolated values from properties outside the shock front, which we have shown in Sections 6.1 and 6.2 to be well represented by equilibrium values. Because of this, we consider Eq. (19) to give a reliable estimate for the surface excess entropy production. Results are compared with the other three methods in Figure 15.

The Gibbs excess method (GEM)
In the GEM, the surface excess entropy production is determined from Eqs. (31) - (33). All quantities in these equations, except the surface temperature T s , are determined by extrapolating properties from the bulk phases adjacent to the wave front. We established in Sections 6.1 and 6.2 that these properties are given by their equilibrium values in the present case. The surface temperature is, however, given by Eq. (30). This equation involves the excess properties ρ s u and ρ s s determined from the entire profiles, including the non-equilibrium properties in the wave front, with the use of Eq. (10). We will return to the question of how σ s s is affected by the uncertainty in T s in the following subsection and for the moment use T s as determined from the available data. Figure 12 shows a plot of ρ s u vs. ρ s s . The relationship is linear with slope T s * = 3.9, which is between the upstream and downstream temperatures.
A plot of the local values of σ q and σ j (Eqs. (32) and (33)) is shown in Fig. 13 as function of x * for t * = 600. The surface excess entropy production is the difference across the shock front as indicated by the double arrow in the figure. The σ q is very small except in the shock front, where ∂T /∂x is large. Ahead of the shock, σ q is zero because the system is at equilibrium there with J q = 0. Behind the shock, J q ≈ 0, because ∂T /∂x is small there, cf. Fig. 6. The difference [σ q ] − of the extrapolated values is therefore practically zero. The blip in the front region is due to the fact that T d < T s < T u , but this is of no importance to the value of [σ q ] − . By comparison, σ j is everywhere large and [σ j ] − is significant. Unlike J q , the mass flux depends on the frame of reference, and j in Eq. (33) is given with the surface as reference. So is also the kinetic energy term in the parenthesis in Eq. (33). The local values of σ j must therefore not be confused with the local entropy production, it is the difference across the wave, indicated as the double arrow in Figure 13, which is relevant for σ s . The individual terms in the parenthesis of Eq. (33) are shown in Fig. 14 for the strong shock. The viscous pressure term varies little over the shock front and the difference between the extrapolated values is practically zero. The kinetic energy term includes the center-of-mass velocity relative to the shock wave velocity. This relative velocity is larger upstream than downstream, so the difference defined by the bracket is positive. Both h and s increase when the shock wave passes. The viscous term is negligible. The mass flux is constant across the shock front because mass is conserved, and therefore equal to the upstream value, j = −ρv s . In total, the term [σ j ] − is positive. Hence, for the propagating shock examined in this work, the overall picture is that kinetic energy is converted to enthalpy. A minor amount of the wave's energy produces entropy across the shock front, leading to a slow retardation of the wave.

How sensitive are the results to uncertainties in the estimated surface temperature?
In the GEM, the surface excess entropy production is determined according to Eqs. (31) - (33). All quantities in these equations, except the surface temperature T s , are determined by extrapolating properties from the bulk phases adjacent to the wave front. We established in Sections 6.1 and 6.2 that these properties are given by their equilibrium values in the present case. The surface temperature is, however, given by Eq. (30). This equation involves the excess properties ρ s u and ρ s s determined from the entire profiles, including the non-equilibrium properties in the wave front, with the use of Eq. (10). We will therefore now estimate how much the uncertainty introduced by using the equilibrium entropy instead of the non-equilibrium entropy in Eq. (10) affects the value of the surface excess entropy production.
For shock waves near steady state and J d q ≈ 0, Eq. (36) holds. In this case, σ s s is independent of T s . Alternatively, we may consider the dominant contribution to σ s s , σ j − , and take the derivative with respect to T s : The results in Section 6.3.3 show that the bracket on the right hand side of Eq. (57) is small, which means that σ j − is rather insensitive to errors in T s . As an example, inserting numerical values for t * = 600 shows that [σ j ] − changes by ±0.1% for a ±10% change in T s . Our conclusion is that the surface excess entropy production is very insensitive to the value of T s and that the Gibbs excess method is an accurate method in the present case.

Entropy production and blast wave decay
A comparison between the surface excess entropy production computed from the four methods employed in this work is shown in Figure 15. The four methods are consistent for both Mach numbers, which adds confidence to the assumptions made. The entropy production decreases with time as expected as the wave moves away from the blast, looses energy and slows down, cf. Figure 5.
The four methods differ in the ways the sources are computed, but they give the same entropy production. The time derivative of the entropy density and the space derivative of the entropy flux used in the BBM are both large in the front region and of opposite sign (cf Figure 10a), and the local entropy production is a small difference between relatively large numbers. The BBM is therefore sensitive to errors in these quantities. The SBM and the GEM depend on the time derivative of the excess surface entropy density, which varies little with time. The main contributions in these methods are jumps in extrapolated quantities determined from bulk properties adjacent to the surface, which are robustly determined from the equilibrium equation of state. The excess entropy production depends strongly on the Mach number, with approximately a factor 10 increase in the produced entropy when the Mach number increases from 2 to 6, or approximately the square of the ratio between the Mach numbers, M strong /M weak 2 . This difference is also reflected in the retardation of the shock wave, the strong shock slows down much faster than the weak shock. A fitted exponential function, σ s s = σ 0 exp(−αt ), to the values for σ s s gave the parameters σ 0 = 0.034 and α = 0.0018 for the weak shock and σ 0 = 0.38 and α = 0.0015 for the strong shock. This means that the relative decay is approximately the same for the two Mach numbers, but the intensities differ by a factor of approximately ten.
We have also included results from the assumption used by Velasco and Uribe, viz. that the Gibbs equation is locally valid in each control volume (marked "LIT" in Figure 15) [30]. These results are systematically lower the other three, albeit not by very much. The difference may be an indication that although the local equilibrium assumption is good, it may lead to systematic errors in the computed entropy production.
A key property in analyses of shock waves is the peak overpressure, i.e. the maximum pressure in the shock wave minus the ambient pressure in front of the shock [11]. The peak overpressure, ∆P, is given by Jones [51] where p 0 is the ambient pressure ahead of the shock, C p and C v are the heat capacities at constant pressure and volume, respectively, and M is the Mach number. The results from Eq. (58) are compared with NEMD results in Figure 16. For this comparison, the ambient pressure and the Mach number were taken from the NEMD simulations and the heat capacities were determined by separate MD simulations. The agreement is good for both the weak and the strong shock. Jones also gave a relation between the overpressure and the blast energy, which for a plane wave reads in our notation: where R 0 is a characteristic distance related to the blast energy [51]. The exponent m is equal to 1 in the limit of M → ∞ and 1/2 for M → 1 for a plane wave. We found m = 0.87 for the strong wave and m = 0.69 for the weak wave. Using our results together with the limiting values, we found the empirical relation m ≈ 1−0.82(1/M )+0.32(1/M ) 2 . Finally, we note that the distance R 0 is proportional to the blast energy [51]. Using R 0 as an adjustable parameter when fitting Eq. (59) to our NEMD data, we found the ratio (R 0 ) strong /(R 0 ) weak = 21.9, in excellent agreement with the ratio between the blast energies in the two cases, which was 22.0.

Conclusions
In this work, we have applied non-equilibrium thermodynamics for surfaces [34] and analyzed the entropy production in two shock waves using four different methods. We have developed a new tool, the "Gibbs excess method" (GEM) and compared it with three other methods. In the "bulk balance method" (BBM), the entropy balance was integrated over the thickness of the shock wave. The LIT method is based on the assumption of local equilibrium in the shock wave front, the local version of the Gibbs equation, and integration of the local entropy production over the shock thickness [30]. In the "surface balance method" (SBM), we used the concept of Gibbs equimolar surface combined with the entropy balance equation across the surface. In the GEM, we took the SBM one step further by using the Gibbs equation for surfaces and derived new expressions for the surface excess entropy production, performed a detailed analysis of the energy conversions in the shock wave front, and found a very simple approximate expression for the entropy production. The most significant difference between the four methods is that BBM and LIT assume local equilibrium everywhere in the fluid, including the shock front, whereas the SBM and the GEM use surface properties in equilibrium. The SBM and GEM are therefore more robust and may be easier to apply. Two plane blast waves were simulated with non-equilibrium molecular dynamics (NEMD) in a Lennard-Jones/spline system with 524,288 particles. Prior to the blast, the system was equilibrated at T * = 1.0 and ρ * = 0.01. The two shocks propagated at almost steady state with Mach numbers approximately 2 and 6. We found the typical difference in the x-and ycomponents of the kinetic temperature, but based on analyses of the particle speeds, potential energy, and pressure, we concluded that the conditions for using non-equilibrium thermodynamics were well satisfied.
The four methods were based on different approximations and used the NEMD data in different ways, but the surface excess entropy productions were in excellent agreement. We found a small deviation from local equilibrium in the front region of the strong shock, but this is of no importance in the GEM, which uses extrapolated equilibrium data from the adjacent bulk regions. From this observation and verifications of some of the assumptions, we conclude that the results are reliable. For the GEM, we found that the differences across the surface in the measurable heat flux and the viscous pressure gave negligible contributions to the entropy production. This is in contrast to the LIT method, in which these are the only two sources to the entropy production. The GEM provides new, detailed information about energy conversions in shock waves. In short, most of the wave's kinetic energy is converted reversibly to enthalpy. A smaller fraction of the waves total energy was dissipated, which led to a weak retardation of the wave.
In principle, the BBM makes no assumption of local equilibrium, but lacking data for the non-equilibrium entropy in the front region, we had to use equilibrium data here. The SBM may be more robust that the BBM because the surface excess entropy density is almost constant with time and its time derivative is almost zero. The LIT gives results with little statistical noise, but may suffer from the lack of local equilibrium in the front region.
As the waves were almost at steady states, the Rankine-Hugoniot (RH) conditions were found to describe the waves well. The peak overpressure determined from the RH conditions agreed very well with the NEMD data. The shock-wave thickness was found to agree with theoretical estimates and experimental data [9] and simulations [29].
The combination of different theories and NEMD data presented here gives new tools to study shock waves. In particular, the GEM is a robust method that relies on equilibrium data adjacent to the shock wave front. For waves close to steady state, a good appoximative value for the surface excess entropy production can be found in a very simple way as given by Eq. (36), viz. as the mass flux in the surface frame of reference times the difference in specific entropy across the surface. and Data Storage in Norway and by Department of Chemistry at The Norwegian University of Science and Technology -NTNU. We are grateful to Olav Galteland for providing data for the shear viscosity using LAMMPS with the SLLOD algorithm. 10

A Equation of state for the Lennard-Jones/spline gas
We consider here a one-component Lennard-Jones/spline (LJ/s) fluid at low density with the purpose to find an expression for its entropy and internal energy. The entropy, S, can be derived from the Helmholtz energy, A, as where T is temperature and V is volume. The Helmholtz energy can be found by integrating where n = N /V is the number density of particles and k B is Boltzmann's constant. Using the polynomial fit in inverse temperature given in ref [41], the first three virial coefficients are with the coefficients given in Table A Note that the corresponding densities are obtained as ρ s = ns and ρ u = nu.