Jean-Bernard Maillet
email: jean-bernard.maillet@cea.fr

Emeric Bourasseau

Laurent Soulard

Jean Clerouin

Gabriel Stoltz

Jean Clérouin

Constant entropy sampling and release waves of shock compressions

Keywords: numbers: 65.20.De, 62.50.Ef, 02.70.Ns

I. INTRODUCTION

The exploration of the thermodynamic behavior of materials under extreme conditions usually follows two paths corresponding to two experimental devices: isothermal compressions and shock compressions. Isothermal compressions are performed with diamond anvil cell (DAC) techniques, and are used to compress materials up to very high pressures, although with limited temperatures. On the other hand, shock compression experiments investigate the high-pressure/high temperature regions through the propagation of dynamic shock waves in the system. Nevertheless, the thermodynamic domain available using shock experiments remains limited to the so-called Hugoniot curve, which is, by definition, the collection of thermodynamic states which can be reached from a system at fixed initial conditions, with shocks of increasing strengths. Such a constraint arises from the fact that the system is assumed to satisfy Euler's fluid equations (i.e inviscid Navier-Stokes equations), and physically meaningful shocks therefore have to fulfill the Rankine-Hugoniot conditions, which relate the thermodynamic parameters of the fluid at rest, and the thermodynamic parameters of the shocked material. Another constraint is that shock waves are adiabatic, therefore leading to very large temperature increases in the material, which limits its compressibility. The equations of state (EOS) used to predict the material's behavior at the extreme conditions encountered are often simple extrapolations of EOS fitted on available data, i.e. shock data and DAC data. It then appeared interesting to enlarge the experimental domain of investigation of materials behavior using dynamic compression set-ups, and particularly isentropic compressions. Several experimental set-ups allow to load a pressure ramp in a material. The first one is the high pulsed power (of which the sandia Z machine and the High Explosive Pulsed Power 1,2 at LANL are good examples). The second one consists in using an impactor with a varying density along one direction, as proposed initially at the AIP -SWCM conference. [START_REF] Lyzenga | proceedings of the AIP conference Shock Compression and Condensed Matter[END_REF] A successful technique is to stack slices of different materials, leading to the so-called PILLOW impactors at Sandia, [START_REF] Barker | proceedings of the AIP conference Shock Compression and Condensed Matter[END_REF] MIVAR impactors in France, [START_REF] Perez | [END_REF] and more recently the FGM (Functionaly Graded Materials) impactors at LLNL, 6 allowing a real design of a thermodynamic path as a succession of shock and release waves. The last one concerns experiments of Barnes' type where the compression is the consequence of the isentropic release of another material, as for example detonation products. 7,8 Experiments involving isentropic compression are of great interest to reach high compres-sion states, or in geophysic applications to reach states representative of the earth's core.

Experiments involving a precise evaluation of release waves in materials need also general numerical methods to compute the states reached by isentropic pocesses.

Up to now, simulations involving shock processes are rather well developed due to the simplicity of the Hugoniot equations. Those studies are performed within the framework of statistical physics, see Refs. 9,10 for reference textbooks on molecular simulation, and Ref. 11 for reference works on nonequilibrium simulation of shock waves. Any state lying on the Hugoniot can be reached from the reference state by searching for a given compression the temperature for which the pressure and the total energy of the system satisfies the Hugoniot relation. The search can be implemented in very efficient manners. 12,13,14 Those methods are now adapted to classical molecular dynamics and Monte Carlo, as well as quantum molecular dynamics. 15 Such an easy method does not exist for isentropic processes. In this paper, we present or recall several equilibrium methods which allow to follow isentropic paths, both for classical or quantum molecular dynamics simulations. We contrast these methods in terms of their precisions, rigor and computational requirements. We compare the results obtained from equilibrium simulations with release waves observed in nonequilibrium molecular dynamics.

The comparison between equilibrium and nonequilibrium methods therefore measures how isentropic the expansion of the system is. It is expected that release waves of a perfect non-viscous fluid are isentropic. For simple monoatomic fluids such as argon, it is often assumed that the release is isentropic, and viscosity effects are neglected. Our results show that even in this simple case, the release is not strictly isentropic and some corrections have to be taken into account. As a by-product of our study, we also explore more precisely the relationship between the Hugoniot and the isentrope curves, from a numercial viewpoint, but also giving a statistical physics proof of the coincidence of the curves for small compressions (see Appendix B).

Organization of the paper

The paper is organized as follows. In Section II, we present the nonequilibrium method used to simulate rarefaction waves, while some equilibrium methods for constant entropy sampling are recalled in Section III -details of the practical implementaion of the thermo-dynamic integration and precisions on the isentropic integration are given in Appendix A and C respectively). In Section IV, a comparison of numerical results obtained in the case of release waves in argon is performed, and we discuss whether release waves are isentropic.

II. NONEQUILIBRIUM SIMULATIONS

A. Microscopic description of physical systems

We first recall how a system is described at the microscopic level using statistical mechanics. [START_REF] Balian | From Microphysics to Macrophysics[END_REF] Statistical physics is a theory which allows to compute thermodynamic (macroscopic) properties of a system knowing the microscopic interactions between its constitutive elements. Consider a microscopic system composed of N particles, confined in a simulation box

