Context Navigation

← Previous Post
Next Post →

Modified giant profiles (continued from last blog)

I continued to work on generating the modified RGB profile to use in AstroBear, as discussed in the previous blog.

As I had been confused by the mismatch between the modified Lane-Emden (MLE) solution and the MESA solution at the transition radius r=h, I wrote S. Ohlmann to ask about this. He explained that the pressure profile of the MLE solution is not actually used. Instead the equations for dm/dr and dp/dr are integrated again, this time backwards starting from some point on the MESA profile, to obtain dp/dr and p(r) for r<h. Thus, only the MLE profile for \rho(r) is used. (So the final profile for r<h will not be a true polytrope, but that is okay.) This made sense but then at least the MLE profile for dp/dr should almost match the MESA profile at r=h, I realized, because both solutions satisfy hydrostatic equilibrium, after all, and by construction m(h) is almost the same for both solutions. The reason it is almost the same and not exactly the same will be discussed below.

This led me to look at my code and sure enough I discovered that the solution was not properly converging. It was oscillating between two apparent solutions, which led to the error in the dp/dr profile. Therefore I set on fixing the PDE solver. The method I had been using was a shooting method (of my own design). I consulted Numerical Recipes (Press et al.) and found the relevant programs, but eventually realized I would have to start from scratch if I did it that way. I modified my own program and it now converges properly, iterating over 2 free parameters to satisfy the two boundary conditions at r=h. (It could instead be solved as a boundary value problem but that is more complicated and unnecessary.) A useful "trick" is to first run the code in low resolution, which converges in a few seconds, and then you have a good estimate of the parameters which you can use to run in high resolution.

As can be seen from this new version of profile.pdf, the profiles are reasonable and match those of Ohlmann. The MLE profile for p(r) does not precisely match MESA at r=h (as explained last blog there is no reason it should) but is now only off by about 20% for the n=3 polytrope solution with h=0.05R. This seems very reasonable. The only concern I had (which I kept revisiting) is that dp/dr obtained from the hydrostatic equilibrium equation by first integrating \rho(r) to obtain m(r), does not perfectly match MESA at the boundary. After studying various solutions (with n=1.5 or h=0.1R or h=0.01R) I believe this is just because the IDL differentiation and smoothing functions used to obtain dp/dr from the MESA p(r) profile is imprecise and sensitive to the number of points chosen for the smoothing, etc. I don't think this will cause problems in AstroBear. Note that the d\rho/dr profile matches perfectly at the boundary. But this is by design: IDL was used to calculate d\rho/dr which was then used to set the boundary condition that d\rho/dr|_h should equal the MESA value. Note that MESA apparently does not give the option of directly outputting dp/dr.

There is, however, one slight inconsistency in the method is the following (as alluded to above). The value of m(h) of the MESA profile is equal to the point particle mass m_c. However, after solving the MLE, we have m(h)=m_c+m_{MLE}(h), where m_{MLE}(r) is the gas profile obtained by solving the MLE. By integrating backward and resetting m(h)=m_c, a discontinuity is avoided. But this means that the core at r<h is being replaced in AstroBear by a mass m_c+m_{MLE}(h) that is greater than the MESA value of m(h)=m_c. In other words, we have added the extra mass m_{MLE}(h) to the core, so the MESA profile for r>h will no longer be in perfect hydrostatic equilibrium. I discussed this issue with Ohlmann who was aware of it. He says he tried setting the particle mass m_c to some value less than m_{MESA}(h) and then iterating over m_c (yet a third parameter over which to iterate). He said it did not make enough of a difference to worry about. Tomorrow morning I will try a manual iteration and report back here.

As promised I have now tried the iteration procedure mentioned. This was done by modifying the code to include an additional parameter: m_c which is not equal to m_{MESA}(h)=m(h). I obtain m(h) in the usual way from MESA but then instead of setting m_c to equal m(h), I set m_c=f m(h), with f a fraction less than (but close to) unity. I first set f to a reasonable value 0.96 given that for f=1, we had m(h)\simeq 1.04m_{MESA}(h). After a few manual iterations, I was able to get a converged solution such that m(h)=m_c+m_{MLE}(h)=m_{MESA}(h). The value of f turns out to be 0.963 for our fiducial profile (n=3 polytropic index and h=0.05R). I have now added the plots to the end of profile.pdf (Figs. 14-16). I have compared these plots against those where the iteration over m_c is not performed (Figs. 2-5). The profiles are almost indistinguishable. However, the new method including this manual iteration is more self-consistent. It is not difficult to implement nor does it require significantly extra time. Therefore I see no reason why it should not be used when constructing the profiles to be used in AstroBear. Including this last step, not included by Ohlmann, will probably not typically make an important difference. Nevertheless, it is an improvement that I think is worth noting in the paper.

astrobear

Context Navigation

Modified giant profiles (continued from last blog)

Comments

Download in other formats: