© The Authors 2025

1 Introduction

Understanding ionisation balance – including the population composition of different charge states – as well as the average degree of ionisation and the equation of state (EOS) of dense carbon plasmas, is important for astrophysics (Dufour et al. 2007; Fontaine et al. 2008; Giammichele et al. 2018), inertial confinement fusion (ICF) (Hurricane et al. 2023), and many other related research fields. White dwarfs of the hot DQ class are characterised by a highly dense carbon and oxygen core surrounded by an envelope primarily composed of carbon (Dufour et al. 2007; Fontaine et al. 2008; Giammichele et al. 2018). The physical properties of EOS and opacity are the direct inputs of star structure and evolution and therefore affect the inferred surface abundances and pulsation properties of this type of white dwarfs (Dufour et al. 2007; Fontaine et al. 2008; Giammichele et al. 2018). White dwarf stars are excellent laboratories for studying the physics of matters under extreme dense conditions. However, accurately modelling pulsating stars remains difficult even at a level of precision asteroseismology, due to insufficient understanding of the microphysics (Winget & Kepler 2008) of the stars. Even for the most studied star, the Sun, the predicted internal structure disagrees with helioseismic observations (Bahcall et al. 2004). This discrepancy needs to be resolved using more accurate opacity (Basu & Antia 2008; Serenelli et al. 2009). However, the accuracy of opacity depends sensitively on the accuracy of the population composition of different charge states in the plasmas.

The EOS is closely related to the ionisation balance. However, accurate determination of the ionisation balance of dense plasmas is challenging both experimentally and theoretically (Kritcher et al. 2020; Kraus et al. 2016; Zeng et al. 2023). In ICF experiments, hydrodynamic instability is a key factor in achieving ignition (Abu-Shawareb et al. 2022), and its growth has been demonstrated to be sensitively dependent on the accuracy of EOS (Clark et al. 2010; Hammel et al. 2010). Most recently, the relationship between pressure and density along the principal shock Hugoniot of hydrocarbon was experimentally measured to a high precision within 5%, which is relevant for the EOS at a pressure range of 100–450 Mbar (Kritcher et al. 2020). However, these experiments were limited to mass densities below or around 10 g cm⁻³, although experiments on ionisation balance are extended to extremely dense Be plasmas at a mass density of 55 g cm⁻³ using the X-ray scattering technique (Döppner et al. 2023). Moreover, an experimental platform for the National Ignition Facility that uses spherically converging shock waves was developed for absolute EOS measurements along the principal Hugoniot (Döppner et al. 2018).

Theoretically, a realistic description of ionisation balance for extremely dense plasmas is difficult due to the complex interactions among the particles. Thus, dense carbon plasmas have attracted a great deal of interest in this area of study over the past few years. As current predictions of the ionisation degree by widely used analytical models have been challenged by a series of experiments, an approach of density functional theory molecular dynamics (DFT-MD) simulations has been proposed to calculate the ionisation degree of carbon plasmas directly from dynamic electrical conductivity using the Thomas-Reiche-Kuhn sum rule at Gbar pressures (Bethkenhagen et al. 2020). Exploration of mean ionisation states was carried out using Kubo-Greenwood conductivity based on the average atom (AA) model (Callow et al. 2023). Considering the inconsistencies between the predictions obtained by the AA models and by the electronic transport properties (Bethkenhagen et al. 2020), Clerouin et al. (2022) proposed two definitions of ionisation based on the plasma structure as given by the pair distribution function and on the total pressure. They show the equivalence of these two definitions on dense carbon plasmas along the 100 eV isotherm and the 10 g cm⁻³ isochor. Ovechkin et al. (2023) investigated plasma ionisation balance and average ion charges in the chemical picture based on the AA model. Gawne et al. (2024) sought a consistent and unambiguous definition for the ionisation state using two potential definitions, with one based on counting the number of continuum electrons and the other on the optical conductivity by finite-temperature DFT calculations. For the ionisation state of carbon in CH plasmas, differences of up to 10% are found between the two methods. The physical effects of Pauli blocking (PB) of bound-states on the ionisation potential depression (IPD) and the average degree of ionisation were studied using a many-body quantum statistical approach for dense carbon plasmas at extreme densities up to 40 g cm⁻³ (Röpke et al. 2019). Recently, an approach based on first-principles quantum Monte Carlo simulations was proposed to investigate the IPD and Fermi barrier in warm dense matter (Bonitz & Kordts 2025).

The EOS of carbon was investigated theoretically in a density range of 3–20 g cm⁻³ by combining the DFT, the liquid free energy model, and the path integral quantum Monte Carlo calculations (Benedict et al. 2014). These investigations into dense carbon plasmas are likewise important for carbon-hydrogen mixtures (Zhang et al. 2018; Hollands et al. 2020). The EOS, transport, and optical properties of dense polystyrene for practical applications can be found in a recent review paper (Hu et al. 2018). Determining the exact plasma conditions at which the transition of K-shell delocalisation occurs, as well as the magnitude of K-shell modifications preceding it, remains challenging both experimentally and theoretically (Döppner et al. 2023). Later, Bethkenhagen et al. (2023) studied the EOS, Hugoniot, pair distribution functions, and diffusion coefficients for carbon plasmas near the principal Hugoniot at densities spanning from 8 g cm⁻³ to 16 g cm⁻³ and temperatures ranging from 100 kK to 10 MK, by employing the Kohn-Sham DFT-MD method.