D = [0, L x ] × [0, L y ] × [0, L z ]. The volume of the domain is denoted by V = |D| = L x L y L z .
The system is characterized by the positions q = (q 1 , . . . , q N ) and momenta p = (p 1 , . . . , p N ) of the particles, which have masses m i , and interact through a potential energy function U.

The phase-space Ω is the collection of all possible microscopic configurations (q, p) of the system.

The central quantity describing the system is the Hamiltonian

H(q, p) = N i=1 p 2 i 2m i + U(q 1 , . . . , q N ), (1)
which gives the energy of a given microscopic configuration (q, p). Average thermodynamic properties of the system can be computed as averages of functions of the microscopic variables O(q, p) (the so-called observables) with respect to the canonical measure at a temperature T , for a given simulation box:

O = O V,T = Ω O(q, p) π V,T (q, p) dq dp. (2)
The canonical measure associated with the Hamiltonian (1) weights microscopic states according to their energies using a Boltzmann weight:

16 π V,T (q, p) = 1 Z V,T e -βH(q,p) , β -1 = k B T, (3)
where k B is the Boltzmann constant, and the partition function Z V,T is a normalization factor so that (3) is indeed a probability measure:

Z V,T = Ω e -βH(q,p) dq dp. (4)
We have indicated explicitly the dependence of the canonical measure (3) and the partition function (4) on the temperature T and the volume V since these will be the parameters allowed to vary in the sequel.

B. Nonequilibrium simulation of release waves

Similarly to what has been proposed for the simulation of shock waves, 11 isentropic compressions or releases can be simulated directly using Non-Equilibrium Molecular Dynamics (NEMD). A straightforward numerical set-up to this end is simply to throw a low speed piston towards the sample (creating a weak shock), and then accelerating the piston in time.

Except this external forcing, the system evolves according to the standard hamiltonian dy-

namics      qi = p i m i , ṗi = -∇ q i U(q), (5)
which is integrated in time with the Verlet scheme. 17 A linear compression ramp would be obtained in the case where the acceleration is constant in time.

To obtain release waves, a shock wave can be loaded in a sample; when this shock wave is reflected when interacting with a free surface, it transforms into a release wave, supposedly isentropic. In this study, we start the release from an equilibrated state obtained from a preliminary canonical simulation, using three dimensional periodic boundary conditions.

When the system is equilibrated, the periodic boundary conditions are removed in the x direction. Two release waves are then created at the two free surfaces, and they propagate in opposite directions towards the center of the box. This process is illustrated in Figure 1.

From the simulation data presented in Figure 1, profiles of thermodynamic quantities (average densities, (kinetic) temperatures and pressures) can be extracted and averaged over thin slices. Moreover, the two release waves being symmetric, their related profiles can be averaged. A superposition of the profiles, taken at different times but projected back in the same thermodynamic diagram, is then obtained and averaged over, leading to a single profile.

The accuracy of the computation increases with the system size: an increase in the size of the transverse directions decreases the uncertainties on the slice averages (thanks to a thermodynamic limit), and an increase in the longitudinal direction allows to accumulate more profiles in time, therefore reducing statistical errors.

III. EQUILIBRIUM METHODS FOR CONSTANT ENTROPY SAMPLING

We present in this section three methods to compute the collection of all states (in terms of their temperature, volume and pressure) which have the same entropy as some reference state. These methods therefore allow to draw a curve in the (V, T ) diagram (or in the (P, T ) or (P, V ) diagrams), called the isentrope, and will be used as benchmark methods in Section IV to check whether release waves computed by NEMD simulations are indeed isentropic or not. We emphasize that, altough presented for the computation of isentropic releases, all the methods described in this section may also be used to determine isentropic compressions. We also recall a fourth method, used to obtain the entropy of a system once the entropy of some reference state (such as the perfect gas) is fixed. The computational cost of the latter method as well as its low accuracy for dense states prevented us from applying it to enough points to obtain an entire isentrope curve, and we therefore limited its use to a consistency check on the results obtained with the other methods.

A. Thermodynamic integration

The entropy of the system varies when the simulation conditions are changed. Here, we consider that the states visited by the release wave are a succession of local equilibrium states, which can be described within the canonical ensemble as given by statistical physics.

Therefore, the state of the system is defined by two parameters, its volume (equivalently, the density) and its temperature.

1. Variables indexing the variations.

Consider a general transformation in which both the volume accessible to the system (equivalently, the density) and the temperature are varied. We restrict ourselves to variations of the domain in one spatial direction only, to model the anisotropic behavior of release waves.

Assuming that the state of the system at rest can be described by some cubic simulation box with periodic boundary conditions, the volume under compression may be indexed by a variable λ 1 , so that the associated simulation domain

D(λ 1 ) = [0, (1 + λ 1 )L x ] × [0, L] 2 has a volume V (λ 1 ) = (1 + λ 1 )L x L 2 .
Notice that we consider L x = L since we may start from a uniaxially compressed state. The temperature variations are indexed by a parameter λ 2 :

