The GYRE Stellar Oscillation Code

GYRE is a stellar oscillation code. Given an input stellar model, GYRE calculates the eigenfrequencies and eigenfunctions for the normal oscillation modes of the model. These data can be put to a variety of uses; the most common is to compare them against observed oscillation frequencies of a star, allowing constraints on the star’s fundamental parameters (mass, radius, etc.) to be established — the discipline of asteroseismology.

GYRE also supports other, related kinds of calculation. One example is evaluating the response of a star to tidal disturbances produced by an orbiting companion; because this is an instance of forced stellar oscillations, similar numerical techniques can be brought to bear.

Preliminaries

Intended Audience

This manual is aimed at a broad audience — whether you’re a GYRE novice or a seasoned veteran, it provides you with the information you’ll need to get the most out of GYRE. However, it does presume some experience with Unix command-line environments, and likewise some basic familiarity with the subject of stellar oscillations. If you need the former, then the Internet is your oyster; and for the latter, we recommend the following online resources:

• Jørgen Christensen-Dalsgaard’s Lecture Notes on Stellar Oscillations;

• Gerald Handler’s Asteroseismology article.

Obtaining GYRE

The source code for GYRE is hosted in the https://github.com/rhdtownsend/gyre git repository on GitHub. GYRE is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, version 3.

Citing GYRE

If you use GYRE in your research, please cite one or more of the relevant ‘instrument’ papers:

• Townsend & Teitler (2013) describes the basic operation of the code;

• Townsend et al. (2018) outlines capabilities for non-adiabatic oscillation calculations;

• Goldstein & Townsend (2020) describes the contour method for finding non-adiabatic modes;

• Sun et al. (2023) introduces the gyre_tides frontend for evaluating tidal responses.

If you find yourself using GYRE on a regular basis, you might consider contributing to the project to ensure its long-term success. Options include

• contributing code to the project (e.g., via GitHub pull requests), to extend GYRE’s capabilities;

• contributing documentation and tutorials to the project, to make GYRE more user-friendly;

• inviting the GYRE team to be co-authors on relevant papers;

• inviting the GYRE team to be co-investigators on relevant grant applications.

Development Team

GYRE remains under active development by the following team:

• Rich Townsend (University of Wisconsin-Madison); project leader

• Warrick Ball (University of Birmingham)

• Earl Bellinger (MPIA Garching)

• Zhao Guo (Cambridge University)

• Mathias Michielsen (KU-Leuven)

• Joel Ong (University of Hawaii-Manoa)

• Jarret Rosenberg (University of Wisconsin-Madison)

• Meng Sun (Northwestern University)

• Vincent Vanlaer (KU-Leuven)

Former developers include:

• Jacqueline Goldstein (MIT)

• Aaron Lopez

Also, the following people have made valuable contributions toward testing GYRE:

• Siemen Burssens (KU Leuven)

• Timothy Van Reeth (KU Leven)

Related Links

• The GYRE discussion forums, the place to post feature requests and bug reports (don’t send emails!).

• The MESA Software Development Kit (SDK), which provides the compilers and supporting libraries needed to build GYRE.

• The MESA Stellar Evolution Code, which can generate stellar models readable by GYRE.

Acknowledgments

GYRE has been developed with financial support from the following grants:

• NSF awards AST-0908688, AST-0904607, ACI-1339606, ACI-1663696, and AST-1716436;

• NASA awards NNX14AB55G, NNX16AB97G, and 80NSSC20K0515.

GYRE has also benefited greatly from contributions (code, bug reports, feature requests) from the academic community. Thanks, folks!

Quick Start

To get started with GYRE, follow these five simple steps:

• install the MESA Software Development Kit;

• download the GYRE source code;

• unpack the source code using the tar utility;

• set the GYRE_DIR environment variable to point to the newly created source directory;

• compile GYRE using the command make -C $GYRE_DIR.

For a more in-depth installation guide that covers alternative use-cases, refer to the Installation chapter. If the code doesn’t compile properly, consult the Troubleshooting chapter. Otherwise, proceed to the next chapter where you’ll put together your first GYRE calculation.

Example Walkthrough

This chapter provides a walkthrough of a example GYRE project, to illustrate the typical steps involved. For this example, we’ll use gyre (the frontend focused on stellar oscillations) to find eigenfrequencies and eigenfunctions of dipole and quadrupole gravity modes for a MESA model of slowly pulsating B (SPB) star.

Making a Place to Work

When starting a new project, it’s a good idea to create a dedicated work directory to contain the various input and output files that gyre operates on. These commands will make a new directory beneath your home directory with the name work, and then set this directory as the current working directory:

mkdir ~/work
cd ~/work

Grabbing a Stellar Model

The next step is to grab the stellar model. There are a number of example models provided in the $GYRE_DIR/models directory; the following commands will copy a MESA model for a \(5\,\Msun\) SPB star into your work directory:

cp $GYRE_DIR/models/mesa/spb/spb.mesa .

Assembling a Namelist File

Now comes the fun part: assembling an input file containing the various parameters which control a gyre run. Using a text editor, create the file gyre.in in your work directory with the following content cut-and-pasted in:

&constants
/

&model
  model_type = 'EVOL'  ! Obtain stellar structure from an evolutionary model
  file = 'spb.mesa'    ! File name of the evolutionary model
  file_format = 'MESA' ! File format of the evolutionary model
/

&mode
  l = 1 ! Harmonic degree
/

&mode
  l = 2 ! Harmonic degree
/

&osc
  outer_bound = 'VACUUM' ! Assume the density vanishes at the stellar surface
/

&rot
/

&num
  diff_scheme = 'COLLOC_GL4' ! 4th-order collocation scheme for difference equations
/

&scan
  grid_type = 'INVERSE' ! Scan grid uniform in inverse frequency
  freq_min = 0.5        ! Minimum frequency to scan from
  freq_max = 1.5        ! Maximum frequency to scan to
  n_freq = 100          ! Number of frequency points in scan
/

&grid
  w_osc = 10 ! Oscillatory region weight parameter
  w_exp = 2  ! Exponential region weight parameter
  w_ctr = 10 ! Central region weight parameter
/


&ad_output
  summary_file = 'summary.h5'                         ! File name for summary file
  summary_item_list = 'l,n_pg,freq,freq_units,E_norm' ! Items to appear in summary file
  detail_template = 'detail.l%l.n%n.h5'        	      ! File name template for detail files
  detail_item_list = 'l,n_pg,omega,x,xi_r,
                      xi_h,c_1,As,V_2,Gamma_1' 	      ! Items to appear in detail files
  freq_units = 'CYC_PER_DAY'                   	      ! Units of freq output items
/

&nad_output
/

This file is in namelist format, containing multiple namelist groups. Detailed information on the groups can be found in the Namelist Input Files chapter; for now, let’s just focus on some of the more-important aspects of the file above:

• the &constants namelist group is empty, telling gyre to use default values for fundamental constants;

• the &model namelist group tells gyre to read an evolutionary model, in MESA format, from the file spb.mesa;

• the two &mode namelist groups tells gyre to search first for dipole (\(\ell=1\)) and then quadrupole (\(\ell=2\)) modes;

• the &osc namelist group tells gyre to assume, when setting up the outer boundary conditions in the oscillation equations, that the density vanishes at the stellar surface;

• the &scan namelist group tells gyre to scan a region of dimensionless angular frequency space typically occupied by gravity modes;

• the &grid namelist group tells gyre how to refine the model spatial grid;

• the &ad_output namelist group tells gyre what adiabatic data to write to which output files; summary data to the file summary.h5, and individual mode data to files having the prefix mode.;

• the &nad_output namelist group is empty, telling gyre not to write any non-adiabatic data.

Running gyre

With the hard work done, it’s now trivial to run gyre:

$GYRE_DIR/bin/gyre gyre.in

As the frontend runs (on multiple cores, if you have a multi-core machine; see the FAQ for more details), it will print lots of data to the screen. Let’s break down this output, chunk by chunk.

First, gyre prints out its version number, tells us (in OpenMP threads) how many cores it is running on, and indicates which file it is reading parameters from (here, gyre.in):

gyre [master]
-------------

OpenMP Threads   : 4
Input filename   : gyre.in

Next, gyre loads the stellar model from the file spb.mesa. This model comprises 1814 points and extends from the surface all the way to the center (which is why gyre decides not to add a central point).

Model Init
----------

Reading from MESA file
   File name spb.mesa
   File version 1.00
   Read 1814 points
   No need to add central point

gyre then prepares to search for modes with harmonic degree \(\ell=1\) and azimuthal order \(m=0\) (not specified in gyre.in, but assumed by default), by building a frequency grid and a spatial grid:

Mode Search
-----------

Mode parameters
   l : 1
   m : 0

Building frequency grid (REAL axis)
   added scan interval :  0.5000E+00 ->  0.1500E+01 (100 points, INVERSE)

Building spatial grid
   Scaffold grid from model
   Refined 0 subinterval(s) in iteration 1
   Final grid has 1 segment(s) and 1814 point(s):
      Segment 1 : x range 0.0000 -> 1.0000 (1 -> 1814)

(The concepts of spatial and frequency grids are explored in greater detail in the Numerical Methods and Understanding Grids chapters). Next, gyre attempts to bracket roots of the discriminant function (again, see the Numerical Methods chapter) by searching for changes in its sign:

Starting search (adiabatic)

Evaluating discriminant
  Time elapsed :     0.886 s

Finally, for each bracket found gyre uses a root solver to converge to the eigenfrequency. Each row of output here corresponds to a mode that gyre has successfully found:

Root Solving
   l    m    n_pg    n_p    n_g       Re(omega)       Im(omega)        chi n_iter
   1    0      -9      0      9  0.50907836E+00  0.00000000E+00 0.6678E-14      7
   1    0      -8      0      8  0.58398491E+00  0.00000000E+00 0.1352E-13      6
   1    0      -7      0      7  0.66078111E+00  0.00000000E+00 0.1666E-13      7
   1    0      -6      0      6  0.73734087E+00  0.00000000E+00 0.3141E-13      6
   1    0      -5      0      5  0.89820448E+00  0.00000000E+00 0.1363E-13      6
   1    0      -4      0      4  0.11322842E+01  0.00000000E+00 0.6270E-13      7
   1    0      -3      0      3  0.13377876E+01  0.00000000E+00 0.6789E-13      6
  Time elapsed :      0.382 s

The columns appearing are as follows:

l

harmonic degree \(\ell\)

m

azimuthal order \(m\)

n_pg

radial order \(n\) (in the Eckart-Osaki-Scuflaire-Takata scheme)

n_p

acoustic-wave winding number \(n_{\rm p}\)

n_g

gravity-wave winding number \(n_{\rm g}\)

Re(omega)

real part of dimensionless eigenfrequency \(\omega\)

Im(omega)

imaginary part of dimensionless eigenfrequency \(\omega\) (zero here because we’ve performed an adiabatic calculation)

chi

convergence parameter

n_iter

number of iterations required for convergence

These values are printed to screen primarily to give an idea of gyre’s progress. Some things to watch out for:

• The convergence parameter chi, defined as the ratio of discriminant values before and after the root finding, should small (on the order of 1E-9 to 1E-15). If it is significantly larger than this, the mode may not be properly converged; and if it is significantly smaller than this, there may be numerical issues with the discretization scheme.

• The number of iterations n_iter should be moderate; values above 20 or so indicate that gyre is having problems converging.

• The mode radial order n_pg should be monotonic-increasing. Departures from this behavior can happen for a number of reasons, that are discussed in the Troubleshooting chapter.

After processing the dipole modes, gyre repeats the search steps for the quadrupole modes. Once the overall run is complete, a number of output files are written:

• A summary file with the name summary.h5

• For each mode found, a detail file with the name detail.lL.nN.h5, where L and N are the harmonic degree and radial order of the mode, respectively.

The Output Files chapter discusses how to read and analyze these files.

Frontends

This chapter summarizes the frontends provided by GYRE — the executables programs that users run to perform calculations. Although the Example Walkthrough chapter focuses on the gyre frontend, there are others available focused on different kinds of task.

gyre

The gyre frontend calculates the free-oscillation modes of a stellar model. The general flow of execution is outlined in the chart to the right. After reading the namelist input file and the model, gyre loops over &mode namelist groups, processing each in turn.

For a given group, gyre searches over a range of oscillation frequencies for modes with a specific harmonic degree \(\ell\) and azimuthal order \(m\). With each mode found, the eigenfrequency, eigenfunctions and other data are optionally written to a detail file. At the end of the run, response data from all modes found (across all &mode groups) are optionally written to a summary file.

The table below lists which namelist groups, and in what number, should appear in namelist input files for gyre.

Description	Namelist group name	Count
Constants	&constants	1
Grid Parameters	&grid	\(\geq 1\)[1]
Mode Parameters	&mode	\(\geq 1\)
Stellar Model Parameters	&model	1
Numerical Parameters	&num	\(\geq 1\)[1]
Oscillation Parameters	&osc	\(\geq 1\)[1]
Output Parameters	&ad_output	1
	&nad_output	1
Rotation Parameters	&rot	\(\geq 1\)[1]
Frequency Scan Parameters	&scan	\(\geq 1\)

Footnotes

[1] (1,2,3,4)

While the input file can contain one or more of the indicated namelist group, only the last (tag-matching) one is used.

gyre_tides

The gyre_tides frontend calculates the response of a stellar model to tidal forcing by a orbiting point-mass companion. The general flow of execution is outlined in the chart to the right. After reading the namelist input file and the model, gyre_tides loops over &tide namelist groups, processing each in turn.

For a given group, gyre_tides solves for the response of the star to the superposition of partial tidal potentials \(\PhiTlmk\) (see the Tidal Effects section). The response wavefunctions and other data associated with an individual partial potential are optionally written to a detail file. At the end of the run, response data from all partial potentials (across all &tide groups) are optionally written to a summary file.

The table below lists which namelist groups, and in what number, should appear in namelist input files for gyre_tides.

Description	Namelist group name	Number
Constants	&constants	1
Grid Parameters	&grid	\(\geq 1\)[1]
Stellar Model Parameters	&model	1
Numerical Parameters	&num	\(\geq 1\)[1]
Orbital Parameters	&orbit	\(\geq 1\)[1]
Oscillation Parameters	&osc	\(\geq 1\)[1]
Output Parameters	&tides_output	1
Rotation Parameters	&rot	\(\geq 1\)[1]
Tidal Parameters	&tide	\(\geq 1\)

Footnotes

[1] (1,2,3,4,5)

While the input file can contain one or more of the indicated namelist group, only the last (tag-matching) one is used.

Numerical Methods

This chapter explains the numerical methods used by the gyre frontend to solve the oscillation equations (similar approaches are followed in other frontends). Although it aims to be user-friendly, gyre is nevertheless a complex piece of software; thus, getting it to produce the ‘best’ results requires some degree of insight into the algorithms it uses to calculate mode eigenfrequencies and eigenfunctions.

The Stretched String Problem

We’ll start our discussion of numerical methods by considering the problem of finding normal-mode eigenfrequencies and eigenfunctions for waves on a stretched string clamped at both ends. Let the string have mass per unit length \(\rho\) and tension \(T\); then, the wave equation describing the transverse string displacement \(y(x,t)\) at spatial position \(x\) and time \(t\) is

\[\npderiv{y}{x}{2} = \frac{1}{c^{2}} \npderiv{y}{t}{2},\]

with \(c \equiv (T/\rho)^{1/2}\). If the string is clamped at \(x=0\) and \(x=L\), then the wave equation together with the boundary conditions

\[y(0,t) = 0 \qquad y(L,t) = 0\]

comprises a two-point boundary value problem (BVP).

Analytic Solution

The stretched-string BVP is straightforward to solve analytically. General solutions of the wave equation take the form of traveling waves,

\[y(x,t) = A \exp [\ii (k x - \sigma t) ],\]

where \(A\) an arbitrary constant, and the frequency \(\sigma\) and wavenumber \(k\) are linked by the dispersion relation

\[\sigma^{2} = c^{2} k^{2}.\]

The phase velocity of these waves is \(\sigma/k = \pm c\).

To satisfy the boundary condition at \(x=0\), we combine traveling-wave solutions with opposite-sign wavenumbers

\[y(x,t) = A \exp [\ii (k x - \sigma t) ] - A \exp [\ii (- k x - \sigma t) ] = B \sin(k x) \exp ( - \ii \sigma t),\]

where \(B = 2A\). For the boundary condition at \(x=L\) to be satisfied simultaneously,

\[\sin(k L) = 0,\]

and so

\[k L = n \pi\]

where \(n\) is a non-zero integer (we exclude \(n=0\) because it corresponds to the trivial solution \(y(x,t)=0\)). Combining this with the dispersion relation, we find that the normal-mode eigenfrequencies of the stretched-string BVP are

()\[\sigma = n \frac{\pi c}{L},\]

and the corresponding eigenfunctions are

()\[y_{n}(x,t) = B \sin \left( \frac{n \pi x}{L} \right) \exp ( - \ii \sigma t).\]

The index \(n\) uniquely labels the modes, and \(y_{n}(x,t)\) exhibits \(n-1\) nodes in the open interval \(x \in (0,L)\).

Separation

Now let’s see how we might go about solving the stretched-string BVP numerically. We begin by performing a separation of variables on the wave equation, assuming trial solutions of the form

()\[y(x;t) = \tilde{y}(x) \, \exp (-\ii \sigma t),\]

where \(\tilde{y}(x)\) is a function of \(x\) alone. Then, the wave equation reduces to an ordinary differential equation (ODE) for \(\tilde{y}\),

\[\nderiv{\tilde{y}}{x}{2} = - \frac{\sigma^{2}}{c^{2}} \tilde{y}.\]

Discretization

To solve the ODE, we discretize it to establish a set of difference equations. The discretization involves transforming the continuous function \(\tilde{y}(x)\) into a finite set of \(N\) values \(\{\tilde{y}_{1},\tilde{y}_{2},\ldots,\tilde{y}_{N}\}\), representing the function sampled on the discrete spatial grid \(\{x_{1},x_{2},\ldots,x_{N}\}\).

For simplicity let’s assume the grid is uniform, so that

\[x_{j+1} - x_{j} \equiv \Delta x = \frac{L}{N-1} \qquad (1 \leq j \leq N-1).\]

Then, the second derivative of \(\tilde{y}\) can be approximated to second order in \(\Delta x\) as

\[\left. \nderiv{\tilde{y}}{x}{2} \right|_{x=x_{j}} \approx \frac{\tilde{y}_{j+1} - 2 \tilde{y}_{j} + \tilde{y}_{j-1}}{\Delta x^{2}} \qquad (2 \leq j \leq N-1).\]

This allows us to replace the ODE with \(N-2\) difference equations

\[\frac{\tilde{y}_{j+1} - 2 \tilde{y}_{j} + \tilde{y}_{j-1}}{\Delta x^{2}} = - \frac{\sigma^{2}}{c^{2}} \tilde{y}_{j} \qquad (2 \leq j \leq N-1).\]

Together with the two boundary conditions

\[\tilde{y}_{1} = 0 \qquad \tilde{y}_{N} = 0,\]

we thus have a linear system of \(N\) algebraic equations and \(N\) unknowns.

Linear System

To find solutions to the linear system, we first write it in matrix form as

()\[\mS \vu = \mathbf{0},\]

where \(\vu\) is the vector with components

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

and the ‘system matrix’ \(\mS\) is an \(N \times N\) tridiagonal matrix with components

\[\begin{split}\mS = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 1 & \sigma^{2} \tau^{2} - 2 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & \sigma^{2} \tau^{2} - 2 & 1 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \end{pmatrix}.\end{split}\]

Here we’ve introduced

\[\tau \equiv \frac{\Delta x} c\]

as the sound crossing time of a single cell.

Equation (4) is a homogeneous linear system, meaning that it has non-trivial solutions \(\vu\) only when the determinant of \(\mS\) vanishes. With this in mind, we formulate the characteristic equation for the BVP,

\[\Dfunc(\sigma) = 0\]

where \(\Dfunc(\sigma) \equiv \det(\mS)\) is a discriminant function whose roots are the characteristic frequencies (eigenfrequencies) of the stretched-string BVP.

Plot of the discriminant function \(\Dfunc(\sigma)\) as a function of the frequency \(\sigma\), for the stretched-string BVP with \(N=50\). The orange dots highlight where \(\Dfunc=0\). The function has been scaled so that \(\Dfunc(0) = 1\). (Source)

Fig. 1 plots the discriminant function for the BVP discretized on a spatial grid of \(N=50\) points. The roots (zeros) of the function are highlighted by the orange markers; they fall very close to the values \(\sigma = \pi c/L, 2 \pi c/L, \ldots\) predicted by the analytic solutions.

Scanning for Eigenfrequencies

While Fig. 1 is useful for visualizing \(\Dfunc\), it’s not the best way to find eigenfrequencies. Instead, we can rely on well-established techniques for isolating and refining roots of monovariate functions.

First, we evaluate a set of \(M\) values \(\{\Dfunc_{1},\Dfunc_{2},\ldots,\Dfunc_{M}\}\), representing the discriminant function sampled on the discrete frequency grid \(\{\sigma_{1},\sigma_{2},\ldots,\sigma_{M}\}\). Then, we scan through these data looking for sign changes between adjacent discriminant values. If \(\Dfunc_{i} \Dfunc_{i+1} < 0\), we know that a root of the discriminant function must lie in the interval \((\sigma_{i},\sigma_{i+1})\) — we have bracketed a root. Fig. 2 illustrates the process of bracket scanning for a frequency grid comprising \(M=32\) points, distributed uniformly in \(\sigma\) across the same range as plotted in Fig. 1. This figure highlights five brackets containing the five roots identified previously.

Plot of the discriminant values \(\{\Dfunc\}\) on the discrete frequency grid \(\{\sigma\}\) (distributed uniformly in \(\sigma\)), for the stretched-string BVP with \(N=50\) and \(M=32\). The orange-haloed segments highlight adjacent points that bracket a root \(\Dfunc=0\). (Source)

Once a bracket is established for a given root, it can be narrowed through a process of iterative refinement until the root is converged upon. There are a variety of well-known root-finding algorithms that perform this refinement; the bisection method is conceptually the simplest, but approaches such as Brent’s method can be much more efficient. For the brackets plotted in Fig. 2, Table 1 compares the eigenfrequencies found using Python’s scipy.optimize.brentq() function, against the analytic values predicted by eqn. (1).

Numerical and analytic eigenfrequencies, in units of \(\pi c/L\), for the stretched-string BVP with \(N=50\). (Source)

n	numerical	analytic
1	0.999829	1.000000
2	1.998630	2.000000
3	2.995378	3.000000
4	3.989047	4.000000
5	4.978618	5.000000

Eigenfunction Reconstruction

For each of the eigenfrequencies found, we reconstruct the corresponding eigenfunction by solving the linear system (4). Because \(\det(\mS)\) is now zero, this system is guaranteed to have a non-trivial solution. The solution vector \(\vu\) resides in the null space of \(\mS\), and we can use standard numerical techniques (e.g., singular value decomposition) to evaluate it. Then, the \(j\)’th element of \(\vu\) corresponds to the eigenfunction sampled at the \(j\)’th spatial grid point:

\[(\vu)_{j} = \tilde{y}_{j} \equiv \tilde{y}(x_{j})\]

Plot of the eigenfunctions \(\tilde{y}\) as a function of spatial coordinate \(x\), for the first three modes of the stretched-string BVP with \(N=50\). The discrete points show the numerical functions, and the solid lines the corresponding analytic functions. In all cases, the eigenfunctions have been normalized to have a maximum \(|\tilde{y}|\) of unity. (Source)

Fig. 3 plots the eigenfunctions found in this way for the first three modes (\(n=1,2,3\)) given in Table 1. Also shown are the corresponding analytic solutions given by eqn. (2). The agreement between the two is good.

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 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.

Separation

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]\]

(this is taken from the Separated Equations section). In addition to the same sinusoidal time dependence as in eqn. (3), a spherical harmonic term \(Y^{m}_{\ell}\) appears because we are separating in three (spherical) spatial coordinates rather than one.

Discretization

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

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

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

Second, while the above equation system can 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_{j+1} = \mY_{j+1;j} \, \vty_{j},\]