Herein, we extend the study of the ionisation balance of carbon plasma from a density of ~10 g cm⁻³ (Zeng et al. 2020a) to a regime of K-shell delocalisation at extremely dense conditions with a pressure reaching a few 100 Gbar. Based on the minimisation of the Helmholtz free energy of a dense plasma, a unified and self-consistent theoretical formalism for ionisation balance and EOS is constructed, in which all the physical quantities required to determine the population distribution of different charge states in the ionisation balance are obtained self-consistently using the same plasma screening potential. The non-ideal effects on the atomic properties and on the ionisation balance are taken into account by properly considering the interactions between the electrons and ions, and among the electrons themselves (Zeng et al. 2023; Huang et al. 2024).

A few past models share some common features with the present plasma model discussed here. Potekhin et al. (2005) (PMC05) investigated the EOS for partially ionised carbon in a wide range of densities using free energy minimisation within the framework of the chemical picture of plasmas. In their approach, each electronic configuration was treated with its own neutrality volume separately, as in the study (Fu et al. 2018). The long-range Coulomb contributions were then added separately, using previous results for the free energy of fully ionised plasmas. Starrett & Saumon (2012) proposed a fully variational AA model with the DFT method and the hypernetted chain (HNC) approximation for the ion-ion correlations. The differences between the present approach and that used in PMC05 are as follows: PMC05 considered the long-range Coulomb contributions separately, using previous results for the free energy of fully ionised plasmas, whereas in the present work, the Coulomb interactions are included within the HNC approximation. The main difference between our work and that of Starrett & Saumon (2012) is that they used a fully variational AA model with the DFT method, while we specified the ionic stages using the Saha equation.

2 Theoretical method

The free-energy minimisation method is a universal tool for determining the ionisation balance and producing the EOS for plasmas. Here, the minimisation procedure is based on the ‘chemical picture’ in which one identifies particle clusters (atoms, ions, and molecules) as distinct particle species within the plasma (Ebeling et al. 1977; Mihalas et al. 1988). If molecule formation can be ignored in a moderately high-temperature plasma, then only the electrons and ions of different charge states are treated explicitly. In this chemical picture, the free energy of a plasma volume V at temperature T consists of the following: the translational part of the ions (F₁), the part associated with the internal excitation of the species with excited states (F₂), the free energy of the ideal free electrons (F₃), and the free energy from Coulomb interactions among all charged particles (ions and electrons) in the plasma (F₄). This is given as $$F(T,V,\left\{N_i\right\})=F_1+F_2+F_3+F_4,$$(1)

where $$F_1=-kT\sum _i{N}_i(\frac{3}{2}\ln T+\ln V-\ln N_i+\frac{3}{2}\ln \frac{2\pi m_{i}k}{h^2}+1),$$(2) $$F_2=\sum _i{N}_i(E_{1,i}-kT\ln Z_i),$$(3) $$F_3=-kT{N}_e[\frac{3}{2}\frac{F_{3/2}(\eta )}{F_{1/2}(\eta )}-\eta ].$$(4)

Here, the occupation number {N_i} refers to the number of particles in charge state i contained within volume V in thermodynamic equilibrium; E_1,i and Z_i [see the detailed expression in Eq. (10)] are the ground state energy and the internal partition function of the charge state i, respectively; k is the Boltzmann’s constant; m_i is the ionic mass of the charge state i; and N_e is the population number of free electrons. Further details on F₁, F₂, and F₃ can be found in the referenced literature (Mihalas et al. 1988; Graboske et al. 1969). In the following expression, F_n(η) is the standard Fermi-Dirac integral of order n, and η is the degeneracy parameter, which is related to the number of electrons by $$F_{1/2}(\eta )=\frac{\sqrt{\pi }}{G_e}\frac{N_e}{V{T}^{3/2}},$$(5)

where $$G_e={\left(\frac{2\pi m_{e}k}{h^2}\right)}^{3/2},$$(6)

and m_e is the rest mass of an electron.

The most challenging part is the calculation of F₄, where the Coulomb interactions are influenced by complex plasma screening and other intricate physical effects. The free energy of Coulomb interactions among the charged particles (i.e. ions and free electrons) in the plasma can be determined by the radial distribution function $$\begin{array}{rl}{F}_4=& \sum _{i<j}\frac{N{N}_i}{2V}\int d\mathit{r}\frac{q_i{q}_j{e}^2}{r}{h}_{ii}(r)[g_{ij}(r)-1]+\frac{{\left(N_e\right)}^2}{2V}\int d\mathit{r}\frac{e^2}{r}\\ & \times h_{ee}(r)[g_{ee}(r)-1]-\sum _i\frac{N_e{N}_i}{2V}\int d\mathit{r}\frac{q_i{e}^2}{r}{h}_{ie}(r)[g_{ie}(r)-1],\end{array}$$(7)

where N is the total ionic number; q_i refers to the charge number of the charge state i in the unit of electron charge; and g_ij(r), g_ee(r), and g_ie(r) represent the radial distribution function between the charge state i and j, the free electrons, and the charge state i and the free electrons, respectively. The plasma screening factors h_ii(r), h_ie(r), and h_ee(r) describe the screening of Coulombic interaction between the charged particles in the plasma. In the AA model, details on obtaining these radial distribution functions can be found in the reference (Hou et al. 2015). The free energy of Coulomb interactions between the electrons, and between the electrons and ions, is incorporated through the non-ideal characteristics (NIC) correction in the free electron partition function calculation below. Therefore, these contributions should not be counted again in the present calculation.

