• Users
  1. Proposal for Jason's student?
  2. Laurence Sabin

• WT
  1. New figures added in the paper

Mach number Vs b
Wind/Grid tracer ratio PDF with tracers

2. Tracers and Gaussian 2 Fit for density PDF Redid the Gaussian 2 fit of density PDF with tracers. Tried Gaussian 2 fit with simple test data to understand the GS2 fit parameters, mainly the relations of the sigmas of GS2 and the sigmas of individual component. While

there's no obvious relations between the sigmas, the sample data & fit shows clearly two peaks which matches the individual component peaks. This can't be found in the WT data with tracers. So interpret the Gaussian 2 fit for the WT density PDF as Grid & Wind material probably won't be proper.

redid tracers figure
original figure without tracers
Test Gaussian 2 fit with simple data

• Working on stellar fallback in ​WASP-12.

• Also have some analysis of WASP wind asymptoting. Took ram and thermal pressures of wind and set equal to disk pressure at required density at various radii - see Mathematica notebook​. For hydrogen extinction, need \kappa_R \lesssim 100 for HI at T = 10^4 K, \rho \approx 10^{-14} \frac{g}{cm^3}. This temperature is right on the border of H- and ff/bf opacity, so hard to say with certainty, but given the relatively low density it should be fine. (​link)

• ​Paper is in decent shape. Still needs some editing and a sentence here or there (figures and a few references, as well), but meat is done.

One of the next steps identified in the last post was to try increasing the ambient pressure to help to prevent the near-surface layers of the star from expanding outward. There are two problems with the (modified) MESA profile near the surface r=R:
1) Due to a sudden drop in pressure as r\rightarrow R, the pressure scale height drops drastically, quickly becoming unresolvable.
2) The pressure profile is not smooth at r=R; i.e. \partial P/\partial r is discontinuous across the surface. This could potentially lead to numerical problems.

We comment on (2) at the end of this post; here we focus on (1). Problem (1) can possibly be addressed by increasing the ambient pressure. This would have the affect of moving the effective surface of the star (P=P_{\mathrm{ambient}}) inward. It is important to be aware how the minimum pressure scale height would be affected. ​Here I have plotted the pressure vs radius (top) and pressure scale height vs radius (bottom) for both modified and unmodified profiles. Three horizontal lines in the bottom plot show the chosen scale height cutoff, while the horizontal lines in the top plot show the corresponding pressure cutoffs. For example, if we want the minimum scale height to be 1R_\odot, then we should set the ambient pressure to 2.83\times10^5dyne/cm^2. This would cause the pressure to equal the ambient pressure at r=45.46R_\odot rather than at r=R=48.06R_\odot.

Given that plot, it seems reasonable to try imposing ambient pressures of 1\times10^5, 1\times10^6 and 1\times10^7 dyne/cm^2. (Recall that our resolution for 256^3 is about 0.5R_\odot.) This was implemented by simply setting the pressure equal to the maximum value of the profile and the specified ambient pressure at each radius. The movies are shown here, with initial and final frames compared below each movie. All units are CGS:
A) P_\mathrm{ambient} set implicitly by AstroBear to a low value: ​3d density ​1d density ​3d pressure ​1d pressure ​1d
​3d density start and end comparison ​1d density start and end comparison ​3d pressure start and end comparison ​1d pressure start and end comparison

B) P_\mathrm{ambient}=1\times10^5 dyne/cm^2: ​3d density ​1d density ​3d pressure ​1d pressure ​1d
​3d density start and end comparison ​1d density start and end comparison ​3d pressure start and end comparison ​1d pressure start and end comparison

C) P_\mathrm{ambient}=1\times10^6 dyne/cm^2: ​3d density ​1d density ​3d pressure ​1d pressure ​1d
​3d density start and end comparison ​1d density start and end comparison ​3d pressure start and end comparison ​1d pressure start and end comparison

