Limits on Mode Coherence in Pulsating DA White Dwarfs Due to a Nonstatic Convection Zone

In the linear, adiabatic theory of stellar pulsations, modes are considered to be perfectly sinusoidal in time. This results in a theoretical eigenmode spectrum with arbitrarily thin, delta function peaks as a function of frequency. In the linear, nonadiabatic theory, modes are allowed to gain or lose energy with their environment, leading to nonzero growth/damping rates, and these rates may be related to the finite widths of peaks in their power spectra. In particular, the widths of modes in solar-like pulsators can be linked to their linear damping rates (e.g., Kumar & Goldreich 1989; Houdek & Dupret 2015).

In the white dwarf (WD) regime, we have direct observational evidence that some modes in pulsating white dwarfs can show a high degree of phase coherence, in a few cases spanning decades. For instance, more than 40 yr of observations have established that the phase of the 215 s mode in the pulsating hydrogen-atmosphere white dwarf (DAV) G117-B15A is coherent over this timescale. Thus, we can show that in this DAV the pulsation period, P, is changing (taking into account the proper motion) at the extremely slow rate of $\dot{P}=(3.57\pm 0.82)\times {10}^{-15}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ (Kepler et al. 2005). The extreme sensitivity of such $\dot{P}$ measurements in this and other white dwarfs has allowed them to be used as testbeds for physical processes that could affect their cooling, such as the interaction of hypothetical dark matter particles (e.g., Bischoff-Kim et al. 2008; Isern et al. 2008; Córsico et al. 2012, 2016).

Another use of the observed stability of these modes is searching for planetary signals in the delayed and advanced light arrival times due to reflex orbital motion of the white dwarf (Mullally et al. 2003, 2008; Hermes et al. 2010; Winget et al. 2015), although no planets have been positively identified orbiting WDs with this technique to date. We note that while the coherent modes that are useful for these studies are mostly found in DAVs near the hot edge of the instability strip, striking seasonal changes in the Fourier spectra of cooler DAVs are commonly observed (e.g., Kleinman et al. 1998). The transition from pulsators with stable Fourier spectra near the blue edge to those showing frequent amplitude changes near the red edge occurs somewhere in the middle of the observed DAV instability strip, in the vicinity of ${T}_{\mathrm{eff}}$ ∼ 11,500 K.

While a handful of coherent modes have been studied in DAVs, no systematic study of mode coherence in a large sample of these stars has been made from the ground. This situation has improved greatly with the launch of the Kepler spacecraft. During its original mission and the follow-on K2 mission, it has obtained nearly continuous time series data, often exceeding 75 days, of a large number of pulsating white dwarf stars. Recently, Hermes et al. (2017) published comprehensive data on the first 27 DAVs studied by Kepler. One of their central results was that longer-period modes (P ≳ 800 s) were observed to have larger Fourier widths than shorter-period modes (P ≲ 800 s), essentially dividing the modes into two populations; this result can hold even for different modes in the same star. We present an explanation for this phenomenon based on whether the mode propagates to the base of the nonstationary convection zone. We also show that this leads to damping of the modes and that the effect is larger near the red edge of the DAV instability strip.

The extended length of observations (≳75 days for most stars) in the Kepler and K2 data sets results in a 1/T resolution in the power spectra of $\lt 0.14\,\mu \mathrm{Hz};$ this sets the observable lower limit for the width of a peak in the Fourier transform. For the first time this enables the measurement of the widths of a large number of modes in many stars that are above this threshold (Bell et al. 2015). For their sample of 27 DAVs, Hermes et al. (2017) obtained follow-up spectroscopy with the WHT and SOAR telescopes; these spectra were fit using the techniques and models of Tremblay et al. (2011) to obtain values of ${T}_{\mathrm{eff}}$ and $\mathrm{log}\,g$ for each star.

Hermes et al. (2017) find that the Fourier width of modes (HWHM, the half width at half maximum of a Lorentzian fit), is a strong function of the mode period, P. To summarize, they find that (1) modes with $\mathrm{HWHM}\gt 0.3$ μHz have P ≳ 800 s and (2) modes with P ≲ 800 s have $\mathrm{HWHM}\lt 0.3$ μHz. This is illustrated in Figure 1, in which we have plotted mode width versus period for the linearly independent periods found in the sample of DAVs from Hermes et al. (2017).⁵ Modes from stars with ${T}_{\mathrm{eff}}$ > 11,500 K are shown as blue squares, while those from stars with ${T}_{\mathrm{eff}}$ < 11,500 K are shown as red circles. The fact that most of the modes with $\mathrm{HWHM}\gt 0.3$ μHz are from the cooler population is not an independent piece of information since longer-period modes are typically found only in the cooler DAVs (Mukadam et al. 2006).

