I continued working toward getting a giant (RGB or AGB) onto the grid.

Recall that when trying to translate the MESA stellar profile to the grid, the compact core gets washed out due to the grid not being fine enough to adequately resolve the local scale height. Last time, I showed that adding a sink particle with the mass of the unresolved RGB core can help to slow down the (unphysical) expansion of the star. But the star is still unstable. This is at least partly due to the fact that the inputted stellar profile has not been modified to take into account the modification to the potential due to the replacement of the central part of the star by a point particle. As explained in Ohlmann 2016 (phd thesis), the asymptotic behaviour of the pressure gradient dp/dr as r tends to 0 must be correctly reproduced (dp/dr \propto r as r\rightarrow0) to avoid unphysical sound waves eminating from the center. This is achieved in Ohlmann 2016 by truncating the MESA profile at some radius r = h (the higher the resolution, the smaller h may be chosen to be) and replacing the inner part of the profile with the solution to a MODIFIED Lane-Emden equation that takes into account not only the potential of the gas (as in the standard LE equation), but also that of the central point particle. Doing this involves a few steps:

1) Choose a truncation radius h and obtain the interior mass m(h) from the MESA profile. This is the mass of the central point particle m_c = m(h).

2) Modify the potential to account for the central point particle using a softening to prevent it from blowing up as r\rightarrow0. Here we follow Ohlmann 2016 (see also Ohlmann et al. 2016, arXiv:1612.00008).

3) For simplicity, approximate the stellar profile for r < h as a polytrope \rho = \rho_0\theta^n, with the polytropic index n a parameter that Ohlmann sets to n = 3. The pressure p = K\rho^{1+1/n}, with K constant, so for n = 3, we obtain 1+1/n = 4/3, as for a relativistic white dwarf. Assuming mass continuity and hydrostatic equilibrium leads to the modified Lane-Emden equation derived in Ohlmann.

4) Now solve the modified Lane-Emden (MLE) equation by first breaking the equation up into a system of two coupled first order PDEs and then solving the system numerically. I have done this using a 3rd order Runge-Kutta "time"-stepping routine in fortran (Brandenburg 2003). For the simplest case of m_c = 0, the MLE equation reduces to the standard LE equation. I've tested my code for this case against a simpler, 1st order ​code downloaded from the web.

Solving the MLE equation is actually much more challenging than solving the LE equation. This is because aside from n, there are now two additional parameters: the central gas density \rho_0, and the characteristic radius \alpha (K can be written in terms of \alpha). Basically, these parameters must be iterated over until a solution is found which satisfies the BOUNDARY CONDITIONS at the transition radius r = h. (Ohlmann mentions using a non-linear root-finder to accomplish this step, but does not provide details.) I have used a simple "home-built" recipe for the iteration that leads to a converged (rapidly enough) solution.

Following Ohlmann, we use the following BCs:
i) \rho(h) matched to MESA profile
ii) d\rho/dr|_h matched to MESA profile

6) This approximate (polytrope) solution for the modified profile for r < h can then be combined with the MESA profile, so that for r > h, one retains the more realistic MESA profile. Now the modified profile can be compared with that of Ohlmann.

Refer to the file ​profile_blog1.pdf.

Figure 1 is the relevant figure from Ohlmann et al. 2016 with which results may be compared.

Figure 2 shows the mass profile of the RGB star as obtained from MESA, and the chosen transition radius h = 0.05R, where R is the outermost radius of the star. The central mass is about 0.4 M_{\odot}, in agreement with Ohlmann.

Figure 3a shows the density \rho(r) obtained from MESA (thick blue) and from the MLE equation (thin black). We have adopted h = 0.05R and the polytropic index n = 3.
Figure 3b shows the density gradient d\rho/dr. Note that \rho and d\rho/dr of the MLE solution are equal to the MESA values at r = h. In other words, the boundary conditions have been correctly implemented. Note also that the MLE solution for \rho(r) agrees with Ohlmann's solution from Figure 1.

Figure 4a shows the pressure p(r). Note that the MLE solution for p (black line) DOES NOT AGREE with the MESA solution at r = h. However, in Figure 1, Ohlmann obtained agreement. Have I made a mistake? Actually, there is no reason to expect perfect agreement, since the MESA solution is not an n = 3 polytrope; MESA employs a more complicated equation of state known as OPAL. However, we can RE-NORMALIZE the MLE solution for p(r) to make it match the MESA solution at r = h (red line). This re-normalized solution now agrees with Ohlmann's solution of Figure 1. Ohlmann mentions "integrating" the stellar structure equations to obtain p(r), presumably while using the MLE profile for \rho(r). But the relation between p and \rho has already been specified when solving the MLE equation, with the constant K (or \alpha) set to satisfy the BC on smoothness of \rho across the transition radius. Thus, after re-normalizing, hydrostatic balance will no longer be perfect. On the other hand, the pressure is now continuous across the transition radius.
Figure 4b shows the pressure gradient dp/dr. Again, we must re-normalize at r = h. Doing so brings the solution into close agreement with that of Ohlmann, seen in Figure 1. The re-normalization factor is close to but not precisely the same as for p(r).

Figures 5 and 6 show the profiles for \rho, d\rho/dr, p, and dp/dr for the case n=1.5, while Figures 7 and 8 show the same for n=5. The solutions obtained are quite similar to that obtained for the n=3 case.

Figures 9, 10, and 11 show the case h=0.1R. That is, with double the transition radius (and n=3). The MLE solutions can be compared with the green lines in Figure 1 and show good agreement.

Figures 12 and 13 show the fiducial h=0.05R and n=3, but now with the boundary conditions modified so that the solution to the MLE equation has \rho(h) equal to the MESA value (as before) and p(h) equal to the MESA value (instead of d\rho/dr|_h). We see that this results in a solution for which d\rho/dr and dp/dr do not match the MESA solution at the transition radius r=h, as would be expected.

Most of the relevant fortran, IDL and tex files for the above work including a README file can be found ​here.

Remarks:
Matching all quantities \rho, d\rho/dr, p, and dp/dr to the MESA solution at r=h AND satisfying hydrostatic equilibrium is not possible to do when approximating the profile for r<h as a polytrope. This is because there are only 2 free parameters (\rho_0 and \alpha), whereas there are 4 boundary conditions, so the problem is overdetermined. Thus, either smoothness or hydrostatic equilibrium must be sacrificed. To the best of my understanding at least, Ohlmann et al. 2016 chose to sacrifice perfect hydrostatic equilibrium to make p(r) and dp/dr smooth accross the transition radius.

Some natural next steps would be:
1) Plot the quantity |dp/dr +Gm(r)\rho(r)/r^2|/max(|dp/dr|,|Gm(r)\rho(r)/r^2|) for the modified profile to assess the level of agreement with hydrostatic equilibrium. Do the same for the unmodified MESA profile and for the pure polytrope solution of the MLE equation.
2) Map the modified solution to the AstroBear grid and perform a run with the appropriate sink particle mass (\sim0.4M_\odot) to assess the level of stability.