Generalized chemical model for plasmas with application to the ionization potential depression

On the basis of the generalized Poisson–Boltzmann equation derived from the Bogolyubov chain of equations for the equilibrium distribution functions in the pair correlation approximation, a general expression is proposed for the Helmholtz free energy of a system that contains any number of components and whose particles interact via arbitrary potentials. This opens up an extraordinary opportunity to simultaneously treat a whole range of physical effects including partial ionization, quantum effects of diffraction and electron degeneracy, short- and long-range interactions of charged particles with neutrals, finite size effects, etc. It is shown that all medium constituents are tied together in a single screening matrix, whose determinant and trace determine the excess contribution to the free energy. The approach developed is then applied to the problem of the ionization potential depression (IPD) leading to quite simple analytical expressions, which turn out to be useful for various practical purposes. In particular, for a single ionization from the neutral state the IPD is shown to significantly depend on the ionization degree such that it consists of the difference of charged and neutral contributions for a fully ionized plasma and turns non-zero for an almost neutral medium. On the other hand, for a multiple ionization process finite size effects of atoms and ions are demonstrated to be of great importance and accounted for in order to achieve good agreement with experimental data on the IPD under warm dense matter conditions.

Nowadays, it is well understood that one of the frontiers of modern plasma theory is to consistently describe warm dense matter (WDM) states [1] appearing in inertial confinement fusion researches [2,3], shock compression experiments [4,5], pulsed-power optical and free electron laser irradiation facilities [6,7], planetary interiors and dense cold stars such as brown and white dwarfs [8][9][10]. As is common for any branch of the contemporary physics, there exists a whole range of ab initio modeling techniques, such as the path integral quantum Monte-Carlo method [11,12], the density functional theory [13,14] and the Thomas-Fermi molecular dynamics method [15], that provide quite an adequate insight into the processes taking place in WDM, thereby urging a further theoretical systematization of the experimental and simulation data already available.
The physical properties of thus omnipresent WDM are determined by numerous phenomena, which include the interparticle correlations, partial ionization and degeneracy, quantum effects, substantial occupation numbers of excited states of atoms and ions, etc. Exactly for this reason in general and the absence of small parameters in particular, only a few theoretical approaches have been worked out to accurately predict the static and dynamic properties of WDM under various external conditions [11,16]. As for the thermodynamic characteristics, the most justified approach is brought by the physical picture that fundamentally treats the mixture of electrons and nuclei within the grand-canonical ensemble of quantum statistics. However, its internal mathematical complexity usually results in the activity/fugacity expansions, like in OPAL and ACTEX, that still imply additional approximations, such as diagrammatic techniques and alike, to assure the applicability to a broader temperature-density domain and/or distinct materials. Note that within the physical picture the occupation numbers for compound particles, such as atoms and molecules, remain unknown but are as yet needed to recover, e.g. the radiative opacities measured experimentally [17,18]. In contrast to the physical picture, chemical models assume that the bounded species are all equal members of the system and the central role is played by the Helmholtz free energy, which severely relies on the linear mixing for plasmas, atomic and molecular liquids rather than strict analysis of interparticle correlations [19,20]. A subsequent free energy minimization procedure provides the corresponding occupation numbers together with the thermodynamic quantities, which are then widely used for astrophysical applications. It is worthwhile mentioning that thus engaged chemical models are capable of assuring quite a reasonable agreement with experimental results at a much lower computational cost.
Below, we deal with a medium of arbitrary composition with unspecified interparticle interaction potentials and the whole consideration is aimed at developing a generalized chemical model that goes far beyond the linear mixing rule by accurate derivation of the Helmholtz free energy. It is firmly based on the precise handling of microscopic correlations between all matter components and, eventually, the key position is occupied by the macroscopic potentials taking into account collective events in the medium. It is anticipated that this will pave a new way for the investigation of not only the ionization balance and the thermodynamic quantities, but for the study of the transport coefficients as well.
The system of interest in the following is confined in a volume V and contains an arbitrary number s of particle species that are then denoted by the subscripts a, b, c, . . ., thereby taking the consecutive natural numbers from 1 to s. All medium components are assumed to be in thermal equilibrium specified by some temperature T and each component, say a, is characterized by the particle number density n a , its individual mass m a and electric charge e a , if applicable. Notice that even different excited states of the same sort of particles can be regarded as distinct entities.
Because of high density/pressure and rather low temperature, a distinctive feature of WDM is rather strong interactions between all plasma components and for the purpose of systematic treatment of interparticle correlations, we employ the following generalized Poisson-Boltzmann equation [21,22]: which distinguishes between the true microscopic interaction potential φ ab (r a i , r b j ) and its macroscopic counterpart Φ ab (r a i , r b j ), accounting for the collective events in the medium. Here r a i denotes the radius vector of the ith particle with ∆ i being the corresponding Laplace operator; β = (k B T) −1 stands for the inverse temperature in energy units with k B beeing the Boltzmann constant. Note that the Einstein notation is implied right above and everywhere below stating that the summation has to be performed over the repeated indexes of particle species.
It should be instantly stressed that the generalized Poisson-Boltzmann equation is originally rooted in the renormalization procedure for interparticle interactions in the medium, which is completely valid in the framework of the pair correlation approximation of the Bogolyubov chain of equations for equilibrium distribution functions. It is important that the above defined macroscopic potentials strictly coincide both with the potential of mean force and the average pair potential [23].
By means of the Fourier transformation the set of governing equation (1) immediately takes the linear algebraic form and its solution for the Fourier transform of the macroscopic potentialΦ ab (k) can be suitably written via the Fourier transform of the microscopic potentialφ ab (k) in the following tensor form: where ε ab (k) = δ ab + βn aφab (k) is then called the screening tensor and δ ab denotes the Kronecker delta. Being defined in the space of particle species, the screening tensor ε ab (k) turns fairly crucial for the present consideration and is neatly written in the matrix from as follows: It is rather timely to make a few remarks concerning the structure of the master relations (2) and (3). First of all, the whole method described herein substantially relies on the existence of the Fourier transforms of the microscopic potentials, which is not the case, for example, for hard core or Lennard-Jones types of potentials. Secondly, provided that the microscopic potential is a symmetric tensor with respect to the particle species permutation due to Newton's third law, φ ab (r a i , r b j ) = φ ba (r b j , r a i ), the same holds for the macroscopic potential, And, most importantly for the following, the Fourier transform of the macroscopic potential is a pure rational function of β, whose order in β is exactly equal to the number s of the particle species in the system.
To proceed with the generalized chemical model, the Helmholtz free energy of a system of interacting particles, F tot , is conventionally represented as a sum of the ideal, F id , and excess, F ex , parts, viz. F tot = F id + F ex . The ideal part of the Helmholtz free energy rigorously obeys the linear mixing rule so that it consists of independent contributions from all system components, i.e.
where λ a = (2π ℏ 2 /m a k B T) 1/2 refers to the thermal de Broglie wavelength with ℏ being the reduced Planck constant and σ a denotes either the internal partition function or the corresponding spin factor as it is required by the statistical physics. In particular, both free electrons and protons, being elementary particles, have the spin factor of 2, whereas for the composite particles a summation has to be performed over the energy spectrum of bounded states. Note that for the electron subsystem, the partial degeneracy effects can also be taken into account via the simple ideal Fermi gas relations [24].
On the other hand, the excess part, F ex , of the Helmholtz free energy compromises the linear mixing rule and is finally yielded by the self-consistent chemical model as follows: [25,26] βF ex V = β 2 n a n bφab (0) − 1 16π 3ˆd kˆdβ βn a n bφ ab (k)Φ ab (k).
The excess part of the Helmholtz free energy formally contains an integral over the system inverse temperature β, which literally requires for the direct knowledge of the correlation functions behavior across the coupling regimes. The same formally holds for expression (6), however, the corresponding obstacle is hidden in the macroscopic potentialΦ ab (k), which, as is mentioned above, is a rational function of β. Of course, for an arbitrary number of particle species the integration in (6) is significantly cumbersome to perform analytically but, nevertheless, the final expression for the excess part of the Helmholtz free energy can be demonstrated to have the following form whose validity can be simply proved using the relations (2) and (3) by differentiating the second term over β. Formula (7) shows beyond doubt that the free energy is readily expressed via the determinant and the trace of the screening tensor (4), which explicitly depends on the number densities and the microscopic interaction potentials of all particle species. Thus, it is inferred that the excess part of the free energy cannot be decomposed into independent contributions from the medium components that are all tied together by the same screening tensor ε ab (k). Nevertheless, the linear mixing rule can still be recovered if either the original physical system can be divided into non-interacting subsystems or if the off-diagonal contributions to the tensor determinant can be estimated to be numerically irrelevant.
It is also worth emphasizing that formula (7) for the free energy is extremely advantageous from the viewpoint of practical numerical calculations because the integration over the inverse temperature has been fully performed. Despite the fact that both expressions (6) and (7) have been formally obtained under assumption that the microscopic potentials are all temperature independent, it is as yet straightforward to demonstrate that formula (7) holds for the most general case of temperature dependent microscopic potentials, which turns rather important for quasi-classical treatment of quantum-mechanical effects of diffraction and symmetry.
As a matter of fact, the general outlook to be advocated is that the excess part of the free energy (7) can be applied to a vast variety of systems with a rather fruitful idea in mind of including all essential physics into the microscopic potentials. To demonstrate how the current method works in practice to achieve quite visual analytical results, we choose the long-standing problem of the ionization potential depression (IPD), which is the phenomenon that less amount of energy is needed for bounded electrons in atoms or ions to be liberated into the plasma continuum rather than into the vacuum. In the past decade significant progress has been made in the experimental determination of the IPD with the further focus on comparison with two available analytical expressions [27,28]. The first of the analytical approaches was put forward by Ecker and Kröll (EK) [29] in an attempt to generalize the Saha equation with regard to the plasma chemical potential, whereas the second one was developed by Stewart and Pyatt (SP) [30] to virtually interpolate between the Debye-Hückel (DH) and the ion-shpere limits. In spite of the increased experimental resolution in measuring the IPD, neither one of the above mentioned analytical expressions has been unambiguously corroborated. Indeed, on the one hand, data on the spectral line disappearance seem to support the SP model [31], but, on the other hand, experiments with the K-shell ionization in aluminium are in a position to validate the EK formula [32], whereas the results of x-ray scattering measurements on imploding CH spheres [33] do not look to favor any of the two analytical expressions. Such an unsatisfactory situation resulted in a whole flow of simulations and numerical calculations to account for various physical effects, which include the Hartree-Fock-Slater model [34], the density functional theory [35,36], the grand canonical Monte Carlo simulation [37], the quantum statistical approximation for ionic dynamical correlations [38], the atomic-solid-plasma model for single and multiple core hole states [39], the in-medium Schrödinger equation for incorporating the electron Pauli blocking and degeneracy [40], etc.
The above brief overview instantaneously reveals that the discrepancy between available experimental data and theoretical predictions on the IPD is primarily prescribed to that both the EK and SP relations are solely designed to classically treat the charged particle correlations in the plasma medium, completely ignoring other important contributions. Below we propose more general analytical formulas for the IPD that are, to some extent, capable of incorporating a whole range of physical phenomena. They are to be firmly based on expression (7) for the Helmholtz free energy, whose subsequent minimization is to provide the exact medium composition under appropriate external constraints.
To begin with, we take a widely spread physical system at thermal equilibrium that consists of three particle species, i.e. free electrons, singly ionized ions and neutral atoms at the ground state denoted by the subscripts 1, 2 and 3, respectively. By fixing the total number density n = n 2 + n 3 of heavy particles, defining the ionization degree via α = n 2 /n and assuming the system quasineutrality, the Helmholtz free energy minimization yields the following generalized Saha equation where the IPD I(α) is found as with σ 1 , σ 2 , σ 3 being the appropriate internal partition functions and λ 1 , λ 2 , λ 3 denoting the corresponding thermal de Broglie wavelengths. The IPD (9) together with the excess part of the free energy (7) renders the calculation scheme quite effective to take into account various physical effects. The first instantaneous inference to be made is that the IPD (9) strongly depends on the ionization degree α and, at the same time, naturally embraces the neutral component, which is neglected both in the EK and SP approaches. In order to demonstrate the importance of neutrals, let us first assume that the plasma is almost fully ionized with the ionization degree α ≈ 1. Then the IPD I high for this regime of high ionization can be rewritten as the difference between the charged and neutral contributions: which are defined via simple integrals with the following notations: The green line represents the neutral contribution In; the red one corresponds to the charged contribution I ch ; the blue line provides the total IPD of I high , while the black dashed line and the black dotted lines display the results obtained within the SP model [30] and the EK model [29], respectively.
On the one hand analytical formulas (10)- (14) are both visual and helpful for evaluating the IPD in various types of plasmas, which can be reasonably realized by a suitable choice of the microscopic potentials. At the same time it is straightforward to demonstrate that for the pure Coulomb interaction potentials and complete elimination of the neutral component, the EK relation for the IPD can be easily recovered, which in this case coincides with the results of the DH theory.
From the theoretical point of view final formula (10) is rather expected beforehand since the IPD is actually determined by the difference in the excess free energy before and after ionization process. From the numerical perspective, however, it has still to be demonstrated that the neutral contribution really matters, which is systematically performed below for warm dense hydrogen that is characterized by the Coulomb coupling Γ = e 2 /ak B T and the density r s = a/a B parameters, where e is the elementary charge, a = (3/4π n) 1/3 stands for the mean interparticle spacing, a B = ℏ 2 /m e e 2 and m e denote the Bohr radius and the electron mass. To treat the quantum diffraction effects in mutual interactions between free electrons and protons, we choose the microscopic potentials in the form of the improved Kelbg potentials [41] with the pure electronic part supplemented by the density-dependent symmetry effects as stated in [42]. As for the interactions between the charged components and neutral hydrogen atoms, the corresponding microscopic potentials are selected in the adiabatic dipole approximation to take into account both the short-range ion-core and long-range polarization effects as suggested in [43]. In figure 1 the IPD I high in a hydrogen plasma is plotted in eV as a function of the coupling parameter in the range of Γ = 0.01 ÷ 0.40 at the fixed value of the density parameter r s = 5 and in the same figure both the EK and SP analytical relations are displayed as well. It can be shown that the quantum diffraction normally suppresses the IPD, but the symmetry (degeneracy) effects are responsible for the opposite action to finally dominate in the charged component contribution I ch depicted in figure 1 in red. It is also seen that the neutral component contribution I n , painted in green, cannot be utterly ignored since it reaches almost 10% of the ionization potential (≈13.6 eV) of an isolated hydrogen atom. Additionally, it can be analytically demonstrated that the contribution of neutrals is proportional to a B /r D as compared to the EK model with r D being the Debye radius. Detailed numerical estimations show that the relative contribution of neutrals reaches approximately 0.5 percent when a B /r D approaches 0.22 or, consequently, it really matters in the range of Γ > r 2 s /125. Note that the total IPD I high reveals a more sophisticated behavior on the coupling parameter, exceeding both the EK and SP predictions at very weak couplings and lying in between them at larger values of Γ.
To contrast the above considered case of high ionization with α ≈ 1, we now turn to the situation, in which the hydrogen medium remains almost neutral. On putting α ≈ 0 into the general expressions (9) and (7), the IPD I low for this low ionization state is ultimately obtained in the following form with the notation We therefore predict in equation (15) the existence of the IPD even for a practically neutral environment, whose behavior is largely driven by the interaction between hydrogen atoms. To carry out numerical estimates, presented graphically in figure 2 for the coupling parameter range Γ = 1.4 ÷ 2.6 at some fixed values of the density parameter r s , the microscopic interaction potential of hydrogen atoms has been chosen in the form of the Morse potential [44], which is quite proper for the description of the pressure dissociation of hydrogen molecules. Figure 2 clearly confirms that the IPD grows when the number density of atoms increases and reaches rather substantial magnitudes in comparison with the ionization potential of an isolated hydrogen atom. Note that the specified region of the plasma thermodynamic parameters has been used to assure that the ionization degree stays low and, at the same time, the formation of hydrogen molecules is effectively tackled. It is worth stressing that in figure 2 direct comparison with the EK and SP predictions seems meaningless since they both formally turn zero because the number densities of charged particles virtually vanish.
Hence, it is inevitably concluded from the two limiting cases considered above that for the warm dense plasma domain the neutral component cannot be neglected in the theoretical analysis when considering the ionization of neutral atomic or molecular states. Note that the corresponding contributions are completely absent in both the EK and SP approximations for the IPD.
With the purpose of comparison with available experimental data, an attention is centered on a somewhat different setting when the WDM consists of free electrons and two sorts of positively charged ions with the respective charge numbers Z and Z + 1, correspondingly denoted in the following by the subscripts 1, 2 and 3. Consequently, the focus is shifted to the multiple ionization process when an ion in the charge state Z + 1 is created by the ion in the charge state Z emitting a free electron. As in the above case of the hydrogen plasma, the total number density n = n 2 + n 3 of heavy ions is presumed fixed, the ionization degree is then defined via α = n 2 /n together with assuming the system local quasi-neutrality, such that the IPD for the case of almost total ionization α ≈ 1 is derived from expression (9) as where On the one hand, in the case of pure Coulomb interactions between the plasma constituents the DH limiting form can be correctly restored from formulas (18)- (20). On the other hand, it can be easily understood that in order to study the IPD, it is necessary to take into account the finite size of ions or of the so-called ion-core effect. Indeed, the ionization energy of an isolated ion is roughly proportional to its charge and inversely proportional to its size, whereas the IPD is determined by some characteristic distance, which could be either the average interparticle distance or the Debye screening radius. For a positive experimental  [30], and the DH and the modified EK [29] models, respectively; the filled circles with error bars provide the experimental data of [32,46]. detection, the IPD must be of the order of the ionization potential itself, which means that the size of the ions must be comparable with the mean interparticle spacing or with the Debye screening length, i.e. the dimensions of ions cannot be ignored when considering the surrounding plasma distribution. It is simple to numerically verify that the above speculations are especially true for the WDM conditions, although, only point-like particles are treated in both the EK and SP theoretical estimates, which is definitely one of the possible sources of discrepancy when a straightforward comparison is undertaken with experimental data.
In order to embody the finite size of ions, a soft empty-core potential is employed in the following form [45]: φ ab (r) = e a e b /r [ , where e a , e b are the charges of interacting particles separated by some distance r and d ab refers to the cut-off distance. In particular, we assume that d 11 = 0 for the interaction between electrons, however, the quantum effects are still treated in the way it has been done above for the hydrogen plasma to account for the diffraction effects and electron degeneracy. As for the interaction of electrons with both sorts of ions, it is assumed, for the sake of simplicity, that in the interionic potentials d 12 = d 13 = d, whereas d 22 = d 33 = d 23 = 2d. Figure 3 depicts the numerical results for the IPD in the aluminum plasma with the mass density ρ = 2.7 g cm −3 and temperature T = 200 eV. The finite size effects are displayed for three particular values of the dimensionless size parameterd = δ/a normalized to the mean interionic spacing a = (3/4π n) 1/3 . It is observed that the IPD drops when the dimensionless size parameter grows and gradually covers the whole domain in between the DH and SP analytical expressions, also reproduced in figure 3 together with the EK results. It is also seen that the experimental data of [32,46] are well fitted for the size parameterd = 0.08, which corresponds to d = 1.26 · 10 −11 m. In reality the ionic size must vary with its charge number Z as well and, moreover, the IPD itself turns out to be rather sensitive to the choice of the form of the interionic potential. Thus, figure 3 clearly demonstrates the importance of handling the finite size of ions to smoothly describe the experimental data. Notice that, as evidenced by the present numerical estimations, the quantum effects of diffraction and symmetry play an important role for rather high Z values, reaching almost 10 eV, which is twice as large as the experimental error.
Finally, it is possible to conclude that the generalized chemical model developed herein seems to be suitable for the description of various physical systems since it is capable of simultaneous accounting for various effects such as the partial ionization, the electron degeneracy, the quantum diffraction, the shortand long-range interaction of neutrals as well as the formation of molecules. In addition, the present approach can be easily extended to study not only the IPD but the thermodynamic and transport properties, which is to be implemented elsewhere in due course.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).