Default parameters for relaxation to a radiative core

[warrickball]

A new feature in MESA is a relaxation of a pre-main-sequence model until it forms a radiative core. This is often a more robust point at which to start changing things like atmospheres, rotation rates, etc. But we've steadily had issues where someone stumbled upon this routine, which is currently on by default, behaving unexpectedly. Let's discuss those defaults, which are (in b9a5f96):

    pre_ms_relax_to_start_radiative_core = .true.
    pre_ms_min_steps_before_check_radiative_core = 50
    pre_ms_check_radiative_core_start = 1d-6
    pre_ms_check_radiative_core_stop = 1d-3
    pre_ms_check_radiative_core_Lnuc_div_L_limit = 0.1d0
    pre_ms_check_radiative_core_min_mass = 0.3d0

I'm personally against this being turned on by default. The models where it helps are often already not "vanilla" stellar models and I think this is more valuable as a helpful tool for those users rather than something for relatively plain and simple models.

Over in #340, I started a tangent about changing pre_ms_check_radiative_core_Lnuc_div_L_limit to help that issue but it doesn't obviously do so, even though the current default of 0.1 seems surprisingly large. I need to experiment quite a lot more (especially with 60+ solar mass stars) to decide on a reasonable value.

Changing it to 1d-6 didn't break any tests, so I'm going to try disabling pre_ms_relax_to_start_radiative_core everywhere in my whb/pre-ms-Lnuc-limit branch to see if any of the tests are actually sensitive to this stuff at all, or if it's a feature that isn't actually being tested.

[mjoyceGR]

I support making "off" the default setting for a number of reasons, one being that it is confusing to see information about a radiative core for stars with no radiative cores (or radiation at all) as the first thing being printed in the terminal. This is not a configuration that makes sense for a solar model, and I think it's reasonable to assume that that will be what most PhD advisers or lecturers tell their students to build first.

[fxt44]

i still do not have a clear picture of the "problem". i've run masses between 0.5 and 100 msun starting from the pre-ms in 15140, and do not recall experiencing an issue of a "bad" pre-ms model. please share an inlist and a version number that demonstrates the "problem".

[mjoyceGR]

One problem is explained in detail on the phase_of_evolution discussion https://github.com/MESAHub/mesa/issues/340

To summarize, leaving these on by default causes a solar-like model to exit the relaxation phase at a central temperature of ~10^6.5, when it should start at something like ~10^5. This messes up the pre-main sequence, and sometimes it is necessary to have things like rotation or other hooks "on" well before this point in low-mass models. See, in particular, discussion about using the deuterium birth line as a physically motivated definition of the pre-main sequence as distinct from the post-relaxation "starting phase"

[fxt44]

ok. please share a minimal inlist and version number that shows the problem.

[mjoyceGR]

the commit I'm using is b9a5f96 and inlist below

&star_job

    save_pulse_data_when_terminate = .false.
    !  save_pulse_data_filename = 'pulse_..._seismic.data'

    ! begin with a pre-main sequence model
    create_pre_main_sequence_model = .true.
    pre_ms_relax_to_start_radiative_core = .true. !.true. ! by default, but it should not be!!! 

    ! save a model at the end of the run
    save_model_when_terminate = .false.
    save_model_filename = 'model.mod'

     new_rotation_flag = .false.

    pre_ms_T_c = 3d5 ! default -- any lower doesn't work
    pre_ms_check_radiative_core_min_mass = 0.3d0 ! default is 0.3d0
    pre_ms_check_radiative_core_Lnuc_div_L_limit = 1d-6

    set_uniform_initial_composition = .true.

    initial_h1 = 0.7409!0.7        ! -1
    initial_h2 = 2.75d-5  !2.75e-5    ! Lodders 2009
    initial_he3 = 8.5d-5    ! Lodders 2009 
    initial_he4 = 0.24 !0.28      !! initial Y!-1

      ! + ``AGSS09_zfracs = 6``
    initial_zfracs = 6 !! scaling of the rest of the elements according to someone's solar abundance

    ! display on-screen plots
    pgstar_flag = .false.

    !use jina
    set_rates_preference = .true.
    new_rates_preference = 2

    !network 
    change_net = .true.
    new_net_name= 'pp_extras.net' !'basic.net' !'pp_and_cno_extras.net'  

/ !end of star_job namelist

&eos
/

&kap
    Zbase = 0.0140 

    !opacities with AGSS09 abundances
    kap_file_prefix = 'a09' 
    kap_lowT_prefix = 'lowT_fa05_a09p'
    kap_CO_prefix = 'a09_co'
/

&controls
    atm_option = 'table' !'photosphere_tables' 
    atm_table = 'photosphere'
    !max_backups_in_a_row = 100 !15

    warn_when_large_rel_run_E_err = 99d0 !0.02d0
    star_history_name = 'history.data'
    log_directory = 'LOGS'

!------------------------------------------------------------------
    !initial_z = 0.0140 !0.014
    !initial_y = 0.24

    max_num_profile_models = -1
    history_interval = 1
    terminal_interval = 100
    write_header_frequency = 10
    photo_digits = 5
    photo_interval = 100

    ! starting specifications
    initial_mass = 1.00  ! in Msun units
    min_timestep_limit=0
    max_model_number = -1

    stop_at_phase_PreMS = .false.
    stop_at_phase_ZAMS = .false.
    stop_at_phase_IAMS = .false.
    stop_at_phase_TAMS = .true.

    cool_wind_RGB_scheme = 'Reimers'
    cool_wind_AGB_scheme = 'Blocker'
    RGB_to_AGB_wind_switch = 1d-4
    Reimers_scaling_factor = 0.1
    Blocker_scaling_factor = 0.2
    max_wind = 1d-3

    mixing_length_alpha = 1.9500  
    MLT_option = 'Henyey'

    do_element_diffusion = .false.
    report_solver_progress = .false.
    ! Timestep and grid

    mesh_delta_coeff = 0.8
    time_delta_coeff = 1.0 

    varcontrol_target = 1e-4 
    !max_allowed_nz = 50000

    !to help with convergence
    okay_to_reduce_gradT_excess = .true.
    Pextra_factor = 2.0

!-------------------------------------------
!! convective overshoot, new formalism
!-------------------------------------------
    overshoot_scheme(:)    = 'exponential' ! ``exponential``, ``step``, ``other``
    overshoot_zone_type(:) = 'any' !  ``burn_H``, ``burn_He``, ``burn_Z``, ``nonburn``, ``any``
    overshoot_zone_loc(:)  = 'any' ! ``core``, ``shell``, ``any``
    overshoot_bdy_loc(:)   = 'any' ! ``bottom``, ``top``, ``any``

    overshoot_f(:)  = 0.014d0
    overshoot_f0(:) = 0.004d0

    ! Convergence checking
    use_gold_tolerances = .true.

    max_number_retries = -1 
    max_model_number = 10000   

    !better resolution of the Henyey hook
    delta_lg_XH_cntr_max = -1

/ ! end of controls namelist

[warrickball]

This is a more limited set of changes rooted in @mjoyceGR's example. It crashes for me on b9a5f96 just after the end of the relax_to_radiative part:

&star_job
    create_pre_main_sequence_model = .true.

    !network                                                                                                                                                                                                                                                  
    change_net = .true.
    new_net_name= 'pp_and_cno_extras.net' !'basic.net' !'pp_and_cno_extras.net'                                                                                                                                                                               
/ ! end of star_job namelist                                                                                                                                                                                                                                  

&eos
/ ! end of eos namelist                                                                                                                                                                                                                                       

&kap
  use_Type2_opacities = .true.
  Zbase = 0.02
/ ! end of kap namelist                                                                                                                                                                                                                                       

&controls
  ! see star/defaults/controls.defaults                                                                                                                                                                                                                       

  ! starting specifications                                                                                                                                                                                                                                   
    initial_mass = 0.75 ! in Msun units
    initial_y = 0.28
    initial_z = 0.02

  ! options for energy conservation (see MESA V, Section 3)                                                                                                                                                                                                   
    energy_eqn_option = 'dedt'
    use_gold_tolerances = .true.

    ! stop when the center mass fraction of h1 drops below this limit                                                                                                                                                                                         
    xa_central_lower_limit_species(1) = 'h1'
    xa_central_lower_limit(1) = 0.6
/

My initial objections were less about particular failures and more about whether it's too much of a departure from "vanilla" modelling. The main thing that gets me is when I start a run and ask for a "pre-main-sequence" model, I expect something up the Hayashi track with the temperature everywhere too low for nuclear reactions, and that's not what I get. Maybe that's just bred of familiarity of what MESA used to do. (Change!?)

However!

After a day of rumination and an hour or so of investigation, I'm inclined to leave this overall relax_to_radiative phase on, document it more thoroughly, and leave this issue as a place to experiment with whether we can find defaults that suit more people better. Here's why this is my position.

