Ionic self-diffusion coefficient and shear viscosity of high-Z materials in the hot dense regime

Yong Hou; Yang Jin; Ping Zhang; Dongdong Kang; Cheng Gao; Ronald Redmer; Jianmin Yuan

doi:10.1063/5.0024409

Abstract

High-Z materials exhibit a broad range of variation of the charge state in the hot dense regime, and so ionic structures become complex with increasing density and temperature owing to ionization. Taking high-Z uranium as example, we study its electronic and ionic structures in the hot dense regime by combining an average-atom model with the hypernetted chain approximation. The electronic structure is described by solving the Dirac equation, taking account of relativistic effects, including broadening of the energy levels, and the effect of other ions via correlation functions. On the basis of the electronic distribution around a nucleus, the ion pair potential is constructed using the modified Gordon–Kim model in the frame of temperature-dependent density functional theory. Because of the presence of ion–ion strong coupling, the bridge function is included in the hypernetted chain approximation, which is used to calculate the correlation functions. To take account of the influence on transport properties of the strong correlation of electrons with highly charged ions, we perform both classical and Langevin molecular dynamics simulations to determine ion self-diffusion coefficients and the shear viscosity, using the Green–Kubo relation and an ion–ion pair potential with good convergence. We show that the influence of electron–ion collisions on transport properties becomes more important as the free electron density increases owing to thermal ionization.

I. INTRODUCTION

The hot dense plasma region covers densities ranging from a few hundredths to hundreds of times solid density and temperatures from several tens of electronvolts (eV) up to the keV region. In this regime, the electrons of an atom are ionized owing to thermal and pressure ionization. Such conditions are relevant to inertial confinement fusion (ICF) capsules,¹ compact astrophysical objects such as white dwarfs,^2–4 and laboratory experiments reaching high energy densities.^5–8 To model, for example, the interaction of high-power lasers with solid targets to generate hot dense plasmas (HDPs), hydrodynamic simulations are performed, for which accurate physical properties such as the equation of state,^9–11 ionic transport properties,^12–16 electronic thermal and electrical conductivities,^17–23 opacity,^24–27 and stopping power^28,29 are needed. In particular, electronic and ionic transport properties are very important for modeling the generation of fast electron beams and energy deposition. However, the available experimental data are sparse, and the calculation of the properties of HDPs is also difficult owing to strong ion coupling and electron degeneracy in this regime. In particular, for high-Z materials such as uranium (U), which are often used for radiation shielding materials, the ion charge state can vary over a broad range, depending on density and temperature. The strong electric field provided by highly ionized states also makes this a situation of matter under extreme conditions. With an increase in the ionization degree, the greater fraction of free electrons will enhance Coulomb screening effects, but ion–ion interactions will also become stronger. The competition between ion–ion interactions and screening effects leads to an almost constant but large ion–ion coupling parameter in a wide range of temperatures and densities,^30–32 and so modeling HDPs poses a challenge for quantum-statistical approaches.

Density functional theory (DFT) was extended to finite temperatures by Mermin,³³ who took electron degeneracy into account by averaging the set of single-electron wave functions with a Fermi–Dirac distribution function at given temperatures for the element under consideration. By solving Newton’s equation to describe the ion motion, quantum molecular dynamics (QMD)^34–38 has become a powerful tool for computing the properties of warm and hot dense matter on the basis of pseudopotentials. Although some modified approaches have been developed to treat high-temperature plasmas as well,^39–41 these calculations are very expensive owing to the large number of single-electron wave functions that are required for the results to converge. Alternatively, the Thomas–Fermi approximation can be applied to replace the single-electron wave function, and orbital-free molecular dynamics (OFMD) simulations^42–49 have been developed for HDPs. If the ion–ion pair potential is obtained consistently from the electronic structure calculations, then the hypernetted chain (HNC) approximation can also be used to describe the ionic structure through solution of the Ornstein–Zernike (OZ) equation.^50–54 However, it is difficult to define the ionization state and its interactions, which are closely related to ionic transport properties,^13,55 especially for high-Z elements, where many charge states occur simultaneously. By considering the correlations of other ions and of free electrons, we have combined the average-atom (AA) model^56,57 with the HNC approximation in the AAHNC method to calculate the mean ionization degree, the x-ray Thomson scattering spectrum,⁵⁸ and nonequilibrium states⁵⁹ of aluminum in the warm dense matter regime, where the ionic structure is determined in the HNC approximation and the ion–ion pair potential is obtained from the electron distribution around the nucleus.⁶⁰

To obtain the ionic structures and transport properties of matter under extreme conditions, the bridge function has been incorporated into the AAHNC method, and this AAHNC+Bridge method has been applied to study the high-Z element U in the HDP region. Because of the wider variation of the ionization degree of U, there is a temperature region where the average charge is proportional to the square root of the temperature and where the ion coupling parameter Γ_ii remains almost constant with increasing temperature. This feature of an almost constant Γ_ii is similar to the behavior of the pair distribution function: the location of its first peak is also almost constant over a wide temperature range. Based on the effective pair potentials derived from the convergent electronic structure calculations, molecular dynamics simulations are performed to calculate the ionic self-diffusion coefficient and the shear viscosity from the autocorrelation function⁶¹ according to the Green–Kubo relation. Because of the large number of free electrons that are present in a hot dense U plasma, we consider the effect of dynamic electron–ion collisions on transport properties by a Langevin molecular dynamics simulation. By comparing the results with those of classical molecular dynamics simulations, we show that the impact of dynamic collisions on transport properties becomes more important with increasing ionization degree. Unless otherwise specified, we use atomic units throughout this paper, i.e., we take ℏ = m_e = e = 1, so that distances are measured in units of the bohr radius a_B = 0.529 177 Å = 5.291 77 × 10⁻¹¹ m, and energies in units of hartree, with 1 hartree = 27.211 386 eV.

II. THEORETICAL METHODS

A. Combining the average-atom model with the hypernetted chain approximation

