XRB Paper Data/Progress - Initial Model
Andrew provides us with an initial model in column-depth space. Column-depth is defined by:
dy = - rho * dr
where y is the column-depth and r is the radial coordinate. In a plane-parallel approximation, using column-depth as the independent coordinate, the equation of hydrostatic equilibrium reduces to
P = |g| * y,
where |g| is the magnitude of gravitational acceleration. Below is a plot of the raw initial model in column-depth space.
The "kink" in these curves is where the composition is changing from pure Fe-56 to pure He-4. This composition jump occurs over only one zone and as such we will want to smooth this out before putting it into our Eulerian code. Furthermore, we have found that running these models as-is still requires a lot of computer time to get until anything interesting starts to happen in the burning layer. We therefore will emulate the slow burning leading up to the interesting part by perturbing the temperature at the base of the helium layer and running an adiabatic profile radially outwards. Any perturbations or smoothing should be done in column-depth space because any changes in the density profile will alter the physical extent of the atmosphere via the first equation above. In partiuclar, addition of an adiabatic region should cause the atmosphere to puff up a bit. Furthermore, we impose the restriction that any smoothing/perturbing be done at constant pressure, which is obtained from the second equation above.
Smoothing
To smooth the original model, we add a number of linearly interpolated column-depth points within the He-4 to Fe-56 transition region. We then use a hyperbolic tangent function to fill the data at these new gridpoints. For example, going from the outside in (i.e. y increasing) X(He4) decreases from ONE to ZERO over the transition region which is between y_base and y_base+1. Therefore we fit
X(He4)_i = -amplitude * tanh( (y_i - center) / phi ) + offset
where
amplitude = ( X(He4)_base + X(He4)_base+1 ) / TWO = HALF
center = ( y_base + y_base+1 ) / TWO
offset = ( X(He4)_base - X(He4)_base+1 ) / TWO = HALF
for the points we added to the model, y_base < y_i < y_base+1. So far, phi has been problem specific in terms of getting a smooth-looking transition - it is a bit of trial and error but we have found a value of 6.0e6 to be sufficient for this smoothing. The image below shows the smoothing of X(He4) in column-depth space.
Perturbing
Applying the temperature perturbation is rather straightforward. We simply multiply the temperature at the base by some multiplicative factor and use a Newton-Raphson method with our EOS to find the density that corresponds to the given temperature, composition and pressure. When everything is thermodynamically consistent, we call the EOS one last time to get the entropy at the base, entr_base. We then loop outward (i.e. decreasing y or i) through the zones and do a root find on the 2D Taylor expansion of:
f(rho, T) = pres_i - pres_want
g(rho, T) = entr_i - entr_base
using a Newton-Raphson iterative method. Here pres_i and entr_i are returned from the EOS given input rho and T, and pres_want is the HSE pressure given by the plane-parallel equation above. If in doing this process, our new temperature falls below that of the original data, we stop the adiabatic integration and take the hotter model. Furthermore, we allow for the option of smoothing out the temperature perturbation below (i.e. greater y) the adiabatic region using a tanh profile as described above. The end result looks like the image below.
Converting to Physical Space
To convert the smoothed and perturbed data to physical space, we need to numerically integrate the first equation above. Currently, to do this we are using a 2 point Newton-Cotes method (i.e. the trapezoid rule) to integrate between each pair of grid points. Note that because of the negative sign in the equation, y(r=0) = y_imax = ymax; to account for this we say rmin = 0 = r_imax and we integrate from i=imax down to i=1:
r_i = r_i+1 - HALF * (y_i - y_i+1) * (ONE/rho_i + ONE/rho_i+1)
and r_i becomes a decreasing function of i. Another method of integration would be to fit a cubic spline or Hermite cubic function to the data and integrate that - I'm going to see if this makes much of a difference. Below shows the result of putting the above smoothed and perturbed model into physical space using the above linear approximation to the function in the interval.
Once these models are put into physical space, then are then ready to
be put into HSE on a uniform grid using the init_1d routine.