MULTIDIMENSIONAL MODELING OF TYPE I X-RAY BURSTS. II. TWO-DIMENSIONAL CONVECTION IN A MIXED H/He ACCRETOR

An accreting neutron star can only build up a thin (∼10 m) surface layer of H/He before the immense gravitational acceleration compresses this fuel to the point of ignition. The ensuing thermonuclear runaway is short lived (10–100 s) but releases an enormous flux of X-rays (total energy ∼10⁴⁰ erg)—a transient event we detect and classify as a Type I X-ray burst (XRB; see Lewin & Joss 1981; Lewin et al. 1993; Bildsten 2000; Strohmayer & Bildsten 2010; in 't Zand 2011 for reviews). Once the explosion subsides, the accretion builds up a fresh layer of fuel in a matter of hours to days, and a new outburst occurs. An XRB light curve shows a sharp rise—about an order of magnitude increase during ∼1 s—in the X-ray luminosity, followed by an extended (∼10 s) decay.

Some ultra-compact systems are thought to accrete pure ⁴He (4U 1820-30, for example; Cumming 2003). The most common systems, however, likely accrete a mixture of H/He from an evolved companion star (see, for example, the compilation of bursts in Galloway et al. 2008). Depending on the local accretion rate, the ¹H accreted in these systems may either (1) burn stably to form a pure ⁴He layer, which then experiences a thin-shell instability resulting in an outburst, or (2) become unstable itself in the presence of helium resulting in a mixed outburst (Fujimoto et al. 1981; Fushiki & Lamb 1987; Cumming & Bildsten 2000; Bildsten 2000). Mixed bursts typically have longer light curves due to the waiting points in the weak nuclear reactions (see Strohmayer & Bildsten 2010 for an overview).

Mixed H/He XRBs are important sites of explosive hydrogen burning via the rp process(Wallace & Woosley 1981; Schatz et al. 2001; Parikh et al. 2014). The nuclear physics of the rp-process nuclei is a focus of the U.S. Department of Energy proposed Facility for Rare Isotope Beams. Understanding the conditions that exist in XRBs is critical to accurately modeling the nucleosynthesis, which may then alter the light curve. Furthermore, the subset of XRBs exhibiting so-called Photospheric Radius Expansion burst phenomena—whereby the burst's luminosity is large enough to lift the photosphere to larger radii (lower effective temperature) before settling back down to the neutron star surface—can yield information about neutron star masses and radii (see, for example, Bhattacharyya et al. 2010; Ozel et al. 2010; Steiner et al. 2010).

Most of our theoretical understanding on XRBs comes from one-dimensional studies with stellar evolution codes, assuming spherical symmetry. These one-dimensional calculations are able to roughly reproduce the observed energies, durations, and recurrence timescales for XRBs (Taam 1980; Taam et al. 1993, 1996; Woosley et al. 2004; Fisker et al. 2008). Due to their one-dimensional nature, these simulations can use larger reaction networks than multidimensional studies to predict the nucleosynthetic yields from advanced burning stages, such as the rp process, and explore the nucleosynthesis in detail (Schatz et al. 1999, 2001; Woosley et al. 2004; Fisker et al. 2008; José et al. 2010; Parikh et al. 2013). However, the one-dimensional nature prevents the simulations from directly modeling the convection, and simplified models such as mixing length theory (see, for example, Kippenhahn & Weigert 1992) are needed. Recent multidimensional simulations have questioned the validity of mixing length theory and emphasized the role of turbulence (Meakin & Arnett 2007; Arnett et al. 2009). If the convection is not modeled properly, then the wrong temperature/pressure history for a fluid element will be obtained, affecting the nucleosynthesis and light curve.

Simulations of the vertical structure of reacting flow on neutron stars are rare. Several models of detonations (Fryxell & Woosley 1982; Zingale et al. 2001; Simonenko et al. 2012) have been done, but these calculations sample density regimes that are not typical of an XRB. Lin et al. (2006) used a low Mach number algorithm to model pure helium bursts in two dimensions, following the rise in temperature and watching the development of convection. Cavecchi et al. (2013) used a simplified hydrodynamic model (the vertical direction was treated as hydrostatic) to model flame propagation across the neutron star on length and timescales appropriate to an XRB, but the numerical technique does not allow for detailed understanding of the dynamics at the front, including mixing and turbulence. It is also important to understand whether the convection can bring ashes up to the photosphere (in 't Zand & Weinberg 2010; Bhattacharyya et al. 2010), altering our interpretation of the radiation.

In Paper I (Malone et al. 2011), we explored the convective dynamics of a pure helium XRB using our low Mach number simulation code, Maestro. Our results differed from those of Lin et al. (2006) in that we found that a much higher resolution is needed to resolve the He burning peak and to properly capture the convective dynamics. Here, we extend the low Mach number methodology to the case of mixed H/He bursts with an extended network that captures the hydrogen burning. Our ultimate goal in this series of papers is to evolve the convective region to the point where we can see a nonlinear rise in the temperature and to assess how the convection impacts the nucleosynthesis.

We use the Maestro algorithm as described in Nonaka et al. (2010). Maestro solves the equations of low Mach number hydrodynamics applicable to hydrostatic (HSE) stellar flows—that is, fluid quantities, e.g., density ρ(x, t), are decomposed into radial background, HSE components, ρ₀(r, t), and perturbational components, ρ'(x, t) that govern the dynamics (see Nonaka et al. 2010, for details). The background pressure, p₀(r, t), governs the thermodynamics and is used in place of the total pressure, p = p₀ + p', everywhere except in the momentum equation—an approximation valid for low Mach number flows up to $\mathcal {O}(M^2)$. A key feature of this equation set is that sound waves are filtered from the system, enabling long timescale evolution of highly subsonic flows. This manifests itself through an elliptic constraint on the velocity field that captures the effects of compressibility on a fluid element due to localized heating and stratification. The overall algorithm consists of an advection step, using the piecewise parabolic method; a projection step that enforces the velocity constraint; and reactions, coupled via Strang splitting for second-order accuracy (Almgren et al. 2008). Maestro is publicly available.⁴

For all simulations in this paper, we use the formulation of Maestro that couples the enthalpy equation into the solution—that is, when needed, the temperature is derived from density (ρ), enthalpy (h), and composition (X_k). In computing the fluxes through the interfaces of the cells, we need time-centered interface values of the state variables. There are a variety of quantities we can predict to the edges to construct these interface states. Experience has shown that in the presence of steep composition gradients, the advection algorithm is more stable when we predict the full state to the interfaces rather than the perturbational fluid quantities. In the present work, we predict ρ and X_k to interfaces and compute the interface state as (ρX_k)_edge = ρ_edgeX_{k, edge}, instead of predicting ρ' and averaging ρ₀ to interfaces and constructing the edge state as $(\rho _0 + \rho ^{\prime }_\mathrm{edge}) X_{k,\mathrm{edge}}$. For the enthalpy state, we now predict T to the interfaces, as this is the most sensitive thermodynamic quantity. We then use the equation of state to compute h on the interfaces and compute the enthalpy state as (ρh)_edge = ρ_edgeh(T_edge). This differs from the previous method of predicting (ρh)' to edges, averaging (ρh)₀ to edges, and then constructing the edge state as $(\rho h)_0 + (\rho h)^{\prime }_\mathrm{edge}$. This change is more in line with the original reconstruction described in Almgren et al. (2006), and that used in Zingale et al. (2013), but continues to be a subject of algorithmic investigation.

Additionally, all the models we run here include the additional term in the momentum equation identified by Klein & Pauluis (2012) and Vasil et al. (2013) that improves the energy conservation in the low Mach number limit, as well as the treatment of gravity waves (we explore this for a variety of problems in A. S. Almgren et al. 2013, in preparation). We note, however, that this additional term has little effect in a convective layer, but we include it for completeness.

2.1. Microphysics

We use an approximate network containing 10 species, based on the description of Wallace & Woosley (1981), Appendix C. This network approximates the hot CNO burning, triple-α, plus rp-process breakout and burning up through ⁵⁶Ni. For details, see the Wallace & Woosley (1981) paper, but below we describe some of the features of the approximate network. We note that our implementation follows that description faithfully, but uses updated ReacLib (Cyburt et al. 2010) rates where they exist; Table 1 lists the particular rates used at some stages in the network, along with the version identifier from ReacLib. The stiff system of ordinary differential equations is integrated using the VODE package (Brown et al. 1989).

Table 1. ReacLib Reaction Rates Used in Approximate Network

Reaction	ReacLib Reference Label
Triple-α	fy05
12C(p, γ)13N	ls09
14N(p, γ)15O	im05
14O(, e+ν)14N	wc07
14O(α, p)17F	Ha96c
15O(, e+ν)15N	wc07
15O(α, γ)19Ne	dc11
16O(α, γ)20Ne	co10
16O(p, γ)17F	ia08
17F(, e+ν)17O	wc07
17F(γ, p)16O	ia08
17F(p, γ)18Ne	cb09
18Ne(, e+ν)18F	wc07
18Ne(α, p)21Na	mv09
19Ne(, e+ν)19F	wc07
19Ne(p, γ)20Na	il10
26Si(α, p)29P	ths8
44Ti(α, p)47V	rh10

Download table as:  ASCIITypeset image

Figure 1 shows a schematic of the rp-process network outlined in Wallace & Woosley (1981) and used in this paper. The isotopes labeled in black are those explicitly included in the network, in addition to ¹H and ⁴He; likewise, black arrows indicate reaction rates explicitly calculated, whereas gray arrows denote approximations. For example, in the HCNO cycle, the sequence of reactions ¹⁴O(β⁺) ¹⁴N(p, γ) ¹⁵O is restricted by the slowest rate in the sequence, the β⁺ decay. In the approximate network, ¹⁴O (plus a proton) is converted directly to ¹⁵O at a rate given by the β⁺ decay. The two colored circles in Figure 1 denote locations where the nuclear flow can break out of the HCNO cycle and start up the rp-process path to heavier nuclei. Each of these breakout points involve competition between a β⁺ decay of an isotope of Ne and an (α, p) or (p, γ) reaction with branching ratios given by λ₁ or λ₂, respectively. For example, the chain of reactions ¹⁵O(α, γ) ¹⁹Ne(p, γ) ²⁰Na(p, γ) ²¹Mg(β⁺) ²¹Na(p, γ) ²²Mg is approximated as turning ¹⁵O (plus an α and three protons) directly into ²²Mg at a rate governed by the α capture on ¹⁵O multiplied by the fraction (1 − λ₂) of ¹⁹Ne isotopes that proton capture before they β⁺ decay. At the high temperatures and densities experienced during an XRB, many of the β⁺-decay waiting points of the traditional rp process are bypassed via (p, γ) and/or (α, p) reactions. For the flow in the rp process up through ⁵⁶Ni, the approximate network accounts for these bypassing reactions by assuming all reactions pass through the fastest path. For example, converting ²²Mg to ³⁰S can proceed through two main channels: (1) a series of eight p-captures and several β⁺ decays, or (2) two (α, p) and two (p, γ) reactions. The faster of the two paths will be used by the network. However, each individual path is limited by its slowest rate. For the first path—the β⁺-decay path—the limiting factor is the mean lifetime of the β⁺ decays; the second path is limited by the slower of the two α-capture rates, namely ²⁶Si(α, p)²⁹P. Likewise, burning ³⁰S to ⁵⁶Ni in the approximate network is restricted by a typical β⁺-decay timescale or α capture on ⁴⁴Ti. Finally, ⁴He burning via the triple-α reaction enters the diagram from the bottom left corner.

Zoom In Zoom Out Reset image size

A general stellar equation of state with ideal ions, arbitrarily degenerate/relativistic electrons, and radiation is used (Timmes & Swesty 2000). Paper I showed that including thermal diffusion did not significantly alter the average properties of the convection or the maximum temperature reached during the simulation. Consequently, we neglect thermal diffusion in the calculations of this paper. Malone et al. (2011) also discussed a volume discrepancy factor that is designed to drive the flow back onto the constraint dictated by the equation of state. For the present simulations we run with this off (equivalent to f = 0 in Equation (13) of Malone et al. 2011).

To assess the validity of the approximate network we perform a one-zone, constant pressure, self-heating burn using the conditions similar to those at the base of the accreted layer in an XRB: T = 9.5 × 10⁸ K, ρ = 5.93 × 10⁶ g cm⁻³, X(¹H) = 0.72, X(⁴He) = 0.24, and the remaining abundance split between ¹⁴O and ¹⁵O in a ratio comparable to their respective β⁺-decay times. The solid lines in Figure 2 show the early evolution of the hydrogen/helium abundance (top), specific energy generation rate (middle), and temperature (bottom) for our approximate network. The dotted lines show a comparison of a simple calculation with the same initial conditions but using an adaptive reaction network from the Kepler code-base, which included many nuclei up to Ge. The two networks agree quite well under these conditions, but the approximate network initially tends to be ≲ 10% hotter than the larger Kepler network.

Zoom In Zoom Out Reset image size

We generate parameterized models of the accreted layer using simulations run by the Kepler stellar evolution code as a guide. Creating our own models alleviates the difficulty of reproducing the proper convectively unstable regions in the atmosphere of the one-dimensional model. In particular, mapping from a Lagrangian grid to an Eulerian grid, interpolating data to a constant mesh spacing, and slight differences in the thermodynamics between codes can all lead to differences in thermal/adiabatic gradients, and hence which region is convectively unstable. Additionally, spurious features in the one-dimensional model can cause difficulties while mapping onto an Eulerian grid and maintaining HSE, especially in regions of sharp discontinuities in the Lagrangian code.

Our simplified models and their generation are described in detail in the Appendix. Briefly, our models consist of a 1.4 M_☉, 10 km, isothermal (T = T_ns) "neutron star" substrate composed of heavy elements (X = X_ns ≡ pure ⁵⁶Ni) with a warm accreted layer of mixed fuel composition (X = X_fuel), part of which is convectively unstable with an isentropic gradient. The transition between the neutron star and the accreted layer is governed by the density (ρ_base) and temperature (T_base) at the base of the accreted layer. The composition of the accreted layer was approximated from an initial Kepler XRB model, but is essentially a slightly metal-rich gas compared to solar metallicity—to reflect prior burning—with the ¹⁵O/¹⁴O ratio set approximately by their respective lifetime against β⁺ decay. Table 2 gives the set of parameters that describe the models—refer to the Appendix for the definitions.

Our first attempt at model creation used values of ρ_base (∼7 × 10⁵ g cm⁻³) and T_base (∼7.5 × 10⁸ K) determined from the original Kepler models. When mapped into two dimensions, these models exhibited weaker burning than their one-dimensional counterparts. This leads to a slowly increasing base temperature that eventually peaked around ≲ 10⁹ K, and then remained roughly constant. As we found in Paper I with semi-analytic models, the multidimensional convection generated in these models was much more efficient at carrying heat than the mixing length theory prescription used in the one dimensional Kepler models.

To create more prolific burning, we therefore modified our initial conditions to those given in Table 2. Namely, we artificially increased the density to ρ_base = 2 × 10⁶ g cm⁻³ and the temperature to T_base = 9.5 × 10⁸ K while keeping the composition fixed. Increasing the temperature and density has two main consequences for the burning in the layer: (1) boosting the triple-α reaction (∝ρ²), which is temperature sensitive, relative to the β-limited HCNO cycle, and (2) decreasing the branching ratios λ₁ and λ₂ (see Figure 1), which allows for an increase in breakout from the HCNO cycle. Ultimately, this quicker consumption of ⁴He nuclei early on will slow the late-time evolution of the burning as (α, p) reactions might not be frequent enough to bridge β⁺-decay waiting points. However, our burning does not reach that far into heavy elements so we are not concerned with this aspect here. Figure 3 shows the initial density, temperature and bulk composition profiles for the model used in our calculations.

Zoom In Zoom Out Reset image size

As with our simulations of convection in/on white dwarfs (Zingale et al. 2013; Nonaka et al. 2012), and also in Paper I, we use a low density cutoff and a sponging technique to control the behavior of the flow at the top of the atmosphere. These are designed to suppress unwanted velocities in the region above our accreted layer. The low density cutoff acts to change the behavior of the low Mach number algorithm when the density drops below a certain value. In particular, when the density drops below ρ_anel, the divergence constraint on the velocity field is altered to resemble that of the anelastic approximation, ∇ · (ρ₀U) = 0. There is another low density cutoff, ρ_cutoff, that is the density at which we fix the ambient medium. Both ρ_anel and ρ_cutoff are runtime parameters, which we set to ρ_cutoff = ρ_anel = 10³ g cm⁻³ in the present simulations. The location where the density first drops below the cutoffs is shown as the rightmost thin vertical gray line in Figure 3.

For the sponge, we pick a multiple, f_sp, of the anelastic cutoff to define the density at which the sponge turns on, ρ_sp = f_spρ_anel. The sponge is designed to be in full force when the density drops to ρ_anel. Here, we use a simplified version of the sponging as compared to our previous work,

$$\begin{eqnarray} s = \left\lbrace \begin{array}{@{} ll} \! 1, & \mathrm{if }\ r < r(\rho _\mathrm{sp}) \\ \! \frac{1}{2} (1 - s_\mathrm{min}) \cos \left(\pi \frac{r - r(\rho _\mathrm{sp})}{r(\rho _\mathrm{anel}) - r(\rho _\mathrm{sp})} \right)\\ \quad + \, \frac{1}{2}(1 + s_\mathrm{min}), &\! \mathrm{if }\ r(\rho _\mathrm{sp}) \le r < r(\rho _\mathrm{anel}) \\ s_\mathrm{min}, &\! \mathrm{otherwise,} \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 1 }$$

and after the velocity is advanced, the sponge is applied as

$$\begin{equation} U^{n+1} \rightarrow s U^{n+1}. \end{equation} \tag{ 2 }$$

For all the simulations presented here, we choose f_sp = 25 and s_min = 0.01. The leftmost vertical gray line in Figure 3 marks the location of the sponge start, r(ρ_sp) (i.e., where ρ = f_spρ_anel).

To seed the convection, we add a number, n_vort, of small vortices to the initial conditions at a fixed height, r_vort, near the base of the convective region. The vortices are spaced equally across the domain, and the orientation—clockwise versus counterclockwise—is altered every other vortex. Each vortex is Gaussian in form with the velocity perturbations given as

$$\begin{eqnarray*} \delta u_i(x,r) &=& -(-1)^iA_\mathrm{vort}(r - r_\mathrm{vort})\exp \left(-\frac{d_i^2}{\sigma ^2} \right) \\ \delta v_i(x,r) &=& (-1)^iA_\mathrm{vort}(x - x_{\mathrm{vort},i})\exp \left(-\frac{d_i^2}{\sigma ^2}\right), \end{eqnarray*}$$

where u and v are the horizontal and vertical components of the velocity, respectively, A_vort is the amplitude of the perturbation, (x_{vort, i}, r_vort) is the center of the ith vortex, d_i = [(x − x_{vort, i})² + (r − r_vort)²]^1/2 is the distance from the center of ith vortex, and σ is the size of the vortex. The superposition of all vortices is used to determine the initial velocity field, that is,

$$\begin{equation} (u,v) =\left(\sum _{i=1}^{n_\mathrm{vort}}\delta u_i,\sum _{i=1}^{n_\mathrm{vort}}\delta v_i\right). \end{equation} \tag{ 3 }$$

For the simulations presented here, we choose r_vort = 1475 cm, A_vort = 1 km s⁻¹, σ² = 200 cm², and n_vort = 16 (each x_{vort, i} is determined from the domain width and n_vort). We note that the initial convective region for this model roughly spans the region between 1470 cm ≲ r ≲ 3350 cm.

We perform simulations at three resolutions: 12, 6, and 3 cm zone⁻¹. The 12 and 6 cm runs are done with a uniform grid. For the 3 cm run, we added a single level of refinement to the 6 cm model in a region encompassing the convective zone. We do not dynamically change the grid with time; this allows us to optimize the load balancing of the simulation. Each of these simulations used a two-dimensional grid of size 1536 cm × 4608 cm. A double-wide model at 6 cm zone⁻¹ resolution was also created and mapped onto a grid of 3072 cm × 4608 cm. This fourth model comprises the bulk of the study in this paper, as it was evolved the furthest in time.

All simulations here assumed plane-parallel geometry with a constant gravitational acceleration of g = −2.45 × 10¹⁴ cm s⁻². Periodic boundary conditions were used in the lateral, or x direction. The vertical, or r-direction, upper boundary was an outflow boundary while the lower boundary used slip-wall conditions.

4.1. General Trends and Resolution Dependence

As in our simulations of Paper I, a transient event occurs at the onset of convective overturn. This transient causes relatively large velocity spikes until the system adjusts to a more steady-state convective flow field. Based on our experience in previous studies of convection with Maestro (Nonaka et al. 2012; Zingale et al. 2013), including an initial velocity perturbation, such as that in Equation (3), helps to minimize the transient by giving the system a small "kick," allowing the energy to be more quickly dispersed via advection. This is different from the approach of Paper I where we let solver noise seed convection. After a few turnover times (≲ 10⁻³ s), however, the convective pattern does not remember how it was initiated; therefore, the results presented here are not sensitive to the choice of initial perturbations.

Figure 4 shows the peak Mach number, M, in the domain as a function of time for all four of our models. At very early time, the Mach number peaks to about 0.18 for all simulations because of the transient discussed above. By about t = 0.01 s, the system has relaxed to a steady-state convective flow. While the 12 cm zone⁻¹ run differs markedly from the others, the 3 and 6 cm zone⁻¹ runs appear to be converged. Also, note that the Mach number is generally smaller for the higher-resolution runs. We suspect that the primary source of the differences with resolution is in the reaction source term. The energy generation is strongly peaked near the base of the burning layer, and, as a result, a coarse simulation will underestimate the total energy generation (since the coarse-zone cell center is at a lower temperature than the a fine-zone cell center at the base of the layer). This was also discussed in Paper I in the context of pure triple-α burning, which required a very fine spatial resolution due to its high temperature sensitivity. Also similar to the results of Paper I, the time-averaged peak velocity remains less than 10% of the sound speed, validating the use of a low Mach number approximation method for XRBs.

Zoom In Zoom Out Reset image size

Similarly, Figure 5 shows the peak temperature in the domain as a function of time for all four runs. Again, we see—excluding the 12 cm zone⁻¹ run—the temperature evolution of the models is quite similar, agreeing to within ≲ 5%. In addition to underestimating the energy generation, the coarse model with 12 cm zone⁻¹ resolution also tends to overestimate the amount of convective overshoot. This can also be inferred from Figure 4; the overall larger velocities allow for parcels of fluid on ballistic trajectories to more readily penetrate the convective boundary, thereby increasing the heat flux away from the base of the accreted layer and increasing the rate at which it can cool.

Zoom In Zoom Out Reset image size

Of the three different resolutions we tested, the 6 cm zone⁻¹ simulations offer the best tradeoff between accuracy and computational cost. The narrow-domain 6 cm zone⁻¹ model, however, has a domain width comparable to the initial vertical extent of the convective region. This forces the convective rolls to assume an aspect ratio ≲ 1, which alters the dynamics and allows for the ≲ 30% difference in the time-averaged M between the normal 6 cm zone⁻¹ model and the wide model. The double-wide domain alleviates the aspect ratio forcing for the duration of the simulations presented here. Therefore, we focus on the 6 cm zone⁻¹ wide-domain simulation for the remainder of this paper.

4.2. Convection in the Wide Domain

4.2.1. Early Adjustment

Figure 6 shows the transient adjustment phase discussed in Section 4.1, as seen in three different variables: magnitude of vorticity (left column), Mach number (middle column), and specific energy generation rate (right column). The initial conditions are in the top row; note that the vortices of Equation (3) are just barely discernible on the vorticity color scale. After 3 × 10⁻⁴ s (middle row), the nonlinear interactions of the perturbations give rise to a handful of dominant plumes that provide the large Mach number of the transient startup. Another 3 × 10⁻⁴ s later (t = 6 × 10⁻⁴ s; bottom row) and the convective pattern has fully filled the region that was originally Schwarzschild-unstable in the initial model. We see localized vortices, and where the flow converges near the bottom of the convective layer, we also have a localized increased burning rate.

Zoom In Zoom Out Reset image size

4.2.2. Well-defined Convection Region

Figure 7 shows the same evolution as Figure 6, except at later times (5 × 10⁻³, 5 × 10⁻², 10⁻¹, and 2 × 10⁻¹ s, from top to bottom). As is well known in two-dimensional simulations, smaller vortices tend to merge into larger vortices. Some of the smaller-scale features are smoothed out and the burning becomes (laterally) very uniform within the burning region. It is also clear from Figure 7 that the convective region expands both radially inward and outward, with the extent of the former being much less than that of the of latter.

Zoom In Zoom Out Reset image size

The expansion of the convectively unstable region can also been seen by looking at the Brunt–Väisälä or buoyancy frequency, N, for a stratified atmosphere (see, for example, Chapter 6 in Kippenhahn & Weigert (1992)):

$$\begin{equation} N^2 = -g \left(\frac{d\ln \rho }{dr} - \left.\frac{d\ln \rho }{dr}\right|_{\mathrm{ad}}\right), \end{equation} \tag{ 4 }$$

where the subscript "ad" means along an adiabat, and we are assuming instability to Schwarzschild convection. Equation (4) gives the square of the frequency at which a fluid element will oscillate if displaced from its equilibrium position (∝e^iNt). When N² < 0, then the buoyancy frequency is purely imaginary, giving rise to an instability. Figure 8 shows profiles of the square of the Brunt–Väisälä frequency as calculated from the laterally averaged quantities—i.e., with the gradients in Equation (4) being calculated as

$$\begin{equation*} \frac{d\ln \rho }{dr}\longrightarrow \frac{d\ln \langle \rho \rangle }{dr} \end{equation*}$$

and

$$\begin{equation*} \left.\frac{d\ln \rho }{dr}\right|_{\mathrm{ad}} \longrightarrow \left.\frac{d\ln \rho (\langle p\rangle,\langle T\rangle,\langle X\rangle)}{dr}\right|_{\mathrm{ad}}, \end{equation*}$$

where angle brackets denote laterally averaged quantities. The different curves are, from top to bottom, profiles at t = 5 × 10⁻³, 5 × 10⁻², 10⁻¹, and 2 × 10⁻¹ s, each being offset for clarity. The regions of instability (N² < 0) for each curve have been highlighted red.

Zoom In Zoom Out Reset image size

In an average sense, the upper boundary of the convective region pushes radially outward, consistent with the color map plots of Figure 7. However, at any given instance, the tenuous material just outside the upper convective boundary gives rise to fluctuations in entropy (or N²) that make a sharp distinction of the boundary difficult. Indeed, the plots of N² at the upper boundary have much more noise than at the lower boundary, where there exists a sharp transition between neutron star and atmosphere. The lower boundary is always marked by a large (positive) N². This results in stable displacements of fluid elements that give rise to internal gravity waves, such as the wave-like features at the base of the convective region in the plots of energy generation in Figure 7.

A peculiar feature evident in Figure 8 at late time (bottom curve) is a small region of convective stability around r = 2850 cm that splits the single convective zone into two layers of convection. This can be seen in either the vorticity or Mach number plots of Figure 7 at the same time (last row), which shows the presence of two layers. Indeed, a close examination of the same plots at an even earlier time, t = 10⁻¹ s (third row of Figure 7), hints at layer formation.

4.2.3. Burning and Breakout

The location of layer splitting in Figure 7 is coincident with, but below, a local maximum in the energy generation rate. Looking back at the initial conditions for the energy generation rate in the top row of Figure 6, we see that indeed there is a secondary, weaker peak above the base of the atmosphere. This peak tends to strengthen with time, while the peak of burning at the base of the accreted layer weakens (as seen in Figures 6 and 7). Figure 9 shows the lateral average (±1σ) of the energy generation rate at the same time as the snapshots in Figure 7. By t = 10⁻¹ s, the two peaks have comparable energy output, while at later times what was initially the secondary peak slightly dominates the burning at the base of the accreted layer. The secondary peak also appears to move radially outward in the atmosphere at about the same rate as the expansion of the convective zone ∼60 m s⁻¹.

Zoom In Zoom Out Reset image size

The presence of the secondary peak is a feature strengthened by the approximations in the network used here. In general, in the accreted layer, the density and temperature are decreasing functions of radius. How, then, can a constant composition profile produce a non-monotonic energy generation rate profile? The answer lies in the branching ratio between the rp-process breakout and returning to the HCNO cycle in the approximate network. In particular, for the initial conditions, Figure 10 shows in blue the branching ratio, λ₁, which was defined in Figure 1 as the fraction of ¹⁸Ne that β⁺ decays to ¹⁸F faster than it can α capture. The β⁺-decay rate is—for lack of better physics of T-dependent weak decays—a constant, whereas the α capture has to overcome a Coulomb barrier. This causes λ₁ to increase with decreasing temperature, and we see a sharp transition in the branching ratio around r ≈ 2200 cm. Note that the branching ratio is unity in the neutron star substrate as there is no ⁴He available for capture on ¹⁸Ne. Also shown in Figure 10 in red is the normalized rate at which ¹⁷F (plus a proton) is converted directly to ¹⁵O (plus an α) in the approximate network—this rate is governed by the ¹⁷F that can proton capture to ¹⁸Ne (λ_{p, γ}(¹⁷F)) and the fraction of ¹⁸Ne that β⁺-decay (λ₁). This particular rate coincides precisely with the secondary peak in the energy generation rate, which is shown (normalized) as the dashed line. It is this approximate reaction sequence that is responsible for the secondary peak in energy generation rate seen in Figures 6, 7, and 9, and for the formation of the two-layered convective structure clearly present at late time in Figure 8. We stress that the formation of distinct layers in the convective region is likely not physical, but rather an artifact of the approximate network used in these studies.

Zoom In Zoom Out Reset image size

Figure 11 shows density, temperature, and various composition profiles as a function of radius both at the beginning (dashed lines) and end (solid lines) of the 6 cm zone⁻¹ wide simulation. Note that of the isotopes plotted, only ¹⁵O and ³⁰S had an initial abundance within the range of the plot (see Table 2). The density and temperature profiles in the top panel show just how much the atmosphere has been heated and puffed up in the short duration of the simulation. The outer layers of the atmosphere have been heated substantially from the secondary energy generation peak. Likewise, the composition profiles betray the presence of layering in the atmosphere due to the secondary burning peak. The base of the accreted layer is around r ∼ 1500 cm. The relatively large abundance of, for example, ¹⁵O, ²²Mg, and ³⁰S below this region is due to convective overshoot; similarly, although not shown, a significant amount of ⁵⁶Ni is dredged up into the atmosphere from the neutron star material via the overshoot. In the first convective layer, 1500 cm ≲ r ≲ 2800 cm, the ²²Mg produced from prior burning, is readily being converted to ³⁰S. In the second layer, 2800 cm ≲ r ≲ 4000 cm, where the branching ratio λ₁ favors ¹⁵O production from ¹⁷F, we see a relatively larger abundance of ¹⁵O. There is still a significant abundance of ²²Mg in this region, mainly due to mixing throughout the convective region before it was split into two layers, but also due to some burning of ¹⁴O (via λ₂) and ¹⁶O. The lower temperature and density in the second convective layer means ²²Mg is converted to ³⁰S at a slightly slower rate than in the in the first convective layer.

Zoom In Zoom Out Reset image size

Very little ³⁰S was burned to ⁵⁶Ni. To disentangle the amount of dredge up from the amount produced from reactions, it is instructive to look at the total mass of each species. Figure 12 shows the total mass of each isotope in our network, normalized to that isotope's initial total mass, as a function of time; the inset plot shows a zoom in on the last 10 ms of evolution. There were a few isotopes—e.g., ¹⁷F and ¹⁶O—in our initial model that were out of equilibrium because we selected only the most abundant species from the Kepler models and then set the abundance of the rest of the species to zero. These species were very quickly produced in a transient event at the beginning of the simulation. The ¹⁷F produced in this transient is burned away to ²²Mg somewhat quickly initially, but is also replenished by ¹⁴O(α, p)¹⁷F; the competition of these two reaction chains results in a slowly decreasing ¹⁷F abundance. Likewise, the ¹⁴O converted to ¹⁷F is quickly replenished by burning ¹²C. Indeed, the ¹²C(p, γ)¹³N reaction rate is the most prodigious rate during the simulation, and any ¹²C is nearly instantaneously converted to ¹⁴O. The total amount of ²²Mg peaks around t ∼ 0.09 s, even though a significant amount of ³⁰S has been produced. The total amount of ⁴He and ¹H burned by mass was ≲ 27% and ≲ 1.5%, respectively. A trace amount (<10⁻⁵ of the total mass) of material was lost off the grid due to expansion of the atmosphere.

Zoom In Zoom Out Reset image size

In absolute terms, the total increase in ⁵⁶Ni mass in the atmosphere (r ⩾ 1500 cm)—accounting for mass lost from the grid—was ∼3.99 × 10¹¹ g. The total production of ⁵⁶Ni from burning was ∼2.24 × 10¹⁰ g. The nuclear production of ⁵⁶Ni was confined to the atmosphere, and, therefore, the total amount of ⁵⁶Ni dredge up during the simulation is the difference between the above masses, or ∼3.77 × 10¹¹ g. Alternatively, in our simulations, approximately 94% of the increase of ⁵⁶Ni in the atmosphere is due to dredge up from convective overshoot. We note that the severity of convective overshoot depends on both the resolution and dimensionality of the simulation. As in Paper I, we note that increasing resolution generally tends to decrease the amount of overshoot until some convergence is reached. We suspect that a three-dimensional simulation will also show reduced overshoot compared to studies presented here; we leave the three-dimensional studies for a future paper.

We presented the first multidimensional calculations of convective flow in the context of mixed H/He XRBs using realistic hydrodynamics and an approximate reaction network for H/He burning. We demonstrated convergence of our results and explored the behavior of the nucleosynthesis and convection. The approximate network we utilized assumed that some reaction chains occurred faster than what would happen in a larger network.

In addition, the approximate network utilized here resulted in the development of two layers of convection as a result of the existence of a secondary peak in energy generation rate, which eventually dominated the total energy output. We traced the cause of this secondary peak in energy generation rate to the critical rp-process breakout branching ratio, λ₁, which measures the ratio of the β⁺-decay rate versus α-capture rate of ¹⁸Ne. With our particular initial conditions, this branching ratio sharply transitions between zero and one midway through our atmosphere.

The secondary peak in energy generation rate becomes comparable to the primary peak around t ∼ 0.1 s (see Figure 9). The split of the convective region lags behind the growth of the energy generation rate peak. Indeed, the two-layered structure is not well defined until about t ∼ 0.14 s. A small entropy jump between the two layers prevents material from crossing between each layer. This causes a quenching of the mixing of isotopes and alters the burning. Coincidentally, Figure 5 shows that the rate of increase of the peak temperature changes around the time of layer formation.

We note that we artificially boosted the temperature and density at the base of our accreted layer to give prodigious burning (see the discussion in 3). Simple Kepler calculations show that in order to reproduce a model with characteristics similar to our base density, the initial metallicity of the accreted material has to be turned down significantly, Z ∼ 10⁻⁴, making such a system perhaps rare in nature. Continuing these calculations in one dimension, the full network in Kepler showed a very slight non-monotonicity of the energy generation rate in a region containing composition profiles qualitatively similar to what we see near the development of our secondary peak in the energy generation rate. However, this minor bump in Kepler's energy generation rate was transient and certainly never dominated the burning as we see in our calculations with the approximate network. The development of layered convective regions in our simulations is likely artificial and due to the approximations used in the reaction network.

These calculations are a starting point for more realistic calculations of H/He XRBs. Future work will focus on improving the nuclear physics, moving to three dimensions, and considering larger domains. We note that very little hydrogen and helium were processed during the short duration of our simulations—and similar Kepler calculations show very little change in energy generation rate profiles on this timescale. This suggests that even further (and possibly more accurate) simplifications can be made to the network for use on the short timescales amenable to multidimensional simulations. Indeed, the original approximate network of Wallace & Woosley (1981) was designed (in part) to model the full burst cycle where copious amounts of heavy elements were produced. Improvements to the nuclear physics and approximations therein are work for a future paper.

In addition, one path that will enable larger simulations is to develop a subgrid model for the burning, and using more realistic initial models (perhaps following the methodology from Cumming & Bildsten 2000, or better mapping techniques between Kepler and Maestro models). It seems that we can capture the convective behavior on a medium-resolution grid, but the steep temperature dependence in the reactions requires fine resolution for the energy generation peaks. This was much more extreme in the case of pure He bursts (Paper I) than in the calculations here; however, a potential path forward is to use subgrid resolution for the burning and average the resulting energetics and compositions back to the hydrodynamic grid. Ultimately, we would like to model a laterally propagating burning front with realistic nuclear physics. In the context of Maestro, this will require support for lateral gradients instead of a single base state. This development will be explored in tandem with the follow-up simulations described above.

The authors thank the referee for useful comments that helped clarify some points of the paper. We especially thank Stan Woosley for providing us with data from simple Kepler burns, allowing us to compare our implementation of the approximate network to a full production network. We also thank Alex Heger for sharing some Kepler X-ray burst models, on which we based the parameterized models used here. As always, we thank Frank Timmes for making his equation of state publicly available. The work at Stony Brook was supported by DOE/Office of Nuclear Physics grants no. DE-FG02-06ER41448 and DE-FG02-87ER40317 to Stony Brook. This work was also supported, at UCSC, by the National Science Foundation (AST 0909129), the NASA Theory Program (NNX09AK36G), and especially by the DOE HEP Program through grants No. DE-FC02-06ER41438 and DE-SC00010676. The work at LBNL was supported by the Applied Mathematics Program of the DOE Office of Advance Scientific Computing Research under U.S. Department of Energy contract no. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231.

We define our parameterized models to consist of three regions: an isothermal lower region, representing the underlying neutron star; a ramp-up region to bring us up to the temperature at the base of the accreted layer; and an isentropic region at the surface, which is where the convection will take place.

We specify the temperature in the lower isothermal region, T_ns; the temperature at the base of the accreted layer, T_base; the density at the base of the accreted layer, ρ_base; the width of the transition region, δ; and the lowest temperature in the convection region, T_cutoff. Additionally, we specify the composition in the accreted layer, X_fuel, and the composition in the underlying neutron star, X_ns. For all these calculations, we hold the gravitational acceleration, g, constant. The basic procedure follows the methodology outlined in Zingale et al. (2002).

We set the start of the base of the accreted layer a distance H from the bottom of the domain. The temperature and composition are set to

$$\begin{eqnarray} X_i &=& X_\mathrm{fuel} + \frac{1}{2} (X_\mathrm{fuel} - X_\mathrm{ns}) \left[ 1 + \tanh \left(\frac{r_i - H}{\delta } \right) \right] \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} T_i &=& T_\mathrm{base} + \frac{1}{2} (T_\mathrm{base} - T_\mathrm{ns}) \left[ 1 + \tanh \left(\frac{r_i - H}{\delta } \right) \right]. \end{eqnarray} \tag{ A2 }$$

