Improved seismic model of the pulsating DB white dwarf KIC 08626021 corrected from the effects of neutrino cooling

Asteroseismology is a powerful tool to unravel the chemical composition and stratiﬁcation inside white dwarfs, as recently achieved by Giammichele et al. (2018, Nature, 554, 73) for the pulsating DB star KIC 08626021. However, Timmes et al. (2018, ApJ, 867, L30) pointed out that neglecting the e ﬀ ects of neutrino cooling, such as in the models used in Giammichele et al. study, could signiﬁcantly impact the derived seismic solution and compromise conclusions drawn upon it. In this context, we perform a complete reevaluation of the seismic solution uncovered for KIC 08626021, using improved static models which incorporate more realistic luminosity proﬁles that reﬂect the still signiﬁcant energy losses induced by neutrino emission mechanisms in hot DB white dwarfs. We ﬁnd that including (or neglecting) neutrino cooling for the speciﬁc case of KIC 08626021 induces frequency di ﬀ erences of ∼ 35 µ Hz on average (with variations up to ∼ 84 µ Hz) for the relevant g -modes, that is, similar to the frequency shifts estimated in Timmes et al. study. However, we show that the propagation of these variations into the derived seismic model properties remain limited and mainly trigger changes of the C / O and C / He composition ratio in the intermediate layers of the seismic model, while other important parameters are only slightly a ﬀ ected. In particular, the derived central oxygen mass fraction and extent of the homogeneous inner part of the core are essentially unchanged. Hence, as found by Timmes et al., seismic investigations of hot pulsating DB white dwarfs that rely on parameterized static models should include the non-negligible e ﬀ ects of neutrino cooling to provide more accurate solutions, but all the important conclusions brought by Giammichele et al. from the analysis of KIC 08626021 remain entirely valid.


