# A New Generation of Cool White Dwarf Atmosphere Models. I. Theoretical Framework and Applications to DZ Stars

S. Blouin, Nicole F. Allard, P. Dufour

► To cite this version:

S. Blouin, Nicole F. Allard, P. Dufour. A New Generation of Cool White Dwarf Atmosphere Models. I. Theoretical Framework and Applications to DZ Stars. The Astrophysical Journal, 2018, 863 (2), pp.184. 10.3847/1538-4357/aad4a9. hal-02280924

HAL Id: hal-02280924

https://hal.science/hal-02280924v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A New Generation of Cool White Dwarf Atmosphere Models. I. Theoretical Framework and Applications to DZ Stars

S. Blouin¹, P. Dufour¹, and N. F. Allard²,³

¹Département de Physique, Université de Montréal, Montréal, QC H3C 3J7, Canada; sblouin@astro.umontreal.ca, dufourpa@astro.umontreal.ca

²GEPI, Observatoire de Paris, Université PSL, CNRS, UMR 8111, 61 avenue de l'Observatoire, F-75014 Paris, France

³Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98bis boulevard Arago, F-75014 Paris, France

Received 2018 June 14; revised 2018 July 13; accepted 2018 July 17; published 2018 August 23

## Abstract

The photospheres of the coolest helium-atmosphere white dwarfs are characterized by fluidlike densities. Under those conditions, standard approximations used in model atmosphere codes are no longer appropriate. Unfortunately, the majority of cool He-rich white dwarfs show no spectral features, giving us no opportunities to put more elaborate models to the test. In the few cases where spectral features are observed (such as in cool DQ or DZ stars), current models completely fail to reproduce the spectroscopic data, signaling shortcomings in our theoretical framework. In order to fully trust parameters derived solely from the energy distribution, it is thus important to at least succeed in reproducing the spectra of the few coolest stars exhibiting spectral features, especially since such stars possess even less extreme physical conditions due to the presence of heavy elements. In this paper, we revise every building block of our model atmosphere code in order to eliminate low-density approximations. Our updated white dwarf atmosphere code incorporates state-of-the-art constitutive physics suitable for the conditions found in cool helium-rich stars (DC and DZ white dwarfs). This includes new high-density metal-line profiles, nonideal continuum opacities, an accurate equation of state, and a detailed description of the ionization equilibrium. In particular, we present new ab initio calculations to assess the ionization equilibrium of heavy elements (C, Ca, Fe, Mg, and Na) in a dense helium medium and show how our improved models allow us to achieve better spectral fits for two cool DZ stars, Ross 640 and LP 658-2.

**Key words:** equation of state – opacity – stars: atmospheres – stars: individual (LP 658-2) – stars: individual (Ross 640) – white dwarfs

## 1. Introduction

Pure helium-rich white dwarfs do not show any spectral lines when $T_{\text{eff}} \lesssim 10,000$ K. The same occurs for $T_{\text{eff}} \lesssim 5000$ K in the case of pure hydrogen-rich atmospheres. Together, these featureless white dwarfs are known as DC stars. One is thus forced to rely solely on the shape of the spectral energy distribution to deduce the chemical composition and effective temperature of these white dwarfs (Bergeron et al. 1997, 2001). Although most cool white dwarfs have featureless spectra, some cool helium-rich white dwarfs do show significant spectral features that can be exploited to retrieve additional information on the physical conditions encountered in their atmospheres. Some contain enough hydrogen to show strong H$_2$–He collision-induced absorption (CIA) features, some show C$_2$ Swan bands (DQ and DQpec stars), and others show metal lines (DZ stars). Interestingly, in all cases, models fail to reproduce these spectra. For instance, the CIA is inadequately modeled (e.g., LHS 3250, SDSS J123812.85+350249.1, SDSS J125106.11+440303.0; Gianninas et al. 2015), the C$_2$ bands are distorted (e.g., LHS 290; Kowalski 2010a), and the metal absorption lines often do not have the right strength or shape (e.g., WD 2356–209, Bergeron et al. 2005; Homeier et al. 2005, 2007; LP 658-2, Wolff et al. 2002; Dufour et al. 2007).

For all of these stars, the discrepancies between models and observations can be related to nonideal high-density effects arising at the photosphere, since for cool ($T_{\text{eff}} < 6000$ K) helium-rich white dwarfs, densities reach fluidlike values. At a Rosseland optical depth $\tau_R = 2/3$, density can be as high as $1 \text{ g cm}^{-3}$ (Bergeron et al. 1995; Kowalski 2010b), which corresponds to a fluid where the separation between atoms is

roughly equivalent to the dimension of the atoms themselves. Clearly, under such conditions, interactions between species are no longer negligible, and the ideal gas approximation must be discarded.

The nonideal effects arising from this high density have remained mostly unnoticed for DC stars, since a featureless spectrum provides little opportunity to test the accuracy of atmosphere models. In contrast, cool helium-rich stars with spectral features (i.e., DQpec, DZ, and those with CIA features) provide a real challenge to atmosphere models and an opportunity to test our understanding of the chemistry and physics of warm dense helium.

In this series of papers, we present and apply our new generation of atmosphere models for cool white dwarf stars. In the first paper of the series, we focus on improving our modeling of cool DZ stars. Note that obtaining better fits of these objects is far more than a mere aesthetic whim. Indeed, because they show spectral lines, cool DZ stars represent a unique opportunity to probe the physics and chemistry of cool helium-rich atmospheres. In a way, they allow us to test the models used for DC stars. Once we have proven that our new models are able to reproduce the rich and complex spectra of cool DZ stars, we will be confident that the constitutive physics is accurate and that the models can reliably be used to measure the atmospheric parameters of all DC stars in general.

This paper describes our new model atmosphere code that includes all nonideal effects relevant for the modeling of the atmospheres of cool DZ and DC stars. This updated atmosphere code is based on the one described in Dufour et al. (2007). Building on other published works, as well as on our own new calculations, we have considerably improved the
<table><caption>**Table 1**
Metal-line Profiles Computed Using the Unified Line Shape Theory Described in Allard et al. (1999)</caption><thead><tr><th>Lines</th><th>Source</th></tr></thead><tbody><tr><td>Ca I 4226 Å</td><td>N. F. Allard (2018, private communication)</td></tr><tr><td>Ca II H & K</td><td>Allard & Alekseev (2014)</td></tr><tr><td>Mg I 2852 Å</td><td>N. F. Allard et al. (2018, in preparation)</td></tr><tr><td>Mg II 2795/2802 Å</td><td>Allard et al. (2016a)</td></tr><tr><td>Mgb triplet</td><td>Allard et al. (2016b)</td></tr><tr><td>Na I D doublet</td><td>Allard et al. (2014)</td></tr></tbody></table>

constitutive physics in our code. Section 2 describes the additions made to correctly calculate radiative opacities, and in Section 3, we discuss the improvements related to the equation of state and the chemical equilibrium. Among the new physics added to the chemical equilibrium calculations, we used ab initio techniques to implement a state-of-the-art description of the chemical equilibrium of heavy elements (C, Ca, Fe, Mg, and Na) in the dense atmosphere of cool DZ stars. These calculations are detailed at length in Section 4. In Section 5, we present two applications that show how the improvements included in our models translate in terms of spectroscopic fits. Finally, in Section 6, we summarize our results and outline the upcoming papers of this series.

## 2. Radiative Opacities

In this section, we describe the additions brought to the code of Dufour et al. (2007) regarding the calculation of radiative opacities. This includes improved line profiles, high-density CIA distortion, and continuum opacities corrected for collective interactions.

### 2.1. Line Profiles

In the atmosphere of cool DZ stars, the wings of heavy-element absorption lines are severely broadened by interactions with neutral helium. Hence, Lorentzian profiles poorly reproduce observed spectral features. It is thus an absolute necessity to implement the unified line shape theory described in Allard et al. (1999) to treat such line profiles. We implemented this formalism for the strongest transitions found in cool DZ white dwarfs (see Table 1). In particular, the line profiles described in Allard & Alekseev (2014), Allard et al. (2014, 2016a, 2016b), and N. F. Allard et al. (2018, in preparation) are used to compute the wings, and a conventional Lorentzian profile is assumed for the core of the spectral lines, where the density is low enough for this approximation to hold. To connect the two profiles, we use a hyperbolic tangent function, which allows a smooth transition. It should also be noted that our Ca I 4226 Å profile is still preliminary, as we do not yet have access to the high-quality ab initio potentials required for the computation of this particular line profile. To make up for this lack, we computed our own ab initio potentials through open-shell configuration-interaction singles calculations with the ROCIS module of the ORCA quantum chemistry package<fnref n="4" /> (Neese 2012).

For transitions not listed in Table 1, our code assumes a simple Lorentzian function or quasistatic van der Waals broadening (Walkup et al. 1984; D. Koester 2018, private

communication). Note that the exact treatment of these secondary transitions has a limited impact on our atmospheric determinations.

We show in Figure 1 a comparison of line profiles calculated using the theory of Allard et al. (1999) to those found assuming a Lorentzian profile for temperature and density conditions representative of the photosphere of cool DZ stars. Clearly, under such conditions, the Lorentzian function fails to provide a satisfactory description of the line profiles. It underestimates the strong broadening observed in the more accurate line profiles and does not take into account the distortion and shift observed for many transitions.

### 2.2. Collision-induced Absorption

The calculation of the H₂–He CIA includes the high-density distortion effects described in Blouin et al. (2017). This pressure distortion effect alters the infrared energy distribution of cool DZ stars with hydrogen in their atmospheres and a photospheric density greater than ≈0.1 g cm⁻³ (n<sub>He</sub> = 1.5 × 10²² cm⁻³). Moreover, we have also included the He–He–He CIA using the analytical fits given in Kowalski (2014).

### 2.3. Rayleigh Scattering

In a dense helium medium, collective interactions between atoms lead to a reduction of the Rayleigh-scattering cross section (Iglesias et al. 2002). For the wavelength domain relevant for white dwarf modeling (i.e., in the low-frequency limit), the reduced cross section can be expressed as (Kowalski 2006a; Rohrmann 2018)

$$ \sigma_{\text{Rayleigh}}(\omega) = S(0)\sigma_{\text{Rayleigh}}^{0}(\omega), \quad (1) $$

