Hydrogen, helium and lithium plasmas at high pressures

The equations of state (EoS) and other thermodynamic properties of plasmas of the light elements H, He, and Li, are calculated using inverted fugacity expansions. Fugacity expansions are known as an alternative to density expansions but show often an inferior convergence. If, however, the inversion can be solved, the fugacity representations may be very efficient. In particular, the contributions of deeply bound states are included in the fugacity expansion in a very effective way. The mathematical problems on nonlinearity connected with the inversion of fugacities to densities are reduced to solvable algebraic problems. The inversion of fugacities to densities is solved separately for two density regions: (i) In the low density, non-degenerate region we consider ring contributions describing screening effects and ladder contributions describing bound state formation. (ii) In the high density, degenerate region the electrons are described by the Fermi–Dirac distribution. Hartree–Fock contributions and Pauli blocking have to be taken into account. The ions are considered as classical, strongly correlated subsystem eventually forming a Wigner lattice. We solve the inversion problem for each of the regions. Near the crossing point, the separate solutions are connected to each other, either by smooth concatenation at the crossing point or by Padé approximations.


Introduction
This work, dedicated to David Blaschke's 60th birthday, on studies of hydrogen, helium and lithium plasmas is based on several methods developed together with David Blaschke [1][2][3][4] as well as on other recent works [5]. The light elements have in common, that the lowest lying bound states, the hydrogen atom, the hydrogen molecule, and the helium and lithium atoms consist of 2-4 charged particles. In each case, the outer electrons are only loosely bound with ionization energies of I  Li = 122.42 eV is required. For full thermal ionization of all helium or lithium electrons we need far more than 100 eV, i.e. more than 10 6 K.
An alternative way to reach high ionization degrees is to increase the particle density to regions of n > 10 23 cm −3 where all atomic and molecular bound states are destroyed by screening and Pauli blocking effects 1 (pressure ionization) [8][9][10][11][12][13][14][15]. For an estimate we consider that a bound state of a nuclear charge +Ze with one or two electrons in an s-shell needs space corresponding to a sphere with a radius r B (Z) = /(m e e 2 Z) (we use Gaussian units). This way the highest density of bound states of nuclei with charge Ze still consistent with a close packing is Earlier work to describe high density and high pressure effects was often based on the chemical picture [13]. In contrast to the physical picture, which normally uses perturbation expansions, the chemical picture has the advantage that it is based on a variational principle. On the other hand, the chemical picture leads to complex mathematical problems including coupled nonlinear equations with solutions, which are sometimes not unique and sometimes even unstable. The aim of this work is to describe an alternative and in many cases more simple method based on the algebraic inversion of fugacity expansions. This method is effective for bound states which consist just of a few charged particles, otherwise the order of the algebra becomes too high. For the light elements H, He and Li, we have to deal with polynomials up to the fourth order. Our approach may be of interest, e.g., for plasmas playing a role in nuclear fusion and in stellar atmospheres because they are dominated by light elements.
One of the so far most accurate equation of states (EOS) has been the one derived from the activity expansion method (ACTEX) by DeWitt and Rogers [16,17]. It has been proven to be very reliable for the sun and sun-like stars, consisting mainly of hydrogen and helium. Closely related is the OPAL opacity project of Livermore [18], which is based on a numerically inverted fugacity expansion with Debye screening. The method has been very successful for applications to plasma at the surface of stars [18]. However, the OPAL equation of state is limited with respect to two aspects. Firstly, it is only available in the form of pre-computed tables that are provided by the Lawrence Livermore National Laboratory. Therefore, applications require interpolations, which always lead to loss of accuracy. Secondly, the OPAL equation of state is proprietary and not freely available. Our method is comparable to OPAL. However, it is an analytical inversion based on the solution of polynomials and is open to the public. It is related to our fugacity expansion method developed earlier [19][20][21] and a different approach developed more recently by Alastuey et al. [22]. In this paper, we develop an analytic inversion of fugacity representations based on the grand canonical ensemble [8,19,[23][24][25].
The OPAL equation of state has to be improved investigating the ionization potential depression at high densities where electron degeneracy becomes relevant, see references [15,26] and references given there. As pointed out by the authors in collaboration with David Blaschke [1][2][3]15], the physical properties of dense hydrogen and other light elements is a hot topic which requires the development of specific methods [1][2][3]13,23]. The basic problem was raised in the past by many authors, starting with Wigner, Huntington, Abrikosov, and others. Of special interest is here again the region of extremely high pressures in the Mbar region, where the transition of hydrogen from a dielectric state to a highly conducting phase is observed, which is considered to be a type of Mott transition. Among the main effects influencing strongly the pressure in this region are: (i) the strong degeneracy of electrons including Hartree-Fock effects, (ii) the strong Wigner -like interaction of the ions, (iii) the lowering of the bound states by screening effects, (iv) Pauli blocking effects.
In some of our earlier works we attacked the present problem by a combination of the chemical approach with Padé approximations, calling this PACH (Padé approximations in the chemical picture) [13]. Here we study the algebraic inversion problem separately for high and low densities and connect the solutions to each other in the spirit of a method proposed by Zimmermann [27], either by a smooth concatenation at the crossing point or by Padé-like interpolations.

The method of fugacity inversion and Fermi-Dirac laws
Within a quantum statistical approach, we choose the grand canonical ensemble to describe the plasma. The state of the plasma is characterized by the temperature T or β = (k B T ) −1 and the chemical potentials µ c corresponding to the conserved components, in our case electrons (c = e) and nuclei (c = i). Instead of chemical potentials we use here fugacities defined as with the thermal wave length Λ c . From the partition function, we obtain the grand canonical thermodynamic potential J = −pV , which defines the pressure p = p(T, z i , z e ) in the fugacity representation. Thus, the pressure in the grand canonical ensemble serves as basic quantity for the description of the thermodynamics of plasmas [16,19,24,[28][29][30]. From this follows the expression for the densities (d = c) of the species c, d in the fugacity representation, i.e., For practical applications, we like to characterize the state of the plasma via the densities, besides the temperature. Expressing the fugacities by densities and in this way also the pressure by densities is the inversion problem which is our main task here.

3406
The European Physical Journal Special Topics Fugacity expansions for plasmas were first given by Montroll and Ward including only ring diagrams in 1958 [31], then by Vedenov and Larkin in 1959 [32] and by Abrikosov et al. [33] including also ladder diagrams. In parallel, Kadanoff and Baym [34] presented a general formalism including all quantum effects based on representations within the grand canonical ensemble. Note that the relations given above for a two-component plasma can be generalized to a multiple-component plasma. Within a chemical picture, we can also consider possible ionization stages of the atoms as new components of the plasma. However, the corresponding particle numbers are not conserved quantities of the plasma. Ionization and recombination as reaction processes determine the chemical equilibrium which minimizes the free energy of the plasma. We will not continue with this chemical picture which becomes problematic in dense plasmas but start from the physical picture using only conserved quantities.
In principle, the pressure may be expressed as [8,[23][24][25] where the main contributions, which will be defined later, are: Here the first two contributions connected with the names Fermi-Dirac and Hartree-Fock are well known, at least for the limiting cases. They are stemming from diagrams with no more than one interaction line. The most important diagrams with two or more interaction lines are the contributions stemming from ring and ladder diagrams p ring i and p ladder i . At very high densities, the contributions of rings and ladders are less relevant since Coulomb bonds need some space to be developed and this space is not available at very high densities, where typical distances are below the Bohr radius. Bound states are, from a technical point of view, contributions connected with ladder diagrams, At densities larger than the estimated maximum density for bound states (1), these bound states disappear and we arrive at regions where only contributions of free states and Hartree-Fock contributions still exist.
Following the idea of Zimmermann [27], we use different methods to treat the low density, weakly degenerate plasmas (n e Λ 3 e < 1) and the strongly degenerate plasmas (n e Λ 3 e > 1). In order to illustrate this, we discuss the well known inversion for the EOS of Fermi-Dirac gases [27]. The Fermi-Dirac contributions to the EOS are exactly given by the Fermi integrals [27,35,36]. With spin factor 2 we have 1 2 n FD e Λ 3 e = I 1/2 (α e ); α e = µ e /k B T.
The problem of inversion from the variable z e = 2Λ −3 e e αe to the density n e can be considered as solved in this approximation, with the expression given by Zimmermann proposed to concatenate both pieces at the crossing point y 0 5. We will use here a more smooth transition by means of weight functions decaying with increasing distance from the crossing point. For a smooth interpolation between both pieces we may use two half-side functions w(ay) and w(−ay) related to the tanh(y) function, where y is the degeneracy parameter defined above, and a is a free parameter [37]. Another alternative is to connect both half-sides of the functions by appropriate Padé interpolations [37]. In a similar way, we will proceed later for all other contributions. Most attention will be paid to the high-density part which is evidently less studied by analytical methods so far.
3 Fugacity inversion for bound states

Inversion for hydrogen plasmas
After explaining our approach for an ideal Fermi-Dirac gas as a rather simple example, we consider now the more difficult case, which is our main task here, the inversion of relations (4) with taking into account bound state contributions. Within the general systematics, the lowest order corrections with respect to interaction are the Hartree-Fock terms which can be incorporated in the quasiparticle picture. This will modify our approach given in Section 2 and also the coefficients of the expansions (7). Note that the Hartree term for Coulomb interaction is diverging what is cured if screening is taken into account, but for charge-neutral plasmas the Hartree term is zero. However, this is not applicable for the plasmas we are considering because the interaction is strong, forming bound states. Improving the ideal Fermi-Dirac gas at low densities, we have to consider the expansion with respect to the density, a so-called cluster expansion, where we take into account binary collisions, three-particle collisions, etc., but in all orders of interaction. The following consideration is fosussed on the bound state contribution and is based on the preliminary drastic assumption, that we may omit all contributions from ring diagrams and exchange diagrams and keep only those contributions of ladder diagrams, representing bound states.
In order to demonstrate our method for the treatment of bound states developed in [20] we consider atomic hydrogen first. A simplification is obtained assuming z e = z i = z, what corresponds in a way to charge neutrality and is reasonable for regions where degeneracy is not relevant. We obtain pressure and densities as convergent series with respect to the fugacities where we use the traditional notation with the reduced mass m −1 ab = m −1 a + m −1 b . Note that the definition of the thermal wave length for the relative motion differs from the definition (2) for the center-ofmass motion. The second virial coefficient b ie (T ) contains contributions of bound and scattering states. The double dash denotes the subtraction of the first order of interaction which is already considered in the Hartree-Fock term. In standard approximation, the convergent bound state part of the partition function appearing in the second virial coefficient is given by the Planck-Brillouin-Larkin expression [8,24] Using a dimensionless fugacity Z = z e /n e = z i /n i , the quadratic equation (9) reads with the solution [20] More convenient for calculations are Padé approximations as e.g.
valid for Z(x) > 1/2. For the pressure follows the representation In the chemical picture, the pressure of hydrogen my be written as (the stars denoting "free densities") βp = n * e + n * i = n e + α n i .
This writing suggests, that Z[n i b ie (T )] corresponds to the degree of ionization α, which is the ratio of the free electron density to the total electron density [8,20] (not to be confused with the activity α c used above). We now extend our approach and include the formation of hydrogen molecules. Within a virial expansion, the H 2 molecule appears as a contribution of the fourth virial coefficient. Neglecting weakly bound ion states like H + 2 or H − , we have a leading contribution of the 4th virial coefficient from the H 2 molecule: The indices denote the combinations of particles in the cluster and the prefactors are due to the total number of physically equivalent combinations. Due to n e = n i we have again z e = z i . The inversion requires the solution of a 4th order algebraic problem for Z. We find a fourth order polynomial in Z, with the assumption a 3 0, i.e. neglecting the terms stemming from hydrogen ions, and get βp = n e + n i Z(a 2 , a 4 ); The solution depends on density and temperature through the coefficients x = a 2 (n i , T ) and y = a 4 (n i , T ), which provide the contributions of bound states. For y = 0 we get back to the known solution for atomic hydrogen equation (10). For small arguments x < 1, y < 1, the function Z(x, y) may be expressed via a Taylor series Z(x, y) = 1 − x − y + 2x 2 + 4y 2 + 6xy + ... (19) In the case that the asymptotic behaviour for one variable is finite and the other one goes to infinity, we have a simpler asymptotic solution, as e.g.
Note that this is an estimate, found by using the analytical expressions of the limiting cases and comparison with numerical results. For still more precise values, we may use numerical iterations of the polynomials which may start from the given approximation. The chemical potentials and the pressure are given by k B T µ e = ln(n e Z(a 2 , a 4 ));

Inversion for helium and lithium plasmas
Within the schema developed above we study now the formation of helium atoms in plasmas, which consist of a doubly charged nucleus (α -particle) and two electrons. This means that three-particle interactions play an important role and we need at least the third virial coefficient to describe the helium atom formation in the Rutherford picture. In order to have a guideline we start with a chemical picture, where the essential contributions to the pressure of He-plasmas can be given as (the stars denote again free densities) βp = n * e + n * i + n * He + + n * He , n * He + = n * e n * i K He + (T ), n * He = (n * e ) 2 n * i K He (T ). (24) We assume, that the densities are connected via an ideal mass action law with mass action constants K He + (T ) and K He (T ) for formation of helium ions and atoms, respectively. We identify corresponding terms in the fugacity expansion in the following way [38]: electron fugacity, ion fugacity and the third virial coefficient with free electron density, free ion density and mass action constant K He (T ), respectively, According to equation (3), we then find et expressions for the densities n i , n e : Taking into account the physical meaning of the contributions to n i we may introduce the fractions of free species as was suggested in some earlier work [11,39]. Although there is no full correspondence between the standard chemical view and the present physical picture, we may approximate α 2 z i /n i as corresponding to the degree of double ionization in the case of helium. Due to the relations between the charges and densities n e = 2n i , it follows that z e = 2z i for Z = z i /n i = z e /n e and the cubic equation The inversion requires solving a cubic polynomial. Neglecting the less important ions He + , i.e. a 2 = 0, the quadratic order is missing, so the problem belongs to a normal form with the Cardano solution. In our case the discriminant is positive since b eei (T ) > 0 and we get one real solution The fugacity function Z(y) depends here only on one dimensionless variable a 3 . Since the Cardano solution is a rather complex expression, we construct a simpler Padélike expression starting again from equation (27), based on the limits for small and large arguments of a 3 , respectively, Better approximations may be found numerically.
The chemical potentials and the pressure of helium plasmas are then given as We may express the pressure related to the Boltzmann pressure with the correct limits Finally we study lithium plasmas using the same methods and similar approximations as for hydrogen and helium. Lithium atoms are formed by a nucleus carrying three positive charges and three electrons in the atomic shells. The outer electron is only loosely bound with an ionization energy of I = 5.39 eV. The relevant bound state consists of 4 charged particles and this makes some mathematical similarity to the inversion of the cluster series for H 2 -plasmas since it may be reduced to a fourth order polynomial. We are starting again from a chemical picture, denoting the concentrations in this view by stars. Neglecting singly charged Li + -ions, the essential contributions to the pressure have the form Again we assume for now, that the densities are connected by ideal mass action laws with appropriate mass action constants. Going to the grand canonical ensemble in the physical picture, we assume for the pressure Applying equation (3), expressions for the densities are With z e = 3z i , z i = n i Z we get a fourth order polynomial with respect to Z: and for the pressure in terms of the fugacity function Z(a 2 , a 3 , a 4 ): In order to simplify, we may either neglect the singly charged ions Li + , as proposed above, or fuse it with the energetically close Li-atoms in order to get, in full equivalence to hydrogen, the more simple polynomial 1 = Z + a 2 Z + a 4 Z 4 . With this assumption, the solution may be approximated by the Padé formula in analogy to hydrogen. Let us underline again, that so far we considered only the contributions from ladder diagrams with three or more ranks, for methodological reasons. They are of the orders e 4 , e 6 , e 8 . In the next section, we will include the missing ring contributions which are connected to the screening effects.
4 Combination of screening and bound state effects in dense plasmas at moderate densities

Region of hydrogen atom formation
Taking into account screening effects by including contributions of the ring diagrams has been worked out by Kraeft et al. [24]. For the low density region considered here, we introduce the so-called reduced mass approximation (RMA) and neglect exchange effects [21,24,37] for the low density region. The RMA is based on the assumption that only collisions of opposite charges are relevant, which follows from the asymptotic properties of the so called quantum virial functions Q(ξ ab ) for ξ ab = βe a e b /(4π λ ab ) 1, for details see [21,24,37]. In the same region, the exchange contribution decreases exponentially [24,37]. The well -known result for pressure contributions stemming from ring diagrams is neglecting symmetry effects [8,19,24] where κ g is the Debye screening parameter as it appears in the grand canonical ensemble, i.e. expressed by fugacities instead of densities, and λ = λ ie due to RMA

3412
The European Physical Journal Special Topics here and in the following. Including higher density orders we get [21] x + 3 10 For this so-called quantum Debye-Hückel approximation (QDHA), the inversion leads to screened fugacities [8,21,24] Taking into account that κ depends on z too, and replacing z by a screened fugacitỹ z and κ byκ, the following expression is obtained [21,37] We insert this result into the inversion equations (17)-(18) for the bound states and get for the screened fugacity functionZ =z/n This is an approximation which is certainly not applicable in the case of degenerate electrons. But we can combine the inversion procedures for screening and for the bound states using the results of the previous section. The further procedure depends on the specific plasma. Let us start with hydrogen and include only atomic states. Combining the approximate solution for the screening contribution (43) with the result for the bound state contribution for hydrogen (13) we get for the fugacities and the pressure respectively. The latter expression is closely related to the Saha equation in QDHA [8,24,37]. This confirms the idea, expressed already in the previous section, that the Z-function corresponds to the degree of ionzation α of electrons and ions, see equation (16). Taking into account Z(x) 1 − x + ... we see that this result is compatible with both limits of "no screening" and "no bound states". For the pressure we find relatively simple formulae which contain three nonlinear functions having convenient Padé approximations. Note that so far our results refer only to the region of non-degenerate electrons. There are several ways for including the degeneracy effects of the electrons. The easiest one is the replacement of the first contribution by the Fermi-Dirac pressure. This will effect the pressure only at higher densities, where it will increase rapidly, Here the Planck -Brillouin -Larkin partition function appears, instead of the full two-particle function, because all terms linear in e 2 are already taken into account by the screening procedure and should not appear twice. This effect is also responsible for the convergence of the partition function [8,24]. Since the fugacity series contains more terms than the density series, we may expect that some of the difficulties connected with the convergence of the density expansion, as, e.g., the strong decrease of the pressure going possibly even down to negative values at high density, may be avoided. Indeed, a representation of the pressure using the fugacity expansion shows a more reasonable behavior with increasing density. It is much closer to the behavior obtained in the chemical picture using Saha-type mass action laws [13], as Figure 1 demonstrates. Taking account all contributions up to second order in the fugacity, the pressure is shown for three temperatures. No negative pressures appear, which we would observe in a pure density expansion. This is due to the fact that the fugacity is not fixed, it decreases with increasing value of the partition function and avoids in this way negative pressures. The good convergence of the new expression (51) including the nonlinear functions Z(x) and φ(x) based on the fugacity series (43) is due to the fact that the fugacities of the electrons and the protons are rather small in the bound state region. We notice again that this result comes from an extension of the density expansion including all quadratic terms from the fugacity virial expansions and the full Fermi pressure of the electrons. We see that the overall behavior of the new representation is much better than that of mere density or mere fugacity representations. Therefore we may conclude that the most appropriate description of Coulomb systems is by inverted fugacity expansions. The new contributions allow to control the exponential increase of the second virial coefficient in the bound state region. The physics behind is that the fugacity expansions describe well the saturable forces between bound states in Coulomb systems, but on the other hand the additive long range Coulomb forces and the screening effects are better represented by density series. It is consistent with the Saha theory and in particular also with the quantum Saha-Debye-Hückel theory developed above. Further we may state that the new extended analytical theory is, up to the density n 10 21 cm −3 , in reasonable agreement with earlier calculations based an the chemical picture and quite complex numerical codes [13]. We stress the following structural differences of equation (51) to the original Saha theory: 1. The exponent exp(βI) in the Saha equation is replaced by the PBL-partition function. 2. The first Debye-Hückel screening function, which is linear in κ, is replaced by a nonlinear function φ(x) describing the quantum statistical ring sums in the grand canonical expansion. 3. The linear second virial coefficient is replaced by a nonlinear ladder function stemming from the representation in the grand canonical ensemble.
The grand canonical ensemble plays an essential role in our derivation. Equation (5) corresponding to a second iteration is in some sense incomplete, e.g., we would expect that in higher iterations more κ-terms are replaced and that more terms corresponding to a mass-action law will appear.
In the region of low temperature T < I/k B and non-degenerate plasmas our rather simple formulae give a rather good behavior and describe well the transition from low density to the valley of bound states. We note however that several physical effects as, e.g., plasma phase transitions are not yet described by the present approach.
Evidently this effect appears only taking into account higher order terms or after transition to some chemical picture [30,37,46]. The treatment of phase transitions is, however, not our aim here, knowing that the description of phase transitions in the grand canonical ensembles requires special tools. Anyhow the overall agreement of mixed representations with chemical descriptions is quite reasonable, the deviations increase only at large densities beyond n 10 21 cm −3 , see Figure 1. The largest 3414 The European Physical Journal Special Topics deviations appear near the minimum of the relative pressure for hydrogen at 20000 K (the green line is the present theory and the red line, the lowest of the curves comes from the PACH approximation [13]) and may be interpreted as due to the formation of molecules which have been neglected in this section so far.

Region of hydrogen molecule formation
In order to describe molecule formation we should include contributions from the fourth order virial coefficient by replacing the solution of a second order polynomial Z(x) by the solution of a fourth order polynomial Z(x, y) as described in Section 3.1, equation (18). This leads to closed equations of the shape in the low-density region. Here we introduced the "screened virial coefficients"ã 2 = γ andã 4 = ζ defined by The hydrogen partition functions are defined as σ H (T ) = σ PBL (T ), see equation (11), and where the g ν denote degeneration factors of the energy levels E ν (electronic and ionic) of the H 2 molecule . The subtraction of the first linear terms of the exponential function is introduced to avoid double counting, since all terms up to e 4 were already taken into account in the summation of ring diagrams.
and the corresponding ground state contribution is Investigating Figure 1, we see small shoulders at intermediate densities which may reflect the formation of atoms and minima at n i ∼ 10 23 cm −3 due to the formation of molecules. For still larger densities, all bound states break down due to the dominance of the Fermi pressure which steeply increases. However, in this region the energy level shifts, which are not yet included, may start to play a role. This will be discussed in forthcoming sections. For the temperature T = 20000 K (in green) the comparison with earlier calculations (lowest curve in red), based on the chemical picture including molecules and Padé approximations PACH [13] shows that our new results for the relative pressure are a bit deeper. Let us draw some conclusions: We have seen that density as well as fugacity expansions have both advantages as well as disadvantages: 1) The density expansion describes well the screening effects but it fails to cope with rapidly increasing contributions from the screening terms and from the BPL-partition function. 2) The fugacity expansions are more smooth, since any rapidly increasing term makes the fugacity factors less relevant. Any finite fugacity expression corresponds to an infinite density series including the partition function σ. If σ is getting large, then the fugacity goes to zero what guarantees even at large densities always finite contributions to the pressure. This is true for the screening contributions and for the bound state contributions.
3) The inverted fugacity expansions used here are useful smooth representations of the pressure. This procedure of inverting fugacity series provides Saha-like terms and therefore chemical effects.
In Figure 1 we give a comparison of the present theory for T = 20000, 30000, 50000 K with the numerical results obtained within an advanced chemical picture by minimization of the free energy [13]. At the present moment the agreement is still only qualitative. This shows that our new methods still need further development. This is basically due to missing terms in the interference between screening and ladder diagrams. With the chemical picture we mean here the Saha equation as well as any description which is based on the optimization of expressions for the free energy depending on the density of free charges, atoms and molecules [11,13]. Based on these observations we may recommend the use of the present alternative to the relative complicated chemical descriptions. We may conclude that the mixed expansions derived first in [21] combine the positive features of density and fugacity expansions avoiding most of the negative features. We should mention that the present approach fails in the region of large densities. The reason is, that at higher densities the Hartree-Fock effects and the Wigner -type ion -ion interaction effects provide Helium plasmas at T = 30000, 50000, 80000 K (red, green, blue curves) obtained by inversion of fugacity series as a function of log(ni). In the region of stronger screening and atom formation the relative pressure decreases from 1 to 1/3 corresponding to the binding of three particles in one helium atom. The thin lines in magenta, turquoise, black show the fractions of free helium ions n He ++ /ni at the corresponding temperatures. Right panel: Lithium plasmas. We show that the relative pressure p rel = p/(3nikBT ) at T = 50000, 80000, 120000 K is going down to about 1/4 and forms also a minimum connected with the binding of three electrons and one nucleus in a lithium atom. We included here the Fermi pressure and see an increase at high densities, which however is only qualitatively correct since there are other corrections also needed (see Section 5).
additional negative contributions to the pressure, as we will see in the last section. In order to describe the transition to full ionization in the degenerate region more correctly, these and other additional effects as the bound state shifts have to be taken into account. This leads in particular to a loss of the symmetry between electrons and ions. The bound state shifts and their influence on the disappearance of the bound states with increasing density will be discussed in the next section following references [1,2].

