Lab 2: Calculating Mode Frequencies

Learning Goals

• Learn how to set up and run Gyre
• Visualize oscillations in an échelle diagram and study their eigenfunctions.

If you need a MESA Model

If for some reason you did not complete lab 1, you can download a rotating model LOGS directory here. This model has an initial rotation rate of $\Omega /{\Omega }_{crit}$ = 1d-4 and has output profiles at ${\nu }_{\max }$ = 100.

Gyre namelist

Just like with MESA, GYRE takes an input file that we will call an “inlist”. Again, just like MESA, these are plaintext files which are use the namelist format/functionality of Fortran. The GYRE inlist file is typically named something like gyre.in. For this lab we can and will start with an empty folder.

&model
/
 
&constants
/
 
&mode
/
 
&osc
/
 
&num
/
 
&scan
/
 
&grid
/
 
&rot
/
 
&ad_output
/
 
&nad_output
/

As you can see, the inlist has sections which function similarly to how the inlist sections in MESA inlists work. However, in GYRE’s case, there are many more types of inputs that we can use. References for each of these sections can be found here in the GYRE docs.

So we need to fill out this inlist with the correct information. The following parts of this section will describe each namelist group starting from the top.

Model

Here’s where we tell GYRE what type of model to read in and calculate frequencies for. The docs have a full explanation of all the options.

We will add the following lines here:

&model
    model_type = 'EVOL'
    file='PATH_TO_PROFILE/profile##.data.GYRE'
    file_format = 'MESA'
/

As described by the docs, this will tell GYRE that we are using an external evolutionary model, and that file is called profile##.data.GYRE and it’s using MESA’s default GYRE file format.

Constants

Here we can define various constants found here. It is good practice to use identical values as used in the code generating your model. This group can be left empty for this lab as we are going to use the default values.

Mode

This namegroup defines which modes you want to calculate. You can state the angular degree ($\ell$ ) and the azimuthal order ($m$ ). For each type of mode we will need one extra &mode namegroup. For now we will leave the azimuthal order at its default value ($m=0$ ); it will come up later when we include rotation. We will also assign the modes with a tag to differentiate between them in other namegroups. For this lab using 'radial' for the $\ell =0$ modes and 'non-radial' for all others is enough.

&mode
    l=0
    tag = 'radial'
/

&mode
    l=1
    tag = 'non-radial'
/

&mode
    l=2
    tag = 'non-radial'
/

Oscillation parameters

In this namegroup we can set various options on how to calculate the oscillations themselves. Which physics and assumptions to use, different boundary conditions, how to scale various parameters and so on. All of the options and their default values can be found here.

&osc
    outer_bound = 'JCD'
    variables_set = 'JCD'
/

Numerical Parameters

Here goes everything to do with numerical parameters. Again, we leave nearly everything at default (see here) and only increase the difference scheme from a second-order to a fourth-order Gauss-Legendre collocation.

&num
    diff_scheme = 'COLLOC_GL4'
/

Frequency Scan Parameters

This section tells GYRE in which frequency range it will scan for the oscillation eigenfunctions. The modes that are of interest to us are those closest to ${\nu }_{\max }$ . We suggest to set freq_min and freq_max to ${\nu }_{\max }\pm 3\Delta \nu$ . Look into your history.data file to find the values of ${\nu }_{\max }$ and $\Delta \nu$ for your model with the profile you choose in the “Model” section. (You wrote this down in our google sheet from Lab1, if you need a refresher). We also set the units of our chosen range to μHz (freq_units = 'UHZ')

Next, we define the resolution, the number of points of our scan. Generally, we need the grid spacing to be smaller than the eigenfrequency separation of adjacent modes, across the full range of the grid. There is no absolute rule as how to determine this value beforehand, but a too low value could miss some modes and a high value means a longer runtime. For our radial p-modes n_freq = 200 is enough but because we calculate the oscillations also for red giants we do not only have p-modes but also mixed-modes (for the non-radial modes). These will need a higher value. To be more efficient we therefore define two different &scan namegroups. To define which modes correspond to which scan we set tag_list = 'radial' (the tag we defined in &mode):

&scan
  tag_list = 'radial'
  grid_type = 'LINEAR'  ! Scan grid uniform in frequency
  freq_min =  ###       ! Minimum frequency to scan from
  freq_max =  ###       ! Maximum frequency to scan to  
  n_freq = 200          ! Number of frequency points in scan
  freq_units = 'UHZ'
/

For the non-radial modes g-modes become relevant. These are not equally spaced in frequency like the p-modes but in period. The higher the star evolves up the RGB (${\nu }_{\max }$ decreases) the more tightly spaced they get. Therefore it is recommended to scan the grid not uniform in frequency (grid_type = 'LINEAR') but uniform in period (grid_type = 'INVERSE').

&scan
  tag_list = 'non-radial'
  grid_type = 'INVERSE' ! Scan grid uniform in inverse frequency
  freq_min = ###        ! Minimum frequency to scan from
  freq_max = ###        ! Maximum frequency to scan to
  n_freq = 3000         ! Number of frequency points in scan
  freq_units = 'UHZ'
/

Grid Parameters

In this group you can modify the spatial resolution of the grid generated by Gyre to calculate the eigenmodes. If it is too low the eigenmodes are not resolved and the accuracy of the results suffers. The grid spacing should be smaller than the scale of the smallest significant variation of the eigenfunction. Here you can see how the following weighting parameters are used to define and refine the grid. The defaults are set to 0 but for a appropriate resolution you should choose higher values.

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

Output

Now it’s time to tell gyre what parameters it should save. Similar to profiles and history files we also have two different output types in GYRE. The summary file is as the name implies a summary of all oscillation modes GYRE found. You can include all parameters that describe the mode with a single value e.g. $l$ , $m$ , $n$ , frequency, inertia, … . The name of the file is given by the summary_file and the parameters it should include are given with summary_item_list. You can rename this if you want to run more than one model (e.g. summary_numax500.txt)(Or put it in the same folder as the details as described in the next paragraph). All the options can be found here. For this lab we will also change the file format to a textfile so it becomes human readable instead of the default HDF5 file.

&ad_output
  summary_file = 'summary.txt'                         
  summary_item_list = 'l,m,n_pg,n_p,n_g,freq,E_norm,E_p,E_g'
  summary_file_format = 'TXT'  
...