where σ₀<sub>Rayleigh</sub>(ω) is the ideal gas result (e.g., Dalgarno 1962) and S(0) is the structure factor of the fluid at a wavenumber k = 0. Therefore, to take into account the reduction of the Rayleigh scattering, we simply need to know S(0), which is a function of the temperature and density of the helium fluid. To compute S(0), we use the analytical fit to the Monte Carlo results of Rohrmann (2018).

### 2.4. He⁻ Free-Free Absorption

Iglesias et al. (2002) also showed that the free-free absorption cross section of the negative helium ion is reduced in a dense helium medium. Given that it is the dominant source of opacity in DZ stars, it is important to take this reduction into account. The corrected cross section for He⁻ free-free absorption is given by Iglesias et al. (2002),

$$ \sigma_{\text{ff}}(\omega) = \delta_{\text{ff}}(\omega)\sigma_{\text{ff}}^{0}(\omega), \quad (2) $$

where σ₀<sub>ff</sub>(ω) is the ideal gas result (e.g., John 1994). Here δ<sub>ff</sub>(ω) can be computed as (Iglesias et al. 2002)

$$ \delta_{\text{ff}}(\omega) = \frac{\int_{0}^{\infty} I(k) dk}{\int_{0}^{\infty} I_{0}(k) dk}, \quad (3) $$

where

$$ I(k) = I_0(k) \frac{S(k)}{|\epsilon(\omega)|^2} \quad (4) $$

<footnote n="4">⁴ https://orcaforum.cec.mpg.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Absorption cross section of metal spectral lines. The black lines correspond to the Lorentzian profiles, and the red ones are the profiles obtained with the unified line shape theory of Allard et al. (1999). These line profiles were computed assuming $T = 6000$ K and $n_{\text{He}} = 10^{22} \text{ cm}^{-3}$. Note that the improved line profile for Ca I 4226 Å relies on approximate potentials (see text).</figcaption></figure>

and

$$ I_0(k) = \frac{1}{k} \exp \left[ -\frac{\hbar^2}{2m_e k_B T} \left( \frac{k}{2} - \frac{m_e \omega}{\hbar k} \right)^2 \right] \times \left| \frac{k^2 \mathcal{F}[\phi_{\text{e-He}}(r)]}{4\pi e^2} \right|^2. \quad (5) $$

In the last expressions, $\epsilon(\omega)$ is the dielectric function, $m_e$ and $e$ are the electron mass and charge, $k_B$ is the Boltzmann constant, $\hbar$ is the reduced Planck constant, and $\mathcal{F}[\phi_{\text{e-He}}(r)]$ is the Fourier transform of the electron–helium potential, for which we use the simple form given by Equations (3.5) and (3.6) of Iglesias et al. (2002). From these equations, it follows that two external inputs are needed to compute $\delta_{\text{ff}}(\omega)$: (1) the structure factor $S(k)$ and (2) the index of refraction of helium $n(\omega) = \sqrt{\epsilon(\omega)}$. The details regarding the calculation of the structure factor are given below, while our evaluation of the index of refraction is described in Section 2.5.

To compute $S(k)$, we rely on the classical fluid theory and the Ornstein–Zernike (OZ) equation. To solve the OZ equation, we use the Percus–Yevick closure relation (Percus & Yevick 1958), since it is well-suited for fluids dominated by short-range interactions (i.e., non-coulombic interactions; Hansen & McDonald 2006). The calculations are performed using a

modified version of pyOZ.<fnref n="5" /> Figure 2 compares our $S(0)$ values to the $S(0)$ analytical fit given in Rohrmann (2018). The agreement between both data sets is satisfactory under $\rho = 1 \text{ g cm}^{-3}$ ($n_{\text{He}} = 1.5 \times 10^{23} \text{ cm}^{-3}$) but worsens at higher densities. This disagreement reflects the limitations of the Percus–Yevick closure relation at high densities in a regime where the Monte Carlo calculations of Rohrmann (2018) are more appropriate. Nevertheless, this small discrepancy is of limited importance in the context of the modeling of cool DZ stars, since the photospheric density of our models never exceeds $\approx 1 \text{ g cm}^{-3}$.

## 2.5. Index of Refraction

The index of refraction, which is needed to compute the correction to the He⁻ free-free cross section (Equations (3) and (4)), is obtained from the Lorentz–Lorenz equation,

$$ \frac{n^2 - 1}{n^2 + 2} = A_R \left( \frac{n_{\text{He}} a_0^3}{N_A} \right) + B_R \left( \frac{n_{\text{He}} a_0^3}{N_A} \right)^2 + \mathcal{O}(n_{\text{He}}^3), \quad (6) $$

where $A_R$ and $B_R$ are the first and the second refractivity virial coefficients, $n_{\text{He}}$ is the helium number density, $a_0$ is the Bohr radius, and $N_A$ is the Avogadro constant. Here $A_R$ is

<footnote n="5">⁵ http://pyoz.vrbka.net</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Structure factor at $k = 0$ as a function of density and for different temperatures. The solid lines show the analytical fits obtained by Rohrmann (2018) from Monte Carlo calculations, and the circles show the results we found by solving the OZ equation.</figcaption></figure>

proportional to the atomic polarizability $\alpha(\omega)$ and is given by

$$A_R(\omega) = \frac{4\pi N_A \alpha(\omega)}{3}. \qquad (7)$$

To compute $A_R$, we use the helium polarizability values reported in Masili & Starace (2003). For the second refractivity virial coefficient, we rely on the classical statistical mechanics expression (e.g., Fernández et al. 1999)

$$B_R(\omega, T) = \frac{8N_A^2 \pi^2}{3} \int_0^\infty \Delta\alpha_{\text{ave}}(\omega, r) \exp\left[-\frac{\phi(r)}{k_B T}\right] r^2 dr, \quad (8)$$

where $\Delta\alpha_{\text{ave}}(\omega, r)$ is the interaction-induced isotropic polarizability and $\phi(r)$ is the helium–helium interatomic potential. To compute $\Delta\alpha_{\text{ave}}(\omega, r)$, we turn to the expansion

$$\Delta\alpha_{\text{ave}}(\omega, r) = \Delta\alpha_{\text{ave}}(0, r) + \omega^2\Delta S(-4, r) + \mathcal{O}(\omega^4), \quad (9)$$

where $\Delta\alpha_{\text{ave}}(0, r)$ is given in Hättig et al. (1999) and Maroulis (2000), and the Cauchy moment $\Delta S(-4, r)$ is given in Hättig et al. (1999). Finally, for the interaction potential $\phi(r)$ in Equation (8), we use the effective pair potential of Ross & Young (1986), which is calibrated to fit experimental data for high-density helium.

To validate our analytical model of the index of refraction, we compared its predicted values with the high-pressure experimental measurements of Dewaele et al. (2003). This comparison is shown in Figure 3 and reveals no significant deviation between our values and the laboratory measurements. Additionally, we checked that our index-of-refraction values are virtually identical to those obtained by Rohrmann (2018).

### 3. Equation of State and Chemical Equilibrium

In this section, we describe how the equation-of-state and chemical equilibrium calculations were modified to take high-density nonideal effects into account.

#### 3.1. Equation of State

The total number density and internal energy density in each atmospheric layer are computed using the ab initio equations of

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Index of refraction of helium as a function of density. The line corresponds to the results of our calculations with Equations (6)–(8), and the circles are the laboratory measurements extracted from Dewaele et al. (2003). For both data sets, $T = 300$ K and $\lambda = 6328$ Å.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Density of a helium medium as a function of pressure and temperature. The solid lines show the results found when using the equation of state of Becker et al. (2014), and the dashed lines correspond to the case where the ideal gas law is assumed.</figcaption></figure>

state for hydrogen and helium published by Becker et al. (2014). As in Blouin et al. (2017), we resort to the additive volume rule for mixed H/He compositions. The mass density $\rho(P, T)$ and internal energy density $u(P, T)$ are given by

$$\frac{1}{\rho_{\text{mix}}(P, T)} = \frac{X}{\rho_{\text{H}}(P, T)} + \frac{Y}{\rho_{\text{He}}(P, T)}, \qquad (10)$$

$$u_{\text{mix}}(P, T) = X u_{\text{H}}(P, T) + Y u_{\text{He}}(P, T), \qquad (11)$$

where X and Y are the mass fractions of hydrogen and helium, respectively.

For the densest cool DZ stars, the pressure at the photosphere exceeds $10^{11}$ dyn cm$^{-2}$. Under such conditions, using the ideal gas law can lead to an important overestimation of the density. In fact, as shown in Figure 4, the ideal gas density can be up to a factor of 5 greater than the value found when using the equation of state of Becker et al. (2014). Such a difference can have a significant effect on the computed atmosphere structure and the synthetic spectrum, since most nonideal
effects included in the code (e.g., detailed line profiles, distorted CIA profiles, high-density continuum opacities, and nonideal chemical equilibrium) are parameterized as functions of the density. For instance, using the ideal gas law would lead to an overestimation of the broadening of spectral lines due to an overestimation of the density of perturbing helium atoms.

### 3.2. Chemical Equilibrium

To compute the ionization equilibrium of helium, we rely on the chemical model proposed by Kowalski et al. (2007). Since it does not rely on any free parameter, this ionization equilibrium model is a major improvement over the occupation probability formalism (Hummer & Mihalas 1988; Mihalas et al. 1988) used in most white dwarf atmosphere codes. Compared to models where the ideal Saha equation is assumed, DZ models that include the helium ionization equilibrium of Kowalski et al. (2007) reach slightly lower densities in their deepest layers. This is the result of pressure ionization, which increases the electronic density and, in turn, the opacity. However, this effect is not as important as in metal-free atmospheres, since heavy elements provide the majority of free electrons and therefore govern the atmosphere structure.

We have also included a detailed description of the ionization equilibrium of heavy elements, which is the subject of Section 4.

## 4. Ionization Equilibrium of Heavy Elements

Properly characterizing the ionization equilibrium of heavy elements in the atmosphere of cool DZ stars is important from several perspectives. First, accurate ionization ratios are necessary to obtain the right spectral-line depths. For instance, in the case of a star that shows both Ca II H & K and Ca I 4226 Å in its spectrum, obtaining the right Ca II/Ca I ratio is a prerequisite for simultaneously reproducing all spectral lines. Moreover, in cool DZ stars, heavy elements provide most of the electrons. Therefore, a change in the ionization equilibrium of these trace species can influence other opacity sources (most importantly, He⁻ free-free) and hence the whole structure of the atmosphere.