Zoom In Zoom Out Reset image size

The region of propagation of g modes ("gravity" or "buoyancy" modes) in a star is defined by the region in which ${\omega }^{2}\lt {N}^{2},{L}_{{\ell }}^{2}$, where ω = 2π/P is the angular frequency of the mode, N is the Brunt–Väisälä (buoyancy) frequency, L_ℓ is the Lamb (acoustic) frequency, and ℓ is the spherical degree of the mode. In Figure 2, we show a propagation diagram for a model, computed with MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019). The black horizontal line denotes the region of propagation of a hypothetical mode whose period, ${P}_{\mathrm{cross}}$, is the minimum required for it to propagate to the base of the convection zone.

Zoom In Zoom Out Reset image size

In Figure 3, we show the outer propagation region for a model with the same parameters as EPIC 201806008 ($\mathrm{log}\,g=8.02$, ${T}_{\mathrm{eff}}$ = 10,910 K). The dashed horizontal lines show the region of propagation that the observed modes would have in this model, where the blue dashed lines denote modes observed to have $\mathrm{HWHM}\lt 0.3$ μHz and the red dashed lines modes with $\mathrm{HWHM}\gt 0.3$ μHz. We note that, by using the convection model $\mathrm{ML}2/\alpha =0.42$, where α is the mixing length to pressure scale height ratio (see Böhm & Cassinelli 1971), we can divide the modes into two groups: (1) narrow (blue) modes whose outer turning point is well beneath the convection zone, and (2) wide (red) modes that propagate all the way to the base of the convection zone. In other words, all the modes with $\mathrm{HWHM}\gt 0.3$ μHz can be explained as propagating to the base of the convection zone, whereas all the modes with $\mathrm{HWHM}\lt 0.3$ μHz have an outer turning point safely inside this point. Our hypothesis is that the long-period modes, through their interaction with the convection zone, will have systematically larger Fourier mode widths than the short-period modes. In the following sections we examine this statement more quantitatively.

Zoom In Zoom Out Reset image size

Since the 27 DAVs in Hermes et al. (2017) have different atmospheric parameters, each star will have a different value of ${P}_{\mathrm{cross}}$ (note: ${P}_{\mathrm{cross}}$ is also a function of the ℓ value of each mode). Thus, for each mode in each star we calculate ${P}_{\mathrm{cross}}$ and from this we predict whether each observed mode in the star has a "wide, mottled" peak in the FT or a "narrow" peak. We assume the ML2 convection model of Böhm & Cassinelli (1971) with $\mathrm{ML}2/\alpha =0.76$. This value of α is taken from the 3D simulations of Tremblay et al. (2015) for a model with ${T}_{\mathrm{eff}}$ = 11,500 K and $\mathrm{log}\,g=8.0$.

We plot the results of this in Figure 4, where we have used different symbols for the different ℓ identifications: blue triangles for ℓ = 1 modes, green diamonds for ℓ = 2 modes, and red stars for unidentified modes. Unidentified modes are assumed to have ℓ = 1 for this analysis. Modes with wide peaks (i.e., $\mathrm{HWHM}\gt 0.3$ μHz) are predicted to interact strongly with the convection zone. If this prediction is correct (i.e., if $P\geqslant {P}_{\mathrm{cross}}$ for the mode) then they are plotted with a filled symbol. If the prediction is incorrect, they are plotted with an unfilled symbol. Similarly, modes with narrow peaks (i.e., $\mathrm{HWHM}\lt 0.3$ μHz) are predicted to interact weakly with the convection zone. If this prediction is correct (i.e., if $P\lt {P}_{\mathrm{cross}}$) then they are also plotted with a filled symbol; otherwise, an unfilled symbol is used. Thus, the filled symbols represent the modes whose widths are correctly predicted by our hypothesis. We find that this procedure correctly classifies the observed mode widths with an accuracy of ∼81%. While there is evidence that the value of α is a function of ${T}_{\mathrm{eff}}$ (Provencal et al. 2012; Tremblay et al. 2015), the percentage of correctly predicted modes is roughly constant for α in the range 0.6–0.9.

Zoom In Zoom Out Reset image size

Modes of different spherical degree (ℓ) are also expected to interact with the convection zone differently as a function of period. Since the Lamb frequency of ℓ = 2 modes is larger than it is for ℓ = 1, for a fixed frequency ℓ = 2 modes have an outer turning point that is closer to the stellar surface, and are therefore more likely to strongly interact with the convection zone. If this can lead, in some cases, to not just a broadened peak in the FT but to a partial or complete suppression of the mode, then we would expect the number of ℓ = 2 modes to be suppressed by this mechanism.