Introduction
For more than three decades, white-dwarf asteroseismology has developed upon the promise of revealing the inner structure of these objects, which represent the ultimate stage of evolution for stars that were born less massive than ∼8 M (see Sect. 2.1 of Giammichele et al. 2017a andWinget & Kepler 2008;Althaus et al. 2010 for reviews on this venture). Recently, a breakthrough occurred with the development of a new technique to explore the C/O core chemical stratification inside white dwarfs, mostly independent of stellar evolution calculations (Giammichele et al. , 2017a. This approach found its first application on the pulsating DB white dwarf KIC 08626021 (Giammichele et al. 2018), a star extensively monitored by the Kepler satellite (Østensen et al. 2011;Bischoff-Kim et al. 2014;Zong et al. 2016).
The study of Giammichele et al. (2018) most importantly established that KIC 08626021 has a central homogeneous core that is ∼15% richer in oxygen (in mass) and ∼40% more massive than expected from standard evolutionary models, thus providing strong and challenging constraints to the modeling of physical processes shaping stellar cores during the preceding helium-burning phase (e.g., Fields et al. 2016;De Gerónimo et al. 2017. However, this result is being questioned after Timmes et al. (2018) pointed out that the impact of neutrino cooling on the thermal structure of the star -an effect neglected in the static equilibrium models employed by Giammichele et al. (2018) -is still significant in a hot DB white dwarf like KIC 08626021. These authors found that low-order gmode frequencies could differ by as much as ∼70 µHz (∼30 µHz on average) between hot DB white dwarf models neglecting and including neutrino losses. Such differences exceed the limit set by Giammichele et al. (2018) for potentially unaccounted systematic effects in their estimation of the error budget associated to the seismic solution. This suggests that significant shifts in the derived parameters of KIC 08626021 could occur when accounting for neutrinos and motivates a reanalysis of this star.
In the present Letter, we address this important issue by introducing new improved static models for white-dwarf asteroseismology that incorporate luminosity profiles accounting for residual contraction and energy losses from neutrino emission processes (Sect. 2). In Sect. 3, the impact of this correction on the g-mode frequency spectrum is investigated, specifically for the optimal model uncovered by Giammichele et al. (2018), and is compared to results of Timmes et al. (2018). A reanalysis of KIC 08626021 based on the improved models then follows in Sect. 4, and the revised seismic solution is discussed in Sect. 5.

Seismic models with neutrino cooling
The models used in Giammichele et al. (2018) are static structures whose advantage is, among other potential applications, to provide a flexible way to explore white dwarf internal configurations with asteroseismology. Details about how they are constructed are provided in Giammichele et al. (2016Giammichele et al. ( , 2017a. One assumption of these models is that the luminosity profile is proportional to the mass distribution (L r ∝ M r ; hereafter referred to as the "zeroth-order" approximation 1 ), which has been robustly verified in white dwarfs that have cooled sufficiently, such as DA pulsators found in the T eff = 13 000-11 000 K range. However, for white dwarfs that are still relatively hot, in particular for the helium-atmosphere DB pulsators found in the T eff = 32 000−22 000 K range, cooling from neutrino emission is still important, with a significant impact on the thermal structure and pulsation frequencies, as illustrated by Timmes et al. (2018), and the above assumption on the luminosity profile is no longer a good approximation.
To address this problem, we introduce improved static models that no longer rely on the L r ∝ M r assumption. Instead, the models are constructed from the relation L r = M r · (L * /M * − k ν ), where k ν accounts for all effects (mainly neutrino cooling, but also residual contraction) that would cause the luminosity profile to differ from the L r ∝ M r relation. This quantity, k ν , is interpolated from a 4D grid of luminosity profiles built after the calculation of 2431 white-dwarf cooling sequences (evolutionary models) covering a wide range of values for T eff , mass, Heenvelope thickness, and core composition (C/O relative mass fraction). We note that models forming this grid are "generic" and consist of a core with a homogeneous composition and a pure helium envelope, providing the main, "first-order" correction to the luminosity profile. A "second-order" correction accounting for a particular chemical stratification, for example, as inferred from a seismic model solution, can be obtained if necessary by calibrating the interpolation of k ν relative to a specific cooling sequence (see Appendix A). This approach produces improved static white-dwarf structures that now account for the impact of neutrino cooling (and other less important effects). Figure 1 compares the model identified by Giammichele et al. (2018) as the best seismic representation of KIC 08626021 (assuming L r ∝ M r ) to an improved version of itself (of the same parameters) incorporating the correction (up to second order) for neutrino cooling. As expected, the addition of neutrinos induces a significant change in the luminosity profile and results in a much cooler temperature in the core, with an off-centered maximum. As a check, the luminosity and temperature profiles obtained from an evolutionary model tuned to represent KIC 08626021, that is with the same mass, effective temperature, and chemical stratification, are also plotted in Fig. 1. The comparison demonstrates that our improved static structures correctly implement the effect. The differences in the thermal structure due to neutrino losses also have a slight but non-negligible impact on the Brunt-Vaisälä frequency profile. This is illustrated in the bottom panel of Fig. 1 and the differences found are similar to those shown in Fig. 2 of Timmes et al. (2018).

Impact on the g-mode frequencies
Table 1 compares the computed adiabatic frequencies for the most relevant g-modes corresponding to the two models shown Comparison between the optimal seismic model obtained by Giammichele et al. (2018) for KIC 08626021 which does not incorporate neutrino cooling (black dotted lines), and an equivalent model using the same parameters and chemical structure that includes this effect (red solid curves). Top panel: composition stratification and surface parameters characterizing the two models. Lower panels: from top to bottom, profiles as functions of fractional mass depth (log q) of the temperature (T ), luminosity (L r /L * ), logarithm of the squared Brunt-Väisälä, and = 1 Lamb frequencies (log N 2 and log L 2 1 ), and relative difference of the squared Brunt-Väisälä frequencies of the two models (E(N 2 ), as defined by Timmes et al. 2018). An evolution model matching the properties of KIC 08626021 is also shown (blue dashed lines) for comparison with profiles obtained from the new white-dwarf static models incorporating neutrino cooling.
in Fig. 1 (see also Table A.1). These frequencies were calculated with a modern version of the stellar oscillation code described in Brassard et al. (1992; see also Brassard & Charpinet 2008), and used by Giammichele et al. (2018). Comparisons of the computed modes indicate that systematic shifts toward lower frequencies by ∼35 µHz on average, with variations up to ∼84 µHz for the worst case, occur when neutrino cooling is considered. Hence, we find that this effect is comparable in magnitude to the variations illustrated in Fig. 3 of Timmes et al. (2018). Table 1. Frequencies of the low-radial-order (k), low-degree ( = 1, 2) g-modes computed for the optimal seismic model of KIC 08626021 derived by Giammichele et al. (2018), which does not include neutrino cooling (ν 0 ), and for the equivalent model including this effect (ν 2 ). k ν 0 (µHz) ν 2 (µHz) ν 0 − ν 2 (µHz) Notes. The modes marked by an asterisk are those fitted to the observed frequencies that closely match ν 0 . The last column lists the frequency differences for each mode, between the two models.
As pointed out by these latter authors, such differences may have an impact on the conclusions drawn from the modeling of Giammichele et al. (2018) regarding, for example, the measured core mass, radius, and composition, considering that the frequency drift exceeds the limit for potentially unaccounted systematic effects set for the determination of uncertainties associated to the proposed seismic solution. Timmes et al. (2018) extrapolated the accuracy tests provided in Giammichele et al. (2017a) to roughly estimate an error of ∼12% for the mass, ∼9% for the radius, and ∼3% for the central oxygen mass fraction triggered by the omission of neutrino cooling in the analysis of KIC 08626021. However, a more robust way to evaluate the real impact on the derived properties of KIC 08626021 is to undertake a new seismic analysis of this star using the improved models introduced above.

Improved seismic solution for KIC 08626021
We performed a complete reanalysis of KIC 08626021 with our updated static models that incorporate more realistic luminosity profiles, following exactly the same approach as described in Giammichele et al. (2018; see the main paper and the supplementary material). We refer to this work for a full description of Table 2. Updated properties of KIC 08626021 when accounting for neutrino cooling in the seismic models (this work) compared to the former solution from Giammichele et al. (2018).
Quantity Giammichele et al. (2018) This work Structural comparison between the optimal seismic model of KIC 08626021 derived by Giammichele et al. (2018; in black) and the new optimal seismic solution obtained with improved models that include neutrino cooling (up to second-order approximation; in red).
Top panel: composition stratification for helium, carbon, and oxygen. Lower panel: logarithm of the squared Brunt-Väisälä and = 1, 2 Lamb frequencies (log N 2 and log L 2 ).
the analysis, which we do not repeat here, and to Appendix B for more details about the uncovered solution. Using the new improved models has no negative impact on the achieved quality of the optimal fit, leading to the identification of a very welldefined seismic solution that also reproduces the frequencies of KIC 08626021 to the precision of the observations (∼0.6 nHz) and does not induce changes in the mode identification (see Table B.2). The obtained solution is however shifted relative to the former model uncovered by Giammichele et al. (2018). The most relevant aspects of this updated solution are synthesized in Table 2 and Fig. 2. The first striking result is that the new improved solution for KIC 08626021 remains close to the former model, as far as global stellar parameters (T eff , log g, M * , and R * ) are concerned (see Tables 2 and B.1). Shifts in these parameters remain close to the 1σ-uncertainties quoted by Giammichele et al. (2018). This suggests that the systematic errors introduced by neglecting neutrino cooling do not propagate evenly among the model parameters that are optimized to determine the best seismic solution. Larger variations are indeed revealed when comparing the internal stratifications of the new optimal model with the former one (see Fig. 2 and Fig. B.1), but again unevenly impacting the characteristics of the profiles.
The chemical composition in the intermediate C/O and C/He layers are the two quantities that are the most affected, "absorbing" the main part of the systematic effect induced by neutrinos. We find that the helium mass fraction in the mixed C/He layer, X(He) env , is now 0.38 ± 0.07 instead of 0.18 ± 0.04, while the amount of oxygen (in mass fraction) in the intermediate C/O layer reaches ∼0.6 instead of ∼0.4. Remarkably, the other important characteristics of the chemical profiles remain nearly unchanged, either within or very close to the 1σ-errors quoted by Giammichele et al. (2018). In particular, we note that all chemical transitions are found to be located nearly at the same fractional mass depths, except for a slight outward displacement of core-to-envelope transition (log q 2 ). The inner core for its part is barely affected, with the central oxygen mass fraction, X(O) centre , estimated at 0.84 ± 0.04 (i.e., very close to the former value: 0.86 ± 0.04) and the homogeneous central part extending up to a fractional mass depth of log q ∼ 0.71 (i.e., encompassing a mass of ∼0.45 M ), which is identical to the value of log q ∼ 0.7 (∼0.45 M ) found from the former solution. Consequently, all the conclusions regarding the inner core size and its C/O composition ratio expressed in Giammichele et al. (2018) are unchanged. Timmes et al. (2018) pointed out a shortcoming in using static models under the assumption that the luminosity profile is proportional to the mass distribution (L r ∝ M r ) for precise asteroseismic analysis of pulsating white dwarfs hotter than ∼20 000 K, thus including a DB pulsator like KIC 08626021 analyzed by Giammichele et al. (2018). Energy losses from neutrino production mechanisms are important enough at these temperatures to markedly affect the thermal structure of the star and its luminosity profile, with a non-negligible impact on g-mode frequencies. Taking this criticism into account, we developed new static models for these hot white dwarfs that incorporate more realistic luminosity profiles interpolated from a grid of evolutionary models. These new structures account for the effects of neutrino cooling to a good level of approximation.

Summary and conclusion
KIC 08626021 was reanalyzed using a similar approach as described in Giammichele et al. (2018), but now relying on the new models, which leads to updated values for the global properties of this star (Table 2). These new values come from a significant improvement of the seismic modeling of KIC 08626021 over the former analysis and become our new reference (we also provide the updated chemical stratification with its error estimates in Fig. 3). Despite some changes in the determined internal structure, we find that for most parameters, accounting correctly for neutrino cooling had no strong impact on the derived values. The only significantly affected quantities, which absorb most of the frequency shifts, are the C/O and C/He mass fractions in the intermediate layers of the star (Fig. 2). Other quantities, such as the location of helium transitions in the envelope or the central oxygen mass fraction and extent of the homogeneous central part of the core, are either unchanged or only slightly shifted. Therefore, despite significant effects on calculated frequencies induced by omitting neutrinos and their impact As a final note, we point out that correcting for the effect of neutrinos did not change the finding that KIC 08626021 has a very thin helium-pure layer on top of its envelope. Such a thin layer might be explained by the presence of any mechanism competing against the separation of He, C, and O in the envelopes of DO and hot DB white dwarfs. For instance, Fontaine & Brassard (2005; see also Brassard et al. 2007) postulated the existence of a weak stellar wind of the order of ∼2 × 10 −13 M yr −1 in order to account for the presence of traces of C in T eff ∼ 20 000-30 000 K DB white dwarfs. Such a wind, similar to the solar wind, would be driven by deposition of acoustic energy generated by strong convective motions which characterize the outer layers of these stars. Perhaps a more likely mechanism is convective overshooting, which was first studied by Freytag et al. (1996) in a white dwarf context and more recently by Tremblay et al. (2015) and Cukanovaite et al. (2018) through detailed 3D simulations. Brassard & Fontaine (2015) explicitly used the overshooting prescription of Freytag et al. to demonstrate that element separation is indeed considerably slowed down in cooling hot white dwarfs. They found that only a "thin" layer of pure He (log q 1 ∼ −8) has accumulated at the surface of their evolving model by the time the latter has cooled down to T eff ∼ 31 000 K. This is in line with our seismic inference for KIC 08626021. In addition, the new solution for KIC 08626021 does not change the fact that no C/O/He triple transition at the core boundary is found in the model that best reproduces the pulsation frequencies of this star. If puzzling from an evolutionary perspective, our seismic probing seemingly requires that mode trapping generated by two chemical transitions around log q ∼ −2.4 and log q ∼ −3.5 must be present. This can be obtained with our current models only if a carbon buffer exists between the envelope and the core. Alternately, enforcing a triple transition and searching for a new optimal solution within this additional constraint strongly degrades the quality of the achieved best frequency match. Improving the capacity of such a configuration to reproduce the observed pulsation modes would likely require the presence of an additional structure (of unclear origin) in the envelope to provide the needed mode trapping. We therefore leave that issue open at this stage. Comparison of three identical models in terms of input parameters, but featuring different levels of correction to the luminosity profile. The zeroth-order model assumes L r ∝ M r and corresponds to the seismic structure uncovered by Giammichele et al. (2018, dotted lines in black). The model with the first-order correction (dashed curves in blue) shows profiles determined from the interpolation of k ν in the grid of "generic" evolutionary models. The model with the second-order correction (plain curves in red) shows profiles corrected for the particular chemical stratification of the model considered.
As briefly explained in Sect. 2, the correction applied to the luminosity profile in order to include the dominant effect of neutrino cooling and, to a lesser extent, residual contraction, is obtained by interpolating the quantity k ν through a grid of pre-calculated profiles taken from "generic" cooling (evolutionary) white dwarf models. These models are generic in the sense that they have simple structures made of a homogeneous C/O core and a pure He-envelope, providing a first-order approximation of the luminosity profile. However, since this profile slightly depends on the detailed chemical stratification inside the star, a more accurate so-called second-order approximation approach may be needed in some cases. The latter is obtained by calibrating the luminosity interpolation scheme with an additional parameter, ∆T ν , that shifts the T eff values of the evolutionary models in the grid by the specified amount. The value chosen for ∆T ν is determined by finding the best match to the luminosity profile obtained from evolutionary models calculated with the specific chemical stratification. Figure A.1 compares identical models at the zeroth-order approximation (no correction for neutrinos; L r ∝ M r ), firstorder approximation accounting for the main effects, and second-order approximation correcting for the specific chemical stratification of the model (illustrated in the upper panel of Fig. 1). The calibrated value is found to be ∆T ν = −400 K in that case, and the resulting static model is very close to the corresponding evolutionary structure (as shown in Fig. 1).
The impact of the various levels of approximation on g-mode frequencies is summarized in Table A.1. As expected from the examination of Fig. A.1, the transition from zerothorder to first-order generates the largest frequency shifts (ν 0 −ν 1 ) Table A.1. Frequencies of low-radial order (k) and low-degree ( = 1 and 2) g-modes computed for the three identical models using various orders of correction for the luminosity profile. Notes. The zeroth-order frequencies (ν 0 ) correspond to the model assuming L r ∝ M r (model from Giammichele et al. 2018). The firstorder frequencies (ν 1 ) are determined with the model constructed from the interpolation of k ν in the grid of "canonical" evolutionary models. The second-order frequencies (ν 2 ) are from the model with the luminosity profile corrected for the particular chemical stratification of the model considered. Absolute differences in mode frequencies between the zeroth-and first-order, and between the first-and second-order are also given.
with variations up to +90 µHz (∼ +39 µHz on average). When correcting for the exact chemical stratification, using the secondorder approximation, changes are an order of magnitude smaller relative to the first-order approach (ν 1 − ν 2 ∼ −2.5 µHz on average). Overall this more accurate second-order model offers slightly reduced differences when compared to the zeroth-order model neglecting neutrino cooling. Even if small, in some circumstances such as in high-precision asteroseismic studies, in particular of hot DB white dwarf pulsators, it is appropriate to push the refinement up to this second-order approximation level.  Table 2 from Giammichele et al. (2018), but providing the optimal parameters obtained from various types of models (see text).

Quantity
Optimal value Error estimate ( §) 0th-order ( * ) 1st-order ( ) 2nd-order ( ) 1st-order ( †) Giammichele et al. (2018). ( ) Solution from models constructed with Akima splines (see text). ( †) Solution from models constructed with Fritsch-Carlson splines (see text). ( ‡) Some models use ∆ 2 as an alternative way to specify the transition width at log(q 2 ). publication. A slight modification of the procedure is however needed to carry out a global search that relies on a model incorporating the second-order correction to the luminosity profile (as described in Appendix A). Indeed, the chemical stratification of the star under consideration is not known a priori, as it is determined from the finding of a best-fit seismic model. Consequently, the value of ∆T ν that needs to be applied to correct the luminosity profile is not known initially. This problem is solved iteratively, by first carrying out the analysis using models at the first-order approximation level (∆T ν = 0). The first-order optimal seismic solution obtained from this initial step provides a first approximation of the chemical stratification inside the star, which is then used to calibrate ∆T ν from the computation of an evolutionary sequence using this specific stratification. The optimization is carried out a second time using models now including the second-order approximation based on the ∆T νvalue determined at the previous step, thus producing a secondorder and more accurate optimal seismic solution. This process could be continued further with additional iterations, by successive re-calibrations of the ∆T ν correction based on the chemical stratification of the seismic solution updated at each step, but in practice stopping at the second-order solution is enough, because the ∆T ν correction depends only slightly on the details characterizing the chemical profiles.
In order to illustrate this process, Table B.1 gives the optimal parameters obtained from the full seismic reanalysis of KIC 08626021 using various flavors of the static models. These solutions are also shown in Fig. B.1 which provides a visual comparison of the internal chemical stratification determined from each type of model. The zeroth-order solution corresponds to the best model uncovered by Giammichele et al. (2018) when neglecting neutrino cooling. An initial set of first-and secondorder solutions is obtained with static structures that use the same Akima-spline interpolation scheme to generate the chemical profiles in the core (see Giammichele et al. 2017a). In that respect, these structures are strictly identical to the models used in Giammichele et al. (2018), except for respectively incorporating the first-order and second-order (with a calibrated value ∆T ν = −160 K) correction to the luminosity profile accounting for neutrinos. Since these three model solutions differ only at the level of the correction applied to the luminosity profile, they are representative of the impact triggered solely by adding (or neglecting) neutrino cooling in the determination of the seismic solution for KIC 08626021.
We took the opportunity of this revision of the seismic solution for KIC 08626021 to include our current most advanced improvements in building seismic models for white dwarf stars. We also present an additional set of solutions, calculated up to the second-order correction for neutrino cooling (with a calibrated ∆T ν = −220 K), obtained from a recent update of our static structures that incorporate a new interpolation scheme based on Fritsch-Carlson (instead of Akima) splines (Fritsch & Carlson 1980). These produce smoother chemical profiles in the core, and close comparisons with full evolutionary models show that these structures can more closely match profiles derived from physical processes at work during evolution, which we consider as an improvement. This technical update induces further variations (not due to neutrinos) in the seismic solution uncovered. However, we note from Table B.1 and Fig. B.1 that the differences with former models remain small. Interestingly, all the optimal models incorporating a correction for neutrinos have similar structures, in particular regarding the core chemical stratification. This indicates that the seismic solution is quite robust to approximations associated with technical choices related to how the underlying models are constructed. Furthermore, to some extent, the differences between each type of model reflect how uncertainties of the order of ∼1−3 µHz (i.e., the typical scale of the difference in frequencies between firstand second-order models, and between models based on Akima vs. Fritsch-Carlson spline representations) propagate into errors on the parameters defining the structure of the star. Among these calculations, we consider that our most advanced model, and  Giammichele et al. (2018) that neglects the effects of neutrino cooling. The two models in magenta and blue are respectively the first-order and second-order solutions using the same type of models based on Akima splines (see text and Giammichele et al. 2017a). The model in red, which is our new adopted representation of KIC 08626021, is the second-order solution obtained with models that use the smoother Fritsch-Carlson spline interpolation to generate the profiles defining the core chemical stratification (the first-order solution found with these models is represented in green). All quantities shown in this figure are described in Fig. 2 caption. new reference to represent KIC 08626021, is the solution incorporating the neutrino correction up to second-order obtained from the models constructed with Fritsch-Carlson splines (i.e., the column appearing in boldface-font in Table B.1 and the red profiles in Fig. B.1, also reproduced in Fig. 2).
In the following, we provide some details about the asteroseismic reanalysis of KIC 08626021, focusing only on the new reference solution mentioned above. Figure B.2 is the equivalent of the Extended Data Fig. 3 of Giammichele et al. (2018) that shows the location of the seismic solution in the log g − T eff plane in comparison with the spectroscopic error box. Our updated solution is found to be at a slightly higher effective temperature than before, but still lying within the estimated 1σ-error for the spectroscopic T eff , while the derived log g value is entirely consistent with the spectroscopic measurement. Figure B.4 is the updated version of the Extended Data Fig. 4 in Giammichele et al. (2018), which provides the marginalized distributions for the most interesting quantities. The values given in Table B.1 and Table 2 are from the optimal model (the blue vertical dashed lines in Fig. B.4), while the errors are evaluated from the 1σ-range defined by fitting a Gaussian to these distributions. Finally, Fig. B.3 is the equivalent of Extended Data Fig. 5 of Giammichele et al. (2018), which shows the derived probabilitydistribution functions for the oxygen, carbon, and helium profiles.