Minilab 3: Meridional Circulation

Introduction

So far in our exploration of rotating stellar models, we have seen how instabilities can be modeled in MESA to approximate particular behavior present in 2D-ESTER models. We have also gone further to calculate the Eddington-Sweet circulation velocity in the 1D case and seen some differences against analogous 2D references.

So what does this mean? Of course a 1D model and a 2D model are going to be different, right?

The good news is that we are not consigned to spend hours on HPC clusters for every test case. Instead, given a selective set of more expensive 2D models, we can tune our MESA models. Luckily, we have various levers we can turn to accomplish this. If we can establish an idea of what the velocity field inside the star should be, then we can implement it into MESA with hooks in run_star_extras.f90. In this lab, we will be doing just that. Beginning from a scaffolded starting point, you will implement a custom torque and angular momentum mixing routine in MESA, outputting data that can be compared against the tracks in Mombarg et al., 2023¹ and Mombarg et al., 2024². By tuning the torque, you can mimic the effect of meridional circulation on the behavior and evolution of the star, as given by these 2D ESTER models.

For a discussion on meridional circulation, click the hint below. This information was all covered in more detail during the lecture, so students are encouraged to only visit this discussion as needed, or while waiting for model runs to complete.

Helpful Links

The lab setup can be found here as a mirror of the Google Drive link found below. This folder contains partial solutions (separated by task), full solutions for each test case, the starting point for the lab, pre-made ZAMS models, and a series of pre-computed runs with no meridional circulation (called no_extra_jdot).

The Google Colab script to make the bonus plots can be found here (local copy) .

Instructions

Step 0: Getting Started

This starting point should be a fairly familiar set of files. Each of these files has been -mostly- pre-prepared with the structure we will need to get going. Additionally, throughout each of these starting point files, variables that need to be changed are explicitly marked with “!!!! !!!!”. This is also true if there is some particular section that needs your input. If you have trouble finding a value, feel free to ctrl+f your way around.

The working directory should now be:

At this stage, we are now ready to dive into some inlists!

Step 1: Inlist Project

Let’s start by looking over inlist_project.

This file should be generally quite familiar from Minilabs 1 & 2. Starting from the top of the file in &star_jobs, you will need to turn on pgstar, load from the saved ZAMS model, and set the initial rotation rate. It is also recommended that you turn on the pause_before_terminate flag to ensure that you can view the pgstar plots before the script exits.

Next, in &controls, set the output directory for the logs, turn on the other angular momentum flag, and turn on the other torque flag. For the log directory, use a standard naming convention. It is recommended that this be something like M<mass>_Omega<initial rotation> (ie. M05_Omega0p60), but generally this is up to you. The use_other_torque and use_other_am_mixing flags are used to provide MESA with a custom subroutine that steps into the execution process and changes key behavior about how the torque and angular momentum mixing are calculated. You can explore all the available use_other_<hook> in the MESA docs here.

We have now set the general starting conditions for the MESA model (from ZAMS). Included in these changes is a directive for MESA to look beyond its regular routines for the torque and angular momentum calculations, which we will deal with later.

Step 2: Inlist Pgstar

“If a MESA model runs to TAMS, and no one sees the pgstar plot, did it even happen?”
— Me, circa 2025

The provided inlist_pgstar is the same as the earlier labs for the day with two modifications. First, the bottom left plot (previously showing gravitational darkening) has been updated to plot star_age on the X-axis and log_total_angular_momentum on the Y-axis. Second, the top right plot (previously the mixing profile) has been updated to plot radius on the X-axis and log_dJ_over_J on the Y-axis. log_dJ_over_J is a new variable that we will define later that can be thought of as the contribution of our planned additional torque on the rotation profile of the star.

Step 3: History/Profile Columns List

Now, we need to add some more information to the history and profile outputs of our model. Most of the values have already been added.

Surf_avg_omega, surf_avg_omega_div_omega_crit, surf_avg_v_rot, and center_omega have already been uncommented in history_columns.list. Radius, omega, r_polar, j_rot, and omega_crit have already been uncommented in profile_columns.list.

Step 4: Run Star Extras

