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

# A New Generation of Cool White Dwarf Atmosphere Models. IV. Revisiting the Spectral Evolution of Cool White Dwarfs

S. Blouin¹, P. Dufour¹, C. Thibeault¹, 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, 98 bis boulevard Arago, F-75014 Paris, France

Received 2019 March 21; revised 2019 May 3; accepted 2019 May 5; published 2019 June 14

## Abstract

As a result of competing physical mechanisms, the atmospheric composition of white dwarfs changes throughout their evolution, a process known as spectral evolution. Because of the ambiguity of their atmospheric compositions and the difficulties inherent to the modeling of their dense atmospheres, no consensus exists regarding the spectral evolution of cool white dwarfs ($T_{eff} < 6000$ K). In the previous papers of this series, we presented and observationally validated a new generation of cool white dwarf atmosphere models that include all the necessary constitutive physics to accurately model those objects. Using these new models and a homogeneous sample of 501 cool white dwarfs, we revisit the spectral evolution of cool white dwarfs. Our sample includes all spectroscopically identified white dwarfs cooler than 8300 K for which a parallax is available in Gaia DR2 and photometric observations are available in Pan-STARRS1 and 2MASS. Except for a few cool carbon-polluted objects, our models allow an excellent fit to the spectroscopic and photometric observations of all objects included in our sample. We identify a decrease of the ratio of hydrogen- to helium-rich objects between 7500 and 6250 K, which we interpret as the signature of convective mixing. After this decrease, hydrogen-rich objects become more abundant up to 5000 K. This puzzling increase, reminiscent of the non-DA gap, has yet to be explained. At lower temperatures, below 5000 K, hydrogen-rich white dwarfs become rarer, which rules out the scenario in which the accretion of hydrogen from the interstellar medium dominates the spectral evolution of cool white dwarfs.

**Key words:** opacity – stars: atmospheres – stars: evolution – white dwarfs

**Supporting material:** figure sets, machine-readable tables

## 1. Introduction

White dwarf stars are compact stellar remnants condemned to a slow cooling that extends over billions of years. One of the fundamental properties of white dwarfs is their intense surface gravity ($\log g = 8$ on average). This strong gravitational field implies that heavy elements are expected to sink toward the core of the white dwarf and that light elements should float high up in the atmosphere (Schatzman 1945). If gravitational settling were the only process influencing the chemical composition of white dwarf atmospheres, then all white dwarfs would show pure hydrogen atmospheres (i.e., only the DA and DC spectral types would exist).

However, we know that at least ~20% of white dwarfs have helium-rich atmospheres (Giammichele et al. 2012; Kepler et al. 2015), that many have metal-polluted atmospheres (DQs and DZs), and that some even have carbon-dominated (hot DQs; Dufour et al. 2007b) or oxygen-dominated atmospheres (DS;<fnref n="4" /> Kepler et al. 2016a). These simple observational facts imply that additional mechanisms can alter the atmospheric composition of a white dwarf. For instance, a white dwarf atmosphere can undergo simple diffusion (Muchmore 1984; Paquette et al. 1986; Koester 2009), convective dilution (Fontaine & Wesemael 1987; MacDonald & Vennes 1991; Rolland et al. 2018), convective mixing (Koester 1976; Vauclair & Reisse 1977; Dantona & Mazzitelli 1979; Rolland et al. 2018), convective dredge-up from the core (Pelletier et al. 1986), radiative levitation (Chayer et al. 1995), accretion from

the interstellar medium (Dupuis et al. 1993; Koester & Wilken 2006), and accretion of rocky planetesimals (Graham et al. 1990; Jura 2003; Farihi et al. 2010). As those mechanisms change in importance during the evolution of white dwarfs, their surface compositions also change as a function of decreasing effective temperature. This phenomenon, known as spectral evolution, is one of the most studied aspects of white dwarf evolution.

The most spectacular feature of the spectral evolution of white dwarfs is probably the existence of a so-called DB gap between $T_{eff} \approx 45,000$ and 30,000 K (Liebert 1986; Liebert et al. 1986), where helium-rich white dwarfs are significantly less abundant than at higher or lower temperatures (Eisenstein et al. 2006b). This deficiency has been explained as the result of a hydrogen float-up process at the blue edge of the gap and a convective dilution process below the red edge (Fontaine & Wesemael 1987; MacDonald & Vennes 1991).

Going down the temperature scale, the next striking feature in the spectral evolution of white dwarfs is the increase of the ratio of non-DA to DA objects below $T_{eff} \approx 12,000$ K (Sion 1984; Tremblay & Bergeron 2008). This observation is the consequence of the convective mixing of a thin hydrogen layer with the thick helium envelope underneath that can transform a hydrogen-rich atmosphere into a helium-rich atmosphere (Koester 1976; Vauclair & Reisse 1977; Dantona & Mazzitelli 1979; Rolland et al. 2018). The lower the temperature of the star, the deeper the convection zone extends into the white dwarf. Therefore, objects with a thicker hydrogen layer turn into helium-rich stars later in their evolution (Tassoul et al. 1990; Bergeron et al. 1997).

<footnote n="4">⁴ DS is the proposed classification for oxygen-rich objects (see also Table 8 of Leggett et al. 2018).</footnote>
At lower temperatures, however, our understanding of the spectral evolution of white dwarfs is more limited. In fact, results obtained in the last two decades point in conflicting directions. On the one hand, analyses performed by Bergeron et al. (1997, 2001) and Kilic et al. (2006a, 2010) have revealed that the ratio of hydrogen-rich to helium-rich stars is greatly enhanced between 5000 and 6000 K. Due to its deficiency of helium-rich objects, this temperature range was termed the non-DA gap. The existence of this gap is usually understood as the consequence of processes that transform a significant fraction of helium-rich stars into hydrogen-rich stars (near 6000 K) and that later transform hydrogen-rich objects into helium-rich objects (near 5000 K). So far, no consistent physical explanation for the existence of the non-DA gap has been proposed (Hansen 1999; Malo et al. 1999; Bergeron et al. 2001).

On the other hand, results obtained by Kowalski & Saumon (2006) and Kilic et al. (2009a, 2009b) suggest that the overabundance of hydrogen-rich objects is not limited to the 5000–6000 K range and that it extends all the way to the coolest objects at $T_{\text{eff}} \approx 4000$ K. In fact, they find that almost every single object in their samples is hydrogen-rich for $T_{\text{eff}} < 5000$ K. The overabundance of hydrogen-rich stars between 5000 and 6000 K would therefore not be a gap but rather a portion of a continuous process that irreversibly transforms helium-rich atmosphere into hydrogen-rich atmospheres. Kowalski (2006) has tentatively suggested that this transformation could be due to the accretion of hydrogen from the interstellar medium.

The fact that the atmospheric composition of the coolest white dwarfs is still open to debate is concerning for white dwarf age dating. As white dwarfs cool down monotonically and as their cooling rate can be precisely modeled (Hansen 1999; Fontaine et al. 2001; Renedo et al. 2010), they can be very accurate cosmic clocks (for applications to various stellar populations, see Oswalt et al. 1996; García-Berro et al. 2010; Kalirai 2012; Tremblay et al. 2014). However, in order to use a white dwarf as a cosmochronometer, a precise determination of its atmospheric parameters is required. In particular, it is important to know its atmospheric composition, because helium-rich and hydrogen-rich atmospheres do not have the same opacities and therefore have different cooling rates. For instance, for a white dwarf with $T_{\text{eff}} = 4000$ K and $\log g = 8$, a mistake on the atmospheric composition can lead to an error of $\sim 1$ Gyr on its cooling age (Fontaine et al. 2001, assuming that the white dwarf evolves with a constant atmospheric composition). In this context, knowing the atmospheric composition of the coolest (and thus oldest) white dwarfs becomes an even more pressing problem.

There are two main reasons that explain the discrepancies between studies that conclude the existence of the non-DA gap and those that suggest a continuous increase of the hydrogen-rich to helium-rich ratio. The first one is that the samples of Kowalski & Saumon (2006) and Kilic et al. (2009a, 2009b) are limited to DA and DC white dwarfs. Therefore, they completely ignore the existence of DQs and DZs, which have helium-rich atmospheres (for limits on the hydrogen abundance in the atmospheres of DQs and DZs, see Dufour et al. 2005, 2007a). For this reason, their samples are strongly biased toward hydrogen-rich stars. In particular, we can be confident that helium-rich white dwarfs do exist below 5000 K, because DQs and DZs are found at those temperatures.

The second reason that explains the discrepancy between both sets of studies is related to the fact that both helium-rich

and hydrogen-rich stars become DCs below $\approx 5000$ K. In fact, the thermal energy becomes too small to excite the atomic states that are required to produce hydrogen or helium spectral lines in the visible, and most white dwarfs with $T_{\text{eff}} \lesssim 5000$ K thus show a flat, featureless spectrum. This implies that the chemical composition of the atmospheres of such stars must be derived solely from photometric observations (as no useful information can be retrieved from the spectroscopy). Atmosphere models are used to fit the spectral energy distribution (SED) assuming both hydrogen-rich and helium-rich models, and the best solution indicates the most likely composition of the atmosphere. As this process depends on a detailed fit of the SED, small differences between different sets of atmosphere models can lead to different solutions. The same star can be classified as having a hydrogen-rich atmosphere using one set of models and as having a helium-rich atmosphere using another set of models. This is precisely what happens here. Compared to the models of Bergeron et al. (1995, used in Bergeron et al. 1997, 2001; Kilic et al. 2006a, 2010), the models employed in Kowalski & Saumon (2006) and Kilic et al. (2009a, 2009b) include numerous improvements for the treatment of nonideal effects in the dense atmospheres of cool white dwarfs. These improved models predict that cool helium-rich white dwarfs have SEDs that are very close to those of blackbodies. As this is not observed, cool white dwarfs are almost all classified as hydrogen-rich (Kowalski 2006).

This paper aims to fix the issues that tarnished previous analyses of the spectral evolution of cool white dwarfs. On the theoretical front, our analysis is based on state-of-the-art atmosphere models that take into account all high-density effects relevant for the modeling of the atmospheres of cool white dwarfs (Blouin et al. 2018a, hereafter Paper I). These new models have the advantage of also being applicable to metal-polluted white dwarfs. Not only is this useful for modeling all cool white dwarfs (and not only pure hydrogen and pure helium objects), but it is also an excellent way of observationally validating the input physics of our models. Indeed, as metal-polluted white dwarfs (DQs and DZs) are the only ones to still show spectral features in the visible below $T_{\text{eff}} \approx 5000$ K, they represent a unique opportunity to test our models against observations. In contrast, as there is no challenge in fitting the flat spectrum of DC white dwarfs, they cannot be used to discriminate between good and poor models. Using this insight, we tested our models against some of the most challenging cool metal-polluted white dwarfs (see Paper I and Blouin et al. 2018b, 2019b, hereafter Papers II and III). In every case, we showed that our models can accurately reproduce the complex spectral features observed in those objects, which validates their constitutive physics.

On the observational front, the analysis presented in this paper is based on the largest (501 objects) homogeneous sample of cool white dwarfs ever studied. The homogeneity of our sample is guaranteed by the fact that all our fits are based on data obtained from the same three surveys (Pan-STARRS1, 2MASS, and Gaia DR2). Additionally, no discrimination on the spectral types was made during the selection of the sample, which ensures a good representation of the diversity of the chemical compositions of cool white dwarf atmospheres.

Section 2 details our methodology, including the selection of the sample, the model atmosphere code and the fitting procedures. In Section 3, we present the solutions of our fits and analyze in detail some interesting objects. In particular, we turn our attention
to the so-called peculiar non-DAs—stars for which the photo-
metric and spectroscopic observations suggest conflicting atmo-
spheric compositions (Bergeron et al. 1997)—and to DQpec white
dwarfs—carbon-polluted objects with distorted Swan bands (Hall
& Maxwell 2008). After a careful analysis of the biases that affect
our sample, Section 4 presents our conclusions on the spectral
evolution of cool white dwarfs. Finally, our main findings are
summarized in Section 5.

**2. Methodology**

2.1. Sample Selection