1. Above all, I forgot that this was already present in r15140 and we haven't seen a flood of messages on the mailing list, so it must be working reasonably well.
2. Aside from these crashes (which I'm hoping to look into further), @fxt44 is right that I can't generally find an issue. I have seen some early spikes while catalytic abundances come into equilibrium but I don't think this is any worse than using homogeneous ZAMS models.
3. I think it's fair to argue that if you're interested in evolution before ¹H fuses all the way to ⁴He (which many write off as "pre-main-sequence evolution"), then you probably want to be careful about how the pre-main-sequence evolution is computed. But I do see the usefulness of the extra relaxation.
4. Though pre_ms_check_radiative_core_Lnuc_div_L_limit = 0.1d0 seems large enough that significant abundance changes would have started to happen (the defaults imply 0.9d0 is near the MS), I haven't actually seen ¹H drop by more than about 1d-7 . I think this is partly because the timesteps are still dictated more by the thermal than nuclear timescale, so even if the nuclear luminosity seems to be becoming appreciable, there still hasn't been enough time to transform much.

So I now lean towards adding some more detailed documentation about this, in particular indicating in create_pre_ms_model that relax_to_radiative_core is on by default, and that you probably want to check those controls. The code is already quite loud about what it's doing.

I'll also add a note about the interaction between set_uniform_composition and initial_Z/Y. If starting from a pre-MS model, the former works better if initial_Z/Y are close to the desired uniform composition. (A separate issue is why there isn't a matching relax_uniform_initial_composition option. There is a relax_initial_composition but that reads from a file.)

[fxt44]

thanks meridith.

[fxt44]

i'll assume that reproducing fig 15 of mesa I with a near current version of mesa will suffice. if yes, then steve kawaler's 2013 and marc pinsonneault's 2019 mesa summer school material appear relevant. if no, then i'm misunderstanding something.

[mjoyceGR]

see my post on #340 about the defaults for initial core temperature...it is higher than the temperature we are currently testing against for setting phase_PreMS

[fxt44]

attached is a redo of mesa I's fig 15 with the current head and appended is my inlist. deuterium and lithium-7 burn lines are present. in my experience, it is best to do the initial relaxation with pure vanilla mesa - no overshoot, no mass loss, etc - then add whatever one wants.

! low mass single star pre-main sequence deuterium and lithium depletion demonstration ; see fig 15 of mesa I

&star_job ! see star/defaults/star_job.defaults

! begin with a pre-main sequence model create_pre_main_sequence_model = .true. pre_ms_relax_to_start_radiative_core = .false.

! save a model at the end of the run save_model_when_terminate = .false. save_model_filename = 'final.mod'

! initial composition set_uniform_initial_composition = .true.

initial_h1 = 0.7409 
initial_h2 = 2.75d-5    ! Lodders 2009
initial_he3 = 8.5d-5    ! Lodders 2009 
initial_he4 = 0.24      ! initial Y 

! metal distribution; AGSS09_zfracs = 6 initial_zfracs = 6

! choose network change_net = .true. new_net_name= 'pp_extras.net' !'basic.net' !'pp_and_cno_extras.net'

! display on-screen plots pgstar_flag = .true. save_pgstar_files_when_terminate = .true.

/ ! end of star_job namelist

&eos ! eos options ! see eos/defaults/eos.defaults

/ ! end of eos namelist

&kap ! kap options ! see kap/defaults/kap.defaults

use_Type2_opacities = .true. Zbase = 0.0140

! opacities with AGSS09 abundances kap_file_prefix = 'a09' kap_lowT_prefix = 'lowT_fa05_a09p' kap_CO_prefix = 'a09_co'

/ ! end of kap namelist

&controls ! see star/defaults/controls.defaults

! starting specifications initial_mass = 1.0 ! in Msun units

! when to stop ! stop when the star nears ZAMS (Lnuc/L > 0.99) Lnuc_div_L_zams_limit = 0.99d0 stop_near_zams = .true.

! wind

! atmosphere

! rotation

! element diffusion

! mlt mixing_length_alpha = 1.9500
MLT_option = 'Henyey'

! mixing

! timestep resolution time_delta_coeff = 1.0 max_model_number = 10000

! mass resolution mesh_delta_coeff = 1.0 max_dq = 1d-3 ! default 1e-2, 1/max_dq = minumum number of cells

! solver ! options for energy conservation (see MESA V, Section 3) energy_eqn_option = 'dedt' use_gold_tolerances = .true.

varcontrol_target = 1e-4 
max_allowed_nz = 50000

! output report_solver_progress = .false. profile_interval = 100
max_num_profile_models = 4000 history_interval = 1 terminal_interval = 50 write_header_frequency = 4 ! mod(model_number, write_header_frequency*terminal_interval) = 0 photo_digits = 8 photo_interval = 100

/ ! end of controls namelist

[warrickball]

I've just been experimenting with this amid some new work and I think we can close this for the reasons I raised above. If anyone disagrees, you're welcome to re-open the issue but we need to demonstrate that there is actually an issue here. From my experiments, the parameters for relax_to_radiative_core are already quite well chosen such that the post-relaxation model is around the bottom of the Hayashi track, where the surface temperature starts increasing (the start of the Henyey track), and that this is broadly useful enough for most users to merit the minor inconvenience for people interested in the details of pre-MS evolution.