T (λ 2 ) = (1 + λ 2 δT )T,
for some reference temperature T and a given relative temperature variation δT , the temperature variation being therefore ∆T = T δT . The reference inverse temperature is still

β -1 = k B T .
The particular case where only the temperature is changed (while the volume is kept constant) corresponds to λ 1 constant, while isothermal transformations are characterized by λ 2 remaining constant. Expansions correspond to λ 1 > 0.

Parametrization of the isentrope curve.

The isentrope is the locus of the points in the (λ 1 , λ 2 ) space such that the entropy normalized by the Boltzmann factor

S k B = U -F k B T (6)
is constant, F denoting the free energy of the system, and U its energy. This thermodynamic relation can be converted into an equivalent formula in the framework of statistical physics, which is much more convenient from a computational viewpoint:

U ≡ U(T, V ) = H V,T ,
where the canonical average is defined in Eq. ( 2), and

F ≡ F (T, V ) = -k B T ln Ω e -βH(q,p) dq dp.
We start from some reference state described by the parameters (λ 1 , λ 2 ) = (0, 0). The statistical physics reformulation of the requirement that (6) be constant is then

S(λ 1 , λ 2 ) -S(0, 0) = 1 k B T (λ 2 ) H V (λ 1 ),T (λ 2 ) - 1 k B T (0) H V (0),T (0) + ln     Ω(λ 1 )
e -H(q,p)/k B T (λ 2 ) dq dp

Ω(0)
e -H(q,p)/k B T (0) dq dp

    = 0.
In this expression, the phase-space Ω(λ 1 ) is the collection of all possible microscopic configurations of the system associated with a domain D(λ 1 ) of volume V (λ 1 ).

Numerical implementation

To determine the isentrope curve, we compute the entropy variation along a given path in (λ 1 , λ 2 ) space going through the reference initial state, and search for the point such that the entropy difference with this state is 0. A simple choice is illustrated in Fig. 2. It consists in performing (i) an isothermal rarefaction, going from the initial compressed state (0, 0) to an intermediate state (λ 1 , 0) with λ 1 ≥ 0;

(ii) in a second step, an isochore cooling, going from the intermediate state (λ 1 , 0) to some final state (λ 1 , λ 2 ), resorting to a maximal temperature difference λ 2 ∆T < 0 large enough.

The idea is that, in general, the first part of the transformation increases the entropy of the system (since more space becomes available for the particles), while the entropy decreases in the second part (since the temperature decreases). Of course, more general paths, with joint variations of λ 1 and λ 2 , could be considered.

The energies H V (λ 1 ),T (λ 2 ) are computed using standard sampling strategies, while the remainder in the expression of S(λ 1 , λ 2 ) -S(0, 0), a ratio of partition functions, is estimated using standard techniques for free-energy calculations. This is detailed in Appendix A.

We emphasize that this procedure is time consuming since it requires many equilibrium samplings to obtain one point on the curve. It is however exact (up to statistical errors and discretization errors in the integrals defining A), and can be straightforwardly parallelized since the equilibrium samplings required are independent.

B. Successive hugoniostat simulations

The variations of macroscopic quantities across a shock interface are governed by the Rankine-Hugoniot relations, which relate the jumps of the quantities under investigation (pressure, density, velocities) to the velocity of the shock front. The third Rankine-Hugoniot conservation law for the Euler equation governing the hydrodynamic evolution of the fluid reads (macroscopic quantities are denoted by curly letters)

H = U -U 0 - 1 2 (P + P 0 )(V 0 -V ) = 0. ( 7
)
In this expression, U is the internal energy of the fluid, P its pressure, and V its volume.

The subscript 0 refers to the initial state (the pole), the other quantities are evaluated at a state obtained from some shock compression, after equilibration. The Hugoniot curve corresponds to all the possible states satisfying (7). In practice, the collection of these states may be computed by nonequilibrium simulations with shocks of different strengths, inducing various compressions.

Alternatively, small equilibrium simulations may be used, relying on the statistical physics reformulation of the Hugoniot relation:

H(λ 1 , λ 2 ) -H(0, 0) = H V (λ 1 ),T (λ 2 ) -H V (0),T (0) + λ 1 2 V (0) P xx V (λ 1 ),T (λ 2 ) + P xx V (0),T (0) = 0. ( 8
)
The xx component of the pressure tensor is, for a simulation domain of volume V (λ 1 ),

P xx (q, p) = 1 V (λ 1 ) N i=1 p 2 i,x m i -q i,x ∂ q i,x U(q). (9)
For a given variation of the volume for instance (indexed by λ 1 ), the variation λ 2 T δT of the temperature is sought for, using for instance the techniques described in Refs. 12,13.

The Hugoniot curve does not have a priori any relationship with the isentrope curve.

However, it can be shown that the entropy variation along the Hugoniot curve is negligible up to terms of order three in the volume variable; the Hugoniot and the isentropic curves are osculatory. We present in Appendix B two proofs, the standard proof based on thermodynamic relations, and a new proof fully relying on a statistical physics reformulation.

