On the Nature of Ultracool White Dwarfs: Not so Cool Afterall

A recent analysis of the 100 pc white dwarf sample in the SDSS footprint demonstrated for the first time the existence of a well defined ultracool -- or IR-faint -- white dwarf sequence in the Hertzsprung-Russell diagram. Here we take advantage of this discovery to enlarge the IR-faint white dwarf sample threefold. We expand our selection to the entire Pan-STARRS survey footprint as well as the Montreal White Dwarf Database 100 pc sample, and identify 37 candidates with strong flux deficits in the optical. We present follow-up Gemini optical spectroscopy of 30 of these systems, and confirm all of them as IR-faint white dwarfs. We identify an additional set of 33 objects as candidates based on their colors and magnitudes. We present a detailed model atmosphere analysis of all 70 newly identified IR-faint white dwarfs together with 35 previously known objects reported in the literature. We discuss the physics of model atmospheres and show that the key physical ingredient missing in our previous generation of model atmospheres was the high-density correction to the He-minus free-free absorption coefficient. With new model atmospheres calculated for the purpose of this analysis, we now obtain significantly higher effective temperatures and larger stellar masses for these IR-faint white dwarfs than the Teff and M values reported in previous analyses, thus solving a two decade old problem. In particular, we identify in our sample a group of ultramassive white dwarfs in the Debye cooling phase with stellar parameters never measured before.


INTRODUCTION
The Hertzsprung-Russel (H-R) diagram for white dwarfs obtained by the Gaia mission (Gaia Collaboration et al. 2018) revealed a wealth of information and helped to identify particularly interesting structures. Most noteworthy was the so-called Q-branch, a nearly horizontal sequence composed mostly of massive DA and DQ stars, which has successfully been interpreted as evidence for crystallization (Tremblay et al. 2019), when the release of latent heat and chemical fractionation decrease the cooling rate of a white dwarf, leading to a pile up of objects in a color-magnitude diagram.
Also of interest was the B-sequence in the H-R diagram, a bifurcation between the non-DA and DA white dwarfs in the range 0.0 < (G BP − G RP ) < 0.8 that could not be reproduced by synthetic colors from model atmospheres with pure helium compositions and normal mass. To match the observed sequence, one needed to invoke either problems with the physics of pure helium model atmospheres, or a higher than av-erage mass for the non-DA sequence. Bergeron et al. (2019) indeed showed that cool (T eff 10, 000 K) non-DA white dwarfs had significantly larger than average masses (M ∼ 0.7 − 0.8 M ) when analyzed with the photometric technique using pure helium atmospheres. However, they also convincingly demonstrated that normal masses could be inferred (see their Figure 11) when using helium-rich models containing a small trace of hydrogen (H/He = 10 −5 ; see also Figure 4 of Bédard et al. 2022 for a prediction in the color-magnitude diagram).
Another previously undetected feature in the colormagnitude diagram was reported by Kilic et al. (2020) who identified an almost horizontal branch with M g ∼ 15.5 in the M g vs (g − z) Pan-STARRS color-magnitude diagram (see their Figure 21) composed of white dwarfs with optical and near-infrared flux deficits. This infrared flux deficiency is most likely caused by collisioninduced absorption by molecular hydrogen due to collisions with helium (H 2 -He CIA).  Hollands et al. (2021) At the low temperatures and high densities of cool white dwarf atmospheres, CIA becomes the dominant source of opacity (e.g., Borysow et al. 2001). In pure hydrogen atmosphere white dwarfs, CIA appears at effective temperatures below 4000 K. Hence, white dwarfs that show significant near-infrared flux deficits have traditionally been classified as ultracool white dwarfs. However, cool helium-rich white dwarfs have lower opacities and higher atmospheric pressures, and the infrared flux deficiency due to H 2 -He CIA starts to dominate at higher effective temperatures. Bergeron & Leggett (2002) presented a model atmosphere analysis of the prototype of this class, LHS 3250, and demonstrated that the pure hydrogen atmosphere models fail to reproduce the observed energy distribution. Even though none of their models provided perfect fits, they concluded that LHS 3250 is better explained with a helium-rich composition (see also Gianninas et al. 2015). Fortunately, a few of the white dwarfs with near-infrared flux deficits are DZ white dwarfs. Blouin et al. (2018b) simultaneously fit the metal ab-sorption lines and the CIA in one of these DZ white dwarfs, WD 0801+228 (see also Blouin et al. 2019b), and found an excellent fit to both the spectroscopic and photometric data with a model that has T eff = 4970 K and log H/He = −1.6, confirming a mixed atmosphere composition and a relatively warm temperature. Hence, it is more appropriate to classify these objects as "IRfaint". And as we shall demonstrate below, the term "ultracool" should definitely be abandoned to describe these white dwarfs showing strong CIA absorption.
The fact that the observed IR-faint white dwarf sequence is so tight indicates that these stars likely have similar hydrogen abundances, and that this must be the cooling sequence of white dwarfs with mixed H/He atmospheres. Bergeron et al. (2019) and Bédard et al. (2022) showed that many of the warmer DC stars (above 6000 K) have trace amounts of hydrogen in their atmospheres. Hence, an IR-faint white dwarf sequence due to CIA seems unavoidable, but was detected only recently.
We now have a fantastic opportunity to explore this sequence thanks to large scale photometric and astro-metric surveys. Here we take advantage of the Gaia DR2 and EDR3 astrometry to expand the IR-faint white dwarf sample to the entire Pan-STARRS footprint and the Montreal White Dwarf Database (MWDD) 100 pc sample . We discuss our sample selection in Section 2, and present follow-up optical spectroscopy of 30 candidates in Section 3. We describe the theoretical framework used in our analysis in Section 4, including new model atmospheres and evolutionary models. We then provide model atmosphere fits for all confirmed and candidate IR-faint white dwarfs in Section 5. We discuss the overall sample properties in Section 6, and conclude in Section 7.

