[mesa-users] Rapidly accreting metal-free protostars

Tue Apr 22 03:25:25 EDT 2014

Hi again,

Following on a stab at this problem about 6 months ago (to which Bill 
replied), I'm now part-supervising a Bachelor's student who's been 
assigned to the problem of rapidly accreting supermassive protostars. 
Once again, I'm looking for any community experience that can help us 
here.

Our basic strategy was to start with low accretion rates (< 10^-6 Msun/yr) 
and gradually crank it up to a breaking point.  We seem to be getting 
somewhere and are now focusing on trying to get the code running longer 
and faster.

Attached is an example inlist for a 10Msun star accreting at 10^-3 
Msun/yr.  This run eventually halts with a central hydrogen abundance 
around 0.3, which seemed strange to me.  I'm used to stars being quite 
happy once on the main sequence.

A closer look shows that, a the start of the main sequence, the convective 
core boundary moves outwards.  It then starts getting a bit flaky a 
wobbling around.  At this point, the surface radius starts to show what 
looks like a numerical instability.  Meanwhile, inside, the convective 
core boundary is bouncing around a bit.

My problem in approaching this is I'm not sure how to disentangle causes 
and effects.  Various amounts of skullduggery have made me able to push 
this model much further (a bit past core He exhaustion), but I'm not 
satisfied with that as a solution.  A large amount of it is in limiting 
both how much the timestep grows and how much it can decrease, both during 
ordinary evolution and on retries and backups.  These changes seem to keep 
the jostling of the core convective zone under some sort of control, which 
stabilizes things a bit, though the code has to work much harder to push 
through.

We've tried various selections of overshooting, semiconvection, 
reaction networks, density instead of pressure as the intrinsic variable, 
equations of state, mesh-delta-coefficients, varcontrols, etc. to no 
avail.  (I thought that overshooting would help here, but it doesn't seem 
to.)