D) P_\mathrm{ambient}=1\times10^7 dyne/cm^2: ​3d density ​1d density ​3d pressure ​1d pressure ​1d
​3d density start and end comparison ​1d density start and end comparison ​3d pressure start and end comparison ​1d pressure start and end comparison

Discussion
Density: Comparing run (B) with the fiducial run (A), we see that the density profile ends up deviating more from that at t=0 than in run (A). The situation seems to improve as we increase P_\mathrm{ambient} in run ©. Low amplitude density waves are seen outside the star, but the density profile remains reasonably constant. In run (D), density waves outside the star are of lower amplitude than in ©, but the final density profile has been squeezed slightly in comparison with the initial profile (the star appears slightly smaller). The 1D density comparison shows the presence of an unphysical kink in the profile. Here the density wave is propagating inward, and the profile has not yet stabilized. In fact, none of the density profiles have stabilized by the end of the simulation.

Pressure: Increasing P_\mathrm{ambient} can be seen to reduce the expansion of the star's outer layers. On the other hand, the profile in the inner part of the star gets perturbed as pressure disturbances propagate inward.

We see then that imposing an ambient pressure helps to prevent the outer layers from expanding, but also causes larger inward-propagating perturbations.

The sudden transitions in \rho and P at the stellar surface (Problem (2) above) should also be addressed. Probably some combination of increasing P_\mathrm{ambient} and smoothing the profile to satisfy mass continuity and hydrostatic equilibrium at all radii (even outside the current stellar surface) is needed. The next step is therefore to think about how to smooth the profile near the surface in a reasonable way.

Update 1
It occurred to me that another possible way to get around the problem of the scale height dropping to low values near the surface is to simply truncate the stellar profile at some value. I did this by truncating the profile where the pressure drops below 1\times10^6dyne/cm^2. This happens at about 0.92R (so the star loses 8\% of its radial extent). The movies are shown here, with initial and final frames compared below each movie. All units are CGS:
E) P_\mathrm{ambient} set implicitly by AstroBear to a low value and cut off profile at P=1\times10^6dyne/cm^2: ​3d density ​1d density ​3d pressure ​1d pressure
​3d density start and end comparison ​1d density start and end comparison ​3d pressure start and end comparison ​1d pressure start and end comparison

Outgoing pressure and density waves are visible. Note also that the ambient values of pressure and density set implicitly by AstroBear are larger than in Model (A), approximately by the factor P(r_\mathrm{trunc})/P(R). The results are not very encouraging… and we conclude that truncating the star just below the surface does not seem to offer any signficant benefit.

Update 2
After discussing with Eric, one idea that came up was to try a model like (D) but with the density also held constant where the pressure is held constant. This would avoid having a small (unresolvable density scale height). This makes the situation worse, probably because it is the pressure gradient that determines the force, so might as well keep the density profile as realistic as possible.
F) P_\mathrm{ambient} set to 1\times10^7dyne/cm^2 and \rho_\mathrm{ambient} set to a constant in the same region: ​1d

As explained in the previous blog entry, the next step was to map the modified RGB profile to the AstroBear grid. In this blog entry I report the initial results.

From the statement in Ohlmann+16c that it is necessary to resolve the softening length h by about 10 grid cells, we had estimated a minimum resolution of about 0.5R_\odot or 280^3 cells, for the case h=0.1R, with R the outer radius of the star. For the case h=0.05R, the minimum resolution would be 0.25R_\odot or 560^3 cells.

