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

The fundamental parameters of DA white dwarfs are obtained with the photometric technique (Bergeron et al. 1997). The solid angle π(R/D)² and T_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 $\mathrm{log}\,{\rm{H}}/\mathrm{He}=2$, the differences between a pure hydrogen atmosphere and a helium-dominated atmosphere ($\mathrm{log}\,{\rm{H}}/\mathrm{He}\lt 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.

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 ≲ 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.

Zoom In Zoom Out Reset image size

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 χ². 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 χ² 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)

$$\begin{eqnarray}&&\mathrm{AIC}=2k-2\mathrm{ln}L,\end{eqnarray} \tag{ 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 $\mathrm{ln}L=C-{\chi }^{2}/2$, where C is a constant that is independent from the model and can be ignored. When we fit 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),

$$\begin{eqnarray}&&\mathrm{AICc}=\mathrm{AIC}+\displaystyle \frac{2k(k+1)}{n-k-1},\end{eqnarray} \tag{ 2 }$$

where n is the sample size (in our case the sample size corresponds to the number of bands included in our photometric 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).⁵ 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\lt \mathrm{log}\,{\rm{H}}/\mathrm{He}\lt 0.5$).⁶ This numerical experiment demonstrates that this approach is 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 $\mathrm{log}\,{\rm{H}}/\mathrm{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.⁷ The cases where our AIC approach fails to retrieve a mixed solution are all for objects with hydrogen abundances near $\mathrm{log}\,{\rm{H}}/\mathrm{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,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}(\mathrm{AIC})={\mathrm{AIC}}_{i}-\min \,\mathrm{AIC}.\end{eqnarray} \tag{ 3 }$$

From there, the Akaike weights w_i are computed as

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{\exp \left[-\tfrac{1}{2}{{\rm{\Delta }}}_{i}(\mathrm{AIC})\right]}{{\displaystyle \sum }_{j}\exp \left[-\tfrac{1}{2}{{\rm{\Delta }}}_{j}(\mathrm{AIC})\right]}.\end{eqnarray} \tag{ 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.

For metal-polluted objects, we still rely on a photometric fit (see Section 2.3.1) to determine T_eff and $\mathrm{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).⁸ 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.

In the 6200 K ≤ T_eff ≤ 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α 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α (see Table 4).⁹ 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_eff > 5000 K (that is, objects for which we should see Hα 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 ± 6)% of those 45 DCs, which is perfectly consistent with our actual error rate in this temperature range of 27% (12 out of 45).¹⁰

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 ≈1σ (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.

Zoom In Zoom Out Reset image size

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 χ² (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 ∼ 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.

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_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_eff and $\mathrm{log}\,g$ determinations are uncertain.

To see by how much T_eff and $\mathrm{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 $\mathrm{log}\,{\rm{H}}/\mathrm{He}=-5$. Note that this hydrogen abundance is too small to produce any detectable Hα feature or to give rise to significant H₂−He CIA. As shown in Figure 5, the addition of this trace amount of hydrogen has an important effect on T_eff and $\mathrm{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_eff). This behavior is analogous to the reduction of the photometric T_eff and $\mathrm{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 $\mathrm{log}\,{\rm{H}}/\mathrm{He}=-5$ solutions, and we adjusted the uncertainties so that they now encompass both solutions. The T_eff, $\mathrm{log}\,g$, M, and τ_cool values reported in Table 3 include these corrections.

Zoom In Zoom Out Reset image size

While we rejected the idea that an undetected amount of hydrogen could bring the density down to a level where the α = 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 Δ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 ΔT_e(ρ), 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 α using the observed spectra of distorted Swan bands. After testing several α parameters, we found that the value that leads—for most objects—to the best agreement between models and observations is α = 0.2. This value is actually a compromise, because we found that no fixed α parameter could perfectly reproduce the distorted Swan bands of all DQpec white dwarfs. Nevertheless, as shown in Figure 9, the α = 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\gt 10\,\mathrm{MG}$; Schmidt et al. 1999) that could affect the position of the Swan bands (Liebert et al. 1978; Bues 1991, 1999).

Zoom In Zoom Out Reset image size

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 $\mathrm{log}\,{\rm{H}}/\mathrm{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 ≈ 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.

Zoom In Zoom Out Reset image size

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₂−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₂−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₂−He CIA (they find mixed H/He compositions of $\mathrm{log}\,{\rm{H}}/\mathrm{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α is better reproduced with $\mathrm{log}\,{\rm{H}}/\mathrm{He}=-5.5$. Another puzzling finding is that the near- and mid-infrared energy distribution of LHS 1126 fits a Rayleigh–Jeans spectrum mimicking a T_eff > 10⁵ 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 $\mathrm{log}\,{\rm{C}}/\mathrm{He}=-8.4$. Regarding our fit to the photometry, we find a hydrogen abundance of $\mathrm{log}\,{\rm{H}}/\mathrm{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 Lyα 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.

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 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¹² 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 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 ≤ 5000 K), relatively carbon-rich ($\mathrm{log}\,{\rm{C}}/\mathrm{He}\geqslant -5$) models, and that this effect is limited to wavelengths beyond the K band.

Zoom In Zoom Out Reset image size

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 ρ_H/ρ_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 ρ_H/ρ_tot between 7000 and 6250 K, suggesting that this feature is not due to SDSS biases.

Zoom In Zoom Out Reset image size

Another potential problem with our analysis of ρ_H/ρ_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 ≈0.45 M_⊙ 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 ρ_H/ρ_tot above 5000 K is the same for both samples, so the contamination of our sample by binary systems cannot explain the decrease of ρ_H/ρ_tot in the 7500–6250 K range and the increase in the 6250–5000 K range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another way to look at the evolution of the hydrogen-rich fraction is to look at ρ_H/ρ_tot as a function of the cooling age (Figure 19). While the horizontal axis is distorted compared to Figure 14, the decrease of ρ_H/ρ_tot before 6250 K and the increase down to 5000 K are still clearly visible.¹⁷ This shows that the trends observed in the ρ_H/ρ_tot versus T_eff plot (Figure 14) are representative of the temporal evolution of cool white dwarfs.

Zoom In Zoom Out Reset image size

All things considered, the decrease of the hydrogen-rich fraction from T_eff ≈ 7500 K to T_eff ≈ 6250 K and the subsequent increase down to T_eff ≈ 5000 K appear to be real. The decrease of ρ_H/ρ_tot is consistent with what we expect from convective mixing. The deepening of the convection zone with decreasing effective temperature mixes the superficial 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 ρ_H/ρ_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.

Let us now examine the decrease of ρ_H/ρ_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 ρ_H/ρ_tot between 5000 and 4000 K is barely significant. Third, no significant counterpart for the decrease below 5000 K in the ρ_H/ρ_tot versus T_eff figure is visible if we plot ρ_H/ρ_tot as a function of the cooling age (Figure 19). Therefore, it is not clear that the decrease of ρ_H/ρ_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 (≈ 6 × 10⁻¹⁷ M_⊙ yr⁻¹; Kowalski 2006) is incompatible with the limits on the accretion rates of cool helium-rich DA/DZA white dwarfs (10⁻²⁰–10⁻¹⁷ M_⊙ yr⁻¹; Dufour et al. 2007a; Rolland et al. 2018).