relating the dependent variable vector at adjacent grid points. The \(\neqn \times \neqn\) fundamental solution matrix \(\mY_{j+1,j}\) is evaluated from the value(s) of \(\mA\) within the interval \([x_{j},x_{j+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 \neqn\) matrix representing the \(\nin\) inner boundary conditions, and \(\subout{\mB}\) is a \(\nout \times \neqn\) matrix representing the outer boundary conditions (note that \(\nin + \nout = \neqn\)). Together, the difference equations and boundary conditions comprise a linear system of \(\neqn\,N\) algebraic equations and \(\neqn N\) unknowns.

Linear System

The linear system can be written in the same form (cf. eqn. 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 \neqn 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 Oscillations section.

Limitations of the Numerical Method

The numerical method used both in the stretched string problem and in GYRE generally performs very well; however, it has a couple of failure scenarios that are important to understand. These scenarios arise through poor choices of the spatial grid used to discretize the governing differential equation(s), and/or the frequency grid used to search for roots of the discriminant function.

Insufficient Spatial Resolution

The cost of evaluating the determinant of the system matrix \(\mS\) scales proportionally to the number of grid points \(N\) used for the discretization. Therefore, in the interests of computational efficiency, we want to make \(N\) as small as possible.

However, things go wrong when \(N\) becomes too small. Fig. 4 demonstrates this by plotting the discriminant function for the stretched-string BVP with \(N=7\). Compared against Fig. 1, we see that toward larger \(\sigma\) the roots of the discriminant function become progressively shifted toward lower frequencies; and, above \(\sigma \approx 3.5 \pi c/L\), they disappear altogether.

Plot of the discriminant function \(\Dfunc(\sigma)\) as a function of the frequency \(\sigma\), for the stretched-string BVP with \(N=7\). The orange dots highlight where \(\Dfunc=0\). The function has been scaled so that \(\Dfunc(0) = 1\). (Source)

To understand this behavior, recall that the determinant of an \(N \times N\) matrix can be expressed (via Laplace expansion) as the sum of N terms; and each term itself involves the product of \(N\) matrix elements, picked so that each row/column is used only once in the construction of the term. With these points in mind, we can see from the definition (4) of \(\mS\) that its determinant (i.e., the discriminant function) must be a polynomial in \(\sigma^{2}\) of order \(N-2\); and as such, it can have at most \(N-2\) (in this case, 5) roots. This leads us to important lesson #1:

Attention

The number of points adopted in the discretization limits the number of modes that can be found. With a spatial grid of \(N\) points, there are only \(\sim N\) distinct numerical solutions.

Plot of the eigenfunctions \(\tilde{y}\) as a function of spatial coordinate \(x\), for the first three modes of the stretched-string BVP with \(N=7\). The discrete points show the numerical functions, and the solid lines the corresponding analytic functions. (Source)

Returning to Fig. 4, the shift in eigenfrequency for the modes that are found occurs due to inadequate resolution of the eigenfunctions. We can see this in Fig. 5, which reprises Fig. 3 for \(N=7\). Clearly, the spatial oscillations of the modes are poorly resolved; the \(n=3\) mode, for instance, is sampled with only one point per quarter wavelength. It’s little wonder that the corresponding eigenfrequencies are off. This brings us to important lesson #2 (closely related to #1):

Attention

The spatial resolution adopted in the discretization determines the accuracy of the modes found. A given eigenfrequency will be accurate only when the grid spacing is appreciably smaller than the spatial variation scale of its corresponding eigenfunction.

Insufficient Frequency Resolution

When searching for root brackets, we have to evaluate the discriminant function a total of \(M\) times. Therefore, as with the spatial grid, computational efficiency dictates that we want to make \(M\) as small as possible. Again, however, things go wrong if \(M\) is too small. Fig. 6 reprises Fig. 2, but adopting a much coarser frequency grid with only \(M=5\) points.

Plot of the discriminant values \(\{\Dfunc\}\) on the discrete frequency grid \(\{\sigma\}\) (distributed uniformly in \(\sigma\)), for the stretched-string BVP with \(N=50\) and \(M=5\). The orange-haloed segments highlight adjacent points that bracket a root \(\Dfunc=0\). (Source)

Clearly, a pair of adjacent roots (corresponding to the \(n=3\) and \(n=4\) modes) is missed in the bracketing process, as a direct result of the too-coarse grid.

Even when many points are included in the frequency grid, issues can still arise when the distribution of points doesn’t match the distribution of roots. An example of this is provided in Fig. 7, which reprises Fig. 2 with the same number \(M=32\) of points in the grid, but now distributed uniformly in \(\sigma^{-1}\).

Plot of the discriminant values \(\{\Dfunc\}\) on the discrete frequency grid \(\{\sigma\}\) (distributed uniformly in \(\sigma^{-1}\)), for the stretched-string BVP with \(N=50\) and \(M=32\). The orange-haloed segments highlight adjacent points that bracket a root \(\Dfunc=0\). (Source)

Now it’s the roots corresponding to the \(n=4\) and \(n=5\) mode pair that are missed. As with the case in Fig. 6, the failure ultimately arises because the spacing between adjacent frequency grid points is (in at least some parts of the grid) larger than the spacing between adjacent roots. This can be summarized in important lesson #3:

Attention

The frequency resolution adopted in the root bracketing influences the completeness of the modes found. All modes will be found only when the grid spacing is smaller than the eigenfrequency separation of adjacent modes, across the full range of the grid.

Interpreting Output Files

This chapter demonstrates using Python to read and plot the summary and detail output files written by the GYRE frontends. Further information about these files is provided in the Output Files chapter.

PyGYRE

PyGYRE is a Python package maintained separately from GYRE, that provides a set of routines that greatly simplify the analysis of summary and detail files. Detailed information about PyGYRE can be found in the full documentation; here, we demonstrate how to use it to read and plot the output files from the Example Walkthrough section.

As a preliminary step, you’ll need to install PyGYRE from the Python Package Index (PyPI). This can be done using the pip command,

pip install pygyre

If PyGYRE is already installed, you can upgrade to a more-recent version via

pip install --upgrade pygyre

Analyzing a Summary File

To analyze the summary file written by gyre during the example walkthrough, change into your work directory and fire up your preferred interactive Python environment (e.g., Jupyter). Import PyGYRE and the other modules needed for plotting:

# Import modules

import pygyre as pg
import matplotlib.pyplot as plt
import numpy as np

(you may want to directly cut and paste this code). Next, read the summary file into the variable s:

# Read data from a gyre summary file

s = pg.read_output('summary.h5')

The pygyre.read_output() function is able to read files in both TXT and HDF formats, returning the data in an astropy.table.Table object. To inspect the data on-screen, simply evaluate the table:

# Inspect the data

s

From this, you’ll see that there are three columns in the table, containing the harmonic degree l, radial order n_pg and frequency freq of each mode found during the GYRE run.

Next, plot the frequencies against radial orders via

# Plot the data

plt.figure()

plt.plot(s['n_pg'], s['freq'].real)

plt.xlabel('n_pg')
plt.ylabel('Frequency (cyc/day)')

(the values in the freq column are complex, and we plot the real part). The plot should look something like Fig. 8.

The frequency \(\nu\) of \(\ell=1\) and \(\ell=2\) modes, plotted against their radial order \(\numpg\). (Source)

The straight line connecting the two curves occurs because we are plotting both the dipole and quadrupole modes together. To separate them, the table rows can be grouped by harmonic degree:

# Plot the data, grouped by harmonic degree

plt.figure()

sg = s.group_by('l')

plt.plot(sg.groups[0]['n_pg'], sg.groups[0]['freq'].real, label=r'l=1')
plt.plot(sg.groups[1]['n_pg'], sg.groups[1]['freq'].real, label=r'l=2')

plt.xlabel('n_pg')
plt.ylabel('Frequency (cyc/day)')

plt.legend()

The resulting plot, in Fig. 9, looks much better.

The frequency nu of \(\ell=1\) and \(\ell=2\) modes, grouped by \(\ell\) and plotted against their radial order \(\numpg\). (Source)

Analyzing a Detail File

Now let’s take a look at one of the detail files, for the mode with \(\ell=1\) and \(\numpg=-7\). As with the summary file, pygyre.read_output() can be used to read the file data into an astropy.table.Table object:

# Read data from a GYRE detail file

d = pg.read_output('detail.l1.n-7.h5')

Inspecting the data using

# Inspect the data

d

shows there are 7 columns: the fractional radius x, the radial displacement eigenfunction xi_r, the horizontal displacement eigenfunction xi_h, and 4 further columns storing structure coefficients (see the Detail Files section for descriptions of these data). Plot the two eigenfunctions using the code

# Plot displacement eigenfunctions

plt.figure()

plt.plot(d['x'], d['xi_r'].real, label='xi_r')
plt.plot(d['x'], d['xi_h'].real, label='xi_h')

plt.xlabel('x')

plt.legend()

The radial (\(\txir\)) and horizontal (\(\txih\)) displacement eigenfunctions of the \(\ell=1\), \(n_{\rm pg}=-7\) mode, plotted against the fractional radius \(x\). (Source)

The plot should look something like Fig. 10. From this figure , we see that the radial wavelengths of the eigenfunctions become very short around a fractional radius \(x \approx 0.125\). To figure out why this is, we can take a look at the star’s propagation diagram:

# Evaluate dimensionless characteristic frequencies

l = d.meta['l']
omega = d.meta['omega']

x = d['x']
V = d['V_2']*d['x']**2
As = d['As']
c_1 = d['c_1']
Gamma_1 = d['Gamma_1']

d['N2'] = d['As']/d['c_1']
d['Sl2'] = l*(l+1)*Gamma_1/(V*c_1)

# Plot the propagation diagram

plt.figure()

plt.plot(d['x'], d['N2'], label='N^2')
plt.plot(d['x'], d['Sl2'], label='S_l^2')

plt.axhline(omega.real**2, dashes=(4,2))

plt.xlabel('x')
plt.ylabel('omega^2')

plt.ylim(5e-2, 5e2)
plt.yscale('log')

Note how we access the mode harmonic degree l and dimensionless eigenfrequency omega through the table metadata dict d.meta. The resulting plot (Fig. 11) reveals that the Brunt-Väisälä frequency squared is large around \(x \approx 0.125\); this feature is a consequence of the molecular weight gradient zone outside the star’s convective core, and results in the short radial wavelengths seen there in Fig. 10.

Propagation diagram for the \(5\,\Msun\) model, plotting the squares of the Brunt-Väisälä (\(N^{2}\)) and Lamb (\(S_{\ell}^{2}\)) frequencies versus fractional radius \(x\). The horizontal dashed line shows the frequency squared \(\omega^{2}\) of the \(\ell=1\), \(n_{\rm pg}=-7\) mode shown in Fig. 10. Regions where \(\omega^{2}\) is smaller (greater) than both \(N^{2}\) and \(S_{\ell}^{2}\) are gravity (acoustic) propagation regions; other regions are evanescent. (Source)

Understanding Grids

This chapter describes how the GYRE frontends set up frequency and spatial grids, and discusses some strategies for ensuring these grids are optimal.

Spatial Grids

The various GYRE frontends all discretize their equations on a spatial grid \(\{x_{1},x_{2},\ldots,x_{N}\}\) in the dimensionless radial coordinate \(x \equiv r/R\). The computational cost of a calculation scales with the total number of points \(N\) in this grid, while the grid’s resolution — i.e., the spacing between adjacent points — impacts both the accuracy of solutions, and in the case of the gyre frontend, the number of solutions that can be found. (The Limitations of the Numerical Method section discusses these behaviors in the context of the stretched string BVP).

Scaffold Grid

A fresh spatial grid is constructed for each iteration of the main computation loop (see the flow-charts in the Frontends chapter). This is done under the control of the &grid namelist groups; there must be at least one of these, subject to the tag matching rules (see the Working With Tags chapter). If there is more than one matching &grid namelist group, then the final one is used.

Each grid begins as a scaffold grid, comprising the following points:

• an inner point \(\xin\);

• an outer point \(\xout\);

• the subset of points of the source grid satisfying \(\xin < x < \xout\)

The source grid can be either the input model grid, or a grid read from file; this choice is determined by the scaffold_src parameter of the &grid namelist group. By default, \(\xin\) and \(\xout\) are obtained from the source grid as well (as its inner-most and outer-most point). However, either or both can be overridden using the x_i and x_o parameters.

Iterative Refinement

Scaffold grids are refined via a sequence of iterations. During a given iteration, each subinterval \([x_{j},x_{j+1}]\) is assessed against various criteria (discussed in greater detail below). If any criteria match, then the subinterval is refined by bisection, inserting an additional point at the midpoint

\[x_{j+\half} = \frac{x_{j} + x_{j+1}}{2}.\]

The sequence terminates if no refinements occur during a given iteration, or if the number of completed iterations equals the value specified by the n_iter_max parameter of the &grid namelist group.

Mechanical Criterion

The wave criterion involves a local analysis of the mechanical parts of the oscillation equations, with the goal of improving resolution where the displacement perturbation \(\vxi\) is rapidly varying. Within the subinterval \([x_{j},x_{j+1}]\), the \(y_{1}\) and \(y_{2}\) solutions (see the Dimensionless Formulation section) take the approximate form

\[y_{1,2}(x) \sim \exp [ \chi \, (\ln x - \ln x_{j+\half}) ],\]

where \(\chi\) is one of the two eigenvalues of the mechanical (upper-left) \(2 \times 2\) submatrix of the full Jacobian matrix \(\mA\), evaluated at the midpoint \(x_{j+\half}\).

In propagation zones the imaginary part \(\chi_{\rm i}\) of the eigenvalue gives the local wavenumber in \(\ln x\) space, and \(2\pi \chi_{\rm i}^{-1}\) the corresponding wavelength; while in evanescent zones the real part \(\chi_{\rm r}\) gives the local exponential growth/decay rate, and \(\chi_{\rm r}^{-1}\) the corresponding e-folding length.

Based on this analysis, the criterion for refinement of the subinterval is

\[( \ln x_{j+1} - \ln x_{j} ) \, \max (\wosc |\chi_{\rm i}|, \wexp |\chi_{\rm r}|) > 2 \pi,\]

where \(\wosc\) and \(\wexp\) are user-definable weighting parameters. This causes refinement if the subinterval width (in \(\ln x\) space) exceeds \(\wosc^{-1}\) times the local wavelength, or \(2\pi \wexp^{-1}\) times the local e-folding length.

Because there are two possible values for \(\chi\), the above refinement criterion is applied twice (once for each). Moreover, because \(\chi\) depends implicitly on the oscillation frequency, the criterion is applied for each frequency in the grid \(\{\omega_{1},\omega_{2},\ldots,\omega_{M}\}\) (see the Frequency Grids section).

Thermal Criterion

Similar to the wave criterion discussed above, the thermal criterion involves a local analysis of the energetic parts of the oscillation equation, with the goal of improving resolution where the thermal timescale is very long and perturbations are almost adiabatic. Within the subinterval \([x_{j},x_{j+1}]\), the \(y_{5}\) and \(y_{6}\) perturbation take the approximate form

\[y_{5,6}(x) \sim \exp [ \pm \tau \, (\ln x - \ln x_{j+\half}) ],\]

where \(\pm\tau\) are the eigenvalues of the matrix formed from the energetic (bottom-right) \(2 \times 2\) submatrix of the full Jacobian matrix \(\mA\), evaluated at the midpoint \(x_{j+\half}\).

Based on this analysis, the criterion for refinement of the subinterval is

\[( \ln x_{j+1} - \ln x_{j} ) \, \wthm |\tau| > 1,\]

where \(\wthm\) is a user-definable weighting parameter.

Because \(\tau\) depends implicitly on the oscillation frequency, this criterion is applied for each frequency in the grid \(\{\omega_{1},\omega_{2},\ldots,\omega_{M}\}\).

Structural Criteria

The structural criteria have the goal of improving resolution where the stellar structure coefficients are changing rapidly. For a given coefficient \(C\), the criterion for refinement of the subinterval \([x_{j},x_{j+1}]\) is

\[( \ln x_{j+1} - \ln x_{j} ) \, \wstr \left| \pderiv{\ln C}{\ln x} \right| > 1,\]

where \(\wstr\) is a user-definable weighting parameter. This criterion is applied separately to the \(V_2 \equiv V/x^{2}\), \(U\), \(A^{*}\), \(c_{1}\) and \(\Gamma_{1}\) coefficients (see the Structure Coefficients section).

Central Criteria

All of the above criteria depend on the logarithmic subinterval width \((\ln x_{j+1} - \ln x_{j})\), and cannot be applied to the first subinterval \([x_{1},x_{2}]\) if it extends to the center of the star, \(x = 0\). In such cases, the resolve_ctr parameter of the &grid namelist group determines whether the subinterval is refined. If set to .FALSE., then no refinement occurs; while if set to .TRUE., then the refinement criteria are

\[\chi_{\rm i} > 0\]

or

\[w_{\rm ctr} | \chi_{\rm r} | > 1\]

where \(\chi\) is the eigenvalue from the local analysis (see the Mechanical Criterion section) corresponding to the solution that remains well-behaved at the origin, and \(w_{\rm ctr}\) is a user-definable weighting parameter. The first criterion causes refinement if the subinterval is in a propagation zone, and the second if the solution slope \(|\sderiv{y}{\ln x}| \sim |\chi_{\rm r}|\) exceeds \(w_{\rm ctr}^{-1}\).

Because \(\chi\) depends implicitly on the oscillation frequency, these criteria are applied for each frequency in the grid \(\{\omega_{1},\omega_{2},\ldots,\omega_{M}\}\).

Limiting Controls

A couple of additional controls affect the iterative refinement described above. Refinement of the \([x_{j},x_{j+1}]\) subinterval always occurs if

\[x_{j+1} - x_{j} > \Delta x_{\rm max},\]

and never occurs if

\[x_{j+1} - x_{j} < \Delta x_{\rm min},\]

where both \(\Delta x_{\rm max}\) and \(\Delta x_{\rm min}\) are user-definable.

Namelist Parameters

The full set of parameters supported by the &grid namelist group is listed in the Grid Parameters section. However, the table below summarizes the mapping between the user-definable controls appearing in the expressions above, and the corresponding namelist parameters.

Symbol	Parameter
\(\wosc\)	w_osc
\(\wexp\)	w_exp
\(\wthm\)	w_thm
\(\wstr\)	w_str
\(\wctr\)	w_ctr
\(\Delta x_{\rm max}\)	dx_max
\(\Delta x_{\rm min}\)	dx_min

Recommended Values

While w_exp, w_osc and w_ctr all default to zero, it is highly recommended to use non-zero values for these parameters, to ensure adequate resolution of solutions throughout the star. Reasonable starting choices are w_osc = 10, w_exp = 2 and w_ctr = 10.

Frequency Grids

The gyre frontend evaluates its discriminant function \(\Dfunc(\omega)\) on a grid \(\{\omega_{1},\omega_{2},\ldots,\omega_{M}\}\) in the dimensionless frequency, and scans for changes in the sign of \(\Dfunc(\omega)\) that are indicative of a bracketed root. The computational cost of a calculation scales with the total number of points \(M\) in this grid, while the grid’s resolution — i.e., the spacing between adjacent points — impacts the completeness of the modes found by gyre. (See the Limitations of the Numerical Method section for a discussion of these behaviors).

A fresh frequency grid is constructed for each iteration of the main computation loop (see the flow-chart in the gyre section). This is done under the control of the &scan namelist groups; there must be at least one of these, subject to the tag matching rules (see the Working With Tags chapter). Each &scan group creates a separate frequency grid; these are then combined and the discriminant function is evaluated on the merged grid.

Grid Types

The grid_type parameter of the &scan namelist group controls the overall distribution of points in a frequency grid. There are currently three options:

Linear Grid

When grid_type = 'LINEAR', gyre first evaluates a sequence of dimensionless angular frequencies in the grid reference frame according to the formula

\[\omega^{\rm g}_{i} = \frac{1}{M-1} \left[ (M - i)\, \omega^{\rm g}_{\rm min} + (i - 1) \, \omega^{\rm g}_{\rm max} \right] \qquad i = 1,2,\ldots,M.\]

(here, the superscript ‘g’ indicates that these are frequencies in the grid reference frame). Then, it transforms from the grid frame to the inertial reference frame via

\[\omega_{i} = \omega^{\rm g}_{i} + \Delta \omega\]

where \(\Delta\omega\) is the frequency shift that transforms from the grid frame to the inertial frame. The actual value of this shift depends on the grid_frame parameter:

• When grid_frame = 'INERTIAL', the shift is \(\Delta \omega = 0\); the grid frame and the inertial frame coincide.

• When grid_frame = 'COROT_I', the shift is \(\Delta \omega = m \, \Orot^{\rm i} \sqrt{R^{3}/GM}\), where \(\Orot^{\rm i}\) is the rotation angular frequency at the inner boundary of the spatial grid; the grid frame coincides with the local co-rotating frame at that boundary.

• When grid_frame = 'COROT_O', the shift is \(\Delta \omega = m \, \Orot^{\rm o} \sqrt{R^{3}/GM}\), where \(\Orot^{\rm o}\) is the rotation angular frequency at the outer boundary of the spatial grid; the grid frame coincides with the local co-rotating frame at that boundary.

The range spanned by the frequency grid, in the grid frame, is set by \(\omega^{\rm g}_{\rm min}\) and \(\omega^{\rm g}_{\rm max}\). These are evaluated via

\[\omega^{\rm g}_{\rm min} = \frac{f_{\rm min}}{\widehat{f}_{\rm min}} + \delta \omega - \Delta \omega, \qquad \qquad \omega^{\rm g}_{\rm max} = \frac{f_{\rm max}}{\widehat{f}_{\rm max}} + \delta \omega - \Delta \omega,\]

where \(f_{\rm min,max}\) are user-definable, \(\widehat{f}_{\rm min,max}\) will be discussed below in the Frequency Units section, and \(\delta\omega\) is the frequency shift that transforms from the frame in which \(f_{\rm min,max}\) are defined to the inertial frame. The actual value of this shift depends on the freq_frame parameter, which behaves analogously to the grid_frame parameter discussed above.

Inverse Grid

When grid_type = 'INVERSE', gyre first evaluates a sequence of dimensionless angular frequencies in the grid reference frame according to the formula

\[\omega^{\rm g}_{i} = (M-1) \left[ \frac{(M - i)}{\omega^{\rm g}_{\rm min}} + \frac{(i - 1)}{\omega^{\rm g}_{\rm max}} \right]^{-1} \qquad i = 1,2,\ldots,M.\]

The grid creation then proceeds as described above in the Linear Grid section.

File Grid

When grid_type = 'FILE', gyre first reads a sequence of dimensioned frequencies \(\{f_{1},f_{2},\ldots,f_{M}\}\) from an external file named by the grid_file parameter. This file is a single-column ASCII table; the number of points \(M\) is determined implicitly from the number of lines in the file. Then, it transforms these frequencies via

\[\omega_{i} = \frac{f_{j}}{\widehat{f}} + \delta \omega,\]

where \(\widehat{f}\) will be discussed below in the Frequency Units section, and \(\delta\omega\) is the frequency shift that transforms from the frame in which \(f\) is defined to the inertial frame. The actual value of this shift depends on the freq_frame parameter, which behaves analogously to the grid_frame parameter discussed above.

Frequency Units

In the expressions above, terms of the form \(f/\widehat{f}\) are used to transform a dimensioned frequency \(f\) into a dimensionless one \(\omega\). The scale factor \(\widehat{f}\) depends on the freq_units parameter. Thus, for example, if freq_units = 'UHZ', then \(f\) is treated as a linear frequency expressed in \({\rm \mu Hz}\), and the scale factor is set by

\[\widehat{f} = \sqrt{\frac{GM}{R^{3}}} \frac{1}{2\pi\,{\rm \mu Hz}}\]

(the factor of \(2\pi\) comes from the transformation between linear and angular frequency).

The full set of values supported by the freq_units parameter is listed in the Frequency Scan Parameters section.

Namelist Parameters

The full set of parameters supported by the &scan namelist group is listed in the Frequency Scan Parameters section. However, the table below summarizes the mapping between the user-definable controls appearing in the expressions above, and the corresponding namelist parameters:

Symbol	Parameter
\(f_{\rm min}\)	freq_min
\(f_{\rm max}\)	freq_max
\(M\)	n_freq

Recommended Values

The default values freq_min=1, freq_max=10, n_freq=10, together with grid_type='LINEAR' are sufficient to find some modes — although unlikely the modes that you want. Choosing good values for these parameters requires some degree of judgment, but here are some suggestions:

• The number of points in the frequency grid should be a factor of 2–3 larger than the number of modes you expect gyre will find. This is to ensure that the frequency spacing of the grid is everywhere smaller than the anticipated eigenfrequency spacing between adjacent modes (see the Limitations of the Numerical Method section for further discussion).

• The distribution of points in the frequency grid should follow anticipated distribution of mode frequencies; this again is to ensure adequate frequency resolution. For p modes, which tend toward a uniform frequency spacing in the asymptotic limit of large radial order, you should chose grid_type = 'LINEAR'; likewise, for g modes, which tend toward a uniform period spacing in the asymptotic limit, you should choose grid_type = 'INVERSE'.

• When modeling rotating stars, you should choose grid_frame = 'COROT_I' or grid_frame = 'COROT_O', because the asymptotic behaviors mentioned above apply in the co-rotating reference frame rather than the inertial one.

Working With Tags

This chapter uses a worked example to demonstrate tags — a simple yet powerful feature allowing a much greater degree of control over GYRE calculations.

Example Tag Usage

Consider applying gyre to calculate the eigenfrequencies of a red giant branch (RGB) stellar model. Because non-radial p-modes in the convective envelope couple with high-order g-modes in the radiative core, the frequency spacing of the non-radial modes is much smaller than that of the radial modes. In such cases, we ideally want to use a coarse frequency scan for the radial modes and a fine frequency scan for the non-radial modes.

The following input file, which is designed to work with the \(2\,\Msun\) RGB model in $GYRE_DIR/models/mesa/rgb/rgb.mesa, achieves this goal using tags:

&constants
/

&model
  model_type = 'EVOL'  ! Obtain stellar structure from an evolutionary model
  file = 'rgb.mesa'    ! File name of the evolutionary model
  file_format = 'MESA' ! File format of the evolutionary model
/

&mode
  l = 0          ! Harmonic degree
  tag = 'radial' ! Matching tag
/

&mode
  l = 1              ! Harmonic degree
  tag = 'non-radial' ! Matching tag
/

&mode
  l = 2              ! Harmonic degree
  tag = 'non-radial' ! Matching tag
/

&osc
  outer_bound = 'VACUUM' ! Assume the density vanishes at the stellar surface
/

&rot
/

&num
  diff_scheme = 'COLLOC_GL4' ! 4th-order collocation scheme for difference equations
/

&scan
  grid_type = 'LINEAR' ! Scan grid uniform in frequency
  freq_min = 41        ! Minimum frequency to scan from
  freq_max = 43        ! Maximum frequency to scan to
  freq_units = 'UHZ'   ! Frequency units
  n_freq = 10          ! Number of frequency points in scan
  tag_list = 'radial'  ! Pair only with 'radial' &mode groups
/

&scan
  grid_type = 'LINEAR' 	  ! Scan grid uniform in frequency
  freq_min = 41        	  ! Minimum frequency to scan from
  freq_max = 43           ! Maximum frequency to scan to
  freq_units = 'UHZ'      ! Frequency units
  n_freq = 100            ! Number of frequency points in scan
  tag_list = 'non-radial' ! Pair only with 'non-radial' &mode groups
/

&grid
  w_osc = 10  ! Oscillatory region weight parameter
  w_exp = 2   ! Exponential region weight parameter
  w_ctr = 100 ! Central region weight parameter
/

&ad_output
  summary_file = 'summary.h5'                  ! File name for summary file
  summary_item_list = 'l,n_pg,freq,freq_units' ! Items to appear in summary file
  freq_units = 'UHZ'                           ! Units of freq output items
/

&nad_output
/

Observe that each &mode namelist groups has a tag parameter. When processing a given &mode, gyre pairs it up with other namelist groups that match one of the following criteria:

• The namelist group doesn’t have a tag_list parameter;

• The namelist does have a tag_list parameter, and the parameter value (a comma-separated list) contains the tag value defined in the &mode group.

In the example given above, the &osc namelist group doesn’t have a tag_list parameter; therefore, it is paired with all three &mode namelist groups, irrespective of their tag values. However, the two &scan namelist groups each have tag_list parameters. In the first group the radial tag appears, and so this group is paired with the first &mode namelist group (i.e., the \(\ell=0\) mode). Likewise, in the second group the non-radial tag appears, and so this group is paired with the second and third &mode namelist groups (i.e., the \(\ell=1\) and \(\ell=2\) modes).

Tag Rules

In addition to the matching criteria given above, there are a couple of rules that must be obeyed by tags:

• Tag names can’t contain commas (however, they can be otherwise arbitrary);

• If a &mode namelist group doesn’t have a tag parameter, then only namelists without a tag_list parameter will be paired with it;

• The &constants, &model, &ad_output and &nad_output namelist groups don’t support tags.

Advanced Usage

This chapter covers a number of topics that fall under the umbrella of ‘advanced’ GYRE usage.

Non-Adiabatic Oscillations

This section discusses how to undertake non-adiabatic oscillation calculations using the gyre frontend. Asteroseismic studies typically rely on adiabatic calculations, because the frequencies of oscillation modes are the primary focus. However, for heat-driven modes the linear growth or damping rates can also be of interest — and evaluating these requires that non-adiabatic effects are included in the oscillation equations.

Note

Not all types of stellar mode include the necessary data (e.g., thermodynamic coefficients, opacity partial derivatives) to undertake non-adiabatic calculations. The Model Capabilities section summarizes this information.

Overview

To include non-adiabatic effects gyre augments the linearized mass, momentum and Poisson equations with the linearized heat and radiative diffusion equations (see the Linearized Equations section for full details). With these additions, the equations and their solutions become complex quantities. The assumed time dependence for perturbations is \(\propto \exp (-\ii \sigma t)\); therefore, the real part \(\sigmar\) and imaginary part \(\sigmai\) of the eigenfrequency are related to the mode period \(\Pi\) and growth e-folding time \(\tau\), respectively, via

\[\Pi = \frac{2\pi}{\sigmar}, \qquad \tau = \frac{1}{\sigmai}.\]

Solving the non-adiabatic equations proceeds using the same general approach as in the adiabatic case, by searching for the roots of a discriminant function \(\Dfunc(\omega)\) (see the Numerical Methods chapter for more details). However, a challenge is that there is no simple way to bracket roots in the complex plane. Instead, gyre must generate initial trial roots that are close to the true roots, and then refine them iteratively. Currently, gyre offers three methods for establishing the trial roots.

Adiabatic Method

The adiabatic method involves adopting the (real) roots found from adiabatic calculations as the initial trial roots for the non-adiabatic problem. This works well as long as the adiabatic and non-adiabatic roots lie close together in the complex plane — typically, when the oscillation modes are only weakly non-adiabatic, with \(|\sigmai/\sigmar| \ll 1\).

To perform non-adiabatic calculations with the adiabatic method, set the following parameters in the &osc namelist group:

• nonadiabatic=.TRUE.

• adiabatic=.TRUE.[1]

and the following parameters in the &num namelist group:

• ad_search='BRACKET'[1]

• nad_search='AD'[1]

You may also wish to use the following setting in the &num namelist group:

• diff_scheme='MAGNUS_GL2'

This tells gyre to evaluate the finite-difference equations using the 2nd order Magnus scheme; experience suggests that this gives the most reliable convergence for the root refinement.

An example of the adiabatic method in action can be found in the $GYRE_DIR/test/nad/mesa/bcep/gyre.in namelist input file, which is set up to find \(\ell=0,\ldots,3\) modes of a \(20\,\Msun\) \(\beta\) Cephei model using the adiabatic method. The important parts are as follows:

&osc
  nonadiabatic = .TRUE.
/

&num
  diff_scheme = 'MAGNUS_GL2'
  restrict_roots = .FALSE.
/

&scan
  grid_type = 'LINEAR'
  freq_min = 3.0
  freq_max = 10.0
  n_freq = 50
/

Note the nonadiabatic parameter in the &osc namelist group, and the diff_scheme parameter in the &num namelist group. The restrict_roots=.FALSE. setting, also in the &num namelist group, tells gyre not to reject any modes that have \(\sigmar\) outside the frequency range specified by the &scan namelist group; this ensures that modes whose non-adiabatic frequencies fall just outside the frequency grid are still found.

Minmod Method

The minmod method involves evaluating the discriminant function along the real-\(\omega\) axis, and then adopting local minima in its modulus \(|\Dfunc|\) as the initial trial roots for the non-adiabatic problem. The method is described in full in Goldstein & Townsend (2020); as shown there, it does not perform significantly better than the adiabatic method, and is included in gyre for the sake of completeness.

To perform non-adiabatic calculations with the minmod method, set the following parameters in the &osc namelist group:

• nonadiabatic=.TRUE.

• adiabatic=.FALSE.[2]

and the following parameters in the &num namelist group:

• nad_search='MINMOD'

As with the adiabatic method, you may also wish to use the following setting in the &num namelist group:

• diff_scheme='MAGNUS_GL2'

An example of the minmod method in action can be found in the $GYRE_DIR/test/nad/mesa/bcep-minmod/gyre.in namelist input file, which is equivalent to $GYRE_DIR/test/nad/mesa/bcep/gyre.in but using the minmod method. The important parts are as follows:

&osc
  adiabatic = .FALSE.
  nonadiabatic = .TRUE.
/

&num
  diff_scheme = 'MAGNUS_GL2'
  nad_search = 'MINMOD'
  restrict_roots = .FALSE.
/

&scan
  grid_type = 'LINEAR'
  freq_min = 3.0
  freq_max = 10.0
  n_freq = 250
/

Note the additional nad_search='MINMOD' parameter in the &num namelist group, which stipulates that the minmod method should be used.

Contour Method

The contour method involves evaluating the discriminant function on a grid in the complex-\(\omega\) plane, and then adopting intersections between the real zero-contours \(\Dfuncr=0\), and the corresponding imaginary ones \(\Dfunci=0\), as the initial trial roots for the non-adiabatic problem. The method is described in full in Goldstein & Townsend (2020); it is very effective even for strongly non-adiabatic modes with \(|\sigmai/\sigmar| \sim 1\), although there is an increased computational cost (see here for one strategy for mitigating this cost).

To perform non-adiabatic calculations with the contour method, set the following parameters in the &osc namelist group:

• nonadiabatic=.TRUE.

• adiabatic=.FALSE.[2]

and the following parameters in the &num namelist group:

• nad_search='CONTOUR'

You must also ensure that at least one &scan namelist group with axis='REAL' is present, and likewise at least one with axis='IMAG'. Together, these groups define the real and imaginary axes of the discriminant grid in the complex-\(\omega\) plane. As a rule of thumb, the resolution along the imaginary axis should be comparable to that along the real axis; this ensures that the contour-tracing algorithm behaves well.

Finally, as with the adiabatic method, you may also wish to use the following setting in the &num namelist group:

• diff_scheme='MAGNUS_GL2'

Note

Because g modes are spaced uniformly in period (in the asymptotic limit of large radial order), it would seem sensible to set grid_type='INVERSE' in the &scan namelist group(s) that correspond to the real axis (i.e., axis='REAL'). However, this typically results in a mismatch between the resolution of the real and imaginary axes, and the contour method doesn’t perform well. A fix for this issue will be forthcoming in a future release of GYRE, but in the meantime it’s probably best to avoid the contour method for g modes.

An example of the minmod method in action can be found in the $GYRE_DIR/test/nad/mesa/bcep-contour/gyre.in namelist input file, which is equivalent to $GYRE_DIR/test/nad/mesa/bcep/gyre.in but using the minmod method. The important parts are as follows:

&osc
  adiabatic = .FALSE.
  nonadiabatic = .TRUE.
/

&num
  diff_scheme = 'MAGNUS_GL2'
  restrict_roots = .FALSE.
  nad_search = 'CONTOUR'
/

&scan
  axis = 'REAL'
  grid_type = 'LINEAR'
  freq_min = 3.0
  freq_max = 10.0
  n_freq = 50
/

&scan
  axis = 'IMAG'
  grid_type = 'LINEAR'
  freq_min = -0.28
  freq_max = 0.28
  n_freq = 5
/

Note the additional nad_search='CONTOUR' parameter in the &num namelist group, which stipulates that the contour method should be used; and, the fact that there are now two &scan namelist groups, one with axis='REAL' and the other with axis='IMAG'.

Footnotes

[1] (1,2,3)

This is the default setting; you don’t need to include it explicitly.

[2] (1,2)

This is optional; leave it out if you want gyre to perform adiabatic calculations as well.

Tidal Forcing

This section discusses how to evaluate the stellar response (fluid displacements and perturbations) to tidal forcing, using the gyre_tides frontend. The response data can be used to calculate the secular rates-of-change of orbital elements, or to synthesize a light curve for a tidally distorted star.

Overview

As discussed in the Tidal Effects section, the tidal gravitational potential (16) of an orbiting companion can be expressed as a superposition of partial potentials \(\PhiTlmk\). For a given &tide namelist group appearing in the namelist input file, gyre_tides solves for the response of the star to each term in the superposition.

Truncating the Sums

Although the sums appearing in eqn. (16) are formally infinite, the terms with large harmonic degree \(\ell\) and/or orbital harmonic \(k\) typically produce a negligible response. gyre_tides offers a couple of approaches for truncating the sums by dropping these terms. The simplest is to set limits on the maximum values of the indices, through the l_max, k_min and k_max parameters of the &tide namelist group (if desired, minimum values can also be set using the corresponding l_min and k_min parameters).

A slightly more sophisticated approach is to set these parameters to large-ish values (say, a few hundred), and then also set one or both of the y_T_thresh_abs and y_T_thresh_rel parameters. These establish a threshold on the magnitude of

\[\yT \equiv \frac{\tPhiTlmk}{GM/R}\]

for a given tidal partial potential \(\tPhiTlmk\) (see eqn. 18) to be included in calculations; if \(|\yT|\) does not meet this threshold, it is ignored.

Optimizing Grids

During the iterative refinement process used in setting up spatial grids, the refinement criteria are evaluated for every partial tide under consideration. If the co-rotating forcing frequency

\[\sigmac \equiv k \Oorb - m \Orot\]

associated with a specific partial tidal potential is small compared to the dynamical frequency of the star, many levels of refinement will occur. While this is exactly what one wants in oscillation calculations (because low-frequency modes have short spatial wavelengths), it often isn’t necessary in tidal calculations because the response of a star to low-frequency forcing is the long-wavelength equilibrium tide.

One way of preventing over-refinement due to low-frequency partial potentials is to set the omega_c_thresh parameter in the &tide namelist group. This establishes a threshold on the dimensionless frequency \(\omegac \equiv \sigmac \, \sqrt{R^{3}/GM}\); partial potentials with \(|\omegac|\) below this threshold are treated as static tides (\(\omegac=0\)), and are not considered during the iterative refinement process.

An alternative approach is to avoid iterative refinement altogether, instead obtaining the spatial grid from an external file (see the 'FILE' option of the scaffold_src parameter, in the Grid Parameters section). This is the most flexible approach, but creating a grid that will adequately resolve the response to each partial potential requires some care.

Output Files

gyre_tides writes response data to summary and detail files. One detail file is created for each partial potential evaluated, and the summary file collects together global data for all partial potentials across all &tide namelist groups. The id output item can be used to determine which group a given response belongs to.

The following Python code demonstrates how the summary data might be used to evaluate the secular rates-of-change of orbital semi-major axis, eccentricity, and argument of periastron, and the stellar torque. The expression for e_dot mirrors eqn. (23) of Sun et al. (2023), and for J_dot eqn. (25) ibid.

import pygyre as pg
import numpy as np

# Read summary file from gyre_tides

s = pg.read_output('summary.h5')

# Extract the first set of responses

sg = s.group_by('id').groups[0]

Omega_orb = sg['Omega_orb']               
R_a = sg['R_a']
q = sg['q']

eps_T = R_a**3*q
    
l = sg['l']
m = sg['m']
k = sg['k']

cbar = sg['cbar']
        
Fbar = -0.5*sg['eul_Phi_ref']/(cbar*eps_T)

x = sg['x_ref']

Gbar_1 = sg['Gbar_1']
Gbar_2 = sg['Gbar_2']
Gbar_3 = sg['Gbar_3']
Gbar_4 = sg['Gbar_4']

# Evaluate secular rates-of-change

kap = np.empty(len(l))

for i in range(len(kap)):
    if k[i] == 0:
        if m[i] == 0:
            kap[i] = 0.5
        elif m[i] > 0 and m[i] <= l[i]:
            kap[i] = 1.
        else:
            kap[i] = 0.
    elif k[i] > 0:
        kap[i] = 1.
    else:
        kap[i] = 0.

# Semi-major axis (units of R per dynamical timescale)

a_dot = np.sum(4. * Omega_orb * q / R_a * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_2)

# Eccentricity (units of per dynamical timescale)

e_dot = np.sum(4. * Omega_orb * q * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_3)

# Argument of periastron (units of radians per dynamical timescale)

pom_dot = np.sum(4. * Omega_orb * q * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.real * Gbar_1)

# Angular momentum (units of GM^2/R)

J_dot = np.sun(4. * q**2 * R_a *
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_4)

Including Rotation

This section discusses how to include the effects of rotation in GYRE calculations. See the Rotation Effects section for further details of GYRE’s rotation treatment.

Setting the Rotation Rate

There are two different ways to define the rotation angular frequency \(\Orot\), via parameters in the &rot namelist group.

• If Omega_rot_source = 'MODEL', then differential rotation is assumed with a spatially varying \(\Orot\) obtained from the stellar model. If the model doesn’t have this capability (see the Model Capabilities section), then \(\Orot\) is set to zero throughout the star.

• If Omega_rot_source = 'UNIFORM', then uniform rotation is assumed with a spatially constant \(\Orot\) set by the Omega_rot and Omega_rot_units parameters.

Incorporating Doppler Effects

The Doppler Shift effect is incorporated automatically whenever calculations are performed with non-zero \(\Orot\) and mode azimuthal order \(m\). To disable this effect, set m= 0 in the &mode namelist group.

Incorporating Coriolis Effects

Incorporating the effects of the Coriolis force can be done using a perturbative treatment or a non-perturbative treatment. In the former case the effects are be applied as a post-calculation correction to non-rotating eigenfrequencies (see the domega_rot output item in the Summary Files and Detail Files sections). In the latter case, the traditional approximation of rotation (TAR) can be enabled by setting coriolis_method='TAR' in the &rot namelist group.

The TAR solution family is controlled by the rossby parameter of the &rot namelist group; set to .TRUE. for the Rossby family, and to .FALSE. for the gravito-acoustic family.

Frequently Asked Questions

This is a list of Frequently Asked Questions about GYRE. Suggestions for new entries are always most welcome!

How Do I…

…obtain the GYRE source code?

See the Download section.

…compile GYRE?

See the Compile section.

…compile a statically linked version of GYRE?

Set the STATIC environment variable to yes prior to compilation. Note that this currently only works on Linux platforms, and when then CRMATH environment variable is set to no.

…run GYRE on multiple cores?

GYRE can take advantage of multiple processors in a shared-memory (multicore) computer through its use of OpenMP. This functionality should be enabled by default, but you can nevertheless force it by setting the OMP environment variable to yes prior to compilation. Then, set the OMP_NUM_THREADS environment variable to the number of threads you want to use.

…run GYRE on a cluster?

GYRE can take advantage of multiple nodes in a computer cluster through its use of MPI. To enable this functionality, set the MPI environment variable to yes prior to compilation. Note that you’ll need to use a version of the MESA Software Development Kit built with MPI support for your cluster’s specific networking technology (e.g., Infiniband); contact the GYRE team for more details.

…pronounce GYRE?

With a soft ‘g’ rhyming with ‘wire’, like this.

…cite GYRE?

See the Citing GYRE section.

…get support for a problem I’m having?

Post a message to one of the GYRE discussion forums.

…contact the GYRE team?

Send an email to Rich Townsend.

…access the documentation for older releases?

Click on the ‘v:…’ dropdown menu at the bottom of the left-hand panel. Note that this menu is only available when viewing the documentation on Read the Docs; it isn’t available in a local build within the $GYRE_DIR/doc/sphinx directory. Also, the legacy documentation is a work-in-progress, and remains incomplete.

What Does…

…’Failed during deflate narrowout-of-domain frequency’ mean?

This is an indication that GYRE’s root solver wandered out of bounds when trying to find a complex root of the discriminant function. Try running with a different choice of diff_scheme parameter (MAGNUS_GL2 seems to be the most robust), and/or using contour method instead (see the Non-Adiabatic Oscillations chapter).

Why Does…

…the error ‘Illegal Instruction’ arise on MacOS when running with large grid sizes?

This behavior is typically caused by overflow of the OpenMP stack (for more info see here). Try setting the OMP_STACKSIZE environment variable to 500K or 1M.

Troubleshooting

This chapter discusses various problems that can arise during normal GYRE operation, and steps that can be taken to resolve them.

Missing Modes

For adiabatic oscillation calculations using gyre, the radial order \(\numpg\) of modes found should be monotonic-increasing[1]. Departures from this behavior can occur for a number of reasons.

Insufficient Frequency Resolution

If the frequency grid has insufficient resolution, then gyre can skip modes during the bracketing phase, as discussed in the Limitations of the Numerical Method section. The signature of insufficient frequency resolution is that an even number of consecutive modes is missed — most often, an adjacent pair of modes.

To fix this problem, first check that the distribution of points in the frequency grids matches (approximately) the expected distribution of mode eigenfrequencies:

• In the asymptotic limit of large radial order, p modes are uniformly distributed in frequency (see, e.g., Aerts et al., 2010). Hence, to search for these modes set grid_type='LINEAR' in the &scan namelist group(s).

• Likewise, in the asymptotic limit of large radial order, g modes are uniformly distributed in period. Hence, to search for these modes set grid_type='INVERSE' in the &scan namelist group(s).

• For rotating stars, the asymptotic behaviors mentioned apply in the co-rotating reference frame, not in the inertial reference frame. So, be sure to also set grid_frame ='COROT_I'|'COROT_O' in the &scan namelist group.

Next, try increasing the number of points in the frequency grids, simply by increasing the n_freq parameter in the &scan namelist group(s).

Tip

A good rule of thumb is that n_freq should be around 5 times larger than the number of modes expected to be found.

Insufficient Spatial Resolution

If the spatial grid has insufficient resolution, then certain modes can simply be absent from the (finite) set of distinct numerical solutions, as discussed in the Limitations of the Numerical Method section. The signature of insufficient spatial resolution is that modes that are found have radial orders comparable to the number of grid points \(N\) in the grid; and that the eigenfunctions of these modes are barely resolved (cf. Fig. 5).

To fix this problem, first check that the w_osc, w_exp and w_ctr weighting parameters in the &grid namelist group are set to reasonable values (see the Recommended Values section). If that doesn’t improve things, try gradually increasing both w_osc and w_ctr.

Non-adiabatic Effects

When undertaking non-adiabatic calculations, modes can be mis-classified or completely missed. The former situation arises because the expectation of monotonic-increasing \(\numpg\) formally applies only to adiabatic oscillations; while it can also work reasonably well for weakly non-adiabatic cases, there are no guarantees. If mis-classification does occur, then it must be fixed manually by determining which adiabatic mode the problematic non-adiabatic mode corresponds to.

Missing modes occur for a different reason: if a mode has a large growth rate, then the usual adiabatic method for establishing initial trial roots can fail to find it. In such cases, the alternative contour method performs very well.

Footnotes

[1]

The sole exception is \(\ell=1\) modes, where \(\numpg=0\) is skipped due to the way the Takata 2006b classification scheme is set up.

Duplicated Modes

Sometimes two oscillation modes with the same \(\numpg\) are found during a gyre calculation. This violates the expectation that \(\numpg\) be monotonic-increasing, and can happen for a few reasons.

Bad Stellar Model

If the input stellar doesn’t conserve mass properly, then one or more bogus (unphysical) modes can appear with the same \(\numpg\) as an existing mode, but a significantly different frequency. Such modes can also be spotted because their radial order is very different from the adjacent-in-frequency modes.

Bogus modes arise because the input stellar model doesn’t conserve mass. GYRE` assumes that the density \(\rho\) and interior mass \(M_{r}\) are related by equation 5. Given that there are many different ways to discretize this equation, there is a certain amount of numerical ‘slop’ that arises when going from discrete form assumed in the stellar evolution code that generated the model (e.g., MESA) through to the discrete form implicitly assumed by GYRE. This slop can sometimes produce bogus modes, and the fix is to re-create the model with a higher spatial resolution.

Another possibility arises if the input stellar model contains density discontinuities that aren’t properly marked using double points. Often, these discontinuities can be diagnosed by plotting the \(I_0\) or \(I_1\) first integrals (for radial and dipole modes, respectively) as a function of radius. These first integrals should vanish everywhere (as shown by Takata 2006a), but will typically show abrupt jumps to non-zero values at the location of unmarked discontinuities. The fix is to re-create the model with double points inserted as necessary; in the case of MESA models, this can be achieved using the add_double_points_to_pulse_data parameter of the &controls namelist group.

Non-adiabatic Effects

When non-adiabatic effects cause a mode to be mis-classified (as discussed in the Missing Modes section), often the incorrect \(\numpg\) value duplicates that of another mode. As before, a mis-classification must be fixed manually by determining which adiabatic mode the problematic non-adiabatic mode corresponds to.

Long Runtimes

gyre

Long runtimes with the gyre frontend occur when the spatial grid and/or frequency grid contain many points. The execution time to process a single &mode namelist group can be approximated by

\[\tau \approx C_{\rm b} N M + C_{\rm s} N N_{\rm m},\]

where \(N\) is the number of spatial points, \(M\) is the number of frequency points, \(N_{\rm m}\) is the number of modes found, and \(C_{\rm b}\) and \(C_{\rm s}\) are constants. The first (\(C_{\rm b}\)) term represents the time take to bracket roots of the discriminant function, and the second (\(C_{\rm s}\)) the time taken to solve for these roots (see the Numerical Methods chapter for details).

The key to ensuring reasonable runtimes lies in judicious choice of parameters in the &scan namelist group(s). The n_freq parameter obviously has an impact on \(\tau\), as it directly sets \(M\). However, the freq_min and freq_max parameters also influence \(\tau\), due to the way the spatial grid is constructed. If the frequency scan includes parts of the star’s oscillation spectrum containing modes with very large radial orders (whether p modes or g modes), then GYRE’s iterative refinement algorithm will insert many grid points in order to resolve the modes’ wavefunctions. This can ultimately lead to huge \(N\) and very long runtimes.

gyre_tides

The narrative is similar with the gyre_tides frontend, although there are no frequency grids involved. The execution time to process a single &orbit namelist group can be approximated by

\[\tau \approx C_{\rm t} N\]

where \(C_{\rm t}\) is a constant.

Installation

This chapter discusses GYRE installation in detail. If you just want to get up and running, have a look at the Quick Start chapter.

Pre-Requisites

To compile and run GYRE, you’ll need the following software components:

• A modern (2003+) Fortran compiler

• The BLAS linear algebra library

• The LAPACK linear algebra library

• The LAPACK95 Fortran 95 interface to LAPACK

• The HDF5 data management library

• The crlibm correctly rounded math library

• The crmath Fortran 2003 interface to crlibm

• An OpenMP-aware version of the ODEPACK differential equation library (optional)

On Linux and MacOS platforms, these components are bundled together in the MESA Software Development Kit (SDK), which can be downloaded from the MESA SDK homepage. Using this SDK is strongly recommended.

Building GYRE

Download

Download the GYRE source code, and unpack it from the command line using the tar utility:

tar xf gyre-7.1.tar.gz

Set the GYRE_DIR environment variable with the path to the newly created source directory; this can be achieved, e.g., using the realpath command[1]:

export GYRE_DIR_DIR=$(realpath gyre-7.1)

Compile

Compile GYRE using the make utility:

make -j -C $GYRE_DIR

(the -j flags tells make to use multiple cores, speeding up the build).

Test

To check that GYRE has compiled correctly and gives reasonable results, you can run the calculation test suite via the command

make -C $GYRE_DIR test

The initial output from the tests should look something like this:

TEST numerics (OpenMP)...
 ...succeeded
TEST numerics (band matrix)...
 ...succeeded
TEST numerics (*_DELTA frequency units)...
 ...succeeded
TEST numerics (rotation, Doppler shift)...
 ...succeeded
TEST numerics (rotation, traditional approximation)...
 ...succeeded

If things go awry, consult the Troubleshooting chapter.

Custom Builds

Custom builds of GYRE can be created by setting certain environment variables, and/or variables in the file $GYRE_DIR/src/build/Makefile, to the value yes. The following variables are currently supported:

DEBUG

Enable debugging mode (default no)

OMP

Enable OpenMP parallelization (default yes)

MPI

Enable MPI parallelization (default no)

DOUBLE_PRECISION

Use double precision floating point arithmetic (default yes)

CRMATH

Use correctly rounded math functions (default yes)

IEEE

Use Fortran IEEE floating point features (default no)

FPE

Enable floating point exception checks (default yes)

HDF5

Include HDF5 support (default yes)

EXPERIMENTAL

Enable experimental features (default no)

If a variable is not set, then its default value is assumed.

Git Access

Sometimes, you’ll want to try out new features in GYRE that haven’t yet made it into a formal release. In such cases, you can check out GYRE directly from the https://github.com/rhdtownsend/gyre git repository on GitHub:

git clone https://github.com/rhdtownsend/gyre.git

However, a word of caution: GYRE is under constant development, and features in the main (master) branch can change without warning.

footnote

[1]

The realpath command is included in the GNU CoreUtils package. Mac OS users can install CoreUtils using MacPorts or Homebrew.

Namelist Input Files

This chapter describes the various groups that can appear in the input files read by the GYRE frontends. These files are in Fortran’s namelist format, a simple text-based format containing one or more namelist groups. Each group begins with the line &name (where name is the name of the group); a list of parameter-value pairs then follows, and the group ends with a slash /.

Constants

The &constants namelist group defines various physical constants, as follows:

G_GRAVITY

Gravitational constant \(G\)

C_LIGHT

Speed of light in vacuo \(c\)

A_RADIATION

Radiation constant \(a\)

M_SUN

Solar mass \(\Msun\)

R_SUN

Solar radius \(\Rsun\)

L_SUN

Solar luminosity \(\Lsun\)

GYRE_DIR

Top-level GYRE directory; overrides the GYRE_DIR environment variable

All of these constants are in cgs units (where applicable), and the default values are defined in $GYRE_DIR/src/common/gyre_constants.fpp.

Grid Parameters

The &grid namelist group defines the parameters used to set up the spatial grid, as follows:

scaffold_src (default 'MODEL')

Source for scaffold grid; one of:

• 'MODEL' : Obtained from the stellar model

• 'FILE' : Read from a file (see the file and file_format parameters)

x_i (default based on model grid)

Inner boundary coordinate of calculation grid

x_o (default based on model grid)

Outer boundary coordinate of calculation grid

w_osc (default 0)

Oscillatory weighting parameter \(w_{\rm osc}\)

w_exp (default 0)

Exponential weighting parameter \(w_{\rm exp}\)

w_ctr (default 0)

Center weighting parameter \(w_{\rm ctr}\)

w_thm (default 0)

Thermal weighting parameter \(w_{\rm thm}\)

w_str (default 0)

Structural weighting parameter \(w_{\rm str}\)

dx_min (default SQRT(EPSILON(1._WP)))

Minimum spacing of grid points

dx_max (default HUGE(0._WP))

Maximum spacing of grid points

n_iter_max (default 32)

Maximum number of refinement iterations

resolve_ctr (default .TRUE.)

Flag to resolve central evanescent region

file (default '')

Name of file containing scaffold grid data (when scaffold_src='FILE')

file_format (default '')

Format of file containing scaffold grid data (when scaffold_src='FILE'); one of:

• 'TEXT': text file with one abscissa value per line

• 'DETAIL': detail file with abscissa values provided in x dataset

tag_list (default '', which matches all)

Comma-separated list of &mode tags to match

See the Spatial Grids section for further details, in particular a discussion of how the weighting (w_*) parameters work.

Stellar Model Parameters

The &model namelist group defines stellar model parameters, as follows:

model_type

Type of model to use; one of:

• 'HOM' : Homogeneous compressible model

• 'POLY' : Polytropic model read from external file

• 'EVOL' : Evolutionary model read from external file

file

Name of file (when model_type= 'POLY'|'EVOL')

file_format

Format of file (when model_type= 'EVOL'); one of

• 'AMDL' : AMDL-format binary file

• 'B3' : B3-format HDF5 file

• 'FAMDL' : FAMDL-format text file

• 'FGONG' : FGONG-format text file

• 'GSM' : GSM-format HDF5 file

• 'LOSC' : LOSC-format text file

• 'MESA' : MESA GYRE-format text file

• 'OSC' : OSC-format text file

• 'WDEC' : WDEC-format text file

data_format (default '', indicates auto-select)

Fortran format specifier for data read from OSC-, FGONG- and FAMDL-format files

deriv_type (default 'MONO')

Cubic interpolation derivatives type (when model_type='POLY'|'EVOL'); one of

• 'NATURAL' : Natural (spline) derivatives

• 'FINDIFF' : Finite-difference derivatives

• 'MONO' : Monotonized derivatives (default)

Gamma_1 (default 5/3)

First adiabatic exponent (when model_type='HOM')

grid_type (default 'UNI')

Model grid type (when model_type='HOM'); one of

• 'UNI' : Uniform spacing

• 'GEO' : Geometric spacing

• 'LOG' : Logarithmic spacing

n (default 10)

Number of points in model grid (when model_type='HOM')

s (default 1)

Skewness parameter for model grid (when model_type='HOM' and grid_type='GEO'|'LOG')

x_i (default 0)

Inner boundary coordinate of model grid (when model_type='HOM')

x_o (default 1)

Outer boundary coordinate of model grid (when model_type='HOM')

dx_snap (default 0)

Threshold for snapping model points together, when model_type is 'EVOL'. If a pair of points are separated by less than dx_snap, they are snapped together.

add_center (default .TRUE.)

Flag to add a center point to the model (when model_type='EVOL'|'POLY'). If a point does not already exist at the origin, then one is added

repair_As (default .FALSE.)

Flag to repair inaccuracies in the dimensionless Brunt-Väisälä frequency at density discontinuities

Mode Parameters

The &mode namelist group defines mode parameters, as follows:

l (default 0)

Harmonic degree \(\ell\)

m (default 0)

Azimuthal order \(m\)

n_pg_min (default -HUGE)

Filter for minimum radial order (applies only to adiabatic calculations)

n_pg_max (default +HUGE)

Filter for maximum radial order (applies only to adiabatic calculations)

tag

Tag for controlling selection of other parameters

Numerical Parameters

The &num namelist group defines numerical method parameters; the input file can contain one or more, but only the last (tag-matching) one is used. Allowable fields are:

diff_scheme (default 'COLLOC_GL2')

Difference equation scheme; one of:

• 'COLLOC_GL2' : Second-order Gauss-Legendre collocation

• 'COLLOC_GL4' : Fourth-order Gauss-Legendre collocation

• 'COLLOC_GL6' : Sixth-order Gauss-Legendre collocation

• 'MAGNUS_GL2' : Second-order Gauss-Legendre Magnus

• 'MAGNUS_GL4' : Fourth-order Gauss-Legendre Magnus

• 'MAGNUS_GL6' : Sixth-order Gauss-Legendre Magnus

• 'MIRK' : Fourth-order mono-implicit Runge-Kutta (experimental)

• 'TRAPZ' : Trapezoidal, with the prescription by Sugimoto (1970) for non-adiabatic cases

r_root_solver (default 'BRENT')

Root solver for real arithmetic; one of:

• 'BRENT' : Brent’s method

c_root_solver (default 'RIDDERS')

Root solver for complex arithmetic; one of

• 'RIDDERS' : Complex Ridders’ method

• 'SECANT' : Secant method

• 'SIMPLEX' : Simplex method

n_iter_max (default 50)

Maximum number of iterations in root-finding algorithm

matrix_type (default 'BLOCK’)

Storage type of system matrix; one of

• 'BAND' : Band-structured

• 'BLOCK' : Block-structured

deflate_roots (default .TRUE.)

Flag to use root deflation, which can avoid the same eigenfrequency being found multiple times

restrict_roots (default .TRUE.)

Flag to check each roots found lies within the bounds of the frequency scan

ad_search (default 'BRACKET')

Initial search method for adiabatic calculations; one of

• 'BRACKET' : Bracket sign changes in the discriminant function

nad_search (default 'AD')

Initial search method for non-adiabatic calculations; one of

• 'AD' : Use adiabatic eigenfrequencies

• 'MINMOD' : Find minima in the modulus of the discriminant function, along the real-\(\omega\) axis

• 'CONTOUR' : Find intersections between real and imaginary zero-contours of the discriminant function

See the Non-Adiabatic Oscillations chapter for more details about these search methods.

tag_list (default '', which matches all)

Comma-separated list of &mode tags to match

Orbital Parameters

The &orbit namelist group defines orbital parameters, as follows:

Omega_orb (default 1)

Orbital angular frequency

Omega_orb_units (default 'NULL')

Units of Omega_orb; one of:

• 'NONE' : Dimensionless angular frequency

• 'HZ' : Linear frequency in Hz[1]

• 'UHZ' : Linear frequency in \(\mu\)Hz[1]

• 'RAD_PER_SEC' : Angular frequency in radians per second[1]

• 'CYC_PER_DAY' : Linear frequency in cycles per day[1]

• 'CRITICAL' : Fraction of the Roche critical rate[1]

q (default 1)

Ratio of secondary mass to primary mass

e (default 0)

Orbital eccentricity

tag_list (default '', which matches all)

Comma-separated list of &tide tags to match

Footnotes

[1] (1,2,3,4,5)

This option is available only for stellar models with D capability

Oscillation Parameters

The &osc namelist group defines oscillation parameters, as follows:

inner_bound (default 'REGULAR')

Inner boundary conditions; one of:

• 'REGULAR' : Regularity-enforcing (only valid when inner grid point is at \(x = 0\))

• 'ZERO_R' : Zero radial displacement (only valid when inner grid point is at \(x \ne 0\))

• 'ZERO_H' : Zero horizontal displacement (only valid when inner grid point is at \(x \ne 0\))

outer_bound (default 'VACUUM')

Outer boundary conditions; one of:

• 'VACUUM' : Vanishing surface density

• 'DZIEM' : Formulation following Dziembowski (1971)

• 'UNNO' : Formulation following Unno et al. (1989)

• 'JCD' : Formulation following Jørgen Christensen-Dalsgaard (ADIPLS)

• 'ISOTHERMAL' : Formulation based on local dispersion analysis for isothermal atmosphere

• 'GAMMA' : Vanishing displacement and derivative at outer boundary, intended for use with \(\gamma\) modes (isolated g modes; see Ong & Basu, 2020)

outer_bound_cutoff (default '')

Outer boundary conditions to use when evaluating cutoff frequencies (see freq_units); same options as outer_bound, and if left blank then takes its value from outer_bound

outer_bound_branch (default 'E_NEG')

Dispersion relation solution branch to use for outer boundary conditions (when outer_bound='UNNO'|'JCD'|'ISOTHERMAL'); one of

• 'E_NEG' : Outward-decaying energy density

• 'E_POS' : Outward-growing energy density

• 'F_NEG' : Outward energy flux

• 'F_POS' : Inward energy flux

• 'V_NEG' : Outward phase velocity

• 'V_POS' : Inward phase velocity

variables_set (default 'GYRE')

Dependent variables in oscillation equations; one of:

• 'GYRE' : GYRE formulation, as described in the Dimensionless Formulation section

• 'DZIEM' : Formulation following Dziembowski (1971)

• 'JCD' : Formulation following Jørgen Christensen-Dalsgaard (ADIPLS)

• 'MIX' : Mixed formulation ('JCD' for \(y_{3,4}\), 'DZIEM' for \(y_{1,2}\))

• 'LAGP' : Lagrangian pressure perturbation formulation

alpha_grv (default 1.)

Scaling factor for gravitational potential perturbations (see the \(\alphagrv\) entry in the Physics Switches section)

alpha_thm (default 1.)

Scaling factor for the thermal timescale (see the \(\alphathm\) entry in the Physics Switches section)

alpha_hfl (default 1.)

Scaling factor for horizontal flux perturbations (see the \(\alphahfl\) entry in the Physics Switches section)

alpha_gam (default 1.)

Scaling factor for g-mode isolation (see the \(\alphagam\) term in entry in the Physics Switches section)

alpha_pi (default 1.)

Scaling factor for p-mode isolation (see the \(\alphapi\) term in entry in the Physics Switches section)

alpha_kar (default 1.)

Scaling factor for opacity density partial derivative (see the \(\alphakar\) entry in the Physics Switches section)

alpha_kat (default 1.)

Scaling factor for opacity temperature partial derivative (see the \(\alphakat\) entry in the Physics Switches section)

alpha_rht (default 0.)

Scaling factor for time-dependent term in radiative heat equation (see the \(\alpharht\) entry in the Physics Switches section)

alpha_trb (default 0.)

Scaling factor for the turbulent mixing length (see the \(\alphatrb\) entry in the Physics Switches section)

inertia_norm (default 'BOTH')

Inertia normalization factor; one of

• 'RADIAL' : Radial amplitude squared, \(|\xi_{\rm r}|^{2}\), evaluated at x_ref

• 'HORIZ' : Horizontal amplitude squared, \(|\lambda| |\xi_{\rm h}|^{2}\), evaluated at x_ref

• 'BOTH' : Overall amplitude squared, \(|\xi_{\rm r}|^{2} + |\lambda| |\xi_{\rm h}|^{2}\), evaluated at x_ref

time_factor (default 'OSC')

Time-dependence factor in pulsation equations; one of:

• 'OSC' : Oscillatory, \(\propto \exp(-{\rm i} \sigma t)\)

• 'EXP' : Exponential, \(\propto \exp(-\sigma t)\)

conv_scheme (default 'FROZEN_PESNELL_1')

Scheme for treating convection; one of:

• 'FROZEN_PESNELL_1' : Freeze convective heating altogether; case 1 described by Pesnell (1990)

• 'FROZEN_PESNELL_4' : Freeze Lagrangian perturbation of convective luminosity; case 4 described by Pesnell (1990)

zeta_scheme (default 'KAWALER')

Scheme for evaluating dimensionless frequency weight function \(\sderiv{\zeta}{x}\) and integral eigenfrequency \(\omega_{\rm int}\); one of:

• 'PESNELL' : Evaluate using eqn. (A5) of Pesnell (1987)

• 'KAWALER' : Evaluate using eqn. (7) of Kawaler et al. (1985), as corrected by Townsend & Kawaler (2023)

• 'KAWALER_GRAV': Evaluate using the g-mode part in eqn. (7) of Kawaler et al. (1985)

• 'DUPRET' : Evaluate using eqn. (1.71) of Dupret (2002, PhD thesis)

deps_scheme (default 'MODEL')

Scheme for calculating nuclear energy generation partials \(\epsnucrho\) and \(\epsnucT\); one of:

• 'MODEL' : Use values from model

• 'FILE' : Use complex (phase-lagged) values from separate file

deps_file (default '')

Name of epsilon partial derivatives file (when deps_scheme='FILE')

deps_file_format (default 'WOLF')

Format of epsilon partial derivative file (when deps_scheme='FILE'); one of:

• 'WOLF' : Format used in preparation of Wolf et al. (2018)

x_ref (default 1 or outer grid point, whichever is smaller)

Reference fractional radius for photosphere, normalizations etc.

x_atm (default -1, implying outer grid point)

Fractional radius for convection-zone crossover point of \(\pi/\gamma\) modes (isolated p and g modes; see Ong & Basu, 2020)

adiabatic (default .TRUE.)

Flag to perform adiabatic calculations

nonadiabatic (default .FALSE.)

Flag to perform non-adiabatic calculations

quasiad_eigfuncs (default .FALSE.)

Flag to calculate quasi-adiabatic entropy/luminosity eigenfunctions during adiabatic calculations

reduce_order (default .TRUE.)

Flag to reduce the order of the adiabatic radial-pulsation equations from 4 to 2

tag_list (default '', which matches all)

Comma-separated list of &mode tags to match

Output Parameters

The &ad_output, &nad_output and &tides_output namelist groups determine the output produced at the end of a run (the first two for the adiabatic and non-adiabatic calculation stages of gyre; the third for gyre_tides). Parameters are as follows:

summary_file (default '')

Name of summary file

summary_file_format (default 'HDF')

Format of summary file; one of

• 'HDF' : HDF5 file

• 'TXT' : Text file

summary_item_list (default 'l,n_pg,omega,freq')

Comma-separated list of output items to write to summary file; see the Summary Files section for possible choices

summary_filter_list (default '')

Comma-separated list of filter criteria for summary files; see the Output Filters section for possible choices

detail_template (default '')

Name template of detail files. Names are generated using the following pattern substitutions:

• '%ID' : Unique mode index, formatted in fixed-width field

• '%id' : Same as '%ID', but formatted in variable-width field

• '%L' : Harmonic degree \(\ell\), formatted in fixed-width field

• '%l' : Same as '%L', but formatted in variable-width field

• '%M' : Azimuthal order \(m\), formatted in fixed-width field

• '%m' : Same as '%M', but formatted in variable-width field

• '%N' : Radial order \(n_{\rm pg}\), formatted in fixed-width field

• '%n' : Same as '%N', but formatted in variable-width field

• '%P' : Acoustic wave winding number \(n_{\rm p}\), formatted in fixed-width field

• '%p' : Same as '%P', but formatted in variable-width field

• '%G' : Gravity wave winding number \(n_{\rm g}\), formatted in fixed-width field

• '%g' : Same as '%G', but formatted in variable-width field

detail_file_format (default 'HDF')

Format of detail files; one of

• 'HDF' : HDF5 file

• 'TXT' : Text file

detail_item_list (default 'l,n_pg,omega,freq,x,xi_r,xi_h')

Comma-separated list of output items to write to detail files; see the Detail Files section for possible choices

detail_filter_list (default '')

Comma-separated list of filter criteria for detail files; see the Output Filters section for possible choices

freq_units (default NONE)

Units of freq output item; one of:

• 'NONE' : Dimensionless angular frequency

• 'HZ' : Linear frequency in Hz[1]

• 'UHZ' : Linear frequency in \(\mu\)Hz[1]

• 'RAD_PER_SEC' : Angular frequency in radians per second[1]

• 'CYC_PER_DAY' : Linear frequency in cycles per day[1]

• 'ACOUSTIC_DELTA' : Fraction of the asymptotic acoustic large frequency separation \(\Delta \nu\)

• 'GRAVITY_DELTA' : Fraction of the asymptotic inverse gravity period separation \((\Delta P)^{-1}\)

• 'UPPER_DELTA' : Greater of \(\Delta \nu\) and \((\Delta P)^{-1}\)

• 'LOWER_DELTA' : Lesser of \(\Delta \nu\) and \((\Delta P)^{-1}\)

• 'ACOUSTIC_CUTOFF' : Fraction of the acoustic cutoff frequency[1]

• 'GRAVITY_CUTOFF' : Fraction of the gravity cutoff frequency[1]

• 'ROSSBY_I' : Fraction of Rossby frequency at inner boundary

• 'ROSSBY_O' : Fraction of Rossby frequency at outer boundary

freq_frame (default INERTIAL)

Frame of freq output item; one of:

• 'INERTIAL' : Inertial frame

• 'COROT_I' : Co-rotating frame at inner boundary

• 'COROT_O' : Co-rotating frame at outer boundary

label (default '')

Textual label to add to all output files

Footnotes

[1] (1,2,3,4,5,6)

This option is available only for stellar models with D capability

Rotation Parameters

The &rot namelist group defines rotational parameters, as follows:

coriolis_method (default 'NULL')

Method used to treat the Coriolis force; one of:

• 'NULL' : Neglect the Coriolis force

• 'TAR' : Use the traditional approximation of rotation

Omega_rot_source (default 'MODEL')

Source for rotational angular frequency \(\Orot\); one of:

• 'MODEL' : Differential rotation, with a spatially varying \(\Orot\) obtained from the stellar model

• 'UNIFORM' : Uniform rotation, with a spatially constant \(\Orot\) set by the Omega_rot and Omega_rot_units parameters

Omega_rot (default 0)

Rotation angular frequency (when Omega_rot_source='UNIFORM')

Omega_rot_units (default 'NULL')

Units of Omega_rot (when Omega_rot_source='UNIFORM'); one of:

• 'NONE' : Dimensionless angular frequency

• 'HZ' : Linear frequency in Hz[1]

• 'UHZ' : Linear frequency in \(\mu\)Hz[1]

• 'RAD_PER_SEC' : Angular frequency in radians per second[1]

• 'CYC_PER_DAY' : Linear frequency in cycles per day[1]

• 'CRITICAL' : Fraction of the Roche critical rate[1]

rossby (default .FALSE.)

Flag to use Rossby solution family in TAR (when coriolis_method='TAR')

complex_lambda (default .FALSE.)

Flag to use complex arithmetic when evaluating the TAR angular eigenvalue \(\lambda\) (when coriolis_method='TAR')

tag_list (default '', which matches all)

Comma-separated list of &mode tags to match

Footnotes

[1] (1,2,3,4,5)

This option is available only for stellar models with D capability

Frequency Scan Parameters

The &scan namelist group defines frequency grid parameters, as follows:

grid_type (default 'LINEAR')

Distribution of frequency points; one of:

• 'LINEAR' : Uniform in frequency

• 'INVERSE' : Uniform in inverse frequency

• 'FILE' : Read from file

grid_frame (default 'INERTIAL')

Reference frame in which grid_type applies; one of:

• 'INERTIAL' : Inertial frame

• 'COROT_I' : Co-rotating frame at inner boundary

• 'COROT_O' : Co-rotating frame at outer boundary

freq_min (default 1)

Minimum frequency, when grid_type is 'LINEAR' or 'INVERSE'

freq_max (default 10)

Maximum frequency, when grid_type is 'LINEAR' or 'INVERSE'

n_freq (default 10)

Number of frequency points, when grid_type is 'LINEAR' or 'INVERSE'

freq_units (default NONE)

Units of freq_min and freq_max, when grid_type is 'LINEAR' or 'INVERSE'; units of read frequencies when grid_type is 'FILE'

• 'NONE' : Dimensionless angular frequency

• 'HZ' : Linear frequency in Hz[1]

• 'UHZ' : Linear frequency in \(\mu\)Hz[1]

• 'RAD_PER_SEC' : Angular frequency in radians per second[1]

• 'CYC_PER_DAY' : Linear frequency in cycles per day[1]

• 'ACOUSTIC_DELTA' : Fraction of the asymptotic acoustic large frequency separation \(\Delta \nu\)

• 'GRAVITY_DELTA' : Fraction of the asymptotic inverse gravity period separation \((\Delta P)^{-1}\)

• 'UPPER_DELTA' : Greater of \(\Delta \nu\) and \((\Delta P)^{-1}\)

• 'LOWER_DELTA' : Lesser of \(\Delta \nu\) and \((\Delta P)^{-1}\)

• 'ACOUSTIC_CUTOFF' : fraction of the acoustic cutoff frequency[1]

• 'GRAVITY_CUTOFF' : fraction of the gravity cutoff frequency[1]

• 'ROSSBY_I' : fraction of Rossby frequency (see eqn. 15) at inner boundary

• 'ROSSBY_O' : fraction of Rossby frequency (see eqn. 15) at outer boundary

freq_min_units (default '')

Units of freq_min; same options as freq_units and overrides it if set

freq_max_units (default '')

Units of freq_max; same options as freq_units and overrides it if set

freq_frame (default 'INERTIAL')

Reference frame in which freq_min and freq_max are defined, when grid_type is 'LINEAR' or 'INVERSE'; one of:

• 'INERTIAL' : Inertial frame

• 'COROT_I' : Co-rotating frame at inner boundary

• 'COROT_O' : Co-rotating frame at outer boundary

file

File to read frequencies from, when grid_type is 'FILE'

axis (default 'REAL’)

Axis that scan applies to; one of

• 'REAL' : Real axis

• 'IMAG' : Imaginary axis

tag_list (default '', which matches all)

Comma-separated list of &mode tags to match

An input file can contain one or more &scan namelist group; the points defined by each (tag-matching) group are merged together to build the frequency grid. See the Frequency Grids section for further details.

Footnotes

[1] (1,2,3,4,5,6)

This option is available only for stellar models with D capability

Tidal Parameters

The &tide namelist group defines tidal parameters, as follows:

y_T_thresh_abs (default 0.)

Absolute threshold on dimensionless tidal potential \(y_{\rm T}\) for a partial tide to contribute

y_T_thresh_rel (default 0.)

Relative threshold on dimensionless tidal potential \(y_{\rm T}\) for a partial tide to contribute

omega_c_thresh (default 0.)

Threshold on dimensionless co-rotating frequency \(\omega_{\rm c}\) for a partial tide to be treated as dynamic (rather than static)

alpha_frq (default 1.)

Scaling parameter \(\alphafrq\) for tidal forcing frequency

l_min (default 2)

Minimum harmonic degree \(\ell\) in spatial expansion of tidal potential

l_max (default 2)

Maximum harmonic degree \(\ell\) in spatial expansion of tidal potential

m_min (default -HUGE)

Minimum azimuthal order \(m\) in spatial expansion of tidal potential

m_max (default HUGE)

Maximum azimuthal order \(m\) in spatial expansion of tidal potential

k_min (default -10)

Minimum orbital harmonic \(k\) in temporal expansion of tidal potential

k_max (default 10)

Maximum orbital harmonic \(k\) in temporal expansion of tidal potential

tag

Tag for controlling selection of other parameters

Output Files

This chapter describes the summary and detail output files written by the GYRE frontends.

Summary Files

Summary files collect together global properties, such as eigenfrequencies and radial orders, of all solutions (modes, tidal responses, etc.) found during a run. The specific data written to a summary file are controlled by the summary_item_list parameters of the &ad_output and &nad_output namelist groups (gyre adiabatic and non-adiabatic calculations, respectively) and the &tides_output namelist group (gyre_tides calculations). These parameters specify the items to be written, via a comma-separated list.

The following subsections describe the items that may appear in summary_item_list, grouped together by functional area. For each item, the corresponding math symbol is given (if there is one), together with the datatype, and a brief description. Units (where applicable) are indicated in brackets [].

Solution Data

Item	Symbol	Datatype	Description
n_row	\(N_{\rm row}\)	integer	number of rows in summary file, each corresponding to a mode found (gyre) or a tidal response evaluated (gyre_tides)
n	\(N\)	integer(n_row)	number of spatial grid points
omega	\(\omega\)	complex(n_row)	dimensionless eigenfrequency

Observables

Item	Symbol	Datatype	Description
freq	—	complex(n_row)	dimensioned frequency; units and reference frame controlled by freq_units and freq_frame parameters
freq_units	—	string	freq_units parameter
freq_frame	—	string	freq_frame parameter
f_T	\(f_{T}\)	real(n_row)	Effective temperature perturbation amplitude; evaluated using eqn. 5 of Dupret et al. (2003)
f_g	\(f_{\rm g}\)	real(n_row)	Effective gravity perturbation amplitude; evaluated using eqn. 6 of Dupret et al. (2003)
psi_T	\(\psi_{T}\)	real(n_row)	Effective temperature perturbation phase; evaluated using eqn. 5 of Dupret et al. (2003)
psi_g	\(\psi_{\rm g}\)	real(n_row)	Effective gravity perturbation phase; evaluated using eqn. 6 of Dupret et al. (2003)

Classification & Validation

Item	Symbol	Datatype	Description
id	—	integer(n_row)	unique mode index
l	\(\ell\)	integer(n_row)	harmonic degree
l_i	\(\ell_{\rm i}\)	complex(n_row)	effective harmonic degree at inner boundary
m	\(m\)	integer(n_row)	azimuthal order
n_p	\(\nump\)	integer(n_row)	acoustic-wave winding number
n_g	\(\numg\)	integer(n_row)	gravity-wave winding number
n_pg	\(\numpg\)	integer(n_row)	radial order within the Eckart-Scuflaire-Osaki-Takata scheme (see Takata, 2006b)
omega_int	\(\omega_{\rm int}\)	complex(n_row)	dimensionless eigenfrequency; evaluated as \(\omega_{\rm int} = \sqrt{\zeta/E}\)
zeta	\(\zeta\)	complex(n_row)	integral of \(\sderiv{\zeta}{x}\) with respect to \(x\)

Perturbations

Item	Symbol	Datatype	Description
x_ref	\(x_{\rm ref}\)	real	fractional radius of reference location
xi_r_ref	\(\txi_{r,{\rm ref}}\)	complex(n_row)	radial displacement perturbation at reference location [\(R\)]
eul_Phi_ref	\(\tPhi'_{\rm ref}\)	complex(n_row)	Eulerian potential perturbation at reference location [\(GM/R\)]
deul_Phi_ref	\((\sderiv{\tPhi'}{x})_{\rm ref}\)	complex(n_row)	Eulerian potential gradient perturbation at reference location [\(GM/R^{2}\)]
lag_S_ref	\(\delta\tS_{\rm ref}\)	complex(n_row)	Lagrangian specific entropy perturbation at reference location [\(R\)]
lag_L_ref	\(\delta\tL_{\rm R,ref}\)	complex(n_row)	Lagrangian radiative luminosity perturbation at reference location [\(L\)]

Energetics & Transport

Item	Symbol	Datatype	Description
eta[1]	\(\eta\)	real(n_row)	normalized growth rate \(\eta\); evaluated using expression in text of page 1186 of Stellingwerf (1978)
E	\(E\)	real(n_row)	mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\)
E_p	\(E_{\rm p}\)	real(n_row)	acoustic mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\) where \(\varpi=1\)
E_g	\(E_{\rm g}\)	real(n_row)	gravity mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\) in regions where \(\varpi=-1\)
E_norm	\(E_{\rm norm}\)	real(n_row)	normalized inertia; evaluation controlled by inertia_norm parameter
E_ratio	—	real(n_row)	ratio of mode inertia outside reference location, to total inertia
H	\(H\)	real(n_row)	mode energy [\(G M^{2}/R\)]; evaluated as \(\frac{1}{2} \omega^{2} E\)
W[1]	\(W\)	real(n_row)	mode work [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{W}{x}\)
W_eps[1]	\(W_{\epsilon}\)	real(n_row)	mode work [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{W_{\epsilon}}{x}\)
tau_ss	\(\tau_{\rm ss}\)	real(n_row)	steady-state torque [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{\tau_{\rm ss}}{x}\)
tau_tr	\(\tau_{\rm tr}\)	real(n_row)	steady-state torque [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{\tau_{\rm tr}}{x}\)

Rotation

Item	Symbol	Datatype	Description
Omega_rot_ref	\(\Omega_{\rm rot,ref}\)	real(n_row)	rotation angular frequency at reference location[\(\sqrt{GM/R^{3}}\)]
domega_rot	\(\Delta \omega\)	real(n_row)	dimensionless first-order rotational splitting; evaluated using eqn. 3.355 of Aerts et al. (2010)
dfreq_rot	—	real(n_row)	dimensioned first-order rotational splitting; units and reference frame controlled by freq_units and freq_frame parameters
beta	\(\beta\)	real(n_row)	rotation splitting coefficient; evaluated by integrating \(\sderiv{\beta}{x}\)

Stellar Structure

Item	Symbol	Datatype	Description
M_star[2]	\(M\)	real(n_row)	stellar mass [\(\gram\)]
R_star[2]	\(R\)	real(n_row)	stellar radius [\(\cm\)]
L_star[2]	\(L\)	real(n_row)	stellar luminosity [\(\erg\,\second^{-1}\)]
Delta_p	\(\Delta \nu\)	real(n_row)	asymptotic p-mode large frequency separation [\(\sqrt{GM/R^{3}}\)]
Delta_g	\((\Delta P)^{-1}\)	real(n_row)	asymptotic g-mode inverse period separation [\(\sqrt{GM/R^{3}}\)]

Tidal Response

Note that these items are available only when using gyre_tides.

Item	Symbol	Datatype	Description
k	\(k\)	integer(n_row)	Fourier harmonic
eul_Psi_ref	\(\tPsi'_{\rm ref}\)	complex(n_row)	Eulerian total potential perturbation at reference location [\(GM/R\)]
Phi_T_ref	\(\tPhi_{\rm T, ref}\)	real(n_row)	tidal potential at reference location [\(GM/R\)]
Omega_orb	\(\Oorb\)	real(n_row)	orbital angular frequency; units and reference frame controlled by freq_units and freq_frame parameters
q	\(q\)	real(n_row)	ratio of secondary mass to primary mass
e	\(e\)	real(n_row)	orbital eccentricity
R_a	\(R/a\)	real(n_row)	ratio of primary radius to orbital semi-major axis
cbar	\(\cbar_{\ell,m,k}\)	real(n_row)	tidal expansion coefficient; see eqn. A1 of Sun et al. (2023)
Gbar_1	\(\Gbar^{(1)}_{\ell,m,k}\)	real(n_row)	secular orbital evolution coefficient; equivalent to \(G^{(1)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_2	\(\Gbar^{(2)}_{\ell,m,k}\)	real(n_row)	secular orbital evolution coefficient; equivalent to \(G^{(2)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_3	\(\Gbar^{(3)}_{\ell,m,k}\)	real(n_row)	secular orbital evolution coefficient; equivalent to \(G^{(3)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_4	\(\Gbar^{(4)}_{\ell,m,k}\)	real(n_row)	secular orbital evolution coefficient; equivalent to \(G^{(4)}_{\ell,m,-k}\) (see Willems et al., 2003)

Footnotes

[1] (1,2,3)

This item is available only for stellar models with N capability

[2] (1,2,3)

This item is available only for stellar models with D capability

Detail Files

Detail files store spatial quantities, such as eigenfunctions and differential inertias, for an individual solution (mode, tidal response, etc.) found during a run. The specific data written to detail files are controlled by the detail_item_list parameters of the &ad_output and &nad_output namelist groups (gyre adiabatic and non-adiabatic calculations, respectively) and the &tides_output namelist group (gyre_tides calculations). These parameters specify the items to be written, via a comma-separated list.

The following subsections describe the items that may appear in detail_item_list, grouped together by functional area. For each item, the corresponding math symbol is given (if there is one), together with the datatype, and a brief description. Units (where applicable) are indicated in brackets [].

Solution Data

Item	Symbol	Datatype	Description
n	\(N\)	integer	number of spatial grid points
omega	\(\omega\)	complex	dimensionless eigenfrequency (gyre) or forcing frequency (gyre_tides)
x	\(x\)	real(n)	independent variable \(x = r/R\)
dx_min	\(\Delta x_{\rm min}\)	real	minimum spacing of spatial grid
dx_max	\(\Delta x_{\rm max}\)	real	maximum spacing of spatial grid
dx_rms	\(\Delta x_{\rm rms}\)	real	root-mean-square spacing of spatial grid
x_ref	\(x_{\rm ref}\)	real	fractional radius of reference location
y_1	\(y_{1}\)	complex(n)	dependent variable
y_2	\(y_{2}\)	complex(n)	dependent variable
y_3	\(y_{3}\)	complex(n)	dependent variable
y_4	\(y_{4}\)	complex(n)	dependent variable
y_5	\(y_{5}\)	complex(n)	dependent variable
y_6	\(y_{6}\)	complex(n)	dependent variable

The definitions of the dependent variables \(\{y_{1},\ldots,y_{6}\}\) are provided in the Oscillation Equations chapter.

Observables

Item	Symbol	Datatype	Description
freq	—	complex	dimensioned frequency; units and reference frame controlled by freq_units and freq_frame parameters
freq_units	—	string	freq_units parameter
freq_frame	—	string	freq_frame parameter
f_T	\(f_{T}\)	real	Effective temperature perturbation amplitude; evaluated using eqn. 5 of Dupret et al. (2003)
f_g	\(f_{\rm g}\)	real	Effective gravity perturbation amplitude; evaluated using eqn. 6 of Dupret et al. (2003)
psi_T	\(\psi_{T}\)	real	Effective temperature perturbation phase; evaluated using eqn. 5 of Dupret et al. (2003)
f_g	\(\psi_{\rm g}\)	real	Effective gravity perturbation phase; evaluated using eqn. 6 of Dupret et al. (2003)

Classification & Validation

Item	Symbol	Datatype	Description
id	—	integer	unique mode index
l	\(\ell\)	integer	harmonic degree
l_i	\(\ell_{\rm i}\)	complex	effective harmonic degree at inner boundary
m	\(m\)	integer	azimuthal order
n_p	\(\nump\)	integer	acoustic-wave winding number
n_g	\(\numg\)	integer	gravity-wave winding number
n_pg	\(\numpg\)	integer	radial order within the Eckart-Scuflaire-Osaki-Takata scheme (see Takata, 2006b)
omega_int	\(\omega_{\rm int}\)	complex	dimensionless eigenfrequency; evaluated as \(\omega_{\rm int} = \sqrt{\zeta/E}\)
dzeta_dx	\(\sderiv{\zeta}{x}\)	complex(n)	dimensionless frequency weight function; controlled by zeta_scheme parameter
zeta	\(\zeta\)	complex	integral of \(\sderiv{\zeta}{x}\) with respect to \(x\)
Yt_1	\(\mathcal{Y}_{1}\)	complex(n)	primary eigenfunction for Takata classification; evaluated using a rescaled eqn. 69 of Takata (2006b)
Yt_2	\(\mathcal{Y}_{2}\)	complex(n)	secondary eigenfunction for Takata classification; evaluated using a rescaled eqn. 70 of Takata (2006b)
I_0	\(I_{0}\)	complex(n)	first integral for radial modes; evaluated using eqn. 42 of Takata (2006a)
I_1	\(I_{1}\)	complex(n)	first integral for dipole modes; evaluated using eqn. 43 of Takata (2006a)
prop_type	\(\varpi\)	integer(n)	propagation type; \(\varpi = 1\) in acoustic-wave regions, \(\varpi=-1\) in gravity-wave regions, and \(\varpi=0\) in evanescent regions

Perturbations

Item	Symbol	Datatype	Description
xi_r_ref	\(\txi_{r,{\rm ref}}\)	complex	radial displacement perturbation at reference location [\(R\)]
xi_h_ref	\(\txi_{\rm h,ref}\)	complex	radial displacement perturbation at reference location [\(R\)]
eul_Phi_ref	\(\tPhi'_{\rm ref}\)	complex	Eulerian potential perturbation at reference location [\(GM/R\)]
deul_Phi_ref	\((\sderiv{\tPhi'}{x})_{\rm ref}\)	complex	Eulerian potential gradient perturbation at reference location [\(GM/R^{2}\)]
lag_S_ref	\(\delta\tS_{\rm ref}\)	complex	Lagrangian specific entropy perturbation at reference location [\(R\)]
lag_L_ref	\(\delta\tL_{\rm R,ref}\)	complex	Lagrangian radiative luminosity perturbation at reference location [\(L\)]
xi_r	\(\txir\)	complex(n)	radial displacement perturbation [\(R\)]
xi_h	\(\txih\)	complex(n)	radial displacement perturbation [\(R\)]
eul_Phi	\(\tPhi'\)	complex(n)	Eulerian potential perturbation [\(GM/R\)]
deul_Phi	\(\sderiv{\tPhi'}{x}\)	complex(n)	Eulerian potential gradient perturbation [\(GM/R^{2}\)]
lag_S	\(\delta\tS\)	complex(n)	Lagrangian specific entropy perturbation [\(\cP\)]
lag_L	\(\delta\tLrad\)	complex(n)	Lagrangian radiative luminosity perturbation [\(L\)]
eul_P	\(\tP'\)	complex(n)	Eulerian total pressure perturbation [\(P\)]
eul_rho	\(\trho'\)	complex(n)	Eulerian density perturbation [\(\rho\)]
eul_T	\(\tT'\)	complex(n)	Eulerian temperature perturbation [\(T\)]
lag_P	\(\delta\tP\)	complex(n)	Lagrangian total pressure perturbation [\(P\)]
lag_rho	\(\delta\trho\)	complex(n)	Lagrangian density perturbation [\(\rho\)]
lag_T	\(\delta\tT\)	complex(n)	Lagrangian temperature perturbation [\(T\)]

Energetics & Transport

Item	Symbol	Datatype	Description
eta	\(\eta\)	real	normalized growth rate \(\eta\); evaluated using expression in text of page 1186 of Stellingwerf (1978)
E	\(E\)	real	mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\)
E_p	\(E_{\rm p}\)	real	acoustic mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\) where \(\varpi=1\)
E_g	\(E_{\rm g}\)	real	gravity mode inertia [\(M R^{2}\)]; evaluated by integrating \(\sderiv{E}{x}\) in regions where \(\varpi=-1\)
E_norm	\(E_{\rm norm}\)	real	normalized inertia; evaluation controlled by inertia_norm parameter
E_ratio		real	ratio of mode inertia outside reference location, to total inertia
H	\(H\)	real	mode energy [\(G M^{2}/R\)]; evaluated as \(\frac{1}{2} \omega^{2} E\)
W	\(W\)	real	mode work [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{W}{x}\)
W_eps	\(W_{\epsilon}\)	real	mode work [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{W_{\epsilon}}{x}\)
tau_ss	\(\tau_{\rm ss}\)	real	steady-state torque [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{\tau_{\rm ss}}{x}\)
tau_tr	\(\tau_{\rm tr}\)	real	steady-state torque [\(G M^{2}/R\)]; evaluated by integrating \(\sderiv{\tau_{\rm tr}}{x}\)
dE_dx	\(\sderiv{E}{x}\)	real(n)	differential inertia [\(M R^{2}\)]; evaluated using eqn. 3.139 of Aerts et al. (2010)
dW_dx[1]	\(\sderiv{W}{x}\)	real(n)	differential work [\(GM^{2}/R\)]; evaluated using eqn. 25.9 of Unno et al. (1989)
dW_eps_dx[1]	\(\sderiv{W_{\epsilon}}{x}\)	real(n)	differential nuclear work [\(GM^{2}/R\)]; evaluated using eqn. 25.9 of Unno et al. (1989)
dtau_ss_dx	\(\sderiv{\tau_{\rm ss}}{x}\)	real(n)	steady-state differential torque [G M^{2}/R]
dtau_tr_dx	\(\sderiv{\tau_{\rm tr}}{x}\)	real(n)	transient differential torque [G M^{2}/R]
alpha_0	\(\alpha_{0}\)	real(n)	excitation coefficient; evaluated using eqn. 26.10 of Unno et al. (1989)
alpha_1	\(\alpha_{1}\)	real(n)	excitation coefficient; evaluated using eqn. 26.12 of Unno et al. (1989)

Rotation

Item	Symbol	Datatype	Description
Omega_rot_ref	\(\Omega_{\rm rot,ref}\)	real	rotation angular frequency at reference location[\(\sqrt{GM/R^{3}}\)]
Omega_rot	\(\Orot\)	real(n)	rotation angular frequency [\(\sqrt{GM/R^{3}}\)]
domega_rot	\(\Delta \omega\)	real	dimensionless first-order rotational splitting; evaluated using eqn. 3.355 of Aerts et al. (2010)
dfreq_rot	—	real	dimensioned first-order rotational splitting; units and reference frame controlled by freq_units and freq_frame parameters
beta	\(\beta\)	real	rotation splitting coefficient; evaluated by integrating \(\sderiv{\beta}{x}\)
dbeta_dx	\(\sderiv{\beta}{x}\)	complex(n)	unnormalized rotation splitting kernel; evaluated using eqn. 3.357 of Aerts et al. (2010)
lambda	\(\lambda\)	complex(n)	tidal equation eigenvalue

Stellar Structure

Item	Symbol	Datatype	Description
M_star[2]	\(M\)	real	stellar mass [\(\gram\)]
R_star[2]	\(R\)	real	stellar radius [\(\cm\)]
L_star[2]	\(L\)	real	stellar luminosity [\(\erg\,\second^{-1}\)]
Delta_p	\(\Delta \nu\)	real	asymptotic p-mode large frequency separation [\(\sqrt{GM/R^{3}}\)]
Delta_g	\((\Delta P)^{-1}\)	real	asymptotic g-mode inverse period separation [\(\sqrt{GM/R^{3}}\)]
V_2	\(V_2\)	real(n)	structure coefficient; defined in Structure Coefficients section
As	\(A^{*}\)	real(n)	structure coefficient; defined in Structure Coefficients section
U	\(U\)	real(n)	structure coefficient; defined in Structure Coefficients section
c_1	\(c_{1}\)	real(n)	structure coefficient; defined in Structure Coefficients section
Gamma_1	\(\Gammi\)	real(n)	adiabatic exponent; defined in Linearized Equations section
nabla[1]	\(\nabla\)	real(n)	temperature gradient; defined in Structure Coefficients section Dimensionless Formulation section
nabla_ad[1]	\(\nabad\)	real(n)	adiabatic temperature gradient; defined in Linearized Equations section
dnabla_ad[1]	\(\dnabad\)	real(n)	derivative of adiabatic temperature gradient
upsilon_T[1]	\(\upsT\)	real(n)	thermodynamic coefficient; defined in Linearized Equations section
c_lum[1]	\(\clum\)	real(n)	structure coefficient; defined in Structure Coefficients section
c_rad[1]	\(\crad\)	real(n)	structure coefficient; defined in Structure Coefficients section
c_thn[1]	\(\cthn\)	real(n)	structure coefficient; defined in Structure Coefficients section
c_thk[1]	\(\cthk\)	real(n)	structure coefficient; defined in Structure Coefficients section
c_eps[1]	\(\ceps\)	real(n)	structure coefficient; defined in Structure Coefficients section
kap_rho[1]	\(\kaprho\)	real(n)	opacity partial; defined in Linearized Equations section
kap_T[1]	\(\kapT\)	real(n)	opacity partial; defined in Linearized Equations section
eps_rho[1]	\(\epsnucrho\)	real(n)	nuclear energy generation partial; defined in Linearized Equations section
eps_T[1]	\(\epsnucT\)	real(n)	nuclear energy generation partial; defined in Linearized Equations section
M_r[2]	\(M_r\)	real(n)	interior mass [\(\gram\)]
P[2]	\(P\)	real(n)	total pressure [\(\barye\)]
rho[2]	\(\rho\)	real(n)	density [\(\gram\,\cm^{-3}\)]
T[2]	\(T\)	real(n)	temperature [\(\kelvin\)]

Tidal Response

Note that these items are available only when using gyre_tides.

Item	Symbol	Datatype	Description
k	\(k\)	integer	Fourier harmonic
eul_Psi_ref	\(\tPsi'_{\rm ref}\)	complex	Eulerian total potential perturbation at reference location [\(GM/R\)]
Phi_T_ref	\(\tPhi_{\rm T, ref}\)	real	tidal potential at reference location [\(GM/R\)]
eul_Psi	\(\tPsi'\)	complex(n)	Eulerian total potential perturbation [\(GM/R\)]
Phi_T	\(\tPhi_{{\rm T}}\)	real(n)	tidal potential [\(GM/R\)]
Omega_orb	\(\Oorb\)	real	orbital angular frequency; units and reference frame controlled by freq_units and freq_frame parameters
q	\(q\)	real	ratio of secondary mass to primary mass
e	\(e\)	real	orbital eccentricity
R_a	\(R/a\)	real	ratio of primary radius to orbital semi-major axis
cbar	\(\cbar_{\ell,m,k}\)	real	tidal expansion coefficient; see eqn. A1 of Sun et al. (2023)
Gbar_1	\(\Gbar^{(1)}_{\ell,m,k}\)	real	secular orbital evolution coefficient; equivalent to \(G^{(1)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_2	\(\Gbar^{(2)}_{\ell,m,k}\)	real	secular orbital evolution coefficient; equivalent to \(G^{(2)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_3	\(\Gbar^{(3)}_{\ell,m,k}\)	real	secular orbital evolution coefficient; equivalent to \(G^{(3)}_{\ell,m,-k}\) (see Willems et al., 2003)
Gbar_4	\(\Gbar^{(4)}_{\ell,m,k}\)	real	secular orbital evolution coefficient; equivalent to \(G^{(4)}_{\ell,m,-k}\) (see Willems et al., 2003)

Footnotes

[1] (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)

This option is available only for stellar models with N capability

[2] (1,2,3,4,5,6,7)

This option is available only for stellar models with D capability

File Formats

The format of summary and detail files depends on the value of the summary_file_format and detail_file_format parameters in the &ad_output and &nad_output namelist groups (see the Output Parameters section). Possible choices are:

• 'HDF' : A binary format based on the HDF5 format

• 'TXT' : A text format modeled after MESA’s profile file format

For both formats, the data stored in the files come in two flavors — scalars (a single value) and arrays (a sequence of values). Files in either format can be read in Python using the pygyre.read_output() function from PyGYRE (see the Interpreting Output Files chapter for examples).

HDF Format

HDF-format output files adhere to the following conventions:

• All data objects are attached to the root HDF5 group (/)

• Attributes are used to store scalar data

• Datasets are used to store array data

• Real values are written with type H5T_IEEE_F64LE when GYRE is compiled in double precision (the default), and type H5T_IEEE_F32LE otherwise

• Integer values are written with type H5T_STD_I32LE

• Complex values are written as a compound type, composed of a real component re and an imaginary component im; the types of these components are the same as for real values

TXT Format

TXT-format files adhere to the following conventions:

• The first three lines contain the scalar data:

  • The first line contains the column numbers for the scalar data, starting at 1

  • The second line contains the column labels for the scalar data

  • The third line contains the actual scalar data values

• The subsequent lines contain the array data:

  • The fourth line contains the column numbers for the array data, starting at 1

  • The fifth line contains the the column labels for the array data

  • The sixth and subsequent lines contain the actual array data (one line per array element)

• Complex values are written as two columns, with the first column containing the real component and the second the imaginary component

Output Filters

TBD

Stellar Models

This chapter documents the different types of stellar models that can be used with GYRE frontends to specify the equilibrium stellar configuration.

Evolutionary Models

Setting the model_type parameter of the &model namelist group to 'EVOL' tells the frontend to read the equilibrium stellar model from a file created by a stellar evolution code (e.g., MESA).

Supported Formats

The format of the model file is specified by the file_format parameter of the &model namelist group (see the Stellar Model Parameters section). Possible choices are summarized in the table below.

file_format	Description
'AMDL'	Binary file describing an evolutionary model in AMDL format, as reverse engineered from the ADIPLS stellar oscillation code (Christensen-Dalsgaard, 2008)
'B3'	HDF5 file describing an evolutionary model in B3 format. This format is for testing purposes only, and will eventually be superseded and/or removed
'FAMDL'	Text file describing an evolutionary model in FAMDL format, as specified in the CoRoT/ESTA File Formats document
'FGONG'	Text file describing an evolutionary model in FGONG format, as specified in the updated FGONG Format document
'GSM'	HDF5 file describing an evolutionary model in GYRE Stellar Model (GSM) format, as specified in the GSM File Format section
'MESA'	Text file describing an evolutionary model in MESA format, as specified in the MESA File Format section
'LOSC'	Text file describing an evolutionary model in the revised LOSC format
'OSC'	Text file describing an evolutionary model in OSC format, as specified in the CoRoT/ESTA File Formats document)
'WDEC'	Text file describing an evolutionary model in WDEC format, as specified in Bischoff-Kim & Montgomery (2018)

Interpolation

Cubic spline interpolation is used to evaluate data between model grid points. The deriv_type parameter in the &model namelist group controls how the spline derivatives are set up.

Double Points

If a model contains a pair of adjacent points with the same radial coordinate \(r\), this pair is treated as a double point representing a discontinuity in the density and some other thermodynamic quantities (but not the pressure or temperature). GYRE does not attempt to interpolate across double points, but instead handles them properly when solving equations through the use of internal boundary conditions.

Polytropic Models

Setting the model_type parameter of the &model namelist group to 'POLY' tells the frontend to read the equilibrium stellar model from a POLY file created by the build_poly tool. See the Composite Polytropes and Building POLY Models appendices for further details.

Homogeneous Models

Setting the model_type parameter of the &model namelist group to 'HOM' tells the frontend to create a homogeneous (uniform density) stellar model, equivalent to a polytrope with index \(n=0\). Because the structure of these model can be computed analytically, there is no need to read from an external file.

The Gamma_1 parameter of the &model namelist group controls the first adiabatic index of the model, while the n, s and grid_type parameters control the model grid. See the Stellar Model Parameters section for further details.

Model Capabilities

Which data items are included in a given stellar model dictates what sorts of calculation can be performed on that model by the frontends. To this end, the different model types and file formats can be classified according to their capabilities (labeled using a single letter):

N

The model supports non-adiabatic calculations.

D

The model supports dimensioned results.

R

The model supports differential rotational.

The table below summarizes the different capabilities of each model-type and file-format combination.

Model Type	File Format	N	D	R
EVOL	AMDL		X
EVOL	B3	X	X
EVOL	FAMDL		X
EVOL	FGONG		X
EVOL	GSM	X	X	X
EVOL	LOSC		X
EVOL	MESA	X	X	X
EVOL	OSC	X	X	X
EVOL	WDEC		X
POLY	POLY
HOM	—

MESA File Format

Files in MESA format store ASCII text data describing a MESA stellar model (note that MESA itself refers to these files as ‘GYRE-format’ files). To create one of these files in MESA, set the pulse_data_format parameter of the &controls namelist group to the value 'GYRE'.

There are a number of versions of the MESA format, distinguished by the initial header line:

Version 0.01

The first line of version-0.01 MESA-format files is a header with the following columns:

Column	Symbol	Datatype	Definition
1	\(N\)	integer	number of grid points
2	\(M\)	real	stellar mass [\(\gram\)]
3	\(R\)	real	photospheric radius [\(\cm\)]
4	\(L\)	real	photospheric luminosity [\(\erg\,\second^{-1}\)]

The subsequent \(N\) lines contain the model data, one line per grid point extending from the center to the surface, with the following columns:

Column	Symbol	Datatype	Definition
1	\(k\)	integer	grid point index (\(k=1,\ldots,N\))
2	\(r\)	real	radial coordinate [\(\cm\)]
3	\(\frac{M_{r}}{M-M_{r}}\)	real	transformed interior mass
4	\(L_{r}\)	real	interior luminosity [\(\erg\,\second^{-1}\)]
5	\(P\)	real	total pressure [\(\barye\)]
6	\(T\)	real	temperature [\(\kelvin\)]
7	\(\rho\)	real	density [\(\gram\,\cm^{-3}\)]
8	\(\nabla\)	real	temperature gradient
9	\(N^{2}\)	real	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
10	\(\cV\)	real	specific heat at constant volume [\(\erg\,\gram^{-1}\,\kelvin^{-1}\))
11	\(\cP\)	real	specific heat at constant pressure [\(\erg\,\gram^{-1}\,\kelvin^{-1}\))
12	\(\chi_{T}\)	real	equation-of-state partial \((\spderiv{\ln P}{\ln T})_{\rho}\)
13	\(\chi_{\rho}\)	real	equation-of-state partial \((\spderiv{\ln P}{\ln \rho})_{T}\)
14	\(\kappa\)	real	opacity [\(\cm^{2}\,\gram^{-1}\)]
15	\(\kapT\)	real	opacity partial
16	\(\kaprho\)	real	opacity partial
17	\(\epsilon\)	real	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
18	\(\epsnuc\,\epsnucT\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
19	\(\epsnuc\,\epsnucrho\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]

Version 0.19

The first line of version-0.19 MESA-format files is a header with the following columns:

Column	Symbol	Datatype	Definition
1	\(N\)	integer	number of grid points
2	\(M\)	real	stellar mass [\(\gram\)]
3	\(R\)	real	photospheric radius [\(\cm\)]
4	\(L\)	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
5	19	integer	version number

The subsequent \(N\) lines contain the model data, one line per grid point extending from the center to the surface, with the following columns:

Column	Symbol	Datatype	Definition
1	\(k\)	integer	grid point index (\(k=1,\ldots,N\))
2	\(r\)	real	radial coordinate [\(\cm\)]
3	\(\frac{M_{r}}{M-M_{r}}\)	real	transformed interior mass
4	\(L_{r}\)	real	interior luminosity [\(\erg\,\second^{-1}\)]
5	\(P\)	real	total pressure [\(\barye\)]
6	\(T\)	real	temperature [\(\kelvin\)]
7	\(\rho\)	real	density [\(\gram\,\cm^{-3}\)]
8	\(\nabla\)	real	temperature gradient
9	\(N^{2}\)	real	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
10	\(\Gamma_{1}\)	real	adiabatic exponent
11	\(\nabla_{\rm ad}\)	real	adiabatic temperature gradient
12	\(\upsT\)	real	thermodynamic coefficient
13	\(\kappa\)	real	opacity [\(\cm^{2}\,\gram^{-1}\)]
14	\(\kapT\)	real	opacity partial
15	\(\kaprho\)	real	opacity partial
16	\(\epsilon\)	real	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
17	\(\epsnuc\,\epsnucT\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
18	\(\epsnuc\,\epsnucrho\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
19	\(\Orot\)	real	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.00

The first line of version-1.00 MESA-format files is a header with the following columns:

Column	Symbol	Datatype	Definition
1	\(N\)	integer	number of grid points
2	\(M\)	real	stellar mass [\(\gram\)]
3	\(R\)	real	photospheric radius [\(\cm\)]
4	\(L\)	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
5	100	integer	version number

The subsequent \(N\) lines contain the model data, one line per grid point extending from the center to the surface, with the following columns:

Column	Symbol	Datatype	Definition
1	\(k\)	integer	grid point index (\(k=1,\ldots,N\))
2	\(r\)	real	radial coordinate [\(\cm\)]
3	\(M_{r}\)	real	interior mass [\(\gram\)]
4	\(L_{r}\)	real	interior luminosity [\(\erg\,\second^{-1}\)]
5	\(P\)	real	total pressure [\(\barye\)]
6	\(T\)	real	temperature [\(\kelvin\)]
7	\(\rho\)	real	density [\(\gram\,\cm^{-3}\)]
8	\(\nabla\)	real	dimensionless temperature gradient
9	\(N^{2}\)	real	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
10	\(\Gamma_{1}\)	real	adiabatic exponent
11	\(\nabla_{\rm ad}\)	real	adiabatic temperature gradient
12	\(\upsT\)	real	thermodynamic coefficient
13	\(\kappa\)	real	opacity [\(\cm^{2}\,\gram^{-1}\)]
14	\(\kappa\,\kapT\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
15	\(\kappa\,\kaprho\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
16	\(\epsilon\)	real	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
17	\(\epsnuc\,\epsnucT\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
18	\(\epsnuc\,\epsnucrho\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
19	\(\Orot\)	real	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.01

The first line of version-1.01 MESA-format files is a header with the following columns:

Column	Symbol	Datatype	Definition
1	\(N\)	integer	number of grid points
2	\(M\)	real	stellar mass [\(\gram\)]
3	\(R\)	real	photospheric radius [\(\cm\)]
4	\(L\)	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
5	101	integer	version number

The subsequent \(N\) lines contain the model data, one line per grid point extending from the center to the surface, with the following columns:

Column	Symbol	Datatype	Definition
1	\(k\)	integer	grid point index (\(k=1,\ldots,N\))
2	\(r\)	real	radial coordinate [\(\cm\)]
3	\(M_{r}\)	real	interior mass [\(\gram\)]
4	\(L_{r}\)	real	interior luminosity [\(\erg\,\second^{-1}\)]
5	\(P\)	real	total pressure [\(\barye\)]
6	\(T\)	real	temperature [\(\kelvin\)]
7	\(\rho\)	real	density [\(\gram\,\cm^{-3}\)]
8	\(\nabla\)	real	dimensionless temperature gradient
9	\(N^{2}\)	real	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
10	\(\Gamma_{1}\)	real	adiabatic exponent
11	\(\nabla_{\rm ad}\)	real	adiabatic temperature gradient
12	\(\upsT\)	real	thermodynamic coefficient
13	\(\kappa\)	real	opacity [\(\cm^{2}\,\gram^{-1}\)]
14	\(\kappa\,\kapT\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
15	\(\kappa\,\kaprho\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
16	\(\epsnuc\)	real	nuclear energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
17	\(\epsnuc\,\epsnucT\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
18	\(\epsnuc\,\epsnucrho\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
19	\(\Orot\)	real	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.20

The first line of version-1.20 MESA-format files is a header with the following columns:

Column	Symbol	Datatype	Definition
1	\(N\)	integer	number of grid points
2	\(M\)	real	stellar mass [\(\gram\)]
3	\(R\)	real	photospheric radius [\(\cm\)]
4	\(L\)	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
5	120	integer	version number

The subsequent \(N\) lines contain the model data, one line per grid point extending from the center to the surface, with the following columns:

Column	Symbol	Datatype	Definition
1	\(k\)	integer	grid point index (\(k=1,\ldots,N\))
2	\(r\)	real	radial coordinate [\(\cm\)]
3	\(M_{r}\)	real	interior mass [\(\gram\)]
4	\(L_{r}\)	real	interior luminosity [\(\erg\,\second^{-1}\)]
5	\(P\)	real	total pressure [\(\barye\)]
6	\(T\)	real	temperature [\(\kelvin\)]
7	\(\rho\)	real	density [\(\gram\,\cm^{-3}\)]
8	\(\nabla\)	real	dimensionless temperature gradient
9	\(N^{2}\)	real	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
10	\(\Gamma_{1}\)	real	adiabatic exponent
11	\(\nabla_{\rm ad}\)	real	adiabatic temperature gradient
12	\(\upsT\)	real	thermodynamic coefficient
13	\(\kappa\)	real	opacity [\(\cm^{2}\,\gram^{-1}\)]
14	\(\kappa\,\kapT\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
15	\(\kappa\,\kaprho\)	real	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
16	\(\epsnuc\)	real	nuclear energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
17	\(\epsnuc\,\epsnucT\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
18	\(\epsnuc\,\epsnucrho\)	real	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
19	\(\epsgrav\)	real	gravothermal energy release rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
20	\(\Orot\)	real	rotation angular frequency [\(\radian\,\second^{-1}\)]

GSM File Format

Files in GSM (GYRE Stellar Model) format store HDF5 data describing a stellar model. This format is intended as a portable, storage-efficient alternative to the MESA File Format. To create one of these files in MESA (from revision 21.12.1 onward), set the pulse_data_format parameter of the &controls namelist group to the value 'GSM'.

The GSM format adheres to the following conventions:

• All data objects are attached to the root HDF5 group (/)

• Real values are written with type H5T_IEEE_F64LE when GYRE is compiled in double precision (the default), and type H5T_IEEE_F32LE otherwise

• Integer values are written with type H5T_STD_I32LE

There are a number of versions of the GSM format, distinguished by the version attribute in the root HDF5 group:

Version 0.00

Data items in the root HDF5 group of version-0.00 GSM-format files are as follows:

Item	Symbol	Object type	Data type	Definition
n	\(N\)	attribute	integer	number of grid points
M_star	\(M\)	attribute	real	stellar mass [\(\gram\)]
R_star	\(R\)	attribute	real	photospheric radius [\(\cm\)]
L_star	\(L\)	attribute	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
r	\(r\)	dataset	real (n)	radial coordinate [\(\cm\)]
w	\(\frac{M_{r}}{M-M_{r}}\)	dataset	real (n)	transformed interior mass
L_r	\(L_{r}\)	dataset	real (n)	interior luminosity [\(\erg\,\second^{-1}\)]
p	\(P\)	dataset	real (n)	total pressure [\(\barye\)]
rho	\(\rho\)	dataset	real (n)	density [\(\gram\,\cm^{-3}\)]
T	\(T\)	dataset	real (n)	temperature [\(\kelvin\)]
N2	\(N^{2}\)	dataset	real (n)	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
Gamma_1	\(\Gamma_{1}\)	dataset	real (n)	adiabatic exponent
nabla_ad	\(\nabla_{\rm ad}\)	dataset	real (n)	adiabatic temperature gradient
delta	\(\upsT\)	dataset	real (n)	thermodynamic coefficient
nabla	\(\nabla\)	dataset	real (n)	temperature gradient
kappa	\(\kappa\)	dataset	real (n)	opacity [\(\cm^{2}\,\gram^{-1}\)]
kappa_T	\(\kapT\)	dataset	real (n)	opacity partial
kappa_rho	\(\kaprho\)	dataset	real (n)	opacity partial
epsilon	\(\epsilon\)	dataset	real (n)	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
epsilon_T	\(\epsnuc\,\epsnucT\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
epsilon_rho	\(\epsnuc\,\epsnucrho\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
Omega_rot	\(\Orot\)	dataset	real (n)	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.00

Data items in the root HDF5 group of version-1.00 GSM-format files are as follows:

Item	Symbol	Object type	Data type	Definition
n	\(N\)	attribute	integer	number of grid points
version	—	attribute	integer	100
M_star	\(M\)	attribute	real	stellar mass [\(\gram\)]
R_star	\(R\)	attribute	real	photospheric radius [\(\cm\)]
L_star	\(L\)	attribute	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
r	\(r\)	dataset	real (n)	radial coordinate [\(\cm\)]
M_r	\(M_r\)	dataset	real (n)	interior mass [\(\gram\)]
L_r	\(L_{r}\)	dataset	real (n)	interior luminosity [\(\erg\,\second^{-1}\)]
P	\(P\)	dataset	real (n)	total pressure [\(\barye\)]
rho	\(\rho\)	dataset	real (n)	density [\(\gram\,\cm^{-3}\)]
T	\(T\)	dataset	real (n)	temperature [\(\kelvin\)]
N2	\(N^{2}\)	dataset	real (n)	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
Gamma_1	\(\Gamma_{1}\)	dataset	real (n)	adiabatic exponent
nabla_ad	\(\nabla_{\rm ad}\)	dataset	real (n)	adiabatic temperature gradient
delta	\(\upsT\)	dataset	real (n)	thermodynamic coefficient
nabla	\(\nabla\)	dataset	real (n)	temperature gradient
kap	\(\kappa\)	dataset	real (n)	opacity [\(\cm^{2}\,\gram^{-1}\)]
kap_T	\(\kappa\,\kapT\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
kap_rho	\(\kappa\,\kaprho\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
eps	\(\epsilon\)	dataset	real (n)	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_T	\(\epsnuc\,\epsnucT\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_rho	\(\epsnuc\,\epsnucrho\)	dataset	real (n)	nuclear energy generation rate partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
Omega_rot	\(\Orot\)	dataset	real (n)	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.10

Data items in the root HDF5 group of version-1.10 GSM-format files are as follows:

Item	Symbol	Object type	Data type	Definition
n	\(n\)	attribute	integer	number of grid points
version	—	attribute	integer	110
M_star	\(M\)	attribute	real	stellar mass [\(\gram\)]
R_star	\(R\)	attribute	real	photospheric radius [\(\cm\)]
L_star	\(L\)	attribute	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
r	\(r\)	dataset	real (n)	radial coordinate [\(\cm\)]
M_r	\(M_r\)	dataset	real (n)	interior mass [\(\gram\)]
L_r	\(L_{r}\)	dataset	real (n)	interior luminosity [\(\erg\,\second^{-1}\)]
P	\(P\)	dataset	real (n)	total pressure [\(\barye\)]
rho	\(\rho\)	dataset	real (n)	density [\(\gram\,\cm^{-3}\)]
T	\(T\)	dataset	real (n)	temperature [\(\kelvin\)]
N2	\(N^{2}\)	dataset	real (n)	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
Gamma_1	\(\Gamma_{1}\)	dataset	real (n)	adiabatic exponent
nabla_ad	\(\nabla_{\rm ad}\)	dataset	real (n)	adiabatic temperature gradient
delta	\(\upsT\)	dataset	real (n)	thermodynamic coefficient
nabla	\(\nabla\)	dataset	real (n)	dimensionless temperature gradient
kap	\(\kappa\)	dataset	real (n)	opacity [\(\cm^{2}\,\gram^{-1}\)]
kap_kap_T	\(\kappa\,\kapT\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
kap_kap_rho	\(\kappa\,\kaprho\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
eps	\(\epsilon\)	dataset	real (n)	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_eps_T	\(\epsnuc\,\epsnucT\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_eps_rho	\(\epsnuc\,\epsnucrho\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
Omega_rot	\(\Orot\)	dataset	real (n)	rotation angular frequency [\(\radian\,\second^{-1}\)]

Version 1.20

Data items in the root HDF5 group of version-1.20 GSM-format files are as follows:

Item	Symbol	Object type	Data type	Definition
n	\(n\)	attribute	integer	number of grid points
version	—	attribute	integer	120
M_star	\(M\)	attribute	real	stellar mass [\(\gram\)]
R_star	\(R\)	attribute	real	photospheric radius [\(\cm\)]
L_star	\(L\)	attribute	real	photospheric luminosity [\(\erg\,\second^{-1}\)]
r	\(r\)	dataset	real (n)	radial coordinate [\(\cm\)]
M_r	\(M_r\)	dataset	real (n)	interior mass [\(\gram\)]
L_r	\(L_{r}\)	dataset	real (n)	interior luminosity [\(\erg\,\second^{-1}\)]
P	\(P\)	dataset	real (n)	total pressure [\(\barye\)]
rho	\(\rho\)	dataset	real (n)	density [\(\gram\,\cm^{-3}\)]
T	\(T\)	dataset	real (n)	temperature [\(\kelvin\)]
N2	\(N^{2}\)	dataset	real (n)	Brunt-Väisälä frequency squared [\(\second^{-2}\)]
Gamma_1	\(\Gamma_{1}\)	dataset	real (n)	adiabatic exponent
nabla_ad	\(\nabla_{\rm ad}\)	dataset	real (n)	adiabatic temperature gradient
delta	\(\upsT\)	dataset	real (n)	thermodynamic coefficient
nabla	\(\nabla\)	dataset	real (n)	dimensionless temperature gradient
kap	\(\kappa\)	dataset	real (n)	opacity [\(\cm^{2}\,\gram^{-1}\)]
kap_kap_T	\(\kappa\,\kapT\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
kap_kap_rho	\(\kappa\,\kaprho\)	dataset	real (n)	opacity partial [\(\cm^{2}\,\gram^{-1}\)]
eps	\(\epsilon\)	dataset	real (n)	total energy generation rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_eps_T	\(\epsnuc\,\epsnucT\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_eps_rho	\(\epsnuc\,\epsnucrho\)	dataset	real (n)	nuclear energy generation partial [\(\erg\,\second^{-1}\,\gram^{-1}\)]
eps_grav	\(\epsgrav\)	dataset	real (n)	gravothermal energy release rate [\(\erg\,\second^{-1}\,\gram^{-1}\)]
Omega_rot	\(\Orot\)	dataset	real (n)	rotation angular frequency [\(\radian\,\second^{-1}\)]

POLY File Format

Files in POLY format store HDF5 data describing a composite polytrope model. This format adheres to the following conventions:

• All data objects are attached to the root HDF5 group (/)

• Real values are written with type H5T_IEEE_F64LE when GYRE is compiled in double precision (the default), and type H5T_IEEE_F32LE otherwise

• Integer values are written with type H5T_STD_I32LE

Data items in the root HDF5 group are as follows:

Data Item	Variable	Object type	Data type	Definition
n	\(N\)	attribute	integer	number of grid points
n_r	\(\nreg\)	attribute	integer	number of regions
n_poly	\(n_{i}\)	attribute	real (n_r)	polytropic indices of regions
z_b	\(z_{i-1/2}\)	attribute	real (n_r-1)	radial coordinates of region boundaries
Delta_b	\(\Delta_{i-1/2}\)	attribute	real (n_r-1)	log density jump of region boundaries
Gamma_1	\(\Gamma_{1}\)	attribute	real	first adiabatic exponent
z	\(z\)	dataset	real (n)	polytropic radial coordinate
theta	\(\theta\)	dataset	real (n)	Lane-Emden variable
dtheta	\(\theta'\)	dataset	real (n)	Derivative of Lane-Emden variable

Oscillation Equations

This chapter outlines how the oscillation equations solved by the GYRE frontends are obtained from the basic equations of stellar structure.

Fluid Equations

The starting point is the fluid equations, comprising the conservation laws for mass

\[\pderiv{\rho}{t} + \cdot \nabla \left( \rho \vv \right) = 0\]

and momentum

\[\rho \left( \pderiv{}{t} + \vv \cdot \nabla \right) \vv = -\nabla P - \rho \nabla \Phi;\]

the heat equation

\[\rho T \left( \pderiv{}{t} + \vv \cdot \nabla \right) S = \rho \epsnuc - \nabla \cdot (\vFrad + \vFcon);\]

and Poisson’s equation

\[\nabla^{2} \Phi = 4 \pi G \rho.\]

Here, \(\rho\), \(P\), \(T\), \(S\) and \(\vv\) are the fluid density, pressure, temperature, specific entropy and velocity; \(\Phi\) is the self-gravitational potential; \(\epsnuc\) is the specific nuclear energy generation rate; and \(\vFrad\) and \(\vFcon\) are the radiative and convective energy fluxes. An explicit expression for the radiative flux is provided by the radiative diffusion equation,

\[\vFrad = - \frac{c}{3\kappa\rho} \nabla (a T^{4}),\]

where \(\kappa\) is the opacity and \(a\) the radiation constant.

The fluid equations are augmented by the thermodynamic relationships between the four state variables (\(P\), \(T\), \(\rho\) and \(S\)). Only two of these are required to uniquely specify the state (we assume that the composition remains fixed over an oscillation cycle). In GYRE, \(P\) and \(S\) are adopted as these primary variables[1], and the other two are presumed to be derivable from them:

\[\rho = \rho(P, S), \qquad T = T(P, S).\]

The nuclear energy generation rate and opacity are likewise presumed to be functions of the pressure and entropy:

\[\epsnuc = \epsnuc(P, S), \qquad \kappa = \kappa(P, S).\]

Footnotes

[1]

This may seem like a strange choice, but it simplifies the switch between adiabatic and non-adiabatic calculations

Equilibrium State

In a static equilibrium state the fluid velocity \(\vv\) vanishes. The momentum equation then becomes the hydrostatic equilibrium equation

\[\nabla P = - \rho \nabla \Phi.\]

Assuming the equilibrium is spherically symmetric, this simplifies to

\[\deriv{P}{r} = - \rho \deriv{\Phi}{r}.\]

Poisson’s equation can be integrated once to yield

\[\deriv{\Phi}{r} = \frac{G}{r^{2}} \int 4 \pi \rho r^{2} \, \diff{r} = \frac{G M_{r}}{r^{2}},\]

where the second equality introduces the interior mass

()\[M_{r} \equiv \int 4 \pi \rho r^{2} \, \diff{r}.\]

The hydrostatic equilibrium equation thus becomes

\[\deriv{P}{r} = - \rho \frac{G M_{r}}{r^{2}}.\]

The heat equation in the equilibrium state is

\[\rho T \pderiv{S}{t} = \rho \epsnuc - \nabla \cdot (\vFrad + \vFcon).\]

If the star is in thermal equilibrium then the left-hand side vanishes, and the nuclear heating rate balances the flux divergence term. Again assuming spherical symmetry, this is written

\[\deriv{}{r} \left( \Lrad + \Lcon \right) = 4 \pi r^{2} \rho \epsnuc,\]

where

\[\Lrad \equiv 4 \pi r^{2} \Fradr, \qquad \Lcon \equiv 4 \pi r^{2} \Fconr\]

are the radiative and convective luminosities, respectively.

Linearized Equations

Applying an Eulerian (fixed position, denoted by a prime) perturbation to the mass and momentum conservation equations, they linearize about the static equilibrium state as

\[\rho' + \nabla \cdot ( \rho \vv' ) = 0,\]

()\[\rho \pderiv{\vv'}{t} = - \nabla P' - \rho' \nabla \Phi - \rho \nabla \Phi'.\]

(in these expressions, the absence of a prime denotes an equilibrium quantity). Likewise, Poisson’s equation becomes

\[\nabla^{2} \Phi' = 4 \pi G \rho'\]

Applying a Lagrangian (fixed mass element, denoted by a \(\delta\)) perturbation to the heat equation, it linearizes about the equilibrium state as

\[T \pderiv{\delta S}{t} = \delta \epsnuc - \delta \left( \frac{1}{\rho} \nabla \cdot \vFrad \right),\]

where the heating term \(\delta (\rho^{-1} \nabla \cdot \vFcon)\) has been dropped[1] due to the continued lack of a workable theory for pulsation-convection coupling. Likewise applying a Lagrangian perturbation to the radiative diffusion equation,

\[\delta \vFrad = \left( 4 \frac{\delta T}{T} - \frac{\delta \rho}{\rho} - \frac{\delta \kappa}{\kappa} \right) \vFrad + \frac{\delta(\nabla \ln T)}{\sderiv{\ln T}{r}} \Fradr.\]

The thermodynamic relations linearize to

\[\frac{\delta \rho}{\rho} = \frac{1}{\Gammi} \frac{\delta P}{P} - \upsT \frac{\delta S}{\cP}, \qquad \frac{\delta T}{T} = \nabad \frac{\delta P}{P} + \frac{\delta S}{\cP},\]

and the perturbations to the nuclear energy generation rate and opacity can be expressed as

\[\frac{\delta \epsnuc}{\epsnuc} = \epsnucad \frac{\delta P}{P} + \epsnucS \frac{\delta S}{\cP}, \qquad \frac{\delta \kappa}{\kappa} = \kapad \frac{\delta P}{P} + \kapS \frac{\delta S}{\cP}.\]

In these expressions, Eulerian and Lagrangian perturbations to any scalar quantity \(f\) are related via

\[\frac{\delta f}{f} = \frac{f'}{f} + \frac{\xir}{r} \deriv{\ln f}{\ln r}.\]

Moreover, the thermodynamic partial derivatives are defined as

\[\Gammi = \left( \pderiv{\ln P}{\ln \rho} \right)_{S}, \quad \upsT = \left( \pderiv{\ln \rho}{\ln T} \right)_{P}, \quad \cP = \left( \pderiv{S}{\ln T} \right)_{P}, \quad \nabad = \left( \pderiv{\ln T}{\ln P} \right)_{S},\]

and the nuclear and opacity partials are

\[\epsnucad = \left( \pderiv{\ln \epsnuc}{\ln P} \right)_{\rm ad}, \quad \epsnucS = \cP \left( \pderiv{\ln \epsnuc}{S} \right)_{P}, \quad \kapad = \left( \pderiv{\ln \kappa}{\ln P} \right)_{\rm ad}, \quad \kapS = \cP \left( \pderiv{\ln \kappa}{S} \right)_{P}.\]

The latter can be calculated from corresponding density and temperature partials via

\[\begin{split}\begin{gathered} \kapad = \frac{\kaprho}{\Gammi} + \nabad \kapT, \qquad \kapS = -\upsT \kaprho + \kapT, \\ \epsnucad = \frac{\epsnucrho}{\Gammi} + \nabad \epsnucT, \qquad \epsnucS = -\upsT \epsnucrho + \epsnucT. \end{gathered}\end{split}\]

Footnotes

[1]

This is known as a frozen convection approximation. GYRE offers multiple ways to freeze convection; see the Convection Effects section for further details.

Separated Equations

With a separation of variables in spherical-polar coordinates \((r,\theta,\phi)\), and assuming an oscillatory time (\(t\)) dependence with angular frequency \(\sigma\), solutions to the linearized equations can be expressed as

()\[\begin{split}\begin{aligned} \xir(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txir(r) \, Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii \sigma t) \right], \\ \xit(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txih(r) \, \pderiv{}{\theta} Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii \sigma t) \right], \\ \xip(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txih(r) \, \frac{\ii m}{\sin\theta} Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii \sigma t) \right], \\ f'(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \tf'(r) \, Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii \sigma t) \right]. \end{aligned}\end{split}\]

Here, \(\xir\), \(\xit\) and \(\xip\) are the radial, polar and azimuthal components of the displacement perturbation vector \(\vxi\); \(Y^{m}_{\ell}\) is the spherical harmonic with harmonic degree \(\ell\) and azimuthal order \(m\); and again \(f\) stands for any perturbable scalar. The displacement perturbation vector is related to the velocity perturbation via

\[\vv' = \pderiv{\vxi}{t}.\]

Substituting the above solution forms into the linearized equations, the mechanical (mass and momentum conservation) equations become

\[\trho' + \frac{1}{r^{2}} \deriv{}{r} \left( \rho r^{2} \txir \right) - \frac{\ell(\ell+1)}{r} \rho \txih = 0,\]

\[-\sigma^{2} \rho \txir = - \deriv{\tP'}{r} - \trho' \deriv{\Phi}{r} - \rho \deriv{\tPhi'}{r},\]

\[-\sigma^{2} \rho r \txih = - \tP' - \rho \tPhi'.\]

Likewise, Poisson’s equation becomes

\[\frac{1}{r^{2}} \deriv{}{r} \left( r^{2} \deriv{\tPhi'}{r} \right) - \frac{\ell(\ell+1)}{r^{2}} \txih = 4 \pi G \trho'\]

and the heat equation becomes

\[-\ii \sigma T \delta \tS = \delta \tepsnuc - \deriv{\delta \tLrad}{M_{r}} + \frac{\ell(\ell+1)}{\sderiv{\ln T}{r}} \frac{\Fradr}{\rho} \frac{\tT'}{T} + \ell(\ell + 1) \frac{\txih}{r} \deriv{\Lrad}{M_{r}},\]

where

\[\delta \tLrad \equiv 4 \pi r^{2} \left( \delta \tFradr + 2 \frac{\txir}{r} \Fradr \right)\]

is the Lagrangian perturbation to the radiative luminosity. The radial part of the radiative diffusion equation becomes

\[\delta\tFradr = \left[ 4 \frac{\delta \tT}{T} - \frac{\delta\trho}{\rho} - \frac{\delta\tkappa}{\kappa} + \frac{\sderiv{(\delta \tT/T)}{\ln r}}{\sderiv{\ln T}{\ln r}} \right] \Fradr.\]

Finally, the thermodynamic, nuclear and opacity relations become

\[\frac{\delta \trho}{\rho} = \frac{1}{\Gammi} \frac{\delta \tP}{P} - \upsT \frac{\delta \tS}{\cP}, \qquad \frac{\delta \tT}{T} = \nabla_{\rm ad} \frac{\delta \tP}{P} + \frac{\delta \tS}{\cP},\]

\[\frac{\delta \tepsnuc}{\epsnuc} = \epsnucad \frac{\delta \tP}{P} + \epsnucS \frac{\delta \tS}{\cP}, \qquad \frac{\delta \tkappa}{\kappa} = \kapad \frac{\delta \tP}{P} + \kapS \frac{\delta \tS}{\cP}.\]

Boundary Conditions

To form a closed system, the separated equations are augmented by algebraic relations at the inner and outer boundaries of the computational domain, and possibly at interior points as well.

Inner Boundary

When the inner boundary is placed at the stellar origin (\(r=0\)), the requirement that solutions remain finite leads to the set of regularity conditions

\[\begin{split}\begin{aligned} \txir - \ell \txih = 0, \\ \ell \tPhi' - r \deriv{\tPhi'}{r} = 0, \\ \delta \tS = 0. \end{aligned}\end{split}\]

Sometimes it’s desirable that the inner boundary is instead placed at \(r > 0\) — for instance, to excise the stellar core from the oscillation calculations. Then, there is much more flexibility in the choice of inner boundary condition. Possible options include setting \(\txir = 0\) or \(\txih=0\) instead of the first equation above.

Outer Boundary

The outer boundary typically corresponds to the stellar surface. Under the assumption that the density vanishes on and above this surface, the gravitational potential must match onto a solution to Laplace’s equation that remains finite at infinity, leading to the potential boundary condition

\[(\ell + 1) \tPhi' + r \deriv{\tPhi'}{r} = 0.\]

Likewise, the assumption that there is no external pressure acting on the star (consistent with the vanishing surface density) gives the momentum boundary condition

\[\delta \tP = 0.\]

Finally, the thermal boundary condition can be derived from the equation

\[T^{4}(\tau) = \frac{4\sigma}{3} \Fradr \left( \tau + \frac{2}{3} \right)\]

describing the thermal structure of an atmosphere under the combined plane-parallel, grey, Eddington, local thermodynamic equilibrium and radiative equilibrium approximations. Here, \(\tau\) is the vertical optical depth and \(\sigsb\) the Stefan-Boltzmann constant. Setting \(\tau=0\) (again, consistent with the vanishing surface density) and perturbing this equation yields the desired boundary condition

\[4 \frac{\delta \tT}{T} = \frac{\delta \tFradr}{\Fradr}.\]

Complications arise when realistic stellar models are considered, because these typically extend only out to an optical depth \(\tau=2/3\) (or some similar value) corresponding to the photosphere. In such cases the density does not vanish at the nominal stellar surface, and the outer boundary conditions must be modified to account for the effects of the atmosphere region lying above the surface. Many stellar oscillation codes, including GYRE, can adopt more sophisticated formulations for the momentum boundary condition — for instance, based on the assumption that the outer atmosphere has an isothermal stratification. However, the atmospheric effects on the potential and thermal boundary conditions are usually neglected.

Internal Boundaries

Internal boundaries arise when the density and other related quantities show a radial discontinuity within the star. Across such a discontinuity \(\txir\), \(\delta \tP\) and \(\delta \tFradr\) must remain continuous[1]. Internal boundary conditions on other perturbations follow from integrating the separated equations across the discontinuity, resulting in

\[\begin{split}\begin{aligned} \tP^{\prime +} - \tP^{\prime -} &= \deriv{\Phi}{r} \left( \rho^{+} - \rho^{-} \right) \txir, \\ \left. \deriv{\tPhi'}{r} \right|^{+} - \left. \deriv{\tPhi'}{r} \right|^{-} &= - 4 \pi G \left( \rho^{+} - \rho^{-} \right) \txir, \\ \tT^{\prime +} - \tT^{\prime -} &= 0. \end{aligned}\end{split}\]

Here, + (-) superscripts indicate quantities evaluated on the inner (outer) side of the discontinuity.

Footnotes

[1]

This is to ensure that the fluid doesn’t ‘tear’, and that pressure forces and radiative heating remain finite.

Dimensionless Formulation

To improve numerical stability, GYRE solves the separated equations and boundary conditions by recasting them into a dimensionless form that traces its roots back to Dziembowski (1971).

Variables

The independent variable is the fractional radius \(x \equiv r/R\) (with \(R\) the stellar radius), and the dependent variables \(\{y_{1},y_{2},\ldots,y_{6}\}\) are

()\[\begin{split}\begin{aligned} y_{1} &= x^{2 - \ell}\, \frac{\txir}{r}, \\ y_{2} &= x^{2-\ell}\, \frac{\tP'}{\rho g r}, \\ y_{3} &= x^{2-\ell}\, \frac{\tPhi'}{gr}, \\ y_{4} &= x^{2-\ell}\, \frac{1}{g} \deriv{\tPhi'}{r}, \\ y_{5} &= x^{2-\ell}\, \frac{\delta \tS}{c_{p}}, \\ y_{6} &= x^{-1-\ell}\, \frac{\delta \tLrad}{L}. \end{aligned}\end{split}\]

Oscillation Equations

The dimensionless oscillation equations are

()\[\begin{split}\begin{aligned} x \deriv{y_{1}}{x} &= \left( \frac{V}{\Gammi} - 1 - \ell \right) y_{1} + \left( \frac{\ell(\ell+1)}{c_{1} \omega^{2}} - \alphagam \frac{V}{\Gammi} \right) y_{2} + \alphagrv \frac{\ell(\ell+1)}{c_{1} \omega^{2}} y_{3} + \upsT \, y_{5}, \\ % x \deriv{y_{2}}{x} &= \left( c_{1} \omega^{2} - \fpigam \As \right) y_{1} + \left( 3 - U + \As - \ell \right) y_{2} - \alphagrv y_{4} + \upsT \, y_{5}, \\ % x \deriv{y_{3}}{x} &= \alphagrv \left( 3 - U - \ell \right) y_{3} + \alphagrv y_{4} \\ % x \deriv{y_{4}}{x} &= \alphagrv \As U y_{1} + \alphagrv \frac{V}{\Gammi} U y_{2} + \alphagrv \ell(\ell+1) y_{3} - \alphagrv (U + \ell - 2) y_{4} - \alphagrv \upsT \, U y_{5}, \\ % x \deriv{y_{5}}{x} &= \frac{V}{\frht} \left[ \nabad (U - c_{1}\omega^{2}) - 4 (\nabad - \nabla) + \ckapad V \nabla + \cdif \right] y_{1} + \mbox{} \\ & \frac{V}{\frht} \left[ \frac{\ell(\ell+1)}{c_{1} \omega^{2}} (\nabad - \nabla) - \ckapad V \nabla - \cdif \right] y_{2} + \mbox{} \\ & \alphagrv \frac{V}{\frht} \left[ \frac{\ell(\ell+1)}{c_{1} \omega^{2}} (\nabad - \nabla) \right] y_{3} + \alphagrv \frac{V \nabad}{\frht} y_{4} + \mbox{} \\ & \left[ \frac{V \nabla}{\frht} (4 \frht - \ckapS) + \dfrht + 2 - \ell \right] y_{5} - \frac{V \nabla}{\frht \crad} y_{6} \\ % x \deriv{y_{6}}{x} &= \left[ \alphahfl \ell(\ell+1) \left( \frac{\nabad}{\nabla} - 1 \right) \crad - V \cepsad - \alphaegv c_{egv} \nabad V \right] y_{1} + \mbox{} \\ & \left[ V \cepsad - \ell(\ell+1) \crad \left( \alphahfl \frac{\nabad}{\nabla} - \frac{3 + \dcrad}{c_{1}\omega^{2}} \right) + \alphaegv c_{egv} \nabad V \right] y_{2} + \mbox{} \\ & \alphagrv \left[ \ell(\ell+1) \crad \frac{3 + \dcrad}{c_{1}\omega^{2}} \right] y_{3} + \mbox{} \\ & \left[ \cepsS - \alphahfl \frac{\ell(\ell+1)\crad}{\nabla V} + \ii \alphathm \omega \cthk + \alphaegv c_{egv} \right] y_{5} - \left[ 1 + \ell \right] y_{6}, \end{aligned}\end{split}\]

where the dimensionless oscillation frequency is introduced as

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

(with \(M\) the stellar mass). These differential equations are derived from the separated equations, with the insertion of ‘switch’ terms (denoted \(\alpha\)) that allow certain pieces of physics to be altered. See the Physics Switches section for more details.

For non-radial adiabatic calculations, the last two equations above are set aside and the \(y_{5}\) terms dropped from the first four equations. For radial adiabatic calculations with reduce_order=.TRUE. (see the Oscillation Parameters section), the last four equations are set aside and the first two replaced by

\[\begin{split}\begin{aligned} x \deriv{y_{1}}{x} &= \left( \frac{V}{\Gammi} - 1 \right) y_{1} - \frac{V}{\Gamma_{1}} y_{2}, \\ % x \deriv{y_{2}}{x} &= \left( c_{1} \omega^{2} + U - \As \right) y_{1} + \left( 3 - U + \As \right) y_{2}. \end{aligned}\end{split}\]

Boundary Conditions

Inner Boundary

When inner_bound='REGULAR', GYRE applies regularity-enforcing conditions at the inner boundary:

\[\begin{split}\begin{aligned} c_{1} \omega^{2} y_{1} - \ell y_{2} - \alphagrv \ell y_{3} &= 0, \\ \alphagrv \ell y_{3} - (2\alphagrv - 1) y_{4} &= 0, \\ y_{5} &= 0. \end{aligned}\end{split}\]

(these are the dimensionless equivalents to the expressions appearing in the Boundary Conditions section).

When inner_bound='ZERO_R', the first and second conditions above are replaced with zero radial displacement conditions,

\[\begin{split}\begin{aligned} y_{1} &= 0, \\ y_{4} &= 0. \end{aligned}\end{split}\]

Likewise, when inner_bound='ZERO_H', the first and second conditions are replaced with zero horizontal displacement conditions,

\[\begin{split}\begin{aligned} y_{2} - y_{3} &= 0, \\ y_{4} &= 0. \end{aligned}\end{split}\]

Outer Boundary

When outer_bound='VACUUM', GYRE applies the outer boundary conditions

\[\begin{split}\begin{aligned} y_{1} - y_{2} &= 0 \\ \alphagrv U y_{1} + (\alphagrv \ell + 1) y_{3} + \alphagrv y_{4} &= 0 \\ (2 - 4\nabad V) y_{1} + 4 \nabad V y_{2} + 4 \frht y_{5} - y_{6} &= 0 \end{aligned}\end{split}\]

(these are the dimensionless equivalents to the expressions appearing in the Boundary Conditions section).

When outer_bound='DZIEM', the first condition above is replaced by the Dziembowski (1971) outer boundary condition,

\[\left\{ 1 + V^{-1} \left[ \frac{\ell(\ell+1)}{c_{1} \omega^{2}} - 4 - c_{1} \omega^{2} \right] \right\} y_{1} - y_{2} = 0.\]

When outer_bound='UNNO' or 'JCD', the first condition is replaced by the (possibly-leaky) outer boundary conditions described by Unno et al. (1989) and Christensen-Dalsgaard (2008), respectively. When outer_bound='ISOTHERMAL', the first condition is replaced by a (possibly-leaky) outer boundary condition derived from a local dispersion analysis of waves in an isothermal atmosphere.

Finally, when outer_bound='GAMMA', the first condition is replaced by the outer momentum boundary condition described by Ong & Basu (2020).

Internal Boundaries

Across density discontinuities, GYRE applies the boundary conditions

\[\begin{split}\begin{aligned} U^{+} y_{2}^{+} - U^{-} y_{2}^{-} &= y_{1} (U^{+} - U^{-}) \\ y_{4}^{+} - y_{4}^{-} &= -y_{1} (U^{+} - U^{-}) \\ y_{5}^{+} - y_{5}^{-} &= - V^{+} \nabad^{+} (y_{2}^{+} - y_{1}) + V^{-} \nabad^{-} (y_{2}^{-} - y_{1}) \end{aligned}\end{split}\]

(these are the dimensionless equivalents to the expressions appearing in the Boundary Conditions section). Here, + (-) superscripts indicate quantities evaluated on the inner (outer) side of the discontinuity. \(y_{1}\), \(y_{3}\) and \(y_{6}\) remain continuous across discontinuities, and therefore don’t need these superscripts.

Structure Coefficients

The various stellar structure coefficients appearing in the dimensionless oscillation equations and boundary conditions are defined as follows:

\[\begin{split}\begin{gathered} V = -\deriv{\ln P}{\ln r} \qquad V_{2} = x^{-2} V \qquad \As = \frac{1}{\Gamma_{1}} \deriv{\ln P}{\ln r} - \deriv{\ln \rho}{\ln r} \qquad U = \deriv{\ln M_{r}}{\ln r} \\ % c_1 = \frac{r^{3}}{R^{3}} \frac{M}{M_{r}} \qquad \fpigam = \begin{cases} \alphapi & \As > 0, x < x_{\rm atm} \\ \alphagam & \As > 0, x > x_{\rm atm} \\ 1 & \text{otherwise} \end{cases}\\ % \nabla = \deriv{\ln T}{\ln P} \qquad \clum = x^{-3} \frac{\Lrad+\Lcon}{L} \qquad \crad = x^{-3} \frac{\Lrad}{L} \qquad \dcrad = \deriv{\ln \crad}{\ln r} \\ % \frht = 1 - \alpharht \frac{\ii \omega \cthn}{4} \qquad \dfrht = - \alpharht \frac{\ii \omega \cthn \dcthn}{4 \frht} \\ % \ckapad = \frac{\alphakar \kaprho}{\Gamma_{1}} + \nabad \alphakat \kapT \qquad \ckapS = - \upsT \alphakar \kaprho + \alphakat \kapT \\ % \ceps = x^{-3} \frac{4\pi r^{3} \rho \epsnuc}{L} \qquad \cepsad = \ceps \epsnucad \qquad \cepsS = \ceps \epsnucS \\ % \cdif = - 4 \nabad V \nabla + \nabad \left(V + \deriv{\ln \nabad}{\ln x} \right) \\ % \cthn = \frac{\cP}{a c \kappa T^{3}} \sqrt{\frac{GM}{R^{3}}} \qquad \dcthn = \deriv{\ln \cthn}{\ln r} \\ % \cthk = x^{-3} \frac{4\pi r^{3} \cP T \rho}{L} \sqrt{\frac{GM}{R^{3}}} \cegv = x^{-3} \frac{4\pi r^{3} \rho \epsgrav}{L} \end{gathered}\end{split}\]

Physics Switches

GYRE offers the capability to adjust the oscillation equations through a number of physics switches, controlled by parameters in the &osc namelist group (see the Oscillation Parameters section). The table below summarizes the mapping between the switches appearing in the expressions above, and the corresponding namelist parameters.

Symbol	Parameter	Description
\(\alphagrv\)	alpha_grv	Scaling factor for gravitational potential perturbations. Set to 1 for normal behavior, and to 0 for the Cowling (1941) approximation
\(\alphathm\)	alpha_thm	Scaling factor for local thermal timescale. Set to 1 for normal behavior, to 0 for the non-adiabatic reversible (NAR) approximation (see Gautschy et al., 1990), and to a large value to approach the adiabatic limit
\(\alphahfl\)	alpha_hfl	Scaling factor for horizontal flux perturbations. Set to 1 for normal behavior, and to 0 for the non-adiabatic radial flux (NARF) approximation (see Townsend, 2003b)
\(\alphagam\)	alpha_gam	Scaling factor for g-mode isolation. Set to 1 for normal behavior, and to 0 to isolate g modes as described by Ong & Basu (2020)
\(\alphapi\)	alpha_pi	Scaling factor for p-mode isolation. Set to 1 for normal behavior, and to 0 to isolate p modes as described by Ong & Basu (2020)
\(\alphakar\)	alpha_kar	Scaling factor for opacity density partial derivative. Set to 1 for normal behavior, and to 0 to suppress the density part of the \(\kappa\) mechanism
\(\alphakat\)	alpha_kat	Scaling factor for opacity temperature partial derivative. Set to 1 for normal behavior, and to 0 to suppress the temperature part of the \(\kappa\) mechanism
\(\alpharht\)	alpha_rht	Scaling factor for time-dependent term in the radiative heat equation (see Unno & Spiegel, 1966). Set to 1 to include this term (Unno calls this the Eddington approximation), and to 0 to ignore the term
\(\alphatrb\)	alpha_trb	Scaling factor for the turbulent mixing length. Set to the convective mixing length to include the turbulent damping term (see the Convection Effects section), and to 0 to ignore the term

Rotation Effects

The oscillation equations presented in the preceding sections are formulated for a non-rotating star. The corresponding equations for a rotating star are significantly more complicated, and a complete treatment of rotation lies beyond the scope of GYRE. However, GYRE can include two important effects arising from rotation.

Doppler Shift

A lowest-order effect of rotation arises in the Doppler shift from transforming between the inertial reference frame and the local co-rotating reference frame. To incorporate this effect in the separated equations, all instances of the inertial-frame frequency \(\sigma\) are replaced by the co-rotating frequency

()\[\sigmac \equiv \sigma - m \Orot,\]

where \(m\) is the azimuthal order of the mode and \(\Orot\) is the rotation angular frequency. GYRE assumes shellular rotation (see, e.g., Meynet & Maeder, 1997), and so the latter can in principle be a function of radial coordinate \(r\). The corresponding modifications to the dimensionless formulation involve replacing the dimensionless inertial-frame frequency \(\omega\) with the dimensionless co-rotating frequency

\[\omegac \equiv \omega - m \Orot \sqrt{\frac{R^{3}}{GM}}.\]

Perturbative Coriolis Force Treatment

Another lowest-order effect of rotation arises from the Coriolis force. For slow rotation, this effect can be determined through a perturbation expansion technique (see, e.g., section 19.2 of Unno et al., 1989). To first order in \(\Orot\), the frequency of a mode is shifted by the amount

\[\Delta \sigma = m \int_{0}^{R} \Orot \, \deriv{\beta}{r} \diff{r},\]

where the rotation splitting kernel is

\[\deriv{\beta}{r} = \frac{\left\{ \txir^{2} + [\ell(\ell+1) - 1] \txih^{2} - 2 \txir \txih \right\} \rho r^{2}} {\int_{0}^{R} \left\{ \txir^{2} + \ell(\ell+1) \txih^{2} \right\} \rho r^{2} \diff{r}}\]

In this latter expression, the eigenfunctions \(\txir\) and \(\txih\) are evaluated from solutions to the oscillation equations without rotation. Therefore, the expression above for \(\Delta \sigma\) can be applied as a post-calculation correction to non-rotating eigenfrequencies.

Non-Perturbative Coriolis Force Treatment

The perturbation expansion technique above breaks down when \(\Orot/\sigmac \gtrsim 1\). To deal with such cases, the gyre frontend [1] can incorporate a non-perturbative treatment of the Coriolis force based on the ‘traditional approximation of rotation’ (TAR). The TAR was first introduced by Eckart (1960; Hydrodynamics of Oceans and Atmospheres) and has since been used extensively within the pulsation community (see, e.g., Bildsten et al., 1996; Lee & Saio, 1997; Townsend, 2003a; Bouabid et al., 2013; Townsend, 2020).

Within the TAR, the solution forms given in eqn. (7) are replaced by

()\[\begin{split}\begin{aligned} \xir(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txir(r) \, \houghr(\theta) \, \exp(\ii m \phi -\ii \sigma t) \right], \\ \xit(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txih(r) \, \frac{\hought(\theta)}{\sin\theta} \, \exp(\ii m \phi -\ii \sigma t) \right], \\ \xip(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \txih(r) \, \frac{\houghp(\theta)}{\ii \sin\theta} \, \exp(\ii m \phi -\ii \sigma t) \right], \\ f'(r,\theta,\phi;t) &= \operatorname{Re} \left[ \sqrt{4\pi} \, \tf'(r) \, \houghr(\theta) \, \exp(\ii m \phi -\ii \sigma t) \right] \end{aligned}\end{split}\]

Here, the Hough functions \(\houghr\), \(\hought\) and \(\houghp\) are the eigenfunctions obtained by solving Laplace’s tidal equations (TEs), a second-order system of differential equations and boundary conditions in the polar (\(\theta\)) coordinate (see Townsend 2020). Together with the associated eigenvalue \(\lambda\), they depend on the harmonic degree \(\ell\)[2] and azimuthal order \(m\), and the spin parameter

\[q \equiv \frac{2 \Orot}{\sigmac}.\]

Solution Families

Solutions to the TEs can be grouped into two families based on the behavior of the eigenfunctions and eigenvalue in the limit \(\Orot \rightarrow 0\). For the gravito-acoustic family,

()\[\begin{split}\left. \begin{aligned} \houghr(\theta) \ \rightarrow & \ Y^{m}_{\ell}(\theta,0) \\ \hought(\theta) \ \rightarrow & \ \sin\theta \pderiv{}{\theta} Y^{m}_{\ell}(\theta,0) \\ \houghp(\theta) \ \rightarrow & \ - m Y^{m}_{\ell}(\theta,0) \end{aligned} \right\} \quad \text{as } \Orot \rightarrow 0.\end{split}\]

and \(\lambda \rightarrow \ell(\ell+1)\). With these expressions, the solution forms (12) reduce to those given in eqn. (7).

Conversely, for the Rossby family

()\[\begin{split}\left. \begin{aligned} \houghr(\theta) \ \rightarrow & \ 0 \\ \hought(\theta) \ \rightarrow & \ m Y^{m}_{\ell}(\theta,0) \\ \houghp(\theta) \ \rightarrow & \ - \sin\theta \pderiv{}{\theta} Y^{m}_{\ell}(\theta,0) \end{aligned} \right\} \quad \text{as } \Orot \rightarrow 0.\end{split}\]

and \(\lambda \rightarrow 0\). Moreover, Rossby-mode eigenfrequencies also show the limiting behavior

()\[\sigmac = \frac{2 m \Orot}{\ell(\ell+1)} \quad \text{as } \Orot \rightarrow 0,\]

which is independent of the stellar structure.

Implementing the TAR

To implement the TAR in the separated equations and boundary conditions, all instances of the term \(\ell(\ell+1)\) are replaced by the TE eigenvalue \(\lambda\). Then, all instances of the harmonic degree \(\ell\) are replaced by \(\elle\), an effective harmonic degree found by solving

\[\elle(\elle+1) = \lambda.\]

Similar steps are taken in the dimensionless formulation, but in the definitions of the dependent variables \(\{y_{1},y_{2},\ldots,y_{6}\}\), \(\ell\) is replaced by \(\elli\), the effective harmonic degree evaluated at the inner boundary.

Footnotes

[1]

Currently the TAR cannot be used with the gyre_tides frontend, because it doesn’t play well with forcing by the tidal potential \(\PhiT\).

[2]

The harmonic degree isn’t formally a ‘good’ quantum number in the TAR; however, it can still be used to identify Hough functions by considering their behavior in the limit \(\Orot \rightarrow 0\), as given in eqns. (13) and (14).

Convection Effects

The oscillation equations presented in the preceding sections neglect the thermal and mechanical effects of convection. GYRE provides functionality for controlling how the thermal effects are suppressed, and how the mechanical effects can be included in a limited way.

Frozen Convection

In the derivation of the linearized equations, a term \(\delta (\rho^{-1} \nabla \cdot \vFcon)\) is dropped from the perturbed heat equation. This is known as a frozen convection approximation, and is grounded in the assumption that the energy transport by convection remains unaffected affected by the pulsation. There’s more than one way to freeze convection; Pesnell (1990) presents a systematic review of different approaches. GYRE currently implements a subset of these:

• Pesnell’s case 1, neglecting \(\delta (\rho^{-1} \nabla \cdot \vFcon)\) in the perturbed heat equation.

• Pesnell’s case 4, neglecting \(\delta \Lcon\) (the Lagrangian perturbation to the convective luminosity) in the perturbed heat equation.

For further details, see the conv_scheme parameter in the Oscillation Parameters section.

Turbulent Damping

The Reynolds number in stars is very large, and thus convection tends to be turbulent. Following the treatment by Savonije & Witte (2002), GYRE can partially incorporate the mechanical effects of this turbulence by adding a term

\[f_{r,{\rm visc}} = - \frac{1}{r^{2}} \pderiv{}{r} \left( \rho \nu r^{2} \pderiv{v'_{r}}{r} \right)\]

to the radial component of the linearized momentum equation (6), representing the viscous force arising from radial fluid motion. Because this term depends on \(v'_{r}\), it is phase-shifted by a quarter cycle relative to the other terms in the equation, and acts like a drag force that damps oscillations. The turbulent viscous coefficient \(\nu\) is evaluated as

\[\nu = \frac{(\alphatrb H_{P})^{2}}{\tconv} \left[ 1 + \tconv \frac{\sigma}{2\pi} \right]^{-1},\]

where \(H_{P}\) is the pressure scale height, \(\alphatrb\) is the turbulent mixing length (in units of \(H_{P}\)), and \(\tconv\) the convective turnover timescale. This expression is adapted from equation (18) of Savonije & Witte (2002), with an exponent \(s=1\).

In GYRE \(\alphatrb\) is implemented as a switch (see the Physics Switches section). A reasonable choice is to set this parameter equal to the MLT mixing length parameter \(\alpha_{\rm MLT}\) of the stellar model.

Tidal Effects

To simulate the effects of tidal forcing by a companion, the gyre_tides frontend solves a modified form of the linearized momentum equation (6), namely

\[\rho \pderiv{\vv'}{t} = - \nabla P' - \rho' \nabla P - \rho \nabla \Phi' - \rho \nabla \PhiT.\]

The final term on the right-hand side represents the external force arising from the tidal gravitational potential \(\PhiT\).

Tidal Potential

The tidal potential can be expressed via the superposition

()\[\PhiT = \sum_{\ell=2}^{\infty} \sum_{m=-\ell}^{\ell} \sum_{k=-\infty}^{\infty} \PhiTlmk.\]

of partial tidal potentials defined by

\[\PhiTlmk \equiv - \epsT \, \frac{GM}{R} \, \cbar_{\ell,m,k} \left( \frac{r}{R} \right)^{\ell} Y^{m}_{\ell}(\theta, \phi) \, \exp(- \ii k \Oorb t).\]

(the summation over \(\ell\) and \(m\) comes from a multipolar space expansion of the potential, and the summation over \(k\) from a Fourier time expansion). Here,

\[\epsT = \left( \frac{R}{a} \right)^{3} = \frac{\Oorb R^{3}}{GM} \frac{q}{1+q}\]

quantifies the overall strength of the tidal forcing, in terms of the companion’s mass \(q M\), semi-major axis \(a\) and orbital angular frequency \(\Oorb\). These expressions, and the definition of the tidal expansion coefficients \(\cbar_{\ell,m,k}\), are presented in greater detail in Sun et al. (2023).

Separated Equations

Because the tidal potential (16) superposes many different spherical harmonics, the solution forms (7) must be replaced by the more-general expressions

()\[\begin{split}\begin{aligned} \xir(r,\theta,\phi;t) &= \sum_{\ell,m,k} \txirlmk(r) \, Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii k \Oorb t), \\ \xit(r,\theta,\phi;t) &= \sum_{\ell,m,k} \txihlmk(r) \, \pderiv{}{\theta} Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii k \Oorb t), \\ \xip(r,\theta,\phi;t) &= \sum_{\ell,m,k} \txihlmk(r) \, \frac{\ii m}{\sin\theta} Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii k \Oorb t), \\ f'(r,\theta,\phi;t) &= \sum_{\ell,m,k} \tflmk'(r) \, Y^{m}_{\ell}(\theta,\phi) \, \exp(-\ii k \Oorb t) \end{aligned}\end{split}\]

(the notation for the sums has been abbreviated). Substituting these solution forms into the linearized equations, and taking advantage of the orthonormality of the spherical harmonics, leads to a fully separated set of differential equations for each combination of \(\ell\), \(m\) and \(k\). A given set resembles the regular separated equations, with just a couple changes:

• The perturbation \(\tPhi'\) is replaced by \(\tPsi' \equiv \tPhi' + \tPhiT\), representing the total (self + tidal) gravitational potential perturbation.

• Rather than being an eigenvalue parameter, the oscillation frequency is set by \(\sigma = k \Oorb\), representing the forcing frequency of the partial tidal potential in an inertial frame.

The latter change means that the dimensionless frequency (10) becomes

\[\omega = \alphafrq \, k \Oorb \sqrt{\frac{R^{3}}{GM}},\]

where \(\alphafrq\) is an additional term introduced to allow tuning of the tidal forcing frequency (see the alpha_frq parameter in the Tidal Parameters section).

Boundary Conditions

The boundary conditions accompanying the separated equations for a given \(\{\ell,m,k\}\) combination resemble those presented previously, except that the outer potential boundary condition becomes

\[(\ell + 1) \tPsi' + r \deriv{\tPsi'}{r} = (2 \ell + 1) \tPhiTlmk,\]

where

()\[\tPhiTlmk \equiv - \epsT \, \frac{GM}{R} \, \cbar_{\ell,m,k} \left( \frac{r}{R} \right)^{\ell}.\]

describes the radial dependence of the partial tidal potential.

Composite Polytropes

Composite polytropes are an extension of standard polytrope models for stellar structure, to allow for discontinuities in the density \(\rho\) and/or the polytropic index \(n\). This appendix lays out the mathematical formalism[1] underpinning them; it is intended to complement the Building POLY Models appendix, which describes how composite polytrope models can be built using the build_poly executable.

Equation of State

Consider a composite polytrope composed of \(\nreg\) regions extending from the origin out to the stellar surface. In the \(i\)’th region (\(1 \leq i \leq \nreg\)), the pressure \(P\) and density \(\rho\) are related by the polytropic equation-of-state

()\[\frac{P_{i}}{P_{i,0}} = \left( \frac{\rho_{i}}{\rho_{i,0}} \right)^{(n_{i} + 1)/n_{i}}\]

where the normalizing pressure \(P_{i,0}\) and density \(\rho_{i,0}\), together with the polytropic index \(n_{i}\), are constant across the region but may change from one region to the next. At the \(\nreg-1\) boundaries between adjacent regions, the pressure and interior mass \(M_{r}\) are required to be continuous, but the density may jump.

Structure Equations

Lane-Emden Equation

In the \(i\)’th region, a composite polytrope satisfies the equation of hydrostatic equilibrium

\[-\frac{1}{\rho_{i}} \deriv{P_{i}}{r} = \deriv{\Phi_{i}}{r}\]

Substituting in the polytropic equation-of-state 19 yields

\[\frac{(n_{i}+1) P_{i,0}}{\rho_{i,0}^{1+1/n_{i}}} \deriv{}{r} \left( \rho_{i}^{1/n_{i}} \right) = - \deriv{\Phi_{i}}{r},\]

which can then be integrated with respect to \(r\) to give

\[\frac{(n_{i}+1)P_{i,0}}{\Phi_{i,0} \, \rho_{i,0}} \left( \frac{\rho_{i}^{1/n_{i}}}{\rho_{i,0}^{1/n_{i}}} - 1 \right) = \left( 1 - \frac{\Phi_{i}}{\Phi_{i,0}} \right).\]

Here, the constants of integration have been chosen so that \(\Phi_{i} = \Phi_{i,0}\) when \(\rho_{i} = \rho_{i,0}\). Rearranging, the density follows as

\[\rho_{i} = \rho_{i,0} \, \theta_{i}^{n_{i}},\]

where the polytropic dependent variable is introduced as

\[\theta_{i} = \left[ \frac{\Phi_{i,0} \, \rho_{i,0}}{(n_{i} + 1) P_{i,0}} \left( 1 - \frac{\Phi_{i}}{\Phi_{i,0}} \right) + 1 \right].\]

With these expressions, Poisson’s equation

\[\frac{1}{r^{2}} \deriv{}{r} \left( r^{2} \deriv{P_{i}}{r} \right) = 4 \pi G \rho_{i}\]

is recast as

\[\frac{1}{r^{2}} \deriv{}{r} \left( r^{2} \deriv{\theta_{i}}{r} \right) = - \frac{1}{A_{i}} \theta_{i}^{n_{i}},\]

where

\[A_{i} \equiv \frac{(n_{i} + 1) P_{i,0}}{4 \pi G \rho_{i,0}^{2}}.\]

A change of variables to the polytropic independent variable \(z \equiv A_{1}^{-1/2} r\) results in the dimensionless form

()\[\frac{1}{z^{2}} \deriv{}{z} \left( z^{2} \deriv{\theta_{i}}{z} \right) = - B_{i} \theta_{i}^{n_{i}},\]

where \(B_{i} \equiv A_{1}/A_{i}\). This can be regarded as a generalization of the usual Lane-Emden equation to composite polytropes.

Continuity Relations

At the boundary between adjacent regions, the pressure and interior mass must be continuous. If \(z_{i-1/2}\) denotes the coordinate of the boundary between the \(i-1\) and \(i\) regions, then these continuity relations are expressed as

\[\begin{split}\left. \begin{gathered} B_{i} = \frac{n_{i-1} + 1}{n_{i} + 1} \frac{\theta_{i}^{n_{i}+1}}{\theta_{i-1}^{n_{i-1}+1}} \frac{\rho_{i,0}^{2}}{\rho_{i-1,0}^{2}} \, B_{i-1}, \\ \theta'_{i} = \frac{n_{i-1} + 1}{n_{i} + 1} \frac{\theta_{i-1}^{n_{i-1}+1}}{\theta_{i}^{n_{i}+1}} \frac{\rho_{i,0}}{\rho_{i-1,0}} \, \theta'_{i-1}, \end{gathered} \right\} \quad \text{at} \ z = z_{i-1/2}\end{split}\]

respectively.

Solution Method

Specification

The structure of a composite polytrope is specified completely by

• a set of \(\nreg\) polytropic indices \(n_{i}\)

• a set of \(\nreg-1\) boundary coordinates \(z_{i-1/2}\)

• a set of \(\nreg\) density jumps \(\Delta_{i-1/2} \equiv \ln [\rho_{i}(z_{i-1/2})/\rho_{i-1}(z_{i-1/2}]\)

Although the normalizing densities \(\rho_{i,0}\) have so far been left unspecified, it’s convenient to choose them as the density at the beginning of their respective regions.

Solution

The structure equations may be solved as an initial value problem. In the first region (\(i=1\)) this IVP involves integrating the Lane-Emden equation 20 from the center \(z=0\) to the first boundary \(z=z_{3/2}\), with the initial conditions

\[\begin{split}\left. \begin{gathered} \theta_{i} = 1, \\ \theta'_{i} = 0, \\ B_{1} = 1, \\ t_{1} = 1 \end{gathered} \right\} \quad \text{at}\ z=0\end{split}\]

(here, \(t_{i} \equiv \rho_{i,0}/\rho_{1,0}\)).

The IVP in the intermediate regions (\(i = 2,\ldots,\nreg-1\)) involves integrating from \(z=z_{i-1/2}\) to \(z=z_{i+1/2}\), with initial conditions established from the preceding region via

\[\begin{split}\left. \begin{gathered} \theta_{i} = 1, \\ \theta'_{i} = \frac{n_{i-1} + 1}{n_{i} + 1} \frac{\theta_{i-1}^{n_{i-1}+1}}{\theta_{i}^{n_{i}+1}} \frac{t_{i}}{t_{i-1}} \, \theta'_{i-1}, \\ B_{i} = \frac{n_{i-1} + 1}{n_{i} + 1} \frac{\theta_{i}^{n_{i}+1}}{\theta_{i-1}^{n_{i-1}+1}} \frac{t_{i}^{2}}{t_{i-1}^{2}} \, B_{i-1}, \\ \ln t_{i} = \ln t_{i-1} + n_{i-1} \ln \theta_{i-1} - n_{i} \ln \theta_{i} + \Delta_{i-1/2}. \end{gathered} \right\} \quad \text{at}\ z=z_{i-1/2}\end{split}\]

The IVP in the final region (\(i=\nreg\)) involves integrating from \(z_{\nreg-1/2}\) until \(\theta_{\nreg} = 0\). This point defines the stellar surface, \(z=z_{\rm s}\). For some choices of \(n_{i}\), \(z_{i-1/2}\) and/or \(\Delta_{i-1/2}\), the point \(\theta=0\) can arise in an earlier region \(i = \nreg_{\rm t} < \nreg\); in such cases, the model specification must be truncated to \(\nreg_{\rm t}\) regions.

Physical Variables

Once the Lane-Emden equation 20 has been solved, the density in each region can be evaluated by

\[\rho_{i} = \rho_{1,0} \, t_{i} \, \theta_{i}^{n_{i}}.\]

The pressure then follows from the equation-of-state 19 as

\[P_{i} = P_{1,0} \, \frac{n_{1}+1}{n_{i}+1} \, \frac{t_{i}^{2}}{B_{i}} \, \theta_{i}^{n_{i}+1}.\]

The interior mass \(m\) is evaluated by introducing the auxiliary quantity \(\mu\), which is defined in the first region by

\[\mu_{1}(z) = - z^{2} \theta'_{1} (z),\]

and in subsequent regions by

\[\mu_{i}(z) = \mu_{i-1}(z_{i-1/2}) - \frac{t_{i}}{B_{i}} \left[ z^{2} \theta'_{i} (z) - z_{i-1/2}^{2} \theta'_{i} (z_{i-1/2}) \right].\]

The interior mass then follows as

\[M_{r} = M \frac{\mu_{i}}{\mu_{\rm s}},\]

where \(\mu_{\rm s} \equiv \mu_{\nreg}(z_{\rm s})\).

Structure Coefficients

The structure coefficients for composite polytropic models are evaluated using

\[\begin{split}\begin{gathered} V_{2} = -(n_{i} + 1) \frac{z_{\rm s}^{2}}{z} \frac{\theta'_{i}}{\theta_{i}}, \qquad \As = V_{2} \frac{z^{2}}{z_{\rm s}^{2}}\left( \frac{n_{i}}{n_{i} + 1} - \frac{1}{\Gamma_{1}} \right), \\ U = \frac{t_{i} z^{3}}{\mu_{i}} \theta_{i}^{n_{i}}, \qquad c_1 = \frac{z^{3}}{z_{\rm s}^{3}} \frac{\mu_{\rm s}}{\mu_{i}}, \end{gathered}\end{split}\]

where \(\mu\) is the auxiliary mass variable introduced above.

Footnotes

[1]

The formalism here is based on Mixed polytrope with density discontinuities (Christensen-Dalsgaard, 2015, unpublished), with extensions to allow for constant-density regions.

Building POLY Models

This appendix describes the build_poly executable, which builds a composite polytropic stellar model and writes it to a file in the POLY format.

Installation

The build_poly executable is automatically compiled when GYRE is built, and installed in the $GYRE_DIR/bin directory (see the main Installation chapter).

Example Walkthrough: Simple Polytrope

As the first example of build_poly in action, let’s build a simple (i.e., single-region) \(n=3\) polytrope, that for instance describes the structure of a radiation-pressure dominated, fully convective star.

Assembling a Namelist File

First, let’s assemble a namelist file containing the various parameters which control a build_poly run. Using a text editor, create the file build_poly.simple.in with the following content cut-and-pasted in:

&poly
	n_poly = 3.0  ! Polytropic index of single region
/

&num
	dz = 1E-2     ! Radial spacing of points
	toler = 1E-10 ! Tolerance of integrator
/

&out
	file = 'poly.simple.h5' ! Name of output file
/

Detailed information on the namelist groups expected in build_poly input files can be found in the Input Files section. Here, let’s briefly narrate the parameters appearing in the file above:

• In the &poly namelist group, the n_poly parameter sets the polytropic index.

• In the &num namelist group, the dz parameter sets the radial spacing of points, and the toler parameter sets the tolerance of the numerical integrator.

• In the &output namelist group, the file parameter sets the name of the output file.

Running build_poly

To run build_poly, use the command

$GYRE_DIR/bin/build_poly build_poly.simple.in

There is no screen output produced during the run, but at the end the poly.simple.h5 will be written to disk. This file, which is in POLY format, can be used as the input stellar model in a GYRE calculation; but it can also be explored in Python (see Fig. 12) using the read_model function from PyGYRE.

Plot of the Lane-Emden solution variable \(\theta\), density \(\rho\), pressure \(P\) and interior mass \(M_{r}\) as a function of radial coordinate, for the simple polytrope. (Source)

Example Walkthrough: Composite Polytrope

As the second example of build_poly in action, let’s build a two-region composite polytrope. The polytropic index is \(n=3\) in the inner region, and \(n=1.5\) in the outer region. At the boundary between the regions, located at radial coordinate \(z=1.4\), the logarithmic density jump is \(\Delta = -0.5\).

Assembling a Namelist File

Using a text editor, create the file build_poly.composite.in with the following content cut-and-pasted in:

&poly
	n_r = 2           ! Number of regions
	n_poly = 3.0, 1.5 ! Polytropic indices of regions
        z_b = 1.4         ! Radial coordinate of region boundaries
        Delta_b = -0.5    ! Logarithmic density jump at region boundaries
/

&num
	dz = 1E-2     ! Radial spacing of points
	toler = 1E-10 ! Tolerance of integrator
/

&out
	file = 'poly.composite.h5' ! Name of output file
/

Again, detailed information on the namelist groups expected in build_poly input files can be found in the Input Files section. Here, let’s briefly narrate the parameters appearing in the file above:

• In the &poly namelist group, the n_r parameter sets the number of regions; the n_poly parameter sets the polytropic indices in the two regions; the z_b sets the radial coordinate of the boundary between the regions; and the Delta_b sets the density jump at this boundary.

• In the &num namelist group, the dz parameter sets the radial spacing of points, and the toler parameter sets the tolerance of the numerical integrator.

• In the &output namelist group, the file parameter sets the name of the output file.

Running build_poly

As before, to run build_poly use the command

$GYRE_DIR/bin/build_poly build_poly.composite.in

There is no screen output produced during the run, but at the end the poly.composite.h5 will be written to disk. This file, which is in POLY format, can be used as the input stellar model in a GYRE calculation; but it can also be explored in Python (see Fig. 13) using the read_model function from PyGYRE.

Plot of the Lane-Emden solution variable \(\theta\), density \(\rho\), pressure \(P\) and interior mass \(M_{r}\) as a function of radial coordinate, for the composite polytrope. Note the density discontinuity, and the associated discontinuities in the gradients of the pressure and interior mass. (Source)

Input Files

The build_poly executable reads parameters from an input file that defines a number of Fortran namelist groups, as described below.

Polytrope Parameters

The &poly namelist group defines polytrope parameters; the input file can contain only one. Allowable parameters are:

n_r (default 1)

Number of regions

n_poly (default 0)

Comma-separated list of length n_r, specifying polytropic indices for regions

z_b

Comma-separated list of length n_r-1, specifying radial coordinates of boundaries between regions

Delta_b

Comma-separated list of length n_r-1, specifying logarithmic density jumps at boundaries between regions

Gamma_1 (default 5./3.)

First adiabatic exponent

Numerical Parameters

The &num namelist group defines numerical parameters; the input file can contain only one. Allowable parameters are:

dz (default 1E-2)

Spacing of grid points in polytropic radial coordinate \(z\)

toler (default 1E-10)

Relative and absolute tolerance of Lane-Emden integrations

Output Parameters

The &out namelist group defines output parameters; the input file can contain only one. Allowable parameters are:

file

Name of POLY-format file to write to

Evaluating Tidal Eigenvalues

This appendix describes the eval_lambda executable, which evaluates the eigenvalue \(\lambda\) appearing in Laplace’s tidal equations (see the Rotation Effects section). This executable is used for the calculations presented in Townsend (2020).

Installation

eval_lambda is automatically compiled when GYRE is built, and installed in the $GYRE_DIR/bin directory (see the main Installation chapter).

Running

Unlike most other GYRE executables, the parameters for eval_lambda are supplied directly on the command line, with the syntax

./eval_lambda l m q_min q_max n_q log_q rossby filename

This evaluates \(\lambda\) for harmonic degree \(\ell\) and azimuthal order \(m\) on a grid \(\{q_{1},q_{2},\ldots,q_{N}\}\) in the spin parameter, writing the results to the file filename. If the flag log_q has the value T then the grid is logarithmically spaced:

\[q_{i} = 10^{(1 - w_{i}) \log q_{\rm min} + w_{i} \log q_{\rm max}},\]

where

\[w_{i} \equiv \frac{i-1}{N-1}.\]

Alternatively, if log_q has the value F, then the grid is linearly spaced:

\[q_{i} = (1 - w_{i}) q_{\rm min} + w_{i} q_{\rm max}.\]

As a special case, when \(n_{q}=1\), \(q_{\rm min}\) and \(q_{\rm max}\) must match, and the single \(q\) point is set to equal them.

If the flag rossby has the value T, then the Rossby \(\lambda\) family is evaluated; otherwise, the gravito-acoustic family is evaluated.

The table below summarizes the mapping between the user-definable controls appearing in the expressions above, and the corresponding command-line parameters:

Symbol	Parameter
\(\ell\)	l
\(m\)	m
\(q_{\rm min}\)	q_min
\(q_{\rm max}\)	q_max
\(N\)	n_q

Interpreting Output

The output file created by eval_lambda is in GYRE’s HDF Format, with the following data:

l (integer scalar)

Harmonic degree \(\ell\)

k (integer scalar)

Meridional order \(k\) (see Townsend, 2003a)

m (integer scalar)

Azimuthal order \(m\)

rossby (logical scalar)

Rossby family flag

q (real array)

Spin parameter \(q\)

lambda (real array)

Eigenvalue \(\lambda\)