Dynamics of Laterally Propagating Flames in X-Ray Bursts. I. Burning Front Structure

X-ray bursts (XRBs) are thermonuclear explosions in an accreted H or He layer on the surface of a neutron star (see Galloway & Keek 2017 for a review). Observations of bursts can assist in constraining the properties of the underlying neutron star, helping to illuminate the nuclear equation of state (Steiner et al. 2010; Özel & Freire 2016). Extensive observations of brightness oscillations during the rise (the initial phase of the burst where the observed flux rapidly increases) have provided evidence that the burning begins in a localized region and spreads over the surface of the neutron star (Bhattacharyya & Strohmayer 2006, 2007; Chakraborty & Bhattacharyya 2014).

One-dimensional studies of XRBs have been very successful in predicting the light curves and recurrence times (see, e.g., Woosley et al. 2004). These assume spherical symmetry and thus cannot capture the effects of localized burning spreading across the neutron star. These studies have also been used to explore the sensitivity of the burst observables to accretion and reaction rates (Cyburt et al. 2010; José et al. 2010; Lampe et al. 2016), and to model individual bursts (Johnston et al. 2019).

Multidimensional simulations of burning on a neutron star are more difficult, with both the temporal and spatial scales presenting challenges (see Zingale et al. 2018 for an overview). For the spatial scales, we need to resolve the reaction zone, ${ \mathcal O }(10\,\mathrm{cm})$ or smaller, the scale height of the atmosphere, ${ \mathcal O }(500\,\mathrm{cm})$, and the Rossby scale where the Coriolis force balances the lateral pressure gradient, ${ \mathcal O }({10}^{5}\,\mathrm{cm})$ (Spitkovsky et al. 2002). For the temporal scales, capturing the rise, ${ \mathcal O }(1\,{\rm{s}})$, and the decay of the burst, ${ \mathcal O }(10\,{\rm{s}})$, as well as the accretion period between bursts, ${ \mathcal O }({10}^{4}\,{\rm{s}})$, is currently beyond the ability of multidimensional hydro codes. Nevertheless, significant progress has been made in understanding the multidimensional nature of XRBs, through various approximations.

Laterally propagating detonations were modeled by Fryxell & Woosley (1982) and Zingale et al. (2001). However, since it is difficult to detonate helium at the densities found in normal XRBs, and even harder to detonate hydrogen because of the waiting times for weak reactions, these would only occur at very high densities. This means that detonations may only really be applicable to superbursts (where carbon is the reactant) (Weinberg et al. 2006; Weinberg & Bildsten 2007).

Global multidimensional studies were performed by Spitkovsky et al. (2002), where it was demonstrated that the Coriolis force plays an important role in confining the burning as it spreads across the neutron star surface. These calculations used the shallow-water approximation, so the vertical details of the atmosphere's structure were not captured. Their model showed that the horizontal pressure gradient between the ash and fuel can be important in accelerating the burning front.

Small-domain studies of convective burning in XRBs preceding flame development have been done in two dimensions (Lin et al. 2006; Malone et al. 2011, 2014) and three dimensions (Zingale et al. 2015). These calculations used low Mach number methods, which approximate the hydrodynamics equations to filter sound waves, enabling large timesteps and efficient modeling of subsonic convection. While these calculations could not support the lateral differences needed for flame spreading, they can help us to understand the role that convection plays in distributing the initial burning products vertically throughout the neutron star atmosphere as well as the nature of any turbulence the burning front might encounter as it propagates through the atmosphere.

The first vertically resolved simulations of lateral deflagrations were obtained by Cavecchi et al. (2013), who showed how the Coriolis force creates a geometrical configuration that increases the flame speed set by conduction by a factor of $\sim {L}_{R}/H$, where L_R is the Rossby radius and H is the scale height of the burning layer. The effect of changing Coriolis confinement across the surface on the flame propagation was explored in Cavecchi et al. (2015), while Cavecchi et al. (2016) showed how magnetic field tension, opposing the Coriolis force, can either speed up or slow down the flame by changing the horizontal extent of the flame front. Finally, Cavecchi & Spitkovsky (2019) studied the effects in 3D of the baroclinic instability at the flame front, measuring flames up to 10 times faster than in the 2D case.

The computational expense of the problem necessitated several approximations to make the simulations of Cavecchi et al. (2013, 2015, 2016) and Cavecchi & Spitkovsky (2019) tractable. The first of these approximations was the assumption of vertical hydrostatic equilibrium, which allowed them to use a less refined grid and filter out vertically propagating sound waves that would otherwise have limited the timestep. This enabled them to reach the simulation times (on the order of seconds) needed to capture the rise of the burst, but they could not fully resolve vortical fluid motions driven by effects such as turbulence and convection. These can facilitate mixing and help to accelerate the propagation of the flame front. They also included only the triple-α process in their implementation of nuclear burning, and used a fixed, constant opacity across the entire domain for each simulation run (although they explored the dependence of the flame propagation speed on this opacity value). Our simulations relax the assumption of vertical hydrostatic equilibrium, and implement a more complete set of nuclear reactions and opacities. However, as discussed later, this comes at the expense of a different set of approximations.

For deflagrations, we either need to resolve the structure of the reaction zone or use a flame model. Flame models usually assume that the flame structure is thin compared to the size of the system (see, e.g., Röpke et al. 2007 for applications to SNe Ia). For XRBs, however, the flame thickness is comparable to the scale height of the atmosphere, so we cannot use these approximations. Accurate models of flames in XRBs therefore require that we resolve the thermal width, which is ${ \mathcal O }(10\,\mathrm{cm})$ for helium flames (Timmes 2000).

The goal of this study is to understand what numerical and physical approximations are required to perform a full hydrodynamical, multidimensional simulation of flame propagation through the atmosphere of a neutron star. For simplicity in this first set of calculations, we will use a pure helium composition. These studies complement the prior multidimensional studies described above in helping us to build a picture of the dynamics of X-ray bursts.

All simulations are performed with the Castro hydrodynamics code (Almgren et al. 2010; see also Zingale et al. 2018 for a recent description). We evolve the system of fully compressible Euler equations for reacting flow:

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {X}_{k})}{\partial t}=-{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}}{X}_{k})+\rho {\dot{\omega }}_{k}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial (\rho {\boldsymbol{U}})}{\partial t} & = & -{\rm{\nabla }}\cdot (\rho {\boldsymbol{UU}})-{\rm{\nabla }}p+\rho {\boldsymbol{g}}\\ & & -2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{U}}-\rho {\boldsymbol{\Omega }}\times ({\boldsymbol{\Omega }}\times {\boldsymbol{r}})\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial (\rho E)}{\partial t}=-{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}}E+p{\boldsymbol{U}})+{\rm{\nabla }}\cdot {k}_{\mathrm{th}}{\rm{\nabla }}T+\rho \dot{\epsilon }\\ +\,\rho {\boldsymbol{U}}\cdot {\boldsymbol{g}}-\rho ({\boldsymbol{\Omega }}\cdot {\boldsymbol{r}})({\boldsymbol{\Omega }}\cdot {\boldsymbol{U}})+\rho | {\boldsymbol{\Omega }}{| }^{2}({\boldsymbol{U}}\cdot {\boldsymbol{r}}).\end{array}\end{eqnarray} \tag{ 3 }$$

Here, ρ is the mass density, ${\boldsymbol{U}}$ is the velocity, p is the pressure, and E is the specific total energy, which is related to the specific internal energy as $e=E-| {\boldsymbol{U}}{| }^{2}/2$. The forcing in the momentum equation includes gravity, described by gravitational acceleration ${\boldsymbol{g}}$, and rotational forces, described by angular velocity ${\boldsymbol{\Omega }}$, with ${\boldsymbol{r}}$ being the position vector from the origin. Species are described by mass fractions, X_k (such that ${\sum }_{k}{X}_{k}=1$), and creation rates, $\dot{\omega }$, and are related to the total specific energy generation rate, $\dot{\epsilon }$. The total mass conservation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}})\end{eqnarray} \tag{ 4 }$$

implies ${\sum }_{k}{\dot{\omega }}_{k}=0$. An equation of state of the form $p=p(\rho ,e,{X}_{k})$ completes the thermodynamic description of the system. Thermal diffusion is described by a thermal conductivity ${k}_{\mathrm{th}}$ and temperature T.

Castro uses an unsplit piecewise parabolic method (PPM) with characteristic tracing for solving the hydrodynamics (Colella & Woodward 1984; Miller & Colella 2002), generalized to an arbitrary equation of state (Zingale & Katz 2015). Reactions are incorporated via Strang splitting (Strang 1968), giving a method that is overall second-order accurate in space and time. Castro uses the AMReX adaptive mesh refinement library (Zhang et al. 2019) to manage a hierarchy of grids at different resolutions.

Since the neutron star rotates, we work in a corotating frame, taking the angular velocity ${\boldsymbol{\Omega }}$ to be constant. Furthermore, for the 2D simulations presented here, we work in axisymmetric coordinates, but we advect a third component of velocity, coming out of the simulation plane, that participates in the Coriolis force (sometimes described as a 2.5D simulation). We will take r to be the cylindrical radial coordinate with corresponding velocity u, z to be the cylindrical vertical coordinate with corresponding velocity v, and the θ to be the cylindrical azimuthal coordinate with corresponding velocity w. As r and z define the simulation plane, we adopt the nonstandard ordering (r, z, θ), for which a right-handed coordinate system has positive w pointing out of the page. We will take ${\boldsymbol{\Omega }}={{\rm{\Omega }}}_{0}\hat{{\boldsymbol{z}}}$ for the angular rotation rate, and ${\boldsymbol{g}}=-g\hat{{\boldsymbol{z}}}$ for the gravitational acceleration, with g constant. With these choices, the Coriolis force is:

$$\begin{eqnarray}&&-2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{U}}=-2\rho \left({{\rm{\Omega }}}_{0}w\hat{{\boldsymbol{r}}}-{{\rm{\Omega }}}_{0}u\hat{{\boldsymbol{\theta }}}\right).\end{eqnarray} \tag{ 5 }$$

We will neglect the centrifugal force—with our plane-parallel geometry, this will act only in the lateral direction, and is not expected to greatly affect the dynamics. Carrying the Coriolis force allows us to capture the geostrophic balance that sets up via lateral hydrostatic equilibrium (Spitkovsky et al. 2002).

Writing the momentum equation in terms of the u, v, and w components, and neglecting the centrifugal force, we have:

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho u)}{\partial t}+{\rm{\nabla }}\cdot (\rho u{\boldsymbol{U}})+\displaystyle \frac{\partial p}{\partial r}=-2\rho {{\rm{\Omega }}}_{0}w\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho v)}{\partial t}+{\rm{\nabla }}\cdot (\rho v{\boldsymbol{U}})+\displaystyle \frac{\partial p}{\partial z}=-\rho g\end{eqnarray} \tag{ 7 }$$

where we cancel $\partial p/\partial \theta$ because there are no variations in the azimuthal direction. This allows us to recast the w-velocity equation as a simple advection equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial w}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}w=2{{\rm{\Omega }}}_{0}u.\end{eqnarray} \tag{ 9 }$$

In our geometry, the flame will propagate from left to right, so u will be positive and the Coriolis force results in $w\gt 0$ (out of the simulation plane). The algorithmic implementation of rotation in Castro is described in Katz et al. (2016).

We use a general stellar equation of state with nuclei (treated as an ideal gas), photons, and degenerate/relativistic electrons, as described in Timmes & Swesty (2000). To model reactions, we use a 13-isotope alpha chain network derived from the aprox13 network (Timmes 2019). For one of our runs, we use the smaller 7-isotope network described in Timmes et al. (2000). We integrate the network using the VODE integration package (Brown et al. 1989), and our implementation is provided in the StarKiller Microphysics source (The StarKiller Microphysics Development Team et al. 2019). We note that we do not explicitly model viscosity. The reactions will provide the small-scale cutoff to the instabilities and turbulence at the flame front. We also do not include species diffusion—astrophysical flames tend to have large Lewis numbers, so this is not expected to be important (Timmes & Woosley 1992). Finally, we use the thermal conductivities described in Timmes (2000).