The free energy F(T, V, {N_i}) is minimized with respect to all variations permitted by the stoichiometric relations for the ionisation equilibrium process from charge state i to a higher state i + 1 plus an electron. This yields $$\frac{\mathrm{\partial }F}{\mathrm{\partial }{N}_i}-\frac{\mathrm{\partial }F}{\mathrm{\partial }{N}_{i+1}}-\frac{\mathrm{\partial }F}{\mathrm{\partial }{N}_e}=0.$$(8)

In this context, the population number N_i+1 of charge state i + 1 in a phase space d³rd³p at a spatial position r and momentum p is determined by $$\frac{N_{i+1}^{rp}}{N_i}=\frac{Z_{i+1}}{Z_i}\frac{1+X_{i+1}}{1+X_i}\frac{1}{N_e}\frac{2d^3\mathit{r}{d}^3\mathit{p}}{h^3}{e}^{-\frac{{\varphi }_i-\mathrm{\Delta }{\varphi }_i}{kT}}\times \frac{e^{\frac{{\mu }_i}{kT}}}{e^{\frac{ϵ+[V_i(r)+V_{scr}(r)]+{\mu }_i}{kT}}+1},$$(9)

where Z_i is the partition function of charge state i, N_e is the population number of free electrons, ϕ_i is the ionisation potential (IP) of isolated charge state i, Δϕ_i is the IPD for ions belonging to the charge state i due to the plasma environment, μ_i is the chemical potential of charge state i in the plasma (determined by requiring that the integrated free electron number equals the number of ionised electrons for charge state i), and ε is the kinetic energy of the free electron. In the above expression, h is the Planck constant, and k is the Boltzmann constant.

The internal partition function Z_i of the charge state i is determined by the energies E_j,i of all possible bound quantum states in the dense plasma, relative to the ground state energy E_1,i, and their statistical weights g_j,i, as follows: $$Z_i=\sum _j{g}_{j,i}{e}^{-[E_{j,i}-E_{1,i}]/kT}.$$(10)

The term X_i in Eq. (9) originates from the inter-ionic correlation (Fu et al. 2018), given by $$X_i=\frac{1}{{\sigma }_{ij}}\sum _j\frac{N_j}{V}\int (e^{\frac{V_{ij}(r)}{kT}}-1)4\pi r^{2}dr,$$(11)

where σ_ij=2 if i = j, and σ_ij=1 if i ≠ j. V_ij(r) is the pair potential of ions i and j. The factor (1 + X_i+1)/(1 + X_i) in Eq. (9) reflects the effect of ion coupling on the Saha equation (Fu et al. 2018). The total potential includes the single-electron effective potential V_i(r) of the isolated charge state i, without the influence of the plasma environment, and the screening potential V_scr(r) arising from the plasma environment. By integrating over the phase space, we obtain the modified Saha equation, $$\frac{N_{i+1}{N}_e}{N_i}=\frac{b_i{Z}_e{Z}_{i+1}}{Z_i}\frac{1+X_{i+1}}{1+X_i}{e}^{-\frac{{\varphi }_i-\mathrm{\Delta }{\varphi }_i}{kT}},$$(12)

where Z_e is the partition function of the free electrons which are assumed to be uniform and ideal, and thus equivalent for all charge states: $$Z_e=2{\left(\frac{2\pi m_{e}kT}{h^2}\right)}^{3/2},$$(13)

In dense plasmas, the partition functions of the free electrons belonging to different charge states vary. The expression b_iZ_e represents the partition function of the free electrons associated with charge state i, where b_i is given by $$\begin{array}{rl}{b}_i=& \frac{1}{{\int }_0^{R_0}4\pi r^{2}dr}{\int }_0^{R_0}dr{\int }_0^{\mathrm{\infty }}dϵ\frac{8\pi }{h^3}{\left(2m_e\right)}^{\frac{3}{2}}{r}^2{ϵ}^{\frac{1}{2}}\\ & \times \frac{e^{\frac{{\mu }_i}{kT}}}{e^{\frac{ϵ+\mid V_i(r)+V_{scr}(r)+{\mu }_i}{kT}}+1}.\end{array}$$(14)

In the above equation, R₀ denotes the radius of a spherical volume containing the ion of charge state i, determined by the matter density as $R_0={\left(\frac{3}{4\pi n_{\text{ion}}}\right)}^{\frac{1}{3}}$, with n_ion being the ion number density. By including the NIC effects on the free electron partition functions, the standard Saha equation is generalised to a dense plasma regime. The above Saha equation can be solved self-consistently together with the constraint of particle and charge conservation: $$N=\sum _i{N}_i,$$(15) $$N_e=\sum _i{q}_i{N}_i.$$(16)

Once the free energy F and the charge state distribution (CSD) {N_i} are determined by solving the modified Saha equation, the pressure can be obtained as $$p={\frac{\mathrm{\partial }F}{\mathrm{\partial }V}|}_{T,\left\{N_i\right\}}$$(17)

by direct differentiation of the free energy F. Because the above equations satisfy the Maxwell relations, thermodynamic consistency of these variables is guaranteed (Mihalas et al. 1988). Since both the modified Saha equation and the pressure are derived from the minimisation of the free energy, the present theoretical formalism can, in principle, be applied to plasmas up at arbitrary high densities, provided that the NIC affecting the IPDs, the free energy, and other physical quantities related to the ionisation balance equation are properly taken into account. However, considering all these NIC effects is a challenging task.

Most previous studies omitted certain physical effects, such as treating free electrons in plasmas as an ideal gas. This made it difficult to apply their theoretical methods to extremely dense plasmas. However, some studies, such as Potekhin et al. (2005), treated free electrons consistently with the effective one-electron potential of the ions.

