Explaining Deviations from the Scaling Relationship of the Large Frequency Separation

In the first-order asymptotic theory of p-modes, the large frequency separation Δν emerges from the phase integral quantization condition

$$\begin{eqnarray}&&{\int }_{{r}_{1}}^{{r}_{2}}{k}_{r}(r,\omega ){dr}=\pi \left(n+\kappa \right)\end{eqnarray} \tag{ 4 }$$

for integer n, where k_r is the WKB wave number, r₁ and r₂ are the turning points of the integral (usually where k_r vanishes), and κ is a function of frequency. Here κ depends only on the functional behavior of k_r² at the turning points. In particular, it can be shown that $\kappa =-\tfrac{1}{2}$ is constant if ${k}_{r}^{2}\propto r-{r}_{t}$ at both turning points (Gough 2007; Aerts et al. 2010). To a first approximation, this expression yields the Duvall law,

$$\begin{eqnarray}&&{\int }_{{r}_{t}}^{R}\displaystyle \frac{{dr}}{{c}_{s}}\sqrt{1-\displaystyle \frac{{S}_{l}^{2}}{{\omega }^{2}}}=\pi \left(\displaystyle \frac{n+\kappa }{\omega }\right),\end{eqnarray} \tag{ 5 }$$

where c_s is the sound speed, ${S}_{l}^{2}=\tfrac{l(l+1)}{{c}_{s}^{2}{r}^{2}}$ is the Lamb frequency, and $\omega =2\pi \nu$ is the angular frequency. Again, the inner turning point of this integral is defined by where the integrand vanishes (or r = 0 for radial modes), and the outer turning point R is the same as that used to evaluate the sound-travel time.

Naturally, higher-order approximations are possible. For instance, Tassoul (1994) derived Equation (1) from a second-order expansion of the asteroseismic equations of motion, using linear combinations of spherical Bessel functions and their derivatives as ansatzen for the eigenfunctions, with a similar phase integral quantization condition emerging for the arguments of these functions. This was in turn taken to fourth order by Roxburgh & Vorontsov (1994) with a similar approach.

On the other hand, the accuracy of the WKB expression can also be improved with more detailed asymptotic analysis (i.e., not setting terms to zero prematurely before performing the WKB analysis). For instance, a more accurate description is afforded by Deubner & Gough (1984),

$$\begin{eqnarray}&&{\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\sqrt{1-\displaystyle \frac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}-\displaystyle \frac{{S}_{l}^{2}}{{\omega }^{2}}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right)}=\pi \left(\displaystyle \frac{n+\kappa ^{\prime} }{\omega }\right),\end{eqnarray} \tag{ 6 }$$

which emerges from just such a more detailed analysis (with still more detail provided in Gough 1993). The integrand in this expression reduces asymptotically to that in Equation (5) far into the interior of the convective envelope of solar-like main-sequence stars, where ${\omega }_{\mathrm{ac}}\ll \omega$ for most p-modes of interest. Here N² is the Brunt–Väisälä or buoyancy frequency, and κ' is another function of frequency (different, in general, from κ in Equation (5)). However, κ' also only depends on the behavior of the integrand in Equation (6) at the turning points.

The formulation of Equation (6) differs from the standard form of the Duvall law in terms of both its explicit dependence on an acoustic cutoff frequency

$$\begin{eqnarray}&&{\omega }_{{\rm{a}}{\rm{c}}}^{2}=\displaystyle \frac{{c}_{s}^{2}}{4{H}^{2}}\left(1-2\displaystyle \frac{dH}{dr}\right)\end{eqnarray} \tag{ 7 }$$

involving $H=-{\left(\tfrac{d\mathrm{ln}\rho }{{dr}}\right)}^{-1}$ (the density scale height) in the integrand and the locations of the turning points (as the integrands vanish at different places). In particular, the outer turning point is essentially determined by the acoustic cutoff frequency, rather than being defined by the boundary conditions of the eigenvalue problem. The precise form of the acoustic cutoff frequency and other quantities involved in the oscillation equations will depend on the choice of dynamical variable used to perform the WKB analysis. Nonetheless, the resulting eigenvalue equation can typically be put into the form of Equation (6), although different expressions for ${\omega }_{\mathrm{ac}}$, S_l², and N² may have to be used in place of those described here. An illustrative example can be found in Gough (1993).