We defined the selection criteria of our sample so that
the maximum number of white dwarfs could be included, while
ensuring that the observations used for the photometric fits
are homogeneous. Six criteria must be satisfied for an object to
be part of our sample.

1. It must have a parallax measurement from the Gaia DR2 (Gaia Collaboration 2016, 2018). The parallax is important not only for measuring the mass of the white dwarf, but also for identifying unresolved binary systems.

2. Each object must have *grizy* photometry from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1; Chambers et al. 2016).

3. For each object, there must be at least *J* photometry (ideally *H* and *K* also) in the Two Micron All Sky Survey (2MASS).

4. A spectrum must be available for each object, so that the appropriate atmospheric composition can be assumed in our models. In particular, the presence of $H\alpha$ is useful for identifying hydrogen-rich atmospheres and metal lines, and $C_2$ Swan bands are used to fix the amount of heavy elements in the model atmospheres of DZ and DQ white dwarfs, respectively.

5. Stars that are part of a known unresolved binary system were rejected.

6. Only objects cooler than 8300 K were retained. This upper limit is motivated by our wish to verify the existence of the non-DA gap. To do so, we need to compare the abundance of hydrogen-rich stars above and below the blue edge of this gap at $T_{\text{eff}} \approx 6000$ K. The 8300 K value was chosen as it allows to compute the fraction of hydrogen-rich stars as a function of effective temperature for 500 K bins up to a bin centered at 8000 K.

We relied on the Montreal White Dwarf Database (Dufour et al. 2017) to identify objects that match these criteria. In total, we found 503 stars that satisfy all six criteria. The bulk of these objects are part of the Limoges et al. (2015) sample (292 objects) and the rest comes from a number of other studies (Bergeron et al. 1997, 2001, 2005; Putney 1997; Kawka & Vennes 2006, 2012; Subasavage et al. 2007, 2008, 2009, 2017; Kilic et al. 2010; Giammichele et al. 2012; Sayres et al. 2012; Kleinman et al. 2013; Gianninas et al. 2015; Kepler et al. 2015, 2016b). Among the 503 objects initially selected, two of them (GJ 1221 and GJ 1228) were rejected because no appropriate models are available at the moment to properly model their atmospheres. These two white dwarfs, classified as DXP stars, are characterized by broad, unidentified absorption features and a very large magnetic field ($B > 100$ MG; Putney & Jordan 1995; Berdyugin & Piirola 1999). Therefore, our final

sample contains 501 objects, which are listed in Tables 1
(astrometric data) and 2 (photometric data).

Due to the relatively small limiting magnitude of 2MASS,
requiring each object to be in 2MASS considerably reduces the
size of our sample. Nevertheless, this constraint is important
because observations in the infrared are necessary to detect
collision-induced absorption (CIA) from molecular hydrogen,
which is very useful to constrain the hydrogen abundance. This
is a particularly important parameter to determine if we want to
get an accurate picture of the spectral evolution of cool DC
white dwarfs.

2.2. Atmosphere Models

Our model atmosphere code is described at length in Paper I.
It is based on the code described in Bergeron et al. (1995) and
Dufour et al. (2005, 2007a), but also includes a number of
nonideal high-density effects that arise in the dense atmo-
spheres of cool white dwarfs. In particular, continuum opacities
are corrected for collective interactions (Iglesias et al. 2002;
Rohrmann 2018), an ab initio equation of state for hydrogen
and helium is assumed (Becker et al. 2014), the pressure
ionization of helium is modeled using the ab initio calculations
of Kowalski et al. (2007), an accurate description of the
pressure broadening of $Ly\alpha$ is included (Kowalski &
Saumon 2006), the $H_2$–He CIA profiles are corrected for
many-body collisions (Blouin et al. 2017), and CIA from He–
He–He interactions is included (Kowalski 2014). Moreover,
our code includes significant improvements relevant for the
accurate modeling of cool metal-polluted atmospheres, such as
accurate line profiles for important heavy element lines
(Paper III; Allard & Alekseev 2014; Allard et al. 2014,
2016a, 2016b, 2018; Blouin et al. 2019a), a state-of-the-art
treatment of the nonideal ionization equilibrium of C, Ca, Fe,
Mg, and Na (Paper I), and the distortion of the $C_2$ Swan bands
in DQpec stars (Kowalski 2010). Note also that we now
include the $C_2$ Swan band opacity using a line-by-line approach
with the line list of Hornkohl et al. (2005). This approach
noticeably improves our fit to the shape of the Swan bands
compared to the just overlapping line approximation (Zeidler-
K.T. & Koester 1982) implemented in the models of Dufour
et al. (2005).

2.3. Fitting Procedures

2.3.1. DAs

The fundamental parameters of DA white dwarfs are obtained with the photometric technique (Bergeron et al. 1997). The solid angle $\pi(R/D)^2$ and $T_{\text{eff}}$ are found by fitting the model fluxes to the observed SED with the Levenberg-Marquardt algorithm. Because $D$ is known from the *Gaia* parallax measurement, the white dwarf radius $R$ can be computed from the solid angle. The mass and the surface gravity of the white dwarf are then found with the evolutionary models of Fontaine et al. (2001). Note that for all DA stars in our sample we assume a pure hydrogen atmosphere. This assumption is very common in the literature (e.g., Bergeron et al. 2001; Gianninas et al. 2011; Limoges et al. 2015) and is supported by the comparison between synthetic and observed Balmer lines profiles. For cool white dwarfs, Balmer lines become shallower and wider when the H/He ratio is decreased. While it is not possible to distinguish between Balmer lines created in a pure hydrogen atmosphere and an atmosphere with $\log H/He = 2$,
<table><caption>**Table 1**

Astrometry of Objects Included in Our Sample</caption><thead><tr><th>PSO<sup>a</sup></th><th>MWDD Object ID</th><th>Gaia Object ID</th><th>R.A. (J2015.5)</th><th>Decl. (J2015.5)</th><th>π (mas)</th></tr></thead><tbody><tr><td>J000224.474+635745.512</td><td>2MASS J00022257+6357443</td><td>431635455820288128</td><td>0.6033724</td><td>63.9627281</td><td>38.080(0.079)</td></tr><tr><td>J000410.478-034008.751</td><td>PHL 2595</td><td>2447889401738675072</td><td>1.0439384</td><td>-3.6691769</td><td>21.076(0.100)</td></tr><tr><td>J000720.831+123018.402</td><td>WD 0004+122</td><td>2766234439302571904</td><td>1.8371872</td><td>12.5048777</td><td>57.308(0.114)</td></tr><tr><td>J000728.917+340342.227</td><td>NLTT 301</td><td>2875903332533220992</td><td>1.8705266</td><td>34.0618834</td><td>29.336(0.099)</td></tr><tr><td>J000754.487+394731.294</td><td>LP 240-30</td><td>383108338321272448</td><td>1.9774884</td><td>39.7919426</td><td>29.038(0.064)</td></tr><tr><td>J000935.139+310840.332</td><td>EGGR 1</td><td>2861792754354276352</td><td>2.3959624</td><td>31.1440912</td><td>18.638(0.077)</td></tr><tr><td>J001122.399+424038.082</td><td>GD 5</td><td>384636109728592768</td><td>2.8432990</td><td>42.6770566</td><td>42.714(0.045)</td></tr><tr><td>J001214.188+502514.381</td><td>GJ 1004</td><td>395234439752169344</td><td>3.0583199</td><td>50.4200515</td><td>91.983(0.029)</td></tr><tr><td>J001412.352-131109.057</td><td>GJ 3016</td><td>2418116963320446720</td><td>3.5508822</td><td>-13.1867135</td><td>53.812(0.070)</td></tr><tr><td>J001737.914-051650.624</td><td>LP 644-81</td><td>2443826805441462656</td><td>4.4082578</td><td>-5.28121904</td><td>21.405(0.153)</td></tr></tbody></table>

**Note.**

<sup>a</sup> Pan-STARRS object name.

(This table is available in its entirety in machine-readable form.)

<table><caption>**Table 2**

Spectral Type and Photometric Data</caption><thead><tr><th>PSO</th><th>Spectral Type</th><th>g</th><th>r</th><th>i</th><th>z</th><th>y</th><th>J</th><th>H</th><th>K</th></tr></thead><tbody><tr><td>J000224.474+635745.512</td><td>DC</td><td>17.63</td><td>16.99</td><td>16.73</td><td>16.63</td><td>16.57</td><td>15.80</td><td>15.57</td><td>15.51</td></tr><tr><td>J000410.478-034008.751</td><td>DA</td><td>16.91</td><td>16.75</td><td>16.73</td><td>16.76</td><td>16.79</td><td>16.12</td><td>15.89</td><td>15.44</td></tr><tr><td>J000720.831+123018.402</td><td>DC</td><td>16.77</td><td>16.25</td><td>16.02</td><td>15.93</td><td>15.92</td><td>15.08</td><td>15.10</td><td>14.90</td></tr><tr><td>J000728.917+340342.227</td><td>DC</td><td>17.61</td><td>17.23</td><td>17.12</td><td>17.05</td><td>17.08</td><td>16.39</td><td>16.25</td><td>15.83</td></tr><tr><td>J000754.487+394731.294</td><td>DC</td><td>17.25</td><td>16.53</td><td>16.23</td><td>16.12</td><td>16.06</td><td>15.18</td><td>14.85</td><td>14.65</td></tr><tr><td>J000935.139+310840.332</td><td>DC</td><td>16.81</td><td>16.74</td><td>16.85</td><td>16.94</td><td>17.01</td><td>16.45</td><td>16.19</td><td>...</td></tr><tr><td>J001122.399+424038.082</td><td>DA</td><td>15.42</td><td>15.27</td><td>15.23</td><td>15.28</td><td>15.31</td><td>14.54</td><td>14.35</td><td>14.39</td></tr><tr><td>J001214.188+502514.381</td><td>DAH</td><td>14.52</td><td>14.27</td><td>14.20</td><td>14.21</td><td>14.22</td><td>13.49</td><td>13.25</td><td>13.19</td></tr><tr><td>J001412.352-131109.057</td><td>DAH</td><td>16.12</td><td>15.76</td><td>15.62</td><td>15.59</td><td>15.60</td><td>14.81</td><td>14.55</td><td>14.63</td></tr><tr><td>J001737.914-051650.624</td><td>DA</td><td>17.92</td><td>17.52</td><td>17.35</td><td>17.32</td><td>17.27</td><td>16.48</td><td>16.35</td><td>...</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

the differences between a pure hydrogen atmosphere and a helium-dominated atmosphere (log H/He < 0) are obvious. As our sample is largely based on the sample of Limoges et al. (2015) and as their spectroscopic analysis shows no evidence of helium-dominated DAs, we conclude that we can safely assume that all DAs in our sample are hydrogen-rich. The question of whether they simply have a hydrogen-dominated atmosphere or a pure hydrogen atmosphere remains, but as we cannot distinguish between both possibilities, we assume a pure hydrogen composition for simplicity. Note that this assumption does not interfere with the main purpose of this work, because we only need to know which stars are hydrogen-rich and which ones are helium-rich.

### 2.3.2. DCs

The fitting procedure for DC stars is identical to the procedure described above for DAs, except that the atmospheric composition is a priori unknown. If the effective temperature is high enough that we should see Balmer lines if the atmosphere were hydrogen-rich, then we assume a helium-rich composition. If the temperature is too low to produce detectable Balmer lines ($T_{eff} \lesssim 5000$ K), then the atmospheric composition determination is based on the photometry. Finding the atmospheric composition of a DC white dwarf solely from the photometry can be difficult. That being said, differences

between SEDs produced by stars with different compositions (see Figure 1) can generally be exploited to confidently establish the atmospheric composition of such objects.

Three distinct photometric fits are performed for each cool DC white dwarf in our sample: one assuming a pure hydrogen atmosphere, one assuming a pure helium atmosphere and one assuming a mixed H/He atmosphere, where the H/He ratio is a free parameter that is adjusted to the observations. To identify the best of the three solutions, we do not simply pick the solution with the smallest $\chi^2$. This approach is misguided when comparing solutions obtained with different numbers of free parameters. In fact, it is almost always possible to obtain a smaller $\chi^2$ by adding one more free parameter to the fit (the H/He ratio), but the gain obtained with this additional parameter is not always statistically significant. Instead, we use the Akaike information criterion (AIC). This criterion penalizes a model according to the number of free parameters it contains, thus enforcing Occam's razor. The best model is the one with the smallest AIC, which is given by (Akaike 1974; Burnham & Anderson 2002)