Hermes et al. (2017) classify all of the statistically significant peaks in the FTs as either ℓ = 1, ℓ = 2, or undetermined. They were able to identify the ℓ value of 87 out of 201 independent modes in these DAVs (more than 40%), based on common frequency patterns in the Fourier transforms (e.g., the rotational splitting of multiplets and the period spacing of different radial overtones).

In Figure 5, we use "kernel density estimation" (kde) to estimate the density of modes as a function of the heights of the Lorentzians that were fit to them in the power spectrum, which we call their amplitude (Pedregosa et al. 2011). Only modes containing peaks with less than a 0.1% false alarm probability (FAP) and whose frequencies are not linear combinations of other mode frequencies are considered. As expected, the ℓ = 1 modes (blue curve) outnumber ℓ = 2 modes (pink curve). We next wish to test the hypothesis that there are as many ℓ = 2 as ℓ = 1 modes, but that the difference in observed numbers is solely due to geometric cancellation reducing their observed amplitudes below detection limits.

Zoom In Zoom Out Reset image size

To test this, we divide the observed ℓ = 1 amplitudes by a factor of 2.4 to simulate the larger geometric cancellation of ℓ = 2 modes observed in the Kepler passband. The result is shown as the green dashed curve in Figure 5. We see that the number of ℓ = 2 modes obtained in this way exceeds those observed by a factor of more than 2. This could indicate that a mechanism in addition to pure geometric cancellation is limiting their number. Of course, without a way of identifying the unidentified modes that comprise ∼60% of this sample, it is not possible to conclusively demonstrate this.

Finally, we note that after scaling the observed ℓ = 1 amplitudes by the appropriate factors for ℓ = 3 and 4 (∼26 and 19, respectively), no modes are found to be above the FAP > 0.1% detection threshold. Thus, these data have nothing to say on the possible presence or amplitude distribution of these higher-ℓ modes.

There are two main ways that a mode's interaction with the outer convection zone could lead to a loss of its coherence. First, if a mode has a low enough frequency, its outer turning point will be the base of the convection zone. Since the base of the convection zone rises and falls with the surface temperature perturbations of the pulsations, the mode will sometimes have to travel a greater (or lesser) distance before being reflected, and in so doing it will acquire an extra phase Δϕ (see Figure 6). A second effect that may play an even larger role is the Doppler shift of the wave caused by reflection from the base of the moving convection zone (see Figure 7). Essentially, the reflection causes a frequency shift of the reflected wave, which again leads to a phase mismatch when the wave next returns to the outer turning point.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Montgomery et al. (2015), we showed that the extra phase acquired by a mode with quantum numbers {n, ℓ} due to a change in the size of its propagation cavity could be approximately calculated from the difference in period of this mode in two models (ΔP) that differ only in the depth of their convection zones; this extra phase is given by

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{cav}}=2\pi n\left(\displaystyle \frac{{\rm{\Delta }}P}{P}\right).\end{eqnarray} \tag{ 1 }$$

We also calculate this phase difference by directly integrating the asymptotic expression for the radial wavenumber. From Gough (1993) we have that the asymptotic radial wavenumber K(r) is given by

$$\begin{eqnarray}&&{K}^{2}(r)=\displaystyle \frac{{\omega }^{2}-{\omega }_{c}^{2}}{{c}^{2}}-\displaystyle \frac{{L}^{2}}{{r}^{2}}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where N is the Brunt–Väisälä frequency, ${L}^{2}\equiv {\ell }({\ell }+1)$, c is the sound speed, ω is the angular frequency of the mode, and ω_c is the acoustic cutoff frequency. The phase difference is then

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{\phi }_{\mathrm{cav}} & = & 2\,{\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})\\ & \equiv & 2\left[{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} -{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}0}}{K}_{0}(r^{\prime} ){dr}^{\prime} \right],\end{array}\end{eqnarray} \tag{ 3 }$$

where r₀ is a fiducial lower radius that remains fixed in mass and ${r}_{\mathrm{tp}}$ is the outer turning point of the mode; we note that ${r}_{\mathrm{tp}}$ is the actual radius at which reflection occurs (i.e., $K({r}_{{r}_{\mathrm{tp}}})=0$). The second term in brackets on the rhs of Equation (3) is the equilibrium value of the phase, and the first term is its instantaneous value. The factor of 2 is to account for the fact that the wave propagates from r₀ to ${r}_{\mathrm{tp}}$, and then back to r₀. For long-period modes the reflection point is very near the base of the convection zone, but for shorter-period modes ${r}_{\mathrm{tp}}$ is deeper in the model.

