The Impact of White Dwarf Luminosity Profiles on Oscillation Frequencies

KIC 08626021 is a pulsating DB white dwarf of considerable recent interest, and first of its class to be extensively monitored by Kepler for its pulsation properties. Fitting the observed oscillation frequencies of KIC 08626021 to a model can yield insights into its otherwise-hidden internal structure. Template-based white dwarf models choose a luminosity profile where the luminosity is proportional to the enclosed mass, $L_r \propto M_r$, independent of the effective temperature $T_{\rm eff}$. Evolutionary models of young white dwarfs with $T_{\rm eff} \gtrsim$ 25,000 K suggest neutrino emission gives rise to luminosity profiles with $L_r$ $\not\propto$ $M_r$. We explore this contrast by comparing the oscillation frequencies between two nearly identical white dwarf models: one with an enforced $L_r \propto M_r$ luminosity profile and the other with a luminosity profile determined by the star's previous evolution history. We find the low order g-mode frequencies differ by up to $\simeq$ 70 $\mu$Hz over the range of Kepler observations for KIC 08626021. This suggests that by neglecting the proper thermal structure of the star (e.g., accounting for the effect of plasmon neutrino losses), the model frequencies calculated by using an $L_r \propto M_r$ profile may have uncorrected, effectively-random errors at the level of tens of $\mu$Hz. A mean frequency difference of 30 $\mu$Hz, based on linearly extrapolating published results, suggests a template model uncertainty in the fit precision of $\simeq$ 12% in white dwarf mass, $\simeq$ 9% in the radius, and $\simeq$ 3% in the central oxygen mass fraction.


