Lab 3: The Hertzsprung Progression

Lab written by: Andy Santarelli (lead TA), Mathijs Vanrespaille (lead TA), Lynn Buchele, Sofia Mesini, and Ebraheem Farag (lecturer)

Background

In this lab you will take one of your saved Cepheid models from Lab 1, use a nonlinear pulsation setup to kick it into motion, and inspect the resulting waveform. Your goal is to identify where the bump appears in the cycle and combine your result with the rest of the class to reconstruct the Hertzsprung progression.

Many classical Cepheids show a distinctive “bump” in their waveform. The figure below shows a few examples of folded Cepheid light curves observed by the OGLE study.

The location of that bump changes with the pulsation period. In the standard single mode picture, this is related to a near 2:1 resonance between the second overtone and the fundamental mode, so a useful quantity to keep in mind is

The linear nonadiabatic (LNA) plot below shows the overtone period ratios as a function of the fundamental mode period, $P_0$ . It shows why this ratio is useful: the second overtone ratio, $P_2/P_0$ , passes through the resonance range near the periods where bump Cepheids become most morphologically interesting.

As the stellar structure changes across the instability strip, the bump shifts from the descending branch, through the middle of the cycle, and onto the rising branch. In this lab, you will see that progression directly in nonlinear MESA models.

Science Goals

1. Evolve a selected Cepheid model into a finite amplitude nonlinear pulsation.
2. Measure the pulsation period and decide whether the model develops a bump.
3. Compare the bump location with the rest of the class sample to reconstruct the Hertzsprung progression.

Lab Goals

1. Choose a starting .mod file from your Lab 1 output or the shared sample.
2. Run the nonlinear setup until the waveform is developed enough to classify.
3. Record the period, bump location, and comments needed for the class comparison.

MESA Goals

1. Load a saved Lab 1 .mod file into the nonlinear MESA star TDC setup.
2. Use a GYRE kick to seed the fundamental radial mode and start the pulsation on lab timescales.
3. Use the history output and PGSTAR diagnostics to inspect the light curve, radius curve, velocity curve, and bump morphology.

Lab Directions

For this lab you will start from a good Cepheid model, run the nonlinear setup, inspect the waveform shape, and decide where the bump appears.

Setting up the work directory

Download lab3_work_dir.zip from this Google Drive, unzip it into some empty directory, and cd into the resulting lab3_work_dir/ directory. You’ll see that it already contains the inlists you will need. However, we need to provide TDC with a starting model to initialize the nonlinear run.

We want to use a Lab 1 .mod file, ideally the one you identified during Lab 2, whose GYRE fundamental mode is unstable and is growing faster than the second overtone: F_growth > 0 and F_growth > O2_growth. Ideally, the fundamental mode also grows faster than the first overtone, F_growth > O1_growth, and the model is closer to the red edge of the instability strip. If you cannot find a model with F_growth > O1_growth, use a redder model with F_growth > 0 and F_growth > O2_growth. Models closer to the blue edge are more likely to switch into an overtone during the nonlinear run.

Copy your Lab 1 .mod files into the Lab 3 work directory:

cp -r /path/to/your/lab1/mod_dir/ .

Alternatively, you can download the models from the Lab 1 mod file solutions, which are grouped by mass. If you do not have a useful saved model ready to run, use the Lab 3 nonlinear start model files, which contains the selected .mod files used for the shared Lab 3 sample.

Task 1: Choose a Starting Model

Take a look inside mod_dir/. These are the saved stellar structures that Lab 3 can use.

The filenames are written in the form

modelNumber_currentMass_effectiveTemperature_luminosity.mod

If you are using the Lab 3 nonlinear start model files instead of your own Lab 1 mod_dir/, the selected files have a longer name with the period and initial mass prepended. In those filenames, the final four underscore separated fields before .mod are still the Lab 1 model number, current mass, effective temperature, and luminosity.

Question: What would make an interesting model to simulate in further detail using TDC?

Task 2: Edit the Lab 3 Inlist

Open inlist_pulses in your editor.

Find the line

