Structural and transport properties of ammonia along the principal Hugoniot

We investigate, via quantum molecular dynamics simulations, the structural and transport properties of ammonia along the principal Hugoniot for temperatures up to 10 eV and densities up to 2.6 g/cm3. With the analysis of the molecular dynamics trajectories by use of the bond auto-correlation function, we identify three distinct pressure-temperature regions for local chemical structures of ammonia. We derive the diffusivity and viscosity of strong correlated ammonia with high accuracy through fitting the velocity and stress-tensor autocorrelation functions with complex functional form which includes structures and multiple time scales. The statistical error of the transport properties is estimated. It is shown that the diffusivity and viscosity behave in a distinctly different manner at these three regimes and thus present complex features. In the molecular fluid regime, the hydrogen atoms have almost the similar diffusivity as nitrogen and the viscosity is dominated by the kinetic contribution. When entering into the mixture regime, the transport behavior of the system remarkably changes due to the stronger ionic coupling, and the viscosity is determined to decrease gradually and achieve minimum at about 2.0 g/cm3 on the Hugoniot. In the plasma regime, the hydrogen atoms diffuse at least twice as fast as the nitrogen atoms.

As ammonia is one of the major constituent of the "hot ice" layers in Uranus and Neptune 1,2 , its properties such as diffusivity, viscosity and structure have considerable impact on the internal evolution and magnetic field of these giant planets. For example, the transport character of warm dense ammonia is essential for understanding the composition of giant planets, and also the magnetohydrodynamic models of these planets depend on the knowledge of diffusion and viscosity of the fluid 3 . In this weakly hydrogen-bonded liquid, chemical processes are expected to take place at deep planetary conditions (several 100 GPa and several 1000 K), which could induce the change of structure and transport properties of ammonia. However, the chemistry and transport properties of ammonia at extreme conditions still remain elusive up to date.
Much of the experimental and theoretical research have advanced the determination and understanding of physical properties of ammonia under extreme conditions. On experimental side, some dynamic techniques, such as standard explosive technique and two-stage light-gas gun, have been used to derive the Hugoniot of ammonia up to 39 GPa first, and later to a higher pressure of 64 GPa 4,5 . Furthermore, Radousky et al. measured the shock temperature at the pressures of 48 and 61 GPa for diagnosing the physics occurring in ammonia 6 . Nonmetal-to-Metal transition for ammonia was found in the electrical conductivity measurements 7,8 . Recently, the phase diagram up to 60 GPa and 2500 K has been studied by several static experiments with diamond anvil cells [9][10][11][12] . However, experiments were restricted to certain areas of the phase diagram due to the limitations of technique. More detailed and broad information at extreme conditions come from theoretical calculations. Ab initio molecular dynamics simulations have predicted the phase transition to a protonic conductor for ammonia above 60 GPa and 1200 K 13 . First-principles calculations have shown that ammonia molecules chemically dissociate to N 2 and H 2 above approximately 7 GPa and 900 K, which was thought to be driven by the entropy of mixing term in the free energy formulation 11 . In addition, ab initio molecular dynamics simulations have been used to examine the equation of state (EOS) and phase diagram covering pressures up to 330 GPa and a temperature range from 500 K to 10000 K 14 . We have recently investigated the thermophysical properties of shocked liquid ammonia up to the pressure of 1.3 TPa and temperature 120000 K using quantum molecular dynamics simulations. It was found that the liquid ammonia undergoes three structural regimes along the principal Hugoniot, accompanying the transformation to a metallic state 15 . Later, Mulford et al. predicted the ammonia EOS based on Cheetah code 16 .
In this paper, we study the structure and transport properties of shocked ammonia under extensive pressure-temperature regimes with temperatures up to 10 eV and densities 2.6 g/cm 3 systematically. Specially, the structural transformation is first studied by the pair distribution function and bond auto-correlation function (BACF), and then the changes in the diffusion coefficients and viscosity of the system at different structural regimes are discussed in detail.