INTRODUCTION
White dwarfs (WDs) are the final evolutionary state of stars whose zero age main-sequence (ZAMS) mass is 8 M (Liebert 1980;Fontaine et al. 2001;Hansen 2004), which for a Salpeter initial mass function is 98% of stars in the Milky Way (e.g., Salpeter 1955;Scalo 1986;Maschberger 2013). The interiors of WDs encapsulate their stellar evolution history, especially the nuclear reactions and mixing that take place during the helium burning stage (Metcalfe et al. 2002;Metcalfe 2005;Fields et al. 2016;De Gerónimo et al. 2017 and the initial cooling that takes place when the WD is newly born. The pulsation properties of variable WDs are sensitive to their mechanical and thermal structure, and hence asteroseismology offers the potential to probe the interior structure and prior evolution history (Kawaler et al. 1985;Brassard et al. 1992;Fontaine & Brassard 2008;Winget & Kepler 2008;Aerts et al. 2010;Althaus et al. 2010;Romero et al. 2012Romero et al. , 2017. KIC 08626021 is a pulsating, He I line dominated WD belonging to the DBV and V777 Her classes (Winget et al. 1982) and the first to be extensively monitored by Kepler for its pulsation properties (Østensen et al. 2011). KIC 08626021 shows a frequency spectrum composed of non-radial, low order g-modes, which are sensitive to the interior stratifications of the WD. Bischoff-Kim & Metcalfe (2011) identified seven oscillation modes from 36 months of Kepler photometric data, some with triplet and doublet structures, to identify the spherical harmonic and m of several modes (also see Zong et al. 2016). Fitting these modes to ab initio WD models (e.g., WDEC, Bischoff-Kim & Montgomery 2018), they found an effective temperature T eff = 29, 650 K, mass M = 0.55 M and evidence for a thin He layer. Giammichele et al. (2016Giammichele et al. ( , 2017a pioneered new techniques to fit observed pulsation frequencies by using flexible, parameterized WD template models, and Giammichele et al. (2018) combined eight oscillation modes with a template model to infer KIC 08626021 has a large oxygen-dominated interior region. Comparison of the density, temperature, mass fraction, and luminosity profiles between the original and relaxed MESA models. The relaxed model is slightly hotter in the core, but density profiles are very similar. The mass fraction profile is identical by construction. The relaxed model has the imposed Lr ∝ Mr profile, which is significantly different than the luminosity profile of the original model.
Template models typically use a WD luminosity profile L r ∝ M r , which usually assumes the WD has largely forgotten its previous evolution history. On the other hand, stellar evolution models of young WDs with T eff 25, 000 K suggest neutrino emission dominates the energy loss budget for average-mass carbon-oxygen (CO) WDs, which yields luminosity profiles with L r ∝ M r (e.g., Vila 1966;Kutter & Savedoff 1969;Winget et al. 2004;Bischoff-Kim & Montgomery 2018). In this Letter we explore the difference this causes in the low order g-mode oscillation frequencies by comparing two nearly identical WD models: one model has an enforced L r ∝ M r luminosity profile and the other model has a luminosity profile determined by the star's evolution history.
In Section 2 we present a MESA model aiming towards KIC 08626021. In Section 3 we relax this model to have a L r ∝ M r luminosity profile while keeping other characteristics unchanged. In Section 4 we compare the GYRE oscillation frequencies and weight functions of the two models, in Section 5 we discuss cooling of the WD MESA model, and we discuss the implications of our findings in Section 6.

AN EVOLUTION MODEL AIMING AT KIC 08626021
We use release 10398 of the MESA software instrument (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018, together with the published set of controls from Farmer et al. (2016) and Fields et al. (2018), to evolve a 2.10 M , Z = 0.02 metallicity model from the zero-age main sequence (ZAMS) through 12 thermal pulses on the asymptotic giant branch (AGB) to a CO WD. The model includes rotation (solid body Ω/Ω crit = 1.9 × 10 −4 applied at ZAMS), wind mass loss (Reimers on the red giant branch with a scaling factor of 0.1 and Blöcker on the AGB with a scaling factor of 0.5), and a 49-isotope reaction network. After winds have reduced the hydrogen envelope mass to 0.01 M , we strip all of the remaining hydrogen from the surface to form a young DB WD, which we then cool until L = 0.137 L , matching the luminosity inferred for KIC 08626021 (Giammichele et al. 2018). Element diffusion is enabled during WD cooling, resulting in a pure He atmosphere and smooth interior composition transitions. The final WD model has a mass M = 0.56 M , radius R = 0.014 R , effective temperature T eff = 29, 765 K, surface gravity log g = 7.90, angular momentum J = 1/380 J , central density ρ c = 2.8 × 10 6 g cm −3 , central temperature T c = 4.8 × 10 7 K, central 16 O mass fraction X 16 = 0.74, central 22 Ne mass fraction X 22 = 0.02, and a location where the core transitions from being 16 O dominated to 12 C dominated of M trans = 0.28 M . Other than X 16 , X 22 , and M trans , this model shares many of the scalar properties with those derived for KIC 08626021 (Giammichele et al. 2018). This WD model is referred to below as "original model". Our inlist, profile, and history files are available at http://mesastar.org.

IMPOSING A LUMINOSITY PROFILE
Rather than evolve a stellar model from the pre-main sequence to a WD, or evolve a hot initially polytropic model to a WD (e.g., Bischoff-Kim et al. 2014), it can be convenient to assume the WD has forgotten its previous evolution history by assigning a specific profile. An example is prescribing the luminosity profile L r ∝ M r .
To facilitate comparing the pulsation properties of an evolutionary WD model with a L r ∝ M r WD model, we modify the original model to achieve a L r ∝ M r profile. The energy conservation equation guarantees that a model in thermal equilibrium ( grav = 0) with a constant energy generation rate nuc will satisfy L r ∝ M r (with nuc setting the constant of proportionality). Thus, we replace the usual nuclear energy generation calculations with a constant nuc throughout the model. We run MESA until this model relaxes to hydrostatic and thermal equilibrium; during this relaxation process, we allow no changes to the abundance profile (due, e.g., to diffusion or burning). We then compare T eff against the effective temperature of the original model, and iterate on nuc until the two match. This process typically leads to a surface gravity different than that of the original model. Therefore, we add a second iteration where we adjust the WD mass M while holding abundance profiles fixed as a function of fractional mass coordinate M r /M. Iterating on both M and nuc simultaneously, we obtain a model with L r ∝ M r that agrees with the T eff and log g of the original model to better than 0.001 %. This model, with nuc = 0.471 erg g −1 s −1 and M = 0.564 M , can be regarded as a 'spectroscopic twin' to the original model, because it shares the same effective temperature, surface gravity, and abundances. Our inlist for creating this model is available at http://mesastar.org. Figure 1 compares the density, temperature, abundance, and luminosity profiles of the original and relaxed models. The density profiles are very similar, while the abundance profiles are identical. The relaxed model follows the desired L r ∝ M r profile, with a larger central temperature (and temperature gradient) than the original model. Figure 2 shows the propagation diagram (e.g., Unno et al. 1989) for dipole ( = 1) modes of the original and relaxed models. The upper panel plots the square of the Lamb and Brunt-Väisälä frequencies as a function of fractional mass M r /M, for the two models. The lower panel shows the relative difference between these critical frequencies. While the models exhibit almost identical S 2 , the relaxed model has a larger N 2 in its core than the original model by up to ≈ 14%; while N 2 is smaller in the surface layers by up to ≈ 4%.

FREQUENCY DIFFERENCES
Assuming the magnitude of the temperature gradient in the interior to be much smaller than the adiabatic gradient, the Brunt-Väisälä frequency reduces to (e.g., Cox 1980; Bildsten    where Γ 1 is the first adiabatic index, χ ρ is the density exponent [∂(lnP)/∂(lnT )] ρ,µI , H is the pressure scale height, k B is the Boltzmann constant, m p is the mass of the proton, and µ I is the mean molecular weight of the ions. There is a degeneracy between the T and µ I profiles in their effect on the Brunt-Väisälä frequency, N 2 ∝ T /µ I . An inaccurate interior temperature profile can then directly impact inferences about the interior composition profile.
The Brunt-Väisälä frequency differences in Figure 2 translate into corresponding differences in the g-mode frequencies. We demonstrate this in Fig. 3, which plots the difference between the = 1 and = 2 g-mode frequencies of the relaxed model (ν relax ) and the original model (ν orig ), as a function of ν orig . Frequencies are calculated using release 5.2 of the GYRE software instrument (Townsend & Teitler 2013;Townsend et al. 2018), and selected modes are labeled by their radial orderñ in the Takata (2006) extension to the standard Eckart-Osaki-Scuflaire classification scheme described e.g. by Unno et al. (1989). The figure reveals frequency differences ranging from ≈ −20 µHz up to ≈ 70 µHz. There is no obvious pattern to the frequency differences, from one mode to the next, although the scatter appears to reduce toward larger values of |ñ|.
The frequency differences demonstrated here can be regarded as a measure of the error introduced by assuming L r ∝ M r during seismic modeling. Therefore, although Giammichele et al. (2018) report that their model frequencies match the observed frequencies to better than ≈ 0.6 nHz, the true error will likely be orders of magnitude larger -and may have an impact on the conclusions regarding core mass, radius, composition etc. that they draw from their modeling. Figure 4 compares the weight functions of the original and relaxed models for pairs of adjacent radial order = 1 and = 2 adiabatic modes The two pairs are chosen so that one mode shows a large frequency difference between the original and relaxed models, but its radial order neighbor shows a small frequency difference. Following Kawaler et al. (1985), the weight function is dζ dr = [C(y, r) + N(y, r) + G(y, r)]ρr 2 r=R r=0 T (y, r)ρr 2 dr , where C(y, r) contains the Lamb frequency, N(y, r) varies with the Brunt-Väisälä frequency, G(y, r) involves the gravitational eigenfunctions, T (y, r) is proportional to the kinetic energy density, and y = (y 1 , y 2 , y 3 , y 4 ) are the Dziembowski (1971) variables. The frequency of an adiabatic mode is then The weight function for the original and relaxed models is dominated by the N(y, r) term except for the surface layers, and the change in the weight function due to a change in N 2 is That is, the change in the weight function in going from the original to the relaxed model is given by the weight function of the original model times the fractional change in N 2 . The lower panel of Figure 2 shows δN 2 /N 2 is positive for M r /M 0.6 and negative by a smaller amount for M r /M 0.6, again ignoring the surface layers. M r /M 0.6 corresponds to the location where the C-O profiles cross. Equation 5 predicts these changes in N 2 will increase the weight function in the inner region (M r /M 0.6) of the relaxed model, and decrease the weight function in the outer region (M r /M 0.6) of the relaxed model by a smaller amount. The left column of Figure 4 verifies the expectations from Equation 5. The C-O crossover occurs at r/R 0.5. The amplitude of the weight functions at the crossover is relatively small for the = 1,ñ = −5 mode but larger for the = 2,ñ = −7 mode. Below the crossover, the amplitude of the weight function in the relaxed model is larger than in the original model. Above the crossover, the amplitude in the relaxed model is smaller. In addition, the change in the amplitude below the crossover is larger than the change in the amplitude above the crossover. The net effect of the weight function redistribution on either side of the crossover causes a shift in ζ, the area under the weight function curves, towards larger mode frequencies in the relaxed model (56.6 µHz for the = 1,ñ = −5 mode and 71.0 µHz for the = 2,ñ = −7 mode).
The right column of Figure 4 shows the same behavior; larger amplitudes below the C-O crossover, smaller amplitude decreases above the crossover. However, the weight functions are concentrated toward the surface layers where N 2 does not significantly change (see Figure 2). Hence, these mode frequencies are relatively unaffected in transitioning from the original to the relaxed model (2.3 µHz for the = 1,ñ = −6 mode and 0.3 µHz for the = 2,ñ = −8 mode). Figures 2 and 4 encapsulate two additional messages. First, the modes which best probe the interior, those whose weight functions are large in the interior, are also the modes most affected by the change in the thermal structure in transitioning from the L r ∝ M r original model to the L r ∝ M r relaxed model. Second, we do not find a large switch between a mode being confined to the core to being confined to the envelope when transitioning from the original to the relaxed model. That is, the physics changing the mode frequencies is the thermal profile rather than mode trapping.

COOLING OF THE EVOLUTION MODEL
When does L r ∝ M r hold in our original model? Figure 5 shows the evolution of the original model as the WD cools down. Plasmon neutrino emission dominates the energy loss budget for average-mass CO WDs with T eff 25, 000 K (e.g., Vila 1966;Kutter & Savedoff 1969;Winget et al. 2004;Bischoff-Kim & Montgomery 2018). The lower end of this range overlaps the observed T eff for the DBV and V777 Her classes, of which KIC 08626021 is a member. Neutrino losses explain why the T eff = 29, 765 K luminosity profile is negative (inward heat flux) from the center to M r /M 0.6.
As the evolution of the original model continues, photons leaving the WD surface begin to dominate the cooling as the electrons transition to a strongly degenerate plasma (van Horn 1971). Energy transport in the interior is dominated by conduction, driven primarily by electron-ion scattering. Energy transport in the outer layers is dominated by radiation or convection associated with the partial ionization of the most abundant element near the surface (e.g., Winget & Kepler 2008;Althaus et al. 2010). For DBV WDs, the partial ionization of He occurs around T eff 30, 000 K, leading to convection and pulsations in relatively hot WDs. Figure 5 shows the relation L r ∝ M r is approximately satisfied in the interior of our CO WD model only after T eff 20, 000 K. These lower temperatures are associated with DAV WDs, where partial ionization in their hydrogen atmospheres leads to the onset of convection and pulsations. This suggests that the approximation L r ∝ M r may be more reliable for asteroseismic studies of this cooler class of objects.

CONCLUSIONS
We have generated a pair of WD models that are 'spectroscopic twins', having the same effective temperature, surface gravity, and abundances. One model has an enforced L r ∝ M r luminosity distribution and the other model has a luminosity distribution given by the star's previous evolution. The loworder g-mode oscillation frequencies of the two models differ by up to 70 µHz, but in an uneven manner. This result suggests that by neglecting the proper thermal structure of the star (e.g., accounting for the effect of plasmon neutrino losses), the model frequencies calculated by using an L r ∝ M r model may have uncorrected random errors as large as 70 µHz. To aid interpretation of this frequency difference, Table 9 of Giammichele et al. (2017b) lists the uncertainty in the derived WD mass, radius, and central oxygen mass fraction for frequency differences of 0.001, 0.01, and 10 µHz. For example, a 10 µHz frequency difference translates into an uncertainty in the fit precision of 4% in the WD mass, 3% in the WD radius, and 1% in the central 16 O mass fraction. Their Table 9 shows a factor of 10 increase in the frequency difference causes about an order of magnitude increase in the uncertainty of the derived WD properties, with a fitting trend of larger masses, smaller radii, and smaller 16 O mass fractions for larger frequency differences. Assuming this trend holds for larger frequency differences, then a linear extrapolation to a mean frequency difference of 30 µHz betwen the original L r ∝ M r and relaxed L r ∝ M r models suggests an uncertainty in the template model fit precision of 12% in WD mass, 9% in the WD radius, and 3% in the central 16 O mass fraction. Alternatively, a frequency difference of 60 µHz for modes that are especially sensitive to the core composition translates to a 6% uncertainty in the central 16 O mass fraction. Finally, Figure 1 shows a 10-20% difference in the core temperature between the original L r ∝ M r and relaxed L r ∝ M r models. If the scaling of Equation 2 is roughly correct, there may be a 10-20% difference in the derived µ I , which may translate into larger uncertainties in the C/O mass fractions.
We encourage future generations of T eff 20, 000 K WD template models to consider using luminosity profiles informed by evolution models or to include a luminosity profile as part of the template model fitting process.