RGB star with sink particle

1) Understand how MESA stellar profiles translate to the grid.
It is important to make sure that the input from MESA is being gridded correctly, so compare with Ohlmann+16c ​https://arxiv.org/abs/1612.00008v1. See last pages of ​this updated presentation. We see that the core gets cut out due to lack of resolution at the center, as expected. Though the units used are different from the Ohlman plots, I tried to lign them up. It can be seen that the initial pressure, density and specific internal energy profiles are in agreement with Ohlmann, as expected. This confirms the following:
(i) (Slide 5 & 6): The MESA density and pressure profiles that we put in are basically what we get out, except for the very center, which gets washed out due to lack of resolution (the resolution is slightly larger than 1 R_sun, which is where the profiles become flat). These profiles also agree with those of Ohlmann;
(ii) (Slide 7): A cutoff radius of 1 R_sun corresponds to a missing central core mass of about 0.4 Msun (which is the core mass used by Ohlmann);
(iii) (Slide 8): The profile for the specific internal energy looks very similar to Ohlmann's profile (though the units are different). This confirms that our code is using the same (ideal gas) equation of state as Ohlmann.

Now that this sanity check has been done, we can try introducing a sink particle at the center.

2) Introduce sink particle to replace the core. Here I've tried sink particles of the following masses:
a) 0 Msun,
b) 0.2 Msun,
c) 0.4 Msun (closest to the actual core mass),
d) 0.8 Msun,
e) 1.7 Msun,
f) 3.4 Msun,
g) 10 Msun.

Movies of density (both 2D slices and 1D profiles along x) are ​here. Simulations have resolution 1283 in a box length 140 R_sun.

Note that the star is not stable for any of the runs, but the 0.4 Msun sink particle does at least seem to help to stabilize the star compared with the case of no sink particle (see remarks below).

For the 0 Msun case (no sink particle) I've checked that the initial local sound-crossing time t_s(r) = r/c_s(r) is comparable to the time it takes for the outgoing pressure disturbance to propagate to radius r, for a few different values of r. Actually, t_s(r) is a few times larger than this expansion time, but within a factor of order unity, as expected. A plot of the initial sound-crossing time for the 0 Msun case is available ​here and can be compared with the expansion times in the movies showing density mentioned above. Movies of the 1D profile of the sound speed for the 0 Msun and 0.4 Msun cases are available ​here.

Another useful quantity is the fractional difference from hydrostatic equilibrium:

hef = | |grad P| - |rho g| |/max(|grad P|,|rho g|).

This has been plotted in the following ​fractional HE movies.

Remarks:
1) A sink particle with mass equal to the core mass helps to preserve the high density at the center, and slows the expansion due to hydrostatic imbalance, but also creates outgoing sound waves eminating from the center;
2) The presence of the sink particle somehow introduces unphysical 4-fold symmetry, followed by complete asymmetry. This seems to be a numerical problem that needs to be addressed;
3) Though it is not computationally optimal for this problem, it would make life simpler to change back to cgs units as the computational units here (rather than length units of R_sun, as in the current simulations). This would allow easier comparison with Ohlmann results and prevent unit conversion errors for these trial runs. I will implement this change in the next set of simulations.

Next step is to solve Ohlmann's modified Lane-Emden equation inside the truncation (cutoff) radius. The choice of truncation radius determines the core mass = m(r_truncation). Then match this solution onto the MESA profile at the truncation radius, as done by Ohlmann (density should be smooth across the transition at r_truncation). I will write a code to do this. Hopefully this will result in a more stable star.

Comments

No comments.