$$AIC = 2k - 2 \ln L, \qquad (1)$$

where k is the number of degrees of freedom in the model and L is the maximum value of the likelihood function. Note that for least-squares model fitting the likelihood function is simply given by $\ln L = C - \chi^2/2$, where C is a constant that is independent from the model and can be ignored. When we fit
<figure><img src="image_2.png" /><figcaption>Figure 1. Synthetic Pan-STARRS and 2MASS photometry obtained with our model atmosphere code. Each panel corresponds to a different effective temperature, and each symbol represents a different atmospheric composition (see legend). In all cases, a surface gravity of log g = 8 is assumed.</figcaption></figure>

photometric data, we are only fitting a small number of data points. In such circumstances, it is necessary to use an AIC that is corrected for small data samples (Hurvich & Tsai 1989),

$$AICc = AIC + \frac{2k(k + 1)}{n - k - 1}, \quad (2)$$

where n is the sample size (in our case the sample size corresponds to the number of bands included in our photo-metric fit).

To our knowledge, this is the first time that such a detailed statistical analysis is applied to the problem of the determination of the atmospheric composition of DC stars. Previous analyses either only used mixed H/He models when CIA was obviously present in the SED (Bergeron et al. 1997, 2001; Kilic et al. 2010) or directly compared mixed and pure solutions without rigorously computing the significance of the χ² reduction following the addition of a third free parameter (Kowalski 2006; Kilic et al. 2009a).<fnref n="5" /> Both approaches are problematic. The first one could miss mixed objects that do not

show strong CIA, and the second one is biased toward mixed solutions (e.g., ≈50% of DCs analyzed by Kilic et al. 2009a are classified as mixed objects). By accurately quantifying the weight of evidence in favor of each solution, our methodology should not be influenced by those biases.

We performed a small numerical experiment that clearly illustrates the danger of using an approach where the H/He ratio is directly considered as a free parameter (i.e., by fitting every DC star with a grid of mixed H/He compositions ranging from pure hydrogen to pure helium and adjusting the H/He ratio as a free parameter). We generated a set of 100 “synthetic” pure hydrogen stars with effective temperatures ranging from 4000 to 6000 K. To do so, we used the synthetic photometry obtained from our models to which we added a Gaussian noise to mimic the observational uncertainties (we assumed a typical 3% uncertainty for the visible photometry and 5% for the infrared). For 45% of these pure hydrogen objects, the fit obtained by adjusting the H/He ratio in order to obtain the smallest χ² corresponds to a mixed H/He solution (i.e., -5 < log H/He < 0.5).<fnref n="6" /> This numerical experiment demonstrates that this approach is

<footnote n="5">⁵ A third approach, proposed by Saumon et al. (2014), consists of plotting χ² as a function of the H/He abundance ratio and checking if χ²(H/He) has a significant minimum. Note, however, that no quantitative criterion was used to define what makes a χ²(H/He) minimum statistically significant.</footnote>

<footnote n="6">⁶ Note that pure helium models are less susceptible to this bias, as the same experiment with pure helium models shows that 10% are (incorrectly) classified as mixed objects.</footnote>
<table><caption>**Table 3**
Atmospheric Parameters</caption><thead><tr><th>PSO</th><th><i>T</i><sub>eff</sub> (K)</th><th>log <i>g</i></th><th><i>M</i> (<i>M</i><sub>☉</sub>)</th><th>τ<sub>cool</sub> (Gyr)<sup>a</sup></th><th>log H/He</th><th>log C/He</th><th>log Ca/He</th></tr></thead><tbody><tr><td>J000224.474+635745.512</td><td>4565(50)</td><td>7.758(0.024)</td><td>0.432(0.013)</td><td>4.45</td><td>He</td><td>...</td><td>...</td></tr><tr><td>J000410.478-034008.751</td><td>6955(60)</td><td>7.931(0.023)</td><td>0.551(0.013)</td><td>1.40</td><td>H</td><td>...</td><td>...</td></tr><tr><td>J000720.831+123018.402</td><td>4885(45)</td><td>8.090(0.021)</td><td>0.625(0.014)</td><td>6.45</td><td>He</td><td>...</td><td>...</td></tr><tr><td>J000728.917+340342.227</td><td>5545(50)</td><td>8.153(0.021)</td><td>0.666(0.013)</td><td>5.02</td><td>He</td><td>...</td><td>...</td></tr><tr><td>J000754.487+394731.294</td><td>4680(25)</td><td>6.768(0.039)</td><td>0.081(0.016)</td><td>1.40</td><td>H</td><td>...</td><td>...</td></tr><tr><td>J000935.139+310840.332</td><td>7960(290)</td><td>8.053(0.100)</td><td>0.609(0.062)</td><td>1.21</td><td>He</td><td>...</td><td>...</td></tr><tr><td>J001122.399+424038.082</td><td>6990(50)</td><td>7.975(0.017)</td><td>0.577(0.010)</td><td>1.47</td><td>H</td><td>...</td><td>...</td></tr><tr><td>J001214.188+502514.381</td><td>6445(40)</td><td>8.247(0.012)</td><td>0.746(0.008)</td><td>3.24</td><td>H</td><td>...</td><td>...</td></tr><tr><td>J001412.352-131109.057</td><td>5855(35)</td><td>8.217(0.016)</td><td>0.724(0.010)</td><td>4.10</td><td>H</td><td>...</td><td>...</td></tr><tr><td>J001737.914-051650.624</td><td>5630(40)</td><td>7.918(0.026)</td><td>0.537(0.015)</td><td>2.40</td><td>H</td><td>...</td><td>...</td></tr></tbody></table>

**Note.**

a Cooling times are computed using the evolutionary models of Fontaine et al. (2001).

(This table is available in its entirety in machine-readable form.)

heavily biased toward finding mixed solutions, even for stars that actually have a pure composition. In contrast, if we use our AIC approach, we find the correct (pure hydrogen) solution in all cases.

We also performed an experiment to make sure that our AIC approach is not biased toward pure compositions. To do so, we fitted a set of 100 “synthetic” mixed H/He white dwarfs with effective temperatures ranging from 4000 to 6000 K and hydrogen abundances between log H/He = −5 and 0. As in our previous experiment, we injected a Gaussian noise to mimic the observational uncertainties. We find that our AIC approach retrieves a mixed solution for 92% of these objects.<fnref n="7" /> The cases where our AIC approach fails to retrieve a mixed solution are all for objects with hydrogen abundances near log H/He = 0, where the SED becomes very similar to that of a pure hydrogen object. Therefore, the bias of the AIC approach against mixed objects is very weak compared to the bias against pure compositions exhibited by the method where the H/He ratio is directly considered as a free parameter, which justifies our methodology.

The AIC can also be used to estimate how likely it is that our atmospheric composition determinations are accurate (Burnham & Anderson 2002). For each model *i* (i.e., for the pure hydrogen, the pure helium, and the mixed H/He models) we can compute the difference in AIC with respect to the best model,

$$ \Delta_i(\text{AIC}) = \text{AIC}_i - \min \text{ AIC}. \quad (3) $$

From there, the Akaike weights $w_i$ are computed as

$$ w_i = \frac{\exp\left[-\frac{1}{2}\Delta_i(\text{AIC})\right]}{\sum_j \exp\left[-\frac{1}{2}\Delta_j(\text{AIC})\right]}. \quad (4) $$

Each weight $w_i$ can be interpreted as the probability that model *i* is the best model. Therefore, we can now quantify the degree of confidence of our atmospheric composition determinations. In particular, given our hydrogen-rich and helium-rich fits, we can compute the odds that a DC star is hydrogen-rich versus

helium-rich. Some examples of the usefulness of the Akaike weights are given in Section 3.1.

### 2.3.3. DQs and DZs

For metal-polluted objects, we still rely on a photometric fit (see Section 2.3.1) to determine $T_{\text{eff}}$ and log g, but we also use the observed spectrum to find the abundances of heavy elements. Our approach is identical to that of Dufour et al. (2005, 2007a). Once a photometric solution is found, we adjust the metal abundances to obtain a good fit to the spectroscopy. As the abundances derived from this spectroscopic fit are usually different from our initial guess, we repeat the whole procedure (including the photometric fit) until a consistent solution is found. For DQ and DQpec objects, we use the C₂ Swan bands to constrain the C/He abundance ratio. For DZs, we use the Ca II H & K doublet and the Ca I resonance line to find Ca/He. The Fe/He, Mg/He, and Na/He ratios are also adjusted when the relevant spectral lines are visible. Otherwise, the abundance ratios of all heavy elements are scaled to the abundance of Ca to match the abundance ratios of CI chondrites (Lodders 2003).<fnref n="8" /> Finally, unless there is direct evidence of the presence of hydrogen in the atmosphere (i.e., Balmer lines, CH bands, or H₂-He CIA), a hydrogen-free atmosphere is assumed for all DQ and DZ white dwarfs.

## 3. Results

Following the fitting procedures described in Section 2.3, we obtained the atmospheric parameters of all 501 objects included in our sample. They are listed in Table 3, which is published in its entirety in machine-readable format. Our statistical analysis of this sample and its implication for the spectral evolution of cool white dwarfs is presented in Section 4. The rest of this section is devoted to a more in-depth analysis of a number of noteworthy objects.

### 3.1. DAs and DCs

Figure 2 displays some examples of photometric fits to DA and DC white dwarfs (the complete set of photometric fits is available

<footnote n="7">⁷ While we find a mixed solution in 92% of cases, we find the correct H/He ratio in 68% of cases. This difference is due to the degeneracy between the SEDs of mixed objects with a high H/He ratio and a low H/He ratio (see Section 3.1 of Paper II).</footnote>

<footnote n="8">⁸ There is an exception to this rule. For Fe, we set the default abundance to half that predicted by the chondritic Fe/Ca number ratio. This is a rough way of taking the faster diffusion of Fe (Paquette et al. 1986; Koester 2009; Hollands et al. 2017) into account when Fe lines are not detected.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Examples of photometric fits to DC (left column) and DA white dwarfs (right column). For DCs, pure hydrogen, pure helium, and mixed solutions are shown (the best fit is indicated by the small arrow next to the legend). Note that the mixed H/He solution is only shown if it corresponds to a different solution than the pure hydrogen and pure helium solutions. The complete set of photometric fits (459 objects) is available in the online journal.</figcaption></figure>

(The complete figure set (46 images) is available).

in the figure set). For DCs, both pure hydrogen and pure helium solutions are shown. The mixed H/He solution is only shown if the solution found when the H/He ratio was adjusted as a free parameter does not correspond to a pure hydrogen or a pure helium solution. For each object, the composition that leads to the smallest Akaike weight is designated as being the most likely composition (Section 2.3.2).

Note that the Akaike weights are particularly useful when trying to assess our degree of confidence in the determination of the atmospheric composition of DC stars. For instance, for a star like WD 1033+714 (Figure 2), the Akaike weights are $w_H = 0.93$ for the pure hydrogen solution, $w_{He} = 0.07$ for the pure helium solution, and $w_{mixed} = 0.00$ for the mixed H/He solution. This indicates that we can be quite confident that it is a hydrogen-rich object. Similarly, the Akaike weights strongly suggest that WD 1346+121 has a pure helium atmosphere ($w_{He} = 0.94$; see Figure 2). Note that if we had simply used the $\chi^2$ to identify the atmospheric composition of WD 1346+121, the mixed model would have been chosen because its $\chi^2$ is smaller than that of the pure hydrogen and pure helium solutions. By adding a penalty for the additional free parameter,

the Akaike weights reverse this verdict and assign a low probability to the mixed solution. Figure 3 shows an example where large observational errors lead to a much more uncertain composition. The Akaike weights for the mixed and pure hydrogen solutions are virtually equal, which reflects our inability to confidently identify the atmospheric composition of this object.

Another interesting aspect of Figure 2 is the presence of two $T_{eff} < 5000$ K DC white dwarfs for which we find a pure helium solution. Note that the Akaike weights for the pure helium solutions are $\ge 0.9$ in both cases, which suggest that we can be quite confident in those solutions. In fact, out of the 74 DC stars cooler than 5000 K in our sample, we find pure helium solutions for 16 of them. We also find that 3 of those 74 DC stars have a mixed H/He atmosphere dominated by helium, which implies that $\approx 25\%$ of DCs cooler than 5000 K have helium-rich atmospheres. This result is in strong disagreement with the conclusions of Kowalski (2006), Kowalski & Saumon (2006), and Kilic et al. (2009a, 2009b), according to whom helium-rich DC white dwarfs cooler than 5000 K are nearly nonexistent. The implications of this finding
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Photometric fits to the DC white dwarf WD 1434+437. It is hard to tell the atmospheric composition of this object.</figcaption></figure>

for the spectral evolution of cool white dwarfs will be discussed in Section 4.

### 3.1.1. Peculiar Non-DAs

In the 6200 K $\le T_{\text{eff}} \le 7600$ K temperature range, Bergeron et al. (1997, 2001) have identified DC stars for which the photometric fit yields a pure hydrogen solution. Given their surface temperatures, those objects are expected to show $H\alpha$ in their spectra if they have a hydrogen-rich atmosphere. The pure hydrogen photometric fits are therefore in contradiction with the spectroscopic observations, which led Bergeron et al. (1997, 2001) to designate those stars as peculiar non-DAs. Furthermore, they interpreted these peculiar non-DAs as forming a distinct physical group and, because of their proximity to the blue edge of the non-DA gap, they speculated that they are objects whose atmospheres are about to become hydrogen-rich.

Our analysis also reveals the presence of 13 objects for which the pure hydrogen composition inferred from the photometry is incompatible with the absence of $H\alpha$ (see Table 4).<fnref n="9" /> Interestingly, out of the seven peculiar non-DAs in Bergeron et al. (1997, 2001) that are part of our sample, only two (PSO J091554.370+532508.067 and PSO J104153.917 +141545.969) are found to have an SED that is better represented by a hydrogen-rich atmosphere. In principle, if peculiar non-DAs did form a physical group, we should find that objects identified as peculiar non-DAs by Bergeron et al. (1997, 2001) are also flagged as peculiar non-DAs in our analysis. These conflicting results could be due to the fact that we rely on different photometric systems. In particular, we note that our 2MASS infrared photometry is of lesser quality than the JHK photometry used in Bergeron et al. (1997, 2001). To evaluate the impact of the photometric system, we performed photometric fits for the nine peculiar non-DAs identified in Bergeron et al. (1997, 2001) using their photometric data, but with our atmosphere model grid. We found that the photometric fits favor a pure hydrogen composition for five of those objects and a pure helium composition for the remaining four.

Note that we reach the same conclusions if we use the trigonometric parallaxes reported in Bergeron et al. (1997, 2001) instead of the more accurate Gaia parallaxes.

The fact that the choice of photometric system and small differences in model grids can tip the photometric solution to another atmospheric composition raises some doubt on the interpretation that peculiar non-DAs form a distinct physical group. Could the existence of those objects simply be explained by random errors in the observations? After all, the difference between a pure hydrogen and a pure helium solution is not always statistically significant (i.e., the Akaike weights can be close to each other), which inevitably leads to a number of misclassifications. Given the Akaike weights of our best solution for every DC with $T_{\text{eff}} > 5000$ K (that is, objects for which we should see $H\alpha$ if they had a hydrogen-rich atmosphere), we can actually compute how many peculiar non-DAs we expect due to the uncertainties of our photometric fits. For the 45 DCs in this temperature range, we find an average Akaike weight of 0.76 for the best solutions. This means that we expect that the composition inferred from the photometry will be incorrect for $(24 \pm 6)\%$ of those 45 DCs, which is perfectly consistent with our actual error rate in this temperature range of 27% (12 out of 45).<fnref n="10" />

While the number of peculiar non-DAs in our sample seems consistent with random errors, we note that Bergeron et al. (1997, 2001) identified systematic trends in the observed SEDs of peculiar non-DAs that supported the interpretation that they form a distinct physical group. In particular, they found that the pure helium models systematically fail to reproduce the observed *B* and *I* photometry by underestimating the *B* flux and overestimating the *I* flux. The four clearest examples of this behavior are given in Figure 17 of Bergeron et al. (1997) and Figure 11 of Bergeron et al. (2001). In Figure 4, we show our fits to those four objects. While we do observe that our pure helium models slightly underestimate the flux in the *B* band and overestimate the flux in the *I* band, the differences are smaller than $\approx 1\sigma$ (except for WD 2011+065) and might not be significant. Moreover, for WD 0000–345 and WD 1039+145, the Akaike weights of the pure helium solution are too high to confidently rule out this solution, and, in the case of WD 0423+120, we even find that the pure helium solution is a better fit to the data than the pure hydrogen solution.

<table><caption>**Table 4**
Peculiar Non-DAs in Our Sample</caption><thead><tr><td>PSO</td><td>MWDD Object ID</td><td>T<sub>eff</sub> (K)</td></tr></thead><tbody><tr><td>J021148.605+711911.081</td><td>NLTT 7194</td><td>4790</td></tr><tr><td>J023759.111+163809.334</td><td>LP 410–67</td><td>5250</td></tr><tr><td>J095120.241+190009.463</td><td>WD 0948+192</td><td>5460</td></tr><tr><td>J084901.422+443934.262</td><td>2MASS J08490170+4439355</td><td>5810</td></tr><tr><td>J234612.394+115849.053</td><td>PM J23462+1158</td><td>5935</td></tr><tr><td>J081411.267+484526.630</td><td>WD 0810+489</td><td>6515</td></tr><tr><td>J142747.971+053230.563</td><td>WD 1425+057</td><td>6725</td></tr><tr><td>J141143.197+220644.803</td><td>WD 1409+223</td><td>6785</td></tr><tr><td>J104718.292+000717.423</td><td>WD 1044+003</td><td>6965</td></tr><tr><td>J104153.917+141545.969</td><td>V* CY Leo</td><td>7140</td></tr><tr><td>J105734.525–073122.157</td><td>LAWD 34</td><td>7155</td></tr><tr><td>J091554.370+532508.067</td><td>EGGR 250</td><td>7170</td></tr><tr><td>J122619.639+183634.295</td><td>WD 1223+188</td><td>7465</td></tr></tbody></table>

<footnote n="9">⁹ Note that in all 13 cases the best mixed solution corresponds either to a pure hydrogen composition or to a mixed H/He composition that should give rise to a detectable $H\alpha$ feature.</footnote>

<footnote n="10">¹⁰ For comparison, Bergeron et al. (2001) have a 33% error rate for DCs between 5000 and 8300 K (i.e., the maximum surface temperature of stars included in our sample).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Photometric fits to four peculiar non-DAs identified in Bergeron et al. (1997, 2001). These fits were obtained using their photometric data and trigonometric parallaxes (as they did not have any parallax measurement for WD 0423+120, we assumed log g = 8 for this object).</figcaption></figure>

As an additional test, we fitted the nine peculiar non-DAs of Bergeron et al. (1997, 2001)—all stars for which they found a better match to pure hydrogen than to pure helium models—using their observations and their atmosphere models. For the pure helium solutions, we find an average Akaike weight of 0.44, meaning that the quality of the pure hydrogen and pure helium solutions is very similar and that the chemical composition determination is quite uncertain. Moreover, for three of those objects, we actually find a smaller $\chi^2$ (and a larger Akaike weight) for the pure helium solution than for the pure hydrogen solution, suggesting that they are not peculiar non-DAs after all.

Our analysis makes a strong case for a simple interpretation according to which the existence of peculiar non-DAs is the unavoidable consequence of the intrinsic uncertainty of the photometric technique. That being said, the fact that two objects—PSO J091554.370+532508.067 (WD 0912+536) and PSO J104153.917+141545.969 (WD 1039+145)—are flagged as peculiar non-DAs in both our analysis and the analysis of Bergeron et al. (1997, 2001) prevents us from ruling out the idea that some peculiar non-DAs are part of a distinct physical group. Given that the same conclusion is reached using both different models and different observations, it appears that those two objects are indeed peculiar. Interestingly, WD 0912+536 harbors a very strong magnetic field ($B \sim 100$ MG; Angel et al. 1972; Angel 1978), possibly hinting at a connection between peculiar non-DAs and magnetism (however, no circular polarization was detected for WD 1039+145; Angel et al. 1981). In the same vein, we note that Bergeron et al. (1997) have suggested the existence of a relation between magnetic fields and the chemical evolution of cool white dwarfs.

### 3.1.2. Hydrogen Traces in Helium-rich DC White Dwarfs

Based on a comprehensive analysis of white dwarfs discovered by the Gaia mission, Bergeron et al. (2019) recently

claimed that pure helium white dwarfs below $T_{\text{eff}} = 11,000$ K are extremely rare. In a nutshell, additional electron donors (metals or hydrogen) have to be included in the atmosphere models of helium-rich objects in order to obtain reasonable masses. Formally, only upper limits on the hydrogen abundance of DC stars can be determined. As the presence of hydrogen traces affects the atmospheric parameters obtained from a photometric fit, our inability to accurately measure the H/He abundance ratio of those objects implies that our $T_{\text{eff}}$ and log g determinations are uncertain.

To see by how much $T_{\text{eff}}$ and log g can be affected by trace amounts of hydrogen, we studied a subsample made up of all DC stars for which we found a pure helium composition. For each of those objects, we performed a second photometric fit assuming a hydrogen abundance of log H/He = -5. Note that this hydrogen abundance is too small to produce any detectable $H\alpha$ feature or to give rise to significant H$_2$–He CIA. As shown in Figure 5, the addition of this trace amount of hydrogen has an important effect on $T_{\text{eff}}$ and log g, particularly at high temperatures. This implies significant uncertainties on the atmospheric parameters of helium-rich DC white dwarfs. By adding free electrons to the atmosphere, small amounts of hydrogen reduce the determined effective temperature (by 400 K at worst and by 120 K on average) and the surface gravity (to keep the luminosity constant, the star has to be inflated to compensate for the decrease of $T_{\text{eff}}$). This behavior is analogous to the reduction of the photometric $T_{\text{eff}}$ and log g following the addition of carbon—another electron donor—in DQ model atmospheres (see Figure 8 of Dufour et al. 2005). To account for this uncertainty, we shifted our helium-rich solutions for DC white dwarfs halfway between the pure helium and the log H/He = -5 solutions, and we adjusted the uncertainties so that they now encompass both solutions. The $T_{\text{eff}}$, log g, M, and $\tau_{\text{cool}}$ values reported in Table 3 include these corrections.
<figure><img src="image_6.png" /><figcaption>Figure 5. Surface gravity of helium-rich DC white dwarfs in our sample as a function of their effective temperatures. The open circles correspond to our pure helium solutions, and the filled circles represent the solutions found if log H/He = -5 is assumed. Typical uncertainties on log g and T<sub>eff</sub> are shown in the top-left corner.</figcaption></figure>

## 3.2. DZs

Figure 6 shows examples of photometric and spectroscopic fits to DZ white dwarfs. The complete set of fits is available in the online journal, except for LP 658-2, Ross 640, SDSS J080440.63 +223948.6 and WD 2251-070 (PSO J055509.987-041037.276, PSO J162824.563+364623.925, PSO J080440.637+223945.828, and PSO J225355.149-064701.663), for which we directly use the values reported in Papers I and II and in Blouin et al. (2019a). As in Papers I, II, and III, our models are in excellent agreement with the observations. Our solutions are consistent across all wavelengths, and they properly reproduce the observed spectral lines. In particular, for the three examples shown in Figure 6, we obtain good fits to the resonance line of Ca I at 4226 Å and to the Ca II H & K doublet. This suggests that the pressure and temperature structure of our DZ models are accurate, as these profiles are sensitive to the physical conditions in the line-forming regions of the atmosphere (Allard & Alekseev 2014).

## 3.3. DQs

