Lab 1: 1D rotating stellar models

Introduction

In the most simple case of a non-rotating single star, the only forces acting on a mass element are pressure and gravity. The resulting spherical symmetry implies that the stellar structure equations depend only on one spatial coordinate (mass or radius) and time. However, rotation breaks the spherical symmetry consequently affecting stellar surface parameters, abundances and even its course of evolution (see e.g. Palacios 2013). To properly account for all the effects associated with rotation we need a multidimensional approach to solve the stellar structure equations. Today, we will compute 1D MESA rotating models and compare them with 2D ESTER models (Mombarg et al, 2023; 2024).

MESA includes rotational effects by applying corrections to the equilibrium stellar structure equations to account for the effect of the centrifugal force and by solving two additional equations: the mixing of chemical elements and the angular momentum (AM) transport. The latter determines the evolution of the angular velocity $(\Omega )$ with time which is what we will focus on this lab 1.

The AM transport is mediated through diffusive and advective processes. In MESA, the AM transport is only included by diffusion using the following equation,

where $i$ is a shell specific moment of inertia, and ${\nu }_{\mathrm{AM}}$ is the turbulent viscosity which determines the efficiency of the AM transport. The first term on the right-hand side accounts for the diffusion transport and the second term accounts for contraction and expansion of the shells at constant mass.

In this lab 1 we will incorporate rotation in a 10 solar mass model using MESA. We will look at the effect of rotation in the global surface parameters and on the internal rotation profile by including different AM transport mechanisms.

Download our template

Create a 10 solar mass rotating main-sequence model

We will start by creating a rotating main-sequence model for a 10 solar mass star. To save some computation time we already computed the pre-main sequence (PMS) track beforehand and saved the model once it reached the zero-age-main-sequence (ZAMS). You will use this ZAMS model as the starting point of all your runs. The particular ZAMS model we provide matches the 2D ESTER models composition, which we will use later.

inlist_project

In this lab we only need to modify the &star_job and &controls sections of the inlist_project. You can find all the commands available in MESA for the &star_job section here and for the &controls section here.

star_job

Let’s start by modifying the &star_job section of the inlist_project. First, you need to indicate that you do not want to create a PMS model. Then, you need to load the precomputed model_ZAMS in your work directory to start the computation from that point.

create_pre_main_sequence_model = .false.
load_saved_model = .true.
load_model_filename = 'model_ZAMS'

The degree of deformation from the spherical symmetry due to rotation depends on the ratio of the rotation rate and the critical rotation of the star ($\Omega /{\Omega }_{\mathrm{crit}}$ ). The critical rotation is reached when the centrifugal force equals the gravitational force in the equatorial plane.

For this run, we will create a slow rotating model with $(\Omega /{\Omega }_{\mathrm{crit}}{)}_{\mathrm{initial}}$ = 0.2. Since our PMS model did not include rotation, you need to include additional steps for the numerical scheme to relax to the set initial rotation rate when including rotation in the ZAMS.

new_omega_div_omega_crit = 0.2d0
set_near_zams_omega_div_omega_crit_steps = 20

controls

Next, in the &controls section add a stopping criterion when the central hydrogen abundance drops below $10^{-3}$ , i.e. at the end of the main-sequence stage.

xa_central_lower_limit_species(1) = 'h1'
xa_central_lower_limit(1) = 1d-3

As we introduced in the beginning of this lab, the turbulent viscosity coefficient ${\nu }_{\mathrm{AM}}$ determines the efficiency of the transport of AM by diffusion. A very high value of ${\nu }_{\mathrm{AM}}$ induces very efficient AM transport and results in near solid body rotation rate as it is the case in convective regions.

For the purpose of this lab we will include an arbitrary viscosity coefficient ${\nu }_{\mathrm{AM}}=10^5{\mathrm{cm}}^2/s$ that is uniform throughout the star.

set_uniform_am_nu_non_rot = .true. 
uniform_am_nu_non_rot = 1d5 !cm^2/s

Lastly, modify the output LOGS directory name to specify the physics you included in your model, for example LOGS_0.2_nuvisc_1d5. Do not forget to change the name of the LOGS directory in each run otherwise the output files will be overwritten.

log_directory = 'LOGS_<omega/omega_crit>_nuvisc_<am_nu_value>'

history_columns.list

Now that we have included the rotation parameters in our inlist let’s add the rotation variables to the history.data file so that we can track how those parameters change along evolution. To do so, you just need to modify the file history_columns.list located in your work directory.

• surf_avg_omega
• surf_avg_omega_div_omega_crit
• grav_dark_L_polar
• grav_dark_Teff_polar
• grav_dark_L_equatorial
• grav_dark_Teff_equatorial

profile_columns.list

Let’s also add the rotation profile information to the profile.data files by modifying the file profile_columns.list in your work directory.

• omega
• omega_div_omega_crit

MESA run

Now that we have included the relevant physics in our inlist let’s start the computation. Do not forget to always clean the executable files, compile the code and, run the executable file generated. The run will stop automatically after it reaches the stopping criterion (this might take some minutes < 5min).

You should be looking at the following plots:

• Top left: evolution of surface rotation rate divided by the critical rotation rate (left axis) and the ratio between the near core and surface rotation rate (right axis);
• Two small panels bottom left: HR diagram as seen from different angles. The left panel shows the intrinsic variables, the right panel shows the polar and equatorial variables;
• Top right: the mixing panel (already seen on Monday);
• Two panels bottom right: the rotation rate against mass and radius.