With respect to the above formalism, Christensen-Dalsgaard & Perez Hernandez (1992) instead define a phase function α so that

$$\begin{eqnarray}&&{\int }_{{r}_{t}}^{R}\displaystyle \frac{{dr}}{{c}_{s}}\sqrt{1-\displaystyle \frac{{S}_{l}^{2}}{{\omega }^{2}}}=\pi \left(\displaystyle \frac{n+\alpha }{\omega }\right),\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&\alpha (\omega )=\displaystyle \frac{\omega }{\pi }\left({T}_{1}-{T}_{2}\right)+\kappa ^{\prime} ,\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{T}_{1} & = & {\int }_{{r}_{t}}^{R}\displaystyle \frac{{dr}}{{c}_{s}}\sqrt{1-\displaystyle \frac{{S}_{l}^{2}}{{\omega }^{2}}},\\ {T}_{2} & = & {\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\sqrt{1-\displaystyle \frac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}-\displaystyle \frac{{S}_{l}^{2}}{{\omega }^{2}}\left(1-\displaystyle \frac{{N}^{2}}{{\omega }^{2}}\right)}\end{array}\end{eqnarray} \tag{ 10 }$$

being integrals with dimensions of time appearing on the left-hand side of Equations (5) and (6). This permits the use of Equation (1) with α in place of κ. Whereas the phase functions κ and κ' defined previously depend only locally on the asymptotic properties of the WKB integrands in the neighborhoods of their classical turning points, α as defined in this manner also depends on the global properties of the entire acoustic mode cavity.

We take the first-order expression (Equation (6)) at face value and assume that the quantities N² and S_l² have no parametric dependence on the frequency, and further that all quantities on the right-hand side are continuous functions of the frequency ν. We proceed to expand finite differences as Taylor series. In particular, we consider finite differences of the terms in Equation (6), multiplied throughout by $\omega /\pi =2\nu$ and evaluated at frequencies of adjacent n for the same l,

$$\begin{eqnarray}&&\begin{array}{l}2{\nu }_{n+1,l}{T}_{2}({\nu }_{n+1,l})-2{\nu }_{n-1,l}{T}_{2}({\nu }_{n-1,l})\\ \quad \quad \sim 2+\kappa ^{\prime} \left({\nu }_{n+1,l}\right)-\kappa ^{\prime} \left({\nu }_{n-1,l}\right)\\ \quad \Longrightarrow \,4\displaystyle \frac{d\nu {T}_{2}}{d\nu }{\rm{\Delta }}\nu \sim 2+2\displaystyle \frac{d\kappa ^{\prime} }{d\nu }{\rm{\Delta }}\nu ,\end{array}\end{eqnarray} \tag{ 11 }$$

with the quantities on the last line being evaluated at the frequency $\nu =\tfrac{1}{2}({\nu }_{n+1,l}+{\nu }_{n-1,l})$ so that the error term is ${ \mathcal O }({\rm{\Delta }}{\nu }^{3})$, and where T₂ is the left-hand side of Equation (6). This yields the expression

$$\begin{eqnarray}&&{\rm{\Delta }}\nu \sim {\left(2\displaystyle \frac{d\nu {T}_{2}}{d\nu }-\displaystyle \frac{d\kappa ^{\prime} }{d\nu }\right)}^{-1}.\end{eqnarray} \tag{ 12 }$$

