Shock compressibility of iron calculated in the framework of quantum-statistical models with different ionic parts

Quantum-statistical calculations of shock compressibility of iron are performed. Electronic part of thermodynamic functions is calculated in the framework of three quantum-statistical approaches: the Thomas-Fermi, the Thomas-Fermi with quantum and exchange corrections and the Hartree-Fock-Slater models. The influence of ionic part of thermodynamic functions is taken into account separately with using three models: the ideal gas, the one-component plasma and the charged hard spheres models. The results of calculations are presented in the pressure range from 1 to 107 GPa for samples with initially densities 7.85, 4.31 and 2.27 g/cm3. Calculated Hugoniots are compared with available experimental data.


Introduction
Iron is a widespread component in various types of structural materials, which are exposed to intensive thermal and mechanical loads [1]. The important technical issue is the prediction of behavior of these materials under extreme conditions, in particular, which are realized at the strong shock waves [2]. At terapaschal pressures, there is lack of reliable experimental data on the equation of state for iron, so different theoretical models are used to estimate thermodynamic properties of the metal [3][4][5].
A way to make realistic calculations of equations of state for various substances over a wide range of temperatures and densities is the application of a quantum-statistical approach [6]. It is based on the Born-Oppenheimer approximation. In the framework of this approximation with respect to the electrons, the atomic nuclei are considered as fixed force centers. When nuclei are displaced, electrons rearrange themselves rapidly, so one can consider the equilibrium state of the electron system for a fixed position of nuclei. Also it is assumed that the equilibrium state of electron system (with accounting for electron-electron and electron-ion interactions) corresponds to the most probable distribution of electrons over their energies. The potential for this distribution is related to electron density by means of Poisson equation.
Additionally, the cell-based approach is used. It provides for the consideration only a single nuclear cell instead of the system with infinitely large number of atoms. For this purpose, one finds the average atomic potential in a certain domain around an arbitrary nucleus, whose position is taken as the origin of coordinates. The averaging of atomic potential are carried out

Equations of electronic part models
The only spherical cell (with radius r 0 ) is considered in the framework of used models. The volume of this cell is assumed equal to the volume, which is attributable to an average of one atom in a substance (Wigner-Seitz cell). The radius of the average atom cell r 0 (in the atomic system of units, in which e = m = = 1) is expressed in terms of the density of matter ρ and the atomic mass A (in atomic mass units m u ) by the formula Here a 0 is the Bohr radius, N A is the Avagadro number. Spherically symmetric averaged atomic potential V (r) satisfies the Poisson equation in spherical coordinates: together with the boundary conditions rV (r)| r=0 = Z, V (r 0 ) = 0, where Z is the nuclear charge, and ρ e (r) is the average electron density corresponding to the averaged atomic potential V (r). In addition, there is the charge neutrality condition for the atomic cell:

The Thomas-Fermi model
It is possible to describe thermodynamics of matter at high densities or high temperatures successfully by using the TF model [9]. In this model, all electrons belong to the continuous spectrum. The equation for TF electron density in the cell is as follows: Here θ is the electron temperature (in units of the Hartree energy E H ), µ is the chemical potential. The value y 0 is responsible for the boundary of the continuous spectrum. For the TF model, this value is zero, y 0 = 0. For the pressure P e and the internal energy E e , following expressions are used [9]: Here I k (x) is the Fermi-Dirac integral: where E eTF0 is the energy of the unexcited electrons of the isolated atom. Its value is determined by condition E eTF = 0 at θ = 0, ρ = 0.

The Thomas-Fermi model with quantum and exchange corrections
Taking into account exchange and quantum second-order corrections in , Kirzhnits implemented the transition from TF to TFC [10]. The expression for the electron density in the TFC model is taken as follows: Here Φ(r) is determined in spherical coordinate system: Atomic TFC potential is related to the TF potential though the correction ζ(r): Value of ζ(r) is obtained from the Poisson equation: To derive the thermodynamic functions of the TFC model, it is necessary to calculate the appropriate functions of the TF model: Values of corrections ∆P eTFC and ∆E eTFC should be small in comparison with values P eTF (5) and E eTF (6). In [10,18], expressions for corrections to the thermodynamic functions were derived, in particular, for pressure and energy: where next notation is used: The constant ∆E eTFC0 is determined so that the value of correction for energy equals zero for infinitely rarefied matter at zero temperature.

