Lab 3: Asteroseismology

Authors: Niall Miller (lead TA), Eliza Frankel, Joey Mombarg - Lecturer: Yaguang Li — MESA Summer School 2026, Tetons, Wyoming

In Labs 1 and 2 we used rotation periods and CMD positions to constrain stellar ages. In this lab we add a third technique: asteroseismology. We will compute two seismic observables – the large frequency separation $\Delta \nu$ and the small frequency separation $\delta {\nu }_{02}$ – and ask whether they constrain ages better than what we measured in Lab 2. We will also build a map of which stars are actually accessible to real space missions.

All of the code that computes the seismic quantities is already implemented in run_star_extras.f90. You do not need to modify it. Your job is to configure the run, interpret the outputs, and combine results across the group.

The large frequency separation $\Delta \nu$

The large frequency separation is the average spacing between p-modes of the same angular degree $\ell$ and consecutive radial order $n$ . It is related to the sound crossing time of the star:

Delta_nu_int = 0.
do k = 2, s% nz, 1
   dr = s% rmid(k-1) - s% rmid(k)
   Delta_nu_int = Delta_nu_int + dr/s% csound(k)
end do
Delta_nu_int = 1./(2.*Delta_nu_int)

The loop accumulates $dr/c_s$ , building up the total sound travel time. The final line takes the reciprocal and divides by 2 to give $\Delta \nu$ in Hz. Because $\Delta \nu \propto \sqrt{M/R^3}$ , it is a measure of the mean stellar density – it decreases as the star expands during main-sequence evolution.

The small frequency separation $\delta {\nu }_{02}$

The small frequency separation is the offset between $\ell =0$ and $\ell =2$ modes of similar frequency. Unlike $\Delta \nu$ , it is sensitive to the gradient of the sound speed in the stellar core, and therefore directly to the central hydrogen abundance:

delta_nu02_int = 0.
nu_max = s% nu_max * 1d-6
do k = 2, s% nz, 1
   dc = s% csound(k-1) - s% csound(k)
   delta_nu02_int = delta_nu02_int + dc / s% rmid(k)
end do
delta_nu02_int = delta_nu02_int - s% csound(1)/s% r(1)
delta_nu02_int = -2 * Delta_nu_int * delta_nu02_int / (3.1415926535**2. * nu_max)