The good agreement between the Hugoniot and the isentrope for small compressions and/or expansions can be used to compute the isentropic curve as a succession of weak shocks or weak releases, this approximation getting more accurate as the shock compressions are weakened. The only parameter left in this method is the relative volume change δV /V = λ n+1 1 λ n 1 during the instantaneous compressions or releases. We used the Hugoniostat method 12,13 to compute a sequence of states (λ n 1 , λ n 2 ) such that H(λ n+1

1 , λ n+1 2 ) = H(λ n 1 , λ n 2 )
, the corresponding thermodynamic properties at these states being obtained as a by-product of the simulation.

C. Isentropic integration

Another way to perform thermodynamic integration along an isentropic path has been proposed by Desjarlais. [START_REF] Desjarlais | Workshop on accelerator-driven warm dense matter[END_REF] The method relies on the equilibrium evaluation of ∂P ∂U (see Eq. ( 11) below). It could be applied to a system where the pressure is not isotropic upon replacing the pressure observable by the xx component of the pressure tensor.

The total differential of the entropy can be written as:

dS = ∂S ∂T V dT + ∂S ∂V T dV. (10)
For constant entropy processes,

∂S ∂T V = 1 T ∂U ∂T V = - ∂P ∂T V , ∂S ∂V T = ∂P ∂T V ,
so that, along the isentrope,

dT T = - ∂P ∂T V ∂U ∂T V dV = - ∂P ∂U V dV.
This equation can be integrated as

T 2 T 1 = exp - V 2 V 1 ∂P ∂U V dV , (11)
giving the temperature T 2 at which the system at volume V 2 has the same entropy as the system in the reference state (T 1 , P 1 ). This formula is evaluated in practice by discretizing the integral appearing in the exponential, and approximating the integrand using standard canonical sampling procedues. We refer to Appendix C for more precisions.

D. Evaluation of the entropy based on the chemical potential

This technique, which can be used only for systems in a fluid phase, follows the classical methodology of computing the free energy F of a system starting from the thermodynamic relation: 19

F = U -T S = Nµ -PV, ( 12
)
where µ is the chemical potential, defined in the canonical ensemble as

µ = ∂F ∂N . (13)
In the case of a canonical simulation, all thermodynamic quantities are functions of the volume and the temperature, so that

S = U(T, V ) + Nµ(T, V ) -P(T, V )V T . (14)
This expression allows to compute the absolute entropy of the system provided the chemical potential is known, 20 the average pressure and energy being computed using standard sampling techniques. The chemical potential is estimated using the Widom insertion method.

IV. NUMERICAL RESULTS FOR RELEASE WAVES

We compare in this section the results for the different techniques presented in Sections II and III, for a release in a Lennard-Jones system (argon). The aim is to assess whether the release is indeed isentropic, and also to demonstrate that approximate equilibrium computations for small systems (successive Hugoniostat, isentropic integration) can approximate the isentrope curve obtained from the more rigourous and costly thermodynamic integration technique.

A. Numerical parameters 1. Initial state.

We consider argon in an initial shocked state, located on the Hugoniot curve for a compression such that L x = cL with c = 0.65, and corresponding to T = 1758 K and P = 1.7 × 10 10 Pa. At these thermodynamic conditions, the system is in a liquid state. The interactions within noble gas atoms are well-described by a Lennard-Jones potential:

V (q 1 , . . . , q N ) = 1≤i<j≤N v(|q i -q j |), v(r) = 4ε σ r 12 - σ r 6 .
In the case of Argon, ε/k B = 120 K and σ = 3.405 Å. The cut-off radius for the Lennard-Jones interaction is here r cut = 2.5σ.

2. Nonequilibrium simulations.

In order to reach this initial state before performing the NEMD release, a preliminary hugoniostat simulation is run for a system of 50×50×500 unit cells, using periodic boundary conditions. Then, the boundary conditions in the longitudinal direction are removed, and the system evolves according to the Hamiltonian dynamics. Profiles of thermodynamic quantities are computed every 0.25 ps for the post-processing procedure described at the end of Section II B.

Equilibrium simulations

Equilibrium computations have been performed with a system composed of N = 4000 atoms, starting in a FCC crystal geometry before melting, using periodic boundary conditions in all directions.

a. Thermodynamic integration. As shown in Section III A, the search of states having the same entropy can be performed using thermodynamic integration, which amounts to performing many equilibrium simulations. The canonical sampling for a given set of parameters (λ 1 , λ 2 ) is done with a Langevin dynamics for N steps = 2 17 time steps, with ∆t = 2 × 10 -15 s, and a friction coefficient γ = 10 13 s -1 .

First, the entropy variation along the isothermal release is computed, with canonical samplings along the path (0, 0) → (0, λ 1 ) with λ 1 = 0.54 (using M + 1 = 15 states). Then, for each compression of interest, the isochore cooling is performed using temperature steps ∆T = -25 K for expansions λ 1 ≤ 0.25, and ∆T = -50 K for states λ 1 ≥ 0.25 (these paths can be restated in terms of λ 2 ∈ [0, 1] upon considering a temperature modification ∆T depending on the compression). The numerical integration for computing the value of A is finally performed using the trapezoidal rule.

Error estimates on the canonical samplings are obtained with block averaging. 21 In all the cases considered, the statistical error (as measured using the 95% confidence interval associated with the variance computed from block-averaging) is inferior to 1%. Therefore, the entropy difference is computed within 1% errors. For fixed λ 1 , the state λ 2 such that S(λ 1 , λ 2 ) is constant is then known with an error depending on the local value of the partial derivative of S with respect to λ 2 . This error can immediately be reformulated as an error on the estimated temperature. The error on the computed pressure is the error arising from the error on the state λ 2 , plus the sampling error. It is found to be at most 2%. d. Entropy evaluation. The test particle insertion method used to evaluate the chemical potential requires many more iterations than the other equilibrium techniques. In the same framework as for isentropic integration (see Appendix C), N steps = 5 × 10 8 iterations were needed to obtain a satisfactory convergence. The statistical error on the calculated entropy (as measured using the 95% confidence interval associated with the variance computed from block-averaging) is estimated to be inferior to 1.2 %.

B. Discussion of the numerical results

Release waves are presented in Figures 345in three different diagrams, (P, ρ),(P, T ) and (T, ρ). It can clearly be seen that the results coming from the three equilibrium techniques of isentropic simulations are very close. This shows that provided the relative volume change parameter is carefully chosen in either the successive hugoniostat method or the isentropic integration, the propagating error remains at a low level; these methods can then be as accurate as the more rigorous and costly method of thermodynamic integration. Moreover, evaluating the chemical potential, we have computed absolute value of the entropy at three densities, ρ = 2780 kg.m -3 , ρ = 2190 kg.m -3 , and ρ = 1806 kg.m -3 . The corresponding values, 83.2±0.95 J.mol -1 , 83.1±0.42 J.mol -1 , 83.2±0.18 J.mol -1 , confirm that the entropy is indeed constant (within the error bars) on the calculated curve, validating once again the different methods.

The comparison with the results of the expansion of the liquid using nonequilibrium MD is also fruitful. The overall agreement is fair enough, which means that release waves are indeed almost isentropic. However, it can be noticed that the temperature is not predicted correctly.

While the different curves look very similar in a (P, ρ) diagrams, some discrepancies appear in the (T, P ) diagram, which are even more obvious in the (T, ρ) diagram. Indeed, for the latter diagram, the observed temperatures around the final density are greater than the error bars. The thermodynamic path followed by the system during its release exhibits systematically a higher temperature than the one of an isentropic process. This means that the release of a monoatomic liquid is not strictly isentropic, as is sometimes expected or assumed.

Recall however that a release is expected to be isentropic only for non-turbulent flows of non-viscous fluids. In the case considered here, the fluid has a finite, non zero viscosity, and therefore dissipates energy under the form of heat. As a consequence, the temperature should be higher than for an isentropic release. A tentative of evaluation of this effect is presented below. On the Hugoniot curve, 33 viscous effects can be introduced in the Hugoniot relation 23 by means of the "viscous pressure" π as

1 2 π(V 0 -V) = U -U 0 - 1 2 (P + P 0 )(V 0 -V ), ( 15
)
where π is defined as

π = -ν du dx , (16)
ν being the fluid viscosity and du dx the velocity gradient. Taking the viscosity of the argon fluid at T = 700 K and P = 1 GPa (the most extremes conditions of available thermodynamic tables), and considering an average velocity gradient (taken during the fluid release), we find a temperature elevation of a few Kelvins. Considering that the pressure is much higher in our simulation, and therefore that the viscosity should be also greater, the actual temperature increase due to the finite viscosity should rather be of the order of a few tens of Kelvins, which is consistent with what can be observed in our numerical results.

Finally, a purely nonequilibrium effect has been observed during the NEMD simulations, that also leads to a temperature increase in the system. Indeed, gradients of thermodynamic and kinematic quantities are large at the first stages of the release, when the hot and dense material is in contact with void. The thermodynamic path followed by the system at those early stages of the simulation does not correspond to the thermodynamic path followed when the release has reached its self-similarity regime. Some equilibration time is needed for some steady-state regime to be reached. We evaluated this time to be around 5 picoseconds.

V. CONCLUSION

We have presented or recalled several equilibrium methods to compute isentropic processes in the high pressure regime, either for compressions or releases. These methods, although very different in nature, lead to similar results when applied to the release of a monoatomic liquid.

We have then compared release waves computed with these equilibrium methods with the nonequilibrium simulation of the release process. The results show that the release is almost, but not strictly isentropic, the system's temperature being systematically greater than the one of the isentropic process. This is the consequence of two effects. First, the fluid actually has a finite viscosity and therefore dissipates heat, leading to a temperature increase.

To our knowledge, this is the first time that this effect has been quantified rigorously using molecular dynamics simulations. Moreover, the thermodynamic path followed by the system during its release takes some time to reach a converged profile. We anticipate that these effects will be enhanced in the case of a more complex fluid, for example in the case of a the release of detonation product. Therefore, the assumption that release waves are isentropic should be carefully verified in each case.

Reformulation of the problem in a fixed geometry

From a computational viewpoint, it is more convenient to work with a fixed simulation domain. For instance, the unperturbed domain V (0) may be used to fix the geometry of the system. The volume variations are then rephrased as variations in the interaction scale between the particles in the direction of compression or release. In the same vein, the temperature may be kept constant, upon rescaling the interactions strengh by a factor depending on the temperature variation. Introducing the rescaled potential energy for a configuration q = (x, y, z):

U λ 1 ,λ 2 (q) = 1 1 + λ 2 δT U((1 + λ 1 )x, y, z).
and the associated Hamiltonian

H λ 1 ,λ 2 (q, p) = U λ 1 (q) + 1 1 + λ 2 δT N i=1 p 2 i 2m i ,
canonical averages for a volume V (λ 1 ) at a temperature T (λ 2 ) can be reformulated as canonical averages in terms of the rescaled Hamiltonian H λ 1 ,λ 2 at the reference state at volume V (0) and temperature T (0). More precisely,

H V (λ 1 ),T (λ 2 ) = 3N 2 k B T (λ 2 ) + (1 + λ 2 δT ) U λ 1 ,λ 2 λ 1 ,λ 2 ,
where

f λ 1 ,λ 2 = Ω(0)
f (q, p) e -H λ 1 ,λ 2 (q,p)/k B T (0) dq dp

Ω(0)
e -H λ 1 ,λ 2 (q,p)/k B T (0) dq dp .

It is then easily seen that

S(λ 1 , λ 2 ) -S(0, 0) = 3N 2 ln(1 + λ 2 δT ) + N ln(1 + λ 1 ) +β [ U λ 1 ,λ 2 λ 1 ,λ 2 -U 0,0 0,0 ] (A1) +A(λ 1 , λ 2 ), with A(λ 1 , λ 2 ) = ln     V (0) N e -βU λ 1 ,λ 2 (q) dq V (0) N e -βU 0,0 (q) dq     .
In the above expression of the entropy difference, the first line is the ideal gas contribution to the entropy difference. As a consistency check, we can verify that the entropy increases when the volume or the temperature is increased, as expected. The terms on the second and third lines in Eq. (A1) are the "excess" contributions associated with the potential interaction energy.

Numerical evaluation of the different terms

To estimate S, two quantities are required:

(i) averages • λ 1 ,λ 2 with respect to the Hamiltonian H λ 1 ,λ 2 are computed using standard sampling techniques such as a Langevin dynamics at an inverse temperature β, implemented using the so-called BBK algorithm. 24 Of course, many other sampling techniques could be used to estimate this canonical average, in particular Nosé-Hoover dynamics 25,26 or Metropolis-Hastings schemes 27,28 (see Ref. 29 for a mathematical review on sampling methods in the context of molecular simulation);

(ii) the term A(λ 1 , λ 2 ) requires more care in its estimation. Since this term is a ratio of partition functions, standard techniques used for the computation of free energy differences may be used. We resorted to thermodynamic integration, 30 in which case the function is rewritten as the integral of some canonical averages:

A(λ 1 , λ 2 ) = λ 1 0 ∂A ∂ λ 1 (x, 0) dx + λ 2 0 ∂A ∂ λ 2 (λ 1 , x) dx, with ∂A ∂λ 2 (λ 1 , λ 2 ) = β δT 1 + λ 2 δT U λ 1 ,λ 2 λ 1 ,λ 2 , (A2) and ∂A ∂λ 1 (λ 1 , λ 2 ) = x • ∇ x U((1 + λ 1 )x, y, z) λ 1 ,λ 2 . (A3)
In conclusion, the numerical procedure consists in first estimating the derivatives of the function A and the average potential energy, for as many points as required on the thermodynamic path chosen. Approximations of S can then be obtained thanks to (A1), after a numerical integration to obtain A. The entropy difference along the path is then plotted, and fixing the volume change λ 1 , the temperature variation is chosen such that the entropy difference is 0. This determines λ 2 as a function of λ 1 .

APPENDIX B: RELATIONSHIP BETWEEN THE HUGONIOT AND THE ISENTROPE CURVES AT THE POLE

We present in this Appendix two proofs of the fact that the isentrope curve and the Hugoniot agree at order 3 in the volume change. The first one is a standard thermodynamic proof, but the second one, based on statistical physics relations, is new to the best of our knowledge.

Standard thermodynamic proof

For the sake of completeness, we reproduce here the proof of Ref. 31. From the thermodynamic relation T dS = dU +P dV and from the Hugoniot relation U = U 0 + 1 2 (P +P 0 )(V 0 -V ), the entropy variation along the Hugoniot curve can be computed.

d 2 P Hug dV 2 + 1 2 (V 0 -V ) d 3 P Hug dV 3 . (B3)
At the initial state (denoted with a subscript 0), that is, in the limit V → V 0 , it holds:

dS Hug dV 0 = 0, d 2 S Hug dV 2 0 = 0, T 0 d 3 S Hug dV 3 0 = - 1 2
d 2 P dV 2 0 = 0.
The entropy variation along the Hugoniot curve is therefore of order three in volume, and the Hugoniot and the isentropic curves are osculatory.

A statistical physics derivation

Without loss of generality (and for notational simplicity), we may set H(0, 0) = S(0, 0) = 0 since we are only interested in differences of F and S.

a. Some useful relations.