load_model_filename = 'mod_dir/YOUR_MODEL.mod'

For your first run:

• update load_model_filename so it points to the model you copied into your local mod_dir/

Task 3: Choose and set an initial kick

It can take a very long time for a MESA TDC model to start pulsating “naturally”. Therefore, we get the pulsation going by giving the model an initial kick. To that end, the eigenfunction of the fundamental pulsation mode, which expresses how the star is perturbed by the pulsation, is computed using GYRE and scaled to an initial kick velocity that you choose.

The closer this kick is to the final pulsational radial velocity, the faster the bump in the light curve will develop. From the figure below, read off a reasonable initial kick for your chosen model.

Now add this value into your inlist_pulses. Question: Can you find which variable stores the initial kick?

    x_ctrl(6) = 10d0 ! initial vsurf (kms)

The kick is also a mode selection choice. Here we seed the fundamental radial mode because the class goal is the fundamental mode Hertzsprung progression. For a production calculation, especially near regions where multiple modes are unstable, you would test that the final state does not depend strongly on the exact kick amplitude or on which mode was initially excited.

Task 4: Compile and Run the Model

First compile the work directory:

./clean
./mk

If the compilation succeeds, start the nonlinear run:

./rn

Task 5: Watch the Diagnostics

The PGSTAR panels are by far the most important diagnostics for this lab. For most students, they already summarise nearly everything that matters for deciding whether the pulsation is developing well and whether a bump is visible.

If you want to know where the supporting output is written, the run also saves files in:

• png_pulsation/
• LOGS_pulsation/
• photos/

Task 6: Decide Whether the Run Is Good Enough

For the purpose of this lab, the run is useful once you can see that the kick has produced a coherent pulsation and the amplitude is either:

• still clearly growing
• or close to a repeating finite amplitude cycle

In the drop down, we provide some information on what standard should be applied to research quality nonlinear runs. Feel free to read this after starting your model.

Signs that the run is doing the right thing:

• growth is positive for at least part of the run
• delta_R, delta_logL, or delta_Mag are no longer consistent with numerical noise
• the light, radius, or velocity curves begin to repeat from cycle to cycle

You can see an example of healthy, developed pulsation below.

The four panels labeled with a red number are the most relevant. They show:

1. Hertzsprung-Russell diagram. We expect an ellipsoidal path for a purely sinusoidal waveform without a bump. In this example however, you can see a downward kink around $\log T_{\mathrm{eff}}=3.7$ .
2. effective temperature variation over time in blue and luminosity variation in solar luminosity over time, also called the light curve, in yellow. This example is clearly not purely sinusoidal as the peak is too thin and there’s a subtle bump just before the peak. Any bumps are usually relatively easy to identify in this panel. You can see other examples of these light curves in the figure in the introduction or at the OGLE catalog.
3. radial variation and surface radial velocity over time in yellow and blue, respectively. Once again, neither curve is sinusoidal, although the bump is not clearly visible here.
4. radial velocity profile against optical depth as a proxy of depth from the surface. The initial kick should be plainly visible here after a few hundred time steps.

Signs that you should stop and rethink:

• the run exits immediately
• the model never develops a clear periodic signal
• the waveform looks obviously pathological rather than pulsational

./re

./re photo_name

Profile_Panels3_xaxis_reversed = .false.

Task 7: Inspect the Waveform

Now look at the waveform and decide where the bump appears in the cycle.

Start with the PGSTAR plot. For almost everyone, this will be the easiest and clearest place to identify the bump.

If you want to check what you are seeing, you can also look at:

• the saved plots in png_pulsation/
• the time series in LOGS_pulsation/history.data

If possible, inspect more than one diagnostic. The bump is often easiest to see in a light curve in luminosity or delta_Mag, but the radius and surface velocity curves can help you decide whether a feature is real.

Use the following simple classification:

• descending branch: the bump appears after maximum light while the curve is falling
• middle: the bump is near the middle of the cycle and not clearly tied to either branch
• rising branch: the bump appears before maximum light while the curve is rising
• no clear bump: use this only if the waveform is genuinely ambiguous