As previously mentioned, the Doppler shift caused by reflection from the base of the moving convection zone may lead to even larger phase shifts (see Figure 7). Essentially, the frequency shift due to reflection leads to a change in the wavelength of the wave (Montgomery et al. 2019). Thus, in traveling down to the inner turning point and back the mode no longer travels an exact integer number of wavelengths, arriving back at the outer turning point with a phase shift. Formally, the frequency perturbation is given by

$$\begin{eqnarray}&&\omega ^{\prime} =\omega \left(1-{v}_{\mathrm{CZ}}/{v}_{\mathrm{gr}}\right),\end{eqnarray} \tag{ 4 }$$

or, setting $\omega ^{\prime} =\omega +{\rm{\Delta }}\omega$,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\omega }{\omega }=-\displaystyle \frac{{v}_{\mathrm{CZ}}}{{v}_{\mathrm{gr}}},\end{eqnarray} \tag{ 5 }$$

where ω is the frequency of the incident wave, v_CZ is the velocity of the base of the convection zone, and v_gr is the group velocity of the wave at the base of the convection zone, given by ${v}_{\mathrm{gr}}\equiv \partial \omega /\partial K$.

Equation (5) is nontrivial to evaluate since v_gr formally goes to zero near the base of the convection zone. To remove this difficulty, we imagine the reflection occurring from a surface of constant phase well beneath the outer turning point, ${r}_{\mathrm{tp}}$, and we calculate the velocity of this surface, v_surf. Using asymptotic formulae, we write this phase between the radii r and ${r}_{\mathrm{tp}}$ as

$$\begin{eqnarray}&&\phi \equiv {\int }_{r}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&={\int }_{r}^{{r}_{0}}K(r^{\prime} ){dr}^{\prime} +{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} ,\end{eqnarray} \tag{ 7 }$$

where r is the radius of the surface of constant phase, K(r) is given by Equation (2), and r₀ is a fiducial radius between r and ${r}_{\mathrm{tp}}$ whose value is held constant.

Setting ∂ϕ/∂t = 0 in Equation (7) and solving for ∂r/∂t yields

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial t}=\displaystyle \frac{1}{K(r)}\displaystyle \frac{\partial }{\partial t}{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} .\end{eqnarray} \tag{ 8 }$$

Since $\partial r/\partial t$ represents the velocity of the surface of constant phase, we set v_surf = ∂r/∂t. Using this in place of v_CZ in Equation (5), we find

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\omega }{\omega }=-\displaystyle \frac{{v}_{\mathrm{surf}}}{{v}_{\mathrm{gr}}}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&=\,-\displaystyle \frac{1}{{v}_{\mathrm{gr}}K(r)}\displaystyle \frac{\partial }{\partial t}{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} .\end{eqnarray} \tag{ 10 }$$

Since we take the point r to be safely beneath the outer turning point, we can use the asymptotic expression for the group velocity of g-modes, ${v}_{\mathrm{gr}}\approx -\omega /K(r)$. Substituting this into Equation (10) yields

$$\begin{eqnarray}&&{\rm{\Delta }}\omega =\displaystyle \frac{\partial }{\partial t}{\int }_{{r}_{0}}^{{r}_{\mathrm{tp}}}K(r^{\prime} ){dr}^{\prime} .\end{eqnarray} \tag{ 11 }$$

We recognize the integral on the rhs of Equation (11) as the mode phase between the point r₀ and the outer turning point. Since the time derivative only acts on the time-dependent variations of this quantity, we can rewrite it using ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ as defined in Equation (3). Thus,

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}\omega & = & \displaystyle \frac{\partial }{\partial t}\left[{\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})\right]\\ & \approx & \langle \omega \rangle \,{\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}}),\end{array}\end{eqnarray} \tag{ 12 }$$

where we have replaced the time derivative with the factor $\langle \omega \rangle =2\pi /\langle P\rangle$, where $\langle P\rangle$ is a characteristic pulsation timescale for the ensemble of excited modes in the star (e.g., shorter for hotter stars, longer for cooler stars). Based on the results of Mukadam et al. (2006), we take $\langle P\rangle =300$ s for ${T}_{\mathrm{eff}}={\rm{12,000}}$ K, $\langle P\rangle =800\,{\rm{s}}$ for ${T}_{\mathrm{eff}}$ = 11,500 K, and $\langle P\rangle =1000\,{\rm{s}}$ for ${T}_{\mathrm{eff}}$ ≤ 11,000 K.