All the relevant data required to solve the modified Saha equation are obtained by solving the Dirac equation for each charge state (atomic units are used here unless otherwise specified), as follows: $$\begin{array}{r}[\sum _{i=1}^{M}c\alpha •p_i+(\beta -1)c^2+V_i^{nuc}(r)+\sum _{i<j}^M\frac{1}{r_{ij}}+V_{scr}(r)]\psi (\mathit{r})\\ =E\psi (\mathit{r}),\end{array}$$(18)

where M denotes the number of bound electrons in the respective ion, α and β are the Pauli matrices, and V_scr(r) is the plasma screening potential. The rest energy c² is subtracted out in the above equation. The plasma screening potential V_scr(r) is determined by the radial distribution function of the charge density of the free electrons ρ(r) (Son et al. 2014; Zeng et al. 2020a,b) $$V_{scr}(r)=4\pi [\frac{1}{r}{\int }_0^r{r}_1+{\int }_r^{R_0}]r_1\rho \left(r_1\right)\mathrm{d}{r}_1-\frac{3}{2}{[\frac{3}{\pi }\rho (r)]}^{1/3}.$$(19)

The free-electron charge density distribution function follows the Fermi-Dirac distribution, $$\rho (r)=\frac{1}{{\pi }^2}{\int }_{k_0(r)}^{\mathrm{\infty }}\frac{k^{2}dk}{e^{\frac{\sqrt{k^2{c}^2+c^4}-c^2+V_{scr}(r)+{\mu }_i}{kT}}+1},$$(20)

where c denotes the speed of light in vacuum, p is the free-electron momentum, and $p_0(r)=\frac{\sqrt{-2V(r)c^2+V(r{)}^2}}{c}$. If V_scr(r) is omitted in the above Dirac equation, then it corresponds to the atomic structure and properties of the isolated ion. The IP of charge state i is determined by the energy difference between the respective ground levels of charge states i and i + 1. By solving the Dirac equation of charge state i of both the isolated ion (without including V_scr(r)) and the screened ion embedded in the plasma (including V_scr(r)), the IPD, Δϕ_i, of the charge state i – required for solving the Saha equation – can be obtained as the difference between the IPs of the isolated ion and the screened ion. Further details on solving the Dirac equation can be found in the flexible atomic code (Gu 2008).

The basis states of the atomic wave functions are constructed from the products of N one-electron Dirac spinors, $${\phi }_{n\kappa m}=\frac{1}{r}(\frac{P_{n\kappa }(r){\chi }_{\kappa m}(\theta ,\psi ,\sigma )}{i{Q}_{n\kappa }(r){\chi }_{-\kappa m}(\theta ,\psi ,\sigma )}),$$(21)

where P_nk(r) and Q_nk(r) are the large and small components of the radial wave functions, χ_κm(θ, ψ, σ) is a two-component spherical spinor, and n, κ, and m denote the principal, the relativistic angular, and the magnetic quantum numbers of the electron orbital, respectively. The large and small components of the radial wave function satisfy the coupled equations, $$\begin{array}{rl}& (\frac{d}{dr}+\frac{\kappa }{r})P_{n\kappa }=\alpha ({\epsilon }_{n\kappa }-V(r)+\frac{2}{{\alpha }^2})Q_{n\kappa }\\ & (\frac{d}{dr}-\frac{\kappa }{r})Q_{n\kappa }=\alpha (-{\epsilon }_{n\kappa }+V(r))P_{n\kappa },\end{array}$$(22)

where α denotes the fine-structure constant, V(r) represents the self-consistent potential including contributions from both the bound electrons and the plasma screening potential V_scr(r) from the free electrons, and ε_nκ is the energy eigenvalue of the electron orbital relative to the potential at the ion sphere boundary.

The radial wave function is solved inside the ion sphere with the following boundary conditions for the radial wave function of a bound orbital (Hou et al. 2006): $$P_{nk}(0)=Q_{nk}(0)=0$$(23)

and $$P_{n\kappa }\left(R_0\right)=0.$$(24)

In the above theoretical formalism, we used a non-relativistic expression of the thermodynamic free energy in practical calculations. However, the relativistic Dirac equation must be used to describe the atomic structure and properties in order to obtain the physical quantities required for calculating the free energy, such as the internal partition functions. This is because the relativistic effects cannot be ignored for the inner-shell orbitals of medium and heavy elements. The use of the relativistic FermiDirac distribution is justified for the same reason as the use of the Dirac equation. The local kinetic energy of a free electron close to the nucleus can reach tens of thousand eV. Therefore, relativistic effects should be included.

3 Results and discussion

In Fig. 1, we show the prediction of IPDs (black lines, left y axis) and IPs (blue lines, right y axis) for the charge states of C²⁺ through C⁵⁺ as a function of mass density in carbon plasmas at temperatures of 50, 200, and 500 eV. Results are presented for both our theoretical model and for the analytical models of the Stewart–Pyatt (SP) (Stewart & Pyatt 1966) and the EckerKröll (EK) (Ecker & Kröll 1963) at a temperature of 200 eV. As shown, all the theories predict an evident lowering of IPs – even at a lower density for all charge states. The IPDs essentially show an increasing effect with increasing mass density. However, they are also temperature dependent, and their values slightly decrease as the temperature increases. Among the theories, the EK model predicts the highest IPD, while the SP model predicts the lowest IPD. Our results lie between the two, demonstrating a better agreement with experimental values, as also reported in previous publications (Zeng et al. 2020a,b; Huang et al. 2024).

