Implementing sub-Kelvin-Helmholtz mass loss in other wind routine
In the previous task of MINILAB3 you evolved a star that losses mass trough winds at a constant rate. However, in reality, things are often more complicated, and the mass loss rate generally depends on the stellar parameters and its external environment.
In the last part of MINILAB2, you computed the Kelvin-Helmholtz mass loss rate as an additional history column. Here we will learn how to implement a fraction of this mass loss rate through other_wind_routine to represent a sub-adiabatic mass loss rate (meaning -- a mass loss rate which the star can adjust its structure to accommodate on a thermal timescale while maintaining its entropy / not undergoing a runaway expansion). A nice expression for this is taking some small fraction of the Kelvin-Helmholtz mass loss rate.
This sub-adiabatic mass loss rate (i.e. MKH >> Mdot ∝ MKH) is often used as a simple representation of the rate at which a star in a binary system could donate mass to its companion / accrete mass from its companion during a dynamically stable mass transfer episode (note that this is true for non-compact accretors, in the case of compact objects the accretion rate could go all the way up to the Eddington limit).
We will again implement this mass loss from the end of core helium burning to the end of core carbon burning, which is a reasonable time for a star to undergo mass transfer, and is also a relevant timeline for some other interesting mass loss processes, such as pulsationally driven superwinds (see e.g. Yoon & Cantiello 2010), arbitrary "instabilities" (see, e.g. Gofman+2020), or potentially additional episodes of eruptive mass loss (see e.g. Cheng+2024) depending on the stellar parameters at the end of core He burning.
We also emphasize that the process of implementing a mass loss rate tied to the Kelvin-Helmholtz timescale (which is a function of M, R, and L) is very similar to the process of implementing many other mass loss rate prescriptions which depend on the star's mass, radius, and luminosity (see, e.g., recent RSG wind rates from Decin+2024 following the work of Beasor+2020,2023, or the slightly-more-complicated double power law rates of Vink+2023 following the work of Yang+2023). If you have extra time at the end of the lab and want those valuable imaginary MESA bonus points, feel free to try to implement your own preferred supergiant wind prescription! Just be careful that it doesn't approach the thermal/Kelvin-Helmholtz mass loss rate, otherwise the star physically won't be able to respond (see the Bonus task below)...
Reminder: The global Kelvin-Helmholtz mass loss rate, which is an estimate of the maximum mass loss rate a star can suffer and still be able to thermally adjust is
If the equation isn't rendering properly, here it is again, and we apologize for the ugly and inconvenient formatting. This equation is too important not to write down in multiple ways:
MdotKH = M/tKH ≈ 6.6e-7 (M/Msun)^-1 (R/Rsun)(L/Lsun) [Msun yr^-1]
1 . Instead of the constant value of w as implemented in Task 1, now, using the above equation (Eq.3.1) and the input variables of other_wind_routine (hint: what's contained in the parentheses when you define the subroutine ?), set
w = ! the output of the above equation in units of solar masses per year
Make sure to check your units! Note that
use const_def
is globally declared at the beginning of the run_star_extras.f90 file, which means that you can use the solar constants that are already stored in MESA in your routine, rather than inserting their values manually. For a list of all the constants see: $MESA_DIR/const/public/const_def.f90.
Hint 1:
Note that all the stellar variables that appear in the equation for the global Kelvin-Helmholtz mass loss rate are conveniently sent as input variables to other_wind_routine!
Partial Solution:
At this stage, w should look like
w = (6.7d-7)*((Msurf/Msun)**(-1))*(Rsurf/Rsun)*(Lsurf/Lsun)
2. Recompile to make sure there are no errors thus far.
3. Set a tunable mass loss fraction constant in inlist_to_end_core_c_burn as follows:
We can set parameters in the &controls section of our inlist to access them later as part of the star structure instead of hard-coding them in the run_star_extras.f90 file. To do that, MESA provides 4 arrays of different types: x_ctrl (for double), x_integer_ctrl (for integer), x_logical_ctrl (for boolean) and x_character_ctrl (for character). This could be convenient in cases we want to make multiple runs automatically changing that value of a pre-factor and explore its effect on the results or our simulations.
3a. Let us set control of double precision in our inlist. Open inlist_to_end_core_c_burn and add
x_ctrl(1) = 0 !replace 0 with your chosen fraction of the KH mass loss rate
to the &controls section, where the value is set to the fraction of the Kelvin Helmholtz mass loss rate you would like to apply to the star. Don’t forget to write your name next to the chosen fraction in the Day2_Massive_lab3 tab of the google spreadsheet.
Note! Runs with higher values of x_ctrl(1) tend to take longer times, so only choose a high value of x_ctrl(1) if you have a strong computer.
3b. Then, in your other_wind_routine in src/run_star_extras.f90, multiply the expression for the Kelvin-Helmholtz mass loss rate you implemented as w by s% x_ctrl(1), where s% is the star pointer, pointing to x_ctrl(1). This is now your scaling factor!
Partial Solution:
At this stage, w should look like
w = s% x_ctrl(1)*(6.7d-7)*((Msurf/Msun)**(-1))*(Rsurf/Rsun)*(Lsurf/Lsun)
Now, try to recompile ( ./mk )
Uh oh! If you did everything right, then something went wrong. Read the error messages in the terminal... what happened?
4. Setting star variables - The recompilation failed due to a missing declaration of the pointer type of s. To be able to use x_ctrl(1) in our other_wind_routine, let's declare the type of the pointer s and initialize it.
a. Declaring a pointer to star_info - Copy the line
type (star_info), pointer :: s
after the declaration of the input variables (id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z).
b. Initializing the pointer - Copy the line
call star_ptr(id, s, ierr)
and add it after the declaration of the output variables (w, ierr) in your other_wind_routine.
5. Re-compile and run your simulation!!
6. Write your results in the Day2_Massive_lab3 tab of the google spreadsheet. Compare your results with your group-mates - how much mass was lost in each case? how did it affect the final radius of the star? What fraction of the Kelvin-Helmholtz mass-loss rate needed to strip the envelope during core carbon burning?
Full Solution:
This is how your inlist_to_end_core_c_burn should look like at this point:
&star_job
show_log_description_at_start = .false.
! You added to start the run by loading a saved model
load_saved_model = .true.
load_model_filename = 'after_core_he_burn.mod'
change_initial_net = .true.
new_net_name = 'co_burn_plus.net'
save_model_when_terminate = .true.
save_model_filename = 'after_core_c_burn.mod'
required_termination_code_string = 'xa_central_lower_limit'
/ ! end of star_job namelist
&eos
/ ! end of eos namelist
&kap
/ ! end of kap namelist
&controls
! You added so you can later access it in run_star_extras.f90
x_ctrl(1) = 0 !replace 0 with your chosen fraction of the KH mass loss rate
xa_central_lower_limit_species(1) = 'c12'
xa_central_lower_limit(1) = 1d-3
! limit max_model_number as part of test_suite
max_model_number = 10000
!max_number_retries = 37
! wind
! atmosphere
! rotation
! mlt
alpha_semiconvection = 0
thermohaline_coeff = 0
! mixing
! timesteps
! mesh
! solver
! output
terminal_show_age_units = 'years'
terminal_show_timestep_units = 'days'
terminal_show_log_dt = .false.
terminal_show_log_age = .false.
! You added to easily find the history values relevant for the lab
num_trace_history_values = 8
trace_history_value_name(1) = 'log_R'
trace_history_value_name(2) = 'star_mass'
trace_history_value_name(3) = 'envelope_mass'
trace_history_value_name(4) = 'total_mass_h1'
trace_history_value_name(5) = 'surface_h1'
trace_history_value_name(6) = 'surface_he4'
trace_history_value_name(7) = 'log_L'
trace_history_value_name(8) = 'log_Teff'
!photo_interval = 10
!profile_interval = 2
!history_interval = 1
!terminal_interval = 1
x_integer_ctrl(1) = 5 ! Inlist part number
/ ! end of controls namelist
&pgstar
! Profile_Panels3_xmin = 0 ! 2.5 ! -101d0
! Profile_Panels3_xmax = -101d0 !
! TRho_Profile_xmin = 2.5
! TRho_Profile_xmax = 7.5
! TRho_Profile_ymin = 8.4
! TRho_Profile_ymax = 9.2
/ ! end of pgstar namelist
This is how your run_star_extras.f90 how should look like at this stage:
module run_star_extras
...
subroutine extras_controls(id, ierr)
integer, intent(in) :: id
integer, intent(out) :: ierr
type (star_info), pointer :: s
ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
s% extras_startup => extras_startup
s% extras_check_model => extras_check_model
s% extras_finish_step => extras_finish_step
s% extras_after_evolve => extras_after_evolve
s% how_many_extra_history_columns => how_many_extra_history_columns
s% data_for_extra_history_columns => data_for_extra_history_columns
s% how_many_extra_profile_columns => how_many_extra_profile_columns
s% data_for_extra_profile_columns => data_for_extra_profile_columns
s% other_wind => other_wind_routine ! You added as a pointer to other_wind
end subroutine subroutine extras_controls
...
! You added to control to make the star lose mass at a constant rate
subroutine other_wind_routine(id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z, w, ierr)
use star_def
integer, intent(in) :: id
real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z !surface values (cgs)
real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0)
integer, intent(out) :: ierr
! Sub KH mass loss rate
w = s% x_ctrl(1)*(6.7d-7)*((Msurf/Msun)**(-1))*(Rsurf/Rsun)*(Lsurf/Lsun)
ierr = 0
end subroutine other_wind_routine
...
end module run_star_extras
Bonus task: Running a stellar model with Kelvin-Helmholtz mass loss
Now, set x_ctrl(1)=1 in your inlist and re-run the simulation. To try and understand the outcome of your simulation by adding the following print statements to other_wind_routine in your run_star_extras.f90 file
print *, "The value of Msurf is", (Msurf/Msun)
print *, "The value of Rsurf is", (Rsurf/Rsun)
print *, "The value of Lsurf is", (Lsurf/Lsun)
print *, "The value of w is", w
Don't forget to recompile!
2. Look at the results of the prints in the terminal, and explain the physical reason for the results of this simulation. Does it crash? When and why (or why not)?
Solution:
This is how your run_star_extras.f90 how should look like at this stage:
module run_star_extras
...
subroutine extras_controls(id, ierr)
integer, intent(in) :: id
integer, intent(out) :: ierr
type (star_info), pointer :: s
ierr = 0
call star_ptr(id, s, ierr)
if (ierr /= 0) return
s% extras_startup => extras_startup
s% extras_check_model => extras_check_model
s% extras_finish_step => extras_finish_step
s% extras_after_evolve => extras_after_evolve
s% how_many_extra_history_columns => how_many_extra_history_columns
s% data_for_extra_history_columns => data_for_extra_history_columns
s% how_many_extra_profile_columns => how_many_extra_profile_columns
s% data_for_extra_profile_columns => data_for_extra_profile_columns
s% other_wind => other_wind_routine ! You added as a pointer to other_wind
end subroutine subroutine extras_controls
...
! You added to control to make the star lose mass at a constant rate
subroutine other_wind_routine(id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z, w, ierr)
use star_def
integer, intent(in) :: id
real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z !surface values (cgs)
real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0)
integer, intent(out) :: ierr
call star_ptr(id, s, ierr)
! Sub KH mass loss rate
w = s% x_ctrl(1)*(6.7d-7)*((Msurf/Msun)**(-1))*(Rsurf/Rsun)*(Lsurf/Lsun)
print *, "The value of Msurf is", (Msurf/Msun)
print *, "The value of Rsurf is", (Rsurf/Rsun)
print *, "The value of Lsurf is", (Lsurf/Lsun)
print *, "The value of w is", w
ierr = 0
end subroutine other_wind_routine
...
end module run_star_extras