Dispersion Relations in Asteroseismology: a neat trick

(carrying on from the last time)

For the last five years, I have worked on asteroseismology — the analysis and interpretation of oscillations of stars. My initial introduction to this topic was via a draft version of the book my current advisor was in the process of writing, which was targeted at astronomers looking to enter the field. Since I did not actually have a formal background in astronomy at the time, I found it quite blasé in its approach to various approximations made in operationalising many of the theoretical aspects of the problem. Unfortunately, this turned out to be quite counterproductive: I have wasted many hours trying to figure out why certain approximations were made, and others were not, and why particular choices may or may not work. It turns out that not everyone does this, and so interrogating the foundations of certain pedagogical elisions can lead to novel research problems — which in my opinion should not have been low-hanging fruit for a graduate student like myself in the first place.

For context, the qualitative phenomena of asteroseismology are constructed around an approximate dispersion relation of the form

where

• $k_r$ and $k_h$ are the radial and tangential components of the wavevector (working here in spherical coordinates): in particular, $k_h^2 = {l (l + 1) \over r^2}$ for normal modes of degree $l$.
• $c_s^2 = {\Gamma_1 P \over \rho}$ is the squared adiabatic sound speed: here $\Gamma_1 = \left.{\partial \log P \over \partial \log \rho}\right|_S$ is a property of the equation of state.
• $N^2$ is the squared Brunt-Väisälä frequency, about the pedagogy of which I have complained quite enough. For our purposes we will use the identities $N^2 = g \left({1 \over \Gamma_1 P}{\mathrm d P \over \mathrm d r} - {\mathrm d \log \rho \over \mathrm d r}\right) \equiv -g\left({g \over c_s^2} + {\mathrm d \log \rho \over \mathrm d r}\right).$

In particular for high frequencies ($\omega^2 > c_s^2 k_h^2 \gg N^2$) one approximately recovers the familiar dispersion for acoustic waves, $k_r^2 + k_h^2 = \omega^2 / c_s^2$, while for low frequencies ($\omega^2 < N^2 \ll c_s^2 k_h^2$) we obtain buoyancy waves: $|k|^2 \propto k_h^2$.

The derivation of this specific dispersion relation has always seemed a little iffy to me, since most pedagogical treatments of it demand certain terms to be neglected with (what seems to me to be) not very much justification. More frustratingly, there are actual waves that defy description by this specific dispersion relation; in typical expositions these are handled as ad-hoc special cases, which I find quite unsatisfactory. On a technical level, I also find this unsatisfactory because it is unclear how one might generalise when the Cowling approximation is relaxed — that is to say, when we do not neglect perturbations to the gravitational potential. Nonetheless, this seems to be good enough for astronomers in general; most of the technical literature on the theoretical aspects of the problem was written by a handful of people between the 80's to the early 90's. While this was very competently executed (for a technical audience), much of the subsequent pedagogical exposition has merely recapitulated the same few (significantly simplified) points over and over, albeit with increasing levels of polish. This being said, I really do think there is significant room here for improvement in the current pedagogy. Consider these a sketch for lecture notes which I might make in the future.

A Particular Linear System

Waves in a spherically symmetric self-gravitating fluid under hydrostatic equilibrium obey this system of coupled differential equations:

The quantities $(\xi_r, P', \Phi')$ are the actual dynamical variables of the system, while the last line is an algebraic constraint relating tangential part of the Lagrangian fluid perturbation wavefunction to these dynamical variables. So far, none of this is novel: this is pretty standard material in textbooks about stellar oscillations. Note that since the third line is a second-order differential equation in $\Phi'$, we can treat $\mathrm d \Phi' \over \mathrm d r$ as an additional dynamical quantity $\Psi$, thereby yielding a set of four coupled first-order linear ODEs: the first three lines of Eq. 2, and the constraint $\Psi = {\mathrm d \Phi' \over \mathrm d r}$.

I now choose a new set of dynamical variables $(y_1, y_2, y_3, y_4)$, which are related to $(\xi_r, P', \Phi', \Psi)$ only by judicious choice of integrating factors. I do this so that each line of Eq. 2 after this change of variables contains terms proportional to either a derivative of a each dynamical quantity, or to that dynamical quantity itself, but never both. In particular:

• We choose $y_1 = e^{h_1} \xi_r$ so that $\mathrm {d \over d r}(e^{h_1} \xi_r) = \left({1 \over r^2}{\mathrm d \over \mathrm d r} (r^2 \xi_r) - {g \over c_s^2}\xi_r\right)e^{h_1} \equiv e^{h_1}\left({\mathrm d \over \mathrm d r}\xi_r + \xi_r h_1'\right)$; this gives us $e^{h_1} = r^2 \exp\left[-\int {g \over c_s^2} \mathrm d r\right]$.
• We choose $y_2 = e^{h_2} {P' \over \rho}$ so that ${\mathrm d \over \mathrm d r}(e^{h_2} {P' \over \rho}) = \left({1 \over \rho}{\mathrm d P' \over \mathrm d r} + {g \over \rho c_s^2}P'\right)e^{h_2}$; this gives us $h_2 = -\int {N^2 \over g} \mathrm d r$.
• From the third line of Eq. 2 it is clear that we should choose $y_3 = r^2 \Psi \equiv e^{h_3} \Psi$.
• We have $y_4 = \Phi'$, needing no integrating factor (i.e. $h_4 = 0$).

Note that the choices of $y_1$ and $y_2$ are identical to the dynamical variables $(\hat{\xi}, \hat{\eta})$ used by Unno et al. (1989) in the Cowling approximation. The choices of $y_3, y_4$ are fairly natural. To illustrate why we chose such a set of variables, we rewrite the system of ODEs in explicit linear form:

In particular, this is a system of linear differential equations of the form $\mathbf{E}^{-1} {\mathrm {d} \over \mathrm d r} \mathbf{y} = \mathbf{A E}^{-1} \mathbf{y}$, where $\mathbf{E}$ is a gauge transformation relating these dynamical quantities $y_i$ to our original set of variables $(\xi_r, P', \Phi', \Psi)\ \hat{=}\ \mathbf{x}$ as $\mathbf{y} = \mathbf{E x}$. Recall that for a system of ODEs of the form $\dot{\mathbf{x}} = \mathbf{Lx}$, such a gauge transformation yields $\dot{\mathbf{y}} = \left(\mathbf{E L E}^{-1} - \mathbf{E} \dot{\mathbf{E}^{-1}}\right)\mathbf{y}$. We see now that we have chosen $\mathbf{E}$ so that the diagonal elements of the matrix $\mathbf{A}$ are all zero.

To turn this into a dispersion relation, we must now introduce some physical restrictions on the solutions to these ODEs. Let us suppose that all of these $y_i$ are associated with the same semiclassical radial wavenumber $k_r$, in that $y_i \sim e^{i \int k_r \mathrm d r} \iff {\mathrm d \over \mathrm d r} y_i \sim i k_r y_i$. This kind of approximation holds well for very large $k_r$, in the same regime as the JWKB approximation (although perhaps a little less rigorously): either the $y_i$ are highly oscillatory (real $k_r$), or they are decaying/increasing very rapidly (imaginary $k_r$). Then we find that locally, $k_r$ must be an eigenvalue of the matrix $\mathbf{A}$ (which is a property of this kind of semiclassical approximation). These eigenvalues may be found as the roots of the characteristic equation

where I have eliminated the degree $l$ in favour of the tangential wavenumber $k_h$. Since $\mathbf{A}$ is a $4\times4$ matrix, it must be a quartic polynomial in $k_r$. Explicitly, we find after some algebra that

Had we restricted our attention to the Cowling approximation (thereby getting rid of $y_3$ and $y_4$), we would have obtained only the upper left $2\times2$ submatrix. Its determinant is in braces on the left-hand-side of Eq. 3; setting it to zero yields precisely the Cowling-approximation dispersion relation of Eq. 1. The terms on the right-hand-side of Eq. 3 are all proportional to $4\pi G \rho$; this also appears on the left-hand-side as a small correction to the acoustic cutoff frequency from accounting for gravitational perturbations, thereby leaving the Cowling approximation.

Note that we now have imaginary coefficients in this dispersion relation (on the right-hand-side). This implies that the roots $k_r$ are generally speaking complex. The interpretation of this is that this dispersion relation might not allow waves of constant amplitude, except where either $N^2$ vanishes, which is the case in the convective regions of the star, or where $k_r \ll {k_h^2 g \over \omega}$, which is the case where $r \to 0$. These happen to correspond to the propagation regions of classical p- and g-modes anyway.

The usual way in which this works, for real $k_r$, is that we first fix a frequency $\omega$, and then use this dispersion relation (which implicitly defines $k_r$ as a function of $\omega$ and other quantities in the stellar structure) to solve for some wavefunction that satisfies the requisite boundary conditions on the dynamical variables. This turns out to be an overdetermined problem, where such solutions only exist for certain $\omega$: the normal modes of oscillation. Techniques for this kind of problem have been covered quite extensively in existing textbooks.

In this case, we can also go the other way around: we can first choose purely imaginary $k_r$ so that the left-hand-side vanishes: in particular we choose $k_r = -i k_h$, so that the dynamical quantities all fall off exponentially as we head inwards away from the surface. This only satisfies the dispersion relation for specific frequencies $\omega_0$, such that the terms on the right-hand-side also vanish. We find that

this is precisely the dispersion relation for $f$-modes, which do not emerge from the usual asymptotic analysis! The conjugate choice for $k_r$ — with the wavefunction decaying exponentially towards the surface — yields $\omega_0^2 = - g k_h$, which does not correspond to physical oscillations. Note, however, that strictly speaking Eq. 4 cannot be used to find normal modes: on the left-hand-side we have a constant, while on the right-hand-side, $g k_h \in \mathcal{O}(1)$ only when the stratification is such that ${\partial m(r) \over \partial r} \sim r^2$. This is approximately true in the near-surface layers, so this approximation becomes increasingly good at high $k_h$ (i.e. high $l$), where the wavefunctions, which then go as $y_i \sim \exp\left[k_h r\right]$, die off increasingly rapidly away from the surface.

Concluding remarks

In the usual asymptotic analysis, the standard procedure is to get rid of factors of $1 \over r$, $g \over c_s^2$ and so forth, with the excuse of such approximations holding only at large radii (so that ${1 \over r} \to 0$) or near the core (so that $g \to 0$). Alternatively, terms are made to vanish by demanding from the very beginning that $k_r \gg {g \over c_s^2}$, for example. However, this has always struck me as being somewhat inconsistent — for example, why do we get rid of terms proportional to $1 \over r$, but retain $k_h^2 \propto {1 \over r^2}$? By contrast, we have played no such games here. That is to say, no such asymptotic foolery is actually necessary to obtain expressions like Eq. 1 and Eq. 3!

The price we pay is that this kind of analysis is strictly speaking not uniquely defined; had we chosen a different set of dynamical variables, we would have obtained a different dispersion relation (because the elements of the matrix $\mathbf{A}$ would have been different). While we have chosen dynamical variables with a well-defined procedure (eliminating the diagonal entries of $\mathbf{A}$), there is no a priori reason why we should have done so — although I would argue that this procedure is still better-defined and mathematically more defensible than the standard pedagogy.

For second-order systems, we can reduce the second-order ODE for each variable to canonical Schrödinger form (which gives us a unique set of dynamical coordinates); this is a necesssary first step for subsequent JWKB analysis. This is why there is a uniquely-defined acoustic cutoff frequency under the Cowling approximation. However, I am not aware of a general procedure for higher-order systems of ODEs.