If you made it this far and you still have some time, very well done! In this bonus exercise we will compute the mass loss rates associated with the global Kelvin-Helmholtz timescale, ṀKH, and the local thermal timescale, Ṁth, namely Eq.(4) and Eq.(5). The aim of the exercise is an order of magnitude comparison of the two rates, as these represent a way to understand which limit you can impose on the rate Ṁ during a sudden mass loss episode for the star to be able to thermally adjust globally or locally, respectively.
NB If you think that you don't have that much time left, but you're still curious about how the order of magnitude of the global Kelvin-Helmholtz mass loss rate ṀKH compares to that of the local thermal mass loss rate Ṁth, at any point you can jump to section Sec.(III) of this page down here, and just follow the instructions there.
Let us remind ourselves how the Kelvin-Helmholtz mass loss rate is defined:
As you can see, the quantities that appear in the equation are global properties of your star. Consequently, you can modify the run_star_extras.f90 to compute an additional history column. This is nothing different than what you have done for the Kelvin-Helmholtz timescale.
Go to the routine data_for_extra_history_columns and find the right spot:
subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
...
! TO DO
! BONUS exercise. Kelvin-Helmholtz mass loss rate M_KH
! computation
! names(2) = "M_KH"
! vals(2) = 0d0
end subroutine data_for_extra_history_columns
Uncomment the appropriate lines and write down the calculation of M_KH. Don't forget to check the how_many_extra_history_columns routine as well and increase the count of extra columns to the appropriate number.
Hint: Mind the units, again
Usually, the most illustrative units with which astronomers express a rate of mass loss from a star is M☉ yr^(-1), namely solar masses per year. So you might want to divide your calculation by a Msun factor. And be mindful of the units with which you expressed tKH too; if you have it calculated in years already, your task now is quite easy.
Solutions
integer function how_many_extra_history_columns(id)
...
! TO DO
! Add a column to the count
how_many_extra_history_columns = 2
end function how_many_extra_history_columns
subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
...
! TO DO
! BONUS exercise. Kelvin-Helmholtz mass loss rate M_KH
! computation
names(2) = "M_KH"
vals(2) = -s% m(1)/Msun/vals(1)
end subroutine data_for_extra_history_columns
This is not mandatory, but we suggest you to try and compile to see if everything makes sense so far:
# Clean and compile
$ ./clean && ./mk
Remember that the minimum mass loss rate Ṁth required to strip away the outer layers of a star down to this mass coordinate m faster than these layers can thermally readjust is defined as
Since you will have to make use of a structural property of your star, i.e. the local thermal timescale at each mass coordinate, you will have to operate into the data_for_extra_profile_columns routine in the run_star_extras.f90.
Go to the routine data_for_extra_profile_columns and find the right spot:
subroutine data_for_extra_profile_columns(id, n, names, vals, ierr)
...
! TO DO
! BONUS exercise. Computation of the local thermal mass loss rate
! correspondent to the local thermal timescale
names(2) = "local_tml"
vals(:,2) = 0d0
end subroutine data_for_extra_profile_columns
You can quickly check the how_many_extra_profile_columns routine as well, but it should already have the right number of columns so nothing to be done there.
Hint: Mind the units, again
Same hint as above here: you might want to divide your calculation by a Msun factor. And be mindful of the units with which you expressed tKH too; if you have it calculated in years already, your task now is quite easy.
Solution
subroutine data_for_extra_profile_columns(id, n, names, vals, ierr)
...
! TO DO
! BONUS exercise. Computation of the local thermal mass loss rate
! correspondent to the local thermal timescale
names(2) = "local_tml"
vals(:,2) = 0d0
do k=1, s% nz
vals(k,2) = -(s% m(1)-s% m(k))/vals(k,1)/Msun
end do
end subroutine data_for_extra_profile_columns
Let us get to the gist of this EXTRA exercise. Whether you managed to implement both the Kelvin-Helmholtz and the local thermal mass loss rates, or just one of them, or you just jumped here right from Exercise 2 because you're curious, we got you covered: you can download the definitive solutions to Minilab 2 from the folder 3_BONUS_Mass_loss_rates, substitute them to the ones you have in your work directory and proceed with reading the instructions down here. You did great regardless!
If you managed to implement the mass loss rates ṀKH and Ṁth yourself, you can go ahead to compiling and running a simulation (you don't have to modify mixing_length_alpha nor the inlists). If you didn't manage, or your compilation failed, just use the run_star_extras.f90 and the inlist_pgstar_minilab2 provided in the download button up here, and then compile and run.
# Clean and compile
$ ./clean && ./mk
$ ./rn
If you did everything correctly, your pgstar window should now look different than before: in the History panel, you should see an additional plot showing the evolution of the Kelvin-Helmholtz mass loss rate as a function of model_number. Additionally, in the Profile panels, you should have a plot showing the behavior of the local thermal mass loss rate log local_tml as a function of the mass coordinate logxq. An example of what you should see is here below:
Example pgstar window for 16 M☉ model at the end of core Helium burning; this is with the History panel plot of the Kelvin-Helmholtz mass loss rate incorporated, as well as with the local thermal mass loss rate Profile panel additional plot.
During the run, inspect the window and try to think about the following points:
Look at the profile plot for the local thermal mass loss rate log local_tml. Why this shape? It seems that there is an appreciable rise towards the surface, and then a plateau...
Hint: The rise is right there!
When you were inspecting the behavior of the local thermal timescale along your star's profile, you were seeing something like a decrease towards the surface, as the mass coordinate increased, and the appreciable bump in correspondence of the superadiabatic zone...
Look now at the history plot for the Kelvin-Helmholtz mass loss rate log M_KH. What do you think of its shape? Does it make sense, knowing that ṀKH=-M/tKH? What do you think about its order of magnitude?
Compare the order of magnitude of the Kelvin-Helmholtz mass loss rate log M_KH with the order of magnitude of the local thermal mass loss rate log local_tml throughout the structure of the star. In particular, take a close look into the order of magnitude of log local_tml right above the superadiabatic zone, namely, approximately the plateau value. Is this value higher or lower than that of the global Kelvin-Helmholtz mass loss rate?
Let's say you would want to put a limit, for the star to thermally adjust with the right pace, on the rate Ṁ during a sudden mass loss episode (for example, during the mass transfer from the envelope of a donor star to an accretor in a binary system). Would you impose the rate of mass loss to be limited by the Kelvin-Helmholtz rate, which describes the capability of the star as a whole to react to the mass loss? Or would you say that Ṁ needs to be limited to the local thermal mass loss rate close to the surface of the star, for it to be able to respond to the mass loss in a self-regulating manner? Which condition would be more restrictive?
Interesting MESA thing for the future you
In this MESA Lab we TAs have already done the basic setup for the run_star_extras.f90 file for the students. But keep in mind that, when you normally copy the $MESA_DIR/star/work directory, the run_star_extras.f90 file is pretty empty and you would have to copy the contents of the $MESA_DIR/star/job/standard_run_star_extras.inc file and replace it with the line include 'standard_run_star_extras.inc' .