In Figure 7, we show examples of photometric and spectroscopic fits to three DQs (the complete set of fits is available in the online journal). Clearly, our best solutions are in good agreement with the observations. Things get more complicated if we look at cooler carbon-polluted atmospheres. The cooling sequence of DQs stops at $T_{eff} \approx 6000$ K (Dufour et al. 2005; Koester & Knist 2006) to give way to DQpec white dwarfs. The spectra of those objects are vaguely similar to those of DQs: they show Swan-like absorption bands that are blueshifted by a few hundred angstroms (Bergeron et al. 1994; Schmidt et al. 1995; Hall & Maxwell 2008). Using density functional theory (DFT) calculations, Kowalski (2010) has convincingly shown that the physical explanation for the DQ → DQpec transition is the pressure-driven distortion of the C₂ Swan bands. As a carbon-polluted white dwarf cools down, its atmosphere becomes denser. Near $T_{eff} = 6000$ K, the helium density becomes high enough to affect the electronic levels of the C₂ molecule, which leads to an increase of the electronic transition energy $T_e$ of the Swan bands. It is this increase that

explains the shift of the C₂ Swan bands toward lower wavelengths. According to the ab initio calculations of Kowalski (2010), the shift of $T_e$ is linear with respect to density,

$$ \Delta T_e (\text{eV}) \approx \alpha \rho (\text{g cm}^{-1}), \quad (5) $$

with $\alpha = 1.6$. However, as pointed out by Kowalski (2010), such a shift overestimates the distortion of Swan bands in the spectra of DQpec white dwarfs. There are two ways to solve this problem: either we change the slope of Equation (5), or we find a way to decrease the photospheric density of DQpec models.

Kowalski (2010) proposed that the photospheric density of DQpec stars could be lowered if hydrogen is added to their atmospheres. In fact, in their analysis of LHS 290, they show a fit where they compensate for the strong $\Delta T_e = 1.6\rho$ shift by supposing an ad hoc log H/He = -2.2 hydrogen abundance. We are very skeptical that this is the correct way of reconciling the DFT calculations with the observations, as such a high amount of H would lead to the formation of a significant quantity of CH. We added the CH rovibrational bands to our code using the Kurucz line lists,<fnref n="11" /> which rely on data from Masseron et al. (2014). Figure 8 shows synthetic spectra computed for parameters very similar to those of LHS 290 ($T_{eff} = 6000$ K and log C/He = -6.0) and with different H/He abundance ratios. Clearly, a hydrogen abundance as low as log H/He = -4.0 is sufficient to produce an unmissable CH G band near 4300 Å, and the log H/He = -2.2 value can be safely rejected. For the models of Figure 8, the pure helium and the log H/He = -4.0 models have photospheric densities of 0.20 g cm⁻³ and 0.16 g cm⁻³ (at $\tau_R = 2/3$), respectively. This slight reduction of the density is not sufficient to compensate for the too important shift implied by the $\alpha = 1.6$ parameter. Therefore, we conclude that adding ad hoc amounts of hydrogen in the atmospheres of DQpec stars is not the solution to the discrepancy between the DFT calculations and the observations.

### 3.3.1. Calibration of $\alpha$

While we rejected the idea that an undetected amount of hydrogen could bring the density down to a level where the $\alpha = 1.6$ value is compatible with the observations, it is still possible that a missing piece of input physics in our models leads to an overestimation of the photospheric density in cool carbon-polluted white dwarfs. However, we believe that it is more likely that the current implementation of the shift of C₂ Swan bands is incorrect. In particular, we still use the standard Swan band spectrum, which we merely shift at every atmospheric layer with the $\Delta T_e$ value given by Equation (5). Until the absorption of C₂ in dense helium is modeled, which is a much more computationally intensive task than computing $\Delta T_e(\rho)$, this approximation will remain unjustified.

In the meantime, for the purpose of this paper, it is sufficient to resort to a more approximative way of fitting DQpec objects, and we calibrate $\alpha$ using the observed spectra of distorted Swan bands. After testing several $\alpha$ parameters, we found that the value that leads—for most objects—to the best agreement between models and observations is $\alpha = 0.2$. This value is actually a compromise, because we found that no fixed $\alpha$

<footnote n="11">¹¹ http://kurucz.harvard.edu</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Examples of photometric and spectroscopic fits for three DZ white dwarfs. Each row represents one object. For the spectroscopic fits, the observations are shown in black, and the best fit is represented by the red spectrum. The complete set of fits (16 objects) is available in the online journal.
(The complete figure set (4 images) is available).</figcaption></figure>

parameter could perfectly reproduce the distorted Swan bands of all DQpec white dwarfs. Nevertheless, as shown in Figure 9, the $\alpha = 0.2$ value allows a good fit to the Swan bands of carbon-polluted white dwarfs across a large temperature range. The most important disagreement is for WD 1008+290, where the bands could be more rounded. Note, however, that additional mechanisms might play a role, because WD 1008+290 has an intense magnetic field ($B > 10$ MG; Schmidt et al. 1999) that could affect the position of the Swan bands (Liebert et al. 1978; Bues 1991, 1999).

### 3.3.2. Problematic Objects

In our analysis of DQ and DQpec objects, we encountered five cases for which a completely satisfying solution could not be found. Here, we detail the challenges posed by those objects and speculate about possible ways of resolving the discrepancies between models and observations.

GJ 1086—GJ 1086 (G99-37) is one of the only two known white dwarfs to show a CH G band (the other one being BPM 27606). Using a hydrogen abundance of log H/He = -4.3, we are able to achieve a reasonable fit to the G band (Figure 10). However, our solution overestimates the distortion of the Swan

bands (the C₂ bands of our fit are too rounded compared to the observations). This problem suggests that GJ 1086 contains more hydrogen than assumed, because adding hydrogen would reduce the distortion of the Swan bands by decreasing the photospheric density. However, the H/He ratio is already constrained by the strength of the CH G band, and thus we cannot find any solution that simultaneously fits all spectral features. One way to fix this problem would be to change the C₂ and/or CH dissociation equilibrium. After all, pressure effects are known to affect the H/H₂ ratio in dense hydrogen (Vorberger et al. 2007; Holst et al. 2008) and in dense helium-rich mediums (Kowalski 2006), so it is likely that something similar occurs with C₂ and CH in the dense atmospheres of cool carbon-polluted white dwarfs. If the nonideal effects on the dissociation equilibria are such that they reduce the CH/C₂ ratio, they could explain (at least in part) the discrepancy described above. Furthermore, it is also worth noting that GJ 1086 is a magnetic white dwarf (Angel & Landstreet 1974), with $B \approx 7$ MG (Berdyugina et al. 2007; Vornanen et al. 2010). As already stated, the impact of such strong fields on the opacities of carbon-polluted atmospheres remains unclear. Note also that convection might be suppressed (Tremblay et al. 2015; Gentile Fusillo et al. 2018), which would significantly affect the atmosphere structure.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Similar to Figure 6, but for DQ/DQpec white dwarfs. Note that a shift parameter of $\alpha = 0.2$ was assumed for the C$_2$ Swan bands (see Section 3.3.1). The complete set of fits (22 objects) is available in the online journal.
(The complete figure set (5 images) is available).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Synthetic spectra computed with different H/He abundance ratios and for atmospheric parameters similar to those of LHS 290. Note that a shift parameter of $\alpha = 0.2$ was assumed for the C$_2$ Swan bands (see Section 3.3.1).</figcaption></figure>

WD 1235+422—WD 1235+422 is characterized by strong Swan bands that are very well represented by our models (Figure 10). However, our fit to the near-infrared photometry is less satisfying. While our overestimation of the flux in the J and H bands might be due to observational errors, it is also possible that we are missing an absorption source in the infrared. Given the effective temperature of WD 1235+422, H$_2$–He CIA cannot explain this flux depletion. In Section 3.3.3, we explore the possibility that the infrared flux depletion of WD 1235+422 is due to C$_2$–He CIA. It is also worth noting that WD 1235+422 is another magnetic white dwarf (Vornanen et al. 2013).

LHS 1126—Wickramasinghe et al. (1982) were the first to identify the strong near-infrared flux deficit that characterizes the SED of LHS 1126 (GJ 2012). Bergeron et al. (1994) and Giammichele et al. (2012) have explained this flux deficit as being due to H$_2$–He CIA (they find mixed H/He compositions of log H/He = 0.1 and -1.2, respectively). However, their solutions are not compatible with the ultraviolet observations of LHS 1126, as Wolff et al. (2002) found that Ly$\alpha$ is better reproduced with log H/He = -5.5. Another puzzling finding
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Spectroscopic fits to the Swan bands of four DQ/DQpec stars. These four stars represent a sequence of decreasing temperature, which also corresponds to a sequence of increasing photospheric density and of increasing distortion of the Swan bands by high-pressure effects. The red spectra (solid lines) correspond to our best solutions assuming $\alpha = 0.2$, and the blue ones (dashed lines) correspond to the case where the Swan bands are not distorted by high-pressure effects (i.e., $\alpha = 0$). Note that the blue spectra were computed by assuming the same atmospheric parameters as those found by fitting the observational data with a grid of $\alpha = 0.2$ models.</figcaption></figure>

is that the near- and mid-infrared energy distribution of LHS 1126 fits a Rayleigh-Jeans spectrum mimicking a $T_{\text{eff}} > 10^5$ K blackbody (Kilic et al. 2006b). As with GJ 1086 and WD 1235+422, we also have to worry about the usual suspect—magnetism—as spectropolarimetric measurements of LHS 1126 cannot rule out the presence of a $B < 3$ MG magnetic field (Schmidt et al. 1995).

As in all previous analyses of LHS 1126, our fit is far from satisfactory. We interpreted the small depression at ≈5000 Å as being due to distorted Swan bands and, based on the depth of this feature, we concluded that $\log C/He = -8.4$. Regarding our fit to the photometry, we find a hydrogen abundance of $\log H/He = -1.2$ if we adjust the H/He ratio to the near-infrared photometry. As explained above, this solution must be rejected because of the $\text{Ly}\alpha$ analysis of Wolff et al. (2002). The origin of the infrared flux depletion remains unknown, and for this reason, the solution displayed in Figure 10 was obtained assuming a hydrogen-free atmosphere and ignoring the *JHK* bands in our fit to the photometric data. Based on novel ab initio calculations, Kowalski (2014) proposed that He-He-He CIA might explain the infrared flux depletion of LHS 1126. This opacity source is included in our models, but fails to explain the SED of LHS 1126. As noted by Kowalski (2014), this failure may be due to our poor constraints on the ionization equilibrium of helium under high-pressure conditions (Fortov et al. 2003; Kowalski et al. 2007), which controls the photospheric density and therefore the intensity of He-He-He CIA. That being said, our sample contains many helium-rich objects that have denser photospheres than LHS 1126 (e.g., the 16 DC stars cooler than 5000 K for which we find a pure helium composition), and none of them shows a discrepancy similar to that observed here.

WD 1036-204 and WD 1008+290—These two objects, identified in Table 3 as PSO J103855.315-204049.761 and

PSO J101141.465+284550.360, present a very similar problem. For both of them, our models were unable to properly match the photometric observations. In particular, we were forced to ignore the *g* and *r* bands in order to obtain an adequate fit to the other photometric bands. We do not know what explains this discrepancy, but we note that both objects have a very low effective temperature (4530 ± 215 K for WD 1036-204 and 4335 ± 165 K for WD 1008+290) and strong magnetic fields (Schmidt et al. 1999; Jordan & Friedrich 2002) that might affect their structures and opacities.

### 3.3.3. C₂-He CIA

For some DQs in our sample, our fits to the infrared photometry are unsatisfactory (see, for instance, WD 1235+422 in Figure 10). A possible explanation for this mismatch is the omission of an important opacity source in our atmosphere models. Here, we investigate the possibility that C₂-He CIA affects the SEDs of cool carbon-polluted atmospheres. A priori, this hypothesis might seem unlikely because the peak of the C₂-He CIA spectrum is expected to be in the mid-infrared, beyond the *JHK* bandpasses considered here. In fact, the fundamental vibrational band of C₂ is located at 1855 cm⁻¹ (Herzberg 1950), which implies that absorption will be especially important near 5.4 µm. In contrast, H₂-He CIA, which dominates the SED of the coolest mixed H/He white dwarfs (Bergeron & Leggett 2002; Kilic et al. 2012; Gianninas et al. 2015), peaks near 2.3 µm where the fundamental vibrational band of H₂ is located. Nevertheless, it is possible that C₂ overtone bands, which are located at lower wavelengths, are important enough to affect the near-infrared SED of some cool DQpec stars.