Unlike the rest of the nonideal effects added to our atmosphere code, the equilibrium of heavy elements in the dense atmosphere of cool DZ stars has not yet been explored by other investigators using state-of-the-art methods. Therefore, we had to perform our own calculations before implementing this improved constitutive physics in our code. In this section, we first give some theoretical background and describe our strategy to compute the ionization equilibrium (Section 4.1). Then, results from our ab initio calculations are presented in Section 4.2 and applied to white dwarf atmospheres in Section 4.3.

### 4.1. Theoretical Framework

#### 4.1.1. The Chemical Picture

To tackle the problem of the ionization equilibrium of heavy elements in the dense atmosphere of cool white dwarfs, we rely on the chemical picture. In this approach, atoms, ions, and electrons are considered as the basic particles, and their interactions are modeled through interaction potentials. This is not as exact as the physical picture, where nuclei and electrons are the basic particles. However, using the chemical picture has

several advantages. Since this approach is semi-analytical, the results derived from it are more easily applicable in stellar atmosphere codes (especially regarding opacity calculations, where thousands of bound states must be taken into account to include the multitude of observed spectral lines). Moreover, it is easier to identify the contribution of every physical effect and thus gain a better physical insight of the problem at hand (Winisdoerffer & Chabrier 2005).

In the chemical picture, the ionization equilibrium problem is reduced to the minimization of the Helmholtz free energy $F(\{N_i\}, V, T)$ associated with a mixture made of species $\{N_i\}$ in a volume $V$ maintained at temperature $T$ (see, for instance, Fontaine et al. 1977; Magni & Mazzitelli 1979; Hummer & Mihalas 1988; Saumon & Chabrier 1992). The total Helmholtz free energy of a mixture of atoms, ions, and electrons can be expressed as the sum of the ideal free energy of the electron gas $F_e^{id}$, the ideal free energy of every ion from every species $F_{j,k}^{id}$, the contribution from the internal structure of bound species $F_{j,k}^{int}$, and the nonideal contribution related to the interaction between species $F^{nid}$,

$$F = F_e^{id} + \sum_j \sum_k F_{j,k}^{id} + \sum_j \sum_k F_{j,k}^{int} + F^{nid}, \quad (12)$$

where $k$ is an ionization state and $j$ is an atomic species.

Since $F$ must be minimized, $dF = 0$, and the ionization equilibrium of species $J$ between ionization states $K$ and $K + 1$ imposes

$$0 = \left( \frac{\partial F}{\partial N_e} \bigg|_{N_{j,k}, V, T} \right) dN_e + \left( \frac{\partial F}{\partial N_K} \bigg|_{N_e, N_{j,k}=K, V, T} \right) dN_K + \left( \frac{\partial F}{\partial N_{K+1}} \bigg|_{N_e, N_{j,k}=K+1, V, T} \right) dN_{K+1}, \quad (13)$$

which, by definition of the chemical potential, is equivalent to the condition

$$\mu_{J,K} = \mu_{J,K+1} + \mu_e. \quad (14)$$

Neglecting the interaction term $F^{nid}$ in Equation (12) and taking $F_e^{id}$ and $F_{j,k}^{id}$ to be the free energy of an ideal nonrelativistic, nondegenerate gas (Landau & Lifchitz 1980), Equation (14) leads to the well-known Saha equation,

$$\frac{n_{K+1}n_e}{n_K} = \frac{2Q_{K+1}}{Q_K} \left( \frac{2\pi m_e k_B T}{h^2} \right)^{3/2} e^{-I/k_B T}, \quad (15)$$

where $h$ is the Planck constant, $n_i$ are number densities, $Q_i$ are partition functions, and $I$ is the ionization potential.

Now, if we keep the nonideal terms in the free-energy equation, we find a result of the form of Equation (15) but with an effective ionization potential $I + \Delta I$ (Kowalski et al. 2007; Zaghloul 2009),

$$\frac{n_{K+1}n_e}{n_K} = \frac{2Q_{K+1}}{Q_K} \left( \frac{2\pi m_e k_B T}{h^2} \right)^{3/2} e^{-(I+\Delta I)/k_B T}, \quad (16)$$

where

$$\Delta I = \mu_e^{nid} + \mu_{K+1}^{nid} - \mu_K^{nid}. \quad (17)$$
Therefore, to compute the nonideal ionization equilibrium of heavy elements in dense helium-rich fluids, all that is needed is to compute the appropriate $\Delta I$ given by the above equation.

In Equation (17), it is the difference in free energy of many-body systems in thermodynamic equilibrium with different ionization states that is computed. This yields an effective ionization potential applicable to thermodynamic ionization equilibrium calculations. As emphasized by Crowley (2014), this ionization potential is not directly applicable to nonequilibrium processes (e.g., photoionization). These are fast (adiabatic) processes that occur before the surrounding plasma has any time to respond.

### 4.1.2. General Strategy

To compute $\Delta I$, we have to evaluate the nonideal chemical potential of every species involved in the ionization process. The electronic term $\mu_K^{\text{nid}}$ is already available in the literature. Kowalski et al. (2007) performed density functional theory (DFT) calculations to evaluate the excess energy of an electron embedded in a dense helium medium and found values that are in good agreement with existing laboratory measurements (Broomall et al. 1976). These calculations, published as polynomial expansions, were performed for a range of temperatures and densities suitable for our purpose.

While $\mu_{K+1}^{\text{nid}}$ and $\mu_K^{\text{nid}}$ were calculated by Kowalski et al. (2007) in the case of helium ionization, we are not aware of any study where the nonideal chemical potentials were computed for heavy elements surrounded by dense helium. The central task of this section is to compute these chemical potentials in order to obtain $\Delta I$ by virtue of Equation (17).

In the limit of strongly coupled systems, the role of entropy can be neglected for the calculation of thermodynamic equilibrium ionization potential, since the configuration of atoms remains the same before and after the ionization takes place. However, plasmas encountered in white dwarf atmospheres have a finite coupling strength. When an atom is ionized, the medium responds, and additional energy is transferred between the atom and the surrounding particles (Crowley 2014). Therefore, the nonideal chemical potential of a species in ionization state K can be expressed as the sum of two contributions,

$$ \mu_K^{\text{nid}} = E_K^{\text{exc}} + \mu_K^{\text{nid,ent}}, \quad (18) $$

where $E_K^{\text{exc}}$ is the excess of internal energy per particle and $\mu_K^{\text{nid,ent}}$ is the entropic contribution to the nonideal chemical potential. Note that this separation of $\mu_K^{\text{nid}}$ into two distinct components directly follows from the definition of the Helmholtz free energy. As $F = E + TS$ and $\mu_K^{\text{nid}} = (\partial F_K^{\text{nid}} / \partial N_K)|_{N_{K=K}, V, T}$, we can write

$$ \mu_K^{\text{nid}} = \left. \frac{\partial(E_K^{\text{nid}} + TS_K^{\text{nid}})}{\partial N_K} \right|_{N_{K=K}, V, T} = E_K^{\text{exc}} + \mu_K^{\text{nid,ent}}. \quad (19) $$

Our general strategy is summarized in Figure 5. To compute the $\mu_K^{\text{nid,ent}}$ contribution, we follow the work of Kowalski (2006b) and Kowalski et al. (2007) and use the classical fluid theory and OZ equation, as detailed in Section 4.2.1. To retrieve $E_K^{\text{exc}}$, we turn to DFT to compute the excess energy of a metallic ion embedded in a dense helium medium. This approach has the advantage of naturally taking into account many-body interaction terms. Prior to using DFT to compute $E_K^{\text{exc}}$, we use molecular dynamics (MD) simulations to obtain

<figure><img src="image_7.png" /><figcaption>Figure 5. Computational strategy used to retrieve the nonideal chemical potential of ionic species. The dashed arrow indicates a validation step described in Section 4.2.2.</figcaption></figure>

representative atomic configurations, as described in detail in Section 4.2.2.

### 4.1.3. Comparison with Previous Studies

To take into account the nonideal ionization of heavy elements, white dwarf atmosphere models (Koester & Wolff 2000; Wolff et al. 2002; Dufour et al. 2007) typically rely on the Hummer–Mihalas occupation probability formalism (Hummer & Mihalas 1988; Mihalas et al. 1988). In this framework, an occupation probability $w_i$ is assigned to every electronic level of every ion. If the level is unperturbed, $w_i = 1$; if the level is completely destroyed by interparticle interactions, $w_i = 0$. This occupation probability appears in the Boltzmann distribution and multiplies every term of the partition function,

$$ Q_K = \sum_i w_i K g_{iK} \exp \left( -\frac{e_{iK}}{k_B T} \right), \quad (20) $$

where the sum is over all states i of species K, and g is a statistical weight. To compute $w_i$ in the particular case of neutral interactions, Hummer & Mihalas (1988) used the second virial coefficient in the van der Waals equation of state to obtain