For the first term, we observe that

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{d\nu }\left(\nu {\int }_{{r}_{1}}^{{r}_{2}}\sqrt{f(r,\nu )}\,\displaystyle \frac{{dr}}{{c}_{s}}\right)=\displaystyle \frac{d}{d\omega }\left(\omega {\int }_{{r}_{1}}^{{r}_{2}}\sqrt{f(r,\omega )}\,\displaystyle \frac{{dr}}{{c}_{s}}\right)\\ \quad ={\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{f(r,\omega )+\tfrac{1}{2}\omega \tfrac{\partial }{\partial \omega }f(r,\omega )}{\sqrt{f(r,\omega )}}\,\displaystyle \frac{{dr}}{{c}_{s}}\\ \quad +\,\omega \left(\displaystyle \frac{\sqrt{f({r}_{2},\omega )}}{{c}_{s}}\displaystyle \frac{{{dr}}_{2}}{{\rm{d}}\omega }-\displaystyle \frac{\sqrt{f({r}_{1},\omega )}}{{c}_{s}}\displaystyle \frac{{{dr}}_{1}}{d\omega }\right).\end{array}\end{eqnarray} \tag{ 13 }$$

Since the turning points of the integral are either fixed (for ${r}_{1}=0$) or defined to be where the relevant integrand vanishes, the boundary contributions to the above derivative also vanish, leaving us with

$$\begin{eqnarray}&&2{\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\displaystyle \frac{1-\tfrac{{S}_{l}^{2}{N}^{2}}{{\omega }^{4}}}{\sqrt{1-\tfrac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}-\tfrac{{S}_{l}^{2}}{{\omega }^{2}}\left(1-\tfrac{{N}^{2}}{{\omega }^{2}}\right)}}-\displaystyle \frac{d\kappa ^{\prime} }{d\nu }\end{eqnarray} \tag{ 14 }$$

inside the parentheses in Equation (12). For radial modes in particular, we can neglect S_l², leaving us with

$$\begin{eqnarray}&&{\rm{\Delta }}\nu \sim {\left(2{\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\displaystyle \frac{1}{\sqrt{1-\tfrac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}}}-\displaystyle \frac{d\kappa ^{\prime} }{d\nu }\right)}^{-1}.\end{eqnarray} \tag{ 15 }$$

As noted previously, $\kappa ^{\prime} =-\tfrac{1}{2}$ in the ideal WKB scenario of the mode cavity exhibiting asymptotically linear behavior near the classical turning points. On the other hand, if the radicand tends to a constant value near the inner turning point (as would be the case for Sun-like stars, where ${\omega }_{\mathrm{ac}}\ll \omega$ far into the interior), then the behavior of the eigenfunctions near the inner turning point is described by the spherical Bessel function j₀ and $\kappa ^{\prime} =-\tfrac{1}{4}$ (Gough 1993). In either case, assuming that it changes very slowly with respect to frequency, its derivative is much smaller than the integral expression. We thus propose the use of a modified integral estimator for the large frequency separation:

$$\begin{eqnarray}&&{\rm{\Delta }}\nu \sim {\left(2{\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\displaystyle \frac{1}{\sqrt{1-\tfrac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}}}\right)}^{-1}.\end{eqnarray} \tag{ 16 }$$

This expression contains an explicit dependence on the frequency. Observationally, however, we are limited to the study of only modes with frequencies near ν_max, the frequency of maximum acoustic power. Moving forward, we will evaluate this quantity at ν_max.

We note that the denominator of the integrand in Equation (16) vanishes at the outer boundary of the domain of integration, so the integrand becomes singular. Depending on the acoustic cutoff frequency of the stellar model, this may also be the case at the inner boundary. In the ideal WKB regime, this results in the presence of integrable singularities at the end points. We would therefore expect that the evolution of these turning points plays a significant role in how the value of this estimator changes over the course of stellar evolution. Our expression does not contain a prescription for the actual turning points of the asymptotic integral, only that they depend on the formulation of the acoustic cutoff frequency that accompanies the degree of approximation being used.

Moreover, within the domain of integration, we note that the integrand is strictly larger than the $1/{c}_{s}$ that appears as the integrand of the sound-travel time. Since the integrand is also strictly positive, we conclude that the integral is strictly larger than the sound-travel time. Given that the sound-travel time generally overestimates Δν (i.e., yields too small a value of the integral, as in Figure 2), this is suggestive as at least a zeroth-order mark of improvement.

To compare the results obtained from the integral and frequencies, we compute the logarithms of the ratios between the value of ${\rm{\Delta }}{\nu }_{\mathrm{fit}}$ (obtained by a least-squares fit of ${\nu }_{n,l=0}$ against n) versus the values of our integral estimator in Equation (16), which we plot as circles in subsequent figures. We also do this with the log ratios of ${\rm{\Delta }}{\nu }_{\mathrm{fit}}$ versus the values of the sound-travel time estimator $1/2{T}_{0}$ (which we show as triangles), as well as versus the values predicted by the scaling relation (squares). As a more quantitative (but still heuristic) comparison, we compute the root mean square (rms) deviation of these log ratios from zero. An rms deviation of zero means that the estimator exactly coincides with ${\rm{\Delta }}{\nu }_{\mathrm{fit}}$ for all models.

The top panel of Figure 3 shows the values of the log ratios described above for stellar models along the main sequence and on the subgiant branch. It is visually evident that the agreement between the estimator in Equation (16) and the fitted values of Δν is considerably better than both the scaling relation and the sound-travel time. We find that the scaling relation results in an rms log deviation of about 0.010 dex, compared to 0.025 dex for the sound-travel time and 0.002 dex for our integral estimator.

Zoom In Zoom Out Reset image size

Whereas the relative deviations from the scaling relation appear to exhibit curlicues for main-sequence stars of different masses (as in Figure 1), we note that the most obvious deviations between our asymptotic estimator and the fitted value of Δν occur as a sharp spike near the main-sequence turnoff during the onset of shell burning (at the core hydrogen exhaustion/subgiant hook for more massive stars). We examine this more closely in Figure 4.

We show in the top panel of Figure 4 the acoustic cutoff frequency, in units of ν_max, of main-sequence, main-sequence turnoff, and subgiant stellar models with $M=1.0{M}_{\odot }$. Following our above discussion of various values of κ in the WKB regime, we see that for the main-sequence star (blue solid curve), the inner turning point is at r = 0, where ${\nu }_{\mathrm{ac}}/{\nu }_{\max }\ll 1$, and the WKB wave function goes as ${j}_{0}(\omega t)$, where $t\,\sim {\int }_{0}^{r}\sqrt{1-{\left({\omega }_{\mathrm{ac}}/\omega \right)}^{2}}\tfrac{{dr}}{{c}_{s}}$, with the same integrand as in Equation (6). On the other hand, for the subgiant star (green dashed curve), the WKB wave function goes as $\mathrm{Ai}(-| \omega t{| }^{\displaystyle \frac{3}{2}})$ near a first-order classical turning point at $r={r}_{1}$, where ${\nu }_{\mathrm{ac}}({r}_{1})\sim {\nu }_{\max }$ and $t\sim {\int }_{{r}_{1}}^{r}\sqrt{1-{\left({\omega }_{\mathrm{ac}}/\omega \right)}^{2}}\tfrac{{dr}}{{c}_{s}}$, and $\mathrm{Ai}$ is the Airy function.

These represent two distinct asymptotic regimes in which the assumptions underlying the WKB approximation hold well. However, this is not necessarily true during the transition between these two regimes. In particular, for some stellar models near the main-sequence turnoff (orange dotted curve), we have a maximum of the acoustic cutoff frequency with a value very close to ν_max. This violates one of the underlying assumptions of the WKB approximation (viz., that k_r varies slowly except near classical turning points, which requires in this case that $\nu \gg {\nu }_{\mathrm{ac}}$), and neither of the above asymptotic descriptions is a good approximation for the actual mode eigenfunctions.

We noted earlier that the radicand in Equation (16) vanishes where ${\nu }_{\mathrm{ac}}={\nu }_{\max }$. For first-order turning points, this results in integrable singularities at the end points of integration. However, when there is a maximum point where ${\nu }_{\mathrm{ac}}/{\nu }_{\max }\,\sim 1$ (as would be the case with a second-order turning point), either the integrand changes very rapidly (becoming very large) in a region smaller than the wavelength of the eigenfunction ($\max {\nu }_{\mathrm{ac}}/{\nu }_{\max }\lessgtr 1$), or a nonintegrable singularity is introduced into the domain of integration ($\max {\nu }_{\mathrm{ac}}/{\nu }_{\max }=1$), and the convergence of the numerical integration is poor. In these cases, a numerical computation of Equation (16) will be very different from the fitted value of Δν: the orange dotted curve of Figure 4 corresponds to the maximum deviation between Equation (16) and the fitted value of Δν along our $1{M}_{\odot }$ track (as in the orange point in the bottom panel). In principle, the phase function κ, as we have defined it, also becomes discontinuous in the neighborhood of $\nu =\max {\nu }_{\mathrm{ac}}\sim {\nu }_{\max }$, and its derivative also becomes singular. A proper accounting should permit these two singularities to cancel each other out. However, this is difficult to handle numerically.

It is possible that these difficulties can be resolved by a more refined analysis. For instance, strictly speaking, the asymptotic behavior of the WKB wave function near such a second-order stationary point should be given by the parabolic cylinder function (Gough 2007). However, such an analysis is fairly involved and may not as readily yield a simple (and easily computed) expression like Equation (16). Alternatively, one might imagine choosing a set of dynamical variables such that the expression for the cutoff frequency is always singular near the center of the model (Gough 1993), so that this transition point does not emerge. Again, doing so would require a more detailed analysis than we have performed, which we leave beyond the scope of this paper.

3.1.1. RGB and Red Clump Stars

In the lower two panels of Figure 3, we show the log ratios of frequency-separation estimators for ascending red giant branch (RGB; middle panel) and descending RGB and red clump stars (lower panel). Once again, our expression in Equation (16) deviates much less from the fitted Δν than do both the scaling relation prediction and the sound-travel time. However, we note the emergence of terrace-like features in these deviation plots, which ultimately originate from the method by which we compute a fitted value for Δν. In particular, since we have included all modes with frequencies lower than the maximum acoustic cutoff frequency in the least-squares fit, and since the acoustic cutoff frequency changes relative to Δν over the course of stellar evolution, the number of modes used in the fit changes also changes discontinuously over the course of the track, as illustrated in Figure 5. These terraces are therefore a property of our benchmarking methodology, rather than being a feature of our asymptotic estimator.

Zoom In Zoom Out Reset image size

In addition, we observe that while the scaling relation overpredicts the large frequency separation for ascending RGB stars, it slightly underpredicts it for red clump stars. Since the use of the scaling relation is (implicitly) a homology argument (Belkacem et al. 2013), this indicates that the acoustic structure of red clump stars differs significantly from that of red giant stars. On this basis, Miglio et al. (2012) proposed the use of different correction factors for the scaling relation for first-ascent red giant and red clump stars, precisely to account for these structural differences.

In terms of our formulation, we show in Figure 6 the radial dependence of the acoustic cutoff frequency for two stellar models—one red giant and one red clump—with identical masses and radii. While their acoustic structure is very similar in the outer parts of the star, we see that their modified radial mode cavities have different inner turning points. Accordingly, our integral estimator correctly returns different estimates for Δν for these different types of stars.

Zoom In Zoom Out Reset image size

3.1.2. Dependence on Model Atmospheres

The models previously discussed were constructed with Eddington-gray atmospheres. To examine the dependence of our integral estimator on atmospheric conditions, we also constructed models using Krishna–Swamy model atmospheres, with Y₀ and α_MLT calibrated separately. Once again, we show the log ratios between our estimator and the fitted Δν in the top panel of Figure 7. Surprisingly, there does not appear to be any significant difference in the structure of the residuals. We interpret this to indicate that the remaining deviations that persist in either case stem from either an insufficiently high order of approximation or issues with our formulation that become significant only in the interior of the star, rather than the atmosphere.

Zoom In Zoom Out Reset image size

One possible such inadequacy in our formulation could be our numerical treatment of the acoustic cutoff frequency. For example, Gough (1993) constructed a different expression for the acoustic cutoff frequency that takes the sphericity of the star into account, which only becomes significant deep in the stellar interior. This is also, therefore, where we expect the expression that we have used to be the most inaccurate. Since a large contribution to the value of our integral expression comes from singularities (integrable or otherwise) at the turning points, it is possible that the remaining error could potentially be considerably reduced or eliminated entirely by using an expression that is more accurate in the interior.

3.1.3. Dependence on Metallicity

3.1.4. Δν from Nonradial Modes

Finally, if we relax our previous restriction to radial modes, we arrive at a similar integral expression,

$$\begin{eqnarray}&&{\rm{\Delta }}{\nu }_{l}\sim {\left(2{\int }_{{r}_{1}}^{{r}_{2}}\displaystyle \frac{{dr}}{{c}_{s}}\displaystyle \frac{1-\tfrac{{S}_{l}^{2}{N}^{2}}{{\omega }^{4}}}{\sqrt{1-\tfrac{{\omega }_{\mathrm{ac}}^{2}}{{\omega }^{2}}-\tfrac{{S}_{l}^{2}}{{\omega }^{2}}\left(1-\tfrac{{N}^{2}}{{\omega }^{2}}\right)}}\right)}^{-1},\end{eqnarray} \tag{ 17 }$$

for the large separation of modes with $l\ne 0$. To test this expression, we evaluate

$$\begin{eqnarray}&&{\rm{\Delta }}{\nu }_{l}\sim {\nu }_{n+1,l}-{\nu }_{n,l}\end{eqnarray} \tag{ 18 }$$

by least-squares fitting of ${\nu }_{n,l}$ against n for l = 1 and 2 on the same set of stellar models as used in the first panel of Figure 3, again including all modes with frequencies below the maximum acoustic cutoff frequency. Again, we plot the log ratio of these against the corresponding integral estimates (this time from Equation (17)) in Figure 8. We see that for main-sequence stars, our estimator also adequately reproduces the fitted value of Δν_l. However, there is a considerable amount of scatter for stars that have evolved off the main sequence, as well as a clear systematic deviation for l = 2.

Zoom In Zoom Out Reset image size

We attribute the above to difficulties in constructing a well-defined average Δν for nonradial modes with evolved stars. For l = 1, avoided crossings emerge as the star evolves through the subgiant phase, and mixed-mode propagation (i.e., evanescent coupling to a g-mode cavity) becomes possible (Scuflaire 1974; Osaki 1975; Aizenman et al. 1977; Deheuvels & Michel 2011), with further mode splitting emerging when Δν becomes sufficiently small. This g-mode coupling also changes the frequencies of observed l = 2 mixed modes relative to pure p-modes. Mixed modes are described by two quantum numbers, n_p (the radial quantum number n of the p-mode being coupled) and n_g (Unno et al. 1989); in computing Δν_l, we have selected for each value of n_p the frequency of the corresponding mode with the lowest inertia. However, this does not necessarily correspond to the eigenfrequency of a pure p-mode with that value of n_p.

We illustrate this problem in Figure 9, where we plot the ratio of the least-squares fitted values of ${\rm{\Delta }}{\nu }_{l}$ to ${\rm{\Delta }}{\nu }_{0}$ for l = 1, 2. It is clear that this ratio is relatively constant (although not necessarily unity) near the main sequence but shows large variations—both systematic and in terms of scatter between timesteps—for evolved stars. Much of the scatter in Figure 8 is a result of these rapid variations in Δν_l, rather than originating from our modified estimator. While our formulation cannot explain the remaining systematic differences for evolved stars (as in the bottom panel of Figure 8, for l = 2), we nonetheless note that these have qualitatively the same morphology as the relative differences between the sound-travel time and the fitted value of Δν. Again, this likely indicates the limitations of our numerical formulation of the acoustic cutoff frequency; a careful accounting (as in Gough 1993) shows that it, too, depends on the degree l (in a somewhat complicated fashion), which is not included in our numerical integration.

Zoom In Zoom Out Reset image size