Results
To confirm the internal energy and pressure achieve good convergence, we first conduct convergence tests for the particle number, the plane-wave cutoff and k-point Brillouin sampling at several P − T conditions, which include (i) QMD simulations with 64-molecule supercells, the Γ point and a 1000 eV plane-wave cutoff; (ii) QMD simulations with 64-molecule supercells, the Γ point and a 550 eV plane-wave cutoff; (iii) QMD simulations with 27-molecule supercells, the Γ point and a 550 eV plane-wave cutoff; (iv) QMD simulations with 27-molecule supercells, a 2 × 2 × 2 k-point grid and a 550 eV plane-wave cutoff. We take the Hugoniot point of 2.0 g/cm 3 and 6036 K as an example, and present the comparison results in Table 1. It can be seen that the internal energies and pressures are converged to less than 1% for all of these choices of parameters.
We further examine the influence of the exchange-correlation potential (PW91, PBE, and LDA) on the equations of state. In Table 1, the pressures and internal energies for different exchange-correlation potentials are also included. The data from QMD simulations with PW91 functional show good agreement with those from PBE functional, while the data from LDA functional present prominent discrepancies, with the relative deviations between PW91 and LDA results as large as 9.9% and 15.6% for pressure and internal energy, respectively. The large deviations may arise from the inaccurate modeling of the breaking and forming of chemical bonds for LDA functional. And thus we present the impact of the exchange-correlation potential on the pair distribution functions (PDFs) in Fig. 1. The PW91 and PBE functionals yield almost the same PDFs, while the LDA functional predicts reduced peak in g NH (r) and higher peak in g NN (r), which signify the larger fraction of dissociated ammonia molecules and formation of N-N bonding. Furthermore, we look insight into the N-N bonding in ammonia through calculating the coordination numbers of N, which is a weighted integral over the PDFs g(r) of N-N The fraction of the ions bounded to a nitrogen molecule is twice the value of K(r) at the first peak in g(r), which is found at around 1.30Å. The plots are inset into Fig. 1(b). The coordination number of N is 0.136 for the PW91 and PBE functionals, while 0.175 for the LDA functional. Thus the fraction of N-N bonding is 27.2% and 35% for PW91/PBE functional and LDA functional, respectively.
In the following, we examine the chemistry that occurs in ammonia along the principal Hugoniot. The Hugoniot points have been obtained using QMD simulations in our previous studies 15 . Figure 2 presents the BACF of ammonia for three types of bonds: N-N, N-H and H-H for different Hugoniot points. At 1.4 g/cm 3 and 2295 K, the bond lifetimes for N-N and H-H are very short, even close to zero. While N-H bonds are at least an order of magnitude longer. This demonstrates that the system mainly remains as ammonia molecules at lower densities and temperatures of the Hugniot. At the intermediate density and temperature, the bond lifetime for N-H becomes much shorter, smaller than 0.2 ps (time for which the correlation function is equal to 1/e). N-N bonds last about six times longer. It is indicated that ammonia molecules dissociate and small amounts of nitrogen bonds form at this state. The mixture phase could also be seen more obviously from the snapshot of the MD simulation of 1.7 g/cm 3 at 3562 K shown in Fig. 1(d). With increasing the density and temperature to 2.2 g/cm 3 and 19180 K, respectively, all three types of molecules are very short-lived and unstable, with the system entering into a plasma state.
In order to gain a further insight into the structural properties of ammonia in the mixture region, we present in Fig. 3 the PDF and BACF for three types of bonds along an isochore of 1.8 g/cm 3 at different temperatures. For the nitrogen-nitrogen PDF, the first peak is at the position of 1.38 Å, longer than the length of triple bond. This signifies the formation of the chain of nitrogen atoms. Here the expression of "the chain of nitrogen atoms" denotes the formation of single bond of N-N. Based on the bond-length criteria, the N-N bonding could be identified. We have used a cutoff radius of 1.70 Å, which is the maximum value for single N-N bond length, to construct a sphere about each nitrogen and the bonding between two atoms is defined as a series of overlapping spheres. It is found that at the corresponding states of 1.8 g/cm 3 and 1000/2000/3000 K, only the nitrogen chains with the typical length of two occur. In addition, we have calculated the fractions of N-N bond as 11.6%, 11.6% and 16% for these three states, respectively. Meanwhile, quite a few hydrogen molecules appear. In addition, the first minimal value of g NH (r) is slightly larger than zero, indicating the occurrence of dissociations of ammonia molecules. At lower temperature, the lifetimes for all three types of bonds are extremely long, indicative of the mixture state of ammonia. With increasing temperature, ammonia dissociates continuously and give rise to shorter lifetime for N-H bond.  Table 1. The pressure and internal energy at the Hugoniot point of 2.0 g/cm 3 and 6036 K for different particle number (N), the plane-wave cutoff (E cut), k-point Brillouin sampling and exchange-correlation potentials (XC).
As we have reported in ref. 15 , the principal Hugoniot of ammonia presents three characterized segmentations, which results from the structural changes. And thus we further explore the effect of structural changes on the diffusion and viscosity properties. Based on QMD simulations, we derive the time-dependent diffusion and viscosity coefficients of ammonia using Eqs (2-4) and (6). It should be noted that the mutual diffusion coefficients are calculated with the thermodynamic factor estimated by (〈Z*〉 + 〈(Z*) 2 〉)/(〈Z*〉 + 〈Z*〉 2 ). Numerical calculations based on average-atom model are performed to determine the average ion Z* 17 . We tabulate the thermodynamic factor for each Hugoniot point in Table 2. As expected, this factor is larger than one, due to an ambipolar effect, depending on the symmetry of the mixture. Through fitting the QMD results using the functional formula of Eqs (7), (9) and (11), the self-and mutual-diffusion coefficients, and viscosity are obtained with the correlation time between 5 to 100 fs. A total uncertainty less than 20% in mutual diffusion coefficient and viscosity come from the fitting and extrapolation to infinite time, while the error in the self diffusion coefficient is less than 4% since the particle average gives an additional N 1 factor.
We examine the convergence of the transport properties with respect to particle number, energy cutoff, k-points and exchange-correlation potential, as shown in Fig. 4. It can be seen that for different particle numbers and energy cutoff, both the self-diffusion coefficient of N and H are converged well, while for mutual-diffusion coefficient and viscosity, smaller particle number and energy may yields about 7.3% and 10.7% error, respectively. Thus in the whole simulations of this paper we set the particle number as 64 and energy cutoff as 1000 eV. For the different K point sets of '111' and '222' , all of the self and mutual diffusion coefficients converged very well, while the viscosity deviates a little, but still less than 5%. As for the impact of exchange-correlation potential on these transport coefficients, the PW91 functional and PBE functional give the very similar results, while the LDA functional diverge from them.
In Fig. 5, we show the fits to the velocity autocorrelation function (VACF) for ammonia at 1.9 g/cm 3 and 4791 K as well as the numerically integrated results. It is evident that correlated behavior and structures exist in these three kinds of VACF. The functional of Eqs (7) and (9) and their analytic integrals fit the raw data well. The diffusion coefficients coming from the fit to the VACF and its analytic integral are in statistical agreement with each other using the standard error in Eq. (13). We further find that there are different time scales and structures in N, H and N-H. In the fit for N, one type of restoring force could give a better fit via the R 2 criteria, while for H, a secondary restoring force is required, and even for N-H, all of the forces of N-N, H-H, and N-H are needed. The fit using a single exponential are also included (indicated as the solid black line). The single exponential cannot describe the oscillating behavior of the VACF and its integral, but goes to zero quickly. Thus the fitting to VACF using the single exponential yield higher diffusion coefficients than that of the fit using Eq. (7) (or Eq. (9)), with  Table 3. The diffusion error is the sum of both statistical and fit contributions. Figure 6 presents the normalized VACF for N and H for different densities along the principal Hugoniot of ammonia. It is shown that the VACFs for N and H exhibit different structural forms and multiple time scales because of the mass difference. Evidently H diffuses on a shorter time scale. The VACF for N has less structures than H, which means N atoms are less correlated. As the density and temperature increase, both of the VACF for N and H evolve from complex structured form to simple exponential form, which corresponds to the process that the system changes from strong-coupling state to moderate-coupling state.
For the complex mixture, we use the functional of Eq. (11) to fit the stress-stress tensor autocorrelation function (STACF) and its integral for extracting viscosity. In Fig. 7 we show the STACF and its integrated value for the ammonia at 1.9 g/cm 3 and 4791 K. In the fitting all of β i are set to be unity. The fitting line goes through the data smoothly, and decreases to zero beyond 50 fs. Therefore the integrated value in a smaller time window could give the final viscosity value. We find that the fit for the STACF and its analytic integral yield results of η = 1.10 ± 0.21 mPa s and η = 1.19 ± 0.13 mPa s, respectively, which are in statistical agreement with each other. The longer fluctuations of the STACF result in the larger error bar in fitting. In addition, the STACF show evident structure with an increase about 7.8 fs, which could not be well fitted with a single exponential functional. Figure 8 presents the diffusion coefficient and viscosity as functions of density along the principal Hugoniot. For comparison, we replot the principal Hugoniot of temperature as functions of density in Fig. 8(a). As can be seen from Fig. 8, similar as the principal Hugoniots, both the diffusion coefficients and viscosity fall into three segments with their respective features, in accordance with the three structural regions respectively. In the molecular region, the hydrogen atoms have essentially the same diffusivity as the nitrogen atoms and the mutual diffusivity is much lower. In the mixture region, nitrogen atoms have diffusivity coefficients typical of fluids (~1.0 × 10 4 cm 2 s −1 ). In the plasma phase, the hydrogen atoms diffuse at least fourfold as fast as nitrogen. In addition, we note that in the plasma regime the mutual diffusion coefficients calculated using Eq. (5) in terms of self-diffusion coefficient agree well with those from QMD simulations directly. This means that the interactions between nitrogen and hydrogen atoms are very weak.