SAMPLE SELECTION
Since the discovery of the first IR-faint white dwarfs, LHS 1126 and LHS 3250 (Bergeron et al. 1994;Harris et al. 1999), we have been able to find only 35 such stars in the past two decades. Table 1 presents this list. Many of these white dwarfs were found serendipitously by the Sloan Digital Sky Survey (SDSS) spectroscopy (Gates et al. 2004;Harris et al. 2008), which targeted them due to their unusually blue colors. Most show flux deficits in the redder SDSS bands, but a significant fraction show significant flux deficits only in the near-infrared range beyond 1 µm (see for example Kilic et al. 2010;Gianninas et al. 2015). All 35 objects in Table 1 show significant flux deficits compared to the best-fitting pure hydrogen and pure helium atmosphere models. However, this list excludes white dwarfs like SDSS J004506.40-060825.7 (Oppenheimer et al. 2001), SDSS J110217.48+411315.4 , and SDSS J145239.00+452238.3 , where the spectral energy distributions do not show significant flux deficits compared to the pure hydrogen or pure helium models, according to our own analysis of these objects.
The previous surveys for IR-faint white dwarfs have been mostly limited to the SDSS footprint (Gates et al. 2004;Harris et al. 2008;Kilic et al. 2020). To take advantage of the significantly larger sky coverage of the Pan-STARRS 3π survey, we selected all white dwarf candidates with 5σ significant parallax measurements from Gaia DR2 that are also in the Pan-STARRS footprint. We used the same color-magnitude and quality cuts as in Kilic et al. (2020) to select a relatively clean white dwarf sample. Figure 1 shows a color-magnitude diagram for the Pan-STARRS + Gaia DR2 sample with /σ ≥ 5 and within 120 pc. Previously known IR-faint white dwarfs are marked by magenta dots. We also show for reference the same cooling sequences for pure hydrogen, pure helium, and mixed H/He atmosphere white dwarfs as in Figure 1. Color-magnitude diagram of the Pan-STARRS white dwarfs within 120 pc and with 5σ parallax measurements from Gaia DR2. Previously known IR-faint white dwarfs are marked by magenta circles. Newly identified IRfaint white dwarfs that were targeted at Gemini are marked by red circles. Also shown are same cooling sequences as those displayed in Figure 21 of Kilic et al. (2020) for various masses and atmospheric compositions, as indicated in the figure. Figure 21 of Kilic et al. (2020). Because IR-faint white dwarfs have M g ranging from 15 to 17 mag, and Gaia's limiting magnitude is ≈ 21 mag, expanding our search beyond the local ∼100 pc volume is currently not feasible.
We selected candidates in the H-R diagram along the nearly horizontal IR-faint white dwarf sequence identified by Kilic et al. (2020), ±1 mag in M g , and fit their Pan-STARRS photometry and Gaia DR2 parallaxes with pure hydrogen and pure helium models to identify objects with significant flux deficits compared to these models. Some of these candidates appear to be consistent with being massive pure hydrogen atmosphere white dwarfs. However, there are many that clearly bear the signature of IR-faint white dwarfs; they are significantly fainter than expected from pure hydro-gen or pure helium atmosphere models in the redder Pan-STARRS bands.
In addition to the Pan-STARRS footprint, we also searched for white dwarfs with strong flux deficits in the MWDD 100 pc sample using Gaia colors. The MWDD selection is based on Gaia DR2 (Gaia Collaboration et al. 2018), and includes all candidates with 10σ significant parallax ( ), G BP and G RP photometry, and + σ > 10 mas. Non-Gaussian outliers in color and absolute magnitude are removed using the recommendations from Lindegren & Hernández (2018), and a cut in Gaia color and absolute magnitude is used to select the white dwarf candidates. We identify several additional IR-faint white dwarf candidates outside of the Pan-STARRS footprint. However, Gaia photometry is not sufficient to reliably confirm them as IR-faint white dwarfs. Two of these objects have grizy photometry available from the Dark Energy Survey (Dark Energy Survey Collaboration et al. 2016), and clearly show significant flux deficits in the izy bands.
In total, we identify 37 IR-faint white dwarf candidates with strong flux deficits in the optical. We mark these with red circles in Figure 1 and present their Gaia source IDs, astrometry, and Pan-STARRS or Dark Energy Survey grizy photometry in Table 2. Four of these objects are missing from Figure 1 because they are either outside of the Pan-STARRS footprint (the two DES objects) or they lack z-band photometry. Only one of these newly identified candidates, SDSS J003908.33+303538.9, has spectroscopy available in the SDSS, confirming it to be a DC white dwarf. A comparison with Figure 21 of Kilic et al. (2020) shows the substantial increase in the number of IR-faint white dwarf candidates based on our selection.

OBSERVATIONAL DATA
We tested the efficiency of our sample selection using the Fast Turnaround observing mode on Gemini. We obtained optical spectroscopy for twelve targets using the 8-m Gemini North telescope equipped with the Gemini Multi-Object Spectrograph (GMOS) as part of the program GN-2020B-FT-104. We used the B600 grating and a 1 slit, and binned the CCD by 4×4. We initially centered the grating at 5300Å and observed J0414+0309 and J1922+0233 in this setup, which provided spectra over the wavelength range 3730-6920Å, with a resolution of 2.05Å per pixel.
IR-faint white dwarfs have spectral energy distributions that usually peak beyond 5000Å. Our initial observations showed that the signal-to-noise ratio of the spectra in the blue, below 5000Å, was poor and the blue portion of the spectra did not provide any meaningful Figure 2. Gemini GMOS spectroscopy of 30 newly identified IR-faint white dwarfs. Objects are plotted based on their g − r colors, increasing from top left to the bottom right. This sample essentially doubles the number of spectroscopically confirmed IR-faint white dwarfs known. All but one of these objects, J1922+0233, are confirmed to be DC white dwarfs with featureless spectra.
constraints on the spectral types of the first two objects observed. To sample the peak of the spectral energy distribution, we changed our observing setup and shifted the central wavelength to 6350Å for the rest of the targets. This setup provided spectra over the wavelength range 4770-7980Å (J1922+0233, mentioned above, was reobserved with this alternative setup). We reduced the Gemini data using the gmos package under IRAF.
Our Fast Turnaround program was a success, confirming all of the observed targets as IR-faint white dwarfs. We observed 18 additional objects as part of the queue programs GN-2021A-Q-203 and GS-2021A-Q-300 using the same setup as above, with a central wavelength of 6350Å. Three more targets were included in the queue, but did not get observed. In total, we obtained spectra for 30 candidates at Gemini.  Figure 2 shows the spectra for all 30 candidates observed at Gemini 1 . The spectra are organized in increasing g − r color from the top left to the bottom right. All but one of these objects are confirmed to be DC white dwarfs with featureless spectra. The exception is J1922+0233, which is a DZ white dwarf with a strong Na D absorption feature (see also Tremblay et al. 2020).
The Gemini spectra confirm the unusual spectral energy distributions of all 30 targets; they peak between 5000 and 7000Å, and show significant absorption in the red. For comparison, the spectral energy distribution of the prototype IR-faint white dwarf LHS 3250 peaks at around 6000Å. About half of the targets in Figure  2 have spectral energy distributions that are even more extreme than LHS 3250.

Photometric Technique
We use the photometric technique as detailed in Bergeron et al. (2019), and follow the same approach and use the SDSS u, Pan-STARRS grizy photometry, and Gaia EDR3 parallaxes in our analysis. If available, we also supplement these data with near-infrared photometry from the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007) and the VISTA Hemisphere Survey (VHS, McMahon et al. 2021). Because CIA dominates in the near-infrared, the majority of our targets are too faint to be detected in the UKIDSS and VHS. Only three targets have near-infrared J-band photometry available, but no H-or K-band data.
We convert the observed magnitudes into average fluxes using the appropriate zero points, and compare with the average synthetic fluxes calculated from model atmospheres with the appropriate chemical composition. Since our sample is restricted to ∼100 pc, we do not correct for reddening. We fit for the effective temperature and the solid angle, π(R/D) 2 , where R is the radius of the star and D is its distance. Given the distance measurements from Gaia, we constrain the radius of the star directly, and its mass based on evolutionary models. The details of our fitting method are further discussed in Gianninas et al. (2015), Bergeron et al. (2019), and Kilic et al. (2020). We use detailed mass-radius relations for CO-core (X C = X O = 0.5) white dwarfs with masses in the range M = 0.2 − 1.3 M (Bédard et al. 2020). For the purpose of this analysis, we also calculated a few He-core models with STELUM (Bédard et al. 2022). These Hecore sequences were started from artificial static white dwarf models, but the initial structures are quickly "forgotten" and thus have no impact on the mass-radius relation at the low effective temperatures of interest here. However, this means that our He-core sequences do not provide reliable cooling times, which can only be obtained through a detailed modeling of pre-white dwarf evolutionary phases (Althaus et al. 2013;Istrate et al. 2016). Figure 3 shows the evolutionary models for He-core (red) and CO-core (blue) white dwarfs. As expected, a He-core white dwarf is larger and more luminous (which corresponds to a lower surface gravity) than a CO-core white dwarf at a given mass. The dashed line shows the lowest surface gravity (log g = 7) model available in our model atmosphere grid.