As shown in Fig. 1, the 2s electron becomes delocalised in C²⁺ at a mass density of 1.55 and in C³⁺ at 1.68 g cm⁻³, both at a temperature of 50 eV. With increasing temperature, the density at which the transition from a bound electron orbital to delocalisation occurs increases slightly, reaching 1.73 and 1.85 g cm⁻³, at a temperature of 200 eV, respectively. However, the K-shell delocalisation of C⁴⁺ and C⁵⁺ occurs at much higher mass densities of 228 and 256 g cm⁻³, respectively, at a temperature of 200 eV. We later demonstrate that the K-shell delocalisation is crucial for determining the CSD and the EOS.

Accurate prediction of IPDs of ions in dense plasmas is crucial for determining the ionisation balance. In particular, the precise determination of the K-shell delocalisation depends sensitively on the IPDs. The advancement of experimental technology makes it now feasible to quantitatively measure the IPDs of ions at a solid density (Ciricosta et al. 2012, 2016), or even a few times the solid density (Hoarty et al. 2013; Fletcher et al. 2014). An experimental platform was developed to measure the IPD – one of the most essential factors in determining the ionisation degree of plasmas – using spectrally resolved X-ray scattering (Kraus et al. 2016, 2019). Our previous predictions on IPDs based on our model, as mentioned above, agree with these experimental measurements (Zeng et al. 2020a; Huang et al. 2024) and, hence, supports the reliability of our theoretical predictions for dense carbon plasmas as well.

Another factor that influences the CSD and ionisation balance is the NIC effect on the partition functions of the free electrons in dense plasmas. The relative partition functions of the free electrons are shown in Fig. 2 for the charge states C⁴⁺ and C⁵⁺ in dense carbon plasmas at mass densities of 3, 10, 100, and 200 g cm⁻³ at a temperature of 200 eV. As shown, the partition functions of free electrons of different charge states of carbon ions are reduced compared to those in ideal plasmas.

Using the obtained IPDs and the partition functions of free electrons, the Saha ionisation balance equation can be solved with the NIC effect taken into account. As examples, Fig. 3 shows the population fractions of different charge states of carbon ions at mass densities of 5, 100, and 235 g cm⁻³ and at temperatures of 50, 100, 200, and 500 eV. For comparison, the CSDs calculated without including the NIC effect on the partition functions of free electrons are also shown, in order to evaluate this non-ideal effect on the CSDs. At the lowest indicated mass density of 5 g cm⁻³, nearly all the ions are C⁴⁺ at temperatures below 50 eV. In this case, the NIC effect on the CSDs, arising from the partition functions of free electrons, is greatly restrained, resulting in nearly identical values with and without the inclusion of the NIC effect. As the temperature increases to 100 and 200 eV, this non-ideal effect becomes stronger than at 50 eV. However, as the temperature increases further to 500 eV, this non-ideal effect becomes much weaker again, and the most abundant charge state is the bare ion, C⁶⁺. At a mass density of 100 g cm⁻³, similar to 5 g cm⁻³, the NIC effect on the CSDs resulting from the partition functions of free electrons is also negligible at temperatures below 50 eV. However, this non-ideal effect exerts a more noticeable influence on the CSDs at temperatures above 50 eV. At a mass density of 235 g cm⁻³, the charge state of C⁴⁺ disappears as a result of K-shell delocalisation. In this case, the non-ideal effect exerts noticeable influence on the CSDs, even at a higher temperature of 500 eV.

Figure 4 shows the average degree of ionisation of dense carbon plasmas as a function of mass density up to 300 g/cm³ at temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV. Results for the average degree of ionisation are provided with and without considering the NIC effect on the partition functions of free electrons in the plasmas. The dynamical effect of the dense plasma environment on the ionisation balance was taken into account when calculating the average degree of ionisation. This was obtained by accounting for the dynamical broadening effect of the plasma on the energy levels and ionisation threshold of all the relevant charge states (Hou et al. 2006). We find that the NIC effect on the average degree of ionisation is negligible at and below 50 eV, consistent with the CSD predictions shown in Fig. 3. As the plasma temperature increases to 300 eV, this effect becomes increasingly significant. However, at temperatures above 800 eV, the NIC effect gradually weakens again. Studies to date demonstrate that the NIC effect influences the ionisation balance in dense plasmas, and that including this effect leads to improved agreement with experimental measurements (Zeng et al. 2023; Huang et al. 2024).

For all the plasma temperatures considered, as shown in Fig. 4, a pressure-induced ionisation process is predicted for the K-shell electrons of the C⁴⁺ and C⁵⁺ charge states. Moreover, the mass density at which K-shell delocalisation occurs also depends on temperature, generally shifting to higher densities as the plasma temperature increases. Below a mass density of ~180 g cm⁻³, thermal ionisation dominates the ionisation balance of dense plasma across all studied temperatures. However, the average degree of ionisation rises sharply to 5.0 at densities above ~217 g cm⁻³, and further increases to 6.0 at densities above ~246 g cm⁻³. The first steep increase in the average degree of ionisation corresponds to the K-shell delocalisation of C⁴⁺, making C⁵⁺ the dominant ion type. The second steep increase in the average degree of ionisation corresponds to the K-shell delocalisation of C⁵⁺, leading to a plasma consisting entirely of bare ions and an average degree of ionisation of 6.0. Above the K-shell delocalisation density threshold, pressure-induced ionisation plays a more important role than thermal ionisation, particularly at lower temperatures. Below the pressure-induced ionisation threshold, the ionisation balance displays typical features of thermal ionisation: the average degree of ionisation increases with rising plasma temperature and decreases with increasing mass density.