To investigate this issue, we use the methodology presented in Blouin et al. (2017) for H₂-He CIA. More precisely, we use ab initio molecular dynamics (MD) to simulate the evolution of
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Similar to Figure 7, but for three problematic carbon-polluted white dwarfs. Note that our fits to WD 1235+422 and LHS 1126 are based on hydrogen-free models (see text for details).</figcaption></figure>

a C₂ molecule in a dense helium medium. In this framework, atoms move according to classical dynamics and DFT is used at each time step to compute the electronic charge density. As shown in Blouin et al. (2017), this methodology is accurate at both low and high densities, where many-body collisions become important.

The simulations were performed with the CPMD<fnref n="12" /> plane-wave DFT code (Marx & Hutter 2000; Hutter et al. 2008). The DFT calculations were carried out using the PBE exchange-correlation functional (Perdew et al. 1996) and ultrasoft pseudopotentials (Vanderbilt 1990). Each simulation consisted of one C₂ molecule surrounded by 15, 31, or 63 He atoms in a cubic box, the length of which was adjusted to obtain the desired density. The density–temperature space was explored with 35 distinct simulations, ranging from T = 4000 to 8000 K and from ρ = 0.08 to 1.4 g cm⁻³. For each simulation, we computed at every time step the dipole moment resulting from the total electronic charge density and the distribution of all nuclei (Silvestrelli et al. 1998; Berghold et al. 2000). At the end of the MD simulations, we obtained the absorption spectra α (ω) using the Fourier transform of the dipole moment time autocorrelation function (Silvestrelli et al. 1997; Kowalski 2014; Blouin et al. 2017).

We took precautions to make sure that the periodic boundary conditions of the simulation cell do not give rise to any artificial absorption features. To do so, we performed simulations with different box sizes while keeping the density constant by adjusting the number of helium atoms in each simulation. We found that absorption spectra obtained from MD simulations performed in a box of at least 10 au (5.3 Å) are virtually identical to those obtained from simulations performed in larger boxes. Therefore, finite-size effects are negligible as long as the simulation cell is at least 10 au large; all simulations presented here satisfy this criterion. Note that we also verified that our simulations have converged with respect to the MD simulation time by comparing absorption spectra obtained for different trajectory lengths. We found that a 32 ps trajectory is usually sufficient to obtain a satisfactory convergence.

Figure 11 shows the results of our DFT-MD calculations for different densities representative of the photosphere of cool carbon-polluted white dwarfs. At low densities, the general shape of the C₂–He CIA spectrum is very similar to that of H₂–He CIA (for comparison, see Figure 4 of Abel et al. 2012), except that the rotational and vibrational bands are located at higher wavelengths. This is a direct consequence of the higher mass of the C₂ molecule. Also similar to H₂–He CIA is the evolution of the absorption spectrum with increasing density. In

<footnote n="12">¹² http://cpmd.org</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** C₂–He CIA spectra for different densities at *T* = 4000 K. All spectra are divided by the number density of C₂ and He.</figcaption></figure>

particular, the fundamental band becomes less and less important, while the rotational band becomes more prominent and shifts toward lower wavelengths (Blouin et al. 2017). For the purpose of this paper, the main result of our simulations is that the overtone bands of the C₂–He CIA are too weak to significantly affect the near-infrared flux of cool carbon-polluted white dwarfs. In fact, we implemented C₂–He CIA in our model atmosphere code—using analytical fits to our results—and we found that it only has an effect for cool ($T_{eff} \le 5000$ K), relatively carbon-rich (log C/He $\ge -5$) models, and that this effect is limited to wavelengths beyond the K band.

# 4. The Spectral Evolution of Cool White Dwarfs

Now that we have determined the atmospheric composition of each of the 501 objects included in our sample (Table 3), we are ready to revisit the spectral evolution of cool white dwarfs. Before going to the results (Section 4.2), we discuss how we can account for selection biases in our sample (Section 4.1).

## 4.1. Correcting for Biases

The quantity of interest for the study of the spectral evolution of white dwarfs is the fraction of stars that have a hydrogen-rich atmosphere in a volume-complete sample, $\rho_H/\rho_{tot} = \rho_H/ (\rho_H + \rho_{He})$.<fnref n="13" /> Unfortunately, the sample considered in this paper is not volume-complete. We decided not to use a volume-complete sample as it would drastically reduce the number of objects in our sample (the almost complete 20 pc sample contains 106 white dwarfs cooler than 8000 K; Hollands et al. 2018) and significantly increase the statistical errors on $\rho_H/\rho_{tot}$.

Given the incompleteness of our sample, a correction must be applied to relate the number of white dwarfs in our sample ($N_H^{sample}$ and $N_{He}^{sample}$) to the unbiased space density ratio. Remember that the four main selection criteria (Section 2.1) are that each star must have a parallax measurement from the Gaia DR2, grizy photometry from Pan-STARRS1, J photometry from 2MASS, and a spectrum to allow an accurate spectral

<table><caption>**Table 5**
$V^{max}$ Correction for the Ratio of Hydrogen-rich to Helium-rich Objects</caption><thead><tr><th><i>T</i><sub>eff</sub> (K)</th><th><i>V</i><sup>max</sup><sub>H</sub> / <i>V</i><sup>max</sup><sub>He</sub></th></tr></thead><tbody><tr><td>8500</td><td>1.15</td></tr><tr><td>8000</td><td>1.16</td></tr><tr><td>7500</td><td>1.17</td></tr><tr><td>7000</td><td>1.19</td></tr><tr><td>6500</td><td>1.21</td></tr><tr><td>6000</td><td>1.23</td></tr><tr><td>5500</td><td>1.22</td></tr><tr><td>5000</td><td>1.20</td></tr><tr><td>4750</td><td>1.18</td></tr><tr><td>4500</td><td>1.15</td></tr><tr><td>4250</td><td>1.09</td></tr><tr><td>4000</td><td>1.01</td></tr><tr><td>3750</td><td>0.82</td></tr></tbody></table>

classification. By far, the two main limiting factors out of those four criteria are the need for J photometry from 2MASS and the need for spectroscopic observations. Below, we examine how each of these two criteria can induce a bias on $N_H^{sample}/N_{tot}^{sample}$ and how we can correct for those biases to obtain the unbiased space density ratio $\rho_H/\rho_{tot}$.

If we ignore for the moment the criterion on spectroscopic observations, our sample is effectively limited by the 2MASS J band, because Pan-STARRS and Gaia have much deeper fields than 2MASS. Depending on the effective temperature, a hydrogen-rich white dwarf can be more or less luminous in the J band than a helium-rich object with the same effective temperature and surface gravity. This can induce a selection bias in favor of hydrogen-rich or helium-rich objects. To correct for this bias, we use the $V^{max}$ correction (Schmidt 1975), where $V^{max}$ is the volume defined by the maximum distance at which a given object would still appear in our sample. More specifically, the quantity of interest is $V_H^{max}/V_{He}^{max}$, the ratio between the observed volume of hydrogen-rich objects and that of helium-rich objects. Using our model fluxes, we computed the ratio of the luminosities and, from there, we obtained $V_H^{max}/V_{He}^{max}$. The corresponding $V^{max}$ ratios are given in Table 5. Note that we assumed log g = 8 to compute $V_H^{max}/V_{He}^{max}$, which is justified by the very similar mass distributions of hydrogen-rich and helium-rich white dwarfs (Giammichele et al. 2012). The decrease of $V_H^{max}/V_{He}^{max}$ at low temperatures is a direct consequence of the onset of CIA, which significantly enhances the opacity in the J band for hydrogen-rich objects. With the data of Table 5, we can relate the number of hydrogen-rich and helium-rich white dwarfs in our sample to the space density ratio,

$$ \frac{\rho_H}{\rho_H + \rho_{He}} = \left( 1 + \frac{N_{He}^{sample} V_H^{max}}{N_H^{sample} V_{He}^{max}} \right)^{-1} . \quad (6) $$

Let us now turn to the bias induced by the requirement of having spectroscopic observations for every object in our sample. An appreciable fraction of our spectra (158 out of 501) are from the Sloan Digital Sky Survey (SDSS), which is plagued by numerous selection effects. The SDSS consists of first imaging the sky in five bandpasses and then using this photometry to select targets that deserve spectroscopic observations. This target selection process is based on many different selection criteria that overlap with one another (Harris et al. 2003; Eisenstein et al. 2006a). Hence, modeling the selection biases induced by the SDSS is an intractable task.

<footnote n="13">¹³ In this work, an atmosphere with a H/He abundance ratio greater than 1 is considered hydrogen-rich.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** The total length of each bar indicates the number of objects in our sample that fall in each effective temperature bin. The dark blue portion of each bar represents the fraction of those objects that were spectroscopically observed in the SDSS.</figcaption></figure>

Fortunately, we do not need to worry about these selection effects for objects cooler than 6000 K, because for the vast majority of those cool objects we rely on spectroscopic observations from other sources (Figure 12).<fnref n="14" /> However, above 6000 K, SDSS objects are much more abundant and we might be at the mercy of SDSS selection effects. To evaluate the severity of those selection effects, we computed the fraction of hydrogen-rich stars for objects of our sample that were spectroscopically observed by the SDSS and for those that were not. We found absolutely no difference for the 6000–8000 K temperature range as a whole (0.80 ± 0.04 for the SDSS objects and 0.82 ± 0.03 for the rest).<fnref n="15" /> However, there are differences if we look at smaller temperature bins. For the 7000–8000 K bin, we found hydrogen-rich fractions of 0.91 ± 0.04 for the SDSS objects and 0.84 ± 0.05 for the rest, and, for the 6000–7000 K bin, we found 0.68 ± 0.06 and 0.80 ± 0.04. Strictly speaking, these differences between SDSS and non-SDSS objects are not statistically significant. The probability that, given the finite size of our samples, the fraction of hydrogen-rich objects is the same for stars observed by the SDSS and for those not observed by the SDSS (i.e., the *p*-value) is 0.26 for the 7000–8000 K bin and 0.11 for the 6000–7000 K bin. Nevertheless, it is worth keeping in mind that the SDSS might induce a small bias in favor of hydrogen-rich objects between 7000 and 8000 K and a small bias in favor of helium-rich objects between 6000 and 7000 K

## 4.2. Results and Discussion

Figure 13 shows the evolution of the $\rho_H/\rho_{tot}$ ratio as a function of decreasing effective temperature and compares our results to those of previous studies. Down to $T_{eff} = 5000$ K, our results are in good agreement with previous studies—

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Fraction of hydrogen-rich white dwarfs as a function of effective temperature for 1000 K bins. Different colors represent data taken from different studies (see legend). The error bars indicate the $1\sigma$ uncertainty associated with the finite number of objects in each bin. Note that the results of our sample were corrected using Equation (6) and the data of Table 5.</figcaption></figure>

particularly with Limoges et al. (2015), whose results were based on a bigger sample than the six other studies shown in Figure 13. It is only below 5000 K that our results stand out from the rest. For the 4000–5000 K temperature bin, we find $\rho_H/\rho_{tot} = 0.69 \pm 0.05$, while other studies found either a much lower ($\approx 0.4$; Bergeron et al. 1997, 2001; Kilic et al. 2010) or a much higher fraction ($\gtrsim 0.95$; Kowalski 2006; Kilic et al. 2009a; Limoges et al. 2015).