Model Atmospheres
The first grid of model atmospheres we used is that described in the analysis of IR-faint white dwarfs by Gianninas et al. (2015, also used in Kilic et al. 2020, based on the original calculations of Bergeron et al. (1995), but with improvements to the H 2 -He CIA calculations from Jørgensen et al. (2000), including the density correction from Hare & Welsh (1958) 2 . These models assume an ideal gas equation of state. Here we have significantly extended this model grid in log H/He = −5.0, −4.0 (0.5) −1.0 (1.0) +2.0, where the numbers in parentheses indicate the step size, and log g = 7.0 (0.5) 9.0.
As was the case for LHS 3250 analyzed in detail by Bergeron & Leggett (2002), or similar IR-faint white dwarfs reported by Kilic et al. (2020), we measured extremely low temperatures (T eff < ∼ 4000 K) for most objects in our sample, but more importantly, the masses we inferred were excessively small (M ∼ 0.15 to 0.35 M ). The reason for such low masses can be understood by looking at Figure 4 where we show the same colormagnitude diagram as in Figure 1, but with evolutionary sequences for various masses, core compositions, and log H/He ratios. The maximum CIA absorption in these models occurs around log H/He ∼ −2.5, where the models reach their maximum luminosity in the optical. Blouin et al. (2018b) discussed the reason for this maximum: a large H/He ratio means that there is more H 2 to cause CIA, but a very small H/He ratio means the density is higher in a helium-dominated atmosphere, which also leads to stronger CIA. At log H/He ∼ −2.5, however, the predicted colors at 0.6 M fall short of matching the observed sequence of IR-faint white dwarfs in Figure 4 by a full magnitude. Even at a mass of 0.2 M , the predicted CO-core sequence is still not luminous enough. Only with He-core models with a mass as low as 0.15 M are we able to match the observed IR-faint sequence.
Another feature of interest in this figure is the almost perfect overlap of the 0.6 M , CO-core models with log H/He = −1 and −4. Such a degeneracy was also reported by Blouin et al. (2018b) who found two solutions at log H/He = −1.6 and −3.5 for the DZ white dwarf WD 0801+228. Even though the latter solution produced a slightly better fit to the photometry, it was incompatible with the observed metal lines in this star and therefore ruled out. We come back to this degeneracy issue below.
In addition to the extreme stellar parameters we measure using these models, the peak of the energy distribution is always predicted too narrow, even though our solutions provide a reasonable match to the overall distribution, as shown for instance in Figure 20 of Kilic et al. 2020. As concluded by Bergeron & Leggett (2002), although there is little doubt that these IR-faint white dwarfs have helium-rich compositions, the prob-  Figure 1 along with the evolutionary sequences for pure hydrogen atmosphere CO-core white dwarfs with 0.6 M shown in red, while green, blue, and cyan lines show the evolutionary sequences for CO-and He-core white dwarfs with a variety of masses and log H/He ratios using our old grid of model atmospheres (see text). Each sequence is labeled based on its core composition, mass, and atmospheric composition. He/0.15/−2.5 means a He-core white dwarf with M = 0.15 M and log H/He = −2.5. lem with these extreme stellar parameters probably lie in the physics included in our model atmospheres, which is either inadequate or incomplete, the most obvious being the non-ideal effects of the equation-of-state at the high atmospheric pressures that characterize these heliumrich atmospheres.
With this idea in mind, Blouin et al. (2018a) presented a new generation of cool white dwarf atmosphere models that include the improved H 2 -He CIA opacity profiles from Abel et al. (2012) with a high-density correction at λ 1.5 µm based on the ab initio molecular dynamics simulations from Blouin et al. (2017), as well as other physical improvements. As we shall see below, the most relevant of these physical improvements, in the context of our study, is the correction to the He − free-free absorption coefficient by Iglesias et al. (2002). The differences between the old and new models are further dis-cussed in detail by Blouin et al. (2018a). Unfortunately, their model grid is incomplete at very low effective temperatures (T eff < ∼ 3000 K) due to numerical issues (convergence of the temperature structure), and could not be used to analyze our sample of (mostly cool) IR-faint white dwarfs.
In order to provide a qualitative and quantitative comparison between both model grids, we first compared the best-fit parameters in the temperature range where they overlap near T eff ∼ 4800 K. We found small differences in the best-fitting parameters from the two model grids, but the fits are qualitatively the same. We then performed an additional comparison by calculating a small set of models similar to those of Blouin et al. (2018a) but this time at much lower temperatures around T eff = 3000 K and log H/He = −3, and compared those with our old models. In this case, the models showed significant differences, with the peak of the energy distribution in the new models considerably shifted to the red with respect to our old models.
In fact, it was impossible to adjust the parameters of the old models to match even approximately the new models. Given that the old models provide fairly decent fits to the observed spectral energy distributions, at least qualitatively, we are thus forced to conclude that the new models would fail to match the photometric observations of the coolest IR-faint white dwarfs in our sample.
We traced back the problem to the use of the improved H 2 -He CIA opacity profiles from Abel et al. (2012) used in the new models, while our old models rely on the earlier calculations from Jørgensen et al. (2000). A comparison between two models calculated with these two sets of CIA opacities, everything else being equal, is displayed in Figure 5. The energy distribution with the CIA profiles from Jorgensen et al. peaks at much shorter wavelengths, in better agreement with the observations, no matter whether we try to vary the effective temperature of the models relying on the calculations from Abel et al. It is worth mentioning that a large difference at short wavelengths (λ < 2 µm) between these two sets of opacity calculations was noted before by Abel et al. (2012, see their Figure 6). It is difficult to interpret our result any further given that the calculations from Abel et al. represent a significant improvement over Jorgensen et al., and a deeper investigation of this discrepancy is clearly outside the scope of this paper.
Despite this situation, we decided to continue our model comparison but this time by using the CIA profiles from Jørgensen et al. (2000) throughout. We calculated a small grid of models using the same atmosphere code as in Blouin et al. (2018a), except for the Figure 5. Energy distribution of a model from Blouin et al. (2018a), which relies on the improved H2-He CIA opacity profiles from Abel et al. (2012, in blue), compared to a similar model but calculated with the earlier CIA profiles from Jørgensen et al. (2000, in red). Note that in this last case, we also include the density correction from Hare & Welsh (1958).
adopted CIA opacity, and fitted J1238+2633, a fairly warm IR-faint white dwarf in our extended sample (see Section 5.4). The comparison of our fits using both model grids is displayed in Figure 6. While both sets of models fail to provide perfect fits to the observed energy distribution, the most striking feature is that the derived stellar parameters are drastically different. In particular, our old models yield an effective temperature that is almost 1000 K cooler, and a stellar mass that is ∼0.3 M smaller than the solution obtained with this test model grid, entirely consistent with the odd results we obtain from the photometric analysis of our IR-faint white dwarf sample.
In order to understand which of the physical effects is responsible for the observed discrepancy between our solutions, we compare in Figure 7 various spectral energy distributions with different assumptions in the input physics. The black (labeled Bergeron) and green (labeled Blouin) lines correspond respectively to our old models and those of Blouin et al. (2018a). As discussed above, even though these two energy distributions appear qualitatively similar, the differences are much more important at lower temperatures. In the model shown in cyan, we used the same input physics as in Blouin et al., but we rely on the H 2 -He CIA opacity profiles from Jørgensen et al. (2000) instead of Abel et al. (2012); the effect here is significant (see also Figure 5). Then, in the model shown in red, we removed the correction to the He − free-free absorption coefficient from Iglesias et al. (2002); as can be seen here, this correction has a major effect on the predicted energy distribution. And finally, in the model shown in blue, we reverted back to the ideal gas equation of state; in the particular physical regime explored here, the nonideal effects due to the equation of state are small but non-negligible compared to the changes in H 2 -He CIA opacity and He − correction. We can also see that this last model is practically identical to our old model (shown in black). These comparisons indicate that for all models calculated with the H 2 -He CIA opacity profiles from Jørgensen et al. (2000), the most important effect is the inclusion of the correction to the He − free-free absorption coefficient.
This particular behavior can be explained in terms of a competition between the H 2 -He CIA opacity and the He − free-free opacity, as illustrated in Figure 8, where we show the contributions of each opacity source at the photosphere of our test model, with and without the correction to the He − free-free opacity taken into account. The correction has the effect of reducing significantly the relative contribution of the He − free-free opacity with respect to that from the H 2 -He CIA opacity, as can also be appreciated by comparing the cyan and red  (Bergeron) and green (Blouin) lines correspond respectively to our old models (see text) and those of Blouin et al. (2018a), while the remaining three energy distributions are similar to Blouin's models but with the H2-He CIA opacity profiles from Jørgensen et al. (2000) instead of Abel et al. (2012). In addition, we have removed the correction to the He − free-free absorption coefficient from Iglesias et al. (2002) in the model shown in red (no corr), and used the ideal gas equation-of-state in the model shown in blue (old EOS). models in Figure 7. In addition, the overall total opacity is reduced, resulting in larger densities that affect the overall atmospheric structure. Given these results, we have thus decided to recalculate our original model grid described at the beginning of this subsection but by properly including this high-density correction factor to the He − free-free absorption coefficient.
We compare in Figure 9 the fits to J2305+3922, a typical IR-faint white dwarf in our sample, obtained with our old model grid to those achieved with our revised models (a full comparison for our complete sample will be discussed in Section 6.2). While the photometric fit with the old models is a poor match to the measured photometry -similar in quality to the fits displayed in Figure 20 of Kilic et al. (2020) -, the fit to the photometry (and spectroscopy as well) using our revised model grid is significantly superior. More importantly, the stellar parameters (given at the top of the figure) have changed drastically, going from an ultracool (T eff = 3146 K) and extremely low-mass (0.18 M ) IR-faint white dwarf, to a much hotter (4550 K) and massive (0.70 M ) white dwarf; the H/He abundance ratio remains essentially unchanged, however. The increase by more than 1400 K in temperature makes the object much more luminous, thus requiring a smaller stellar radius and larger mass. These revised models will have important consequences on the results of our analysis.
Note also that despite the fact that we now obtain much higher effective temperatures for our objects, well within the range of T eff values calculated by Blouin et al. (2018a), those models would still fail to match the photometry of the objects in our sample with strong IR-flux deficiencies because of the problem related to the use of the H 2 -He CIA opacity profiles from Abel et al. (2012), as discussed above.
We end this section by stressing that even though our revised models represent a significant improvement over our previous model grid, both qualitatively and quantitatively, they represent by no means the best models that could be achieved. First, we have neglected the nonideal effects in the equation of state. Second, instead of using the improved H 2 -He CIA opacity calculations of Abel et al. (2012), we relied on the more approximate calculations from Jørgensen et al. (2000), which for some reason provide a much better description of the observed energy distribution of the coolest IR-faint objects in our sample. Clearly, this deserves further investigation. Figure 10 shows the spectral energy distributions of the 37 newly identified IR-faint white dwarfs with strong Figure 9. Fits to the spectral energy distribution of the IRfaint white dwarf J2305+3922. This object was also spectroscopically confirmed at Gemini, and the spectrum is shown in cyan. Here we show two solutions, one with our original model grid (black), and the other solution using our revised model grid where the correction to the He − free-free absorption coefficient from Iglesias et al. (2002) has been taken into account (red).

IR-Faint White Dwarfs with Strong Flux Deficits
flux deficits in the optical. Error bars show the observed photometry and the red lines show the spectra, normalized at ∼6000Å. Optical spectra are available for 31 of these objects, and in most cases the spectra follow the photometric spectral energy distributions, revealing relatively narrow peaks with significant absorption in the red. Each object is labeled based on its Gaia Source ID, object name based on Gaia DR2 coordinates, and the photometry used in the fitting: ugrizy means SDSS u + Pan-STARRS grizy. Solid black lines and dots show the monochromatic fluxes and the synthetic photometry for the best-fit model, respectively. The parameters of this model are included in each panel. For reasons discussed in Section 6.1, with this particular model grid, we always found a unique and robust solution to all the IR-faint white dwarfs analyzed here and below. In particular, we did not find any degeneracy in our solutions similar to those reported in previous investigations (see Section 4.3).
In contrast with the fits obtained with our old model grid (see, e.g., Bergeron & Leggett 2002, Gianninas et al. 2015, Kilic et al. 2020, those displayed in Figure 10 show a remarkably good agreement with the optical and infrared photometry, as well as with the optical spectroscopy. In those cases where the predicted monochromatic fluxes do not match the observed spectrum, there is also a discrepancy between the spectrum and the measured photometry (see, e.g., J0231+3254, J0440−0414, etc.), suggesting that some of the spectra suffer from flux calibration issues. But in most cases, the agreement is excellent. While the solutions for these IR-faint white dwarfs obtained with our old model grid (not shown here) span a range of T eff = 1950 to 3370 K and M ∼ 0.15 to 0.35 M , the best fits with our revised models yield much higher temperatures well above 4000 K for most objects (see also Section 6.2), and significantly larger masses, with more than half of the objects in our sample in the range M ∼ 0.8 − 1.2 M . We discuss further the global properties of our sample in Section 6.
Since the objects displayed in Figure 10 represent extreme cases of white dwarfs with the strongest observed infrared flux deficiencies, we are forced to conclude that IR-faint white dwarfs are not ultracool afterall. Also, they do not require unusually large radii and low masses either, as was the case with the prototype LHS 3250 analyzed by Bergeron & Leggett (2002). Obviously, the key physical ingredient missing in our previous generation of model atmospheres was the correction to the He − free-free absorption coefficient described in Iglesias et al. (2002).

An IR-Faint DZ White Dwarf
J1922+0233 is the first IR-faint white dwarf discovered with an absorption feature and strong flux deficits in the optical bands (see Figure 10, and also McCleery et al. 2020;Tremblay et al. 2020). It shows a strong Na doublet feature, which can be used as an independent diagnostic of the H/He ratio. However, our attempts to fit both the sodium doublet and the CIA feature in J1922+0233 failed. Figure 11 shows a comparison between the observed sodium feature and a synthetic spectrum calculated using the unified Na profiles described in Blouin et al. (2019a) with the stellar parameters given in Figure 10 (T eff = 4436 K, M = 1.065 M , log H/He = −1.73), and by adjusting the Na abundance (only in the synthetic spectrum calculation) to match the depth of the observed Na D doublet. As can be seen in the figure, the inferred abundance of log Na/He = −9.7 predicts way too much broadening of the Na D doublet as a result of the high photospheric pressure encountered in He-rich atmospheres. The relatively narrow sodium feature in J1922+0233 can only be explained by lowering the density significantly, at least in the line-forming region. For instance, it is possible to achieve a good fit to the sodium feature by assuming a pure hydrogen atmosphere, where the photospheric density is two orders of magnitude lower than in our best-fit model, but in this case the predicted energy distribution is totally inconsistent with the observations, as hydrogen-dominated at- Figure 10. Fits to the spectral energy distributions of 37 new IR-faint white dwarfs with strong flux deficits in the optical. Black lines show the monochromatic fluxes for the best-fit model for each star, and dots show the synthetic photometry of those models in each filter. Thirty of these objects were spectroscopically confirmed at Gemini, and one (J0039+3035) has a spectrum in the SDSS. These spectra are shown in red. The remaining six targets currently lack optical spectroscopy. All 37 targets have spectral energy distributions that peak in the optical and display significant absorption in the near-infrared bands, and are best explained by mixed H/He atmospheres.  mospheres do not produce strong CIA features that are required to fit the photometry.
We know of other cool DZs with mixed H/He atmospheres and Na D features. For example, Blouin et al. (2019a) was able to fit successfully both the CIA and the Na D features in WD J2356−209 with a model that has T eff = 4040 ± 110 K and log H/He = −1.5 ± 0.2. Hence, it is surprising in the case of J1922+0233 to see simultaneously very strong CIA (which demands a low hydrogen abundance) and a narrow Na D feature (which demands a high hydrogen abundance). It is difficult to reconcile those two observations. It is also interesting that only one out of the 37 objects presented here is a DZ, whereas ∼30% of helium atmosphere white dwarfs in the 4000-5000 K range are DZs. Hence, J1922+0233 seems to be rather different than the general population of IR-faint white dwarfs with strong flux deficits in the optical. J1922+0233 may be an exception, rather than the rule, among the IR-faint white dwarf sequence.

Known IR-Faint White Dwarfs
Given our revised model atmospheres, we felt it was worth reanalyzing the previously known IR-faint white dwarfs published in the literature and summarized in Table 1. Our best fits to these 35 IR-faint white dwarfs are Figure 11. Na D region of J1922+0233. The atmospheric parameters are set to the values given in Figure 10, and the Na abundance of the synthetic spectrum is adjusted (log Na/He = −9.7) to match the depth of the observed Na D doublet. We find that the inferred He-rich composition leads to too much broadening of the Na D doublet. presented in Figure 20 (in Appendix). Again, the quality of these fits is excellent in most cases. For J1238+3502, we found that we could achieve a much better fit to this object (also displayed in Figure 20) by dropping the J magnitude reported by Kilic et al. (2010). However, since there is no reason to believe this measurement is erroneous, we provide here both solutions for this star.
Some of the previously known IR-faint white dwarfs in Figure 20 show significantly less absorption in the infrared than those displayed in Figure 10 as a result of their larger hydrogen abundances, and thus smaller H 2 -He CIA opacity (see, e.g., J0041−2221, J0309+0025, J0854+3503, etc.). Otherwise, they are found in the same temperature and mass range as those reported in our spectroscopic sample. Again, we defer the discussion of the global properties of this sample to Section 6. Also worth mentioning is the case of the prototype LHS 3250 (J1653+6253), for which we obtain T eff = 4993 K, M = 1.049 M , and log H/He = −2.74 (see the fit in Figure 20), while Bergeron & Leggett (2002) reported significantly different values of T eff = 3480 K, M = 0.23 M , and log H/He = −4.7, and a considerably worse fit (see their Figure 7). Our improved models have successfully solved this two decade old problem.

Additional IR-Faint White Dwarfs
Our follow-up spectroscopy presented in Section 5.1 specifically targeted objects with strong flux deficits in the Pan-STARRS bands (see Figure 1). However, this sample represents only the tip of the iceberg, and there are likely many other IR-faint white dwarfs hiding in the 100 pc sample. In order to identify those missing white dwarfs with weaker flux deficits in the optical, we first investigated the spectral energy distributions of all targets below the observed white dwarf sequence, and then expanded our search to all white dwarfs in the MWDD 100 pc sample.
We fitted the spectral energy distribution of each target with pure hydrogen and pure helium atmosphere models, and compared the resulting χ 2 to the bestfitting mixed H/He atmosphere model fits. We identified 210 candidates as potential IR-faint white dwarfs. We then cross-correlated this list with the UKIDSS and VHS, and found near-infrared photometry for 71 targets, and updated our model fits for those stars with both optical and near-infrared data. Four of these objects have follow-up spectroscopy available in the literature, and all four are confirmed to be DC white dwarfs by Tremblay et al. (2020). Figure 12 shows our model fits to one of these DC white dwarfs, J2150−0439. Note that because our search for these additional candidates was performed Figure 12. Fits to the spectral energy distribution of a newly identified IR-faint white dwarf with mild flux deficits in the Pan-STARRS bands; note that we rely here on our old model atmospheres. J2150−0439 is a spectroscopically confirmed DC white dwarf within 40 pc . The symbols are the same as in Figure 10. The red dots show the best-fitting pure hydrogen atmosphere model with parameters also highlighted in red. Near-infrared photometry reveals strong flux deficits compared to the pure H atmosphere models, confirming the IR-faint nature of this object.
prior to upgrading our model atmospheres, the fits displayed here are based on our old model atmosphere grid. Error bars show the observed photometry, and the red dots show the best-fitting pure hydrogen atmosphere model. The parameters for this model along with the reduced χ 2 of the fit are presented in red. IR-faint white dwarfs occupy the same region of the color-magnitude diagram as massive white dwarfs; the pure hydrogen atmosphere solution for J2150−0439 requires a mass of 1.2 M . However, this solution is clearly bad, with a reduced χ 2 = 659.
The solid black line and dots show the monochromatic fluxes and the synthetic photometry for the best-fitting mixed H/He atmosphere model, respectively. The bestfitting mixed atmosphere model parameters are included in the figure for comparison. Even though the reduced χ 2 value for our helium-dominated model fit is not great (but recall that we are using our old model grid here), it provides a much better fit than the pure hydrogen model. J2150−0439 shows relatively mild absorption in the optical compared to the IR-faint white dwarfs discussed in the previous section. However, near-infrared photometry clearly shows significant flux deficits compared to the pure hydrogen atmosphere model; it is undoubtedly IR-faint.
We identified 33 additional IR-faint candidates using this procedure, including four confirmed DC white dwarfs, where the available photometry clearly favors an IR-faint classification. These IR-faint white dwarfs are listed in Table 3 and our best fits using our revised model atmospheres are displayed in Figure 21 (in Appendix). One remarkable feature of this particular subsample is that all objects have more normal masses below 0.8 M , the reason of which is discussed in the next section. Also, we uncovered a few IR-faint white dwarfs with extremely low masses below ∼0.3 M (J0357−2606, J0406−0333, J1448+2935, and J2035+4054).
Even though there are other white dwarfs with mild flux deficits in the Pan-STARRS bands, many would require infrared photometry and optical spectroscopy to confirm their nature. In order to avoid mis-classifying objects based on noisy z-or y-band data, we adopted a conservative approach and only included objects where available optical and near-infrared photometry clearly shows a strong flux depression in the red filters. For example, 11 of the 30 candidates without follow-up spectroscopy have near-infrared photometry available, and they are unambiguously IR-faint. Similarly, 14 of these objects have SDSS u-band photometry available, which anchors the spectral energy distribution in the blue and helps identify significant absorption in the red. Followup spectroscopy of all of these objects would be helpful in confirming their nature. 6. DISCUSSION

Color Diagrams
We used Gaia astrometry and Pan-STARRS photometry to identify 37 IR-faint white dwarfs with strong flux deficits in the optical, and 33 additional white dwarfs with milder deficits. We have follow-up spectroscopy available for 30 of them from Gemini and five from the literature. The optical spectra are featureless for all but one of these white dwarfs. The spectra follow the photometric spectral energy distributions with significant absorption in the red, confirming the nature of these targets as IR-faint white dwarfs. Hence, these two samples increase the number of IR-faint white dwarfs from 35 to 105, thus tripling the sample size. Figure 13 shows the color-magnitude diagram of the MWDD 100 pc sample (black dots) along with evolutionary sequences described further below. The previously known IR-faint white dwarfs and the newly identified IR-faint white dwarfs with strong and mild absorption are marked by magenta, red, and green dots, respectively. The four spectroscopically confirmed DC white dwarfs with mild absorption (see Figure 12) are marked by yellow dots. We also show color-color diagrams of the spectroscopically confirmed white dwarfs Figure 13. Color-magnitude diagram of the MWDD 100 pc sample (black dots) along with the previously known IRfaint white dwarfs (magenta dots, Section 5.3) and the newly identified IR-faint white dwarfs with strong (red dots, Section 5.1) and mild (green dots, Section 5.4) CIA absorption. Yellow symbols mark the four spectroscopically confirmed DC white dwarfs with mild absorption. The red and blue curves show the cooling sequences for 0.6 M CO-core white dwarfs with pure H and pure He atmospheres, respectively. The solid cyan curves have the same mass but with log H/He = 0, −1, −2, and −3, from bottom to top, while the solid black curves are with log H/He = −3.5, −4, and in the 100 pc sample in various Pan-STARRS filters in Figure 14.
These two figures clearly show that the 33 additional IR-faint white dwarf candidates with mild optical flux deficits (green and yellow dots) are an extension of the IR-faint white dwarf sequence. The previously known 35 IR-faint white dwarfs (magenta), newly identified 37 with strong flux deficits (red), and the 33 objects with mild CIA in the optical bands (green and yellow) clearly form a sequence in the color-magnitude diagrams and  Figure 14. Pan-STARRS color-color diagrams for the spectroscopically confirmed white dwarfs in the Montreal White Dwarf Database 100 pc sample (black dots) along with the IR-faint white dwarfs discussed here (colored symbols). The symbols are the same as in Figure 13. also in the g − r versus r − i color-color diagram. Even though many of the mild CIA cases are difficult to identify based on the gri photometry alone, these objects are clearly outliers in their i − z and z − y colors. The IR-faint white dwarf sequence around M g ∼ 16 is rather tight, both in the color-magnitude diagrams and the g−r versus r − i color-color diagram, suggesting that these objects probably have similar masses and compositions. We also find a dozen objects or so with luminosities below this tight sequence; we come back to these objects later. Also indicated as a black dotted line in Figure 13 is Gaia's limiting magnitude at G ∼ 21 in our 0.6 M models. This limiting magnitude may partially explain why we observe a lack of continuous IR-faint sequences in the color-magnitude diagram. Even with our detailed analysis of the 100 pc sample, we have only scratched the surface in terms of the IR-faint white dwarfs. Near-infrared observations of 112 cool DC white dwarfs by Kilic et al. (2010) found eight previously unknown IR-faint white dwarfs that show significant absorption in the H-and K-bands or only in the K-band. As the name indicates, given that these objects are relatively faint in the infrared, they are usually missing in the wide-field surveys like the UKIDSS and VHS. There are likely many more IR-faint white dwarfs within 100 pc that can only be identified through followup near-infrared observations. Also shown in Figure 13 are evolutionary sequences, calculated with our new model atmospheres, for 0.6 M CO-core white dwarfs with pure H, pure He, and mixed H/He compositions (red, blue, and cyan/black lines). In the case of mixed H/He atmospheres, we show the effects on the predicted colors of varying the atmospheric composition (log H/He between −5 and 0) as well as the stellar mass (M = 0.2 M and 1.0 M at log H/He = −2). Our new models reach a maximum optical luminosity (M g ) around log H/He = −4.0 to −3.5 at 0.6 M depending on the temperature. Even though our old models displayed in Figure 4 reach their maximum luminosity at roughly the same H/He value, the predicted luminosities are well below the observed sequence of IRfaint white dwarfs, even with a mass as low as 0.2 M with a CO-core. For this reason, our fitting code with the old models always pick the solution where the CIA is strongest in the model grid (see, e.g., Figure 20 of Kilic et al. 2020), precisely because these models underestimate the CIA opacities compared to the observations, as discussed above. In contrast, some of our new models are now more luminous than the observed sequence of IR-faint white dwarfs, even at 0.6 M . This should lead to a larger spread in the best-fitted parameters, in particular the stellar mass and the hydrogen-to-helium abundance ratio. Also, the degeneracy observed in Figure 4 between the 0.6 M , CO-core models with log H/He = −1 and −4 is no longer observed in Figure  13 with our new model grid. This is the reason why we always find a unique solution to all the IR-faint white dwarfs in our sample. Finally, even though the results displayed in Figure 13 illustrate how a change in H/He can be compensated to some extent by a change in mass, this occurs only in this color-magnitude diagram. Indeed, while the luminosity in the optical is controlled by the radius, and thus the mass of the star, the CIA opacity, which dominates in the infrared, is controlled by the H/He ratio. One parameter cannot be simply compensated by the other when the overall energy distribution is considered. Table 4 presents the best-fit parameters for all 105 IR-faint white dwarfs in our sample, including the previously known and newly identified systems. Despite the fact that our revised models show significant improvements over our previous model grid, those models remain approximate, as mentioned above, and the numbers given in Table 4 should be taken with caution. As discussed in Section 5.3, we provide in Table 4 two solutions for J1238+3502, including one with the J magnitude omitted from the fit (see Figure 20). In both cases, the solutions appear way too cool for the inferred mass, which leads to cooling ages that are largely extrapolated.

Global Properties of the Sample
In Figure 15, we compare the physical parameters, T eff , M , and log H/He one against each other, using the same color symbols as in Figure 13. The IR-faint white dwarfs in our sample have a range in T eff between roughly 3000 K and 5600 K, although most objects are clustered around 4600 ± 200 K, reinforcing our conclusion that most IR-faint white dwarfs are not ultracool afterall. The masses span an interval between 0.2 and 1.3 M , as usually found for other mass distributions, although here we have a strong excess of massive (M > 0.8 M ) white dwarfs. Note that all these massive objects consist of previously known IR-faint white dwarfs (magenta) and newly IR-faint white dwarfs with strong CIA absorption (red) identified in our spectroscopic sample. In contrast, those with mild IR absorption (green) have normal masses around ∼ 0.6 M , indicating that more normal mass white dwarfs with their lower photospheric densities produce, in general, less CIA absorption, and they are thus more difficult to identify. This interpretation is consistent with their location in Figure 13, where they lie much closer to the main white dwarf sequence. Incidentally, these white dwarfs with mild CIA absorption tend to have similar hydrogen abundances around log H/He ∼ −2.5. The fact that we find normal masses for these IR-faint white dwarfs gives us confidence in our model atmospheres, at least at these photospheric densities.   The bottom panel of Figure 15 shows that the massive white dwarfs identified in the upper panel also tend to have the smallest hydrogen abundances around log H/He ∼ −3.5. In contrast, the IR-faint white dwarfs with mild CIA absorption and normal masses (green symbols) are characterized by hydrogen abundances around log H/He ∼ −2.5. This trend of log H/He with mass is most likely related to the origin, and thus the progenitors, of these IR-faint white dwarfs. If they represent the outcome of convectively mixed DA stars, our results indicate that more massive white dwarfs have smaller hydrogen abundances at the photosphere as a result of a larger dilution in the stellar envelope. We come back to this point in Section 6.4.
We also observe in Figure 15 IR-faint white dwarfs with equal amounts of hydrogen and helium in their atmospheres (log H/He ∼ 0), spanning a large range in stellar mass from 0.2 to 1.0 M . Two of these objects are the DZ white dwarfs J1824+1213 and J2317+1830, recently analyzed by Hollands et al. (2021). Given that Figure 15. Masses, effective temperatures, and H/He ratios for our IR-faint white dwarf sample. The symbols are the same as in Figure 13. the presence of metals in cool DZ stars is generally interpreted as evidence for accretion from external sources such as the interstellar medium, tidally disrupted asteroids, small planets, or even comets, the large hydrogen abundance measured in these stars may also be the result of accretion from such external sources, possibly containing water (see, e.g., Gentile Fusillo et al. 2017). The large hydrogen abundances in all other objects in our sample most likely have the same origin, even if they show no metal, because heavy elements eventually diffuse at the bottom of the mixed H/He convection zone, while hydrogen always remain within the stellar envelope. Interestingly, we obtain similar stellar parameters for these two DZ white dwarfs if we use the model atmospheres from Blouin et al. (2018a), because in such hydrogen-enriched atmospheres, the dominant source of opacity in the infrared is the H 2 -H 2 CIA opacity (instead of H 2 -He CIA), and both sets of model atmospheres in this case rely on the same opacity calculations.
For completeness, we compare in Figure 16 the masses and effective temperatures obtained with the original and revised models; the changes in log H/He values are not as important (see, e.g., Figure 9) and not discussed any further. As was described previously in the case of J2305+3922 displayed in Figure 9, both the M and T eff values using our revised models have increased significantly and systematically with respect to the values obtained with the original model grid, although the differences are larger in the most massive objects where the high-density correction to the He − free-free absorption coefficient is more important. Note also that with the original model grid, the extremely low masses inferred at 0.15 M represent only upper limits since this corresponds to the lowest mass value of our cooling sequences assuming a He-core (see Figure 3).

Kinematics
Kinematics can provide additional insights into the nature of IR-faint white dwarfs. Figure 17 plots the distribution of Galactic U , V , and W velocity components for the IR-faint white dwarf sample, along with the 1σ (dashed) and 2σ (dotted) velocity ellipsoid values for the thick-disk and halo (Chiba & Beers 2000). The thick disk velocity ellipsoid is (σ U , σ V , σ W ) = (46, 50, 35) km s −1 . Since we do not have any radial velocity constraints on the IR-faint white dwarfs, we assume zero radial velocity. The tangential velocities clearly center on the disk, except for three objects, WD 0343+247 (J0346+2455), LHS 2139 (J0925+0018), and J1824+1213, which are clearly halo stars.
The 2D velocity dispersion of the 105 stars, 60 km s −1 , is practically identical to the Chiba & Beers (2000) thick disk value. This motivated us to draw hypothetical radial velocities from the radial component of the thick disk velocity ellipsoid, 10,000 times for each star. The mean radial velocity of each star remains zero in this approach, but the dispersion (44 km s −1 , on average) provides a quantitative measure of how plausible radial velocities may impact the kinematics. For completeness, we also re-draw the measured proper motions and parallaxes using the full Gaia covariance matrix.
The error bars in Figure 17 show the dispersion of the U V W velocity components assuming the stars have thick disk radial velocities. Note that radial velocities project into different U V W velocity components depending on a star's location on the sky (i.e. at the pole, or anti-center, etc). Thick disk radial velocities inflate the velocity of any thin disk star, of course, but any thin disk star will remain consistent with the disk. It is possible that halo stars may remain hidden in the sample. Still, the overall sample remains clustered around the thick disk velocity ellipsoid, and we conclude that the IR-faint white dwarfs are a thick disk population on the basis of their kinematics. This is consistent with our relatively large white dwarf cooling age estimates presented in Table 4.

On the Origin of IR-faint White Dwarfs
We show in Figure 18 the location of the 105 IRfaint white dwarfs in our sample in a M vs T eff diagram, together with the spectroscopically confirmed white dwarfs in the 100 pc sample and the SDSS footprint from Kilic et al. (2020, see their Figure 18). Note that in Kilic et al., the DQ and DZ white dwarfs have been analyzed using detailed model atmospheres appropriate for these stars, but they are treated simply as non-DAs (white circles) in Figure 18. So even though the Kilic et al. sample is restricted to the SDSS footprint and is thus considerably smaller than the large Gaia sample displayed in Figure 11 of Bergeron et al. (2019), the analysis from Bergeron et al. relies on pure H and pure He model atmospheres only, and their results are therefore considerably less accurate.
At the scale used in Figure 18, the IR-faint white dwarfs are found in a rather narrow range of effective temperatures. The IR-faint objects with normal masses around 0.6 M appear as a natural extension of the Figure 18. Stellar masses as a function of effective temperature for all spectroscopically confirmed white dwarfs in the 100 pc sample and the SDSS footprint from Kilic et al. (2020); DA and non-DA white dwarfs are shown by red and white circles, respectively. Our sample of 105 IR-faint white dwarfs is shown by green circles. Solid curves are theoretical isochrones (white dwarf cooling age only), labeled in units of 10 9 yr, obtained from the cooling sequences of Bédard et al. (2020) with C/O-core compositions, q(He) = 10 −2 , and q(H) = 10 −4 . The lower blue solid curve indicates the onset of crystallization at the center of evolving models, and the upper one indicates the locations where 80% of the total mass has solidified. The dashed curve indicates the onset of convective coupling, while the dotted curve corresponds to the transition between the classical to the quantum regime in the ionic plasma (see text).
DA and non-DA white dwarfs found at higher temperatures. More striking in this figure is the excess of massive white dwarfs in the IR-faint population, which appears even more extreme here. There is also an apparent gap in T eff between the coolest massive DA stars in Figure 18 and the massive IR-faint white dwarfs. We have to be careful not to overinterpret these results, however, because there are various selection effects at play with these samples. First, the sample from Kilic et al. (2020) is restricted to the SDSS footprint, while ours covers the Pan-STARRS 3π survey. Moreover, the objects displayed in Figure 18 do not include white dwarf candidates without spectroscopic confirmation, many of which are still present in the 100 pc sample, and even more so if they are massive, with smaller radii, and thus lower luminosities. Also, as mentioned in Section 2, although we restricted our analysis to a distance of 100 pc, we performed a deeper search of all white dwarf candidates with 5σ significant parallax measurements from Gaia DR2. Hence the most extreme IR-faint objects in our sample, which also happen to be more massive than average, are most likely over represented in Figure 18 with respect to the rest of the white dwarfs displayed here.
Nevertheless, the massive IR-faint white dwarfs in Figure 18 occupy a region the M -T eff plane where few or no objects have ever been reported 3 . Obviously there are none in the Kilic et al. sample. All cool, and massive white dwarfs in this particular region of the M vs T eff diagram have mixed H/He atmospheres and show strong IR-flux deficiencies. If we interpret these objects as the result of convectively mixed DA stars, we encounter one obvious problem related to the measured values of the hydrogen abundance in these massive IRfaint white dwarfs. Figure 19 shows the variation of the H/He ratio as a function of T eff for DA models with 0.6 and 1.0 M and with various hydrogen layer masses (log q(H) = log M H /M ). These envelope calculations are similar to those of Rolland et al. (2018). Upon mixing, the H/He ratio goes from infinity down to a value set by the size of the mixed H/He convection zone after mixing. The H/He ratio then continues to vary with decreasing effective temperature as a result of changes in the size of the mixed H/He convection zone. Our results indicate that in a 1.0 M white dwarf at T eff = 4800 K with a photospheric value of log H/He ∼ −3.5 -typical of the massive IR-faint white dwarfs in our sample -, the total hydrogen mass is log q(H) ∼ −10.4. The DA progenitor with such a small hydrogen layer mass would mix at a temperature of T eff = 9200 K and cool off as a non-DA star. The problem then of course is that most, if not all, massive white dwarfs in Figure 18 are DA stars (see also Figure 11 of Bergeron et al. 2019). Given that there are remaining uncertainties in our model atmosphere calculations, in particular at the high densities encountered in these massive white dwarfs, it is possible that the hydrogen-to-helium abundance ratios are underestimated in our analysis. For the same reasons, it would be precarious at this point to estimate the total amount of hydrogen in individual objects until better models become available.
In contrast, the IR-faint white dwarfs with mild CIA absorption and normal masses (green symbols in Figure  15), and thus lower atmospheric pressure, are characterized by hydrogen abundances around log H/He ∼ −2.5, indicating that their DA progenitors had thicker hydrogen layers, and thus mixed at much lower effective temperatures, around T eff ∼ 7500 K according to the results of Figure 19. And in this case, there are plenty of non-DA stars with normal masses that could be the immediate progenitors of these IR-faint white dwarfs.
Also shown in Figure 18 are the regions in the M -T eff plane where crystallization occurs: the lower solid blue curve marks the region where crystallization starts at the center of the white dwarf, while the upper solid blue curve marks the limit of 80% crystallization of the core resulting from the solidification front moving upward in the star with further cooling. The dashed curve marks the onset of convective coupling, when the convection zone first reaches into the degenerate interior, where all of the thermal energy of the star resides. As discussed in more detail in Bergeron et al. (2019), the onset of convective coupling is weakly mass-dependent, and clearly interacts with the manifestations of crystallization, as illustrated by the behavior of the curved isochrones at high mass and low temperatures in Figure  18. The massive IR-faint white dwarfs in our sample located in this specific range of parameters have undoubt-edly entered the so-called Debye cooling phase, i.e., the rapidly cooling phase resulting from the transition, in the solid phase, from the classical regime where the specific heat of a solid is independent of temperature, to the quantum regime where it goes down from that constant value with decreasing temperature (as indicated by the blue dotted curve in Figure 18). Here, we defined this transition from the classical to the quantum regime by isolating the evolving model where the central temperature becomes equal to 1/20 the central Debye temperature (θ D /T = 20), as defined in Althaus et al. (2007). In the Debye cooling phase, the specific heat decreases quickly with cooling, which depletes rapidly the reservoir of thermal energy and produces the extreme increase of the cooling rate, as observed in Figure 18. It is thus perhaps surprising to have uncovered such a large number of massive IR-faint white dwarfs caught in this rapidly cooling phase.

CONCLUSION
We presented a detailed model atmosphere analysis of 105 IR-faint white dwarfs using the so-called photometric technique, based on Gaia astrometry combined with Pan-STARRS and near-infrared photometry. In particular, we identified 37 new IR-faint white dwarfs with strong flux deficits in the optical, 30 of which have follow-up Gemini optical spectroscopy, and 33 additional white dwarfs with milder deficits. We showed that the extremely low temperatures and stellar masses previously reported in the literature (see, e.g., the case of LHS 3250 analyzed by Bergeron & Leggett 2002) were related to inaccuracies in the model atmospheres. We convincingly demonstrated that the key physical ingredient missing in the previous generations of model atmospheres was the high-density correction to the He − free-free absorption coefficient described in Iglesias et al. (2002). Despite the significant improvements of our new models over our previous model grid, both qualitatively and quantitatively, we also stressed the importance that much work remains to be done on the theoretical front in terms of the equation of state and the H 2 -He CIA opacity calculations.
The two major uncertainties in IR-faint white dwarf models are (1) the atmospheric density, which depends on factors like continuum opacities and the pressure ionization of helium, and (2) the shape of the CIA opacities at a given density (Abel et al. 2012;Blouin et al. 2017). The only way to untangle these uncertainties and improve the physics in our models is to obtain spectra of IR-faint white dwarfs in the infrared. The molecular absorption features that dominate the infrared spectra of these objects have never been resolved. Yet they are predicted to exhibit structures that contain critical information on the physical conditions that characterize the atmospheres of IR-faint white dwarfs.
Recent first-principles calculations predict that the shape of the CIA features and the overall infrared spectral energy distribution are highly sensitive to the atmospheric density (see Figure 10 in Blouin et al. 2017). The James Webb Space Telescope provides an excellent opportunity to obtain such spectra, and for the first time constrain the shape of the CIA opacities in the high density environments of helium-dominated atmosphere white dwarfs. Laboratory experiments may also help constrain the opacity of helium under the conditions prevalent at the photospheres of IR-faint white dwarfs (McWilliams et al. 2015). In this work, we have shown how sensitive the inferred properties of IR-faint white dwarfs are to helium continuum opacities; experimental measurements of dense helium opacities would help confirm the nature of those stars.
Regardless of the problems with the model atmospheres, there are a few things that we can safely conclude: 1-The majority of the IR-faint white dwarfs form a tight sequence around M g = 16 in the color-magnitude diagram. Pure hydrogen atmosphere composition can safely be ruled out based on the location of these objects in the color-magnitude diagrams and the strength of the CIA observed in the optical bands. The observed sequence can be explained as the cooling sequence of mixed H/He atmospheres. The majority of the warm DC white dwarfs with T eff > 6000 K require trace amounts of hydrogen to explain their location in the Gaia color-magnitude diagrams (Bergeron et al. 2019), hence an IR-faint white dwarf sequence due to H 2 -He CIA seems unavoidable.
2-Not all white dwarfs become IR-faint. In fact, IRfaint white dwarfs represent only a small fraction of the overall cool white dwarf population in the solar neighborhood. This indicates that most DC white dwarfs below 5000 K have hydrogen-dominated atmospheres, otherwise there would be significantly more objects that display H 2 -He CIA.
3-There is a lack of continuous IR-faint sequences in the color-magnitude diagram. As discussed above, the majority of the IR-faint white dwarfs are clustered on a sequence at M g ∼ 16, and there are several IR-faint white dwarfs with M g ∼ 17 that may or may not form another sequence. Given Gaia's limiting magnitude of G ∼ 21, our survey volume is limited for these objects (see Figure 13). Now if we take the results we obtain with our revised model atmospheres at face value, we are forced to conclude that the vast majority of IR-faint white dwarfs are not as cool as previously estimated. Most of the IRfaint white dwarfs with extremely strong flux deficits in the optical also have large masses well above 0.8 M , and they occupy a region in the M vs T eff plane where very few objects -if none -have ever been identified. These must necessarily be in the Debye cooling phase.
Some of the progenitors of the IR-faint white dwarfs are likely DA white dwarfs that got convectively mixed. Bédard et al. (2022) computed for the first time the spectral evolution of a helium-rich DO white dwarf turning into a pure hydrogen DA star through the float-up process, and then into a helium-dominated DC white dwarf with trace amounts of hydrogen through convective mixing. This DO-DA-DC transition model predicts a final surface hydrogen abundance of log H/He = −4.0 around 8000 K, and naturally explains the increase in the number of helium atmosphere white dwarfs at low effective temperatures (Blouin et al. 2019b;Cunningham et al. 2020) and the observed bifurcation in the Gaia colormagnitude diagrams (Gaia Collaboration et al. 2018). Bédard et al. (2022) note that the most important parameter for the spectral evolution is the mass fraction of hydrogen in the initial models, with a larger amount of hydrogen resulting in a larger final hydrogen abundance in the atmosphere. Their evolutionary calculations do not include the IR-faint temperature regime, but strongly suggest that convectively mixed DA white dwarfs turn into IR-faint white dwarfs when they are sufficiently cool to show H 2 -He CIA. A more thorough exploration of the parameter space should shed some light on the progenitors of the IR-faint white dwarfs analyzed in this paper. Until such calculations become available, we can nevertheless conclude that the IR-faint white dwarfs with the largest hydrogen abundances in our sample cannot be explained as the result of convectively mixed DA stars alone. External sources of hydrogen must also be invoked, most likely from tidally disrupted asteroids, small planets, and comets.
We would like to thank the anonymous referee for a careful reading of our manuscript and several constructive comments that helped improve this paper. This work was supported in part by the NSF under grant AST-1906379, the NSERC Canada, by the Fund FRQ-NT (Québec), and by NASA under grant 80NSSC22K0479. SB is a Banting Postdoctoral Fellow and a CITA National Fellow, supported by NSERC.
Based on observations obtained at the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).  Table 1. Error bars show the observed photometry. Black lines show the monochromatic fluxes for the best-fit model for each star, and dots show the synthetic photometry of those models in each filter. Figure 21. Additional IR-faint white dwarf candidates identified in our analysis. Three of these objects, J0027+0554, J1004−0506, and J1951+4026, are spectroscopically confirmed DC white dwarfs.