We note that the large time-dependent changes occur near the base of the convection zone, so ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ and therefore Δω should be independent of the reference depth r₀. Numerically, we find that ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ changes by <0.01% as r₀ increases in depth by one mode wavelength. This is necessary since the calculated value of the frequency shift Δω should be independent of the assumed distance from the base of the convection zone at which the calculation is made. For the calculations in Section 10, we arbitrarily choose this fiducial radius to be the point below the outer turning point at which ϕ = 5π/4; this is slightly over half a mode wavelength beneath the outer turning point.

The frequency difference given by Equation (12), as the ray propagates down to the inner turning point and back out to the outer turning point, results in a total accumulated phase change of

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Delta }}{\phi }_{\mathrm{dop}} & = & 2\pi n\left(\displaystyle \frac{{\rm{\Delta }}\omega }{\omega }\right)\\ & = & 2\pi n\,\displaystyle \frac{\langle \omega \rangle }{\omega }\,{\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}}).\end{array}\end{eqnarray*}$$

These phase shifts lead to damping of the mode, as we show in the following section. Finally, using Equation (3) we see that the phase change due to the Doppler shift of the mode frequency can be related to that due to the changing size of the propagation cavity, i.e.,

$$\begin{eqnarray}&&{\rm{\Delta }}{\phi }_{\mathrm{dop}}=n\pi \,\displaystyle \frac{\langle \omega \rangle }{\omega }\,{\rm{\Delta }}{\phi }_{\mathrm{cav}}.\end{eqnarray} \tag{ 13 }$$

Since n ≥ 1 and $\omega \approx \langle \omega \rangle$, we expect that the phase shift (and therefore the damping rate) due to the Doppler shift to dominate that due to the changing cavity size.

We note here that this approach assumes that the traveling waves that comprise the mode propagate freely from the outer turning point down to the inner turning point and back again without encountering any regions of rapid spatial variation. Such regions could produce partial reflection and transmission of the waves, resulting in an unequal distribution of kinetic energy between the core and envelope regions. This "mode trapping" would lead to differential sensitivity of consecutive radial orders to phase shift effects. The fact that we ignore these possible internal reflections means that mode trapping effects will be suppressed in our damping calculations.

For the mode to be completely coherent, it needs to accumulate exactly 2πn radians of phase each time it propagates back and forth in the star. As shown in the previous section, a changing convection zone can upset this condition, leading to destructive interference and a change in amplitude given by

$$\begin{eqnarray}&&\displaystyle \frac{{dA}}{{dt}}=-\displaystyle \frac{A}{{nP}}\left(1-\cos {\rm{\Delta }}\phi \right)\end{eqnarray} \tag{ 14 }$$

(see Montgomery et al. 2015). Assuming that $A\propto {e}^{-\gamma t}$, and averaging γ over values of the phase shift from −Δϕ to +Δϕ gives

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{1}{{nP}}\left(1-\displaystyle \frac{\sin {\rm{\Delta }}\phi }{{\rm{\Delta }}\phi }\right).\end{eqnarray} \tag{ 15 }$$

This is the equation we use to calculate the damping rates of modes given the amplitudes (Δϕ) of their phase shifts.

The WD models used in these calculations were computed with the MESA stellar evolution code (Paxton et al. 2011, 2013, 2015, 2018, 2019), and all the models assumed a mass of 0.6 M_⊙ (the mean mass of WDs in our sample is 0.62 M_⊙). Since the depth of the convection zone is the key parameter in this study, we have used the mixing length calibration of Tremblay et al. (2015) to set this parameter. Specifically, we have used the value of ML2/α_Schwa from their Table 2 for the $\mathrm{log}\,g=8.0$ sequence interpolated to our ${T}_{\mathrm{eff}}$ values of interest.

To simulate the ${T}_{\mathrm{eff}}$ changes due to finite-amplitude pulsations, we compute models with ${T}_{\mathrm{eff}}$ values centered on the equilibrium state, e.g., ${T}_{\mathrm{eff}}$ = 12,000 ± 50 K. Since only the surface layers of the model correspond to these perturbed ${T}_{\mathrm{eff}}$ values, in a previous study (Montgomery et al. 2015) we spliced the outer portion of these models onto the equilibrium model. This was justified by the fact that the models rapidly converge with depth to nearly identical equilibrium structures so that the splicing produces no visible numerical artifacts. For the present study we adopt a different procedure. We first compute the depth of the convection zone with the perturbed ${T}_{\mathrm{eff}}$, say ${T}_{\mathrm{eff}}$ = 12,050 K. We then find the value of ML2/α that reproduces this depth for the equilibrium ${T}_{\mathrm{eff}}$ of 12,000 K. We now identify this model as the "12,050 K" model since it has a convection zone depth that is identical to that of a "true" 12,050 K model. This approach has the advantage that no splicing of models is required.