As a final step before we can run the model, we need to modify src/run_star_extras.f90 to include an additional torque from meridional circulation and increase the viscosity to match the sample 2-D ESTER models.

Step 4.1: Subroutines, A New Hope

Start by looking over the run_star_extras file. Take note of the general form of the code. You will notice that this file is really a collection of labelled subroutines and functions (particularly labelled ‘integer functions’). On a high level, these are just the basic execution steps that MESA uses to solve the model. You can see when they get called in the flowchart below. This same flowchart can be found in the MESA docs here.

In Step 1, we turned on two flags to use other_torque and other_am_mixing. Since we have not reached an AI Singularity (yet), MESA cannot intuit what these custom subroutines are or where they will be. Hence, we need to provide a pointer that says “This new procedure is defined by this other piece of code” in the subroutine extras_controls. The form of this pointer can be seen elsewhere in extras_controls or in the MESA docs here. Do not worry about what these custom subroutines are yet, we will cover that in the next step.

s% <star procedure> => <local subroutine>

s% other_torque => meridional_circulation
s% other_am_mixing => additional_nu

Step 4.2: Fortran Strikes Back

Now that MESA knows where to look, what exactly is going on in these new subroutines? In run_star_extras, scroll down to our custom other torque subroutine, meridional_circulation.

meridional_circulation is declared with the subroutine keyword, meaning we do not expact any output. Instead, this routine will grab the object containing all the star’s information, identified with the pointer s%, and modify values directly. You will also see that this subroutine takes in two values, id and ierr. id is the unique identifier that is tied to each instance of star_info, the object which holds everything about the star. ierr is an integer passed across MESA to keep tabs on the status of each operation. If this value becomes non-zero, it means that an error has occurred.

Next, we have variable declarations. In Fortran, ALL variables must be declared before they are used. This includes arrays, which have the added complication of needing to be allocated as well! The types of these variables should be explicitly provided as well. In fact, because run_star_extras.f90 contains the implicit none statement at the beginning of the file, this explicit declaration is not optional. We will need three additional integers for this subroutine: k, k0, and nz. k will be used as a counter value for the index within the model, k0 will be the index where the radiative envelope starts, and nz will hold onto the total number of zones in the model. These values will be described in more detail later.

<type> :: <variable name>

<type> :: <variable name0>, <variable name1>, <variable nameN>

integer :: <variable name>

integer :: k, k0, nz

Following this section, we reset the ierr integer to 0 (meaning success), grab the relevant pointer to our star object, then check if anything went wrong. This portion of the routine is functionally boilerplate. You will find the same three lines across all the routines that access an instance of the star.

The next two sections set up some important constants for subsequent calculation steps. The first, nz, is the integer number of zones in the model. This value can be accessed directly from the star pointer. The second value, k0, is the index of the zone where the radiative envelope starts. This value is not simply stored in the star pointer. Instead, we have to use a DO loop which can be thought of in the same way as a FOR loop in Python, C, or Java. This loop starts at the center and works it way to the surface, setting k0 to the first index where the mixing_type == minimum_mixing. Minimum_mixing means that the amount of diffusive mixing, D_mix, is being set to the minimum value specified in our inlist_project of 1d3. We should expect this mixing_type to switch from convective mixing to minimum_mixing only once we have entered the envelope, as we leave the convective core and enter the radiative envelope of our star. Therefore, the index where this value first switches is the base of the radiative envelope.

Next, we allocate and calculate three arrays: U_r, mer_comp, and dmer_comp_dr. U_r will be a simple velocity field with a constant scale factor, C, that should improve the comparison to the 2D ESTER models. Note, U_r is a by-eye linear fit to the 2D runs and NOT an actual, physical relation. The equation we will use for this field is:

mer_comp is then the meridional torque produced by that velocity field, given by:

dmer_comp_dr is then differential meridional torque per unit radius:

U_r = (s% r) / (s% r(1)) * 0.0001 ! cm/s
mer_comp = s% rho * pow4(s% r) * (s% omega) * U_r

We also want to store a copy of dmer_comp_dr with the star object in memory. This can be done by passing this value over to an xtraN_array, where N is an integer 1 through 6. These six of these arrays are declared within star_info to hold double precision values and are guaranteed to have a length equal to s% nz. With this array, we then call a weighted smoothing subroutine to decrease the amount of noise and increase how pretty our end plots are. We will not be exploring how the weighted smoothing subroutine works in this lab, but feel free to explore it at a later time after the school, if interested.