This temperature serves as an initial guess in the isentropic region and will be reset there.

We choose the base of the accreted layer as the starting point of integration. We determine the base pressure and entropy from the equation of state

$$\begin{eqnarray} &&p_\mathrm{base} = p(\rho _\mathrm{base}, T_\mathrm{base}, X_\mathrm{fuel}); \quad s_\mathrm{base} = s(\rho _\mathrm{base}, T_\mathrm{base}, X_\mathrm{fuel}). \nonumber\\ \end{eqnarray} \tag{ A3 }$$

We integrate radially outward from the base, using the condition of hydrostatic equilibrium to find the needed pressure, while constraining the thermodynamic conditions to yield constant entropy. Finding the density, ρ_i, and temperature, T_i, in zone i is an iterative procedure.

• 1.  
  Pick an initial guess for ρ_i and T_i.
• 2.  
  Compute the hydrostatic pressure

  $$\begin{equation} p_i^\mathrm{HSE} = p_{i-1} + \frac{1}{2} \Delta r (\rho _i + \rho _{i-1}) g. \end{equation} \tag{ A4 }$$
• 3.  
  Define

  $$\begin{eqnarray} F &\equiv & p_i^\mathrm{HSE} - p_i^\mathrm{EOS} \end{eqnarray} \tag{ A5 }$$

  $$\begin{eqnarray} G &\equiv & s_\mathrm{base} - s_i^\mathrm{EOS} , \end{eqnarray} \tag{ A6 }$$

  where $p_i^\mathrm{EOS} = p(\rho _i, T_i, X_\mathrm{fuel})$ and $s_i^\mathrm{EOS} = s(\rho _i, T_i, X_\mathrm{fuel})$ from the equation of state. We then use the Newton–Raphson method to correct ρ_i and T_i subject to F = G = 0.
• 4.  
  We continue to iterate on this zone until we converge. If the temperature falls below T_cutoff, then we switch to constraining only the hydrostatic pressure, keeping the temperature constant.

This defines the isentropic region above the base of the accreted layer. Beneath the base (including the transition region characterized by width δ), we integrate radially inward, constraining the pressure to that demanded by HSE,

$$\begin{equation} p_i^\mathrm{HSE} = p_{i+1} - \frac{1}{2} \Delta r (\rho _i + \rho _{i+1}) g. \end{equation} \tag{ A7 }$$

Since we are using the prescribed temperature, we only need to zero $F \equiv p_i^\mathrm{HSE} - p_i^\mathrm{EOS}$ using the Newton–Raphson technique.

The code to generate these initial models is provided in the Maestro source release in MAESTRO/Util/initial_models/toy_atm.