Note that I'm still using a "normal" photospheric boundary condition: no 
shocks or gravitational radiation from the infalling gas.  I realize that 
this is fairly unphysical, but given that the code's struggles seem 
unconnected to subtler choices in input physics, I don't *think* that's 
connected to why the code is struggling.  (But maybe something like 
correctly setting the specific entropy of the infalling gas would 
stabilize the situation?  I'm out of my depth on such questions.)

The only physical effect that I can see bearing on the problem is the way 
that nuclear burning proceeds.  For those unfamiliar with massive Pop III 
stars, there is at first no CNO to catalyse the CNO-cycles, so the star 
initially burns only through pp-chains.  But this isn't enough to reach 
equilibrium before the core becomes hot enough to start the triple-alpha 
process.  As a result, a small amount of helium is burned into C (or O) 
and then catalysed into a very low metallicity, in C, N and O (about 
10^-8).  In this star, this sensitivity to the amount of CNO available and 
the rapid mixing in convective zones seems to create a step-like profile 
in both the abundances and the nuclear energy generation rate.  None of 
this is helped by the sudden outward jumps of the convective core 
boundary, which mixes in fresh, unburnt hydrogen.

I've attached 2 profiles from the last model of the run.  The first plot 
shows the total energy generation rate and hydrogen abundance as a 
function of temperature.  The hydrogen abundance is a bit of a mess at 
this point.  The second profile is dq as a function of the mass 
co-ordinate (in Msun).  To me, it looks pretty flaky in the interior.  Is 
this symptomatic, causal, or expected?  We haven't played with things like 
the mesh-spacing function or re-meshing criteria, so I guess that those 
may yet help.

As ever, any and all input is greatly appreciated!  If there's anything 
else that'd be useful, just let me know.  We think we've already been 
able to reproduce a few interesting effects presented by other groups 
(notably Hosokawa et al.), so it'd be great to be able to flesh these 
models out a bit.

Cheers,
Warrick

On Fri, 18 Oct 2013, Warrick Ball wrote:

> Hi MESA users,
>
> I'm interested in trying to reproduce some recent results by Hosokawa, Yorke, 
> Inayoshi, Omukai and Yoshida, where they model the stellar evolution of 
> *very* rapidly accreting (up to 1Msun/yr) Population III (i.e. metal-free) 
> protostars, with initial masses around 1 or 2 Msun.
>
> As a first try, I though I'd try to have MESA model the fiducial case of 
> Hosokawa, Omukai & Yorke (2012).  This is a 2 Msun protostar with an initial 
> radius of 177.8 Rsun, accreting at 0.1 Msun/yr.  I set this up with the 
> attached inlist and gave it a try.  (This is MESA 5456, compiled with the 
> April 2013 version of the MESA SDK, running on Scientific Linux 6.4 
> "Carbon".)
>
> And before I discuss the output, I acknowledge that there are a few 
> differences between my MESA input and the Hosokawa et al. models. Notably, 
> the surface boundary condition is different: Hosokawa et al. implemented a 
> shock boundary condition and also solve for the structure of the material 
> falling onto the protostar.  I also noted that the initial model isn't 
> exactly the same: mine starts with initial radius of 181 Rsun, but that 
> should be close enough.
>
> Attached is a (basic) plot showing the evolution of 3 timescales and the 
> timestep (in years) as a function of the total stellar mass.  The three are 
> the dynamical (t_dyn), Kelvin-Helmholtz (t_KH) and accretion (t_acc) 
> timescales.  Also shown is the accretion timescale divided by the number of 
> points in the models (usually around 600).  The accretion rate is constant at 
> 0.1 Msun/yr, so the age is basically 10 times the mass.  The plot stops 
> roughly at the onset of hydrogen burning, when the central temperature has 
> reached about 120 MK.
>
> In short, my crude setup runs, but it runs *very* slowly.  The timestep gets 
> up to about a month and then declines to a few days, ultimately not much more 
> than the dynamical timescale.  It stutters a lot around the onset of hydrogen 
> burning.  (I haven't plotted that far because this test run ran out of 
> cluster computing time and the next run is still waiting in the queue.)  Does 
> the community have any tips about modelling Pop III and/or rapidly-accreting 
> protostars?  I had a look around the literature for papers with Pop III MESA 
> models but didn't find anything that I thought might steer me toward a faster 
> calculation.
>
> I tried to follow the recent discussion about accretion and, if I understand 
> correctly, I shouldn't get too hopeful because the results are only really 
> reliable as long as Mdot*dt < dm(surf). i.e. the while the mass accreted in a 
> timestep is less than the mass in the outermost shell. My crude average is a 
> start, but the outermost shells are almost certainly less massive (they are 
> in other codes) and thus the constraint is tighter than my plot.
>
> In fact, I probably ought to be careful that the outer-mass-shell accretion 
> timescale isn't becoming shorter than the *dynamical* timescale...
>
> Cheers,
> Warrick
>
>
> Hosokawa, T., Omukai, K., & Yorke, H.W. 2012, ApJ, 756, 93 
> http://adsabs.harvard.edu/abs/2012ApJ...756...93H
>
>
> ------------
> Warrick Ball
> Postdoc, Institut für Astrophysik Göttingen
> wball at astro.physik.uni-goettingen.de
> +49 (0) 551 39 5069

------------
Warrick Ball
Postdoc, Institut für Astrophysik Göttingen
wball at astro.physik.uni-goettingen.de
+49 (0) 551 39 5069
-------------- next part --------------
! iinlist_project 

&star_job

      show_log_description_at_start = .false. 

      mesa_dir = '' ! empty string defaults to environment variable MESA_DIR

      change_lnPgas_flag = .true.

      new_lnPgas_flag = .true.

      eos_file_prefix = 'mesa'

      kappa_file_prefix = 'gs98'

      pgstar_flag = .false. ! .true.

      create_pre_main_sequence_model = .true.

      pre_ms_T_c = 3e5 ! 3e5, helps with convergence, it seems

      initial_model_relax_num_steps = 200

/ !end of star_job

&controls

  ! starting specifications and global/fixed parameters

      !use_Henyey_MLT = .true. ! .false.

      !MLT_option = 'Henyey' ! options are ML1, ML2, Mihalas, Henyey (default)

      which_atm_option = 'simple_photosphere' ! 'Eddington_grey' ! or 'simple_photosphere' or 'Krishna-Swamy' (needs different alpha_MLT)

      mesh_delta_coeff = 1.5 ! 1.0

      varcontrol_target = 1d-2 ! 1d-4

      terminal_interval = 10

      write_header_frequency = 10

  ! model parameters

      initial_mass = 10.0 ! in Msun units

      initial_Z = 0.0

      initial_Y = 0.25

      mixing_length_alpha = 2

  ! when to stop

      max_age = 4.95d8 ! in years

      mass_change = 1d-3

  ! output to files and terminal

      photostep = -1

      write_profiles_flag = .true.

      profile_interval = 100

  ! some numerical controls

      max_iter_for_resid_tol1 = 8 !

      newton_iterations_limit = 20

      refine_solution = .true.

      refine_mtx_solution = .true.

/ ! end of controls

&pgstar

/ ! end of pgstar namelist

-------------- next part --------------
A non-text attachment was scrubbed...
Name: final_profile_eps_nuc_and_X.png
Type: image/png
Size: 28285 bytes
Desc: 
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20140422/debed275/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: final_profile_dq.png
Type: image/png
Size: 51189 bytes
Desc: 
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20140422/debed275/attachment-0001.png>