Exact asymptotic expansions for thermodynamics of the hydrogen gas in the Saha regime
We consider the hydrogen quantum plasma in the Saha regime, where it almost reduces to a partially ionized atomic gas. We briefly review the construction of systematic expansions of thermodynamical functions beyond Saha theory, which describes an ideal mixture of ionized protons, ionized electrons and hydrogen atoms in their ground-state. Thanks to the existence of rigorous results, we first identify the simultaneous low-temperature and low-density limit in which Saha theory becomes asymptotically exact. Then, we argue that the screened cluster representation is well suited for calculating corrections, since that formalism accounts for all screening and recombination phenomena at work in a more tractable way than other many-body methods. We sketch the corresponding diagrammatical analysis, which leads to an exact asymptotic expansion for the equation of state. That scaled low-temperature expansion improves the analytical knowledge of the phase diagram. It also provides reliable numerical values over a rather wide range of temperatures and densities, as confirmed by comparisons to quantum Monte Carlo data.
Obtaining asymptotically exact formulae for the equation of state of quantum Coulomb matter is important, both at a theoretical level and for practical applications. They provide a better understanding of basic phenomena like molecular recombination and screening in the framework of statistical mechanics. Such formulae are free from any intermediate phenomenological modelization and uncontrolled approximation. They provide moreover reliable and accurate data in some range of thermodynamical parameters. Exact expansions are of particular interest for hydrogen described as a gas of quantum protons and electrons interacting via the Coulomb potential. Indeed, thanks to its relative simplicity, analytical calculations can be carried out further than for heavier species. In practice, the corresponding expansions are quite useful since hydrogen is the most abundant element in the universe. Astrophysicists need accurate equations of state over a wide range of temperatures and densities, including the so-called Saha regime where hydrogen reduces to a partially ionized atomic gas.
In the ’s , a first kind of asymptotic expansion was derived for an electron gas at high densities, which behaves as a free Fermi gas in a first approximation. Corrections can be computed in a systematic way in the framework of standard many-body perturbation theory , where the small-expansion parameter is the charge of the electrons. A second kind of asymptotic expansions, namely the familiar virial expansions, where constructed in the ’s , for two- or more component systems, including quantum hydrogen. In such expansions, temperature is fixed at a non-zero value, and numerical densities for species are driven to zero. For hydrogen, we set . At lowest order, the system behaves as an ideal mixture of nuclei and electrons, as rigorously proved in Ref. . Corrections are represented by series involving integer and half-integer powers of the ’s, as well as logarithmic terms. Corresponding calculations have been first performed up to order , by using the effective-potential method introduced by Morita . Further corrections of order have been derived in the ’s within another formalism based on the path integral representation , and retrieved later by Ebeling-Morita’s method .
Above expansions are suited for regimes where the systems are almost fully ionized. They cannot describe the Saha regime for hydrogen, where a finite fraction of electrons and protons recombine into atoms in their groundstate. Then, according to the familiar Saha theory , the system is expected to behave as an ideal mixture of ionized protons, ionized electrons and atoms . The construction of suitable asymptotic expansions requires first to identify, if it exists, a regime of thermodynamical parameters where Saha theory becomes asymptotically exact. As described in Section 2, that regime is obtained by setting , while ratio is kept fixed with temperature-dependent density given by expression (4).
Once the proper limit which defines the Saha regime has been identified, the construction of systematic expansions beyond Saha theory requires a formalism which accounts for all recombination and screening phenomena at work. For that purpose, the screened cluster diagrammatical representation  is particularly adequate, as briefly described in Section 3. The corresponding analysis of all involved graphs, provides the so-called SLT expansion of pressure , i.e.
where Saha pressure is given by formula (5) in section 2. Functions only depend on ratio , while temperature dependent functions decay exponentially fast when vanishes, except for possible multiplicative powers of . Expansion (1) is ordered with respect to the decaying rates of the ’s functions. The first five corrections computed in Ref.  are schematically presented in Section 4. They account for non-ideal phenomena such as plasma polarization, shift in the atomic energy levels, interactions between ionized charges and atoms, and also formation of molecules or ions and . Further corrections decay exponentially faster than , where is the atomic groundstate energy and is the reduced mass for the two-body electron-proton problem. Along a given low-temperature isotherm, we also study the behaviour of the ’s. SLT expansion (1) then appears as a partial infinite resummation of ordinary virial expansions at low densities . Expansion (1) remains valid at intermediate () and large () densities, but breaks down at too large densities because molecular recombination becomes then prominent.
As usual for asymptotic series, keeping only the first few terms of SLT expansion (1) should lead to an accurate equation of state, provided that thermal energy is smaller than Rydberg energy . In Section 5, we give a flavour of numerical calculations based on the truncation of (1) up to term included . Because of the relatively large temperature scale , and of the occurrence of exponentially decaying factors, that truncated equation of state is reliable over a rather wide range of thermodynamic parameters, as confirmed by comparisons to quantum Monte Carlo simulations by Militzer and Ceperley . This allows us to introduce a semi-empirical criterion for the convergence of SLT expansions, which provides the validity domain of the corrected EOS in the temperature-density plane. Our formulae should be particularly useful in physical situations where deviations to Saha theory play an important role even if they remain small. For instance, a very accurate EOS is needed for interpreting recent seismology measurements in the Sun . Our corrected EOS should be useful since, according to usual models, the Sun adiabat lies in the previous validity domain (see Fig. 1). However, notice that an accurate description of that adiabat requires to take into account heavier species like helium, carbon, nitrogen and oxygen. If exact calculations become much more complicated, a simple account of obvious ideal contributions, in particular for less abundant species, should be sufficient for significantly improving a pure-hydrogen EOS.
We stress that our approach does not provide unambiguous definitions of neither free and bound charges, nor ionization rates. As commented in Section 4, this does not cause any trouble as far as thermodynamic quantities are concerned. Other quantities like conductivity or opacity cannot be obtained within the present formalism. Usually, such quantities are computed within the framework of the chemical picture, where atoms and molecules are introduced phenomenologically as preformed entities. A first-principles description is of course possible in principle, but it becomes quite cumbersome since it would require a rigorous introduction of either real-time evolution or coupling to radiation.
2 The hydrogen gas in the Saha regime
2.1 The physical picture
Within the physical picture, a hydrogen gas is viewed as a system of quantum point particles which are either protons or electrons, interacting via the instantaneous Coulomb potential . Protons and electrons have respective charges, masses, and spins, and , and , . In the present non-relativistic limit, the corresponding Hamiltonian for particles reads
where is the species of the th particle and is the Laplacian with respect to its position . The system is enclosed in a box with volume , in contact with a thermostat at temperature and a reservoir of particles that fixes the chemical potentials equal to and for protons and electrons respectively. Because the infinite system maintains local neutrality in any fluid phase, the bulk equilibrium quantities depend in fact solely on the mean
while the difference is not relevant as rigorously proved in ref. .
2.2 Identification of the scaled limit
In the so-called Saha regime, a finite fraction of protons and electrons combine into hydrogen atoms, forming a partially ionized hydrogen gas. That regime is attained when several conditions are met. The temperature must be sufficiently low so that atoms can form, namely . The density must be sufficiently low as well so that atoms maintain their individuality thanks to , where is the mean inter-particle distance and is the Bohr radius. If the density becomes too low, atoms dissociate by entropy, while if it becomes too large, they recombine into molecules . According to those simple considerations, both and must go to zero, in a related way. The precise form of that relation can be inferred from a rigorous analysis in the grand-canonical ensemble devised by Macris and Martin , who extended Fefferman’s work on the atomic phase of the hydrogen plasma . They introduce a scaling limit where the temperature goes to zero, while the average chemical potential of protons and electrons approaches the ground-state energy with a definite slope . Then, they proved that pressure , within that scaling limit, tends to its Saha expression , which describes an ideal mixture of hydrogen atoms, ionized protons and ionized electrons.
In terms of temperature and density, the previous scaling limit can be rephrased as a low temperature expansion at fixed ratio [see eq. (1)], where is the temperature-dependent density
Notice that density vanishes exponentially fast when is sent to zero. This ensures the proper energy-entropy balance which keeps a finite ionization rate that is entirely determined by the fixed ratio . Pressure in units of tends to Saha formula
apart from exponentially small terms when . The respective behaviours of Saha pressure (5), for , and for clearly illustrate that is a cross-over density between full ionization and full recombination.
3 Construction of SLT expansion
3.1 Introduction of a suitable formalism
Corrections to Saha pressure (5) involve interactions between ionized charges and atoms, as well as formation of ions and molecules. Standard many-body theory is not well-suited for taking into account recombination, since that mechanism is not perturbative with respect to the charge. For instance, an infinite number of Feynman ladder graphs must be resummed for describing a single atom .
The effective-potential method, which amounts to introduce a classical equivalent system of point particles with many-body effective interactions, is a priori more efficient for dealing with recombination. Indeed, -body effective interactions are inferred from -body quantum Gibbs factors which do account for recombination of particles at short distances. The contributions of two-body effective interactions can be analyzed within standard methods of classical statistical mechanics. In particular, that feature has been exploited for computing virial expansions up to order , where atomic contributions appear. Ionic and molecular recombination are embedded in three- and four-body effective interactions. Unfortunately, the analysis of the corresponding contributions becomes rather cumbersome, in particular because no standard classical tool is available.
Above drawbacks of both standard many-body theory and effective-potential method, clearly emphasize the need for a formalism more appropriate to deal with Saha regime. The ACTEX method, introduced by Rogers  in the 70’s, is intended to account for the formation of chemical species in the framework of the physical picture. That approach starts from the usual activity expansion of thermodynamical quantities in the grand-canonical ensemble. Despite its rather successful predictions at moderate densities and temperatures, that approach cannot be applied here as it stands, because ACTEX series are not exactly reorganized via a systematic treatment of both recombination and screening. In fact, quantum Gibbs factors involved in activity series do not factorize as products of two-body counterparts, like in the case of classical charges. Resummations, which are crucial for taking into account screening effects, are very hard to handle in the usual quantum activity series.
The difficult task of controlling simultaneously recombination and screening effects in quantum activity series can be accomplished thanks to the Feynman-Kac path integral representation. Within that formalism, the genuine quantum system of point protons and electrons is shown to be equivalent to a classical gas of extended loops . Thermodynamic quantities of hydrogen are then represented by activity series in the world of loops, which can be suitably rearranged as described below.
3.2 The screened cluster representation
Since loops are classical objects with two-body interactions, standard Mayer diagrammatical methods can be applied. In particular, the one-body loop density is represented by a series of Mayer-like graphs in the grand-canonical ensemble. Usual points are replaced by loops, loop fugacities are simply related to particle fugacity , while Mayer bonds are built with the loop-loop interaction potential. Since that potential behaves as the Coulomb interaction at large distances, Mayer graphs are plagued with long-range divergences. Such divergences are systematically removed via chain resummations, which amount to introduce a screened potential . Contrary to the familiar classical Debye potential which decays exponentially fast, that quantum potential decays only as at large distances . However, at low densities, it reduces to its classical Debye counterpart plus small corrections. At the same time, the whole series is reorganized in terms of particle clusters. Eventually, particle density , obtained by integrating loop density over all possible shapes, is exactly rewritten as the following diagrammatical series 
The corresponding graphs are constructed with topological rules close to that of ordinary Mayer graphs, except for some exclusion constraints avoiding double counting. Usual points are now replaced by particle clusters. The statistical weight of a given cluster involves particle fugacities, as well as screened interactions. Two clusters can be connected by a single bond, which is either , , or , where is the screened interaction between those clusters.
Expansion (6) accounts in a fully consistent way for all effects of interactions in the system at finite density and finite temperature. The various phenomena at work are embedded in well-defined graphs. For instance, the first graphs shown in (6) account respectively for a single ionized proton, formation of an atom and of a molecule , and interactions between two atoms. As the scaling limit defines quite diluted conditions, only a few simple graphs in the screened cluster expansion are expected to contribute to the first corrections to Saha theory.
3.3 Behaviour of graphs in the scaling limit
The scaling limit defined in Section 2.2 within the grand-canonical ensemble, can be rephrased as follows. Starting from particle fugacity and temperature , we introduce a new couple of independent thermodynamic parameters defined through the relation equal a constant times . Then, we set at fixed . The behaviours of graphs in representation (6) result from the competition between three mechanisms, which can be roughly described as follows.
Screening : Contributions of bonds , and are controlled by the inverse Debye screening length for ionized protons and ionized electrons with density . They behave as positive or negative powers of , which itself decays exponentially fast as .
Recombination : For a particle-cluster made with protons and electrons, the behaviour of its statistical weight gives raise to a cluster partition function in the vacuum. Each is a truncated trace of Gibbs operator involving only bare Coulomb Hamiltonians, which converges thanks to a truncation inherited from screening by ionized charges . In , contributions of all possible recombined entities made with protons and electrons, are mixed together. Remarkably, the contribution of a given chemical species made with protons and electrons, naturally emerges through Boltzmann factor , where is the groundstate energy of Hamiltonian . That factor increases exponentially fast since .
Entropy : In a given cluster, the presence of particles generates activity powers , which decay exponentially fast as .
The behaviour of a graph in the scaling limit is obtained, roughly speaking, by taking the product of the exponential factors generated by each of the above mechanisms. Then, every graph is found to decay exponentially fast. The leading contributions arise from the first two graphs in representation (6). They do reduce to the ideal terms and predicted by Saha theory. Further corrections decay exponentially faster than in agreement with rigorous bounds . Dividing all terms by , the SLT expansion of reads 
where the first two terms account for contributions from ionized particles and hydrogen atoms in their ground state. The remaining terms involve functions that decay exponentially fast when , while is an integer or half-integer power of (which may be multiplied by integer powers of when ). Expansion (7) is ordered with respect to increasing decay rates of the ’s. The corresponding hierarchy follows from subtle inequalities between the groundstate energies of all Coulomb Hamiltonians . For instance, the molecular contribution, which determines the leading low-temperature behaviour of , indeed decays exponentially fast thanks to : this ensures that molecules are very scarce in the Saha regime compared to atoms , despite they are more stable energetically, i.e. .
4 Equation of state beyond Saha theory
4.1 Scaled low-temperature expansion of pressure
Representation (7) expresses the density in terms of variables and , or equivalently and since there is a one-to-one correspondence between those sets of variables. As the natural thermodynamical parameters are the temperature and the density, it is quite useful to invert the SLT expansion (7) to determine , in order to compute all thermodynamical quantities as functions of and . In the present scaling limit, this can be done in a perturbative way around Saha expression
easily obtained by keeping only the first two terms in (7). Each correction to that form reduces to a product of an algebraic function of times a temperature-dependent function which decays exponentially fast.
The standard thermodynamical relation, which expresses density in terms of the partial derivative of pressure with respect to at fixed , can be rewritten here as
After inserting SLT expansion (7) of into (9), a straightforward integration with respect to provides the SLT expansion of in terms of and . The corresponding expansion (1) of in terms of and , then follows by using the inversion relation determined above. The physical content of first five corrections in (1), as well as the expressions and values of the corresponding decay rates are summarized in the following table
|Correction ()||Physical content||(in eV)|
|1||plasma polarization around ionized charges|
|2||formation of molecules, atom-atom interactions|
|3||atomic excitations, charge-charge interactions|
|4||formation of ions, atom-charge interactions|
|5||fluctuations of plasma polarization|
First correction is equivalent to a modification of Saha ionization equilibrium  derived within Green functions techniques (see also Ref. ), where rate arises from the behaviour in the scaling limit. All further corrections are entirely new, as well as the structure of SLT expansion (1). Beyond their leading behaviours , functions include further corrections which decay exponentially faster. For instance, if the leading behaviour of is controlled by the first atomic excited state, all the other contributions of excited states are incorporated into . Similarly, includes not only the contribution of the molecular groundstate, but also all contributions of molecular excited states. As mentioned above about recombination, the sums of all those contributions are indeed finite. Also, we stress that such contributions of recombined entities, like atoms in or molecules in , are entangled with that of their dissociation products. Thus, purely atomic or molecular contributions cannot be unambiguously defined. This does not any trouble here, since only full contributions embedded in and are relevant for thermodynamics. In other approaches based on the chemical picture, that ambiguity has been the source of many controversies since the introduction of Planck-Larkin formula (see e.g. Ref. , Ref.  and Ref. ). Eventually, notice that the contributions to expansion (1) of more complex entitites, like , or , decay exponentially faster than as detailed in Ref. .
4.2 Low-temperature isotherms
Let us consider now a small fixed temperature , and study the behaviour of various corrections to Saha pressure in (1) when density is varied. The corresponding low () and large ( ) density behaviours are summarized below :
At low densities , the leading correction of order is given by plasma polarization (term ), while at large densities the leading correction of order arises from both molecules and atom-atom interactions (term ). We have checked that the familiar virial expansion at low densities is indeed recovered up to order included, since all terms provide powers higher than . We stress that at too large densities, expansion (1) breaks down because various corrections, in particular those due to molecular recombination, prevail over Saha pressure which grows only as . In fact, a semi-empirical criterion based on that observation allows us to infer a validity domain for the SLT expansion, as described below.
5 Numerical applications and comparisons to Monte Carlo data
The truncated EOS obtained by keeping the first five terms in SLT expansion (1) can be computed numerically. The ’s are easily computed since they reduce to simple algebraic functions of ratio . Numerical values for temperature-dependent functions , and are also readily inferred from explicit analytic expressions. No similar expressions for and are available, since analytical results on the three- and four-body quantum problem are very scarce. Then, we use simple modelizations of those functions which account for their exact low-temperature forms on one hand, and incorporate familiar phenomenological descriptions of ions and and of molecule on the other hand .
Various isotherms corresponding to increasing temperatures have been considered. For temperatures below , the Saha regime defines extremely diluted conditions which do not make physical sense : this explains why Earth or Brown Dwarfs atmospheres only involves molecules . Temperature can be increased up to , which is still small compared to the characteristic temperature scale . At a given temperature, calculations within the truncated EOS are not reliable above some density , for which corrections to Saha pressure become too large, in agreement with previous estimations for . That breakdown is due to molecular recombination below K, and to atom-atom interactions for higher temperatures.
At relatively low temperatures, i.e. below , PIMC calculations  have been performed at rather high densities for which atoms are mainly recombined into molecules . In the Saha regime, statistics in PIMC results are poor because the corresponding densities are too diluted. In fact, under such conditions, our analytical results might serve as testbenchs for simulation methods. For higher temperatures, i.e. above , there exists a density range (see PIMC crosses displayed in Fig. validitydomain), where comparisons between our results and PIMC data  are instructive. A good agreement is observed in some density range, which can be inferred from a semi-empirical criterion : corrections cannot exceed a few per cent of Saha pressure. This defines, at a heuristic level, the validity domain of SLT expansion shown in Fig. validitydomain. In that whole domain, weak-coupling and weak-degeneracy conditions are fullfilled. The tongue structure of the domain between K and K results from the increase of the strength of interactions between ionized charges and atoms. Notice that the whole domain is restricted to rather low densities in general, so high-density phenomena, like the celebrated plasma phase transition, remain beyond the scope of our approach.
Eventually, we emphasize that our numerical calculations will be detailed in a forthcoming paper , where simple representations of functions () will be given, while the corresponding functions can already be found in Ref.. Our results will be compared to PIMC data, and also to phenomenological calculations. Applications to the Sun adiabat should be considered later.
-  M. Gellmann and K.A. Brueckner, Correlation energy of an electron gas at high density, Phys. Rev. 106 :364-368 (1957); E.W. Montroll and J.C. Ward, Quantum statistics of interacting particles: general theory and some remarks on properties of an electron gas, Phys. Fluid 1 :55 (1958)
-  A.L. Fetter and J.D. Walecka, Quantum Theory of Many-Particle Systems (McGraw-Hill, San Francisco, 1971). Ph.A. Martin and F. Rothen, Many-Body Problems and Quantum Field Theory: An Introduction, 2nd ed (Springer, 2004).
-  W. Ebeling, Ann. Phys. Leipz. 19 :104 (1967)
-  J. Lebowitz and R. Pena, Low density form of the free energy of real matter, J. Chem. Phys. 59:1362-1364 (1973)
-  W.D. Kraeft, D. Kremp, W. Ebeling, and G. Ropke, Quantum Statistics of Charged Particle Systems (Plenum Press, New York, 1986)
-  T. Morita, Equation of state of high temperature plasma, Prog. Theor. Phys. 22 :757 (1959)
-  A. Alastuey and A. Perez, Virial expansion of the equation of state of a quantum plasma, Europhys. Lett. 20 :19-24 (1992)
-  T. Kahlbaum, The quantum-diffraction term in the free energy for Coulomb plasma and the effective-potential approach, J. Phys. IV 10 (P5):455 (2000)
-  M. Saha, Philos. Mag. 40:472 (1920)
-  A. Alastuey, V. Ballenegger, F. Cornu and Ph.A. Martin, Screened cluster expansions for partially ionized gases, J. Stat. Phys. 113:455-503 (2003)
-  A. Alastuey, V. Ballenegger, F. Cornu and Ph.A. Martin, Exact results for thermodynamics of the hydrogen plasma: low-temperature expansions beyond Saha theory, J. Stat. Phys. 130:1119-1176 (2008)
-  A. Alastuey, V. Ballenegger and F. Cornu, submitted for publication
-  B. Militzer and D.M. Ceperley, Path integral Monte Carlo simulation of the low-density hydrogen plasma, Phys. Rev. E 63:066404 (2001)
-  W. Däppen, The equation of state for the solar interior, J. Phys. A: Math. Gen. 39: 4441 (2006)
-  E.H. Lieb and J. Lebowitz, The constitution of matter : existence of thermodynamics for systems composed of electrons and nuclei, Adv. Math. 9: 316-398 (1972)
-  N. Macris and Ph.A. Martin, Ionization equilibrium in the proton-electron gas, J. Stat. Phys. 60:619-637 (1990)
-  C. Fefferman, The atomic and molecular nature of matter, Rev. Math. Iberoamericana 1:1-44 (1985)
-  F.J. Rogers, Statistical mechanics of Coulomb gases of arbitrary charge, Phys. Review A 10: 2441 (1974)
-  J. Ginibre, Some applications of functional integration in statistical mechanics, in Statistical mechanics and quantum field theory, C. DeWitt and R. Stora, eds, Les Houches (Gordon-Breach, 1971); F. Cornu, Correlation in quantum plasmas : Resummations in Mayer-like diagrammatics, Phys. Rev. E 53 :4562 (1996); P.A. Martin, Quantum Mayer graphs: applications to Bose and Coulomb gases, Acta Phys. Pol. B 34 :3629 (2003)
-  V. Ballenegger, Ph.A. Martin and A. Alastuey, Quantum Mayer graphs for Coulomb systems and the analog of the Debye potential, J. Stat. Phys. 108 :169-211 (2002)
-  D. Kremp, W.D. Kraeft and A.J.M.D. Lambert, Equation of state and ionization equilibrium for non-ideal plasmas, Physica A 127:72-86 (1984)
-  I.L. Iosilevski and V.K. Gryaznov, Comparative accuracy of thermodynamic description of properties of a gas plasma in the Thomas-Fermi and Saha approximations, High Temp. 19:799 (1981)
A.S. Kaklyugin and G.E. Norman, Thermodynamic functions (analytic expressions) for partially ionized strongly coupled plasmas. Bound and free states contributions, Journal de Physique IV 10(P5):153 (2000)
-  A.N. Starostin and V.C. Roerich, Bound states in non-ideal plasmas : formulation of the partition function and application to the solar interior, Plasma Sources Sci. Technol. 15: 410-415 (2006); Equation of state of weakly non-ideal plasmas and electroneutrality condition, J. Phys. A : Math. Gen. 39: 4431-4439 (2006)