Another constraint to satisfy is to resolve the scale height H=-\frac{P}{\partial P/\partial r}. We estimated that the scale height needs to be resolved by at least few grid cells, to avoid unreasonably large fluctuations in the Mach number. Near to the surface of the star r\sim R, H becomes very very small, and, moreover, dP/dr suffers a discontinuity at r=R, where the star transitions suddenly to the ambient medium. So resolving the scale height at the surface is numerically unfeasible (we return to this issue below). The scale height H may also be too small to resolve near the transition radius r=h. A plot of the pressure scale height H vs. radius r for the RGB profile, showing both original MESA and modified profiles, is available ​here for the case h=0.1R and ​here for the case h=0.05R. Note that the unphysical behaviour of H near the transition radius r=h and noise in H for r>h are probably consequences of how the IDL routine used to differentiate the pressure profile handles the different sampling in radius for r<h (obtained using modified Lane-Emden) and r>h (MESA). This is not expected to cause a problem in AstroBear, which interpolates the inputted profile for the AstroBear grid. For h=0.1R, the minimum value near r=h is H=1.8R_\odot at r=2.4R_\odot, whereas for h=0.05R, the minimum value near r=h is H=0.78R_\odot at r=1.2R_\odot. With a resolution of 0.5R_\odot, we resolve H by at least a few cells for h=0.1R, which may be sufficient. For h=0.05R, and 0.5R_\odot resolution, H is only marginally resolved by \sim 1 cell, which is apparently insufficient (see results below).

We now present plots for the density from each of the following simulations:

A) Fiducial run: resolution 256^3 (0.56R_\odot resolution), modified RGB profile with h=0.1R, box size L=10^{13}cm \approx144R_\odot. ​3d, ​1d
B) Low res run: resolution 128^3 (1.12R_\odot resolution) but otherwise same as (A). ​3d, ​1d
C) High res run: resolution 320^3 (0.45R_\odot resolution) but otherwise same as (A) (400^3 was attempted but protections caused the simulation to terminate). ​3d, ​1d
D) No iteration over core mass run: m_\mathrm{c}=m_\mathrm{MESA}(h), so not fully self-consistent (see previous blogs), otherwise same as (A). ​3d, ​1d
E) Small softening length run: h=0.05R instead of 0.1R, otherwise same as (A). ​3d, ​1d
F) No spline run: Spline softening of potential not employed in AstroBear, otherwise same as (A). ​3d, ​1d
G) Original MESA profile but with spline potential run: No modification to profile, but otherwise same as (A). ​3d, ​1d
H) Direct MESA run: No modification to profile, nor is spline softening employed, otherwise same as (A). ​3d, ​1d

I also made side-by-side movies of the fiducial run (on the left) and an alternative run (on the right), for easy comparison:
i) (A) and (B) ​3d, ​1d
ii) (A) and © ​3d, ​1d
iii) (A) and (D) ​3d, ​1d
iv) (A) and (E) ​3d, ​1d

Discussion
From the above plots we can conclude:
I) The fiducial run (A) is quite stable at the center, as expected, but there is still a slight drop in the central density during the run. This might be alleviated by using higher resolution, but run ©, with slightly higher resolution, shows only a marginal increase in stability.
II) In all the runs, the SURFACE of the star becomes unstable on the local dynamical time, as expected.
III) A resolution of \sim 0.5R_\odot is needed for approximate stability at the center, as expected.
IV) The h=0.05R run (E) is NOT quite stable at the center when using the same resolution as for the fiducial model. This is not surprising because both h and H are smaller in this model, so not as well-resolved.
V) The iteration procedure used for the particle mass m_\textrm{c} to ensure that m(h) is equal to the MESA value makes a small but noticeable difference to the stability and final profile of the star (see the 1d plot comparison in (iii) above). If this step is not included, the density profile develops an unphysical kink at r\sim2.5\times10^{12}\textrm{cm}\sim35R_\odot.

Next steps
Some logical next steps are:
1) Figure out a viable method for reducing and smoothing out the pressure gradient at the surface. This will involve specifying an ambient value and/or an extension of the profiles for r>R. This involves coming up with some different methods, testing them, and comparing the results.
2) Once (1) is implemented, protections should no longer be a problem, so higher resolution (e.g. 512^3) could be tested.
3) It will then be useful to make other plots to study, e.g. the degree to which hydrostatic equilibrium is being satisfied, or the local Mach number of any residual motions.
4) If the star is reasonably stable after these steps, then AMR can be implemented and tested.
5) If damping of residual motions is still needed, it can first be implemented as in Ohlmann+16c.