The detail file on the other hand gives you a detailed description of one mode. We have several files for one run and therefore we give a template of the name instead to avoid overwriting the file for each mode we calculate. E.g. '%l' will be replaced with the angular degree of the corresponding mode and '%n' with its radial order. All options can be found here. A good idea is to make a new folder for all your detail files if you want to run several models and provide the path to it (for example: detail_template = 'details_numax500/detail.l%l.n%n.h5). Be aware that GYRE does not generate any folders for you and you will get an error if it can’t find the folder. Keep an eye out on your storage, one file may not be large but it sums up quite quickly when you have a lot of modes (That’s why we will stick to the HDF5 format). Additionally we will define the frequency units again with freq_units.

...
  detail_template = 'detail.l%l.n%n.h5'               
  detail_item_list = 'l,n_pg,omega,x,xi_r,xi_h,c_1,As,
                      V_2,Gamma_1,rho,P,R_star,M_star' 
  !detail_file_format = 'TXT'  

  freq_units = 'UHZ'                      
/

&nad_output
/

We put everything in the namegroup &ad_output. It tells gyre that we will assume adiabatic conditions. If we want to instead calculate it non-adiabatically we would put it in &nad_output instead.

Rotation

For now we will ignore rotation and leave &rot empty. We will come back to it in Lab3.

But if we calculate rotation we need to assume some approximations. In Lab 3 we will set up GYRE to calculate the effects of the Coriolis force based on the traditional approximation of rotation (TAR). This assumption allows us to neglect the radial component of the Coriolis effects. But this also means we need to be able to separate them into a horizontal and radial component. This is not quite possible if we consider gravitational potential perturbations. We need to use the Cowling Approximation, which sets the gravitational potential perturbation to zero.

alpha_grv = 0 
reduce_order = .false.

Be aware that this is only done to better compare the rotating to non-rotating results. For other science cases it is not necessary to use these approximations.

Running GYRE

$GYRE_DIR/bin/gyre gyre.in

GYRE will output a bunch of numbers describing the grids it generates and the modes it found. When it finishes you can open the summary file and inspect its contents.

Tables are nice but we also want to plot these outputs. Gyre doesn’t come with something like pgstar for MESA and we have to make the plots ourself. We prepared a google colab link for you with various plotting functions. You can open it here (local copy) . You may need to save it into your google drive folder to be able to edit it.

Upload your history.data, profile#.data, profiles.index and summary.txt files (There is a folder symbol on the left under the key symbol. Click on it and drag and drop your files into it). Run the first 2 cells by clicking on the arrow to the right of the cell. It will install mesa_reader and pygyre, which allow you to easily read MESA and Gyre output files. Next, write your profile number in the third cell and run this and all following cells.

The first plot shows a HRD with a red cross marking the model corresponding to your chosen profile.

The second plot shows an échelle diagram, which is a useful tool to display solar-like oscillations. In a power spectrum, the different oscillation modes form a repeating pattern. P-modes with the same angular degree ($\ell$ ) but a different radial order ($n$ ) are separated approximately by the large frequency separation ($\Delta \nu$ ). In the échelle diagram we plot the frequency over frequency mod $\Delta \nu$ , we basically split the power spectrum into parts with length $\Delta \nu$ and stack them on top of each other. This way modes with the same $\ell$ are sorted on top of each other. This pattern gets a little bit more “messy” when we move up the RGB and mixed-modes become relevant, when p-modes couple with g-modes. The g-modes are not arranged like the p-modes. Instead of being equally spaced in frequency, they are equally spaced in period. Thus we can see several modes (for $\ell =1$ and $2$ not for $\ell =0$ . there are no g-modes for $\ell =0$ ) in the same ‘row’ of the echelle diagram. The markersize of the modes corresponds to the p-mode inertia. The bigger the marker the higher the amplitude of the mode on the surface of the star and the easier to observe it.

The last cell called “plot eigenfunctions” will output the names of 2 detail files. Upload them and rerun the cell.

In the new figure we chose two $\ell =1$ modes close to ${\nu }_{\max }$ to look at them in more detail. These two modes are marked with squares in the echelle diagram to the right. One of them has a low p-mode inertia and one a high one. We also added a propagation diagram at the top left to get a better idea of what is happening inside the star and where the modes are propagating. Remember, we made a version of it in Lab1 with pgstar already. It tells you in which areas of the star (x-axis is the radius coordinate) the modes can propagate for which frequencies (y-axis is frequency squared). In the orange section g-modes can propagate and in the blue area p-modes. In the white area (the “evanescent zone”) both types of modes get damped. If this evanescent zone is small enough (like it is in red giants), the p-modes on the surface and the g-modes in the core can couple to each other, forming “mixed modes”. We added a dot-dashed line at ${\nu }_{\max }$ to see at which frequency the modes get excited. We plot the radial and horizontal displacement functions for both modes. As the p-modes are acoustic modes (longitudinal waves) the radial displacement is more relevant for them. For g-modes (transverse waves) the buoyancy is the restoring force and thus the horizontal displacement is more important. Can you identify which plot shows the p-mode dominated and which the g-mode dominated mixed mode? How can you tell?

We generated a small movie to show how the oscillations change over the evolution of the star.

Full Solutions

Here you can find a full solution of Lab2 including the history and profile files needed for plotting.