Low and moderate densities in helium and lithium plasmas
The calculations for helium and lithium plasmas are performed in the same way as for hydrogen except the different definitions of the Z-functions explained in Section 3.2. The result for the relative pressure of a helium plasma is shown in Figure 2. In order to compare with the result for the density of free helium ions, which is given by the Z-function, we show Z(n i , T ) by thin lines.
From the physical point of view we observe the expected behavior. Bound states formation is reflected by shoulders and minima of the relative pressure. The agreement with earlier calculations of the pressure using the mass action law [11,12,39] is, in particular around the minimum, not yet satisfactory.
We summare now the results of this section. We have shown that transitions from full to partial ionization and back to full ionization at high densities occur in the temperature range 10 4 −10 5 K in a density region around n i 10 20 −10 24 cm −3 . The formation of H atoms, also H 2 molecules, He and Li atoms may be described in the physical picture. We included ring and ladder diagrams in higher order approximations as well as certain combinations of these basic diagrams. The formation of the minima is essentially connected with the contributions of inverted ladder diagrams. The principal scheme how the formation of bound states depends on the density at a fixed temperature is shown in Figure 2. We see that the main effects as the formation of atoms and molecules and their destruction at very high densities are correctly described. The agreement with advanced calculations based on the more complicated chemical picture is, at higher densities around the minima and beyond, only qualitatively so far. The inversion of fugacity functions stemming from summations of partial fugacity series leads to quite simple descriptions of the bound state formation without using mass action law equations. At extremely high densities, the present methods fails due to the missing Hartree-Fock and Pauli blocking effects which are relevant at close packing of electrons and ions. This will be discussed next.

