Putting in the MESA equation of state

Work done

1. Examined MESA RGB model (used as initial condition in our simulation) in more detail, including
   1. Temperature profile
   2. Free electron fraction profile
   3. Mean molecular weight profile
   4. Condition for convective instability as a function of radius
   5. Convective time scale as a function of radius

Results

1. Initial conditions
   1. ​RGB Pressure profile from MESA showing total pressure, gas pressure and radiation pressure. We see that gas pressure dominates outside the softening sphere.
   2. ​RGB Temperature profile from MESA. We see that outside the softening sphere the temperature is predicted accurately by assuming the gas to be ideal. Also, the temperature is higher than that obtained by assuming pure fully ionized hydrogen (mu=m_H/2), and lower than that obtained assuming pure fully neutral hydrogen (mu=m_H), except within 1 Rsun from the surface, where it is higher than both. The ionization temperatures of H and He are also plotted for reference.
   3. ​Graph showing gamma1 (adiabatic index) free electron fraction and mean molecular mass profiles from MESA. We see that gamma1 remains below 5/3 and dips below 4/3 near the surface. This behaviour is roughly predicted by an ​analytic model contained in these notes by G. Glatzmaier that assumes pure hydrogen, using the free electron fraction profile from the MESA simulation, shown in black in the plot. We also see that gamma1=(dlnP/dlnrho)_S is quite similar to the polytropic Gamma1=dlnP/dlnrho calculated for the stellar profile, which tells us that the envelope is probably convective. The electron fraction is fairly constant but dips after about r=30Rsun, which coincides with where mu increases. Finally we plot Del=dlnT/dlnP and the critical value of Del=dlnT/dlnP for (i) an ideal gas with gamma1=5/3, which gives (gamma1-1)/gamma1=2/5, (ii) the MESA simulation neglecting the compositional gradient (Schwarzchild criterion) and (iii) the MESA simulation including the compositional gradient (Ledoux criterion). The plots below zoom in on this region to show more clearly what is happening as per convective stability/instability.
   4. ​Graph showing condition for convective instability. This graph should be compared with ​Figure 6 (bottom left) from ​Ohlmann+2017. Our results are consistent with theirs as long as we use the Ledoux condition outputted from MESA directly. If we compute the compositional gradient dln mu/dln P from mu and P we get the red line instead of the blue line. I don't understand why, but anyhow, as we will see in the next graph, both lines are consistent in that they predict convective instability. The black dash-dotted curve would be the relevant curve for the fiducial RGB simulation with adiabatic EoS (gamma1=5/3). We see that for that simulation, the star is predicted to be convectively stable.
   5. ​Graph showing condition for convective instability, plotted in a way that more clearly shows whether profile is stable or unstable. We clearly see that for the Ledoux criterion (blue line), which is the relevant criterion in our case, we are convectively unstable everywhere outside the softening radius. Note that the red line should match the blue line, but doesn't for some reason (see above), but in any case the red line also predicts convective instability. So to summarize, we see that the profile is convectively stable if the adiabatic gamma1=5/3 EoS is assumed, but convectively unstable for the more realistic EoS, consistent with Ohlmann+2017.
   6. ​Graph showing convective speed from MESA model. We see that the typical convective speed is ~0.7 km/s, rising to ~1 km/s near the softening radius and steadily rising to ~4 km/s near the surface.
   7. ​Graph showing convective Mach number from MESA model. This can be compared with Figure 9a from ​Ohlmann+2017.
   8. ​Graph showing convective time scale computed from MESA. The convective time scale is computed either by taking H_P/v_conv where H_P is the pressure scale height, or H_P²/eta_MLT where eta_MLT is the diffusion coefficient from mixing length theory, outputted by MESA, or L_MLT²/eta_MLT where L_MLT is the mixing length outputted by MESA.
2. AstroBEAR results
   1. I ran the simulation up to 2 frames on bluehive. It slowed down by about a factor of 3 compared to the first frame. The reason is probably that the code was struggling to resolve the large gradients near the particles. The density reduced markedly around the RGB core particle between frames 1 and 2, which does not happen in the fiducial run. The velocity with respect to the RGB core particle points outward. Could this be caused by convection? Or is it just that more resolution is needed at the RGB core when the MESA EoS is used? (Note that inside the softening radius, the profile is the reconstructed profile, not the original MESA profile.)
      • ​Zoom-in on RGB core, MESA EoS Run 207
      • ​Zoom-in on RGB core, fiducial Run 143
   2. The gas profile around the secondary is quite different from the fiducial run.
      • ​Zoom-in on secondary, MESA EoS Run 207
      • ​Zoom-in on secondary, fiducial Run 143 Although it would be worth adding more resolution around the secondary, the difference seen is probably physical and likely has to do with the smaller gamma1 near the surface in the MESA EoS model. A smaller gamma1 means that more compression is required to achieve a given pressure, so we might expect greater density in Run 207, which is what is seen.

Summary

• Pressure is dominated by gas pressure.
• Temperature is as expected for an ideal gas.
• Adiabatic gamma = gamma1 = C_P/C_V (with C_P and C_V the specific heats at constant pressure and volume) = (dlnP/dlnrho)_S (at constant entropy), which is 5/3 for a monatomic ideal gas. Here we obtained gamma1 from MESA and find that it is close to the polytropic Gamma1 in the envelope (as expected for a convective envelope) except near the surface.
• We then used gamma1 along with the stellar profile, to compute Del-Del_Ledoux at each point, where Del=dlnT/dlnP and Del_Ledoux=(gamma1-1)/gamma1 +dln mu/dln P. We find the envelope to be convectively unstable outside the softening radius. For the adiabatic EoS with gamma1=5/3 as used in the fiducial model, we find that the envelope is stable to convection.
• Our results for convective instability for the adiabatic gamma1=5/3 model (stable) and also for the MESA tabular EoS model (unstable) agree with those presented in Ohlmann+2017 for their profile, which is very similar to ours.
• The convective time scale varies between a few days near the surface, to ~50-400 days in the bulk of the envelope, to ~10-50 days at the softening radius. Thus, we would expect the surface layers to be somewhat less stable than in the fiducial run during the dynamical plunge-in (about 12.5 days to the first periastron passage in the fiducial run). Futhermore, we would likely see more mixing inside of the particle orbit at late times in the simulation since the simulation time (40 days) is comparable to the convective time near the softening radius.
• More resolution is probably needed at the primary core.