Speed of sound in solid molecular hydrogen-deuterium: Quantum Molecular Dynamics Approximation

Uniformity of the solid layer is one of the critical points for an efficient ignition of the Deuterium-Tritium (DT) target. During the compression process this layer, perturbations grow as the Rayleigh-Taylor instability. Knowing the mechanical properties of this layer and its thermo-mechanical limits is necessary if we want to control or to minimize these instabilities. In this work we have used a simplified approach, replacing the DT ice system with a mixture of hydrogen-deuterium (HD) because beta decay of tritium complicates the analysis in the former case. Through simulation with ab initio methods we have calculated the elastic constants, the bulk modulus and sound velocity for hydrogen isotopes in solid molecular state. In this work we present the results for hydrogen-deuterium mixtures 50%-50%, at 15 K and with a compression which covers the range of 1 to 15 GPa. This system is interesting for study the early stages of the dynamic compression and provides conditions that are close to the manufacture of DT target in inertial confinement fusion. Discontinuities in the curve that have been observed on pure hydrogen, which are associated with phase transitions and the phase hysteresis.


Introduction
The study of the mechanical properties of solid molecular hydrogen (SMH) is of great importance for the development of inertial confinement fusion (ICF) system [1]. In the ICF system, fusion is achieved by compressing a sphere of Deuterium Tritium (DT) ice using high-power lasers [1]. The compression process is modeled to increase ignition efficiency of the target. Our purpose is to model systems composed of H, D and T and mixtures of them. In recent studies we have compared our results with experimental data, obtaining a good agreement, both for the structural changes undergone by the hydrogen and for the properties mechanical bulk modulus [2]. In this way we have demonstrated the reliability of the simulation in simple systems composed of hydrogen, deuterium and tritium. This work allows us to understand the changes in mechanical properties such as bulk modulus and sound velocity for systems with mixtures of H and D, a preliminary step for more detailed study for combinations of DT and HDT. Different studies have been conducted to analyze the mixture of isotopes HD. Phase diagram of HD is very different from having the H 2 and D 2 separately, Cui et al [3] show that the phase diagram has a different behavior in the pure H curve representing the transition from phase I and II for temperatures below 100 K. Several authors demonstrate the "broken symmetry" associated with the restriction of rotation of the rotor of the molecule HD [4,5]. The HD has been insufficiently studied, it becomes important to analyze the temperature and pressure region near to manufacture and early compression DT ice in ICF, for this we choose 15 K as initial temperature and pressure range of 1 to 15 GPa, to thereby representing dynamic compression system.
The method used in this study was applied in the study of various properties of SMH, describing the three phases of solid H pure [6], calculating the bulk modulus tritium [2] and relating changes in the sound velocity curve vs. pressure changes in the structural configuration of the solid [7]. The interest of these simulations is focused on generating parameters describing the behavior of the solid in conditions seldom studied experimentally. The temperature of interest in this article is near to that used in the storage and handling of DT ice spheres in ICF systems [8], analyzing factors that can induce problems in the solid layer of the target in early stage compression. For commercial applications targets can be produced in situ [9] or stored for later use; hence, storing thousands of targets is very relevant for systems of magnetic or inertial. The results obtained in this work also apply to models used to study geophysics and formation of exoplanets and giants planets, such as Jupiter and Saturn [10,11]. Finally, the paper presents the changes that occur in the bulk modulus and sound velocity of the HD solid. We show the methodology in section 2, the results of bulk and sound velocity in the range of 1-15 GPa with 15 K in section 3, and our conclusions and hypotheses for future work in Section 4.

Methodology
For HD simulation molecular solid we use a box with 576 atoms, 288 HD molecules, each molecule is composed of one hydrogen and one deuterium. The concentration of isotopes is variable in our code but the proposed concentration 1:1 for a sample basis and compared with pure H. The SIESTA code [12] is used, based localized orbitals on Density Functional Theory, the complete method has been described in previous papers [2,6,7]. To get the lowest-energy configuration of the HD system we use the method of conjugate gradients, whereby an initial structure that minimize the total energy of the system. Then we apply the Nose-Hoover thermostat and the Parrinello-Rahman barostat the crystal structure is obtained for a determinate pressure and temperature. For an understanding of the dynamic compression experiments, we applied the pressure as ramps. The initial pressure is 1 GPa maintained for 1 ps, with this structure we proceed to calculate the bulk modulus (B) using Hooke's law for this pressure applying deformation of 1% and observing stress in the structure. Then without deformation structure is compressed again applying 2 GPa pressure during the same period 1 ps, creating a new structure for each step of the pressure ramp, 1 GPa as increased variation of P, until the 15 GPa (see figure 1). The simulation is made with SIESTA standard basis double zeta polarized (DZP) and LDA functional, with MeshCutoff 100 Ry. The pseudo-potential used is Troullier-Martin type. As a first approximation use the model of isotropic transverse material. HD is a molecular solid for low pressure in the Phase I rotor associated to molecule is free, the molecular rotation decreases with increasing pressure and rotor is restricted. The C ij allow us to calculate each of the directions of deformation getting information about the speed of sound propagation in different aspects of the material [7].
where σ is the stress calculated by simulation with the deformation ϵ. Calculating the bulk modulus is achieved through equation of the Voigt average [13]: 9B=(C 11 +C 22 +C 33 )+2(C 12 +C 13 +C 23 ) For us C 22 =C 11 , C 12 =C 13 and the sound velocity is, as a first approximation: Vs = √ (B / ρ), where ρ is solid density and Vs is the speed of sound in solid.

Results
The figure 2 shows a change in the magnitude of the bulk modulus and sound velocity. HD bulk modulus is greater than that of H (see figure 2a) and the difference increases with pressure. In the case of the HD sound velocity (see figure 2b) it is lower than that of H in the analyzed range. Variations in pressure below 5 GPa in the graph for bulk modulus and the sound velocity for HD, would be related to the phase transition from solid-solid [7]. There are two discontinuities in red and blue HD curves, one for 2 to 5 GPa and for P> 9 GPa, this last is larger than the first. This can be explained due to variations in the phase diagram of HD, and the same points are not present in the phase change of H [3]. For 1 and 2 ps the bulk modulus is similar between 5 GPa <P <9 GPa, while the rest of pressures the bulk modulus is different, this means that the range between 5 to 9 GPa change in time of application of pressure ramp (here 1ps and 2 ps), does not alter the system configuration and therefore the parameters calculated from this. In the figure 2b the complexity of the curve for HD remains but also appears some differences in the magnitude of the sound velocity compared to that of H, this effect arises because the velocity depends on the solid density and this quantity is larger for HD than for H (50% larger than the original H). This density increase according to the ratio of deuterium, causing the HD phase diagram depends on the isotope ratio in the solid, changing its mechanical properties for a given P and T compared to H.
The discontinuities in the sound velocity curve (see figure 2b) have been associated with changes in the restriction of the molecules rotation in the solid H [7], changes being related to solid-solid phase [14]. For P < 2 GPa H 2 molecule center move outside the lattice point and the rotor is free, these effects have been studied in previous works [7], for 2 < P < 10 GPa H 2 molecule center coincides with lattice point and rotor is free, to P> 12 GPa the H 2 rotor is restricted [6,7]. The difference between 1ps and 2 ps is only appreciable in a range of pressure 5 GPa <P <10 GPa, for the rest of the curve the difference between H and HD does not exist. This shows that the application time of pressure ramp intervenes in the ability of the solid to respond to the change in the phase transition, this can be explained by the concept of hysteresis phase. This phenomenon is not present in experiments static compression [15] but is very important in dynamic compression experiments, especially compression laser ICF systems such as NIF experiments [16], where the pressure profile varies with time. The hysteresis phase has important implications for the process in which the target of DT ice is compressed, and has already been observed in experiments iron dynamic compression [17]. Remains to be studied more in depth this hysteresis of phase phenomenon in compression of SMH isotope, and see the effects on the generation of inhomogeneities in the front of the shock wave [18] that can alter conditions plasma heating.

Conclusions
The method presented in this work allows control of the H-D ratio in the volume, so we can get different relative concentrations and associate them with changes in their mechanical properties. We have shown that our methodology represents very well the experimental results for a wide range of temperature, 300 K [2] and 15 K (see figure 2) and a broad pressure range from 1 GPa to 200 GPa [6,7] ( in this work between 1 and 15 GPa). The H:D ratio in the solid changes the magnitude of the sound velocity, which in turn changes the shape of the sound velocity and bulk modulus curves, producing variations in mechanical properties compared with that of pure H. The difference between 1 to 2 ps allows a deeper understanding of the effect of hysteresis in phase and pressure regions near the phase transition temperature. Getting to know the response of the solid-solid phase transition as a function of time of application P-ramps gives us a map of delays in the response of solid-liquid or solid-plasma transition in dynamic compression experiments.