The derivatives of the function A are useful for comparing the Hugoniot and the isentrope relations. The average xx component of the pressure tensor for the volume V (λ 1 ) and the temperature T (λ 2 ) is obtained by averaging the observable

P xx (q, p) = 1 V (λ 1 ) (Nk B T (λ 2 ) -x • ∇ x U(q)) .
Therefore,

P xx V (λ 1 ),T (λ 2 ) = N (1 + λ 1 )V (0) k B T (λ 2 ) - 1 V (λ 1 ) V (λ 1 ) N x • ∇ x U(q) e -U (q)/k B T (λ 2 ) dq V (λ 1 ) N e -U (q)/k B T (λ 2 ) dq = N βV (0) 1 + λ 2 δT 1 + λ 1 - 1 + λ 1 V (λ 1 ) V (0) N x • ∇ x U((1 + λ 1 )x, y, z) e -βU λ 1 ,λ 2 (q) dq V (0) N e -βU λ 1 ,λ 2
. This shows that, using (A3),

P xx V (λ 1 ),T (λ 2 ) = N βV (0) 1 + λ 2 δT 1 + λ 1 + 1 + λ 2 δT βV (0) ∂A ∂ λ 1 (λ 1 , λ 2 ). b. Hugoniot curve.
With the above computations, it is easily seen that the Hugoniot relation ( 8) can be restated as

βH(λ 1 , λ 2 ) = 3N 2 λ 2 δT + β [(1 + λ 2 δT ) U λ 1 ,λ 2 λ 1 ,λ 2 -U 0,0 0,0 ] + Nλ 1 2 1 + λ 2 δT 1 + λ 1 + 1 + λ 1 2 ∂A ∂ λ 1 (0, 0) + (1 + λ 2 δT ) ∂A ∂ λ 1 (λ 1 , λ 2 ) .
c. Comparison between the Hugoniot and the isentrope.

We now Taylor expand the difference βH(λ 1 , λ 2 ) -S(λ 1 , λ 2 ) up to the third order, i.e.

neglecting a remainder term r(λ 1 , λ 2 ) which is such that |r [START_REF] Lyzenga | proceedings of the AIP conference Shock Compression and Condensed Matter[END_REF] . We denote such remainders by O(λ 3 ) in the sequel. It holds

(λ 1 , λ 2 )| ≤ C(|λ 1 | + λ 2 |)
βH(λ 1 , λ 2 ) -S(λ 1 , λ 2 ) = 3N 2 (λ 2 δT -ln(1 + λ 2 δT )) + N λ 1 2 1 + 1 1 + λ 1 -ln(1 + λ 1 ) +λ 1 λ 2 δT 2 N 1 + λ 1 + ∂A ∂ λ 1 (λ 1 , λ 2 ) + βλ 2 δT U λ 1 ,λ 2 λ 1 ,λ 2 + λ 1 2 ∂A ∂ λ 1 (0, 0) + ∂A ∂ λ 1 (λ 1 , λ 2 ) -A(λ 1 , λ 2 ).
Introducing the notation

A i = ∂A ∂ λ i (0, 0), A ij = ∂ 2 A ∂ λ i ∂ λ j (0, 0),
the Taylor expansions of the function A and its first derivatives at an arbitrary state (λ 1 , λ 2 ) read (using A(0, 0) = 0):

A(λ 1 , λ 2 ) = λ 1 A 1 + λ 2 A 2 + λ 2 1 2 A 11 + λ 1 λ 2 A 12 + λ 2 2 2 A 22 + O(λ 3 ), ∂A ∂ λ i (λ 1 , λ 2 ) = A i + λ 1 A i1 + λ 2 A i2 + O(λ 2 ).
With these Taylor expansions and the relation (A2), it is straightforward to show that

βH(λ 1 , λ 2 ) -S(λ 1 , λ 2 ) = 3N 4 λ 2 2 δT 2 + λ 1 λ 2 δT 2 (N + A 1 ) + λ 2 (1 + λ 2 δT ) ∂A ∂ λ 2 (λ 1 , λ 2 ) = λ 2 δT 2 λ 1 N + A 1 + A 12 δT + λ 2 3N 2 δT + 2A 2 + A 22 δT + O(λ 3 ).
Using (A2), the derivatives of the entropy differences can be computed:

∂S ∂λ 1 (λ 1 , λ 2 ) = N 1 + λ 1 + ∂A ∂λ 1 (λ 1 , λ 2 ) + 1 + λ 2 δT δT ∂ 2 A ∂λ 1 ∂λ 2 (λ 1 , λ 2 ), ∂S ∂λ 2 (λ 1 , λ 2 ) = 3NδT 2(1 + λ 2 δT ) + 2 ∂A ∂λ 2 (λ 1 , λ 2 ) + 1 + λ 2 δT δT ∂ 2 A ∂ 2 λ 2 (λ 1 , λ 2 ).
This shows that

βH(λ 1 , λ 2 ) -S(λ 1 , λ 2 ) = λ 2 δT 2 λ 1 ∂S ∂λ 1 (0, 0) + λ 2 ∂S ∂λ 2 (0, 0) + O(λ 3 ), (B4) so that, since S(λ 1 , λ 2 ) = S(0, 0) + λ 1 ∂S ∂λ 1 (0, 0) + λ 2 ∂S ∂λ 2 (0, 0) + O(λ 2 ),
and S(0, 0) = 0, it holds