For long-period modes whose outer turning points are very near the base of the convection zone, all of our approaches yield very similar results. That is, computing ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ from either the difference in oscillation periods of two neighboring models (Equation (1)) or using Equation (3) to integrate the radial wavenumber in two models yields the same result. For instance, this can be seen in Figure 9(b), in the agreement between the two sets of solid and dotted curves for periods greater than ∼1000 s.

For shorter-period modes that have an outer turning point farther beneath the convection zone, ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ depends on how quickly the differences in model quantities decay with depth. Denoting the differences between the two models (at constant mass coordinate) by δ, δN²/N², δc²/c², and $\delta {\omega }_{c}^{2}/{\omega }_{c}^{2}$ initially decay exponentially with depth. However, instead of decaying to zero, they asymptote to a level of ∼10⁻⁶ (see inset in Figure 8). Since these models are designed to represent the same star at different times, the internal structure of the models (which should not be strongly affected by the changing surface convection zone) should be nearly identical. We force this to be true in our analysis by fitting the model differences with an exponential solution that decays to zero with increasing depth. We then add these differences back to the unperturbed model to create the new perturbed model, and these are the models used in integrating ${\rm{\Delta }}\phi ({r}_{0},{r}_{\mathrm{tp}})$ in Equation (3). This leads to a net decrease in the calculated damping rates of short-period modes, which were dominated by small, unphysical model differences in deeper layers.

Zoom In Zoom Out Reset image size

In an asymptotic treatment, modes with outer turning points sufficiently far beneath the perturbed region will be completely unaffected by changes in the convection zone, and so will have zero calculated damping. However, in a full calculation the modes are global, and while they are evanescent in the region beneath the convection zone, they do still sample it. For long-period modes that propagate near the base of the convection zone, the period difference ΔP between models is dominated by the difference in convection zone depths, so we calculate ΔP as a direct difference of the pulsation periods of corresponding modes in the two models (Montgomery et al. 2015). For short-period modes, such period differences are dominated by numerical artifacts from the deeper layers, so we instead calculate the period difference via a variational approach (see Equation (14.19) of Unno et al. 1989 and Equation (5.80) of Christensen-Dalsgaard 2014), i.e.,

$$\begin{eqnarray}&&{\rm{\Delta }}P=-\displaystyle \frac{\pi }{{\omega }^{3}}\,\displaystyle \frac{\int {{dM}}_{r}\,\delta {N}^{2}\,{\xi }_{r}^{2}}{\int {{dM}}_{r}\left[\,{\xi }_{r}^{2}+{\ell }({\ell }+1)\,{\xi }_{h}^{2}\,\right]},\end{eqnarray} \tag{ 16 }$$

where ξ_r and ξ_h are the radial and horizontal displacements of the spatial eigenfunction, and δN² is the difference in the Brunt–Väisälä frequency between the two models. We then use Equations (1), (12), (13), and (15) to calculate the phase shifts and damping associated with the period change of this mode.

Combining Equations (1), (3), (12), and (13) with Equation (15), we compute finite-amplitude damping rates for ℓ = 1 pulsation modes. In Figures 9(a), (b) we show the damping rates using 0.6 M_⊙ WD models at two different temperatures: ${T}_{\mathrm{eff}}$ = 12,000 K for Figure 9(a), with assumed temperature excursions of ±50 K; and ${T}_{\mathrm{eff}}$ = 11,000 K for Figure 9(b) with assumed temperature excursions of ±100 K. After geometric cancellation, these cases correspond to observed luminosity perturbations of ∼1% and ∼2%, respectively. The blue curves with circular points show the damping rate due to the changing size of the g-mode cavity, while the orange curves with square points give the damping rate due to the Doppler shifting of the reflected wave's frequency. In addition, curves with filled symbols employ a direct integration of the wavenumber (Equation (3)) while those with open symbols use the period difference (Equation (1)) with ΔP calculated as described in Section 9. As can be seen, the two methods agree well for long-period modes, e.g., the blue filled and unfilled symbols in Figure 9(b) having P > 1000 s. The red curve in both plots is an estimate of the linear growth rates of these modes. These estimates are based on calculations by Wu (1998), as shown in Figures 5.5–5.8 of her thesis (see also Figure 6 of Wu & Goldreich 1999 and Figure 5 of Wu & Goldreich 2001).

Zoom In Zoom Out Reset image size

The linear growth rates, by definition, do not depend on the amplitude of the pulsations, while the nonlinear damping mechanisms presented here do.⁶ As expected, we see that the damping due to the Doppler shift of the mode's frequency is much larger than that due to the variation in the size of its g-mode cavity, for both models and sets of modes. In Figure 9(a), this driving is larger than the assumed total damping for all periods that are driven, using either method of calculating the damping, while in Figure 9(b), the driving exceeds the damping for periods less than ∼650 s when the damping is calculated by direct integration of the radial wavenumber (filled points), with the driving exceeding the damping for periods less than ∼500 s if the period difference method is used (unfilled symbols).