The Hartree-Fock-Slater model
Although the electron density of the TF and TFC models is consistent with the atomic potential V (r), the wave functions calculated for this potential yield the electron density, which does not coincide with the original [6]. It is due to the fact that the TF and TFC models do not reflect oscillations of the wave functions. Within the atomic cell, the electron states are divided into three groups: (i) states with a continuous energy spectrum (the corresponding electron density is ρ cont ), (ii) states with a discrete energy spectrum (ρ bound ) and (iii) states of the intermediate group (ρ band ). Accordingly, the total electron density is determined as a sum of three terms: In the HFS model, oscillations of wave functions are effectively taken into account [11]. It is provided by choosing an effective energy boundary of the continuous spectrum ǫ 0 , so that the electrons, which give the main contribution to the oscillating part of the electron density, belong to the discrete spectrum.
For atomic potential V (r), in addition to the Coulomb contribution V c (r), an exchange correction is taking into account as additional term: The Coulomb contribution to the self-consistent potential of the HFS model is calculated similarly to other cell models: The exchange contribution is calculated in the semiclassical approximation. In this paper, for V ex (r), the following interpolation formula [11] is used: The bound electron density is calculated directly in terms of the wave functions: Here ǫ nl and R nl (r) are the eigenvalues and the wave functions corresponding to the principal quantum number n and the azimuthal quantum number l. In the framework of approximation of the average atom, instead of a set of ions in different states, one ion with an average occupation numbers of states N nl is considered. The numbers are obtained according to the Fermi-Dirac distribution: Eigenvalues ǫ nl and wave functions R nl (r) are solutions of the Schrödinger equation: with boundary conditions R nl (0) = 0, R nl (r 0 ) = 0 and the normalization condition To solve the Schrödinger equation (22), the phase method [19] is used. Additionally, for eigenvalues of the wave functions of discrete spectrum electrons, the following condition is applied: where r 1 and r 2 are the only two roots (turning points) of equation p(r) = 0, r 1 < r 2 ≪ r 0 . If this condition is not valid, or there are more than two roots, then this state can be assigned to the intermediate group.
Due to the influence of neighbor atoms the energy spectrum of electrons of intermediate group consists of a number of bands of allowed energies [6]. The band structure is taking into account via periodic boundary conditions [20]. The formulas for electron density and density of states take on the form where k = r 0 | k|, k is the quasi-momentum. Summation is carried out over n (the principal quantum number), l (the azimuthal quantum number), m (the magnetic quantum number). Within the cell, the wave functions are expanded in spherical harmonics l ′ with the corresponding coefficients A nlml ′ . The wave function for electrons of the intermediate group is normalized as usual: In the HFS model for the determination of effective boundary between continuous and discrete spectrum ǫ 0 , the following condition is used: Here the summation is carried out only over the states with energy ǫ nl < ǫ 0 , ǫ nlm (k) < ǫ 0 . The fulfilment of this condition provides the thermodynamic consistency for HFS equations [6]. To calculate the density of electrons of the continuous spectrum ρ cont , equation (4) is used with the value y 0 from (28). Use of the found values of the self-consistent potentials of the HFS model provides calculation of thermodynamic functions, in particular, pressure and internal energy [6]: The parameter E eHFS0 for the internal energy is obtained from the normalization condition E eHFS = 0 at T = 0, ρ = 0.

Technical details of calculations
For all considered models atomic potentials V (r), electron densities ρ e (r) and the radial wave functions R nl (r), R ǫl ′ (r) are calculated on a grid with a constant step in a new variable x = √ r, ∆x = √ r 0 /N . Here N is the number of steps in the grid along the radius of the cell; N = 2000.
The integration of differential equations and computing definite integrals is carried out by the Runge-Kutta method with 4-th order of accuracy with additional using specific Runge-Romberg rule [21]. For the HFS model, states with the value of the principal quantum number n ≤ 5 are classified as the states of the discrete spectrum or the intermediate group. The rest of the states belong to the continuous spectrum. Integration of the quasi-momentum k is carried out on uniform grid with the number of steps N k = 500. The maximum value of l ′ is limited by the value l ′ max = 10. The atomic potential was calculated by the iteration technique with the precision e error = 10 −7 . The convergence of the iterations is determined by the condition where p is the consequent number of iteration step.