$$ w_i = \exp \left[ -\frac{4\pi}{3} \sum_{i'} n_{i'} (r_i + r_{i'})^3 \right], \quad (21) $$

where $n_i$ is the number density of the particles in state i and $r_i$ is the radius of the particles in this state. The interpretation of Equation (21) is straightforward: when a state occupies a volume of the same order as the mean volume allowed per
particle, it is gradually destroyed. Although simple and easy to implement in atmosphere models, we see three important drawbacks with this approach.

1. This formalism is expected to break down above $\approx 0.01 \text{ g cm}^{-3}$ (Hummer & Mihalas 1988), which is insufficient for many cool DZ white dwarfs.

2. The excluded volume effect is only a caricature of the real interaction potential between two neutral particles.

3. There is no theoretical prescription for the radii $r_i$. For instance, for a ground-state He I atom, should $r$ be given by the hydrogenic approximation ($r = n^2a_0/Z_{\text{eff}} = 0.39 \text{ Å}$) or the van der Waals radius (1.40 Å; Bondi 1964)? To address this problem, it is always possible to calibrate the radii to fit the spectral lines observed in white dwarf stars. This was successfully done by Bergeron et al. (1991) for hydrogen, but it would be impracticable for DZ stars, where many ions contribute to the total electronic density.

Our approach aims at answering these three concerns. First, by taking into account many-body interaction terms, it is designed to remain physically exact up to densities of the order of 1 g cm$^{-3}$. Second, the interaction between species is modeled through ab initio calculations that accurately describe the complex behavior of electrons under these high-density conditions. Finally, since we rely only on first-principles physics, our method does not require any free parameters.

### 4.1.4. Approximations

Before moving to the calculation of the nonideal chemical potentials and $\Delta I$ in Section 4.2, we take time to justify three important approximations that we use throughout Section 4.

*Electrons and heavy elements as trace species.* We are interested in helium-rich plasmas, where heavy elements and electrons can be considered as trace species. Hence, we completely neglect the interaction of metallic ions with other metallic ions and with electrons. This approximation is justified by the very low abundance of heavy elements in white dwarf atmospheres. Indeed, to our knowledge, the most metal-rich DZ star mentioned in the literature has an atmosphere with a number density ratio of log Ca/He $\approx -6$ (Ton 345; Wilson et al. 2015).

As a consequence of this approximation, we completely ignore the excess energy resulting from the interaction between charged species. Since electrostatic interactions occur at long range, this approximation deserves some additional justification. To show that electrostatic interactions are negligible, we computed the contribution of electrostatic interactions to the Helmholtz free energy. The latter can be broken down into three components (Chabrier & Potekhin 1998),

$$F^{\text{elec}} = F^{ee} + F^{ii} + F^{ie}, \qquad (22)$$

where $F^{ee}$ is the exchange-correlation contribution from the electron fluid, $F^{ii}$ is the contribution from the one-component ion plasma, and $F^{ie}$ is the electron screening contribution. To evaluate $F^{elec}$, we used the equations reported in Ichimaru et al. (1987) for $F^{ee}$ and those in Chabrier & Potekhin (1998) for $F^{ii}$ and $F^{ie}$. If all electrons originate from singly ionized species, then $F^{elec}$ is a function of only the electronic density $n_e$ and the plasma temperature

$T$. Figure 6 shows $\Delta I^{\text{elec}} = \left( \frac{\partial}{\partial N_e} + \frac{\partial}{\partial N_{j,i+1}} \right) F^{\text{elec}}$ for different

<figure><img src="image_8.png" /><figcaption>Figure 6. Contribution of the electrostatic interaction to the effective ionization potential with respect to the electronic density and temperature. The dashed line indicates the electronic density at $\tau_R = 2/3$ for vMa2, a typical cool DZ star.</figcaption></figure>

$n_e$ and $T$. The dashed line indicates the electronic density at the photosphere ($\tau_R = 2/3$) of vMa2, a typical cool DZ star. At these electronic densities and temperatures, the effect of electrostatic interactions on $\Delta I$ is only a few meV and is therefore negligible compared to the total $\Delta I$ reported later in this paper (which is of the order of a few eV). The charged-particle density is simply too low for electrostatic interactions to have any significant effect.

*Omission of the quantum behavior of ions.* We do not take into account the quantum behavior of ions and atoms. To justify this approximation, we can compute the first quantum correction of the Helmholtz free energy (Wigner 1932), which can be seen as a correction for the overlapping wave functions of nearby particles. For an $m$-component mixture, it can be expressed as (Saumon & Chabrier 1991)

$$F^{\text{quant}} = \frac{\pi \hbar^2}{12kT V} \sum_{a,b}^{m} \frac{N_a N_b}{\mu_{ab}} \int \nabla^2 \phi_{ab}(r) g_{ab}(r) r^2 dr, \quad (23)$$

where $\phi_{ab}(r)$ and $g_{ab}(r)$ are, respectively, the pair potential and pair distribution function between species $a$ and $b$, and $\mu_{ab} = \frac{m_a m_b}{m_a + m_b}$ is the reduced mass of particles $a$ and $b$. The contribution of this term to $\Delta I$ is computed as $\Delta I^{\text{quant}} = \left( \frac{\partial}{\partial N_{j,k+1}} - \frac{\partial}{\partial N_{j,k}} \right) F^{\text{quant}}$. Using the pair distribution functions and pair potentials described in Section 4.2.1, we find that $\Delta I^{\text{quant}}$ remains below 5 meV for all physical conditions relevant for the modeling of the atmosphere of cool DZ stars. As this is well below $E^{\text{exc}}$ and $\mu^{\text{nid,ent}}$, we can safely ignore the quantum behavior of ions.

*The ground-state approximation.* To compute the ionization equilibrium of heavy elements, we assume that every atom is in its electronic ground state. This solely means that we consider all species to be in their ground state when computing the ionization equilibrium. Once the ionization equilibrium is computed, the population of every electronic state can be obtained through the Boltzmann distribution. How good is this approximation? For helium atoms, this approximation is excellent. The first excited state of He I lies at 19.8 eV, so almost all helium atoms are in their fundamental state for the temperature domain in which we are interested ($k_B T < 1$ eV).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of the contributions of the fundamental state and the excited states to the partition function $Q$ at $k_B T = 0.5$ eV for heavy ions found in cool DZ stars. The number at the end of each bar gives the fraction of $Q$ resulting from excited states. This figure was made using the atomic data of the NIST Atomic Spectra Database (Kramida et al. 2015).</figcaption></figure>

For heavy elements, this approximation could be problematic. It is well known that excited states are typically more affected by nonideal effects than the fundamental state (e.g., Hummer & Mihalas 1988). Therefore, since the $\Delta I$ term in Equation (16) only takes into account the destruction of the fundamental state, an error could be introduced in the ionization equilibrium if excited states are affected in a significantly different way and if they account for a large portion of the partition function $Q$.

To investigate the maximum error associated with this approximation, we computed the fundamental-state contribution to the partition function $Q$ for C, Ca, Fe, Mg, and Na. The results are shown in Figure 7 for $k_B T = 0.5$ eV. The worst possible error associated with this approximation will occur if all excited states are destroyed while the fundamental state remains unperturbed (see Equation (20)). This scenario is highly unlikely but provides an easy way of assessing the maximum error. If it is the case, then, as shown in Figure 7, the maximum error on $Q$ is ≈40% (see Fe II). Therefore, in the worst case, the ionization fraction will be wrong by a factor of ≈2.

This maximum error is not a cause of concern for the modeling of the atmosphere of cool DZ stars. First, for all other atomic species (C, Ca, Mg, and Na), $Q$ is far more dominated by the fundamental-state contribution, and the maximal error associated with this approximation is thus much smaller than the value derived for Fe. Second, for the coolest DZ stars, the relative contribution of the fundamental state to the partition function is higher than for their warmer counterparts. Therefore, the ground-state approximation becomes more accurate for the stars for which the departure for the ideal chemical equilibrium is expected to be the most important. Last but not least, for the conditions relevant for the modeling of cool DZ stars, both this work and the formalism of Hummer & Mihalas (1988) predict deviations for the ideal gas equilibrium that are much more important than the aforementioned factor of ≈2 (see, for instance, Figure 15).

## 4.2. Results

In this section, we detail the computations performed to obtain $\Delta I$ for C, Ca, Fe, Mg, and Na. In Sections 4.2.1 and 4.2.2, we describe the computational setup and our intermediate results, and our final results are given in Section 4.2.3. For the sake of clarity, we only refer to Ca in the discussion of Sections 4.2.1 and 4.2.2, although all of the reported calculations were also performed for C, Fe, Mg, and Na.

### 4.2.1. Entropic Contribution

To compute the entropic contribution to the nonideal chemical potential, we first use the OZ equation (and the Percus–Yevick closure relation) to find the radial distribution function $g_{\text{He-Ca}}(r)$ describing the spatial configuration of Ca relative to He atoms. Then, once the radial distribution function $g_{\text{He-Ca}}(r)$ is obtained, $\mu_{\text{Ca}}^{\text{nid}}$ can be obtained through Equations (9) and (12) of Kiselyov & Martynov (1990). From there, we simply subtract the excess energy of Ca (as computed in the OZ framework) to obtain $\mu_{\text{Ca}}^{\text{nid,ent}}$ (Equation (18)).

To compute $g_{\text{He-Ca}}(r)$ with the OZ equation, the pair potentials $\phi_{\text{He-He}}(r)$ and $\phi_{\text{He-Ca}}(r)$ must be specified (in accordance with the approximation detailed in Section 4.1.4, $\phi_{\text{Ca-Ca}}(r) = 0$, since the metal–metal interactions are neglected). For the helium–helium pair potential, we use the effective pair potential of Ross & Young (1986).

As metal–helium pair potentials are not available in the literature for every metallic ion considered in this work, we had to compute ab initio pair potentials between helium and metallic ions. To do so, we used the ORCA quantum chemistry package to obtain the potential energy $\phi_{\text{Ca-He}}$ at various separations,

$$ \phi_{\text{Ca-He}}(r) = E_{\text{Ca-He}}(r) - E_{\text{He}} - E_{\text{Ca}}, \quad (24) $$

where $E_{\text{Ca-He}}(r)$ is the total energy for a separation $r$ and $E_{\text{He}}$ and $E_{\text{Ca}}$ are the computed energies of isolated He and Ca atoms. We rely on the CCSD(T) method (Raghavachari et al. 1989) as implemented in ORCA (Neese et al. 2009; Kollmar & Neese 2010) with the aug-cc-pCVQZ basis sets (Dunning 1989; Kendall et al. 1992; Woon & Dunning 1993). Using the counterpoise method (Boys & Bernardi 1970), we verified that the basis set superposition error is small enough ($<2$ meV) to be neglected for our purpose.

In the particular case of Ca, a few interaction potentials can be found in the literature for the Ca I–He I (Partridge et al. 2001; Lovallo & Klobukowski 2004) and Ca II–He I interactions (Czuchaj et al. 1996; Allard & Alekseev 2014). We used the values reported by these authors to validate our computational setup. This comparison, which reveals no significant differences, is shown in Figure 8.

The main limitation of these pair potentials is that they were obtained in the infinite-dilution limit (i.e., Ca interacts with only one He atom). Therefore, when we use these potentials, we implicitly assume that the total potential is pairwise additive, and an error may be introduced if many-body terms are important. This is the main reason why we resort to the OZ equation only to compute the entropic contribution and not to compute the excess energies. In fact, as described in Section 4.2.2, we turn to DFT to compute excess energies, which guarantees that many-body interaction terms are properly taken into account.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between the pair potentials for the Ca I–He I and Ca II–He I interactions computed in this work and the values reported in Lovallo & Klobukowski (2004), Partridge et al. (2001), Czuchaj et al. (1996), and Allard & Alekseev (2014).</figcaption></figure>

### 4.2.2. Excess Energy Contribution

The excess energy of Ca embedded in a dense helium medium made of N He atoms is given by

$$E_{\text{Ca-He}}^{\text{exc}} = E_{\text{NCa+He}} - E_{\text{NHe}} - E_{\text{Ca}}, \quad (25)$$

where $E_{\text{NCa+He}}$ is the total energy of the system, $E_{\text{NHe}}$ is the energy of the N He atoms, and $E_{\text{Ca}}$ is the computed energy of the isolated Ca atom. This calculation requires two steps. First, we need to find meaningful atomic configurations for the system (i.e., configurations that are representative of the thermodynamic fluctuations undergone by the real system). Then, we can use these configurations to compute the excess energy with Equation (25).

**MD.** To obtain representative atomic configurations of a system consisting of one Ca atom surrounded by N He atoms at a given temperature and density, we turned to classical MD simulations. More precisely, we used LAMMPS<fnref n="6" /> (Plimp- ton 1995) and the pair potentials described in Section 4.2.1. The simulations were performed in a cubic box with periodic boundary conditions. The box size and the number of He atoms included in the simulations were chosen to attain the desired density (additional considerations regarding finite-size effects are discussed in the next paragraph), and the temperature was kept near the target value using a Nosé–Hoover thermostat (Nosé 1984; Hoover 1985). The simulations were run for 5 ns using 0.2 fs time steps. At regular time intervals, the atomic positions were saved, and it is these configurations that we use in the next section to compute the excess energies.

**DFT calculations.** To compute the excess energy of Ca in the atomic configurations extracted from the MD simulations, we used the QUANTUM ESPRESSO<fnref n="7" /> DFT package (Giannozzi et al. 2009) with the PBE exchange-correlation functional (Perdew et al. 1996) and norm-conserving pseudopotentials. For all DFT calculations, we chose a kinetic energy cutoff of 45 Ry (612 eV) and a charge density cutoff of 180 Ry. We checked that this cutoff is enough to achieve a <0.05 eV

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Excess energy of Ca at T = 4000 K for configurations taken at 25 ps intervals from MD trajectories for different helium densities.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Autocorrelation function of the excess energy time series shown in Figure 9.</figcaption></figure>

convergence of the metal excess energy. To remove the electrostatic interaction associated with periodic boundary conditions, we used the Martyna–Tuckerman correction (Martyna & Tuckerman 1999) as implemented in QUANTUM ESPRESSO, which allows us to correct both the total energy and the self-consistent field potential.

Furthermore, to make sure that the finite size of the box does not result in undesired artifacts, we performed simulations using different numbers of helium atoms per simulation box and different box sizes (up to N = 160 helium atoms and a = 30 au). We found that using at least N = 50 helium atoms and a simulation box of at least a = 15 au (7.94 Å) allows a <0.1 eV convergence of the excess energy compared to results obtained at the same density with higher N and a values. This indicates that finite-size artifacts are negligible when these two conditions are met. Hence, all DFT calculations reported in this work were performed with a ≥ 15 au and N ≥ 50.

When computing the excess energy $E_{\text{exc}}$ using configuration snapshots extracted from MD simulations, the results can fluctuate drastically from one configuration to the other. This is shown in Figure 9, where the lines represent the evolution of $E_{\text{exc}}$ from configuration to configuration. In Figure 10, we show

<footnote n="6">⁶ http://lammps.sandia.gov</footnote>

<footnote n="7">⁷ http://quantum-espresso.org</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Standard error of the mean of the Ca excess energy at T = 4000 K with respect to the number of independent configurations used to compute the mean for different helium densities.</figcaption></figure>

the autocorrelation function of the $E_{exc}$ time series,

$$r_k = \frac{\sum_{i=1}^{N-k} (E_{exc}^i - \langle E_{exc} \rangle)(E_{exc}^{i+k} - \langle E_{exc} \rangle)}{\sum_{i=1}^{N} (E_{exc}^i - \langle E_{exc} \rangle)^2} \quad (26)$$

Since the autocorrelation function quickly decays to zero, we conclude that the time elapsed between each configuration snapshot is long enough for the $E_{exc}$ time-series values to be statistically independent. Therefore, we can safely apply the central-limit theorem to compute the standard error of the mean,

$$\sigma_{\langle E_{exc} \rangle} = \frac{\sigma_{E_{exc}}}{\sqrt{N}} \quad (27)$$

Figure 11 shows the evolution of $\sigma_{\langle E_{exc} \rangle}$ with respect to the number of configurations used to compute the mean. For both $\rho = 0.1$ and $\rho = 1.0 \text{ g cm}^{-3}$, we notice the $1/\sqrt{N}$ decay of $\sigma_{\langle E_{exc} \rangle}$. This implies that to improve the error by a factor of 2, the number of configurations needs to be quadrupled. From this analysis, we chose to use 100 configurations for each ($T, \rho$) condition. This value is enough to obtain $\sigma_{\langle E_{exc} \rangle} \lesssim 0.1 \text{ eV}$ for most physical conditions considered in this work, which is an error that we consider acceptable for our purpose.

Validation with ab initio MD. Since our $\phi_{Ca-He}(r)$ potential was calculated in the infinite-dilution limit, one could be worried about the exactitude of the atomic configurations obtained through MD using this potential. To check this point, we computed the excess energy of Ca using configurations extracted from ab initio MD simulations. In this framework, no pair potential is assumed. The electronic density, energy, and forces on ions are recomputed at every time step of the simulation using DFT. This approach is expected to be more exact than the classical MD approach, but its computational cost is larger by orders of magnitude. These calculations were performed using Born-Oppenheimer MD with the CPMD package<fnref n="8" /> (Marx & Hutter 2000; Hutter et al. 2008), with the PBE exchange-correlation functional and ultrasoft pseudopotentials (Vanderbilt 1990). We employed 0.5 fs time steps and

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Excess energy of Ca at T = 5000 K for different helium densities, obtained from configurations extracted from ab initio MD (DFT-MD) and classical MD using the pair potentials described in Section 4.2.1.</figcaption></figure>

an energy cutoff of 35 Ry. As before, we extracted atomic configurations from these simulations and used these configurations to compute the interaction energy of Ca with the surrounding medium through DFT calculations.

Figure 12 compares the results obtained to those found with the classical MD simulations. This comparison shows that there is only a negligible difference between the two approaches, at least below $\rho = 1 \text{ g cm}^{-3}$. We did not perform any comparison at higher densities because of the prohibitive calculation time of such calculations. In any case, densities above $1 \text{ g cm}^{-3}$ are never encountered at the photosphere of cool DZ white dwarfs (Section 4.3). Therefore, we conclude that our infinite-dilution limit potential $\phi_{Ca-He}(r)$ is sufficient to generate the atomic configurations used to compute the excess energy (and it is much faster than resorting to ab initio MD simulations).

### 4.2.3. Ionization Equilibrium

Following the methodology described in the previous sections, we computed $\mu_K^{\text{nld,ent}}$ and $E_K^{\text{exc}}$ for C I/C II, Ca I/Ca II, Fe I/Fe II, Mg I/Mg II, and Na I/Na II. By adding these excess chemical potentials to the electron excess energy, we computed how much the ionization potential is altered at a given density and temperature (Equation (17)). Figure 13, which shows the three contributions to $\Delta I$ (free electron excess energy, variation of $E_K^{\text{exc}}$, and change in $\mu_K^{\text{nld,ent}}$), illustrates this process in the case of Ca.

Figure 14 shows our final results. First, for every ion considered, we notice that $\Delta I \to 0$ when $\rho \to 0$. This is the expected behavior, and it shows that our methodology is consistent with the ideal regime when we push it to low densities. Second, we note that $\Delta I$ is always negative and that its absolute value increases with density. This result means that ionization becomes easier with increasing density, which also corresponds to the expected behavior. Finally, for all elements except Fe, we notice that higher temperatures are associated with slightly larger ionization potential depressions. This result is consistent with the findings of Kowalski et al. (2007), who found a reduction of the band gap of warm dense helium with increasing temperature.

<footnote n="8">⁸ http://cpmd.org</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 13. Contributions added to the reference ionization potential of Ca to obtain its effective ionization potential at various densities (see legend). These results are for T = 4000 K.</figcaption></figure>

To easily implement these nonideal ionization potentials in atmosphere models, we have fitted our results with a simple function of $\rho$ and $T$,

$$ \Delta I(\rho, T) = \min\{0, (a + bT)\rho + c\rho^2\}, \quad (28) $$

where a, b, and c are parameters found using a $\chi^2$ minimization algorithm, $\rho$ is the helium density in g cm$^{-3}$, and $T$ is the temperature in K. This expression allows a satisfactory fit to the data and yields $\Delta I = 0$ at $\rho = 0$. The analytical fits are shown in Figure 14, and the fitting parameters are reported in Table 2. Formally, in order to stay within the limits of our calculations, the use of these analytical expressions should be limited to densities between 0 and 1.5 g cm$^{-3}$ and temperatures between 4000 and 8000 K. Nevertheless, we have verified that Equation (28) can safely be extrapolated to lower (down to 2000 K) or higher (at least up to 10,000 K) temperatures if needed.

### 4.2.4. Comparison with Previous Studies

It is instructive to compare these results with the ionization equilibrium predicted by the Hummer–Mihalas occupation probability formalism, which is widely used in atmosphere codes. Since there is no theoretical prescription for the values of the hard sphere radii used to compute the occupation probabilities (Equation (21)), a somewhat arbitrary choice must be made to perform this comparison. We chose to compute the hard sphere radii with the hydrogenic approximation, as described by Beauchamp (1995). In this approximation, the radius of a species in state i is given by

$$ r_i = \frac{n_i^2 a_0}{Z_i^{\text{eff}}}, \quad (29) $$

where $n_i$ is the principal quantum number of the uppermost electron, $a_0$ is the Bohr radius, and the effective nuclei charge $Z_i^{\text{eff}}$ is given by

$$ Z_i^{\text{eff}} = n_i \sqrt{\frac{I_i}{13.598 \text{ eV}}}, \quad (30) $$

where $I_i$ is the energy needed to ionize an electron from state i. In the Hummer–Mihalas formalism, every term in the partition

function is multiplied by the occupation probability (Equation (20)). If we stick to the ground-state approximation (Section 4.1.4), the occupation probability is the same for every level and can be factored out of the partition-function sum. Hence, the net effect of the Hummer–Mihalas formalism is to multiply the right-hand side of the Saha equation (Equation (15)) by the ratio of occupation probabilities, $w_{ZII}/w_{ZI}$.

Figure 15 compares the multiplicative factors that need to be applied to the right-hand side of the Saha equation for the Ca I/ Ca II ionization equilibrium to account for nonideal effects (i.e., $w_{CaII}/w_{CaI}$ in the case of the Hummer–Mihalas formalism and $e^{-\Delta I/(k_B T)}$ for our ionization model). The most obvious aspect of Figure 15 is that we find a weaker pressure ionization than what is predicted using the Hummer–Mihalas formalism and hard sphere radii computed in the hydrogenic approximation. We checked that this result holds true for C, Fe, Mg, and Na. This conclusion is consistent with the findings of Bergeron et al. (1991) for the ionization equilibrium of hydrogen in cool DA stars. Using the Hummer–Mihalas formalism and a hydrogen radius given by $r_n = n^2 a_0$, they found that the high Balmer lines are predicted to be too weak, indicating that pressure ionization in the Hummer–Mihalas formalism is too strong. They showed that using a smaller radius in the computation of the occupation probabilities, $r_n = 0.5n^2 a_0$, allows better spectral fits.

Unfortunately, we cannot compute the ionization potential depression of H to directly confirm the conclusion of Bergeron et al. (1991). The problem is that the H II–He potential (e.g., Kolos & Peek 1976; Pachucki 2012) has a deep attractive well (since $H^+$ and He can form the He$H^+$ molecule) that prevents proper convergence of the OZ equation solver. The same issue arises if we try to compute the ionization potential of H in an H-rich medium, since the H II–H I potential (e.g., Frost & Musulin 1954) also has an important attractive well ($H^+$ and H can form the $H_2^+$ molecule).

### 4.3. Atmosphere Models

Using the analytical model described in the previous section, we implemented the improved ionization equilibrium of heavy elements in our atmosphere code to investigate how it affects the synthetic spectra of cool DZ stars. Before even examining any spectrum, we can get an idea of the impact of the new nonideal ionization equilibrium by looking at the densities involved in the model atmospheres. Figure 16 shows the density at $\tau_\nu = 2/3$ as a function of $\lambda$ for a few atmosphere models with different effective temperatures and calcium abundances.<fnref n="9" /> This type of figure is useful to identify which densities are probed at different wavelengths. In the previous section, we saw that no important deviation from the ideal ionization equilibrium is expected below 0.1 g cm$^{-3}$ (see Figure 14). From Figure 16, it is clear that the probed densities are below this threshold for Ca/He $\gtrsim 10^{-10}$ and above this threshold for Ca/He $\lesssim 10^{-10}$. Therefore, it should become important to take into account the nonideal ionization equilibrium for cool DZ atmosphere models with Ca/He $\lesssim 10^{-10}$, but it is probably superfluous for models with Ca/He $\gtrsim 10^{-10}$ (note that nonideal effects on the opacities and the equation of state nevertheless remain

<footnote n="9">⁹ In this paper, the abundance of all metallic species, from C to Cu, is scaled to the abundance of Ca to match the abundance ratios of chondrites reported in Lodders (2003).</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Depression of the ionization potential of C, Ca, Fe, Mg, and Na embedded in a dense helium fluid. Circles show the results of our ab initio calculations, and error bars indicate the statistical errors associated with the configuration sampling. The solid lines show the analytical fits found through a $\chi^2$ minimization of Equation (28). Data in red are for T = 8000 K, and data in yellow are for T = 4000 K.</figcaption></figure>

<table><caption>**Table 2**
Fitting Parameters for $\Delta I(\rho, T)$ (Equation (28))</caption><thead><tr><th>Ion</th><th>a<sup>a</sup></th><th>b<sup>b</sup></th><th>c<sup>c</sup></th></tr></thead><tbody><tr><td>C</td><td>1.91782</td><td>-3.24813</td><td>-1.19948</td></tr><tr><td>Ca</td><td>-2.20703</td><td>-0.14431</td><td>0.57494</td></tr><tr><td>Fe</td><td>-2.23142</td><td>0.48427</td><td>0.21301</td></tr><tr><td>Mg</td><td>0.45809</td><td>-0.85522</td><td>-1.01958</td></tr><tr><td>Na</td><td>-0.52305</td><td>-0.62471</td><td>0.04833</td></tr></tbody></table>

**Notes.**

<sup>a</sup> eV g⁻¹ cm³.

<sup>b</sup> 10⁻⁴ eV g⁻¹ K⁻¹ cm³.

<sup>c</sup> eV g⁻² cm⁶.

important in this regime). For intermediate densities (Ca/He ≈ 10⁻¹⁰), using the nonideal ionization equilibrium should result in small changes in the spectral-line wings of the coolest models.

Figure 17 compares synthetic spectra computed with our ionization equilibrium model to spectra computed using the occupation probability formalism and the ideal Saha equilibrium (in each case, the atmosphere model structure and the synthetic spectrum were computed using the same ionization model). This figure focuses on the region between 3500 and 4500 Å, since it contains several Ca, Fe, and Mg absorption lines likely to be affected by the choice of the ionization model.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Multiplicative factor applied to the right-hand side of the Ca I/Ca II Saha equation (Equation (15)) to take nonideal effects into account. The blue line is $w_{CaII}/w_{CaI}$, the result obtained using the Hummer–Mihalas formalism, and the green curve is $e^{-\Delta I/(k_B T)}$, the result obtained with our ionization model.</figcaption></figure>

The first thing to note is that for the high-density models (i.e., those with a low metal abundance and effective temperature), there are important differences between spectra obtained using the ideal Saha equilibrium and our ionization model. These
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Density at an optical depth $\tau_\nu = 2/3$ with respect to $\lambda$. The top panel shows the results for $T_{\text{eff}} = 4000$ K models and the bottom panel for $T_{\text{eff}} = 6000$ K. The Ca abundance is given in the legend, and a surface gravity $\log g = 8$ is assumed.</figcaption></figure>

differences are mostly due to a shift in the continuum associated with the increased electronic density in models that take pressure ionization into account. Next, we notice that the spectra computed using the Hummer–Mihalas formalism are even further from the spectra obtained using the ideal Saha equilibrium than those computed with our ionization model. This is not surprising, since, as seen in Figure 15, the Hummer–Mihalas formalism predicts a very strong pressure ionization. Finally, for the low-density models (i.e., those with a high metal abundance and/or effective temperature), all three sets of spectra are virtually identical, which is consistent with our analysis of Figure 16.

The nonideal chemical equilibrium of heavy elements also has a small impact on the model atmosphere structure. The increased electronic density associated with pressure ionization leads to an increase of the Rosseland mean opacity and therefore a reduction of the pressure at the photosphere. For instance, for $T_{\text{eff}} = 4000$ K, $\log g = 8$, and Ca/He = $10^{-11}$, a model that assumes the ideal Saha equation has a photospheric density of 0.93 g cm$^{-3}$, while an atmosphere structure based on our ionization model has a photospheric density of 0.89 g cm$^{-3}$. Moreover, the occupation probability formalism predicts a density that is still lower ($10^{-8}$ g cm$^{-3}$). Given Figure 15, this result is not surprising: compared to our calculations, the Hummer–Mihalas formalism overestimates the efficiency of pressure ionization.

Our results constitute a physically grounded answer to the question of the importance of pressure ionization in cool DZ stars, which will help to reduce the gap between solutions found with different atmosphere codes. A good example to illustrate this point is vMa2 (WD 0046+051). On one hand, using an ideal treatment of chemical equilibrium, Dufour et al.

(2007) found a solution with $T_{\text{eff}} = (6220 \pm 240)$ K. On the other hand, using the Hummer–Mihalas occupation probability formalism, Wolff et al. (2002) found $T_{\text{eff}} = (5700 \pm 200)$ K. In their analysis, Dufour et al. (2007) showed that the difference between both solutions can largely be explained by the different chemical equilibrium models used in both studies. This uncertainty can be removed by relying on the accurate description of the chemical equilibrium described in the current work.

## 5. Applications

To show how the improved constitutive physics presented in this work translates in terms of better spectroscopic fits, this section presents the analysis of two well-known DZ stars: Ross 640 (WD 1626+368) and LP 658-2 (WD 0552-041). Applications to other objects will be presented in other papers of the series.

Our new analysis of these two objects makes use of Gaia DR2 parallaxes (Prusti et al. 2016; Brown et al. 2018), BVRI and JHK photometry published in Bergeron et al. (2001; see Table 3), optical spectra published in Giammichele et al. (2012), and UV spectra obtained with HST and the Faint Object Spectrograph (FOS; Koester & Wolff 2000; Wolff et al. 2002).

### 5.1. Ross 640

At $T_{\text{eff}} \approx 8000$ K, Ross 640 is technically not a “cool” white dwarf. Since the density at its photosphere is $\approx 0.01$ g cm$^{-3}$ ($n_{\text{He}} = 1.5 \times 10^{21}$ cm$^{-3}$), nonideal effects affecting the equation of state and chemical equilibrium are minimal. However, this density is high enough to induce important differences between Lorentzian profiles and the improved line profiles presented in Section 2.1. This object is therefore the perfect candidate to test our line profiles separately, without the interference of other nonideal effects.

To fit this star, we follow the procedure described in Dufour et al. (2007). In short, we first find $T_{\text{eff}}$ and $\log g$ using the photometric technique described in Bergeron et al. (2001). The photometric measurements are first converted into fluxes using the constants reported in Holberg & Bergeron (2006). Then, these observed fluxes $f_\nu$ are compared to the model fluxes $H_\nu$ to obtain $T_{\text{eff}}$ and the solid angle $\pi(R/D)^2$, where R is the radius of the star and D is its distance to the Earth. These parameters are found using a $\chi^2$ minimization technique relying on the Levenberg–Marquardt algorithm. Since D is known from the parallax measurement, the radius R can be computed from the solid angle. The mass of the star and the corresponding surface gravity $g = GM/R^2$ are then found using the evolutionary models of Fontaine et al. (2001). This $\log g$ value being generally different from our initial guess, we repeat the fitting procedure until all fitting parameters are converged.

Once a consistent solution for $T_{\text{eff}}$ and $\log g$ is obtained from the procedure described in the previous paragraph, we move to the determination of the abundances using spectroscopic observations. We keep $T_{\text{eff}}$ and $\log g$ fixed to the values found using the photometric observations and then fit the Ca/He and H/He ratios by minimizing the $\chi^2$ between our synthetic spectra and the observed spectrum. Since the abundances found with this technique are generally different from those initially used for the photometric fit, the whole fitting procedure is repeated until internal consistency is reached.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Comparison between synthetic spectra computed using the Hummer & Mihalas (1988) formalism (blue), the ionization equilibrium presented in this work (red), and the ideal Saha equation (black). All models were computed assuming log g = 8 and H/He = 0. The effective temperature and metal abundance are indicated above each panel.</figcaption></figure>

<table><caption>**Table 3**
Observational Data</caption><thead><tr><th></th><th>Ross 640</th><th>LP 658-2</th></tr></thead><tbody><tr><td>Parallax (mas)</td><td>62.915 ± 0.022</td><td>155.250 ± 0.029</td></tr><tr><td>B<sup>a</sup></td><td>14.02</td><td>15.49</td></tr><tr><td>V</td><td>13.83</td><td>14.45</td></tr><tr><td>R</td><td>13.75</td><td>13.99</td></tr><tr><td>I</td><td>13.66</td><td>13.54</td></tr><tr><td>J</td><td>13.58</td><td>13.05</td></tr><tr><td>H</td><td>13.57</td><td>12.86</td></tr><tr><td>K</td><td>13.58</td><td>12.78</td></tr></tbody></table>

**Note.**

Although the abundance ratio between the different heavy elements is kept constant during the $\chi^2$ minimization procedure, we manually adjust the abundance ratio of Mg, Fe, and Si to fit the spectral lines labeled in Figure 18. All other heavy elements (from C to Cu) are included in the models, but since we could not use any spectral line to fit their abundances, we simply assume the same abundance ratio with respect to Ca as in chondrites (Lodders 2003).

As shown in Figure 18, our solution is consistent with observations across all wavelengths. Our fitting parameters, given in Table 4, are roughly similar to those found by Dufour et al. (2007), Koester & Wolff (2000), and Zeidler-KT et al. (1986), although they all found a higher effective temperature (8440±320, 8500±200, and 8800 K, respectively). One major improvement compared to previous authors is our fit to the broad Mg II 2795/2802 Å lines. To obtain a good fit, Koester & Wolff (2000) arbitrarily multiplied the van der Waals broadening constant of these lines by 10. No arbitrary constants are needed using our new line profiles, and a consistent abundance is found from both the optical and ultraviolet magnesium lines.

## 5.2. LP 658-2

LP 658-2 is a DZ star that exhibits a weak Ca II H & K doublet. During the last two decades, many authors have tried

to fit this star, but none has reached a consistent solution across all wavelengths. Because they relied on different models and observations, the solutions they found are quite diverse (see Table 5).

First, Bergeron et al. (2001) found that LP 658-2 has a helium-rich atmosphere with $T_{eff} = (5060 \pm 60)$ K. However, their analysis was based on atmosphere models that did not include heavy elements, which strongly influence UV opacities and the temperature profile.

Then, using HST data (FOS), Wolff et al. (2002) extended the analysis of Bergeron et al. (2001) with an investigation of the UV portion of the spectrum of LP 658-2. The large absorption feature observed in the UV was interpreted as strong broadening from the wing of Ly$\alpha$. Keeping the effective temperature fixed at the $T_{eff} = 5060$ K value found by Bergeron et al. (2001), they used this UV absorption feature to fit the hydrogen abundance and found that $H/He = 5 \times 10^{-4}$. However, contrary to other stars in their sample (e.g., LHS 1126 and BPM 4729), they were not able to properly reproduce the shape of this UV absorption feature.

Subsequently, using models that include heavy elements in the atmosphere structure, Dufour et al. (2007) determined a much cooler temperature for LP 658-2 ($T_{eff} = 4270 \pm 70$ K). At this temperature, the photometric data can completely exclude the presence of traces of hydrogen at the level found by Wolff et al. (2002), since H$_2$-He CIA would cause a strong IR flux depletion that is not observed. However, the solution of Dufour et al. (2007) does not explain the UV absorption feature seen in the FOS data, and their spectroscopic solution predicted a large Ca I 4226 Å line, which is completely absent from the observations.

More recently, Giammichele et al. (2012) argued that the narrow H & K lines observed in the spectra of LP 658-2 indicate that it is perhaps a hydrogen-rich star after all. However, although an H-rich composition allowed a better fit to the visible spectrum than that of Dufour et al. (2007), the photometric fit was not as good (and it cannot explain the shape of the UV spectrum).

<footnote><sup>a</sup> There is a 3% uncertainty on all photometric measurements.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Our best solution for Ross 640. The top panel shows our fit to the UV spectrum, the middle panel is our fit to the visible spectrum, and the bottom panel shows our photometric fit to the BVRI and JHK bands.</figcaption></figure>

<table><caption>**Table 4**
Fitting Parameters</caption><thead><tr><th></th><th>Ross 640</th><th>LP 658-2</th></tr></thead><tbody><tr><td>T<sub>eff</sub> (K)</td><td>8070 ± 140</td><td>4430 ± 40</td></tr><tr><td>log g</td><td>7.923 ± 0.008</td><td>7.967 ± 0.022</td></tr><tr><td>log H/He</td><td>-3.5 ± 0.2</td><td>&lt;-5</td></tr><tr><td>log Ca/He</td><td>-9.12 ± 0.05</td><td>-11.38 ± 0.05</td></tr><tr><td>log Fe/He</td><td>-8.44 ± 0.10</td><td>...</td></tr><tr><td>log Mg/He</td><td>-7.40 ± 0.10</td><td>-8.66 ± 0.20</td></tr><tr><td>log Si/He</td><td>-7.90 ± 0.20</td><td>...</td></tr></tbody></table>

<table><caption>**Table 5**
Literature Review of LP 658-2</caption><thead><tr><th>Authors</th><th>T<sub>eff</sub> (K)</th><th>H/He</th></tr></thead><tbody><tr><td>Bergeron et al. (2001)</td><td>5060 ± 60</td><td>He</td></tr><tr><td>Wolff et al. (2002)</td><td>5060 ± 60</td><td>H/He = 5 × 10<sup>-4</sup></td></tr><tr><td>Dufour et al. (2007)</td><td>4270 ± 70</td><td>He</td></tr><tr><td>Giammichele et al. (2012)</td><td>5180 ± 80</td><td>H</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Our best solution for LP 658-2. The top panel shows our fit to the visible spectrum, and the bottom panel displays our fit to the photometric observations and FOS data. The bottom panel also shows two synthetic UV spectra computed without the Mg II 2795/2802 Å and Mg I 2852 Å lines, one without hydrogen (blue) and one with H/He = 10⁻⁵ (green).</figcaption></figure>

Using our improved models, we can now obtain a solution that agrees perfectly with the observations across all wavelengths assuming a helium-rich atmosphere (Figure 19). We can also constrain the amount of hydrogen to H/He < 10⁻⁵, as a higher hydrogen abundance would produce an IR flux depletion that is incompatible with the observations. Given this limit, the shape of the UV continuum can no longer be explained by the wing of Lyα (see the green dash-dotted line in Figure 19). Instead, we find that the absorption in the UV can naturally be explained by the presence of trace amounts of magnesium (absorption from the Mg II 2795/2802 Å and Mg I 2852 Å lines). While there are no lines formally detected, the amount of magnesium needed to reproduce the UV continuum is small enough as to not produce features in the optical spectrum.

Finally, our new models do not predict the strong Ca I 4226 Å line that was predicted using the models of Dufour et al. (2007). This is mainly due to the use of our improved line profiles (Section 2.1), as well as our new nonideal Ca ionization equilibrium calculation (Section 4), the former effect being the most important. Our fitting parameters, given in Table 4, were found using the same fitting procedure as for Ross 640.

## 6. Conclusion

We have developed an updated atmosphere model code that incorporates all the necessary constitutive physics for an accurate description of cool DZ stars. This code includes

* the most important heavy-element line profiles computed using the unified line shape theory of Allard et al. (1999),
* CIA profiles suitable for fluids where the density exceeds 0.1 g cm⁻³,

* He Rayleigh scattering and He⁻ free-free absorption corrected for collective interactions between atoms,

* an ab initio equation of state for H and He, and

* a nonideal chemical equilibrium model for He, C, Ca, Fe, Mg, and Na.

While most of these nonideal effects were implemented using results previously published by various authors, we performed our own calculations to assess the chemical equilibrium of heavy elements.

More precisely, we used the classical theory of fluid and DFT calculations to characterize the ionization equilibrium of C, Ca, Fe, Mg, and Na in a dense helium medium and under the temperature and density conditions found in the atmosphere of cool DZ stars. These calculations show that the effective ionization potential begins to decrease when the density exceeds 0.1 g cm⁻³, reaching a depression of ≈1–2 eV at ρ = 1 g cm⁻³. We provided analytical fits to our data that can be implemented in atmosphere model codes to obtain the effective ionization potential for a given temperature and density.

We computed atmosphere models using this improved description of the ionization of heavy elements and found that under the right conditions (i.e., weakly polluted, low-T<sub>eff</sub> objects), the synthetic spectrum can significantly differ from results obtained using the ideal Saha equation. Moreover, we found that the Hummer–Mihalas formalism—when used in conjunction with hydrogenic hard sphere radii—leads to a much stronger pressure ionization than our model, which indicates an overestimation of pressure ionization. This result is consistent with previous findings based on comparisons between atmosphere models and observed spectra (Bergeron et al. 1991). Finally, we showed how the improved constitutive physics included in our code translates into better spectral fits for Ross 640 and LP 658-2, two cool DZ stars that presented a challenge to previous atmosphere model codes.

In the next papers of this series, we will use our updated models to analyze in detail other cool white dwarfs, in particular WD 2356–209 (a peculiar cool DZ star showing an exceptionally strong Na D feature) and the first cool DZ star to show CIA absorption. We will also analyze the bulk of the known cool white dwarfs taking advantage of the Gaia data and revisit the spectral evolution of these objects.

We wish to thank Piotr M. Kowalski for useful discussions regarding the DFT calculations presented in Section 4. This work was supported in part by NSERC (Canada).

This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This work has made use of the Montreal White Dwarf Database (Dufour et al. 2017).

This work used observations made with the NASA/ESA Hubble Space Telescope and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope

European Coordinating Facility (ST-ECF/ESA), and the Canadian Astronomy Data Centre (CADC/NRC/CSA).

Software: CPMD (Marx & Hutter 2000), LAMMPS (Plimp- ton 1995), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), ORCA (Neese 2012), QUANTUM ESPRESSO (Giannozzi et al. 2009).

## ORCID iDs

S. Blouin https://orcid.org/0000-0002-9632-1436

## References

* Allard, N., & Alekseev, V. 2014, AdSpR, 54, 1248

* Allard, N., Guillon, G., Alekseev, V., & Kielkopf, J. 2016a, A&A, 593, A13

* Allard, N., Leininger, T., Gadéa, F., Brousseau-Couture, V., & Dufour, P. 2016b, A&A, 588, A142

* Allard, N., Royer, A., Kielkopf, J., & Feautrier, N. 1999, PhRvA, 60, 1021

* Allard, N. F., Homeier, D., Guillon, G., Viel, A., & Kielkopf, J. 2014, JPhCS, 548, 012006

* Beauchamp, A. 1995, PhD thesis, Université de Montréal, Montréal, QC

* Becker, A., Lorenzen, W., Fortney, J. J., et al. 2014, ApJS, 215, 21

* Bergeron, P., Leggett, S., & Ruiz, M. T. 2001, ApJS, 133, 413

* Bergeron, P., Ruiz, M., & Leggett, S. 1997, ApJS, 108, 339

* Bergeron, P., Ruiz, M. T., Hamuy, M., et al. 2005, ApJ, 625, 838

* Bergeron, P., Saumon, D., & Wesemael, F. 1995, ApJ, 443, 764

* Bergeron, P., Wesemael, F., & Fontaine, G. 1991, ApJ, 367, 253

* Blouin, S., Kowalski, P., & Dufour, P. 2017, ApJ, 848, 36

* Bondi, A. 1964, JPhCh, 68, 441

* Boys, S. F., & Bernardi, F. D. 1970, MolPh, 19, 553

* Broomall, J. R., Johnson, W. D., & Onn, D. G. 1976, PhRvB, 14, 2819

* Brown, A., Vallenari, A., Prusti, T., et al. 2018, arXiv:1804.09365

* Chabrier, G., & Potekhin, A. Y. 1998, PhRvE, 58, 4941

* Crowley, B. 2014, HEDP, 13, 84

* Czuchaj, E., Rebentrost, F., Stoll, H., & Preuss, H. 1996, CP, 207, 51

* Dalgarno, A. 1962, Spectral Reflectivity of the Earth's Atmosphere III: The Scattering of Light by Atomic Systems, Geophysics Corporation of America Technical Report, 62-28-A

* Dewaele, A., Eggert, J., Loubeyre, P., & Le Toullec, R. 2003, PhRvB, 67, 094112

* Dufour, P., Blouin, S., Coutu, S., et al. 2017, in ASP Conf. Proc. 509, 20th European White Dwarf Workshop, ed. P.-E. Tremblay, B. Gaensicke, & T. Marsh (San Francisco, CA: ASP), 3

* Dufour, P., Bergeron, P., Liebert, J., et al. 2007, ApJ, 663, 1291

* Dunning, T. H. 1989, JChPh, 90, 1007

* Fernández, B., Hättig, C., Koch, H., & Rizzo, A. 1999, JChPh, 110, 2872

* Fontaine, G., Brassard, P., & Bergeron, P. 2001, PASP, 113, 409

* Fontaine, G., Graboske, H., Jr., & Van Horn, H. 1977, ApJS, 35, 293

* Frost, A. A., & Musulin, B. 1954, JChPh, 22, 1017

* Giammichele, N., Bergeron, P., & Dufour, P. 2012, ApJS, 199, 29

* Gianninas, A., Curd, B., Thorstensen, J. R., et al. 2015, MNRAS, 449, 3966

* Giannozzi, P., Baroni, S., Bonini, N., et al. 2009, JPCM, 21, 395502

* Hansen, J.-P., & McDonald, I. R. 2006, Theory of Simple Liquids (3rd ed.; New York: Academic)

* Hättig, C., Larsen, H., Olsen, J., et al. 1999, JChPh, 111, 10099

* Holberg, J., & Bergeron, P. 2006, ApJ, 132, 1221

* Homeier, D., Allard, N., Allard, F., et al. 2005, in ASP Conf. Proc. 334, 14th European Workshop on White Dwarfs, ed. D. Koester & S. Moehler (San Francisco, CA: ASP), 209

* Homeier, D., Allard, N., Johnas, C. M. S., Hauschildt, P. H., & Allard, F. 2007, in ASP Conf. Proc. 372, 15th European Workshop on White Dwarfs, ed. R. Napiwotzki & M. R. Burleigh (San Francisco, CA: ASP), 277

* Hoover, W. G. 1985, PhRvA, 31, 1695

* Hummer, D., & Mihalas, D. 1988, ApJ, 331, 794

* Hunter, J. D. 2007, CSE, 9, 90

* Hutter, J., Alavi, A., Deutsch, T., et al. 2008, CPMD: Car-Parinello Molecular Dynamics, v3.17.1, IBM Corp 1990–2008 and MPI für Festkörperforschung Stuttgart 1997–2001. www.cpmd.org

* Ichimaru, S., Iyetomi, H., & Tanaka, S. 1987, PhRv, 149, 91

* Iglesias, C. A., Rogers, F. J., & Saumon, D. 2002, ApJL, 569, L111

* John, T. 1994, MNRAS, 269, 871

* Kendall, R. A., Dunning, T. H., & Harrison, R. J. 1992, JChPh, 96, 6796

* Kiselyov, O. E., & Martynov, G. A. 1990, JChPh, 93, 1942

* Koester, D., & Wolff, B. 2000, A&A, 357, 587

* Kollmar, C., & Neese, F. 2010, MolPh, 108, 2449
Kolos, W., & Peek, J. 1976, CP, 12, 381

Kowalski, P. 2006a, PhD thesis, Vanderbilt Univ.

Kowalski, P., Mazevet, S., Saumon, D., & Challacombe, M. 2007, PhRvB, 76, 075112

Kowalski, P. M. 2006b, ApJ, 641, 488

Kowalski, P. M. 2010a, A&A, 519, L8

Kowalski, P. M. 2010b, in AIP Conf. Proc. 1273, 18th European White Dwarf Workshop, ed. K. Werner & T. Rauch (Melville, NY: AIP), 424

Kowalski, P. M. 2014, A&A, 566, L8

Kramida, A., Ralchenko, Y., Reader, J., & NIST ASD Team 2015, NIST Atomic Spectra Database, v5.5.6, National Institute of Standards and Technology, Gaithersburg, MD. http://physics.nist.gov/asd

Landau, L., & Lifshitz, E. 1980, Course of Theoretical Physics, Vol. 5 (3rd ed.; Oxford: Pergamon)

Lodders, K. 2003, ApJ, 591, 1220

Lovallo, C. C., & Klobukowski, M. 2004, JChPh, 120, 246

Magni, G., & Mazzitelli, I. 1979, A&A, 72, 134

Maroulis, G. 2000, JPCA, 104, 4772

Martyna, G. J., & Tuckerman, M. E. 1999, JChPh, 110, 2810

Marx, D., & Hutter, J. 2000, Modern Methods and Algorithms of Quantum Chemistry, 1, 141

Masili, M., & Starace, A. F. 2003, PhRvA, 68, 012508

Mihalas, D., Dappen, W., & Hummer, D. 1988, ApJ, 331, 815

Neese, F. 2012, WIREs Compt. Mol. Sci., 2, 73

Neese, F., Hansen, A., Wennmohs, F., & Grimme, S. 2009, AcChR, 42, 641

Nosé, S. 1984, JChPh, 81, 511

Pachucki, K. 2012, PhRvA, 85, 042511

Partridge, H., Stallcop, J. R., & Levin, E. 2001, JChPh, 115, 6471

Percus, J. K., & Yevick, G. J. 1958, PhRv, 110, 1

Perdew, J. P., Burke, K., & Ernzerhof, M. 1996, PhRvL, 77, 3865

Plimpton, S. 1995, JCoPh, 117, 1

Prusti, T., De Bruijne, J., Brown, A. G., et al. 2016, A&A, 595, A1

Raghavachari, K., Trucks, G. W., Pople, J. A., & Head-Gordon, M. 1989, CPL, 157, 479

Rohrmann, R. D. 2018, MNRAS, 473, 457

Ross, M., & Young, D. A. 1986, PLA, 118, 463

Saumon, D., & Chabrier, G. 1991, PhRvA, 44, 5122

Saumon, D., & Chabrier, G. 1992, PhRvA, 46, 2084

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Vanderbilt, D. 1990, PhRvB, 41, 7892

Walkup, R., Stewart, B., & Pritchard, D. 1984, PhRvA, 29, 169

Wigner, E. 1932, PhRv, 40, 749

Wilson, D. J., Gänsicke, B. T., Koester, D., et al. 2015, MNRAS, 451, 3237

Winisdoerffer, C., & Chabrier, G. 2005, PhRvE, 71, 026402

Wolff, B., Koester, D., & Liebert, J. 2002, A&A, 385, 995

Woon, D. E., & Dunning, T. H. 1993, JChPh, 98, 1358

Zaghloul, M. R. 2009, ApJ, 699, 885

Zeidler-KT, E.-M., Weidemann, V., & Koester, D. 1986, A&A, 155, 356