For accurate calculation of electronic structure, the correlation functions of electrons and ions are needed. In addition, effective pair potentials between ions and electrons are based on the electron and ion density distributions, and so the electronic and ionic structures have to be calculated self-consistently.⁵⁸ The AA model is a statistical approach used to calculate the electronic structures of atoms and ions in an average ion sphere for HDPs. The Wigner–Seitz ion-sphere radius is determined from the ion number density n_i as $R_{W S} = {(3 / 4 π n_{i})}^{1 / 3}$ . The single-electron wave function is obtained by solving the Dirac equation: $\begin{align} \frac{d P_{n κ} (r)}{d r} + \frac{κ}{r} P_{n κ} (r) & = \frac{1}{c} [ϵ_{n κ} + c^{2} - V (r)] Q_{n κ} (r), \end{align}$ (1) $\begin{align} \frac{d Q_{n κ} (r)}{d r} - \frac{κ}{r} Q_{n κ} (r) & = - \frac{1}{c} [ϵ_{n κ} - c^{2} - V (r)] P_{n κ} (r), \end{align}$ (2)where P_nκ(r) and Q_nκ(r) are respectively the large and small components of the wave function of relativistic orbital nκ, and c is the speed of light. In the normal AA model, the central symmetric potential V(r) is often divided into three parts, namely, the interaction potential of the nucleus, the interaction potential of the electrons, and the exchange correlation, as shown in the first line of the following equation, which are self-consistently determined by the electron distribution in the ion sphere: $\begin{align} V (r) = & - \frac{Z}{r} + \int \frac{ρ_{b} (r^{'})}{|r - r^{'}|} d^{3} r^{'} + V_{x c} (ρ_{b} (r) + ρ_{e}^{0}) - V_{x c} (ρ_{e}^{0}) \\ - \frac{ρ_{e}^{0}}{β} \int C_{e e} (|r - r^{'}|) h_{i e} (r^{'}) d^{3} r^{'} \\ - \frac{ρ_{i}^{0}}{β} \int C_{i e} (|r - r^{'}|) h_{i i} (r^{'}) d^{3} r^{'} . \end{align}$ (3)To take account of the influence of free electrons and other ions, their statistical distributions are described by the correlation functions of ions and electrons, as shown in the second and third lines of the above equation.^58,62 In Eq. (3), β = 1/k_BT is the inverse temperature, and $ρ_{i}^{0}$ and $ρ_{e}^{0}$ are the uniform densities of ions and electrons, respectively. If the correlation functions are known, then Eqs. (1)–(3) can be solved through the self-consistent method to obtain the wave functions. In local thermodynamic equilibrium, to ensure that the electron density is the same at the boundaries of all ion spheres, the potential at the ion-sphere boundary is usually chosen as the energy reference point. In this way, we define electrons with energies larger than zero as free and those with negative energies as bound. The bound-state contribution to the density is $ρ_{b} (r) = \frac{1}{4 π r^{2}} \sum_{j} \int_{ϵ - Δ ϵ}^{ϵ + Δ ϵ} b_{j} (ϵ) [P_{j}^{2} (r) + Q_{j}^{2} (r)] d ϵ,$ (4)where pressure broadening effects are taken into account by energy level broadening,^56,57 Δϵ is determined by two different boundary conditions in solving the Dirac equations (1) and (2), and b_j(ϵ) is the density of the occupation number of state j expressed in terms of the Fermi–Dirac distribution. For hot plasmas, there are a large number of free electrons owing to thermal and pressure ionization, and we can use a semiclassical statistical approach, the Thomas–Fermi approximation, to compute the distribution of free electrons: $ρ_{f} (r) = \frac{1}{π^{2}} \int_{k_{0} (r)}^{\infty} \frac{k^{2} d k}{\exp {β [\sqrt{k^{2} c^{2} + c^{4}} - c^{2} - V (r) - μ]} + 1},$ (5)where $k_{0} (r) = {[2 V (r) c^{2} + V {(r)}^{2}]}^{1 / 2} / c$ and μ is the chemical potential, which is determined by the condition of electrical neutrality in the ion sphere.

For calculation of the central symmetric potential V(r), the electron–electron, electron–ion, and ion–ion correlation functions of ions and electrons are needed. Here, the HNC approximation with the Ornstein–Zernike closure relation is used to calculate electron–electron and ion–ion correlation functions based on the pair potential. In principle, the ionic interactions depend on the spatial distribution of nuclei and the electron density around the nuclei. In HDPs, in particular for high-Z U, owing to the static Coulomb screening effect of free electrons, only interactions of nearest-neighbor ions are important. Therefore, we consider only the pair potential and apply the modified Gordon–Kim (GK) model;⁶⁰ details of the calculation method are given in the Appendix. Considering the ionic strong coupling, a bridge function is included in the HNC approximation to calculate the correlations: $h_{i i} (r) = \exp [- β V_{i i} (r) + h_{i i} (r) - C_{i i} (r) + B_{i i} (r)] - 1,$ (6)where V_ii represents the interactions of ion i and B_ii(r) is the bridge function. B_ii(r) has a well-defined diagrammatic representation, but it is impossible to compute this infinite series, and so different approximations are applied to describe the bridge function. Here, we use the bridge function of a one-component plasma (OCP)^63,64 to modify the HNC approximation (this is the HNC+Bridge model): $B_{i i} (r) = - α Γ^{β} \exp (- \frac{b_{1}}{b_{0}} r^{2}),$ (7)where $\begin{aligned} b_{0} & = 0.258 - 0.0612 \ln Γ + 0.0123 {(\ln Γ)}^{2} - 1 / Γ, \\ b_{1} & = 0.0269 + 0.0318 \ln Γ + 0.008 14 {(\ln Γ)}^{2}, \end{aligned}$ (8)and α and β are constants determined from the results of molecular dynamics simulations. Because the pair distribution function (PDF) reflects the interaction of ions,⁶⁵ we calculate the PDF at a temperature of 100 eV and different densities (1.893 g/cm³, 9.465 g/cm³, 18.93 g/cm³, and 94.65 g/cm³) using the AAHNC+Bridge and Langevin molecular dynamics (LMD) methods to obtain the values of α and β. The best-fitting results are α = 0.026 and β = 1.3, which are then used for the bridge function (7). All the results are shown in Fig. 1. The difference between the AAHNC and AAHNC+Bridge results increases with increasing density, but the agreement is still reasonable. We conclude that the effects of the bridge function become more pronounced in the regime of strong ionic coupling, as expected.

Figure 1.Pair distribution functions as calculated with the AAHNC (red dashed lines), AAHNC+B (black solid lines), and LMD (blue dot-dashed lines with circles) methods as functions of the ion–ion distance for U at a temperature of 100 eV and densities of 1.893 g/cm³, 9.465 g/cm³, 18.93 g/cm³, and 94.65 g/cm³.

For the free electrons, we use the Deutsch pair potential^50,66 $\begin{align} V_{e e}^{D e u t s c h} (r) = & \frac{1}{r} [1 - \exp (- \frac{r}{λ_{e e}})] \\ + k_{B} T \ln 2 \exp [- \frac{\ln 2}{π} {(\frac{r}{λ_{e e}})}^{2}], \end{align}$ (9)with $λ_{e e} = 1 / \sqrt{k_{B} T}$ . The first term in Eq. (9) arises from quantum effects such as the uncertainty principle. The second term is due to the exchange interaction, which yields a repulsion if electrons have parallel spin; here, the spin-average result is used. The Deutsch pair potential gives Coulomb-like behavior at long distances and a finite value at r = 0. It has been proved to give a reasonable description for high temperatures k_BT ≥ 0.5 Hartree as considered here and λ_ee ≪ r_s, where r_s is the effective electron radius.⁶⁷ Because the electron density distribution around the nucleus has been derived from the self-consistent calculation in the ion sphere, the ion–electron pair correlation function can be defined via the excess free electron density: $h_{i e} (r) = \frac{ρ_{f} (r)}{ρ_{e}^{0}} - 1,$ (10)where $ρ_{e}^{0}$ is the uniform electron density given by $ρ_{e}^{0} = ρ (r_{b})$ .

B. Molecular dynamics simulation and calculation of transport properties

Based on the results of the self-consistent calculations of the electronic and ionic structures using the AAHNC+Bridge method, we have obtained an ion pair potential for hot dense U plasmas that exhibits good convergence. Molecular dynamics simulations can now be performed to compute ionic transport properties. Note that all calculations are based on the Born–Oppenheimer approximation, with electron–ion screening included in the ion–ion interactions. However, owing to ionic motion, the screening cloud around the ion will deform, which leads to an additional electron–ion force known as the Debye–Onsager relaxation effect.⁶⁸ The force produced by electron–ion collisions hinders the motion of the ions, which has not been treated in the simulations. In particular, with increasing temperature and density, more and more electrons become free owing to thermal and pressure ionization, and thus there is a higher frequency of electron–ion collisions. An efficient way to account for nonadiabatic effects is to use LMD^39,40,59 simulations, which consider the motion of ions in a dense electron gas as the motion of Brownian particles, with random electronic collisions leading to a frictional force: $m_{i} \frac{d^{2} r_{i} (t)}{d t^{2}} = F - γ m_{i} \frac{d r_{i} (t)}{d t} + N_{i},$ (11)where m_i, r, and F are the ionic mass, position, and force, respectively. The ion–ion force is derived from the convergent pair potential. N_i is the Gaussian random noise corresponding to the Langevin friction coefficient γ. Quantifying this extra friction in Langevin mechanics is a complex and difficult procedure. For light elements, some efforts have been made to compare LMD results with those obtained using time-dependent DFT (TDDFT), and it has been found that the friction coefficient has a significant influence on the dynamic transport properties in LMD simulations.⁶⁹ In hot dense high-Z plasmas, there are an increased number of free electrons, and so the influence of electron–ion collisions on dynamical properties become more important. To consider the effects of the friction on ionic transport properties, we estimate γ using the Rayleigh model:⁷⁰ $γ = 2 π \frac{Z^{*}}{R_{W S} m_{i}} \sqrt{k_{B} T},$ (12)where Z* is the average ionization degree. Although the Rayleigh model considers only binary elastic collisions and assumes that these are short-ranged with a negligible collision time in a rarefied gas, the model can be extended to study the effects of nonequilibrium fluctuations and of dynamic electron–ion collisions, and the predictions of this more general model can be compared with the results of the original model. As examples here, ionic transport properties are calculated by performing both Langevin and classical molecular dynamics simulations.

According to the Green–Kubo relation, the ionic self-diffusion coefficient D is obtained from the velocity autocorrelation functions:^13,61 $D = \frac{1}{3 N} \int_{0}^{\infty} 〈\sum_{i = 1}^{N} v_{i} (t) \cdot v_{i} (0)〉 d t,$ (13)where v_i(t) is the velocity of ion i at time t, N is the number of simulation particles, and $〈\cdot〉$ denotes the statistical ensemble average. The shear viscosity is computed from the autocorrelation function of the off-diagonal components of the stress tensor:⁶¹ $η = \frac{V}{k_{B} T} \int_{0}^{\infty} 〈P_{12} (0) P_{12} (t)〉 d t,$ (14)where V is the volume of the simulation box. P₁₂(t) represents the five independent off-diagonal components P_xy(t), P_yz(t), P_zx(t), [P_xx(t) − P_yy(t)]/2, and [P_yy(t) − P_zz(t)]/2 at time t, which are calculated as $P_{12} (t) = \frac{1}{V} [\sum_{i} m_{i} v_{i 1} (t) v_{i 2} (t) + \sum_{i} \sum_{j > i} r_{i j 1} (t) F_{i j 2} (t)] .$ (15)

The Stokes–Einstein relation provides a connection between the self-diffusion coefficient and the shear viscosity through the expression $F_{S E} [D, η] = \frac{D η}{k_{B} T n_{i}^{1 / 3}} = C_{S E},$ (16)where C_SE is a constant, often taken in the range from 1/6π to 1/4π, depending on the limits of the slip coefficient from infinity (stick) to zero (slip), respectively.

III. RESULTS AND DISCUSSION

A. Ionic structures based on the AAHNC+Bridge model

For the high-Z element U, there is a wide region of ionic charge states in the HDP regime. We first calculate the average charge of U using the AAHNC+Bridge model and compare the results with those of the INFERNO model¹⁶ as functions of temperature for the solid density ρ₀ = 18.93 g/cm³ and for 5 × ρ₀ = 94.65 g/cm³, as shown in Fig. 2. In a very wide temperature range, from 50 eV to 5000 eV, the resulting average ionization degree is in good agreement with the INFERNO model for these two densities, although small differences occur for temperature in the middle of that range. The INFERNO model also calculates the electronic structure by solving the Dirac equation in a self-consistent-field approach and assumes an all-electron distribution in the ion sphere that is determined by the plasma density. The correlations of electrons and other ions have been included in our calculations of the electron structures, and so slight differences occur in the temperature range from 100 eV to 1000 eV. In Fig. 2, we also give average charges at other densities: 0.1 × ρ₀ = 1.893 g/cm³ and 0.5 × ρ₀ = 9.465 g/cm³. At lower temperatures, the resulting average ionization degree for 5 × ρ₀ is higher than for the other densities because of pressure ionization. Thermal ionization becomes important at high temperatures and makes the average ionization degree at lower densities stronger. In the inset of Fig. 2, the results are shown on a log–log plot, along with the curve of the square root of temperature, $\sqrt{k_{B} T}$ . Except for the results for 5 × ρ₀, the behavior of the average charge is very similar to the $\sqrt{k_{B} T}$ behavior for temperatures lower than 600 eV.

Figure 2.Average charge of U as a function of temperature at solid density ρ₀ = 18.93 g/cm³ (circles), 0.1 × ρ₀ = 1.893 g/cm³ (triangles up), 0.5 × ρ₀ = 9.465 g/cm³ (triangles down), and 5 × ρ₀ = 94.65 g/cm³ (squares). Solid lines with different symbols represent the results for the different densities calculated using the AAHNC+Bridge model. The blue dashed lines with circles (18.93 g/cm³) and squares (94.65 g/cm³) are calculated using the INFERNO model.¹⁶ The inset shows a log–log plot along with the $\sqrt{k_{B} T}$ behavior (brown dot-dashed line).

As shown in Fig. 2, the average ionization degree becomes larger and larger with increasing electron temperature, and so the ion–ion interactions become stronger. Simultaneously, the greater fraction of ionized electrons strengthens the interionic shielding. To study the variations in the ion–ion interactions, the ionic coupling parameter Γ_ii is computed as $Γ_{i i} = \frac{Z^{* 2}}{R_{W S} k_{B} T} .$ (17)Normally, Γ_ii decreases with isochoric heating according to Eq. (17). However, a region of almost constant Γ_ii has been reported in Refs. 30–32 for tungsten plasmas at densities of 19.2 g/cm³, 40 g/cm³, and 100 g/cm³. The average charge is proportional to the square root of the temperature, $Z^{*} \propto \sqrt{k_{B} T}$ , and so the influence of temperature cancels in Eq. (17). This feature becomes more obvious for a hot U plasma because of the wider variation of the average ionization degree. In Fig. 3, we show the coupling parameter Γ_ii as a function of temperature at densities from 0.1 to 5.0 times solid density. We also find a region with an almost constant coupling parameter. In particular, for the solid density ρ₀, the results do not change up to very high temperatures T ≥ 1000 eV. For the lowest density, 0.1 × ρ₀, Γ_ii first increases with temperature from about 20 eV owing to the rapid increase in the average charge, as shown in Fig. 2, i.e., Z* increases more rapidly than $\sqrt{k_{B} T}$ . At temperatures T ≥ 600 eV, the rate of increase of Z* becomes slower than that of $\sqrt{k_{B} T}$ , and Γ_ii begins to decrease. However, for the high density of 5 × ρ₀, Γ_ii decreases with increasing temperature, because the ionization degree increases only slightly, as shown in Fig. 2. To compare with the results of Arnault et al.³¹ and Clérouin et al.,³² we also calculate the reduced density ρ* = ρ/AZ, where A is the mass of the element. The reduced density of solid-density U, $ρ_{0}^{*}$ , is 0.000 87 mol/cm³, and thus $5 \times ρ_{0}^{*} = 0.0043$ mol/cm³. According to Refs. 31 and 32, a critical density $ρ_{0}^{*} = 0.0045$ mol/cm³ exists above which the coupling parameter decreases monotonically with increasing temperature, while there is a plateau below the critical density. The present results are in agreement with those earlier findings. These results also show impressively that the coupling parameter Γ_ii of high-Z plasmas does not depend linearly on temperature but also on details of the ionization process that depend on temperature and density.

Figure 3.Ionic coupling parameter Γ_ii for a hot dense U plasma as a function of temperature for densities in the range (0.1–5.0) × ρ₀.

From Fig. 3, we also find that the ionic coupling parameter Γ_ii is still very large owing to the highly ionized charge of U in hot plasmas, even at low densities 0.1 × ρ₀, which indicates that ionic correlation effects are very strong. This effect can be illustrated by the ion–ion PDF derived from the AAHNC+Bridge model, which is shown in Fig. 4 over the whole temperature range from 40 eV to 4000 eV for (a) 0.1 × ρ₀, (b) ρ₀, and (c) 5 × ρ₀. From the results for solid density ρ₀ in Fig. 4(b), it can be seen that from 200 eV to 1500 eV, there is almost no change in the PDF, and the height of the first peak does not decrease with increasing temperature up to 1500 eV. This behavior is very similar to that of Γ_ii, which is almost constant in that temperature range at solid density. This correlation between the ionic coupling parameter and the PDF is even more obvious for the density 0.1 × ρ₀, as can be seen in Fig. 4(a). As the temperature increases from 40 eV to 600 eV, the first peak of the PDF becomes sharper and sharper. For T ≥ 800 eV, the first peak begins to decrease, which is similar to the trend exhibited by the ionic coupling parameter as shown in Fig. 3: first increasing and then dropping down at 800 eV. For dilute plasmas, for example at 0.1 × ρ₀, the ionic ionization increases faster than linearly with temperature. Simultaneously, the screening effects of free electrons increase slowly, but the ion–ion interaction is very strong owing to the rate of increase of the ionic charge. For solid density, the competition between ionic interactions and electron screening tends to balance, which leads to an almost constant Γ_ii for a wide range of temperature. At higher density, 5 × ρ₀, owing to the limited space between ions, the static screening effects of free electrons become more important, which causes Γ_ii and the first peak of the PDF to fall for all temperatures, which is in good agreement with the conclusions of Refs. 31 and 32.

Figure 4.Pair distribution functions as derived from the AAHNC+Bridge method are shown as functions of the ion–ion distance for U at densities 1.893 g/cm³ (left), 18.93 g/cm³ (middle), and 94.65 g/cm³ (right) for different temperatures.

B. Transport properties using Langevin molecular dynamics simulations

On the basis of the pair potentials obtained from the self-consistent calculations within the AAHNC+Bridge method, molecular dynamics simulations can be performed to compute ionic transport properties. The discussion of the ionic coupling parameter and the pair distribution function of the high-Z element U in Sec. III A has shown that electron screening effects are very important with increasing ionization degree at high densities, and these have been included in the calculation of the pair potentials. The Born–Oppenheimer approximation has been used to separate the interaction between ions and electrons. However, with increasing free electron density, dynamic random collisions between electrons and ions become important, although these are neglected in classical molecular dynamics (CMD) simulations. Therefore, we consider the effect of dynamic electron–ion collisions by introducing a friction coefficient γ in the LMD simulations. In both CMD and LMD simulations, we use 13 500 ions in a cubic box with the NVT ensemble (constant N, V, and T), and apply periodic boundary conditions to simulate the positions and velocities of particles at the boundary. The time step of ionic motion, δt, is determined by the plasma temperature and density: $δ t \sim R_{W S} / T_{i}^{1 / 2}$ . All ions are initially are put in an ideal face-centered cubic lattice, and their velocities are set according to a Maxwellian velocity distribution. Based on the molecular dynamics simulations, we calculate the ionic self-diffusion coefficients and the shear viscosity via the velocity autocorrelation function according to the Green–Kubo relation.

We have used this approach previously to treat a two-temperature aluminum plasma⁵⁹ and have discussed the influence of the charge state on the ionic dynamic structure and the transport properties at a given ion temperature. For high-Z plasmas, there are more complex charge states and ionic structures with increasing temperature, as shown in Figs. 2 and 3. At the same time, the increasing fraction of free electrons makes the effect of electron–ion dynamic collisions on the ionic motion more important. Here, we discuss the electron–ion collision effect via the friction coefficient, which is determined by the average ionization charge and the electron temperature as in Eq. (12); it increases more rapidly than the charge state with temperature. If the Born–Oppenheimer approximation is assumed when calculating the electron structure, the nonadiabatic dynamic effect is neglected. However, in dense plasmas, with increasing temperature, both pressure and thermal ionization produce more and more free electrons, and the ions move in a free-electron gas in a manner similar to Brownian motion. Therefore, nonadiabatic dynamic effects will influence the transport properties,⁶⁹ and this becomes more obvious in the results for hot dense U plasmas, as can be seen in Fig. 5, where the self-diffusion coefficient and shear viscosity are shown as functions of temperature at three densities of 1.893 g/cm³, 18.93 g/cm³, and 94.65 g/cm³. All the results are obtained from LMD and CMD simulations with an effective pair potential derived from the AAHNC+Bridge model. We find that the difference between the LMD and CMD results becomes bigger and bigger with increasing temperature owing to the inclusion of the friction coefficient in the LMD simulations. More electrons of the U atom become ionized at high temperatures, and the effect of electron–ion dynamic collisions becomes more and more important with this increasing ionization degree. This effect is most obvious at the lowest density of 1.893 g/cm³, because of the higher ionization degree compared with that at higher densities (see Fig. 2). We also find that the self-diffusion coefficients and shear viscosities from the LMD simulations are all smaller than those from the CMD simulations. Electron–ion dynamic collisions increase the effective collision cross section and weaken the interaction between ions.

Figure 5.Self-diffusion coefficient D (a) and shear viscosity η (b) of U as functions of temperature at different densities: 1.893 g/cm³ (dashed lines), 18.93 g/cm³ (solid lines), and 94.65 g/cm³ (dot-dashed lines). The CMD (black circles) and LMD (red squares) simulations are based on the effective pair potential from the AAHNC+Bridge calculations. For comparison, the results of OFMD (blue triangles up) simulations and the R–OCP model (green stars)¹⁶ are also shown.

In Fig. 5, we compare our results with those of OFMD and the R–OCP model. Here, R–OCP refers to the classical one-component plasma model with effective charge determined by the INFERNO model.¹⁶ The LMD and CMD results are in good agreement with those of the R–OCP model, because their average charges agree very well (see Fig. 2). However, only the shear viscosity from OFMD agrees with the results of the other models, and big differences occur for the self-diffusion coefficient. In the OFMD model, the Thomas–Fermi approximation is applied to describe the electronic structure, which gives different ionic charges compared with the present model and R–OCP, which influences the self-diffusion coefficient strongly.^16,55

We also check the Stokes–Einstein relation according to Eq. (16) against the data for the self-diffusion coefficient and the shear viscosity. In Fig. 6, the quantity F_SE is shown as a function of temperature at densities 1.893 g/cm³, 9.465 g/cm³, 18.93 g/cm³, and 94.65 g/cm³. For the highest densities and with increasing temperature, F_SE is almost constant and independent of temperature and density except at the highest temperatures. The values computed from the LMD results, represented by filled circles, are larger than the stick value, 1/6π, and smaller than that for slip, 1/4π. This means that the Stokes–Einstein relation is also valid for hot dense high-Z plasmas, and that ions move in the dense free-electron gas very similarly to Brownian motion. However, the values of the ionic transport coefficients derived from CMD simulations, represented by open circles, which are similar to those from the INFERNO-OCP model in Ref. 16, are higher than those from the LMD simulations and always larger than the slip values. Therefore, it is necessary to consider the influence of the dense free-electron gas on the ionic motion. For the lowest density, 1.893 g/cm³, the values of F_SE comply with the Stokes–Einstein relation only in the intermediate temperature region. The other values deviate from the stick and slip values, and so these plasmas cannot be considered as simple fluids, owing to the presence of pressure and thermal ionization, which makes the effective ion-sphere radius very large.

Figure 6.Stokes–Einstein relation F_SE as a function of temperature calculated from the diffusion coefficients and the shear viscosities at densities 1.896 g/cm³ (orange), 9.465 g/cm³ (black), 18.93 g/cm³ (red), and 94.65 g/cm³ (blue). The filled circles represent the results of LMD simulations and the open symbols those of CMD simulations. The gray dashed lines show the constant values of the Stokes–Einstein relation for stick (1/6π) and slip (1/4π) boundary conditions, respectively.

IV. CONCLUSIONS

We have extended the AAHNC model to describe high-Z U plasmas by including the bridge function for the ionic structure (AAHNC+Bridge). The results for the average ionization degree are in good agreement with those of the INFERNO model in a wide temperature and density region. Because the average ionization degree is almost proportional to the square root of the temperature, a plateau feature of the coupling parameter Γ_ii is observed, consistent with previous predictions.^31,32 We have presented results for the ion–ion pair distribution function and have found that the first peak does not change very much over a wide range of temperature. Even for the lowest densities, 1.893 g/cm³ and 9.465 g/cm³, the first peak of the PDF becomes sharper and sharper, which corresponds to an increase in the ionic coupling parameter with increasing temperature. The increased average ionization degree produces more and more free electrons in the plasma. Besides the screening effect of free electrons, dynamic electron–ion collisions are accounted for by a friction force in LMD simulations, where the force is estimated by the Rayleigh model. We have compared the ionic transport properties, the self-diffusion coefficient, and the shear viscosity as calculated by LMD and CMD simulations. We have found that although only short-ranged elastic collisions are considered in the Rayleigh model, the increased number of free electrons strengthens the electron–ion collisions, and their influence on transport properties becomes more and more apparent with increasing temperature.

APPENDIX: MODEL OF THE ION–ION PAIR POTENTIAL

The GK model^71,72 was initially proposed to calculate the interatomic potential in ionic crystals and systems of noble-gas atoms. The total electron density of two interacting atoms is taken as the sum of the two separate atomic densities, which are obtained in the single-center force-field approximation. The total energy of the system includes the Coulomb interaction between all the charges, the electronic kinetic energy, the electron exchange effects, and the correlations: $V (R) = V_{C o u l} (R) + V_{k} (R) + V_{e} (R) + V_{c} (R),$ where R is the distance between the two nuclei, V_Coul(R) is their static Coulomb interaction, V_k(R) is the kinetic energy, and V_e(R) and V_c(R) are the exchange and correlation energies, respectively. The Coulomb energy is calculated by using the electron densities, while the other three terms are obtained from the free-electron-gas expressions according to the density. The interatomic interaction is the total energy minus the energy of the two separate atoms.

In a hot dense plasma, the density of free electrons is almost uniformly distributed and does not change very much when the two ions come closer. Therefore, we can divide the electron density into two parts: a uniformly distributed free-electron sea ρ(r_b) with a density equal to that of the electrons at the ion-sphere boundary, and a density of quasilocalized electrons $ρ_{i}^{l o c} (r)$ representing the dramatic spatial variations of the electronic distribution around the nucleus. The latter quantity is the total density minus the uniformly distributed free-electron density, which is not distorted when the two ions come closer. The total density is $ρ = ρ_{A}^{l o c} + ρ_{B}^{l o c} + ρ (r_{b})$ , where A and B represent the different ions. To keep the electron neutral, the volume of the uniformly distributed free electrons becomes the two truncated spherical volumes when the two ions come closer.⁶⁰

References

[1] T. C. Sangster, V. N. Goncharov, S. P. Regan et al. The National direct-drive inertial confinement fusion program. Nucl. Fusion, 59, 032007(2019).

[2] R. F. Smith, T. S. Duffy. Ultra-high pressure dynamic compression of geological materials. Front. Earth Sci., 7, 23(2019).

[3] R. W. Lemke, M. P. Desjarlais, M. D. Knudson et al. Probing the interiors of the ice giants: Shock compression of water to 700 GPa and 3.8 g/cm³. Phys. Rev. Lett., 108, 091102(2012).

[4] D. C. Swift, J. H. Eggert, D. G. Hicks et al. Mass-radius relationships for exoplanets. Astrophys. J., 744, 59(2012).

[5] W. M. Zhou, H. B. Cai, F. Zhang et al. Enhanced energy coupling for indirect-drive fast-ignition fusion targets. Nat. Phys., 16, 810(2020).

[6] D. D. Ryutov, R. P. Drake, B. A. Remington. Experimental astrophysics with high power lasers and Z pinches. Rev. Mod. Phys., 78, 755(2006).

[7] D. Batani, M. Koenig, J. L. Miquel et al. Development of the PETawatt Aquitaine laser system and new perspectives in physics. Phys. Scr., T161, 014016(2014).

[8] W. J. Hogan, E. M. Campbell. The national ignition facility-applications for inertial fusion energy and high-energy-density science. Plasma Phys. Control. Fusion, 41, B39(1999).

[9] P. Arnault, S. X. Hu, J. A. Gaffney et al. A review of equation-of-state models for inertial confinement fusion materials. High Energy Density Phys., 28, 7-24(2018).

[10] J. Nilsen, M. E. Martin, A. L. Kritcher et al. Understanding the effects of radiative preheat and self-emission from shock heating on equation of state measurement at 100 s of Mbar using spherically converging shock waves in a NIF hohlraum. Matter Radiat. Extremes, 5, 018401(2020).

[11] Q. Zhang, H. Liu, H. Song et al. Validation for equation of state in wide regime: Copper as prototype. Matter Radiat. Extremes, 1, 123(2016).

[12] S. D. Bergeson, C. L. Ellison, S. D. Baalrud et al. Exploring the crossover between high-energy-density plasma and ultracold neutral plasma physics. Phys. Plasmas, 26, 100501(2019).

[13] J. Dai, D. Kang, Y. Hou et al. Equations of state and transport properties of mixtures in the warm dense regime. Phys. Plasmas, 22, 022711(2015).

[14] E. R. Meyer, C. Ticknor, J. D. Kress, L. A. Collins. Effect of correlation on viscosity and diffusion in molecular-dynamics simulations. Phys. Rev. E, 90, 043101(2014).

[15] J. Daligault, S. D. Baalrud. Effective potential theory for transport coefficients across coupling regimes. Phys. Rev. Lett., 110, 235001(2013).

[16] J. D. Kress, D. P. Kilcrease, J. S. Cohen et al. Orbital-free molecular dynamics simulations of transport properties in dense-plasma uranium. High Energy Density Phys., 7, 155-160(2011).

[17] B. B. L. Witte, E. Galtier, L. B. Fletcher et al. Warm dense matter demonstrating non-Drude conductivity from observations of nonlinear plasmon damping. Phys. Rev. Lett., 118, 225001(2017).

[18] P. Sperling, M. French, B. B. L. Witte et al. Observations of non-linear plasmon damping in dense plasmas. Phys. Plasmas, 25, 056901(2018).

[19] M. W. C. Dharma-wardana, D. D. Klug. Isochoric, isobaric, and ultrafast conductivities of aluminum, lithium, and carbon in the warm dense matter regime. Phys. Rev. E, 96, 053206(2017).

[20] C. E. Starrett. Kubo-Greenwood approach to conductivity in dense plasmas with average atom models. High Energy Density Phys., 19, 58-64(2016).

[21] W. R. Johnson, J. Nilsen. Average-atom treatment of relaxation time in x-ray Thomson scattering from warm dense matter. Phys. Rev. E, 93, 033205(2016).

[22] J. C. Pain, G. Dejonghe. Electrical resistivity in warm dense plasmas beyond the average-atom model. Contrib. Plasma Phys., 50, 39-45(2010).

[23] M. Schörner, M. Bethkenhagen, B. B. L. Witte et al. Carbon ionization at gigabar pressures: An ab initio perspective on astrophysical high-density plasmas. Phys. Rev. Res., 2, 023260(2020).

[24] T. Nagayama, J. E. Bailey, G. P. Loisel et al. A higher-than-predicted measurement of iron opacity at solar interior temperatures. Nature, 517, 56-59(2015).

[25] J. E. Bailey, G. P. Loisel, T. Nagayama et al. Systematic study of L-Shell opacity at stellar interior temperatures. Phys. Rev. Lett., 122, 235001(2019).

[26] P. Liu, Y. Hou, C. Gao et al. Transient space localization of electrons ejected from continuum atomic processes in hot dense plasma. Commun. Phys., 1, 95(2018).

[27] F. B. Rosmej, O. Renner. Challenges of x-ray spectroscopy in investigations of matter under extreme conditions. Matter Radiat. Extremes, 4, 024201(2019).

[28] J. A. Frenje, P. E. Grabowski, A. B. Zylstra et al. Measurement of charged-particle stopping in warm dense plasma. Phys. Rev. Lett., 114, 211002(2015).

[29] C. Deutsch. Correlated ion stopping in dense plasmas. Matter Radiat. Extremes, 4, 034201(2019).

[30] G. Robert, J. Clérouin, P. Arnault. Behavior of the coupling parameter under isochoric heating in a high-Z plasma. Phys. Rev. E., 87, 061101(R)(2013).

[31] G. Robert, P. Arnault, J. Clérouin. Thomas-Fermi Z-scaling laws and coupling stabilization for plasmas. Phys. Rev. E., 88, 063106(2013).

[32] J. D. Kress, G. Robert, J. Clérouin, L. A. Collins, P. Arnault. Self-organization in dense plasmas: The Gamma-Plateau. Contrib. Plasma Phys., 55, 159-163(2015).

[33] N. D. Mermin. Thermal properties of the inhomogeneous electron gas. Phys. Rev., 137, A1441(1965).

[34] V. Recoules, S. X. Hu, V. V. Karasiev et al. Interspecies radiative transition in warm and superdense plasma mixtures. Nat. Commun., 11, 1989(2020).

[35] R. Redmer, H. R. Rüter. Ab initio simulations for the ion-ion structure factor of warm dense aluminum. Phys. Rev. Lett., 112, 145007(2014).

[36] F. Lambert, S. Mazevet, F. Bottin et al. Ab initio molecular dynamics simulations of dense boron plasmas up to the semiclassical Thomas-Fermi regime. Phys. Rev. E, 75, 056404(2007).

[37] M. Krack, F. R. Mohamed, T. D. Kühne, M. Parrinello. Efficient and accurate Car-Parrinello-like approach to Born-Oppenheimer molecular dynamics. Phys. Rev. Lett., 98, 066401(2007).

[38] M. P. Desjarlais. Density-functional calculations of the liquid deuterium Hugoniot, reshock, and reverberation timing. Phys. Rev. B, 68, 064204(2003).

[39] Y. Hou, J. Dai, J. Dai, J. Yuan, Y. Hou, J. Yuan. Quantum Langevin molecular dynamic determination of the solar-interior equation of state. Astrophys. J., 721, 1158(2010).

[40] J. Dai, D. Kang, Y. Hou et al. Structure, equation of state, diffusion and viscosity of warm dense Fe under the conditions of a giant planet core. New J. Phys., 15, 045003(2013).

[41] W. Kang, S. Zhang, H. Wang et al. Extended application of Kohn-Sham first-principles molecular dynamics method with plane wave approximation at high energy-From cold materials to hot dense plasmas. Phys. Plasmas, 23, 042707(2016).

[42] G. Zérah, F. Lambert, J. Clérouin. Very-high-temperature molecular dynamics. Phys. Rev. E, 73, 016403(2006).

[43] V. Recoules, F. Lambert. Plastic ablator and hydrodynamic instabilities: A first-principles set of microscopic coefficients. Phys. Rev. E, 86, 026405(2012).

[44] L. Burakovsky, C. Ticknor, J. D. Kress et al. Transport properties of lithium hydride at extreme conditions from orbital-free molecular dynamics. Phys. Rev. E, 87, 023104(2013).

[45] J. Daligault, T. Sjostrom. Fast and accurate quantum molecular dynamics of dense plasmas across temperature regimes. Phys. Rev. Lett., 113, 155006(2014).

[46] T. Sjostrom, J. Daligault. Ionic and electronic transport properties in dense plasmas by orbital-free density functional theory. Phys. Rev. E, 92, 063304(2015).

[47] C. E. Starrett. Thomas-Fermi simulations of dense plasmas without pseudopotentials. Phys. Rev. E, 96, 013206(2017).

[48] C. E. Starrett, D. Saumon. Equation of state of dense plasmas with pseudoatom molecular dynamics. Phys. Rev. E, 93, 063206(2016).

[49] C. E. Starrett, D. Saumon, J. Daligault. Pseudoatom molecular dynamics. Phys. Rev. E, 91, 013104(2015).

[50] R. Redmer, R. Bredow, Th. Bornath, W.-D. Kraeft. Hypernetted chain calculation fo multi-component and non-equilibrium plasmas. Contrib. Plasma Phys., 53, 276(2013).

[51] K. Wünsch, M. Schlanges, D. O. Gericke, P. Hilse. Structure of strongly coupled multicomponent plasmas. Phys. Rev. E, 77, 056404(2008).

[52] Th. Bornath, W.-D. Kraeft, V. Schwarz et al. Hypernetted chain calculations for two component plasmas. Contrib. Plasma Phys., 47, 324(2007).

[53] V. Bezkrovniy, D. Kremp, W.-D. Kraeft, M. Schlanges. Reaction ensemble Monte Carlo technique and hypernetted chain approximation study of dense hydrogen. Phys. Rev. E, 69, 061204(2004).

[54] J. Hansen, M. Baus. Statistical mechanics of simple Coulomb systems. Phys. Rep., 59, 1(1980).

[55] Y. Fu, D. Kang, Y. Hou et al. Multi-charge-state molecular dynamics and self-diffusion coefficient in the warm dense matter regime. Phys. Plasmas, 25, 012701(2018).

[56] Y. Hou, J. Yuan, F. Jin. Influence of the electronic energy level broadening on the ionization of atoms in hot and dense plasmas: An average atom model demonstration. Phys. Plasmas, 13, 093301(2006).

[57] Y. Hou, J. Yuan, F. Jin. Energy level broadening effect on the equation of state of hot dense Al and Au plasma. J. Phys.: Condens. Matter, 19, 425204(2007).

[58] J. Yuan, R. Bredow, R. Redmer, Y. Hou. Average-atom model combined with the hypernetted chain approximation applied to warm dense matter. Phys. Rev. E, 91, 033114(2015).

[59] R. Bredow, Y. Fu, Y. Hou et al. Average-atom model for two-temperature states and ionic transport properties of aluminum in the warm dense matter regime. High Energy Density Phys., 22, 21-26(2017).

[60] J. Yuan, Y. Hou. Alternative ion-ion pair-potential model applied to molecular dynamics simulations of hot and dense plasmas: Al and Fe as examples. Phys. Rev. E, 79, 016402(2009).

[61] D. J. Tildesley, M. P. Allen. Computer Simulation of Liquids(1987).

[62] C. E. Starrett, J. Clérouin, D. Saumon, J. D. Kress. The quantum hypernetted chain model of warm dense matter. High Energy Density Phys., 8, 150(2012).

[63] S. Ichimaru, S. Ogata, H. Iyetomi. Bridge functions and improvement on the hypernetted-chain approximation for classical one-component plasmas. Phys. Rev. A, 46, 1051(1992).

[64] M. S. Murillo, A. Diaw. A dynamic density functional theory approach to diffusion in white dwarfs and neutron star envelopes. Astrophys. J., 829, 16(2016).

[65] J. Clérouin, N. Desbiens, P. Arnault. Parametrization of pair correlation function and static structure factor of the one component plasma across coupling regimes. Phys. Plasmas, 23, 092120(2016).

[66] M. M. Gombert, C. Deutsch, H. Minoo. Classical modelization of symmetry effects in the dense high-temperature electron gas. Phys. Lett. A, 66, 381(1978).

[67] J. P. Hansen, I. R. McDonald. Microscopic simulation of a hydrogen plasmas. Phys. Rev. Lett., 41, 1379(1978).

[68] G. Röpke, A. Esser. Debye-Onsager relaxation effect in fully ionized plasmas. Phys. Rev. E, 58, 2446(1998).

[69] T. G. White, P. Mabey, S. Richardson et al. A strong diffusive ion model in dense ionized matter predicted by Langevin dynamics. Nat. Commun., 8, 14125(2017).

[70] A. V. Plyukhin. Generalized Fokker-Planck equation, Brownian motion, and ergodicity. Phys. Rev. E, 77, 061136(2008).

[71] Y. S. Kim, R. G. Gordon. Theory for the forces between closed-shell atoms and molecules. J. Chem. Phys., 56, 3122(1972).

[72] R. G. Gordon, Y. S. Kim. Theory of binding of ionic crystals: Application to alkali-halide and alkaline-earth-dihalide crystals. Phys. Rev. B, 9, 3548(1974).