Task 2. The local thermal timescale

Congrats, you just modeled a star with a fiducial mixing_length_alpha = 2.0. Now we will try to gain a better understanding on the local properties of our star's convective envelopes by computing the associated local thermal timescale tth of every layer of the star from the mass coordinate m to its surface at mass coordinate M:

where cP(m) is the specific heat at constant pressure at mass m, T(m) is the temperature profile of the star and L is its surface luminosity. I should stress, this is a local property of the star rather than global. Therefore, you would like to instruct MESA to compute a profile quantity, and save it as a custom additional column in your profiles.data output files.

I. Modify the run_star_extras.f90

Find your spot

The way in which you tell MESA to add a profile column is very similar to what you did for the additional history column. Open your run_star_extras.f90 and find the routine data_for_extra_profile_columns. The routine works exactly like the data_for_extra_history_columns, but rather than returning two 1-dimensional arrays called names and vals, it returns a 1-dimensional and a 2-dimensional array with the same names, respectively. You should be seeing something like this:

subroutine data_for_extra_profile_columns(id, n, names, vals, ierr)

...

character (len=maxlen_profile_column_name) :: names(n)

real(dp) :: vals(nz,n)

...

! TO DO

! Computation of local thermal timescale t_th, see Eq.(3)

! of your instructions

names(2) = "local_tt"

vals(:,2) = 0.d0

end subroutine data_for_extra_profile_columns

You see that the definitions of the names and vals arrays are such that you can store n names of variables, and per each n index you can store the value of the variable throughout the entire star, from zone 1 to zone nz. The value of n is, again, set by the routine called how_many_extra_profile_columns, which you can check too. For convenience, we already put there the right count for the extra columns, so don't modify it; however, keep in mind that you would normally have to increase the count as you did for the how_many_extra_history_columns routine.

Compute the local thermal timescale

You see that now the index corresponding to the relevant quantity, local_tt, is zero, vals(:,2) = 0.d0. If you already know what to replace this with, you can try to put in there the equation for tth (yr), Eq. (2). Otherwise, give a look at these hints, and / or the solutions.

Hint 1: How is it called?

You might like to know how to refer to the cP(m) information in the star's info structure s: just s% cp(k) will do. Also, the mass element dm̃ is trivially called dm(k) for zone k. You can find these variable declarations in $MESA_DIR/star_data/public/star_data_step_work.inc. Remember to mind the units too...

Hint 2: How do I "if" in fortran?

In fortran, an if statement has the following structure:

if (2 > 1) then

...

end if

Hint 3: Integrals?!

Your star is divided in a number nz of zones, which is of course discretizing your system. Integrals become then just finite sums, and whether or not this is a good approximation is a problem of good or bad resolution. You will hear about it tomorrow morning. The way to achieve a finite sum in fortran is by means of do loops, the equivalent of for loops. A suitable structure for you can be the following:

! ! Remember that the outermost zone is at k=1 and the

! ! innermost is at k = s% nz

do k=1, s% nz

vals(k,2) =

end do

Solution

subroutine data_for_extra_profile_columns(id, n, names, vals, ierr)

...

! TO DO

! Computation of local thermal timescale t_th, see Eq.(3)

! of your instructions

names(1) = "local_tt"

vals(:,1) = 0.d0

! Remember that the outermost zone is at k=1 and the

! innermost is at k = s% nz

do k=1, s% nz

vals(k,1) = (1/s% L(k))*(s% cp(k)*s% T(k)*s% dm(k))/secyer

if (k > 1) then

vals(k,1) = vals(k-1,1)+(1/s% L(1))*(s% cp(k)*s% T(k)*s% dm(k))/secyer

end if

end do

end subroutine data_for_extra_profile_columns

Try to compile

Modifications to the source code in run_star_extras.f90 for this exercise are done. Try now to clean and recompile MESA, to see if everything you typed makes sense:

# Clean and compile

$ ./clean && ./mk

If the compilation is successful, great work! Don't run your simulation just yet; we still need to set up the inlist_pgstar_minilab2 and pick a different value for αMLT.

II. Visualizing the local thermal timescale

A nice(r) mass coordinate

Today has been all about the envelope of stars. As you know, most of the mass of your star is in the core, where gravity is stronger and density is increasing due to the subsequent burning stages. The higher you move away from the core, the fluffier your envelopes (and their mass distribution) will be. We would like to map the local thermal timescale tth with some details in the outer envelope, therefore we introduce the following mass-coordinate:

You can see that this is really "stretching out" the outer layers (i.e., where the argument of the logarithm is small, m→M. If MESA were not smart, you see that we would have to calculate it as an additional profile column, in the same fashion as for local_tt. Luckily, MESA is in fact very smart and the information on the logarithmic mass coordinate is automatically saved in the profile.data files. Open the profile_columns.list file and look for log10(1-q), where q represents the fractional mass ratio, q = m/M.

You should be seeing something like this:

! profile_columns.list -- determines the contents of star model profiles

...

!xm ! mass exterior to point (Msun units)

!dq ! mass of zone as a fraction of total star mass

logxq ! log10(1-q)

!logxm ! log10(xm)

...

You see that a very eloquent logxq profile quantity is uncommented in the profile_columns.list file: this means that MESA is instructed to save this information in any profile.data file. Uncommenting any other quantity of this file will instruct MESA to add the new quantity too! Keep in mind the name logxq: we will use it in a second.

Modify the inlist_pgstar_minilab2

We are now going to manipulate a little bit the inlist_pgstar_minilab2 in order to show tth along our star's profile in the Profile panel of the pgstar window. Open the inlist in your editor and find something like the snippet here below. You will first have to increase the number of Profile panels from Profile_Panels1_num_panels = 1 to Profile_Panels1_num_panels = 2, since you want to add the tth vs logxq plot:

&pgstar

...

! PROFILE PANELS

! ==============

! ! TO DO

! ! Choose the number of panels you want to show

Profile_Panels1_title = 'Profile'

Profile_Panels1_num_panels = 2

...

/ ! end of pgstar namelist

Immediately below this part, you will find the spot in which you can give the direct instructions to plot your favorite quantity, Profile_Panels1_yaxis_name(2) = 'local_tt', as a function of a custom x-axis coordinate, Profile_Panels1_xaxis_name, which will be applied to and shared by all the plots you will add in the Profile panel. Find the lines here below and uncomment them. Don't forget to edit the control Profile_Panels1_xaxis_name = 'mass' with the name of the wanted x-coordinate, your custom profile column logxq.

&pgstar

...

! ! TO DO

! ! x-axis coordinate "logxq" to better

! ! visualize the envelope properties

Profile_Panels1_xaxis_name = 'mass'

! Profile_Panels1_xaxis_reversed = .true.

! ! TO DO

! ! Local thermal timescale to be visualized

! ! into the second profile panel

! Profile_Panels1_yaxis_name(2) = 'local_tt'

! Profile_Panels1_yaxis_log(2) = .true.

...

/ ! end of pgstar namelist

Notice this instruction: Profile_Panels1_xaxis_reversed = .true.: this makes sure that, once you picked the coordinate logxq, its representation in plots will be as a reversed x-axis, i.e. going from the core (log(1-m/M)=0) to the surface (log(1-m/M)=-∞).

NB: After modifying the inlist_pgstar_minilab2, you don't need to recompile MESA, similarly to when you edit inlist_minilab2. However, inlist_pgstar_minilab2 is rather special: you can edit it while the simulation runs, and it will adjust the window on the fly 😊 Therefore, don't worry if you are not sure about these last modifications. During the run, you will easily be able to verify if something is wrong, and modify it again.

Run an α≠2.0 simulation

Now let's change the value of αMLT from αMLT=2.0 to something different. Specifically, let's now pick either a higher or a lower value. You are basically increasing/decreasing the length over which a convective plume transports flux (since we assume the plume transports energy over a characteristic length scale αMLT times the pressure scale height) -- therefore we will be increasing/decreasing the ability of a convective fluid parcel to lift up (or sink down) energy.

Again, we will make use of the PRELAB_models model bank from the Google Drive folder. The starting model of interest for your simulation is going to be located into a path that follows this notation:
[M]Msun/alpha_[alpha]/after_core_he_burn_[M]Msun_alpha_[alpha].mod

Select a value of αMLT from those available. Coordinate with those at your table to get a good range of values! If your neighbour picked a higher or a lower value for αMLT, you should do the opposite case ☝🏻

Similarly to what you did before, you will have to load a saved model that evolved the star till core Helium depletion with a modified αMLT, i.e. a after_core_he_burn_[M]Msun_alpha_[alpha].mod, where alpha is of course your new value >/< 2.0. Pick one of the available models from the PRELAB_models folder; per each mass, there should be the possibility to pick from multiple different αMLT. Move the model inside the working folder:

# Substitute the picked mass M in place of [M] and the picked αMLT in place of [alpha]

$ cp ../PRELAB_models/[M]Msun/alpha_[alpha]/after_core_he_burn_[M]Msun_alpha_[alpha].mod .

You probably know how it goes now: open your inlist_minilab2 and modify the following namelist items:

load_model_filename = 'after_core_he_burn_[M]Msun_alpha_[alpha].mod'

...

save_model_filename = 'after_core_c_burn_[M]Msun_alpha_[alpha].mod'

...

log_directory = 'LOGS_[M]Msun_alpha_[alpha]'

...

Grid2_file_dir = 'png_[M]Msun_alpha_[alpha]'

Additionally, make sure to change αMLT in your inlist_minilab2 to a higher (lower) value:

mixing_length_alpha = !set to your chosen alpha ≠ 2.0

Done with setting up! Let's run the model and stare a little bit at the pgstar window, which should be slightly different with respect to before: now, all the way to the right, you should be able to see a nice new Profile panel which shows you the logarithm of the local thermal timescale log local_tt as a function of the mass coordinate logxq. Rememeber that its x-axis is reversed, namely it goes from the core (log (1-m/M)=0) to the surface (log (1-m/M)=-∞).

Inspect the profile of log tth in your pgstar window.

As you saw during the theory session, the timescale at which you can radiate the thermal energy from a layer at mass m upwards is decreasing as you approach the surface. This just means that you have less energy budget to radiate at the same surface luminosity L. However, the decrease is not smooth: there's a relatively sharp local decrease. By looking at the plot right above this one, namely the profile of the mixing diffusion coefficient D (cm^2 s^(-1)), can you guess what that bump corresponds to?

Hint: Convection and energy transport...

In general, we expect that there is a part of the sub-surface layers of your star where the stellar interior is unstable to convection, but energy transport via convection is radiatively inefficient. This occurs where the optical depth is low enough that radiation is able to carry significant flux. This results in temperature gradients that are significantly steeper than the adiabatic temperature gradient; sometimes these regions are called superadiabatic regions.

Within these radiatively-inefficient, superadiabatic subsurface layers, the steep temperature gradients result in a relatively sharp decrease in the local thermal timescales with m in the upper parts of these regions, compared to above and below the superadiabatic layers.

One more important note: Superadiabatic regions (in the sense that the temperature gradient is much steeper than the adiabatic temperature gradient) can also exist in the deep interior of the star even at high enough optical depths such that radiative diffusion is unable to carry heat flux out of a convective parcel. In fact, there are multiple notions and criteria for "convective efficiency" in the literature, and this can be a source of much confusion. See Jermyn+2022's Research Note on Convective Efficiency for additional discussion.

Looking again at the mixing diffusion coefficient logD (cm^2 s^(-1)) plot, you can notice a sharp drop and the light blue color indicating that convection is happening where log D is larger than zero. On the other hand, stars on the coolest part of their evolutionary track are expected to have convective envelopes, or at least this is what you would find in textbooks. And indeed, you can see the light blue envelope extending down your star in the Kippenhanhn diagram of your pgstar window. So how is it that you seem to have a radiative envelope towards the surface as log(1-mass/M) goes to ?

Hint: Don't be fooled by the log

The amount of radiative envelope mass on the top of the superadiabatic zone is minuscule, look at the values of log(1-mass/M); you wouldn't be able to resolve it in the pgstar Kippenhahn diagram for sure. In general, you can always expect to find such regions right below the photosphere of a star. Moreover, given the large-scale vigorous convection happening below in the envelope, this radiative region would very likely have substantial fluid motion.

Example Profile panel of the pgstar window for 12 M☉ model at the end of core Helium burning, with local thermal timescale tth vs log(1-m/M) plot included.

Crowdsource and compare

As before, annotate, in the Google spreadsheet of Minilab 2, the values of the final surface luminosity logL(L☉), radius logR(R☉), effective temperature logTeff (K) and Kelvin-Helmholtz timescale tKH (yr). Don't forget to specify also the value of αMLT that you picked.

Compare the values for tKH that you just found with the ones from the run with the fiducial αMLT=2.0. Do you understand what is happening to the surface properties of your star, when you increase (decrease) the efficiency of convective energy transport?

Hint Convection in the surface!

Your star has a large convective envelope, which is all of a sudden is more (or less) efficient in driving out energy radiated from the lower star's layers...

You can find the saved pngs of the pgstar window inside the ./png_[M]Msun_alpha_[alpha] folder in the work directory. Grab the last one and look at the Profile panel plot of tth along your star's profile. Which order(s) of magnitude does the local thermal timescale span? Where do you think it becomes more or less comparable to the global thermal timescale?

Annotate, in the Google spreadsheet, the order of magnitude estimate (from looking at the plot) of tth in correspondence of the top of the convective zone (the value at the kink just after the little bumpy feature; or looking at the plot of D (cm^2 s^(-1)) see where it goes to zero or drops down a lot). For even more imaginary MESA bonus points, you could try to report tth at the top of the convection zone algorithmically as an extra history column.

If you still have time, please check out the Bonus exercises for MINILAB2, which explore the concept of a thermal mass loss rate (the stellar mass divided by the thermal time) -- this will be especially relevant in MINILAB3!

Page updated

Google Sites

Report abuse