We next examine these results both as a function of ℓ and ${T}_{\mathrm{eff}}$. In order to be conservative, in the remainder of this paper we use the damping estimates computed by direct integration of the radial wavenumber (e.g., filled symbols in Figures 9(a), (b)), since these appear to provide a lower limit for these rates. In Figures 10(a)–(d), we show the finite-amplitude total damping rates as computed in the previous section for ℓ = 1, 2, and 3 modes. The different panels give the results for different values of ${T}_{\mathrm{eff}}$; for context, we also plot estimates of the driving rates based on Figures 5.5–5.8 from the thesis of Wu (1998). In Figure 10(a), the range of ℓ = 1 modes assumed to be driven by convective driving is P ∼ 100–500 s (blue dashed curve), which is plausible for a ${T}_{\mathrm{eff}}$ = 12,000 K model. For this case, we see that the damping is essentially insignificant for the driven modes. The same is true for ℓ = 2 and 3 modes.

Zoom In Zoom Out Reset image size

Figure 10(b) is a slightly cooler case (${T}_{\mathrm{eff}}$ ∼ 11,500 K), and the range of linearly driven ℓ = 1 modes is assumed to be P ∼ 100–900 s, with ranges of 100–600 s and 100–450 s for ℓ = 2 and 3, respectively. For the assumed 200 K amplitude, the damping mechanism could now potentially affect the amplitudes of the longest periods that are driven: 800 s for ℓ = 1, and 500 s and 350 s for ℓ = 2 and 3.

Finally, Figures 10(c), (d) show the trends at cooler temperatures. As ${T}_{\mathrm{eff}}$ decreases, the driving decreases relative to the damping, implying smaller overall amplitudes. This has partially been accounted for by a decrease in the assumed amplitude of the pulsations, from ±200 K in Figure 10(b) to ±100 K in Figures 10(c), (d). For ℓ = 1, the period at which the driving equals the damping decreases from 800 s at ${T}_{\mathrm{eff}}$ = 11,500 K to 650 s and 550 s at ${T}_{\mathrm{eff}}$ = 11,000 K and 10,500 K, respectively. A similar trend can be seen for ℓ = 2 and 3 modes. Overall, these results imply that both the amplitudes seen in cooler models and their maximum periods should decrease. While the decrease in overall amplitude near the observed red edge has been previously noted (Mukadam et al. 2006), any observational decrease in maximum period with ${T}_{\mathrm{eff}}$ is less dramatic and has not been conclusively identified.

Due to the nature of the stochastic excitation mechanism, the theoretical damping rates of solar-like pulsators provide a prediction for the observed Fourier widths of the modes. While this connection is less clear for pulsators with modes that are linearly unstable (such as the DAVs), we will attempt to calculate mode widths in the context of the damping mechanism presented above.⁷

Equation (11) provides an estimate of the frequency change produced by a single interaction of a wave with the convection zone. If we interpret this frequency shift as an estimate of the width of the Fourier peak of the mode, we can calculate the expected widths for the pulsation modes in a given stellar model.

We show the results of such a calculation in Figure 11. The dashed magenta curve shows the frequency width for ℓ = 1 modes calculated in this way, while the dotted orange curve shows the width for ℓ = 2 modes. These calculations are based on the equilibrium model of Figure 10(b) (${T}_{\mathrm{eff}}$ = 11,500 K, $\mathrm{log}\,g=8.0$, $\mathrm{ML}2/\alpha =0.76$, with ${T}_{\mathrm{eff}}$ variations of ±200 K). We see that there is a steep rise in the ℓ = 1 predicted widths at ∼800 s, nominally coinciding with the observed rise seen in the data. We also note that the ℓ = 2 modes with P ∼ 700 s would not be expected to have large widths if they were ℓ = 1 modes, but are easily accounted for as ℓ = 2 modes. While certainly not the final word, this calculation provides a partial explanation for the observed rise in mode widths with increasing period.

Zoom In Zoom Out Reset image size

In the preceding sections we have shown that the modes with observed Fourier widths greater than 0.3 μHz have longer periods (P ≳ 800 s), and that in most cases these modes are calculated to propagate all the way to the base of the surface convection zone in pulsating DA WDs. We have proposed mechanisms through which the convection zone can cause a lack of coherence in these modes. We find the dominant mechanism to be the Doppler shift of the mode frequency as it reflects off the time-dependent outer turning point. Preliminary results are that this mechanism is stronger near the observed red edge of the instability strip and may be important for the observed reduction in mode amplitudes there.

The phase shifts are a finite-amplitude effect that naturally lead to damping. For hotter models (${T}_{\mathrm{eff}}$ ≳ 12,000 K), this damping is quite small and likely unimportant. On the other hand, for cooler models (${T}_{\mathrm{eff}}$ < 12,000 K), this damping can be significant for longer-period modes. The DAVs known to have stable pulsations coherent over a time period of years have short-period pulsations (P ≲ 400 s) and are near the blue edge of the instability strip. From Figure 9(a), we see that the proposed damping mechanism is quite small for these modes, which could be part of the reason these modes exhibit such extreme coherence. For cooler models, the damping is still small for the short-period modes, so the possibility remains that cooler pulsators could also harbor very stable pulsations in the form of low-period modes (P < 300 s). These modes could have propagation regions sufficiently distant from the convection zone that they would still have extreme stability. Such modes, if they exist, could be used to asteroseismically trace secular evolution of cooler DAVs or expand the search for planets around white dwarfs with pulsation timing variations. As far as we are aware, no systematic search for such stable modes in cooler DAVs has been made, though some short-period modes do appear stable in cool DAVs over the span of months in Kepler/K2 observations.

We note that the calculations presented here could be improved in two respects. First, the convective driving rates based on the formalism of either Wu & Goldreich (1999) or Dupret et al. (2004), and Van Grootel et al. (2012) could be calculated for the stellar models employed here, leading to a consistent set of driving and nonlinear damping rates. In addition, the calculation of the difference in structure of neighboring models could be improved. Although the differences in their outer layers are dominated by their different convection zone depths, which our current calculations take into account, the differences deeper in the models will be directly due to the pulsation modes themselves. Thus, realistic eigenfunctions should be used to calculate the instantaneous structure of these models, which are then used in the subsequent analysis.

Finally, many DAVs observed with the Kepler and K2 missions undergo outbursts, increases in average brightness of 10%–40% that can typically last from 5 to 15 hr (Bell et al. 2015, 2017; Hermes et al. 2015). While the best-known theory for this process involves a resonant transfer of energy from a driven parent mode to damped daughter modes as an amplitude threshold is reached (Wu & Goldreich 2001; Luan & Goldreich 2018), we speculate that the mechanism we have proposed involving phase shifts of reflected modes could be relevant as well. It is possible that some pulsators reach an amplitude threshold in which there is a slight increase in the damping, and this increased damping leads to a slight heating of the surface layers. This in turn causes a net thinning of the convection zone, leading to larger phase shift mismatches, which leads to further damping, and the cycle reinforces itself. In addition, as the surface convection zone becomes thinner, radiative damping may play an important role in removing energy from the high overtone modes (Luan & Goldreich 2018).

Nonlinear resonant energy transfer also appears to modify the frequency and amplitudes of modes on timescales shorter than expected from secular evolution, and this has been observed in a number of pulsating white dwarfs with long-baseline observations (e.g., Dalessio et al. 2013; Zong et al. 2016). However, the frequency changes of these effects are considerably smaller than those observed in the broadened mode line widths discussed here. Additionally, many of the DAVs with frequency changes incompatible with secular evolution occur in hotter stars with short-period (<300 s) pulsations that should be relatively unaffected by the convection zone (e.g., Hermes et al. 2013).

We present a mechanism that may be relevant for the limitation of pulsation amplitudes in white dwarf stars for modes above a threshold period. As the convection zone changes depth during the pulsation cycle, the condition for coherent reflection of the outgoing traveling wave is slightly violated. In effect, this causes the amplitude of the mode (viewed as the superposition of inward- and outward-propagating components) to decrease, leading to damping. This mechanism should be present at some level in all pulsating WDs, and should be larger near the red edge of the DAV instability strip. In addition, this mechanism could possibly be relevant for other g-mode pulsators with surface convection zones (e.g., Gamma Doradus stars), or even large-amplitude pulsators such as high-amplitude delta Scuti stars (HADs) or RR Lyrae stars.

M.H.M. and D.E.W. acknowledge support from the United States Department of Energy under grant DE-SC0010623 and the NSF grant AST 1707419. M.H.M., D.E.W., and B.H.D. acknowledge support from the Wootton Center for Astrophysical Plasma Properties under the United States Department of Energy collaborative agreement DE-NA0003843. J.J.H. acknowledges support from NASA K2 Cycle 5 grant 80NSSC18K0387 and K2 Cycle 6 grant 80NSSC19K0162. K.J.B. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1903828.

Software: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), Scipy (Virtanen et al. 2019), Numpy (Oliphant 2006), Matplotlib (Hunter 2007).