βH(λ 1 , λ 2 ) -1 + λ 2 δT 2 S(λ 1 , λ 2 ) = O(λ 3 ). ( B5
)
This relation shows immediately that H(λ 1 , λ 2 ) = O(λ 3 ) on the isentrope, and so, the initial slopes of the curves, and their first derivatives, coincide.

The numerical implementation of this method is done as follows. The partial derivative of the pressure with respect to the energy is first computed with a Monte Carlo simulation for the given initial conditions (N, V 1 , T 1 ). The temperature T 2 is then evaluated from Eq. (C1).

The partial derivative is next computed at volume V 2 to predict the next temperature.

Proceeding incrementally, the whole isentrope curve can be constructed.

The numerical results presented in this work have been obtained by performing canonical samplings with a Metropolis algorithm, using the Monte-Carlo Gibbs code. 34 Partial derivatives have been computed in the NVT ensemble. The convergence of simple thermodynamic averages was generally obtained after N steps = 10 7 iterations, but derivative properties (related to the covariance of some observables) required about N steps = 10 8 iterations for a satisfactory convergence. Error estimates on the canonical samplings have been obtained

with block averaging, 21 and the error propagation estimated along the integration scheme has been computed using standard propagation rules. In all the cases considered, the statistical error (as measured using the 95% confidence interval associated with the variance computed from block-averaging) on the predicted temperature on the isentrope curve is inferior to 1.5 %.

  b. Successive Hugoniostat. Successive Hugoniostat simulations have been performed with a Langevin version of the Hugoniostat method (see Eq. (11) in Ref. 22, with the parameters ξ = 10 12 s -1 and ν = 10 12 s -1 ). Trajectories of N steps = 50, 000 timesteps at each compression are considered, with a timestep ∆t = 5 × 10 -16 s. The relative volume change δV /V 0 from one point on the curve to another is set to 0.01. c. Isentropic integration. See Appendix C.

Figure 1 :

 1 Figure 1: (a) (color online) Non-Equilibrium Molecular Dynamics of isentropic release waves. The four pictures represent snapshots of the system during the release process, the expansion proceeding in the longitudinal x direction. Atoms are colored according to their potential energies (scaling corresponding to -1.38 × 10 -20 J for blue up to 7.55 × 10 -20 J for red).

Figure 1 :

 1 Figure 1: (b) (color online) Density profiles taken at different times of the simulation (from blue to red as time increases).

Figure 2 :

 2 Figure2: Path in the (λ 1 , λ 2 ) space used to compute states with the same entropy as the initial state. Each cross represents some equilibrium canonical sampling along the thermodynamic path. First, the isothermal expansion is performed (horizontal line in the diagram), starting from the initial state (0, 0), until the required density is reached. The entropy of the state (λ 1 , 0) is S init + ∆S expansion . Then, an isochore cooling is performed (λ 1 is kept fixed; vertical line in the diagram), until the entropy difference during this process is the opposite of the entropy variation found in the expansion part. The final state (λ 1 , λ 2 ), located at the intersection of the curve ∆S = 0 and the vertical line, has then the same entropy as the initial state.

Figure 3 :

 3 Figure 3: (color online) Isentropic release in a (P, ρ) diagram. Symbols represent resultsfrom equilibrium methods, red diamonds for the successive hugoniostat, blue squares for the thermodynamic integration and yellow triangles for the entropy integration. NEMD results are plotted in green, the width of the so-obtained tube corresponding to the error bars. The arrow indicates the path followed during the release.

Figure 4 :

 4 Figure 4: (color online) Isentropic release in a (T, P ) diagram. The symbols are the same as in Figure 3. Notice that there is slight deviation of the NEMD results for the lowest temperatures.

Figure 5 :

 5 Figure 5: (color online) Isentropic release in a (T, ρ) diagram. The symbols are the

FIG. 1 :

 1 FIG.1:(a)

Since the Hugoniot curve and the isentrope are close enough, we assume that the impact of the viscosity has roughly the same amplitude on the isentrope curve.

The Monte-Carlo Gibbs code is owned by the Institut Francais du Pétrole, the Université Paris-Sud, and the CNRS, and developped in collaboration with the CEA.

APPENDIX C: PRECISIONS ON THE ISENTROPIC INTEGRATION

Several numerical schemes may be used to integrate (11). The simplest one consists in approximating the integral appearing in the exponential factor with a Riemman formula using the value of the integrated function on the left side of the interval:

Of course, higher order integration methods could be used.

It remains to decide how to compute the derivative ∂P ∂U V 1 . Finite differences may be used to this end, but this would require at least two very carefully converged simulations with volumes V 1 ± ∆V . It seems more appealing to compute the partial derivative using standard fluctuations formulas: 14,32

and

where

is the specific heat at constant volume, and the pressure observable for a simulation domain of volume V 1 reads

The partial derivative ∂P/∂U can then be evaluated in a single simulation at (N, V 1 , T 1 ) using (C2)-(C3).