So far, we have been doing a lot of prep work collecting variables. The next sections of the subroutine are all directed at finally applying these new torques to the model. We start at the surface and work our way down to the radiative envelope boundary k0 that we calculated previously. At each zone k we encounter, using our stored copy of dmer_comp_dr(k), we get the rate at which the specific angular momentum is changed as a result of our velocity field at t=0 and add it to whatever the value of extra_jdot was at t=-1.

Once we reach the convective core, we calculate the total torque that we applied to the envelope (saving it to a temporary value, s% xtra(1)). Then, to ensure that angular momentum is conserved across the entire star, we apply a opposing torque to zones in the core. This torque is scaled from the total torque on the envelope and in the opposite direction.

Pay attention to the places that we saved data in the star pointer, s%. Specifically, s% xtra(1) is holding onto the total torque applied to the envelope, a scalar. Meanwhile, s% xtra1_array and s% xtra2_array are storing the values of extra_jdot and dmer_comp_dr as arrays. You will need some of these values later.

You may remember that we also pointed to another subroutine, additional_nu. We will not be exploring how that subroutine works, but feel free to explore it at your own pace after the school. In practice, it is just a few DO loops to increase the viscosity and help bring the 1D model further into agreement with the 2D ESTER model.

Before we continue, we need to check that all the updates made thus far are working as expected.

Step 4.3: Return of the Profile Columns (from Lab 2)

So, we have made the necessary calculations and saved off some variables into the star pointer. Now, we need to make those values presentable in the history and profile columns.

Now that MESA expects some additional values, lets add the history data first.

names(<column number>) = <your column name>
vals(<column number>) = <your value>

names(1) = "total_torque_envelope"
vals(1) = s% xtra(1)

Now, we will need to do some calculations to get the data for the profile columns. Remember, since these are profiles, they are arrays with entries at each zone k. We will be adding two columns called extra_jdot and log_dJ_over_J. extra_jdot is the same value we encountered in meridional_circulation. log_dJ_over_J is an account of how much of the specific angular momentum was due to our additional velocity field. value can be calculated as,

where ${\dot{J}}_{\mathrm{extra}}$ is extra_jdot.

To put these values into columns, we need to make a DO loop across all the zones in the model, saving data in each one. An example of how to do this with a variable named beta is already given in the subroutine, data_for_extra_profile_columns.

names(<column number>) = <your column name>
do k = 1, nz
    vals(k, <column number>) = <your value at k>
end do

log10()
abs()

s% xtra1_array(k)
s% dt
s% j_rot(k)

names(1) = "extra_jdot"
names(2) = "log_dJ_over_J"
do k = 1, nz
    vals(k, 1) = s% xtra1_array(k)
    vals(k, 2) = log10(abs(s% xtra1_array(k) * s% dt / s% j_rot(k)))
end do   

Step 5: Running the model

Let’s see MESA in action.

You should see the rotation begin at your set initial value. Further along, the plot for log_dJ_over_J will display the additional torque being forced along the radius of the star. How does this value relate to the convective zone boundary? How noisy is the plot? Do you see any information about how energy is being transported?

Once the model completes, take note of the following values as printed in the pgstar plot:

• center_omega
• surf_avg_omega
• star_age

You were also provided completed output logs & pgstar plots for each of the mass cases with extra_jdot set to a near-zero value, representing models with no meridional circulation. Importantly, in these runs extra_jdot was explicitly set to 1d-99 at each zone, k , to avoid log_dJ_over_J from becoming undefined.

Pgstar Output with meridional circulation (M08, Omega 0.75)

Pgstar Output without meridional circulation (M08, Omega 0.75)

Congratulations, you have completed Lab 3! Feel free to pursue the following two bonus exercises if you have additional time.

BONUS: Recreating the results from Mombarg et al., 2024²

BONUS: Exploring other velocity field scale factors

So far we have been working through the models with an ad-hoc approximation of $U_r$ with C = 1e-4.

Primary References