Figure 4 illustrates the interplay between pressure-induced K-shell ionisation and thermal ionisation, with the competition between these processes becoming increasingly evident at higher plasma temperatures. At the lowest temperature considered here (i.e. 50 eV), the average degree of ionisation remains nearly constant at 4.0, reflecting the delocalisation of the L-shell electrons from C³⁺ ions. Here, pressure-induced L-shell ionisation dominates the plasma balance. Thermal ionisation can only be observed at the lowest mass densities, below ~2 g cm⁻³. Additionally, an anomalous phenomenon is predicted: the average degree of ionisation increases slightly as the mass density rises up to ~180 g cm⁻³. Such a phenomenon is attributed to the increasing IPDs of carbon ions with rising mass density. At higher temperatures, this anomaly disappears, and the average degree of ionisation decreases with increasing mass density. Below the pressure-induced K-shell ionisation threshold, thermal ionisation plays an increasingly important role in the ionisation balance as mass density increases. Accurately determinating the pressure-induced K-shell ionisation threshold is therefore vital for describing the competition between pressure-induced K-shell ionisation and thermal ionisation.

In Fig. 5, we compare our theoretical results for the average degree of ionisation of dense carbon plasma at a temperature of 100 eV as a function of mass density with other widely used ionisation models in astrophysics and inertial confinement fusion. Among all the theories, the calculation carried out by Potekhin, Massacrier, and Chabrier (i.e. PMC05, performed at a slightly lower temperature of 10⁶ K compared to the 100 eV considered here) (Potekhin et al. 2005) predicts a mass density threshold for pressure-induced K-shell ionisation (~200 g cm⁻³) that is closest to our theoretical result. Above a density of 200 g cm⁻³, our predictions are in reasonable agreement with that by PMC05. However, below this density threshold, there is significant discrepancy: PMC05 predicts that the average degrees of ionisation increases with rising mass density, exhibiting three local minima at mass densities of ~1.1, ~2.2, and ~190 g cm⁻³. In contrast, our prediction yields a decreasing average degree of ionisation as mass density increases. We still do not understand why the PMC05 results exhibit the three local minima. A similar difference is also observed for the DFT-MD results (Bethkenhagen et al. 2020), which closely follow the referred PMC05 (Potekhin et al. 2005) results. There is overall reasonable consistency between our results and those of the neutral pseudoatom (NPA) model (Dharma-wardana 2021), which predicts a slower rise in the average degree of ionisation at densities above ~100 g cm⁻³. At ~300 g cm⁻³, the ATOMIC model (Hakel et al. 2006) predicts a higher threshold of pressure-induced K-shell ionisation than ours. The Beth-Uhlenbeck (BU) calculations – obtained via Eq. (15) in the reference (Röpke et al. 2019) – which take into account the PB effects (short-dashed black line) predict a much lower density threshold of pressure-induced K-shell ionisation at ~20 g cm⁻³.

Other theories either predict a much lower density threshold of pressure-induced K-shell delocalisation, accompanied by much smoother and softer increase in the average degree of ionisation with increasing mass density (e.g. DFT-MD (Bethkenhagen et al. 2020), PURGATORIO (Sterne et al. 2007), the quantum AA model (QAAM) (Faussurier et al. 2021), and the NPA model (Dharma-wardana 2021)), or do not predict any pressure-induced K-shell delocalisation at all (e.g. OPAL (Rogers et al. 1996) and both BU calculations with the SP and EK models (Röpke et al. 2019)). At lower densities, below the threshold for pressure-induced K-shell delocalisation considered here, the predictions for the average degree of ionisation from three different BU models (Röpke et al. 2019), as well as from other such models as PURGATORIO (Sterne et al. 2007), OPAL (Rogers et al. 1996), and ATOMIC (Hakel et al. 2006), qualitatively agree with our results and predict a decreasing average degree of ionisation with increasing mass density. The Thomas-Fermi (TF) model with a Dirac exchange correction (TFD), along with definitions of ionisation based on its effect on plasma structure (Qs) and on total pressure (Qp) (Clerouin et al. 2022), yields an average degree of ionisation below 4.0 near a mass density of 15 g cm⁻³. This indicates that the L-shell electrons remain bound to C³⁺ ions or ions with even lower charge states.

As is well known, EOS is closely related to CSDs and the average degree of ionisation. In Fig. 6, we show pressure as a function of mass density in dense carbon plasmas up to 400 g cm⁻³ at temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV with and without (b_i=1) considering the NIC effect on the partition functions of free electrons in the plasmas. A general increase in plasma pressure with rising mass density and temperature is observed across the plasma conditions studied. Notably, a steeper increase in plasma pressure with rising mass density occurs at ~300 g cm⁻³ for all the temperatures considered here.

This reflects the pressure-induced K-shell delocalisation discussed above. The NIC effect on the partition functions of free electrons exerts less influence on the pressure at lower temperatures (50 and 100 eV) and at the higher temperature of 800 eV, compared to other intermediate temperatures. Furthermore, the NIC effect on pressure becomes negligible at densities above ~300 g cm⁻³, where the K-shell delocalisation occurs. At a temperature of 300 eV, the pressure can be overestimated by ~10% at densities below ~260 g cm⁻³ if the NIC effect on the partition functions of free electrons in the plasmas is not included in the calculation.

Figure 7 compares the pressure as a function of mass density predicted by our theory with results from the first principles molecular dynamics (FPMD) method (Clerouin et al. 2022), the QAAM, and the quantum molecular dynamic (QMD) simulations (Faussurier et al. 2021) at a temperature of 100 eV. There is reasonable agreement between our calculated pressure and that from FPMD, QAAM, and QMD below a mass density of 220 g cm⁻³. Above 220 g cm⁻³, however, our model predicts a higher pressure than the above methods due to the onset of K-shell delocalisation in C⁴⁺ and C⁵⁺. For example, at densities above 220 g cm⁻³, our predicted pressure is 24% higher than that from QMD. Such a significant difference in the EOS of carbon plasma can impact the evolution of white dwarfs, most of which have carbon and oxygen cores, and hot DQ white dwarfs with carbon-dominated atmospheres.

Using the present theoretical formalism, we can better understand the discrepancy between the theoretical and experimental data on the opacity of the dense hydrocarbon (C₉H₁₀) (Kritcher et al. 2020). In the experiment, the C₉H₁₀ samples were mixed with uniform trace amounts of fluorine (1 at.%). In Fig. 8, we show the opacity of dense C₉H₁₀ with 1 at.% fluorine obtained using detailed level accounting (DLA) method, normalised to that of the cold material, at a photon energy of 9 keV and compare it with the experimental results and theoretical data of non-equilibrium modelling using the detailed configuration accounting (DCA) method (Scott & Hansen 2010). Our results show better agreement with the data obtained from experimental measurements (Kritcher et al. 2020) when the NIC effect is properly included in the calculation of ionisation balance.

Fig. 1 Ionisation potential depressions and IPs as a function of mass density in dense carbon plasmas at plasma temperatures of 50, 200, and 500 eV for the charge states C2+ (a), C3+ (b), C4+ (c), and C5+ (d). For each plot, the left y axis represents the IPDs (black lines), and the right y axis represents the IPs (blue lines). The meaning of the black lines in plot (a) for the IPDs and the blue lines in plot (b) for the IPs applies to all plots (a)–(d) for the charge states C2+, C3+, C4+, and C5+. For comparison, predictions obtained from the SP (Stewart & Pyatt 1966) and EK (Ecker & Kröll 1963) analytical models at a temperature of 200 eV are shown. The legends in plot (a) for the IPDs (black lines) and in plot (b) for the IPs (blue lines) apply for all plots for the charge states of C2+-C5+.

Fig. 2 Relative partition functions (bi in the modified Saha equation) of the free electrons for different charge states in dense carbon plasmas at a temperature of 200 eV and the mass densities of 3, 5, 100, and 200 g cm−3. For each plasma condition, the partition functions are scaled relative to that of the uniform electron gas $(Z_e=2{\left(\frac{2\pi m_{e}kT}{h^2}\right)}^{3/2})$.

Fig. 3 Population fractions of different charge states of dense carbon plasmas at mass densities of 5 (a), 100 (b), and 235 g cm−3 (c), shown for temperatures of 50, 100, 200, and 500 eV. For each plasma density and temperature, solid and dashed lines indicate the results with and without the NIC effect on the partition functions of free electrons, respectively. All other effects, including the IPDs and internal partition functions arising from the plasma environment, are treated identically in both sets of results.

Fig. 4 Average degree of ionisation of carbon plasmas as a function of mass density at different plasma temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV. Black lines and red lines with circles refer to the average degree of ionisation with and without the NIC effect on the partition functions of free electrons in the plasmas, respectively.

Fig. 5 Comparison of the average degree of ionisation of carbon plasmas as a function of mass density at a plasma temperature of 100 eV, predicted by various theoretical methods. Shown are results from DFTMD simulation (black dotted line) (Bethkenhagen et al. 2020); BU calculations (Röpke et al. 2019) incorporating two different IPD models – SP (long black-dashed line) and EK (dotted, short-dashed black line) – and including the PB effects (short-dashed black line); PURGATORIO (dotted, long-dashed black line) (Sterne et al. 2007); OPAL (dot-dot-dashed black line) (Rogers et al. 1996); ATOMIC (dot-dashed-dashed black line) (Hakel et al. 2006); TF and TFD models (Clerouin et al. 2022); Qp and Qs ionisation (Clerouin et al. 2022); QAAM (Faussurier et al. 2021); PMC05 calculation at a temperature of 106 K (Potekhin et al. 2005); and the NPA model (Dharma-wardana 2021).

Fig. 6 Pressure as a function of mass density in carbon plasmas at different temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV with and without considering the NIC effect on the partition functions of free electrons in the plasmas.

Fig. 7 Comparison of pressure as a function of mass density for carbon plasmas at a temperature of 100 eV, as predicted by the FPMD method (Clerouin et al. 2022), the QAAM model, and the QMD simulations (Faussurier et al. 2021).

Fig. 8 Opacity of dense C9H10 with 1 at.% fluorine (as a tracer element), normalised to that of cold material, at a photon energy of 9 keV as a function of pressure. The NIC effect on the partition functions of free electrons is included in our DLA method. Plasma conditions are determined at the shock front along the shock Hugoniot measured by Kritcher et al. (2020). Experimental and theoretical data from DCA (Kritcher et al. 2020) are shown for comparison. The upper and lower limits of the experimental data indicate the uncertainty contours of ±1σ in the measurement.

4 Summary and conclusions

Based on the chemical picture of free energy minimisation, the Saha equation approach was modified and extended to apply to an extremely high density plasma regime with pressures of a few 100 Gbar. This allowed for the investigation of the ionisation balance and EOS by properly accounting for NIC effects on free energy, particularly those arising from Coulomb interactions between electrons and ions, and among electrons themselves. The modified equation was then applied to study the ionisation balance and EOS of hot, dense carbon plasmas.