The loop computes $d{c}_s/r$ shell by shell, approximating the radial derivative of the sound speed as a finite difference. The surface boundary term $c_s(R)/R$ is subtracted after the loop. The whole integral is then scaled by $-2\Delta \nu /({\pi }^2{\nu }_{\max })$ to give $\delta {\nu }_{02}$ in Hz. Because the sound speed gradient in the core steepens as hydrogen is depleted, $\delta {\nu }_{02}$ decreases monotonically with age - making it a direct age clock. (https://ui.adsabs.harvard.edu/abs/2005MNRAS.356..671O/abstract)

Step 1 – Setup

Get the working directory. Download the Lab 3 starter and unzip it:

→ Download the Lab 3 starter (Google Drive)

Everything you need is already inside – including the run_star_extras.f90 with the seismic calculations and the ../data/ files – so Lab 3 runs on its own. Work inside this directory; there is nothing to bring over from Labs 1 or 2.

For this lab, all of your configuration lives in a single file: inlist_run. Every namelist – &star_job, &controls, &colors, and &pgstar – sits together in that one file.

This is a change from Lab 2. There, a small top-level inlist acted as a header that pointed MESA off to separate files: the physics in inlist_project and the plotting settings in inlist_pgstar. Here there is no header and no splitting up by namelist – you point the star binary at one self-contained file and that is the whole configuration.

Have a look around the directory before running anything:

cd ~/Downloads              #or wherever the zip file is

unzip Lab3_start.zip        #unzip it

cd Lab3_start/Lab3          #Change directories to it.

ls                          #to see whats here

cat rn                      #lets see what rn does (it copies 'inlist_run' to 'inlist' then calls rn1)

cat rn1                     #it calls star with 'inlist'
                            #so we have established that 'inlist_run' is what we modify. It will be copied to inlist for execution.

run_star_extras.f90 already implements the seismic calculations, so you will not need to edit it.

We are going to crowdsource our science as a group. Open the shared Google Sheet, add your name to column A, and pick a mass from 0.5, 0.6, 0.8, 0.9, 1.0, 1.1, 1.2 $M_{\odot }$ . You can take more than one mass, but try not to duplicate masses that others have already claimed.

Once you have chosen a mass, open inlist_run and set it:

initial_mass = X.X   ! set to your assigned value

you’ll see it set to a default of 0.4 — replace that number with your assigned mass.

Step 2 – Configure the inlist

Look through the below pgstar display. It configures pgstar to show five panels that update in real time:

Before you run anything, you need to set up two things: the pgstar live display and the colors module.

pgstar is MESA’s built-in plotting system – it opens a window that updates every few timesteps so you can watch the star evolve in real time. The colors module computes synthetic photometry from the model atmosphere at each step, which is what gives you the 2MASS magnitudes in the history file.

The namelist below configures five panels:

• HR diagram
• $\Delta \nu$ vs age
• $\delta {\nu }_{02}$ vs age
• Interior abundance profile
• 2MASS magnitude track

Copy this into the your inlist_run, making the &pgstar section:

&pgstar

  file_white_on_black_flag = .false.

  Grid1_win_flag = .true.
  Grid1_win_width = 18
  Grid1_win_aspect_ratio = 0.56

  Grid1_file_flag = .true.
  Grid1_file_dir = 'pgplot'
  Grid1_file_prefix = 'grid_'
  Grid1_file_interval = 10
  Grid1_file_width = 18

  Grid1_num_cols = 9
  Grid1_num_rows = 9
  Grid1_num_plots = 5

  Grid1_xleft = 0.01
  Grid1_xright = 0.99
  Grid1_ybot = 0.02
  Grid1_ytop = 0.98

  ! Panel 1 -- HR diagram (top left)
  Grid1_plot_name(1) = 'HR'
  Grid1_plot_row(1) = 1
  Grid1_plot_rowspan(1) = 4
  Grid1_plot_col(1) = 1
  Grid1_plot_colspan(1) = 3
  Grid1_plot_pad_left(1) = 0.06
  Grid1_plot_pad_right(1) = 0.02
  Grid1_plot_pad_top(1) = 0.05
  Grid1_plot_pad_bot(1) = 0.06
  Grid1_txt_scale_factor(1) = 0.65
  HR_title = 'HR diagram'

  ! Panel 2 -- Delta_nu vs age (top centre)
  Grid1_plot_name(2) = 'History_Track1'
  Grid1_plot_row(2) = 1
  Grid1_plot_rowspan(2) = 4
  Grid1_plot_col(2) = 4
  Grid1_plot_colspan(2) = 3
  Grid1_plot_pad_left(2) = 0.06
  Grid1_plot_pad_right(2) = 0.02
  Grid1_plot_pad_top(2) = 0.05
  Grid1_plot_pad_bot(2) = 0.06
  Grid1_txt_scale_factor(2) = 0.65
  History_Track1_title = 'Large frequency separation'
  History_Track1_xname = 'star_age'
  History_Track1_yname = 'Delta_nu_int'
  History_Track1_xaxis_label = 'Age (yr)'
  History_Track1_yaxis_label = 'Delta-nu (Hz)'
  History_Track1_ymin = 0
  History_Track1_ymax = 5e-4
  History_Track1_reverse_xaxis = .false.



  ! Panel 3 -- delta_nu02 vs age (top right)
  Grid1_plot_name(3) = 'History_Track2'
  Grid1_plot_row(3) = 1
  Grid1_plot_rowspan(3) = 4
  Grid1_plot_col(3) = 7
  Grid1_plot_colspan(3) = 3
  Grid1_plot_pad_left(3) = 0.06
  Grid1_plot_pad_right(3) = 0.04
  Grid1_plot_pad_top(3) = 0.05
  Grid1_plot_pad_bot(3) = 0.06
  Grid1_txt_scale_factor(3) = 0.65
  History_Track2_title = 'Small frequency separation'
  History_Track2_xname = 'star_age'
  History_Track2_yname = 'delta_nu02_int'
  History_Track2_xaxis_label = 'Age (yr)'
  History_Track2_yaxis_label = 'delta-nu02 (Hz)'
  History_Track2_ymin = 0
  History_Track2_ymax = 5e-5


  ! Panel 4 -- Interior composition (bottom left)
  Grid1_plot_name(4) = 'Abundance'
  Grid1_plot_row(4) = 5
  Grid1_plot_rowspan(4) = 4
  Grid1_plot_col(4) = 1
  Grid1_plot_colspan(4) = 5
  Grid1_plot_pad_left(4) = 0.06
  Grid1_plot_pad_right(4) = 0.02
  Grid1_plot_pad_top(4) = 0.04
  Grid1_plot_pad_bot(4) = 0.06
  Grid1_txt_scale_factor(4) = 0.65
  Abundance_title = 'Interior composition'
  Abundance_num_isos_to_show = 4
  Abundance_which_isos_to_show(1) = 'h1'
  Abundance_which_isos_to_show(2) = 'he4'
  Abundance_which_isos_to_show(3) = 'c12'
  Abundance_which_isos_to_show(4) = 'n14'
  Abundance_xaxis_name = 'mass'
  Abundance_log_mass_frac_min = -4.0

  Grid1_plot_name(5) = 'History_Track3'
  Grid1_plot_row(5) = 5
  Grid1_plot_rowspan(5) = 4
  Grid1_plot_col(5) = 6
  Grid1_plot_colspan(5) = 4
  Grid1_plot_pad_left(5) = 0.06
  Grid1_plot_pad_right(5) = 0.04
  Grid1_plot_pad_top(5) = 0.05
  Grid1_plot_pad_bot(5) = 0.06
  Grid1_txt_scale_factor(5) = 0.65
  History_Track3_title = '2MASS mags'
  History_Track3_xname = 'J'
  History_Track3_yname = 'H'
  History_Track3_xaxis_label = 'M_J'
  History_Track3_yaxis_label = 'M_H'

/ ! end of pgstar namelist

Now for the colors module. The namelist below uses paths that are interpreted relative to the directory you launch ./rn from (your Lab 3 working directory). They point to the ../data/ folder that ships with this lab, so they should resolve as they are. Before you paste the namelist in, confirm that they do – for example:

ls ../data/stellar_models     # Kurucz2003all__alpha_00, vega_flam.csv, ...
ls ../data/filters/2MASS      # the 2MASS filter throughputs

If a file is missing, the same data also ships with MESA under $MESA_DIR/data/colors_data – point the path there instead. Once every path checks out, add this &colors namelist to your inlist_run:

&colors

   use_colors = .true.                                              !use colors

   instrument = '../data/filters/2MASS/2MASS'                       !We are assuming that data is in this dir
   stellar_atm = '../data/stellar_models/Kurucz2003all__alpha_00'   !check to see where it actually is. 

   vega_sed = '../data/stellar_models/vega_flam.csv'                !same as above
   mag_system = 'Vega'                                              

   distance = 3.0857d19                                             !10 parsecs in cm -> absolute magnitudes

   make_csv = .true.                                                !Make a csv for each filter
   colors_results_directory = 'SED'                                 !put them in the SED/ directory
   sed_per_model = .false.                                          !overwrite them at every step

/ ! end of colors namelist

The distance = 3.0857d19 sets the distance to 10 parsecs (expressed in cm), which is what puts the magnitudes on the absolute scale you used in Lab 2.

Step 3 – Run the model

./clean 
./mk
./rn

The model will run from the pre-main sequence to Terminal Age Main Sequence. In the pgstar window, $\Delta \nu$ is the top-centre panel (“Large frequency separation”) and $\delta {\nu }_{02}$ is the top-right panel (“Small frequency separation”). As it runs, pay attention to:

• How quickly does $\Delta \nu$ change compared to the HR diagram position?
• How does $\delta {\nu }_{02}$ behave – does it change monotonically?
• What is happening to the interior composition at the same time?

Step 4 – Let’s think about this as the model runs – Age diagnostics: seismic vs CMD

As the model runs, compare what each observable is actually telling you.

$\Delta \nu$ as an age clock. Because $\Delta \nu \propto \sqrt{\bar{\rho }}\propto \sqrt{M/R^3}$ , it tracks the mean density. On the main sequence the star expands slowly, so $\Delta \nu$ decreases – but the rate depends on mass. This makes $\Delta \nu$ a useful density indicator but a coarse age clock on its own, because you need to know $M$ independently to convert density to age.

$\delta {\nu }_{02}$ as an age clock. The small separation probes the sound speed gradient in the core. As hydrogen burns, the mean molecular weight of the core increases, the sound speed drops, and the gradient steepens – so $\delta {\nu }_{02}$ decreases monotonically with central hydrogen abundance. This makes $\delta {\nu }_{02}$ a nearly mass-independent age indicator along the main sequence at fixed $\Delta \nu$ . This is the basis of the Christensen-Dalsgaard (C–D) diagram: plotting $\delta {\nu }_{02}$ versus $\Delta \nu$ produces a grid where lines of constant mass and constant age cross at different angles, allowing both to be read off simultaneously from two observables.

Comparison with Lab 2. In Lab 2 you measured stellar ages from CMD position – the $J-K_s$ colour and absolute $K_s$ magnitude.

Step 5 – Plotting beyond pgstar with Python

Once the run has enough history data, use mesa_reader to reproduce the four key plots. The python_helpers/ directory contains more complete plotting scripts – the code below is a minimal example you can run directly.

#make sure you are in the working labs module. running an 'ls' or 'pwd' should tell you where you are. 
cd Lab3                     # or Lab3_start/Lab3 or wherever your labs are

python                      #open python live and paste the below code into it one by one to see the plots.

import mesa_reader as mr
import matplotlib.pyplot as plt

h = mr.MesaData('LOGS/history.data')

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# HR diagram
axes[0, 0].plot(h.log_Teff, h.log_L)
axes[0, 0].invert_xaxis()
axes[0, 0].set_xlabel(r'$\log\,T_\mathrm{eff}$')
axes[0, 0].set_ylabel(r'$\log\,L/L_\odot$')
axes[0, 0].set_title('HR diagram')

# CMD -- 2MASS Vega-system columns written by the colors module
# Verify the exact column names with: head -7 LOGS/history.data
axes[0, 1].plot(h.J - h.Ks, h.Ks)
axes[0, 1].invert_yaxis()
axes[0, 1].set_xlabel(r'$J - K_s$')
axes[0, 1].set_ylabel(r'$K_s$ (mag)')
axes[0, 1].set_title('CMD')

# Delta_nu vs age
axes[1, 0].plot(h.star_age / 1e9, h.Delta_nu_int * 1e6)
axes[1, 0].set_xlabel('Age (Gyr)')
axes[1, 0].set_ylabel(r'$\Delta\nu$ ($\mu$Hz)')
axes[1, 0].set_title('Large frequency separation')

# delta_nu02 vs age
axes[1, 1].plot(h.star_age / 1e9, h.delta_nu02_int * 1e6)
axes[1, 1].set_xlabel('Age (Gyr)')
axes[1, 1].set_ylabel(r'$\delta\nu_{02}$ ($\mu$Hz)')
axes[1, 1].set_title('Small frequency separation')

plt.tight_layout()
plt.savefig('lab3_history.png', dpi=150)
plt.show()

3D CMD – lift the CMD into 3D using $\Delta \nu$ as the third axis, revealing how the seismic observable evolves along the sequence:

import mesa_reader as mr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

h = mr.MesaData('LOGS/history.data')

fig = plt.figure(figsize=(9, 7))
ax  = fig.add_subplot(111, projection='3d')

ax.scatter(h.J - h.Ks, h.Ks, h.Delta_nu_int * 1e6, s=5, c=h.star_age, cmap='plasma')

ax.set_xlabel(r'$J - K_s$')
ax.set_ylabel(r'$K_s$')
ax.set_zlabel(r'$\Delta\nu$ ($\mu$Hz)')
ax.invert_yaxis()
plt.tight_layout()
plt.show()

Step 6 – Crowd-source the seismic grid

Run the snippet below to extract values at 1, 3, 5, 7, and 9 Gyr, then enter one row per age into the shared spreadsheet:

import mesa_reader as mr
import numpy as np

h = mr.MesaData('LOGS/history.data')

target_ages = [1e9, 3e9, 5e9, 7e9, 9e9]

for age in target_ages:
    idx = np.argmin(np.abs(h.star_age - age))
    print(f"Age: {h.star_age[idx]/1e9:.1f} Gyr  "
          f"Teff: {h.Teff[idx]:.0f} K  "
          f"L/Lsun: {10**h.log_L[idx]:.4f}  "
          f"J-Ks: {h.J[idx] - h.Ks[idx]:.4f}  "
          f"Delta_nu: {h.Delta_nu_int[idx]*1e6:.2f} uHz  "
          f"delta_nu02: {h.delta_nu02_int[idx]*1e6:.2f} uHz")

…and to check the mass from the inlist (we can check the inlist as we earlier showed that the inlist file is made as a copy of the inlist_run file when we do ./rn)

with open('inlist') as f:
    for line in f:
        if 'initial_mass' in line:
            print(line.strip())

awk '
!hdr {
  for (i = 1; i <= NF; i++) if ($i == "star_age") { for (j = 1; j <= NF; j++) c[$j] = j; hdr = 1 }
  next
}
{
  n++
  age[n]=$(c["star_age"]); teff[n]=$(c["Teff"]); logL[n]=$(c["log_L"])
  jmag[n]=$(c["J"]); ksmag[n]=$(c["Ks"])
  dnu[n]=$(c["Delta_nu_int"]); dnu02[n]=$(c["delta_nu02_int"])
}
END {
  m = split("1 3 5 7 9", t, " ")
  for (k = 1; k <= m; k++) {
    target = t[k] * 1e9
    best = 1; bd = 1e99
    for (i = 1; i <= n; i++) { d = age[i] - target; if (d < 0) d = -d; if (d < bd) { bd = d; best = i } }
    i = best
    printf "Age: %.1f Gyr  Teff: %.0f K  L/Lsun: %.4f  J-Ks: %.4f  Delta_nu: %.2f uHz  delta_nu02: %.2f uHz\n", \
           age[i]/1e9, teff[i], 10^logL[i], jmag[i]-ksmag[i], dnu[i]*1e6, dnu02[i]*1e6
  }
}
' LOGS/history.data

grep initial_mass inlist

Add the results to the google sheets file (Lab3 tab). This link : google sheets

Once the full group has contributed, look at the complete grid. How well do $\Delta \nu$ and $\delta {\nu }_{02}$ separate stars of different masses at the same age? How does this compare to the CMD separation from Lab 2?

Age constraints from seismology are only useful for stars we can actually observe oscillating. The next step takes the grid you just built and allows us to probe which of these stars are detectable, and with what instrument.

Step 7 – Who cares if we can’t even measure it?

Get your own copy of the notebook

The analysis lives in a shared Google Colab notebook. Open the link below, then make a personal copy before you do anything else – this is important so your edits don’t overwrite someone else’s work and vice versa.

→ Open the Lab 3 detection map notebook

Once it opens:

File → Save a copy in Drive

or click the “Copy to Drive” button in the toolbar at the top of the page. Either way, a personal copy lands in your Google Drive and the original is untouched. Work from your copy for the rest of the lab.

Download the data

In the shared Google Sheet:

File → Download → Comma Separated Values (.csv)

Save the file – you’ll upload it into the notebook using the file upload button in the first cell.

What the notebook does

The notebook reads the crowd-sourced grid and produces two detection map plots (photometric and RV), coloured by stellar mass with age annotated on each point. Horizontal threshold lines show the minimum detectable amplitude for each instrument and campaign length. Stars above a line are detectable; stars below are not.

The parameters block near the top of the notebook lets you change instrument names, noise floors (${\sigma }_0$ ), and campaign lengths – the threshold lines update when you rerun the cell.

What to look for

The y-axis is logarithmic. The steep drop in amplitude towards red $J-K_s$ (K and M dwarfs) is immediately visible. What masses can we realistically use this technique for? (https://arxiv.org/abs/1103.0702)

Bonus Step – What actually matters for wobbly stars?

Lab 2 showed you that changing the atmospheric boundary condition (atm_T_tau_relation) and the mixing length parameter (${\alpha }_{\mathrm{MLT}}$ ) visibly shifts a star’s track on the HR diagram and CMD https://arxiv.org/pdf/2303.09596. Here you will ask: by how much do those same changes affect seismic observables?

The answer matters because if $\Delta \nu$ and $\delta {\nu }_{02}$ are insensitive to surface physics, they give more trustworthy ages than CMD position – the seismic signal comes from the core, not the atmosphere. This does not mean they are completely immune to surface physics, we are investigating just how robust these observables are.

Setup

Keep your Lab 3 mass. Run a small grid varying one parameter at a time, recording seismic and photometric values at the same target ages like before:

Vary the atmospheric boundary condition (fix ${\alpha }_{\mathrm{MLT}}=1.8$ ):

atm_T_tau_relation = 'Eddington'       ! reference
atm_T_tau_relation = 'solar_Hopf'
atm_T_tau_relation = 'Krishna_Swamy'
atm_T_tau_relation = 'Trampedach_solar'

Vary the mixing length parameter (fix atm_T_tau_relation = 'Eddington'):

mixing_length_alpha = 1.5
mixing_length_alpha = 2.0
mixing_length_alpha = 2.5
mixing_length_alpha = 3.0

Extract the values

Use the same snippet as Step 6, just targeting your chosen ages:

import mesa_reader as mr
import numpy as np

h = mr.MesaData('LOGS/history.data')

target_ages = [1e9, 3e9, 5e9, 7e9, 9e9]

for age in target_ages:
    idx = np.argmin(np.abs(h.star_age - age))
    print(f"Age: {h.star_age[idx]/1e9:.1f} Gyr  "
          f"Teff: {h.Teff[idx]:.0f} K  "
          f"L/Lsun: {10**h.log_L[idx]:.4f}  "
          f"J-Ks: {h.J[idx] - h.Ks[idx]:.4f}  "
          f"Delta_nu: {h.Delta_nu_int[idx]*1e6:.2f} uHz  "
          f"delta_nu02: {h.delta_nu02_int[idx]*1e6:.2f} uHz")

Enter results into the shared spreadsheet

Go to the Lab 3 Bonus tab in the shared spreadsheet. Enter one row per run. Columns D and E are the parameters you changed; columns F–J are the MESA output. Columns K–M calculate automatically.

The four bar charts at the bottom of the sheet show how much Teff, $J-K_s$ , $\Delta \nu$ , and $\delta {\nu }_{02}$ vary across all models. Compare the vertical scale of the seismic panels to the photometric ones.