Equations of ionic part models
Obtained values for electron pressure P e and energy E e (in atomic system units) are used to determine the total pressure P and the specific internal energy E (in SI units) by expressions: where P i = P i (ρ, T ) and E i = E i (ρ, T ) are the pressure and the energy of ionic subsystem, which should be determined separately. The IG model is the simplest approximation, which takes into account the thermal motion of the ions [12]. This approximation is applicable for heated rarefied matter. The pressure (in units of E H /a 3 0 ) and energy (in units of E H ) of the Boltzmann ideal gas is expressed by Unlike the IG case, the OCP model [13] takes into account the interaction of ions. In the framework of this model, ions are considered as point particles with charge Z 0 = (4/3)πr 3 0 ρ e (r 0 ) moving in homogeneous neutralizing surroundings (of compensating charge −Z 0 ). In the OCP approximation, where ∆E i is the correction to the internal energy [13]: The CHS model takes into account the strong repulsion between ions at distance r * , which is equal to the effective radius of the ion core. The value of r * is calculated by formula 4π r * 0 ρ e (r)r 2 dr = Z − Z 0 .
The expression for the pressure is obtained by adding the correction of the CHS model to pressure of the OCP model, and energy of the CHS model is the same as energy of the OCP model [6], namely: where η = (r * /r 0 ) 3 is the packing parameter.

Calculation of Hugoniots
Under compression of a substance in the shock-wave front, the Hugoniot relation holds expressing the law of conservation of energy [22]: where P , E and ρ are the pressure, specific internal energy and density of the shock-compressed matter; P 0 , E 0 and ρ 00 are the relevant parameters in the initial state before the front. At  extremely high pressures and temperatures, the compression ratio σ = ρ/ρ 0 on the Hugoniot ceases to depend on the pressure and reaches the asymptotic value Here γ = C P /C V is the adiabatic index, C P is the specific heat at constant pressure, C V is the specific heat at constant volume. At high temperatures, both electron and ion gases are almost homogeneous and perfect, and because relativistic effects are not taken into account, then γ = 5/3 [22]. In the present paper, the calculations are performed for pressure less than 10 7 GPa, where the compression ratio still depends on pressure.
In the domain close to normal conditions T 0 = 293 K and P 0 = 0.1 MPa, the used quantumstatistical models are not applicable. In particular, for all considered combinations of models of electronic part (E = TF, TFC, HFS) and models of the ionic part (I = IG, OCP, CHS) the calculated value of the normal density of iron ρ 0EI , determined by expression P 0 = P (ρ 0EI , T 0 ), differs from the experimental value ρ 0 = 7.85 g/cm 3 (table 1).
The experimental data on shock compression of iron are shown in figure 1 on coordinate space pressure versus compression ratio for samples with initial densities ρ 00 = 7.85, 4.31 and 2.27 g/cm 3 . In the range of pressures P < 500 GPa, the experimental position of the principle shock adiabat is determined quite reliably. Instead, for the pressure P > 500 GPa, experimental data have relatively large error in determination of the density of shocked samples. This region of parameters corresponds to data mainly from experiments with nuclear explosions.
Using the Hugoniot equation (37), the density can be expressed in terms of mass velocity U and the shock velocity D [22] as By use of expressions (39), it is possible to estimate the error of determination of density and pressure from the errors of determination of kinematic parameters. In figure 1 for data in the range of high pressures, the error of determination of density, which is obtained by using expressions (39), is shown. Note that for the chosen logarithmic scale, the error of determination of pressure is negligible. The choice of the values of the initial specific internal energy E 0 is ambiguous, because the initial state of samples is outside from the range of application of the quantum-statistical approach. In this paper, a value of E 0 , which is substituted in the Hugoniot equation (37), is determined from the condition where T 1EI is determined from the condition P 0 = P (ρ 0 , T 1EI ). The state of matter at T = 0, ρ = 0 is selected as state, where internal energy is equal to zero. Equation (40) relates the initial value of energy with state at normal pressure and density. But this state could be at temperature differing from standard room conditions. Inasmach as the temperature of matter is not explicitly present in the Hugoniot equation (37), this estimate of initial internal energy is convenient for calculations of the shock adiabats (if it can be implemented). So in this work, the equation (40) was used for calculations by HFS model. Estimated values of E 0 for HFS model are −8.00 (with IG model), −6.24 (with OCP model) and −8.63 kJ/g (with CHS model). For TF and TFC models at any non-negative temperatures T 1EI the condition P 0 = P (ρ 0 , T 1EI ) is not satisfied. In particular, for TF and TFC models the point P 0 , ρ 0 is below the cold curve of iron. In such cases, the value of E 0 was taken equal to −7.44 kJ/g [36]. For this value of E 0 , the Hugoniot equation (37) has a solution only at pressures essentially greater than P 0 . In figure 2, for various values of the parameter E 0 , calculated shock adiabats are shown for the HFS and CHS combination. The main difference in the behavior of the Hugoniots causing by variation of E 0 belongs to the region of small pressures, where the quantum-statistical approach is not applicable. Determination of E 0 by the equation (40) gives a qualitatively correct behavior of the shock adiabats in the region of small pressures, in particular, it provides that the shock adiabat passes through the point corresponding initial parameters.