Two factors contribute to this difference. On the one hand, in contrast to Kowalski (2006) and Kilic et al. (2009a, 2009b), our sample includes metal-polluted white dwarfs. As those objects are all helium-rich in the 4000–5000 K temperature bin, they contribute to give rise to a lower $\rho_H/\rho_{tot}$ ratio. On the other hand, our models (Paper I) include the nonideal input physics (in particular the improved Ly$\alpha$ opacities; Kowalski & Saumon 2006) that led Kowalski (2006), Kilic et al. (2009a, 2009b), and Limoges et al. (2015) to conclude that the vast majority of cool DC white dwarfs are hydrogen-rich. This factor explains why we find a higher $\rho_H/\rho_{tot}$ ratio than Bergeron et al. (1997, 2001) and Kilic et al. (2010), who did not include such high-density nonideal effects in their atmosphere models. Note, however, that our conclusions are less extreme than Kowalski (2006), as we find that the SEDs of $\approx 25\%$ of DC white dwarfs cooler than 5000 K are better fitted with helium-rich models (Section 3.1). This is the expected result. We know many DZ white dwarfs cooler than $T_{eff} = 5000$ K that must have a helium-rich atmosphere in order to explain their broad spectral lines (for examples, see Dufour et al. 2007a, Hollands et al. 2017, and Paper III). As the pollution of a white dwarf by rocky debris is independent from the evolution of the white dwarf itself, a significant number of nonpolluted counterparts (i.e., helium-rich DCs) must exist.

From Figure 13, the fraction of hydrogen-rich stars seems almost consistent with an evolution at a fixed $\rho_H/\rho_{tot}$ ratio. However, the picture looks quite different if we use smaller temperature bins. Figure 14 shows the same plot as Figure 13, but, this time, 500 K temperature bins were used. To make sure that no feature is missed by our arbitrary binning choice, we also added a continuous $\rho_H/\rho_{tot}$ versus $T_{eff}$ curve to Figure 14

<footnote n="14">¹⁴ For the most part, objects that are not in the SDSS are from the Limoges et al. (2015) sample. White dwarfs in that sample were selected from reduced proper motion diagrams (Limoges et al. 2013) and are therefore not subject to biases in favor of hydrogen-rich or helium-rich objects.</footnote>

<footnote n="15">¹⁵ The uncertainties are due to the finite size of the samples and are estimated as $\sigma_f = \sqrt{f(1-f)/N}$, where N is the sample size and f is the fraction of objects that are hydrogen-rich. Note that $f \pm \sigma_f$ corresponds to a 68% confidence interval.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Fraction of hydrogen-rich white dwarfs as a function of effective temperature for fixed 500 K bins (in blue) and for a 500 K moving bin (in gray; see the text for details). As in Figure 13, the error bars indicate the 1σ uncertainty associated with the finite number of objects in each bin, and the correction given by Equation (6) was applied.</figcaption></figure>

(in gray). This continuous curve was obtained by computing $\rho_H/\rho_{tot}$ within a 500 K moving bin, which eliminates the arbitrariness of the choice of bin boundaries. Note also that the weight of each star in the calculation of $\rho_H/\rho_{tot}$ was computed using a Gaussian with a standard deviation given by the uncertainty on $T_{eff}$.

Figure 14 reveals a much more complex picture than Figure 13. The hydrogen-rich fraction decreases from $T_{eff} \approx 7500$ K until $T_{eff} \approx 6250$ K, then increases until $T_{eff} \approx 5000$ K, and decreases again below 5000 K. To confirm that these fluctuations of $\rho_H/\rho_{tot}$ are statistically significant, we compared the effective temperature distributions of hydrogen-rich and helium-rich stars in our sample (Figure 15). An Anderson-Darling test shows that hydrogen-rich and helium-rich white dwarfs do not follow the same temperature distribution, as the probability that both temperature distributions belong to the same population is 0.002.<fnref n="16" /> Moreover, note that these conclusions remain unchanged if we take into account the $V_H^{max}/V_{He}^{max}$ corrections given in Table 5. If we randomly add or remove hydrogen-rich stars in each temperature bin (where the number of additions or withdrawals is determined by the $V_H^{max}/V_{He}^{max}$ ratio), we find that the probability that both the hydrogen-rich and helium-rich temperature distributions come from the same distribution is $\approx 0.004$. Therefore, it is highly unlikely that the large increases and decreases of $\rho_H/\rho_{tot}$ with respect to $T_{eff}$ (Figure 14) are due to random fluctuations associated with small number statistics.

### 4.2.1. Behavior above 5000 K

As explained above, SDSS selection effects could lead to a bias in favor of hydrogen-rich objects between 7000 and 8000 K, and a bias in favor of helium-rich objects between 6000 and 7000 K (Section 4.1). To check if the decrease of

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Effective temperature distribution of hydrogen-rich and helium-rich stars in our sample. An Anderson-Darling test shows that the probability that both distributions belong to the same parent population is 0.002.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Fraction of hydrogen-rich white dwarfs as a function of effective temperature for a 500 K moving bin. The results obtained using our whole sample are shown in gray and those obtained by excluding SDSS objects are shown in red.</figcaption></figure>

$\rho_H/\rho_{tot}$ between 7500 and 6250 K is solely due to this bias, we compare the evolution of the hydrogen-rich fraction obtained using all stars in our sample to that obtained from a subsample that excludes all SDSS objects (Figure 16). Although both samples have different behaviors at high temperatures, they both show a clear decrease of $\rho_H/\rho_{tot}$ between 7000 and 6250 K, suggesting that this feature is not due to SDSS biases.

Another potential problem with our analysis of $\rho_H/\rho_{tot}$ is that it includes several objects with a very low mass (Figure 17). Because it is impossible for a star with a mass below $\approx 0.45 M_{\odot}$ to become a white dwarf through single-star evolution within the age of the universe, those objects are likely part of binary systems (Liebert et al. 2005; Rebassa-Mansergas et al. 2011). Figure 18 compares the hydrogen-rich fractions that we find if we include and if we exclude those low-mass objects. Clearly, the evolution of $\rho_H/\rho_{tot}$ above 5000 K is the same for both samples, so the contamination of our sample by binary systems

<footnote n="16">¹⁶ The corresponding probability given by the more usual Kolmogorov-Smirnov test is significantly higher (0.013). However, the Anderson-Darling test is more appropriate, because it is more sensitive than the Kolmogorov-Smirnov test when the differences between both distributions are more important near their extremities (Feigelson & Babu 2012), which is precisely the case here.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Mass of white dwarfs in our sample as a function of their effective temperatures. Typical uncertainties on *M* and *T*<sub>eff</sub> are shown in the top-left corner.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Fraction of hydrogen-rich white dwarfs as a function of effective temperature for a 500 K moving bin. The results obtained using our whole sample are shown in gray and those obtained by excluding objects with $M \le 0.45 M_\odot$ are shown in blue.</figcaption></figure>

cannot explain the decrease of $\rho_{H}/\rho_{tot}$ in the 7500–6250 K range and the increase in the 6250–5000 K range.

Another way to look at the evolution of the hydrogen-rich fraction is to look at $\rho_{H}/\rho_{tot}$ as a function of the cooling age (Figure 19). While the horizontal axis is distorted compared to Figure 14, the decrease of $\rho_{H}/\rho_{tot}$ before 6250 K and the increase down to 5000 K are still clearly visible.<fnref n="17" /> This shows that the trends observed in the $\rho_{H}/\rho_{tot}$ versus $T_{eff}$ plot (Figure 14) are representative of the temporal evolution of cool white dwarfs.

All things considered, the decrease of the hydrogen-rich fraction from $T_{eff} \approx 7500$ K to $T_{eff} \approx 6250$ K and the subsequent increase down to $T_{eff} \approx 5000$ K appear to be real. The decrease of $\rho_{H}/\rho_{tot}$ is consistent with what we expect from convective mixing. The deepening of the convection zone with decreasing effective temperature mixes the superficial

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Fraction of hydrogen-rich white dwarfs as a function of cooling age for a 1 Gyr moving bin (in gray) and for fixed 1 Gyr bins (in blue).</figcaption></figure>

hydrogen layer with the more massive helium layer underneath, which can turn DAs into non-DAs (Tassoul et al. 1990; Bergeron et al. 1997; Rolland et al. 2018). However, the increase of $\rho_{H}/\rho_{tot}$ in the 6250–5000 K range is much more intriguing. This increase is similar to that observed at the blue edge of the non-DA gap (Bergeron et al. 1997, 2001), albeit located at slightly lower temperatures. No satisfactory physical explanation for this increase can be found in the literature (Hansen 1999; Malo et al. 1999; Bergeron et al. 2001), and we do not have any new scenario to propose.

### 4.2.2. Behavior below 5000 K

Let us now examine the decrease of $\rho_{H}/\rho_{tot}$ below 5000 K. The first thing to note is that this feature cannot be affected by SDSS selection biases, because there are very few SDSS objects with $T_{eff} < 5000$ K in our sample (see Figures 12 and 16). Second, we note that the decrease of the hydrogen-rich fraction becomes less obvious if we eliminate the low-mass objects from our sample. For the blue curve of Figure 18, the decline of $\rho_{H}/\rho_{tot}$ between 5000 and 4000 K is barely significant. Third, no significant counterpart for the decrease below 5000 K in the $\rho_{H}/\rho_{tot}$ versus $T_{eff}$ figure is visible if we plot $\rho_{H}/\rho_{tot}$ as a function of the cooling age (Figure 19). Therefore, it is not clear that the decrease of $\rho_{H}/\rho_{tot}$ below 5000 K reflects the temporal evolution of cool white dwarfs. Instead, the finding that helium-rich white dwarfs become more abundant below 5000 K may simply reflect the fact that they cool down faster than their hydrogen-rich counterparts.

In any case, our results rule out the possibility that cool white dwarfs all become hydrogen-rich as a result of hydrogen accretion from the interstellar medium. This scenario was already challenged by the fact that the accretion rate required for this conversion ($\approx 6 \times 10^{-17} M_\odot \text{yr}^{-1}$; Kowalski 2006) is incompatible with the limits on the accretion rates of cool helium-rich DA/DZA white dwarfs ($10^{-20}-10^{-17} M_\odot \text{yr}^{-1}$; Dufour et al. 2007a; Rolland et al. 2018).

## 5. Conclusion

A detailed photometric and spectroscopic analysis of a homogeneous sample of 501 cool white dwarfs was presented.

<footnote n="17">¹⁷ A hydrogen-rich white dwarf with log g = 8 has a cooling age of 2.0 Gyr when $T_{eff}$ = 6250 K and of 5.7 Gyr when $T_{eff}$ = 5000 K (Fontaine et al. 2001).</footnote>
Our analysis was based on a state-of-the-art model atmosphere code that includes all of the nonideal input physics required to accurately model the dense atmospheres of cool white dwarfs. As our models are the first to successfully model the most challenging cool DZ white dwarfs (Papers I, II, and III) and as our analysis makes use of the largest homogeneous sample of cool white dwarfs studied to date, our results are on firmer grounds than previous attempts at constraining the spectral evolution of cool white dwarfs.

A satisfactory fit was found for all 501 white dwarfs studied in this work, except for a handful of carbon-polluted objects. We explored a few avenues to improve our fits to those stars, such as empirically adjusting the pressure shift of the distorted Swan bands and computing the absorption resulting from C₂–He collisions. We also revisited a peculiar class of DC stars (known as peculiar non-DAs) for which spectroscopic and photometric observations suggest conflicting chemical compositions. We reported evidence that is at odds with the interpretation according to which those objects form a real physical group.

We found that hydrogen-rich white dwarfs become less abundant when the effective temperature decreases from 7500 to 6250 K, which is likely explained by convective mixing. From 6250 to 5000 K, the trend is reversed and the fraction of white dwarfs with a hydrogen-rich atmosphere increases. So far, no physical scenario is available to explain this observation. Finally, at lower temperatures, we find that hydrogen-rich white dwarfs become rarer from 5000 to 4000 K. This trend invalidates the scenario according to which accretion of hydrogen from the interstellar medium dominates the spectral evolution of cool white dwarfs.

We wish to thank Pierre Bergeron for many enlightening discussions on the spectral evolution of white dwarfs. We are grateful to Mukremin Kilic for sharing with us his spectrum of SDSS J134118.68+022737.0. We also thank Adela Kawka and Stephane Vennes for allowing us to use their VLT/FORS2 spectra of NLTT 8733, NLTT 14553, and NLTT 57760.

This work was supported in part by the NSERC (Canada) and the Fund FRQNT (Québec). This work has made use of the Montreal White Dwarf Database (Dufour et al. 2017).

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://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.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No.