Discussion
The complex behavior of viscosity is a consequence of the subtle interplay between coupling and kinetic effects. In kinetic regime with small coupling, the mechanism of bodily movement of particle is predominant and, as in a gas, the viscosity rises with temperature and density, reaching a maximum at about 1.6 g/cm 3 and  3 and 6036 K, which results from the similar contribution from the two competitive mechanisms. With the system becoming strongly coupled, the viscosity of ammonia behaves like dense ordinary liquids, demonstrating the importance of caging effect in this region. Figure 9 shows the diffusion and viscosity properties as functions of temperature along two isochores: 1.8 and 2.4 g/cm 3 . In the temperature ranges we consider here, the thermodynamic factors computed with average-atom model almost keep constants of 1.289 and 1.334 for the isochores of 1.8 and 2.4 g/cm 3 , respectively. For these two cases, both self-and mutual diffusion coefficients increase systematically with temperature, except that the nitrogen atoms diffuse with typical values of fluid (~1.0 × 10 4 cm 2 s −1 ) for temperature below 5000 K at 1.8 g/cm 3 . In addition, along the isochore of 2.4 g/cm 3 , the self-diffusion coefficients for H are of the same order of magnitude as the mutual-diffusion coefficient, while the self-diffusion coefficients for N are reduced by almost one order of magnitude. The viscosity behaves much different for these two densities: At 1.8 g/cm 3 the viscosity decreases with temperature, while at 2.4 g/cm 3 the viscosity presents a shallow minimum as temperature increased to 12000 K. That is because the system lye in strong coupling (50 ≤ Γ) and intermediate coupling (10 ≤ Γ ≤ 50) regions at 1.8 and 2.4 g/cm 3 , respectively. As mentioned above, at intermediate coupling the two mechanisms of coupling and kinetic effects contribute with a similar magnitude, and thus leads to a minimum. At stronger coupling the caging effects dominates, and the viscosity exhibits the characterization of liquids.
In summary, we explore the complex behavior of ammonia at extreme conditions. The PDF and BACF suggest that ammonia undergoes three structural changes along principal Hugoniot: pure molecules, mixture, and plasma. This quantitatively demonstrates that profound chemistry process occurs in shocked ammonia. Especially, we note that the shocked nitrogen atoms form the chains with the dissociation of ammonia. Through using a complex fitting functional including structures and multiple time scales, accurate diffusion and viscosity constants of strong correlated ammonia are derived. It is found that the diffusion and viscosity exhibit distinct features in these three structural regions, which result from the competitive mechanisms of coupling and kinetic effects.

Methods
The equilibrium configurations of the ions at different density-temperature states are determined by employing the quantum molecular dynamics (QMD) simulations, which are implemented in the Vienna ab initio simulation package (VASP) 18,19 . The equations of motion for active electrons are solved in the framework of finite-temperature density functional theory (FT-DFT) 20 Both of them are converged to be less than 0.5%, Thus the rcut = 1.2 bohrs (N)/1.1 bohrs (H) PAW potential is enough to reach convergence for the density ranges we considered in this manuscript. We also use the PBE 26 and LDA 27 functionals to study the influence of different exchange-correlation potentials on equation of state and structure, which are discussed in detail in the Results section.
In the simulations, the canonical (NVT) ensembles are set with the fixed number of particles, volume and temperature. The simulated density and temperature ranges from 0.8 to 2.6 g/cm 3 and temperature from 250 to 120000 K, respectively. We fix the plane-wave cutoff at 1000 eV and use only the Γ point to sample the Brillouin zone of a cubic supercell with 256 atoms. A sufficient number of bands are included to ensure the occupation of the highest band to a level of 10 −5 . The simulations last 20000~40000 steps, with a time step varying for different density-temperature ranges: 0.5~1.0 fs for temperatures lower than 6000 K, 0.2~0.3 fs for higher temperature. These equilibrated trajectories provide a consistent set of structure and transport properties of ammonia.
To characterize the structural change in ammonia, we not only calculate the pair distribution function (PDF) for each possible bondings of N-N, H-H, and N-H regularly, but also analyze the MD trajectories through the use of the bond auto-correlation function (BACF), which can shed light on the complex chemistry occurring in ammonia under extreme pressure-temperature conditions. The BACF β(t) is defined as where the bond vector B denotes a bit vector with value of 1 or 0 corresponding to the case for each bond of a certain type within or beyond a cutoff distance respectively. the direction of B vector is along the bond between two atoms. This correlation function describes the average duration of a bond between two atoms. In the present study, we use 1.7 Å for bonds between nitrogen atoms, 1.3 Å for bonds between nitrogen and hydrogen, and 0.9 Å for hydrogen-hydrogen bonds, which are the locations of the first minimum of pair distribution function for each species respectively. However, we note that the decay timescale (or the bond lifetime) is not strongly dependent on the chosen bond cutoff.
The transport properties, such as diffusivity and viscosity, are examined at different density-temperature conditions for ammonia. The self-diffusion coefficient of a particular ion species D α is computed from the integral of the VACF 28,29 , which is i i 0 with V i being the velocity of the ith particle of species α. We also derive the mutual-diffusion coefficient D αβ according to the Green-Kubo relation 30-32 Figure 5. The data from QMD run on ammonia at 1.9 g/cm 3 and 4791 K. (a) VACF for N in ammonia using Eq. (7) with i = 1; (b) Self-diffusion coefficient for N in ammonia using the fit to the analytic integral of Eq. (7); (c) VACF for H in ammonia using Eq. (7) with i = 2; (d) Self-diffusion coefficient for H in ammonia using the fit to the analytic integral of Eq. (7); (e) VACF for N-H using Eq. (9) with i = 3; (d) Mutual-diffusion coefficient for N-H using the fit to the analytic integral of Eq. (9). The green "x" are the numerical data, while the light-red dashed line is the fits using the respective equations. The solid black line is the example of a single exponential fit.
The thermodynamic factor Q accounts for nonideal mixing contributions to the Gibbs free energy with the expression of , . This factor reduces to unity for ideal gas mixtures, but in ionic mixtures goes to (〈Z*〉 + 〈(Z*) 2 〉)/(〈Z*〉 + 〈Z*〉 2 ) with Z* the average ionization degree 33 . x α and N α represent the concentration and particle number of species α, respectively. If the interspecies interactions remain small compared to the intraspecies ones, then the mutual diffusion coefficient takes a very simple form in terms of the concentrations and the self-diffusion coefficients: This relation does not correspond to a standard mixing rule since the self-diffusion coefficients arise from simulations on the full mixture, not on each pure species.
In the stress-stress approach, the viscosity is computed from the STACF of non-diagonal components of the stress tensor (P xy , P yz , P zx , (P xx − P yy )/2, and (P yy − P zz )/2) 34  It should be noted that to obtain similar statistical accuracy, the viscosity and mutual diffusion require much longer trajectories than the self-diffusion coefficient because the single-particle correlations are averaged over the particles and gain significant statistical improvement. To avoid the obstacle of computational cost, one can use empirical fits to the time dependence of the autocorrelation function, which can substantially shorten the length of the trajectory required. Then the long time behavior of the integrated functional form yields the desired transport property. The fitting functional form of simple exponential (or Gaussian) model is good enough to produce accurate results in the weakly coupling limit 36 . However, in strongly coupled regime, the autocorrelation function usually exhibits correlated and structured behavior, and thus the fitting functionals need more physics incorporated 37 . Especially for the dense mixture, there exist multiple restoring forces or multiple frequencies of collection motion acting on a given ion. Therefore, the fitting form with multiple times scales should be used to fit the autocorrelation function. As proposed in ref. 34 , the VACF could be given by with c i = a i /Σ i a i . The differences between the statistical error for Eqs (7) and (9) are that the disappearance of