Task 8: Record Your Result

Add one row for your successful model to the shared class table.

At minimum, record:

• model filename
• initial mass
• log L and log T_eff, if available from Lab 1
• fundamental period
• whether the pulsation was clearly established
• bump classification
• a short note such as clear bump, weak bump, or needed restart

Once the class table starts to fill up, switch to the “Lab 3: Sorted by Period” tab on the google sheets and look for the bump progression across the sample.

Questions for Discussion

As the class table fills in, discuss these questions at your table:

• how does bump location change with period?
• where does the bump move from the descending branch to the rising branch?
• does the class sample support the idea that the morphology is tied to the $P_2/P_0\approx 0.5$ resonance?
• what finite amplitude information does the TDC waveform add beyond the Lab 2 periods and growth rates?
• did any model show behavior that was not cleanly single mode, such as overtone selection, alternating cycles, or a waveform that did not settle?

If You Are Still Waiting on a Run

These TDC runs can take a long time to converge, often more than 10 minutes. If your model is still running and you have some idle time, use that time to do one or more of the following:

• review your Lab 2 notes and make sure your expected period is written down
• look at the shared class table and decide which period range is still undersampled
• inspect the output files you already have and make a preliminary guess about the bump
• compare what you are seeing in PGSTAR with what appears in history.data

Task 9 (BONUS): If You Have Extra Time

If your group finishes the core lab early, here are the most useful next steps:

1. compare your TDC result directly with the linear information you already gathered in Lab 2
2. look more carefully at how the bump appears in luminosity, radius, and velocity together
3. compare your PGSTAR animation with the other students at your table
4. make a movie of your PGSTAR output

You do not need to complete all of these. Pick the next one that feels most useful.

Option 1: Compare Back to Lab 2

If you completed Lab 2, compare your nonlinear result with the linear information you already had for the same model.

Ask yourself:

• did a model with positive linear growth turn into a useful nonlinear pulsator?
• is the nonlinear period similar to the period you expected from the linear analysis?
• did the model you thought would be interesting actually produce a clear bump?

If you also estimated where $P_2/P_0$ is closest to $0.5$ , compare that expectation with the waveform shape you actually see in the TDC run.

Option 2: Compare Different Diagnostics

If you have a clearly pulsating model, compare the bump location in:

• behavior related to luminosity
• radius
• surface velocity

You may find that the bump is easier to identify in one diagnostic than another. Record that in your notes if it helps explain your classification.

Option 3: Compare with Other Students at Your Table

If several people at your table have useful runs, compare them directly:

• do the bumps appear in different parts of the cycle?
• do the PGSTAR animations suggest a smooth progression across period?
• which models develop the clearest bump?

Option 4: Making a movie

Isn’t that animated PGSTAR window neat? Unfortunately, it vanishes once you end the run. Luckily, a bunch of .png files are output by MESA, which can be used to recreate the animated PGSTAR plots. You could either flick through them in an image viewer or combine them into a proper movie. MESA includes tools to make such movies. To do so, run the following in your terminal:

images_to_movie "png_pulsation/*.png" my_Cepheid_movie.mp4

Troubleshooting

Bonus coding task: average the light curve over one cycle

If you would like a more coding focused extension, modify run_star_extras so that it measures a cycle averaged quantity from the nonlinear light curve and compares that average with the corresponding static value from the original model.

You can find a complete set of Lab 3 bonus solutions with these source changes already made. As with the starter, you still need to copy in a .mod file and update load_model_filename.

One possible version of this task is:

1. identify one full pulsation cycle after the model has reached a reasonably repeatable waveform
2. measure a quantity over that cycle, such as luminosity or magnitude
3. compute the time average over the cycle
4. compare that cycle averaged value with the static value from the original stellar model

For example, you might compare:

• cycle averaged luminosity versus the original static luminosity
• cycle averaged magnitude versus the magnitude implied by the static model

Here is one reasonable way to implement this:

Step 1: Find the relevant source files

The most useful files to inspect are:

• src/run_star_extras.f90
• src/run_star_extras_TDC_pulsation.inc
• src/run_star_extras_TDC_pulsation_defs.inc

The existing Lab 3 setup already computes per cycle quantities such as:

• period
• delta_logL
• delta_Mag
• KE_growth_avg

and writes them out through the extra history column machinery. That makes this a natural place to add one or two more derived quantities.

Step 2: Choose what you want to average

Start with one quantity only. Good choices are:

• luminosity
• log_L
• a magnitude like quantity

The simplest first version is to average luminosity over one completed pulsation cycle.

Step 3: Decide what you will compare against

Pick the corresponding static quantity from the original model. For example:

• cycle averaged luminosity compared with the model luminosity before the nonlinear pulsation becomes large
• cycle averaged magnitude compared with the magnitude implied by the static luminosity

You do not need to design a perfect scientific definition here. The point is to compare a static value with the value implied by the nonlinear cycle.

Step 4: Accumulate the average over one cycle

As the run advances, keep track of:

• the value of the quantity you are averaging
• the elapsed time associated with each sample
• the running weighted sum over the current cycle
• the total elapsed time over the current cycle

In other words, your code should conceptually build something like

over one pulsation cycle.

      real(dp) :: cycle_sum_L_dt, cycle_sum_dt
      real(dp) :: cycle_avg_L, cycle_avg_L_minus_static
      real(dp) :: static_L

Step 5: Reset the accumulators at the start of a new cycle

The existing pulsation code already keeps track of completed cycles and period level information. Use that logic to decide when one cycle has ended and the next has begun.

At the end of each completed cycle:

• compute the average
• save the result somewhere
• reset the running sums for the next cycle

         cycle_sum_L_dt = cycle_sum_L_dt + s% L(1) * s% dt
         cycle_sum_dt   = cycle_sum_dt   + s% dt

         if (cycle_sum_dt > 0d0) then
            cycle_avg_L = cycle_sum_L_dt/cycle_sum_dt
            cycle_avg_L_minus_static = cycle_avg_L - static_L
         end if

         cycle_sum_L_dt = 0d0
         cycle_sum_dt   = 0d0

Step 6: Expose the new quantity in the history output

Once you have computed your new cycle averaged quantity, add it to the custom extra history columns so it appears in history.data.

That means updating the part of the code that currently exports values like:

• period
• growth
• delta_R
• delta_Mag

You can either:

• replace one of the less important bonus outputs for your experiment, or
• increase the number of extra history columns and append your new quantity

      TDC_pulsation_how_many_extra_history_columns = 14

         names(i) = 'cycle_avg_L'; vals(i) = cycle_avg_L; i=i+1
         names(i) = 'cycle_avg_L_minus_static'; vals(i) = cycle_avg_L_minus_static; i=i+1

Step 7: Recompile and rerun

Because this bonus task edits source include files, rerun both ./clean and ./mk before continuing. A plain rebuild can miss changes in included source files.

If you want to start and kick a new TDC model after editing the Fortran source:

./clean
./mk
./rn

If you want to continue from a previously saved TDC run after recompiling:

./clean
./mk
./re

Step 8: Compare the nonlinear average with the static value

Once your new quantity appears in the output, compare:

• the cycle averaged value obtained with TDC
• the corresponding static value from the original model

You can do this for one cycle or for several successive cycles if the run is still evolving.

Step 9: Interpret what you find

Write down a short conclusion:

• are the two values nearly the same?
• is there a systematic offset?
• does the offset shrink or grow as the pulsation settles?
• does averaging luminosity directly give a different answer than averaging a magnitude like quantity?

Suggested Reading

• OGLE Atlas of Variable Star Light Curves, classical Cepheids
• OGLE Atlas of Variable Star Light Curves, Type II Cepheids
• Farag et al. 2026, self-consistent nonlinear classical Cepheid pulsations during stellar evolution with MESA
• Smolec 2014, mode selection in pulsating stars
• Bono, Marconi, and Stellingwerf 2000, the Hertzsprung progression
• Marconi et al. 2024, the Hertzsprung progression of classical Cepheids in the Gaia era