All simulations use adaptive mesh refinement to refine on the atmosphere (leaving the space between the top of the atmosphere and upper boundary at low resolution). As seen in Figure 1, we use up to three refinement levels in addition to the base grid, the first one is a factor of 4 finer than the previous and the remaining are each a factor of 2 finer than the previous. Castro subcycles in time, so the finer grids are evolved with a finer timestep than the coarse grids. Occasionally, the timestep chosen at the start of a cycle will violate the CFL condition during the advancement of the finer grids. In this case, we restart the finer grid evolution with a smaller timestep, subcycling within the larger timestep hierarchy. We use a CFL number of 0.8 for our simulations. The base grid for our standard simulations is 768 × 192 zones, and the equivalent finest grid when three refinement levels are added would be 12288 × 3072 zones. Our standard domain is $1.2288\,\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$, corresponding to 10 cm resolution on the finest grid. We only refine the fuel layer in the left half of the domain at the highest resolution (and only down to densities of $2.5\times {10}^{4}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$), since this is where we expect the flame to propagate. At the start of the simulation, 3.4% of the domain is at the finest resolution. This increases to 7.4% by the end of the simulation, because of the increase in the scale height of the atmosphere behind the flame.

Zoom In Zoom Out Reset image size

Thermal diffusion is modeled explicitly, using a predictor-corrector scheme to achieve second-order accuracy. A verification test of the diffusion scheme is shown in the Appendix. The explicit thermal diffusion requires a timestep limiter of the form:

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{\mathrm{diff}}\leqslant \displaystyle \frac{1}{2}\displaystyle \frac{{h}^{2}}{{ \mathcal D }},\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal D }={k}_{\mathrm{th}}/(\rho {c}_{v})$ is the thermal diffusivity and h is the width of a cell (either ${\rm{\Delta }}r$ or ${\rm{\Delta }}z$). The diffusivity increases rapidly at the top of the atmosphere, causing these low density regions to determine the overall timestep for the simulations. Therefore, we disable thermal conduction at low densities where it is not expected to be important.

We use hydrostatic boundary conditions on the lower boundary, using a discretized hydrostatic equilibrium equation of the form:

$$\begin{eqnarray}&&{p}_{i}={p}_{i-1}+\displaystyle \frac{1}{2}{\rm{\Delta }}z({\rho }_{i}+{\rho }_{i-1}){\boldsymbol{g}}\cdot \hat{{\boldsymbol{z}}}\end{eqnarray} \tag{ 11 }$$

and holding the temperature constant in the ghost cells. This is solved together with the equation of state. The velocity is reflected at this boundary. This procedure follows the form described in Zingale et al. (2002). The left boundary is reflecting and the right boundary is a zero-gradient outflow. The top boundary sets the state to simply the conditions in our outer buffer region of the initial model (see below), with the normal velocity set to the larger of zero or the velocity at the top of the domain (this prevents incoming velocities at the top) and the transverse velocities set to zero.

When we begin the simulation, there is a transient phase as the flame gets established. Material that is forced upward will encounter the steep density gradient at the top of the atmosphere and accelerate as it is blown out of the atmosphere. Eventually this material will fall back to the top of the atmosphere. In this paper, we are mostly concerned with the behavior of the flame and not any material that is violently blown out of the top of the atmosphere, so we apply a sponge to this region. This is similar to the method we previously used in Malone et al. (2014), and takes the form of a source term to the momentum and energy equations of the form:

$$\begin{eqnarray}&&{{\boldsymbol{S}}}_{\rho {\boldsymbol{U}}}=\rho {\boldsymbol{U}}\displaystyle \frac{f}{{\rm{\Delta }}t}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{S}_{\rho E}=\rho {\boldsymbol{U}}\cdot {{\boldsymbol{S}}}_{\rho {\boldsymbol{U}}}\end{eqnarray} \tag{ 13 }$$

with the sponge forcing f dependent on the density. We define the sponge shape as:

$$\begin{eqnarray}s=\left\{\begin{array}{cc}0 & \rho \gt {\rho }_{\mathrm{upper}}\\ \displaystyle \frac{1}{2}\left[1-\cos \left(\displaystyle \frac{\pi (\rho -{\rho }_{\mathrm{upper}})}{{\rm{\Delta }}\rho }\right)\right] & {\rho }_{\mathrm{upper}}\geqslant \rho \gt {\rho }_{\mathrm{lower}}\\ 1 & \rho \lt {\rho }_{\mathrm{lower}}\end{array}\right..\end{eqnarray} \tag{ 14 }$$

Here ${\rho }_{\mathrm{upper}}$ and ${\rho }_{\mathrm{lower}}$ are the densities where the sponge transitions to being fully applied. We take ${\rho }_{\mathrm{upper}}={10}^{2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and ${\rho }_{\mathrm{lower}}=1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, with ${\rm{\Delta }}\rho ={\rho }_{\mathrm{upper}}-{\rho }_{\mathrm{lower}}$. The sponge update is done implicitly to get the effective forcing, f:

$$\begin{eqnarray}&&f=-\left[1-\displaystyle \frac{1}{1+\alpha s}\right]\end{eqnarray} \tag{ 15 }$$

with $\alpha ={\rm{\Delta }}t/{\tau }_{\mathrm{sponge}}$. Here ${\tau }_{\mathrm{sponge}}$ is the timescale over which the sponge acts. We take ${\tau }_{\mathrm{sponge}}={10}^{-7}\,{\rm{s}}$. The sponge drives the velocity of the material in the low density regions at the top of the atmosphere to zero. This sponging helps increase our timestep as well.

The speed and thickness of a laminar helium flame are determined by the energy generation rate and conductivity, and scale roughly as

$$\begin{eqnarray}&&{s}_{L}\approx \sqrt{{k}_{\mathrm{th}}\dot{\epsilon }}\qquad {\lambda }_{L}\approx \sqrt{\displaystyle \frac{{k}_{\mathrm{th}}}{\dot{\epsilon }}}\end{eqnarray} \tag{ 16 }$$

(O'Rourke & Bracco 1979; Khokhlov 1993). At the densities we consider in this simulation, the pure He flame speed is quite slow and would require long integration times to see significant evolution of the burning. We therefore consider boosted flames in this first paper, to accelerate the burning and allow us to understand the qualitative effects of laterally propagating flames. To boost the flame while keeping the thickness the same, we can multiply both the burning rate and the conductivity by the same factor. For our standard calculations, we choose 10 for each, to give a 10× faster flame speed. We call this the "10/10" flame. We will also do a simulation with the reactions and conductivity both boosted by 5, the "5/5" flame.

To understand the time and length scales involved in flame propagation, we do a 1D simulation of a laminar flame using our microphysics. Figure 2 shows the flame thermodynamic profile and properties for the 10/10 flame using our conductivities and the aprox13 reaction network. This flame had a density of $2\times {10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and temperature of $5\times {10}^{7}$ K. We observe that this flame speed is about ${10}^{5}\,\mathrm{cm}\,{{\rm{s}}}^{-1}$, the flame width is about 40 cm, and it takes about 3 ms to settle into a sustained flame. Note that this speed is quite small compared to the speeds of $\sim {10}^{6}\,\mathrm{cm}\,{{\rm{s}}}^{-1}$ estimated in Spitkovsky et al. (2002). Table 1 gives the properties of the 10/10 and 5/5 laminar flames. We expect a multidimensional flame to accelerate due to hydrodynamics interactions (wrinkling, turbulence interactions, directed flows feeding fuel into the flame, etc.).

Zoom In Zoom Out Reset image size

We measure the laminar flame width as:

$$\begin{eqnarray}&&{\lambda }_{L}\equiv \displaystyle \frac{{\rm{\Delta }}T}{\max \{| {\rm{\nabla }}T| \}}.\end{eqnarray} \tag{ 17 }$$

Experience with modeling resolved flames suggests that we need a spatial resolution, ${\rm{\Delta }}x$ of ${\lambda }_{L}/{\rm{\Delta }}x\sim 5$ (Bell et al. 2004). These conditions represent the bottom of the He layer. As the density decreases with altitude, the flame thickness increases and the flame speed decreases, so we will easily resolve the flame structure throughout the rest of the atmosphere.

We wish to create an initial atmosphere consisting of a hot "post-flame" region and a cooler atmosphere that the flame will laterally propagate into. We put the hot region at the very left of the domain (the origin of the axisymmetric coordinates). To create these initial conditions, we produce two different hydrostatic models, a "hot" model that will represent the perturbation that drives the flame and a "cool" model that will represent the state ahead of the flame. These will have different scale heights. To create these models, we break the vertical structure of the atmosphere into four layers: (1) the underlying neutron star, (2) a ramp-up to the base of the accreted atmosphere, (3) a fuel layer representing the bulk of the atmosphere where the flame will propagate, and (4) an outer, low density, isothermal buffer above the atmosphere that allows us to capture expansion and explosive dynamics.

The temperature profile in the star and ramp region is given as:

$$\begin{eqnarray}&&T(z)={T}_{\star }+\displaystyle \frac{1}{2}({T}_{\mathrm{hi}}-{T}_{\star })\left[1+\tanh \left(\displaystyle \frac{\tilde{z}}{2{\delta }_{\mathrm{atm}}}\right)\right]\end{eqnarray} \tag{ 18 }$$

with

$$\begin{eqnarray}&&\tilde{z}=z-{H}_{\star }-\displaystyle \frac{3}{2}{\delta }_{\mathrm{atm}}.\end{eqnarray} \tag{ 19 }$$

Here, ${\delta }_{\mathrm{atm}}$ is a characteristic width of the transition ramp, ${T}_{\star }$ is the temperature of the neutron star, and T_hi is the highest temperature in the HSE model—it will represent the base of the fuel layer. We also define a floor temperature, T_lo, which the atmosphere reaches and is held constant at in the outer, isothermal layer.

The species mass fractions use this same profile, switching from a set describing the underlying star, ${{X}_{k}}_{\star }$, and the set for the accreted material, ${{X}_{k}}_{\mathrm{atm}}$, which is used in the isentropic and outer regions. Note that since the profile above is linear in X, if the initial mass fractions sum to one, then the blended mass fractions in the ramp region also sum to one.

We specify the density, ${\rho }_{\mathrm{int}}=\rho (z={H}_{\star })$ as the starting point for the integration of hydrostatic equilibrium. This is just below the ramp-up region—this ensures that regardless of what the peak temperature (T_hi) is, the state beneath the ramp-up region remains unchanged. Therefore, we will still be in lateral equilibrium in the star region. We will denote the density where $T={T}_{\mathrm{hi}}$ as ${\rho }_{\mathrm{fuel}}$.

Creating the model involves specifying ${T}_{\star }$, T_hi, T_lo, ${\rho }_{\mathrm{int}}$, ${H}_{\star }$, ${\delta }_{\mathrm{atm}}$, ${{X}_{k}}_{\star }$, ${{X}_{k}}_{\mathrm{atm}}$, and g. We then integrate outwards from the base of the ramp region ($z={H}_{\star }$), enforcing the discrete form of hydrostatic equilibrium, Equation (11). Integrating upwards, we find p_i and ${\rho }_{i}$ using a Newton–Raphson solver together with the equation of state, with either the temperature specified, ${T}_{i}=T({p}_{i},{\rho }_{i},{\{{X}_{k}\}}_{i})$ (in the neutron star, ramp, and isothermal buffer layers) or constant entropy, ${s}_{i}=s({p}_{i},{\rho }_{i},{\{{X}_{k}\}}_{i})$, in the fuel layer. This follows the procedures described in Zingale et al. (2002). We use a constant temperature for all $z\lt {H}_{\star }+3{\delta }_{\mathrm{atm}}$. Above this, we switch to an isentropic atmosphere until the temperature drops to T_lo, at which point we again keep the temperature constant. The integration of the atmosphere continues until the density falls to a low density cutoff, ${\rho }_{\mathrm{cutoff}}$. The material above this height is taken to have constant density and temperature. The choice of factors in front of ${\delta }_{\mathrm{atm}}$ were designed to make sure the peak T is attained at the desired density of the burning layer. The parameters we use for the model generation are listed in Table 2 and the initial model profiles are shown in Figure 3.

Zoom In Zoom Out Reset image size

Table 2.  Initial Model Parameters

Parameter	Cool	Hot
${T}_{\star }$	108 K
Thi	$2\times {10}^{8}$ K	$1.4\times {10}^{9}$ K
Tlo	$8\times {10}^{6}$ K
${\rho }_{\mathrm{int}}$	$3.43\times {10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$
${\rho }_{\mathrm{fuel}}$ a	$2.36\times {10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$	$1.20\times {10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$
${\rho }_{\mathrm{cutoff}}$	${10}^{-4}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$
${H}_{\star }$	2000 cm
${\delta }_{\mathrm{atm}}$	50 cm
${X}_{\star }({}^{56}\mathrm{Ni})$ b	1.0
${X}_{\mathrm{atm}}({}^{4}\mathrm{He})$ b	1.0
${\boldsymbol{g}}$	$-1.5\times {10}^{14}\,\mathrm{cm}\,{{\rm{s}}}^{-2}\hat{{\boldsymbol{z}}}$

Notes.

^aThis is not an input parameter, but instead is computed during integration. We list it here for reference. ^bAll other species are taken as 0.

Download table as:  ASCIITypeset image

We blend the hot and cold models laterally to produce the perturbation needed to initiate a localized flame, with the hot model at the origin of the axisymmetric geometry. The blending is done as:

$$\begin{eqnarray}&&p(r,z)=f(r){p}_{\mathrm{hot}}(z)+[1-f(r)]{p}_{\mathrm{cool}}(z)\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\rho (r,z)=f(r){\rho }_{\mathrm{hot}}(z)+[1-f(r)]{\rho }_{\mathrm{cool}}(z)\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{X}_{k}(r,z)=f(r){{X}_{k}}_{\mathrm{hot}}(z)+[1-f(r)]{{X}_{k}}_{\mathrm{cool}}(z)\end{eqnarray} \tag{ 22 }$$

with

$$\begin{eqnarray}f(r)=\left\{\begin{array}{ll}1 & r\lt {r}_{\mathrm{pert}}\\ 1-\displaystyle \frac{r-{r}_{\mathrm{pert}}}{{\delta }_{\mathrm{blend}}} & {r}_{\mathrm{pert}}\leqslant r\leqslant {r}_{\mathrm{pert}}+{\delta }_{\mathrm{blend}}\\ 0 & r\gt {r}_{\mathrm{pert}}+{\delta }_{\mathrm{blend}}.\end{array}\right.\end{eqnarray} \tag{ 23 }$$

Since the equation of hydrostatic equilibrium is linear and our blending is a linear combination of two models in hydrostatic equilibrum, the blended model is also in vertical equilibrium initially. We choose ${r}_{\mathrm{pert}}=1.024\times {10}^{4}\,\mathrm{cm}$ and ${\delta }_{\mathrm{blend}}\,=2048\,\mathrm{cm}$. Once the blended model is constructed, we compute $T(r,z)$ and $(\rho e)(r,z)$ from the equation of state,

$$\begin{eqnarray}&&T(r,z)=T(\rho (r,z),p(r,z),{X}_{k}(r,z))\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&(\rho e)(r,z)=\rho (r,z)\cdot e(\rho (r,z),p(r,z),{X}_{k}(r,z)).\end{eqnarray} \tag{ 25 }$$

The initial model is created on a uniform grid at the resolution corresponding to the finest level of refinement. At regions of the atmosphere that are not refined to the finest level, we interpolate density, pressure, and composition on the grid and then obtain the temperature from the EOS.

To assess the sensitivity of the flame propagation to the various approximations we made, we ran a suite of simulations. Table 3 summarizes these simulations. The majority of them used a reaction rate boosting of 10 and a conductivity boosting of 10, which should increase the flame speed by a factor of 10. Most simulations used a resolution of 10 cm and a domain width of a little more than 1 km. We use an artifically high rotation rate of $2000\,\mathrm{Hz}$, which gives a Rossby length of

$$\begin{eqnarray}&&{L}_{R}=\displaystyle \frac{\sqrt{{{gH}}_{0}}}{{\rm{\Omega }}}\sim 3\times {10}^{4}\,\mathrm{cm}\end{eqnarray} \tag{ 26 }$$

using a scale height ${H}_{0}={10}^{3}\,\mathrm{cm}$. This is about one-quarter of the domain width. The run at $20\,\mathrm{cm}$ resolution used one fewer level of refinement. The slower rotating case (1000 Hz) uses a slightly wider domain to accommodate the expected larger Rossby length. We note that the entire simulation framework for these calculations is freely available in the Castro github repository.⁶ In the discussions below, we will use the simulation names defined in the table to refer to specific runs and we will use the 10/10 run as the reference calculation.

Table 3.  Simulation Parameters

Name	Reaction	Conductivity	Fine Grid	Rotation	Domain Size	Network
	Boost	Boost	Resolution	Rate
10/10	10	10	10 cm	2000 Hz	$1.2288\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$	aprox13
5/5	5	5	10 cm	2000 Hz	$1.2288\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$	aprox13
10/10-iso7	10	10	10 cm	2000 Hz	$1.2288\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$	iso7
10/10-lowres	10	10	20 cm	2000 Hz	$1.2288\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$	aprox13
10/10–1000 Hz	10	10	10 cm	1000 Hz	$1.8432\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$	aprox13

Download table as:  ASCIITypeset image

5.1. General Features

Figure 4 shows the time evolution of the 10/10 simulation, focusing on the mean molecular weight,

$$\begin{eqnarray}&&\bar{A}={\left(\sum _{k}\displaystyle \frac{{X}_{k}}{{A}_{k}}\right)}^{-1}.\end{eqnarray} \tag{ 27 }$$

In each frame the buffer of ${}^{56}\mathrm{Ni}$ that serves as the underlying neutron star is seen spanning the bottom of the domain. Above that, the composition begins as pure ${}^{4}\mathrm{He}$, but as the simulation progresses, the flame processes this to heavier nuclei, increasing $\bar{A}$. By about 10 ms, we see the flame front is reasonably well-defined. We see that the bottom of the burning front is lifted off of the base of the atmosphere, greatly increasing the surface area of the burning compared with a perfectly vertical flame front. By $20\,\mathrm{ms}$ the flame has moved out substantially and we are beginning to see a gradient in the composition of the ash, with the heavier nuclei furthest behind the flame. The boosting of the burning likely artificially increases this effect, as we will see in the 5/5 case below.

Zoom In Zoom Out Reset image size

Figure 5 shows the temperature, energy generation rate, and w component of the velocity (the out-of-plane velocity induced by the Coriolis force). This latter field illustrates the hurricane effect set up by the laterally spreading burning front. In the energy generation rate plot, we see that the burning is mostly concentrated toward the bottom of the layer, as expected since the density is greatest there. We see that the peak of the burning has moved off of the symmetry axis, demonstrating that the burning front is propagating to the right. In the temperature plot, we see the effect of our refinement criteria focusing only on the part of the atmosphere where we are most dense, with an artificial change in the temperature at the refinement boundary due in part to the construction of the initial model using the fine grid resolution for HSE. As we will see in the lower resolution case, this does not affect the results.

Zoom In Zoom Out Reset image size

A final feature worth noting is the ash that seems to move along the surface at a higher velocity than the flame, via surface gravity waves. The sponging that we perform is likely damping this to some extent, and the method by which we initialize the flame may induce a larger transient than in nature. Nevertheless, this surface ash is intriguing because it affects the composition of the photosphere ahead of the burning front, potentially changing our interpretation of observations. This is something that will be explored more fully in the future.

Our default resolution puts ∼80 zones vertically in the "cool" model atmosphere height. To understand the effects of resolution, we also performed a run with one fewer refinement level, giving a $20\,\mathrm{cm}$ resolution overall (this is our 10/10-lowres run). Figure 6 shows the fields for this run at 20 ms. The structure is largely the same as the 10/10 run, with largely the same flame shape and position, and the same structure in the energy generation rate. Since we do not have a jump in refinement right below the atmosphere, we do not see the temperature feature from the initial model mapping there, but we do see some cooling in the ${}^{56}\mathrm{Ni}$ region. This is likely a resolution effect, again due to the strong gradient in temperature at the base of the atmosphere. The strong agreement in the low resolution run to the 10/10 run gives us confidence that we are capturing the flame physics properly.

Zoom In Zoom Out Reset image size

5.2. Effects of our Approximations

The results above all used a boosting of 10/10. To see how the results are sensitive to this boosting, we ran a simulation with a reduced boosting of 5/5. This is shown in Figure 7. The results look qualitatively the same—a laterally propagating flame develops that is lifted off of the bottom of the fuel layer. The flame has not moved out as far as in the 10/10 simulation, simply because there is less energy release, but we expect that if we were to run this out twice as long, the flame would have advanced to the position seen in our 10/10 runs. The ash is also not as evolved to as high of an $\bar{A}$. The good agreement in the structure of the flame seen with the lower boosting gives us confidence that the overall aspects of the flame structure and acceleration we are seeing are robust to the approximations we make.

Zoom In Zoom Out Reset image size

We also considered the effect of the network size. Figure 8 shows a comparison of the 10/10 boosting run with the standard 13-isotope aprox13 network and the reduced 7-isotope iso7 network. We see that the flame in the aprox13 case is slightly more advanced and has a higher $\bar{A}$ than the iso7 case. We see in the next section that these two networks give largely the same flame speed. We only ran the simulation with the iso7 reaction network out to 15 ms, and we observe this was sufficient to establish a steady flame speed for comparison to our 10/10 simulation with the aprox13 network. The reduced network size saves a lot of memory, which will be useful when we transition to 3D simulations.

Zoom In Zoom Out Reset image size

The final approximation to explore is our choice of rotation rate. We ran the 10/10–1000 Hz simulation for 12.6 ms (shorter than the 20 ms we used for the other runs). We also used a larger domain, $1.8432\times {10}^{5}\,\mathrm{cm}\times 3.072\times {10}^{4}\,\mathrm{cm}$ to account for the larger Rossby radius. Figure 9 shows the flame structure. Overall it looks much like the faster rotator. In the next section, we explore the effect of rotation on the flame acceleration.

Zoom In Zoom Out Reset image size

5.3. Flame Propagation

To measure the propagation rate of the burning front, we first collapse our nuclear energy generation rate (${\dot{e}}_{\mathrm{nuc}}$) data at each time into a 1D radial profile by averaging over the vertical coordinate. We then take the peak ${\dot{e}}_{\mathrm{nuc}}$ value across all profiles to provide a fixed reference point. We define the position of the flame front to be the location to the right of the hottest part of the flame where ${\dot{e}}_{\mathrm{nuc}}$ first drops to $\lt 0.1 \%$ of the global maximum. This corresponds roughly to the leading edge of the burning region. Averaging over the vertical coordinate helps to reduce sensitivity to localized fluid motions, as does tracking the 0.1% contour rather than a local maximum. As we see in Figure 10, the flame settles into a state of steady propagation after an initial transient period spanning ∼3 ms. The position data here are well fitted by a linear function of time, and the resulting slope gives the velocity of the flame front.

Zoom In Zoom Out Reset image size

Table 4 gives the flame speed measured in each multidimensional simulation. Uncertainties in the velocities are the standard errors of the slopes returned by the least-squares fits, assuming uncertainties in front position on the order of the coarse grid size. The 2D flames propagate at speeds about an order of magnitude faster than their 1D counterparts (Table 1). The increase in flame speed is likely a product of the larger flame surface area and hydrodynamical effects such as turbulence, wrinkling, and convective cycles, which bring cooler fuel from ahead of the front into the hottest part of the burning region.

Table 4.  Flame Speeds Measured in 2D Calculations

Run	sfront (km s−1)
10/10	9.18 ± 0.03
5/5	4.00 ± 0.01
10/10-iso7	7.56 ± 0.02
10/10-lowres	9.33 ± 0.04
10/10–1000 Hz	16.5 ± 0.21

Download table as:  ASCIITypeset image

All of the various approximations have the expected effects. The flame in the 5/5 boost run propagates at about half the speed of the 10/10 flame. This is consistent with the 1D laminar flames, and is the behavior predicted by Equation (16). The 1000 Hz run goes nearly 2 times faster than the 10/10 2000 Hz run, as expected from the inverse relation between rotation rate and the ratio of burning front area to scale height (and as observed in Cavecchi et al. 2013). This confirms that the role of the Coriolis force is to limit the rate of flame spreading by the anticipated geometrical/hydrodynamical effect (${s}_{\mathrm{front}}={s}_{L}* {L}_{R}/H$). Reducing the resolution had minimal impact on the flame speed—the low resolution line is right on top of the fit for the standard run in Figure 10, with their slopes differing by only a few percent. Using a smaller network also produces similar behavior to the standard run, although there is a small reduction in speed owing to less energetic burning.

Cavecchi et al. (2013) ran simulations at rotation rates of 450 Hz exploring a range of opacities. Accounting for the boosting, we would expect to see our 10/10 flame propagate at ∼2.25 times the speed of these 450 Hz runs, but instead the speeds range from 4.5 to 7 times those reported in Cavecchi et al. (2013). The additional boost may come from cyclic fluid motions that would have been underresolved in the 450 Hz runs due to the enforcement of vertical hydrostatic equilibrium (we already saw that hydrodynamic effects can accelerate the burning front in our comparison with the 1D flames). The inclusion of reaction chains beyond triple-α may also increase energy generation and consequently the propagation speed.

5.4. Entrainment, Flow Features

We explore the baroclinic instability by looking at the magnitude of the baroclinicity, calculated as

$$\begin{eqnarray}&&{\boldsymbol{\psi }}=\displaystyle \frac{1}{{\rho }^{2}}{\boldsymbol{\nabla }}p\times {\boldsymbol{\nabla }}\rho ,\end{eqnarray} \tag{ 28 }$$

and shows the misalignment of the local density and pressure gradients (we also explored this in Malone et al. 2014). As we are considering a 2D system, the component of the baroclinicity we consider here (out of the plane) reflects the misalignment of the fields in the plane of the simulation. Figure 11 shows that the baroclinicity peaks along the flame front. This baroclinicity generates vorticity, which in turn entrains material along the surface of the flame (as seen in Cavecchi et al. 2013). In 3D, this same vorticity should perturb the flame front and further increase the flame speed (Cavecchi & Spitkovsky 2019).

Zoom In Zoom Out Reset image size

Vortical motions are also evident in a direct visualization of the velocity field, as seen in Figure 12. We observe turbulence in the wake of the flame and in the unburnt fuel ahead of the front, while a convective cycle sets up near the hottest part of the burning region. The cycle draws the cooler fluid near the interface into the center of the flame and drives hot ash out toward the instability at its surface, helping to facilitate the flame spreading. Convective mixing should still be important in 3D, but we would expect much more complicated flow patterns, with a greater contribution from small-scale features (Zingale et al. 2015).

Zoom In Zoom Out Reset image size

To further demonstrate the relationships between different properties of the flame, in Figure 13 we present some phase plots of the energy generation rate. The phase plot of the energy generation rate as a function of the x- and y-velocities shows that energy is preferentially generated in regions with negative x-velocity. This most likely corresponds to the burning along the underside of the flame front, where fuel has been entrained along the surface of the flame and so is moving in the direction opposite to the direction the flame is propagating in. This is further supported by the phase plot of the energy generation rate as a function of the x-velocity and the density, which shows that the energy generation rate peaks at $\rho \sim 3\times {10}^{5}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, at the base of the flame. The energy generation rate also depends strongly on density.

Zoom In Zoom Out Reset image size

In the first plot, there are "loops" of points of similar ${\dot{e}}_{\mathrm{nuc}}$ in the outer edges of the plot. These are likely to correspond to the vortices that appear within the flame, where the fluid moves in a circular motion in u − v phase space.

Figure 14 shows the energy generation rate as a function of the temperature and $\bar{A}$. The energy generation rate peaks at low $\bar{A}$, where there is a high fraction of unburnt material. The burning raises the temperature of this material, such that the peak temperature coincides with the peak energy generation rate. The material then cools again as the burning converts the fuel into ashes, increasing $\bar{A}$ and reducing the energy generation rate as there is less available fuel. Along the base of the plot we see cool unburnt fluid and ashes. In the center of the plot there are thin "trails" in phase space, which could correspond to less common reaction pathways in the reaction network.

Zoom In Zoom Out Reset image size

We have shown the results of our fully hydrodynamical, multidimensional simulations of flame propagation through the atmosphere of a neutron star. By using the fully compressible hydrodynamics equations, we are able to capture the vertical dynamics of the system. To accurately model the flame propagation, it is necessary to sufficiently resolve the scale height of the atmosphere. It is also important that there is a thin transition between the underlying neutron star and the atmosphere, to ensure that the peak temperature occurs at the correct base density.

In order to satisfy these requirements while allowing the simulations to be computationally feasible, we used several approximations: a higher than normal rotation rate, a boosted flame, a simplified reaction network, and a 2D axisymmetric model for the flow. The first two of these were used to reduce the model's spatial and temporal scales. Using the simulation framework developed here, these both can be relaxed in the future, at the cost of more computer time. The same goes for 2D axisymmetry versus full 3D—the only difference is computer time, and our future calculations will explore the 3D evolution and compare to the 2D simulations presented here. In particular, in 3D we will be able to explore shear instabilities at the flame front. We can also capture the baroclinic instability (Cavecchi & Spitkovsky 2019) and the competition between it and shear. Larger networks are a straightforward change, and already supported in Castro using the pynucastro framework (Willcox & Zingale 2018) and JINA ReacLib rate database (Cyburt et al. 2010).

In addition to extending our models to 3D and relaxing the approximations used in this study, in the future we plan to perform a number of additional further studies. These include investigating the effects of ignition latitude and different initial models. We also plan to model mixed H/He bursts. This will require a different reaction network, and perhaps different resolution requirements. A long-term goal is to include magnetohydrodynamics in our models. Although the magnetic fields of neutron stars exhibiting X-ray bursts are relatively weak, with $B\lesssim {10}^{8}\mbox{--}{10}^{9}\,$ G (Mukherjee et al. 2015; at least compared to, e.g., magnetars, which have $B\lesssim {10}^{15}\,$ G), it has been found by Cavecchi et al. (2016) that even weak magnetic fields could have a nonnegligible effect on the flame propagation, reducing confinement due to the Coriolis force and leading to increased flame speeds. Ultimately, we wish to link our simulations to observed light curves. To do this, we will need to explore radiation transport in order to model how the burst energy propagates through the outer layers of the neutron star atmosphere.

In addition to exploring other aspects of the XRB physics, there are several changes to the algorithm used to model these XRBs we will pursue. First, as shown in Zingale et al. (2019), we have developed a fourth-order (in space and time) method for coupling hydrodynamics and reactions that should greatly improve the accuracy of the simulations. We expect that by using this new high-order algorithm we can drop a level of refinement from the simulations while still accurately modeling the evolution. We have also ported Castro to GPUs, giving an order of magnitude performance boost on nodes with both CPUs and GPUs. Flames without any boosting are currently running using the new GPU-enabled Castro and will be the focus of the next study.

Castro is open-source and freely available at http://github.com/AMReX-Astro/Castro. The problem setups used here are available in the git repo as flame and flame_wave. The work at Stony Brook was supported by DOE/Office of Nuclear Physics grant DE-FG02-87ER40317 and the SciDAC program DOE grant DE-SC0017955. M.Z. acknowledges support from the Simons Foundation. Y.C. was supported by the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Global Fellowship grant agreement No. 703916. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725, awarded through the DOE INCITE program. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

Facilities: NERSC - , OLCF. -

Software: AMReX (Zhang et al. 2019), Castro (Almgren et al. 2010), GCC (https://gcc.gnu.org/), linux (https://www.kernel.org/), matplotlib (Hunter 2007, http://matplotlib.org/), NumPy (Oliphant 2007; van der Walt et al. 2011), python (https://www.python.org/), valgrind (Nethercote & Seward 2007), VODE (Brown et al. 1989), yt (Turk et al. 2011).

To test our implementation of diffusion, we created a unit test in Castro that just uses the diffusion solver. This test uses the standard Gaussian initial conditions—the diffusion of a Gaussian remains Gaussian, with a lower amplitude and wider width. This is run with a constant diffusion coefficient. The point of this test is to demonstrate that the implementation done in a predictor-corrector form is second-order accurate.

We use a 2D axisymmetric geometry, with a reflection boundary on the symmetry axis (left) and Neumann boundaries on all other sides. In this geometry, the solution behaves like a spherical problem, and has the form:

$$\begin{eqnarray}&&T(r,t)={T}_{1}+({T}_{2}-{T}_{1}){\left(\displaystyle \frac{{t}_{0}}{t+{t}_{0}}\right)}^{3/2}{e}^{-{r}^{2}/(4D(t+{t}_{0}))},\end{eqnarray} \tag{ A1 }$$

where t₀ is a small time used to set the initial width of the Gaussian, r is the distance from the origin, and D is the diffusion coefficient (see, e.g., Swesty & Myra 2009). T₁ is the ambient temperature and T₂ is the temperature of the peak of the Gaussian. We use $D=1$, ${T}_{1}=1$, ${T}_{2}=2$, and ${t}_{0}=0.001$.

Table A1 shows the L-inf error against the analytic solution for the diffusion test problem. The tests were run in two fashions: no AMR and AMR with one level of refinement (a jump of 2×). In both cases we see second-order convergence of the error. Figure A1 shows the temperature at the end of the simulation of the $64\times 128$ base grid plus one level of refinement simulation, with the grids shown to illustrate where the refinement jump is.

Zoom In Zoom Out Reset image size

Table A1.  Convergence of Diffusion Test Problem

Base Resolution	L-inf Error	Convergence Ratea
no refinement
32	$5.995\,\times \,{10}^{-3}$	2.04
64	$1.477\,\times \,{10}^{-3}$	2.01
128	$3.664\,\times \,{10}^{-4}$	2.00
256	$9.139\,\times \,{10}^{-5}$
one level of refinement
32	$1.480\,\times \,{10}^{-3}$	2.01
64	$3.685\,\times \,{10}^{-4}$	1.97
128	$9.406\,\times \,{10}^{-5}$	1.94
256	$2.446\,\times \,{10}^{-5}$

Note.

^aConvergence, n, is defined between two successive resolutions as $n={\mathrm{log}}_{2}({e}_{2h}/{e}_{h})$, where e_2h is the error for the coarser resolution run and e_h is the error for the finer resolution run.

Download table as:  ASCIITypeset image