Don’t worry if your run has finished before you grasped the content of the plots. The pgplot of the final model was saved in the /png folder under the name grid1000200.png.

At the end of your run your pgstar plot should look like this.

The effect of rotation on the surface parameters

We will now look at the effect of rotation in the global surface parameters. To do so we need to increase the rotational velocity in our models.

inlist_project

Let’s modify our inlist_project to create another model with the same physics as our previous model but with faster rotation.

new_omega_div_omega_crit = 0.6d0

log_directory = 'LOGS_0.6_nuvisc_1d5'

MESA run

Let’s run MESA again. The run will stop automatically after it reaches the stopping criterion (this might take some minutes < 5min). After the run has finished do not forget to modify the png file name, for e.g. grid1000495_0.6_nuvisc1d5.png.

(Magneto)-hydrodynamic instabilities

In the previous task we took a very simple approach by including a constant ad-hoc viscosity in MESA models. Now, we will take a more physically motivated approach by including (magneto)-hydrodynamic instabilities in our stellar models.

In addition to the effects of rotation on the surface parameters, rotation also triggers several hydrodynamical instabilities that can transport AM and chemical elements in the radiative regions. MESA includes several (magneto)-hydrodynamical (MHD) processes in its AM transport prescription:

• dynamical shear instability (DSI),
• Solberg-Høiland instability (SH),
• secular shear instability (SSI),
• Eddington-Sweet circulation (ES),
• Goldreich-Schubert-Fricke instability (GSF),
• Spruit-Taylor dynamo (ST).

See more details on these physical processes in Heger et al. (2000).

In this exercise we will compute models with different $(\Omega /{\Omega }_{\mathrm{crit}}{)}_{\mathrm{initial}}$ and different rotational-instabilities.

inlist_project

Let’s modify inlist_project according to what you have chosen in the google spreadsheet. First, we need to disable the uniform viscosity in the &controls section of the inlist_project that we used in the previous exercise.

set_uniform_am_nu_non_rot = .false. 
uniform_am_nu_non_rot = 0d0

Next, we will include the AM transport by (magneto)-hydrodynamic processes in our run. In MESA, the viscosity coefficient $({\nu }_{\mathrm{AM}})$ is calculated as a sum of the diffusion coefficients for convection, semi-convection and the (magneto)-hydrodynamical instabilities described above.

D_DSI_factor = 0d0   ! dynamical shear instability
D_SH_factor  = 0d0   ! Solberg-Hoiland
D_SSI_factor = 0d0   ! secular shear instability
D_ES_factor  = 0d0   ! Eddington-Sweet circulation
D_GSF_factor = 0d0   ! Goldreich-Schubert-Fricke
D_ST_factor  = 0d0   ! Spruit-Tayler dynamo

Lastly, you need to modify the $(\Omega /{\Omega }_{\mathrm{crit}}{)}_{\mathrm{initial}}$ according to the value you chose in the Google spreadsheet. Remember that we have already done this step in the previous tasks.

new_omega_div_omega_crit = <VALUE>

log_directory = 'LOGS_<VALUE>_<INSTABILITY>'

MESA run

Finally let’s run MESA! The run will stop automatically after it reaches the stopping criterion (this should take some minutes < 5min). After the run has finished do not forget to modify the png file name, for e.g. grid1000495_0.6_DSI.png.

After your run is finished, open the history.data file in your favourite text editor and look in the last line (which corresponds to the last computed model) for the data required to fill in the google spreadsheet.

Comparison with 2D Ester models

The last task in this lab1 is to compare your MESA tracks with 2D Ester models. In this lab we mainly focus on the evolution of the rotation rate and the HR diagram. In the next labs you will be able to explore and test other important physical quantities.

If you are running into problems with the jupyter notebook, you can send your tutor the history.data file of your runs (rename those files according to your initial rotation rate and instability).

Bonus Task: Critical rotation

Some of the runs in Exercise 2 may have led to rotation rates above critical rotation. In order to automatically stop a computation when this takes place you need to add a customized stopping criterion to the run_star_extras.f90 file located in the src directory.

run_star_extras.f90

You can find out which variables you need by searching around in the $MESA_DIR/star_data/public/star_data_step_work.inc file (in particular amongst the surface variables).

As our PMS model did not include rotation, the numerical scheme requires additional steps to relax to the chosen rotation rate (which you already set in the first tasks). During this relaxation process the initial rotation rate can reach above critical values for some models before it relaxes to the set initial rotation rate. Therefore, only check the stopping condition for models above model_number = 5 (e.g. s% model_number > 5).

integer function extras_check_model

if (s% model_number > 5) then
    if (s%omega_avg_surf/s%omega_crit_avg_surf >= 1d0) then
        extras_check_model = terminate
        write(*, *) 'termination code: reached critical rotation'
        return
        end if
end if

MESA run

You can test if your stopping condition works for other combinations of initial rotation rate and instabilities. You might need to increase the minimum model number (s% model_number > 5) in some cases.

Troubleshooting

If you are running into solver problems due to the high rotation rates (e.g. a lot of retries or termination code message), try increasing the mesh size and decreasing the timestep. The values we provide in the inlist_project ensure a fast run but are not necessary the best to ones to ensure convergence of the numerical scheme.

Vary the parameters mesh_delta_coeff and time_delta_coeff according to the specific physics included in your model and always perform convergence testing in your real science cases to ensure that your runs have indeed converged.