Results of calculations
Hugoniots P (σ) for various combinations of electronic and ionic parts are shown in figures 3-5. Also the isotherm T = 0 is shown. Cold curve was calculated only for electronic subsystem, because considered models of ionic part give zero contribution to the thermodynamics at zero temperature.  The TF model does not pretend to describe the thermodynamic properties of matter at low pressures. In particular, in the TF model, condensed state of matter does not appear at P = 0. Experimental data for the principle shock Hugoniot (ρ 00 = ρ 0 = 7.85 g/cm 3 ) at pressures P < 3 TPa are situated below the cold curve (see figure 3), which results from the TF model as a function of density. So it is possible to evaluate the agreement between the calculation (by the TF model) with experimental data only at P > 3 TPa. Generally, it is possible to assume that the TF model in combination with all considered models of ionic part does not contradict with the experimental data at this pressure range, nevertheless the IG and OCP models provide better agreement than the CHS model. Similarly, for porous samples of iron with initial density 2.27 g/cm 3 , shock Hugoniots calculated by the IG and OCP models pass through (with taking into account the error of its measurement) the experimental points [25] and [35]. For samples with initial density 4.31 g/cm 3 , there are no the data lying above the TF isotherm T = 0.  During the analysis of the Hugoniots within the TFC model, it is important to determine the area boundaries, where the values of corrections to pressure and energy are small in comparison with the values of pressure and energy of the TF model. In table 2, for different models of ionic part and for different initial densities of samples, the values of pressure, where condition ∆P eTFC /P eTF = 0.25 is satisfied for Hugoniots are listed. For all these values of pressure, the condition of smallness of the corrections to the energy ∆E eTFC /E eTF < 0.25 is certainly satisfied. All available experimental data refer to the region of lower pressures, where the condition of the smallness is not valid.
To take into account the difference, which is caused by selection of various models of ionic part, it is possible to obtain a good agreement of principle Hugoniot calculated by the HFS model with the experimental data at pressures more than 100 GPa. Generally, the CHS model of ionic part in combination with the HFS model of electronic part provides for the best agreement of the results of calculations with available experimental data.
The qualitative difference in the behavior of the TFC and HFS models in the region of strong loading is due to the fact that in the framework of the TFC model shell effects are not taken into account [37,38]. Oscillatory behavior of the shock Hugoniots calculated by the HFS model is connected with the ionization of electron shells with increasing temperature.
At pressures P > 10 TPa, it is possible to observe the qualitative difference in the behavior of the Hugoniots calculated by TF and TFC models on one side and HFS on the other (figure 6). In this area of pressure, the HFS curve has several inflection points and crosses the TF and TFC curves many times before reaching the asymptotic value (the so-called oscillating behavior).  Figure 4. The same as in figure 3 but for the TFC model. This is because the more complex HFS model takes into account shell effects, which manifest themselves in the oscillating behavior connected with the stepwise ionization of electron shells with increasing temperature. In the area of influence of these effects on the position of the shock adiabats of iron, there is no reliable experimental data that would confirm these oscillations with sufficient accuracy.

Conclusion
From among the three considered models of electronic part, the HFS model has the widest range of agreement of calculation results with experimental data on shock compression of iron. In particular, for samples with the normal density ρ 00 = 7.85 g/cm 3 , the results of calculation have a good agreement with experimental data in the area from 0.1 to 3 TPa, where the TF and TFC models are not applicable. At pressures greater than 3 TPa, considering the scatter of experimental points and the difference, which gives various ways to take into account ionic contribution, it is possible to assume that all considered models do not contradict with experimental data. But the best agreement is observed for HFS model in combination with CHS model.