From Stretched String to GYRE

The numerical technique demonstrated in the The Stretched String Problem section provides a powerful analog to how GYRE solves the oscillation equations. The full gory details of GYRE’s approach are laid out in Townsend & Teitler (2013); in this section we briefly summarize it, highlighting similarities and differences with the stretched-string problem.


Similar to the stretched-string problem, GYRE begins by separating variables in space and time. For the radial displacement perturbation \(\xir\), trial solutions take the form

\[\xir(r,\theta,\phi;t) = \operatorname{Re} \left[ \sqrt{4\pi} \, \txir(r) \, Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii \sigma t) \right],\]

In addition to the same sinusoidal time dependence as in eqn. (3), an additional spherical harmonic term \(Y^{m}_{\ell}\) appears because we are separating in three (spherical) spatial coordinates rather than one.


As with the stretched-string problem, GYRE then discretizes the governing equations on a spatial grid \(\{x_{1},x_{2},\ldots,x_{N}\}\). However, a couple of important differences arise at this juncture. First, the system of oscillation equations are fourth order (sixth, in the non-adiabatic case). Rather than employing finte-difference approximations to high-order differential operators, GYRE instead decomposes the problem into set of coupled first-order equations. This set is written generically as

\[x \deriv{\vty}{x} = \mA \, \vty\]

where \(\vty\) is a vector of \(\neq\) dependent variables, and \(\mA\) is a \(\neq \times \neq\) Jacobian matrix. In the adiabatic case, \(\neq=4\); in the non-adiabatic case, \(\neq=6\).

Second, while this equation could be discretized using a simple finite-difference approximation to the left-hand side, GYRE offers more-sophisticated approaches with higher orders of accuracy. These include the Magnus schemes described in Townsend & Teitler (2013), and implicit Runge-Kutta schemes mentioned in Townsend et al. (2018). The choice of scheme is set by the diff_scheme parameter of the &num namelist group. The discretization leads to difference equations of the form

\[\vty_{k+1} = \mY_{k+1;k} \, \vty_{k},\]

relating the dependent variable vector at adjacent grid points. The \(\neq \times \neq\) fundamental solution matrix \(\mY_{k+1,k}\) is evaluated from the value(s) of \(\mA\) within the interval \([x_{k},x_{k+1}]\) using the discretization scheme.

There are \(N-1\) of these sets of difference equations. They are augmented with the boundary conditions

\[\subin{\mB} \, \vty_{1} = 0, \qquad\qquad \subout{\mB} \, \vty_{N} = 0,\]

where \(\subin{\mB}\) is a \(\nin \times \neq\) matrix representing the \(\nin\) inner boundary conditions, and \(\subout{\mB}\) is a \(\nout \times \neq\) matrix representing the outer boundary conditions (note that \(\nin + \nout = \neq\)). Together, the difference equations and boundary conditions comprise a linear system of \(\neq\,N\) algebraic equations and \(\neq N\) unknowns.

Linear System

The linear system can be written in the same form (4) as with the stretched-string problem. However, now \(\vu\) is the vector with components

\[\begin{split} \vu = \begin{pmatrix} \vty_{1} \\ \vty_{2} \\ \vdots \\ \vty_{N-1} \\ \vty_{N} \end{pmatrix}\end{split}\]

and the system matrix \(\mS\) is an \(\neq N \times \neq N\) block-staircase matrix with components

\[\begin{split}\mS = \begin{pmatrix} \subin{\mB} & \mz & \cdots & \mz & \mz \\ -\mY_{2;1} & \mI & \cdots & \mz & \mz \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mz & \mz & \cdots & -\mY_{N;N-1} & \mI \\ \mz & \mz & \cdots & \mz & \subout{\mB} \end{pmatrix}.\end{split}\]

As before, the linear system (4) has non-trivial solutions only when the determinant of \(\mS\) vanishes. Thus, GYRE finds eigenvalues of the oscillation equation by solving the characteristic equation

\[\Dfunc(\omega) \equiv \det(\mS) = 0,\]

where the dimensionless frequency

\[\omega \equiv \sqrt{\frac{R^{3}}{GM}} \, \sigma,\]

is the product of the star’s dynamical timescale and the oscillation frequency \(\sigma\). (Internally, GYRE works extensively with such dimensionless quantities, as it improves the stability of the numerical algorithms).

Scanning for Eigenfrequencies

In the adiabatic case, GYRE searches for roots of the discriminant function \(\Dfunc\) using the same bracketing and refinement strategies as the stretched-string problem.

In the non-adiabatic case, a complication is that the discriminant function and the dimensionless frequency are both complex quantities. Solving the characteristic equation in the complex plane is computationally challenging because there is no equivalent to bracketing and refinement. GYRE implements a couple of different approaches to the problem, as discussed in the Non-Adiabatic Calculations chapter.