AST-1238877, the University of Maryland, Eötvös Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

**Software:** CPMD (Marx & Hutter 2000), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011).

## ORCID iDs

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

## References

Abel, M., Frommhold, L., Li, X., & Hunt, K. L. C. 2012, JChPh, 136, 044319

Akaike, H. 1974, ITAC, 19, 716

Allard, N. F., & Alekseev, V. A. 2014, AdSpR, 54, 1248

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

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

Allard, N. F., Kielkopf, J. F., Blouin, S., et al. 2018, A&A, 619, A152

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

Angel, J. R. P. 1978, ARA&A, 16, 487

Angel, J. R. P., Borra, E. F., & Landstreet, J. D. 1981, ApJS, 45, 457

Angel, J. R. P., Illing, R. M. E., & Landstreet, J. D. 1972, ApJL, 175, L85

Angel, J. R. P., & Landstreet, J. D. 1974, ApJ, 191, 457

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

Berdyugin, A. V., & Pirola, V. 1999, A&A, 352, 619

Berdyugina, S. V., Berdyugin, A. V., & Pirola, V. 2007, PhRvL, 99, 091101

Bergeron, P., Dufour, P., Fontaine, G., et al. 2019, ApJ, 876, 67

Bergeron, P., & Leggett, S. K. 2002, ApJ, 580, 1070

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

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

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

Bergeron, P., Ruiz, M.-T., Leggett, S. K., Saumon, D., & Wesemael, F. 1994, ApJ, 423, 456

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

Berghold, G., Mundy, C. J., Romero, A. H., Hutter, J., & Parrinello, M. 2000, PhRvB, 61, 10040

Blouin, S., Allard, N. F., Leininger, T., Gadéa, F. X., & Dufour, P. 2019a, ApJ, 875, 137

Blouin, S., Dufour, P., & Allard, N. F. 2018a, ApJ, 863, 184

Blouin, S., Dufour, P., Allard, N. F., & Kilic, M. 2018b, ApJ, 867, 161

Blouin, S., Dufour, P., Allard, N. F., & Salim, S. 2019b, ApJ, 872, 188

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

Bues, I. 1991, in White Dwarfs, NATO Advanced Science Institutes (ASI) Series C, Vol. 336, ed. G. Vauclair & E. Sion (Dordrecht: Kluwer), 285

Bues, I. 1999, in ASP Conf. Ser. 169, 11th European Workshop on White Dwarfs, ed. S.-E. Solheim & E. G. Meistas (San Francisco, CA: ASP), 240

Burnham, K. P., & Anderson, D. R. 2002, Model Selection and Inference (New York: Springer)

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Chayer, P., Fontaine, G., & Wesemael, F. 1995, ApJS, 99, 189

Dantona, F., & Mazzitelli, I. 1979, A&A, 74, 161

Dufour, P., Bergeron, P., & Fontaine, G. 2005, ApJ, 627, 404

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

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

Dufour, P., Liebert, J., Fontaine, G., & Behara, N. 2007b, Natur, 450, 522

Dupuis, J., Fontaine, G., Pelletier, C., & Wesemael, F. 1993, ApJS, 84, 73

Eisenstein, D. J., Liebert, J., Harris, H. C., et al. 2006a, ApJS, 167, 40

Eisenstein, D. J., Liebert, J., Koester, D., et al. 2006b, AJ, 132, 676

Farihi, J., Barstow, M. A., Redfield, S., Dufour, P., & Hambly, N. C. 2010, MNRAS, 404, 2123

Feigelson, E. D., & Babu, G. J. 2012, Modern Statistical Methods for Astronomy (Cambridge: Cambridge Univ. Press)

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

Fontaine, G., & Wesemael, F. 1987, in IAU Coll. 95. Second Conf. on Faint Blue Stars, ed. A. G. D. Philip, D. S. Hayes, & J. W. Liebert (Schenectady, NY: Davis Press), 319

Fortov, V. E., Ternovoi, V. Y., Zhernokletov, M. V., et al. 2003, JETP, 97, 259

Gaia Collaboration 2016, A&A, 595, A1

Gaia Collaboration 2018, A&A, 616, A1

García-Berro, E., Torres, S., Althaus, L. G., et al. 2010, Natur, 465, 194
Gentile Fusillo, N. P., Tremblay, P.-E., Jordan, S., et al. 2018, MNRAS, 473, 3693

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

Gianninas, A., Bergeron, P., & Ruiz, M. T. 2011, ApJ, 743, 138

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

Graham, J. R., Matthews, K., Neugebauer, G., & Soifer, B. T. 1990, ApJ, 357, 216

Hall, P. B., & Maxwell, A. J. 2008, ApJ, 678, 1292

Hansen, B. M. S. 1999, ApJ, 520, 680

Harris, H. C., Liebert, J., Kleinman, S. J., et al. 2003, AJ, 126, 1023

Herzberg, G. 1950, Molecular Spectra and Molecular Structure. Vol. 1: Spectra of Diatomic Molecules (Princeton, NJ: Van Nostrand)

Hollands, M. A., Koester, D., Alekseev, V., Herbert, E. L., & Gänsicke, B. T. 2017, MNRAS, 467, 4970

Hollands, M. A., Tremblay, P.-E., Gänsicke, B. T., Gentile-Fusillo, N. P., & Toonen, S. 2018, MNRAS, 480, 3942

Holst, B., Redmer, R., & Desjarlais, M. P. 2008, PhRvB, 77, 184201

Hornkohl, J. O., Parigger, C. G., & Nemes, L. 2005, ApOpt, 44, 3686

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

Hurvich, C. M., & Tsai, C.-L. 1989, Biometrika, 76, 297

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

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

Jordan, S., & Friedrich, S. 2002, A&A, 383, 519

Jura, M. 2003, ApJL, 584, L91

Kalirai, J. S. 2012, Natur, 486, 90

Kawka, A., & Vennes, S. 2006, ApJ, 643, 402

Kawka, A., & Vennes, S. 2012, MNRAS, 425, 1394

Kepler, S. O., Koester, D., & Ourique, G. 2016a, Sci, 352, 67

Kepler, S. O., Pelisoli, I., Koester, D., et al. 2015, MNRAS, 446, 4078

Kepler, S. O., Pelisoli, I., Koester, D., et al. 2016b, MNRAS, 455, 3413

Kilic, M., Kowalski, P. M., Reach, W. T., & von Hippel, T. 2009a, ApJ, 696, 2094

Kilic, M., Kowalski, P. M., & von Hippel, T. 2009b, AJ, 138, 102

Kilic, M., Leggett, S. K., Tremblay, P.-E., et al. 2010, ApJS, 190, 77

Kilic, M., Munn, J. A., Harris, H. C., et al. 2006a, AJ, 131, 582

Kilic, M., Thorstensen, J. R., Kowalski, P. M., & Andrews, J. 2012, MNRAS, 423, L132

Kilic, M., von Hippel, T., Mullally, F., et al. 2006b, ApJ, 642, 1051

Kleinman, S. J., Kepler, S. O., Koester, D., et al. 2013, ApJS, 204, 5

Koester, D. 1976, A&A, 52, 415

Koester, D. 2009, A&A, 498, 517

Koester, D., & Knist, S. 2006, A&A, 454, 951

Koester, D., & Wilken, D. 2006, A&A, 453, 1051

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

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

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

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

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

Kowalski, P. M., & Saumon, D. 2006, ApJL, 651, L137

Leggett, S. K., Bergeron, P., Subasavage, J. P., et al. 2018, ApJS, 239, 26

Liebert, J. 1986, in IAU Coll. 87, Hydrogen Deficient Stars and Related Objects, ed. K. Hunger, D. Schoenberner, & N. Kameswara Rao (Dordrecht: D. Reidel), 367

Liebert, J., Angel, J. R. P., Stockman, H. S., & Beaver, E. A. 1978, ApJ, 225, 181

Liebert, J., Bergeron, P., & Holberg, J. B. 2005, ApJS, 156, 47

Liebert, J., Wesemael, F., Hansen, C. J., et al. 1986, ApJ, 309, 241

Limoges, M.-M., Bergeron, P., & Lépine, S. 2015, ApJS, 219, 19

Limoges, M.-M., Lépine, S., & Bergeron, P. 2013, AJ, 145, 136

Lodders, K. 2003, ApJ, 591, 1220

MacDonald, J., & Vennes, S. 1991, ApJ, 371, 719

Malo, A., Wesemael, F., & Bergeron, P. 1999, ApJ, 517, 901

Marx, D., & Hutter, J. 2000, in Modern Methods and Algorithms of Quantum Chemistry, Vol. 1, ed. J. Grotendorst (Jülich: John von Neumann Institute for Computing), 1

Masseron, T., Plez, B., Van Eck, S., et al. 2014, A&A, 571, A47

Muchmore, D. 1984, ApJ, 278, 769

Oswalt, T. D., Smith, J. A., Wood, M. A., & Hintzen, P. 1996, Natur, 382, 692

Paquette, C., Pelletier, C., Fontaine, G., & Michaud, G. 1986, ApJS, 61, 197

Pelletier, C., Fontaine, G., Wesemael, F., Michaud, G., & Wegner, G. 1986, ApJ, 307, 242

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

Putney, A. 1997, ApJS, 112, 527

Putney, A., & Jordan, S. 1995, ApJ, 449, 863

Rebassa-Mansergas, A., Nebot Gómez-Morán, A., Schreiber, M. R., Girven, J., & Gänsicke, B. T. 2011, MNRAS, 413, 1121

Renedo, I., Althaus, L. G., Miller Bertolami, M. M., et al. 2010, ApJ, 717, 183

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

Rolland, B., Bergeron, P., & Fontaine, G. 2018, ApJ, 857, 56

Saumon, D., Holberg, J. B., & Kowalski, P. M. 2014, ApJ, 790, 50

Sayres, C., Subasavage, J. P., Bergeron, P., et al. 2012, AJ, 143, 103

Schatzman, E. 1945, AnAp, 8, 143

Schmidt, G. D., Bergeron, P., & Fegley, B. 1995, ApJ, 443, 274

Schmidt, G. D., Liebert, J., Harris, H. C., Dahn, C. C., & Leggett, S. K. 1999, ApJ, 512, 916

Schmidt, M. 1975, ApJ, 202, 22

Silvestrelli, P. L., Bernasconi, M., & Parrinello, M. 1997, CPL, 277, 478

Silvestrelli, P. L., Marzari, N., Vanderbilt, D., & Parrinello, M. 1998, SSCom, 107, 7

Sion, E. M. 1984, ApJ, 282, 612

Subasavage, J. P., Henry, T. J., Bergeron, P., et al. 2007, AJ, 134, 252

Subasavage, J. P., Henry, T. J., Bergeron, P., Dufour, P., & Hambly, N. C. 2008, AJ, 136, 899

Subasavage, J. P., Jao, W.-C., Henry, T. J., et al. 2009, AJ, 137, 4547

Subasavage, J. P., Jao, W.-C., Henry, T. J., et al. 2017, AJ, 154, 32

Tassoul, M., Fontaine, G., & Winget, D. E. 1990, ApJS, 72, 335

Tremblay, P.-E., & Bergeron, P. 2008, ApJ, 672, 1144

Tremblay, P.-E., Fontaine, G., Freytag, B., et al. 2015, ApJ, 812, 19

Tremblay, P.-E., Kalirai, J. S., Soderblom, D. R., Cignoni, M., & Cummings, J. 2014, ApJ, 791, 92

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

Vanderbilt, D. 1990, PhRvB, 41, 7892

Vauclair, G., & Reisse, C. 1977, A&A, 61, 415

Vorberger, J., Tamblyn, I., Militzer, B., & Bonev, S. A. 2007, PhRvB, 75, 024206

Vornanen, T., Berdyugina, S. V., & Berdyugin, A. 2013, A&A, 557, A38

Vornanen, T., Berdyugina, S. V., Berdyugin, A. V., & Piirola, V. 2010, ApJL, 720, L52

Wickramasinghe, D. T., Allen, D. A., & Bessell, M. S. 1982, MNRAS, 198, 473

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

Zeidler-K. T., E. M., & Koester, D. 1982, A&A, 113, 173