At a temperature of 200 eV, the 2 s electron in C²⁺ and C³⁺ is predicted to become delocalised at mass densities of 1.73 and 1.85 g cm⁻³, respectively, whereas the K-shell delocalisation for C⁴⁺ and C⁵⁺ occurs at much higher mass densities of 228 and 256 g cm⁻³. Accurately predicting the critical density for pressure-driven K-shell delocalisation is crucial for determining the average degree of ionisation and EOS. At temperatures below 50 eV, pressure-driven L-shell delocalisation dominates the plasma equilibrium below the density of 220 g cm⁻³, resulting in an almost constant average degree of ionisation of 4.0 Above that density, pressure-induced K-shell ionisation becomes the dominant process affecting the plasma ionisation balance.

As the plasma temperature increases, a competition arises between pressure-induced K-shell ionisation and thermal ionisation. Which mechanism dominates depends on plasma density and temperature. Below the density threshold for pressure-driven K-shell delocalisation, our predicted average degrees of ionisation for dense carbon plasmas show reasonable agreement with other theoretical models. However, the prediction results show large discrepancies in the critical density at which K-shell electrons become delocalised by pressure. Above ~220 g cm⁻³, our model generally predicts higher pressure than most other theories, primarily due to the onset of K-shell delocalisation.

Acknowledgements

This work was supported by the National Natural Science Foundation of China under Grants No. 12174343, No. 12274384, and No. 12335015.

References

All Figures

Fig. 1 Ionisation potential depressions and IPs as a function of mass density in dense carbon plasmas at plasma temperatures of 50, 200, and 500 eV for the charge states C2+ (a), C3+ (b), C4+ (c), and C5+ (d). For each plot, the left y axis represents the IPDs (black lines), and the right y axis represents the IPs (blue lines). The meaning of the black lines in plot (a) for the IPDs and the blue lines in plot (b) for the IPs applies to all plots (a)–(d) for the charge states C2+, C3+, C4+, and C5+. For comparison, predictions obtained from the SP (Stewart & Pyatt 1966) and EK (Ecker & Kröll 1963) analytical models at a temperature of 200 eV are shown. The legends in plot (a) for the IPDs (black lines) and in plot (b) for the IPs (blue lines) apply for all plots for the charge states of C2+-C5+.

In the text

	Fig. 2 Relative partition functions (bi in the modified Saha equation) of the free electrons for different charge states in dense carbon plasmas at a temperature of 200 eV and the mass densities of 3, 5, 100, and 200 g cm−3. For each plasma condition, the partition functions are scaled relative to that of the uniform electron gas $(Z_e=2{\left(\frac{2\pi m_{e}kT}{h^2}\right)}^{3/2})$.
In the text

Fig. 3 Population fractions of different charge states of dense carbon plasmas at mass densities of 5 (a), 100 (b), and 235 g cm−3 (c), shown for temperatures of 50, 100, 200, and 500 eV. For each plasma density and temperature, solid and dashed lines indicate the results with and without the NIC effect on the partition functions of free electrons, respectively. All other effects, including the IPDs and internal partition functions arising from the plasma environment, are treated identically in both sets of results.

In the text

	Fig. 4 Average degree of ionisation of carbon plasmas as a function of mass density at different plasma temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV. Black lines and red lines with circles refer to the average degree of ionisation with and without the NIC effect on the partition functions of free electrons in the plasmas, respectively.
In the text

Fig. 5 Comparison of the average degree of ionisation of carbon plasmas as a function of mass density at a plasma temperature of 100 eV, predicted by various theoretical methods. Shown are results from DFTMD simulation (black dotted line) (Bethkenhagen et al. 2020); BU calculations (Röpke et al. 2019) incorporating two different IPD models – SP (long black-dashed line) and EK (dotted, short-dashed black line) – and including the PB effects (short-dashed black line); PURGATORIO (dotted, long-dashed black line) (Sterne et al. 2007); OPAL (dot-dot-dashed black line) (Rogers et al. 1996); ATOMIC (dot-dashed-dashed black line) (Hakel et al. 2006); TF and TFD models (Clerouin et al. 2022); Qp and Qs ionisation (Clerouin et al. 2022); QAAM (Faussurier et al. 2021); PMC05 calculation at a temperature of 106 K (Potekhin et al. 2005); and the NPA model (Dharma-wardana 2021).

In the text

	Fig. 6 Pressure as a function of mass density in carbon plasmas at different temperatures of 50, 100, 150, 200, 250, 300, 500, and 800 eV with and without considering the NIC effect on the partition functions of free electrons in the plasmas.
In the text

	Fig. 7 Comparison of pressure as a function of mass density for carbon plasmas at a temperature of 100 eV, as predicted by the FPMD method (Clerouin et al. 2022), the QAAM model, and the QMD simulations (Faussurier et al. 2021).
In the text

Fig. 8 Opacity of dense C9H10 with 1 at.% fluorine (as a tracer element), normalised to that of cold material, at a photon energy of 9 keV as a function of pressure. The NIC effect on the partition functions of free electrons is included in our DLA method. Plasma conditions are determined at the shock front along the shock Hugoniot measured by Kritcher et al. (2020). Experimental and theoretical data from DCA (Kritcher et al. 2020) are shown for comparison. The upper and lower limits of the experimental data indicate the uncertainty contours of ±1σ in the measurement.

In the text