Effective wave equations and energy shifts
The study of level shifts and broadening due to interaction effects is of large significance for the investigation of plasmas [10,23]. According to elementary estimates, the ground state level disappears because of screening effects at r D < 2a B . Later studies give the more precise value of the Debye radius r D < 0.84a B at which the levels disappear. This idea goes back originally to Ecker and Weizel and was worked out by Theimer and Kepple and in more detail in [8,24,25,27]. The cited authors developed the general idea, that the continuum comes down with the density, so that the ionization energies are squeezed. Note that Mott-like effects may be of different origin [23]. The disappearance of the bound states with increasing density was discussed in many works [8][9][10][23][24][25]27]. The appropriate tool is the use of a generalization of the Schrödinger equation, the Bethe -Salpeter equation which is named after Hans Bethe and Edwin Salpeter. It is difficult to trace back the origin of this fundamental equation. Evidently it was first published in 1950 in the final part of a paper by Yoichiro Nambu, but without derivation and then formulated in a more systematic way by Bethe and Salpeter (1951), so one should more correctly say "Nambu-Bethe-Salpeter" equation. This equation was used in particular for the development of a more advanced theory of the energy levels in plasmas [23][24][25]27]. Nowadays this concept is found in fruitful applications in the theory of nuclear bound states [4,[40][41][42] and recently also in the physics of graphene [5].
The bound states of pairs of particles play a major role for the properties of gases and plasmas, they are responsible for the formation of atoms, molecules, clusters etc. Let us first consider the treatment of a pair bound state in the case that the pair is isolated from the surrounding. This way we have a many-body problem for N = 2 particles (e.g. electron and proton) with a central force potential U = U (r), where r = |r 1 − r 2 |. The wave function Ψ (r 1 , σ 1 , r 2 , σ 2 , t) obeys a Schrödinger equation in the 6-dimensional space After separating the center of mass motion, the Schrödinger equation for the relative motion in momentum representation reads The European Physical Journal Special Topics For the electron-proton system, the reduced mass is m −1 ep = m −1 e + m −1 p ≈ m −1 e . Using Rydberg units ( = 1, m e = 1/2, e 2 = 2), the Coulomb interaction in momentum representation reads V (q) = 8π/q 2 . The normalized wave function for the ground state (n = 1) is The corresponding ground state energy for hydrogen is 1 Ry = 13.56 eV, the higher levels depend on two quantum numbers E s,l . If the bound states are imbedded into a plasma, this shifts the levels as pointed out above.

Level shifts and disappearance of the discrete spectrum
An elementary approach is the confined atom model [24], which assumes that the bound state electrons are confined by a spherical box with the radius r 0 given by According to this estimate, the bound state levels merge into the continuum at The physical reason for this shift is the Pauli exclusion principle, i.e. the bound electrons are not admitted to penetrate the space occupied by the atoms in the neighborhood. The elementary estimate given above provides us with the condition for the existence of a ground state This condition, which is less stringent than the Mott conditions, needs certainly improvements. A more strict calculation of the level shift may be given with the quantum statistical approach based on the Nambu-Bethe-Salpeter equation [24,25], introducing concepts such as the self-energy, dynamical screening and the spectral function. For these many-body quantities, special approximations can be performed which reflect different processes in the plasma. In this way, an effective wave equation is obtained [10,24]: We assumed here the adiabatic limit m e /m p 1, the center of mass motion p has been neglected. In general, the plasma Hamiltonian H pl (q) will depend also on p and on the energy, if dynamical and retardation effects are taken into account. The plasma Hamiltonian will shift the energy eigenvalues E n = E 0 n + ∆E n and will modify the wave functions ψ n (p). With increasing density, the influence of the plasma increases and the binding energies may merge into the continuum so that bound states disappear. This breakup of bound states is called Mott effect and has important consequences for the macroscopic properties of the plasma. We focus on the influence of the mean-field contributions to the effective Schrödinger equation of pairs, − f e (p + q)ψ n (p) = E n ψ n (p) .
Here, f e (p) = f e (E p ) is the Fermi distribution. In the case of low density the perturbation due to the plasma Hamiltonian is small, the shift of the energy eigenvalues is obtained with the unperturbed wave functions as Here, inserting the Schrödinger equation, the Pauli blocking term can be rewritten as [1,2] ∆E Pauli A simple expression is found in the low-density limit, where the Fermi distribution with the normalization p f e (p) = n e /2 is concentrated near p = 0. In the zero temperature limit, we have a Fermi sphere with Fermi momentum p F = (3π 2 n e ) 1/3 1. The energy shift of the ground state φ 1 (p), equation (60), results as Here and in the following we are using Rydberg units where the dimensionless density of free electrons is n e in units of a 3 B , and the ground state energy is E 0 1 = −1. In this approximation, the Pauli blocking shift is linear in the density. In the general case of arbitrary temperatures, we have to introduce the Fermi function and the integrals in equation (16) cannot be taken in a simple way, The European Physical Journal Special Topics If we approximate the Fermi distribution by a Boltzmann distribution normalized to the same density, we obtain an analytical expression where the function G(T ) expressing the temperature-dependence is given in a simple approximation [23]. Here the Rydberg units of temperature are given by Ry/k B = 157886 K, i.e. by the ionization temperature of hydrogen. In the asymptotic approximation T 1 (i.e. temperatures below 20000 K where G(x) 1) this leads back to the zero temperature expression for the shift given by equation (17).
Similar expressions can be given for the Fock term. In the low temperature, low density limit we get It compensates partially the Pauli shift so that the total shift is ∆E Fock It is shown in Figure 3, full line, indicating a rather steep shift of the bound state energy. Taking into account the temperature effects, the Fock shift is given by Due to phase space occupation, the bound state energy is shifted and may merge with the continuum of scattering states, indicating the breakup of bound states. Considering in equation (67) the continuum part of the spectrum describing scattering states, only the Fock shift contributes to the energy shift. The lowest energy in the continuum occurs at p = 0 and is shifted by However, the two-particle continuum state can only be created at the Fermi momentum since all states below that are occupied. Thus the continuum of scattering states, where the bound states become free, begins at p F where we have in the zero temperature limit the Fock shift shown also in Figure 3. Extrapolating these low-density results to higher densities, the ground state disappears at a density which corresponds (in Rydberg units) in first approximation to n e 0.015. This leads to an average distance of r 0 2a B and is with respect to the density lower than the Mott criterion. The Mott condition expresses the idea that atoms are destroyed if the mean distance of the electrons is equal or smaller than the Bohr radius. We remember that the elementary estimate based on the confined atom model led us even to a higher limit density r 0 0.4a B . For a variational evaluation of energy shifts see [1]. In order to summarize: In this section the influence of the medium is studied in mean-field approximation, i.e., by Pauli blocking and Fock self energy. The first of the contributions to the shift we denote as Pauli-blocking and the second one as Fock contribution. We have shown above that the Fock term and the Pauli blocking contributions have opposite signs and compensate each other partially. We consider the density region where Pauli-Fock shifts dominate and remember that the Fock term and the Pauli blocking contributions have opposite signs. According to the present estimates (see Fig. 6), the hydrogen levels merge into the (lowered) continuum at a density of about n i a 3 B 0.01. This gives a proton density in H-plasmas of about 10 23 cm −3 . In He or Li plasmas, because of the smaller Bohr radii, the density of fully ionization of all atoms is about one order of magnitude higher. At these high densities, where the bound states disappear, the equation of state needs an entirely different approach, which is the topic of the next section.
6 Equation of state at high densities and pressures 6

.1 Two fluid model and Hartree-Fock-Wigner contributions
In the limit of very high densities, we may think about nuclear densities beyond 10 22 −10 23 cm −3 , we expect that the plasmas of light elements are determined mainly by the Pauli principle and by strong Coulomb repulsion of the nuclei and show therefore some universality. In this region of very high densities, the plasmas of light elements may be described by the so-called two-fluid model, as a composition of an electronic fluid, i.e., a degenerate one-component plasma (OCP) with a positive background, and a strongly interacting ionic fluid with a negative background. In the region of strongly degenerate electrons, the thermodynamics is in good approximation determined by the ideal Fermi-Dirac and by the Hartree-Fock contributions to the thermodynamic functions. Bound states are not relevant due to the fact that atomic shells are no more stable mainly due to the shifts and the Pauli blocking effect discussed above. We neglect continuum correlations in the electron-proton channel which are described by the generalized Beth-Uhlenbeck formula [42]. Correlations in the degenerate electron system are described by the Macke-Gell-Mann-Brueckner result, see [24].
We investigate now the triangular region between n e Λ 3 e ∼ 1 where the electrons become degenerate and the value of the Brueckner parameter Note that we have to generalize here the definition of the coupling parameter Γ e including the degeneration effects [20]. The electrons are degenerate and strongly coupled if Γ e > 1. The ions behave classically and form a strongly coupled classical subsystem. For the pressure, correlations in the ion subsystem lead to the Wigner contributions stemming from short-range order and lattice-like structures. Let us first consider the electronic part. In the limit of high densities the pressure is dominated by the Fermi-Dirac (FD), Hartree-Fock (HF), and Gell-Mann-Brueckner (GB) contributions. At T = 0, they are depending on the parameter r s , see Section 6.4. in [24], where the corresponding contributions to the pressure (in Ryd units) are At the very high densities studied here, the first two contributions dominate and the GB term corresponding to the ring diagrams may be omitted. At finite temperatures the Fermi-Dirac and the Hartree-Fock contributions to the chemical potential and the corresponding pressure are determined by Fermi-integrals and their inversion functions. As a result we have, see, e.g., [24] α FD = βµ FD e , n e Λ 3 e /2 = I 1/2 (α FD e ), The pressure may be calculated either by integrations or obtained directly from the general formula We study now the temperature dependence of the Hartree-Fock contributions to the chemical potential and the pressure following a procedure developed by Zimmermann [27]. First we find α HF = βµ HF e as a function of the dimensionless electron density y = n e Λ 3 e /2 where the limiting cases of small or large density respectively have to be observed. After finding piecewise representations we use the method of concatenating the pieces at the crossing points developed for the ideal Fermi contributions (1 + 0.8543y + 10y 3.666 ) , z HF e = n e exp − e 2 k B T Λ e (y − 0.2069y 2 + 3.103y 4 ) (1 + 0.8543y + 10y 3.666 ) .
By construction of the interpolating procedure these expressions are compatible with the known approximations for low or high degeneracy respectively. For alternative representations including the higher order terms and covering the full range see a review [43]. For the intermediate region E F < T < 12E F a parametrization of the Hartree-Fock contributions was given by Perrot and Dharma-wardana and discussed in [24]. The dependence of the Hartree-Fock functions on the density is demonstrated in Figure 4. We observe here a smooth transition of the expansions in the region of intermediate degeneracy. Numerical Monte-Carlo (MC) results obtained by path integral MC can be found in [44]. The contributions stemming from the ionic fluid is numerically of the same order and may be approximated in the framework of two fluid models by Wigner-type expansions [20,24]. In the density region n i 10 23 − 10 25 cm −3 which we will study here, we may use in zeroth approximation the following limits for pressure and fugacities: We see that the Wigner contributions as well as the Hartree-Fock contributions to the pressure increase strongly with the density. At the same time the fugacities strongly decrease. Both are changing their analytical shape at certain degeneracy parameter, where the non-degenerate behavior goes more or less abrupt over into the degenerate behavior, see Figure 5.
In the higher approximations, the ionic contribution may be approximated at high densities by Wigner type expansions [20,24,45]. In the density region which we will study here we may use for y i = n i (z 2 i e 2 /k B T ) 3 1 an analytical approximation for the ionic pressure stemming from an interpolation of Monte-Carlo calculations [11,12,24,45] In difference to the situation at small densities, the contribution of the bound states is in the high-density region just a perturbation which sometimes may be neglected in the limit of very high densities and pressure. Let us make a final note on this section. Hartree-Fock and Wigner effects are strictly speaking not closely related but are stemming from quite different order contributions in the interaction parameter e 2 . However both contributions are similar with respect to their order in 1/r s and are numerically of similar order. e /2 for the interaction parameter ξ = βe 2 /Λe = 2 (above) and ξ = 8 (below). Right panel: Fugacity coefficient of the ions due to Wigner effects in dependence on the ion interaction strength yi = ni(βe 2 ) 3 for rs = 0 (above) and rs = 1 (below).

Bound state effects at high degeneration
In order to include the bound state contributions which are beyond Hartree-Fock approximations we represent the contributions to the pressure beyond the Hartree-Fock term by two-particle Greens functions where G 2 (1, 2, 1 , 2 ) is the full two-particle Green's function which may be represented by Feynman diagrams. The contributions of the omitted first diagrams have already been calculated below. For the treatment of the light elements up to lithium we have to include all Feynman diagrams including one heavy particle and 1 to 3 electrons. This leads to the formula (93) for the pressure in the grand canonical ensemble which should correspond to the low density representations given in Section 3. We assume in the following that at very high densities / high degeneration the bound state effects are just a correction. Note that beyond the Mott density we still may have relevant scattering state contributions [42]. Here we assume that in a first order the previous formulae are still valid with the difference that the fugacities are determined by the Hartree-Fock and Wigner expressions. This way we get, e.g., for the pressure of hydrogen the estimate (compare with Sections 1-3 and note that the index "high" stands here for "high density", and B for Boltzmann) For helium we have correspondingly and for and the grand-canonical pressure of lithium we get Here the coefficientsb denote the bound state virial contributions with energy levels from the HF solution of the Bethe-Salpeter equation.

Applications to plasmas at extreme densities and pressure
The theory presented in the first sections provides isotherms of the pressure for the non-degenerate region. After crossing the line of degeneracy n e Λ 3 e 1 to the degenerate regions above it, the previous theory is only a rough approximation since degeneracy effects become dominant. Further the dissymmetry of the masses maybe relevant. In this region the formulae given above are just an extrapolation which is exact only with respect to the (sometimes dominant) ideal Fermi contributions. For very dense, strongly degenerate plasmas, bound states do not exist, the plasma behaves like a degenerate non-ideal Coulomb gas. We will study now in more detail the region where no bound states exist any more, according to the solutions of the Nambu-Bethe-Salpeter equations. According to the estimates from the previous section the main contribution to he EOS comes now from Pauli and Fock effects. For hydrogen at high densities we may use equation (96). This approximation and the one at low density p low are crossing in the region of about n i = 10 23 cm −3 (see Fig. 6). A convenient interpolation is obtained here again by using the function w(x) = tanh(x/2)/2 which connects the two pieces in a smooth at the crossing point y 0 18 where y = n e Λ 3 e /2. This way we get in our case with a 0.1 p low/high = p low /(1 + exp(a(y − y 0 ))) + p high /(1 + exp(−a(y − y 0 ))).
The present theory may be considered as a simple extension of the two approaches given above to the whole range of densities. We demonstrate the method of of smooth concatenation of low density and high density approximations using tanh(x) functions in Figure 6. This way we found now for the limit of high density closed expressions for the pressure of hydrogen and other light plasmas expressed in the physical picture of densities. Basically we have developed now by inversion of the fugacity series different formulae for the region of low, i.e. small and moderate, densities including bound states on one side and large densities on the other side neglecting bound states. Both limits are simple composed here, so that some uncertainties in the transition region remain to be discussed.
As an application we calculated the pressure of hydrogen related to the classical ideal pressure for the case of high densities beyond n ∼ 10 23 cm −3 at temperatures 10000−50000 K using the extended Hartree-Fock and Wigner approximations given above, see Figures 7 and 8. The agreement of our new, rather simple high-density approach with the previous ones [13,[46][47][48][49] is satisfactory. A comparison with data from the numerical approach VASP [46][47][48][49] as well as other methods discussed in [46,50] shows that the curves for the relative pressure obtained by our method for the region 10 23 −10 25 cm −3 are about 10 percent too high. On the other hand, our analytical formulae allow the calculation of a density series at given temperature in Fig. 8. Calculations of the high-density pressure of hydrogen plasmas based on the extended Hartree-Fock-Wigner approximation. The left panel shows the relative pressure calculated using our low-density approach for the density region below 10 23 cm −3 at temperatures T = 20000, 30000, 50000 K. Note that the bound states are in the region below n ∼ 10 23 cm −3 described only in certain approximation. At densities above n ∼ 10 23 cm −3 we used the approximations developed in this section. In the right panel we show the pressure related to the classical ideal pressure at the temperatures T = 10000 (in blue) and 20000 K (in red) in high-density approximation. We present also a comparison with a curve (in green) from the PACH theory by Beule et al. [13] for T = 10000 K (PACH -Padé approximations and chemical picture). just a few minutes on a standard personal computer and can easily be extended to other light elements and to mixtures of them like in the atmospheres of stars.
In certain region of the density-temperature plane the pressure at high densities is essentially determined by the Hartree-Fock approximation for the thermodynamic functions and the Hartree-Fock approximation and Pauli blocking for the energy shifts. We assumed here that in the high-density region these shifts provide the most important effects for the destruction of bound state. In particular we contribute here to the theory of hydrogen at high pressures in the region where a Mott transition to full ionization has been predicted and where recent experiments have shown a transition from insulating behavior to metal-like conductivity. In order to understand this transition, several effects have to be taken into account. We concentrated here on so-called Pauli blocking effects expressing the rule that states occupied by atomic electrons cannot be occupied by free electrons with the same spin state. This leads at high electron densities to the destruction of atomic states which need a relatively high amount of phase space. We calculated the energy shifts due to Pauli effects and discuss the Mott effects solving effective Schrödinger equations for strongly correlated systems. A few remarks: (i) For simplicity we neglected the formation of molecules which play a role at temperatures below 20000 K [13].
(ii) In the temperature region which we studied here (10000−50000 K), no first order phase transitions were observed in this approximation. (iii) The approach based on the physical picture which we used in the present section, can be extended to lower temperatures, where H 2 molecules dominate. In the high density region studied here, the Fermi-and the Fock contributions are most essential and the present picture is appropriate and quite simple. In Figure 9 we show the corresponding results for helium and lithium plasmas.

Conclusions
In the present work, we developed the method of inversion of fugacity representations of the pressure with applications to molecular hydrogen plasmas and atomic helium and lithium plasmas at high pressures/densities. We studied the EOS of various plasmas at constant temperature for increasing density. Typically for all, graphs start with a Boltzmann behavior i.e. βp/(n e + n i ) 1 at low densities. With increasing density the pressure related to the Boltzmann pressure decreases first due to the formation of atoms and molecules. Then it goes through a minimum and rises steeply due to Fermi, Hartree-Fock, and Wigner effects. We describe the formation of the atomic and molecular minima without using mass action laws, using instead algebraically inverted fugacity expansions up to four-particle cluster integrals. At very high densities the EOS is determined beside the ideal Fermi pressure by electronic Hartree-Fock and ionic Wigner-DeWitt approximations. This leads to a strong increase of the pressure which is in the high-density limit quite universal. The methods presented here lead to analytical approximations which are easily programmed 3430 The European Physical Journal Special Topics on standard personal computers and provide within seconds of calculations curves for the pressure for hydrogen, helium and lithium. The approach is applicable to any mixture of plasmas of the light elements. The present method does not include the region of extreme, beyond Megabar, pressures, where according to estimates [11], plasma phase transitions may occur. Fugacity expansions are always rather smooth and do not show van der Waals wiggles, so that one may have difficulties to detect phase transitions by standard methods. However the approach is in principle applicable to any mixture of plasmas of the light elements possibly up to carbon and to conditions typical for some astrophysical applications.
Open Access funding enabled and organized by Projekt DEAL. The authors acknowledge fruitful discussions with Dietmar Ebert, Wolf-Dietrich Kraeft, Ronald Redmer, Sergey Trigger and Jan Vorberger. The work of G.R. was partially supported by the MEPhI Academic Excellence Program under contract number 02.a03.21.0005.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.