Dynamics of Laterally Propagating Flames in X-ray Bursts. II. Realistic Burning & Rotation
A. Harpole, N. M. Ford, K. Eiden, M. Zingale, D. E. Willcox, Y. Cavecchi, M. P. Katz
DD RAFT VERSION F EBRUARY
2, 2021Typeset using L A TEX preprint style in AASTeX63
Dynamics of Laterally Propagating Flames in X-ray Bursts. II. Realistic Burning & Rotation
A. H
ARPOLE , N. M. F
ORD , K. E
IDEN ,
3, 1
M. Z
INGALE , D. E. W
ILLCOX , Y. C
AVECCHI , AND
M. P. K
ATZ Dept. of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800 Lawrence Berkeley National Laboratory, Berkeley, CA Dept. of Astronomy, University of California, Berkeley, CA 94720-3411 Universidad Nacional Aut´onoma de M´exico, Instituto de Astronom´ıa, Ciudad Universitaria, CDMX 04510, Mexico NVIDIA Corporation
ABSTRACTWe continue to investigate two-dimensional laterally propagating flames in type I X-raybursts using fully compressible hydrodynamics simulations. In the current study we relaxprevious approximations where we artificially boosted the flames. We now use more physi-cally realistic reaction rates, thermal conductivities, and rotation rates, exploring the effectsof neutron star rotation rate and thermal structure on the flame. We find that at lower rotationrates the flame becomes harder to ignite, whereas at higher rotation rates the nuclear burn-ing is enhanced by increased confinement from the Coriolis force and the flame propagatessteadily. At higher crustal temperatures, the flame moves more quickly and accelerates as itpropagates through the atmosphere. If the temperature is too high, instead of a flame prop-agating across the surface the entire atmosphere burns steadily. All of the software used forthese simulations is freely available.
Keywords:
X-ray bursts (1814), Nucleosynthesis (1131), Hydrodynamical simulations (767),Hydrodynamics (1963), Neutron stars (1108), Open source software (1866),Computational methods (1965) INTRODUCTIONConsiderable evidence suggests that ignition in an X-ray burst (XRB) starts in a localized region andthen spreads across the surface of the neutron star (Bhattacharyya & Strohmayer 2007; Chakraborty &Bhattacharyya 2014). We continue our study of flame spreading through fully compressible hydrodynamicssimulations of the flame. Building on our previous paper (Eiden et al. 2020), we relax the approximationswe made previously (artificially boosting the speed of the flame in order to reduce the computational cost)and explore how the flame properties depend on rotation rate and the thermal structure of the neutron star.This new set of realistic simulations is possible because of the work done to offload our simulation code,
Castro (Almgren et al. 2020), to GPUs, where it runs significantly faster.
Corresponding author: A. [email protected] a r X i v : . [ a s t r o - ph . H E ] J a n H ARPOLE ET AL .We investigate the effect of rotation rate on the flame. With the exception of IGR J17480-2446 (Altami-rano et al. 2010, spinning at
11 Hz ), most observations of XRBs come from sources with rotation rates of −
600 Hz (Bilous & Watts 2019; Galloway et al. 2020). There are a number of factors that could explainthis lack of observations below
200 Hz . It could be that there is some physical process which inhibits theflame ignition and/or spread at lower rotation rates. It could be that bursts at lower rotation rates are smallerin amplitude and therefore more difficult to detect. It could be that it does not have anything to do withthe flame at all, but that neutron stars in accreting low mass X-ray binaries rarely have rotation rates below
200 Hz .Previous studies have found that rotation can have a significant effect on the flame’s propagation. Asthe rotation rate increases, the Coriolis force whips the spreading flame up into a hurricane-like structure(Spitkovsky et al. 2002; Cavecchi et al. 2013). The stronger Coriolis force leads to greater confinement ofthe hot accreted matter, leading to easier ignition of the flame (Cavecchi et al. 2015).The temperature structure of the accreted fuel layer can also affect the flame propagation. Timmes (2000)showed that laminar helium flames have higher speeds when moving into hotter upstream fuel. It has beensuggested that crustal heating may be stronger at lower accretion rates and weaker at higher accretion rates,due to the effect of neutrino losses (Cumming et al. 2006; Johnston et al. 2019). On the other hand, at veryhigh accretion rates the atmosphere is so heated that it simmers in place rather than forming a propagatingflame (Fujimoto et al. 1981; Bildsten 1998; Keek et al. 2009). A shallow heating mechanism of as yetunknown origin has been found necessary to reproduce observed properties of XRBs in 1D simulations(Deibel et al. 2015; Turlione et al. 2015; Keek & Heger 2017). In our models, we keep the crust at aconstant temperature, so by increasing this temperature we can effectively increase the crustal heating,shallow heating and/or mimic the effects of accretion-induced heating.In the following sections, we conduct a series of simulations at various rotation rates and crustal temper-atures to investigate their effects on the flame. We find that at lower rotation rates, the flame itself becomesharder to ignite. At higher rotation rates, nuclear burning is enhanced and the flame propagates steadily.At higher crustal temperatures, burning is greatly enhanced and the flame accelerates as it propagates. Wediscuss the implications that this may have for burst physics and observations. NUMERICAL APPROACHWe use the
Castro hydrodynamics code (Almgren et al. 2010; Almgren et al. 2020) and the simulationframework introduced in Eiden et al. (2020). The current simulations are all performed in a two-dimensionalaxisymmetric geometry. For these axisymmetric simulations, we add an additional geometric source termfrom Bernard-Champmartin et al. (2012) that captures the effects of the divergence of the flux operating onthe azimuthal unit vector. This term is a small correction, but was missing from our previous simulations.The simulation framework initializes a fuel layer in hydrostatic equilibrium, laterally blending a hot modelon the left side of the domain (the coordinate origin) and a cool model on the right. The initial temperaturegradient between the hot and cool fluids drives a laterally propagating flame through the cool fuel. In ouroriginal set of calculations (Eiden et al. 2020), in order to make the simulations computationally feasible weartificially boosted the flame speed by adjusting the conductivity and reaction rate to produce a flame moving5–10 × faster than the nominal laminar flame speed. We also used high rotation rates ( ≥ ) to reducethe lateral lengthscale at which the Coriolis force balances the lateral flame spreading in order to reduce thesize of the simulation domain. The port of Castro to GPUs (Katz et al. 2020) significantly improved itsoverall performance, enabling us to run these new simulations without the previous approximations whilecontinuing to resolve the burning front. For these simulations, we no longer boost the flame speed—the true
ATERAL F LAME D YNAMICS
II 3 z [cm]10 T [ K ] F500F500_2E8F500_3E8F500_4E8 0 500 1000 1500 2000 2500 3000 3500 4000z [cm]10 [ g c m ] Figure 1.
Initial temperature structure (left panel) and density structure (right panel) as a function of height in the“hot” state. conductivities (taken from Timmes 2000) and reaction rates are used. We are also able to use slower, morephysically realistic rotation rates. We continue to use a 13-isotope α -chain to describe the helium burning.The initial model is set up in the same fashion as described in Eiden et al. (2020). In particular, we createa “hot” and “cool” hydrostatic model representing the ash and fuel states and blend the two models laterallyto create a hot region near the origin of the coordinates and a smooth transition to the cooler region at largerradii. The cool initial model is characterized by three temperatures: T star is the isothermal temperature ofthe underlying neutron star, T hi is the temperature at the base of the fuel layer, and T lo is the minimumtemperature of the atmosphere. The atmosphere structure is isentropic as it falls from T hi down to T lo . Forthe hot model, we replace T hi with T hi + δT . In the calculations presented here, we explore the structureof the initial models by varying these parameters. All models have the same peak temperature in the hotmodel, T hi + δT .For the current simulations, we explore a variety of initial rotation rate and temperature conditions for theflame. The main parameters describing the models and the names by which we shall refer to them in thispaper are provided in Table 1. Figure 1 shows the temperature and density structure for our hot models andFigure 2 shows the temperature and density structure for the cool models. SIMULATIONS AND RESULTSWe present six simulations in total, summarized in Table 1. These simulations encompass three differentrotation rates:
250 Hz ,
500 Hz , and , and for the
500 Hz run, four different temperature profiles.In the following subsections, we look at how the flame properties depend on the model parameters. Allsimulations are run in a domain of . × cm × . × cm with a coarse grid of × zones and two levels of refinement (the first level refining the resolution by factor of four, and the secondby a factor of two again). This gives a fine-grid resolution of
20 cm . In these simulations, refinement iscarried out in all zones within the atmosphere with density ρ > . × g cm − . We use an axisymmetriccoordinate system, with the horizontal r -direction pointing along the surface of the star and the vertical z -direction pointing perpendicular to the surface. H ARPOLE ET AL . z [cm]10 T [ K ] F500F500_2E8F500_3E8F500_4E8 0 500 1000 1500 2000 2500 3000 3500 4000z [cm]10 [ g c m ] Figure 2.
Initial temperature structure (left panel) and density structure (right panel) as a function of height in the“cool” state.
Table 1.
Rotation rate and temperature properties of the simulations. In theleft-hand column we list the names we shall use to refer to each simulationthroughout this paper.run Rotation Rate (Hz) δT (K) T hi (K) T star (K) T lo (K) F1000 . × × × F500
500 1 . × × × F500 2E8
500 1 . × × × × F500 3E8
500 1 . × × × × F500 4E8
500 10 × × × F250
250 1 . × × × For some of our analysis, we would like to have a means of estimating the temperature ( T ) and nuclearenergy generation rate ( ˙ e nuc ) in the burning region of each simulation. For this purpose, we define themass-weighted average (cid:104) Q (cid:105) w of some quantity Q to be (cid:104) Q (cid:105) w ≡ (cid:80) c i m ( c i ) Q ( c i ) (cid:80) c i m ( c i ) ; c i ∈ C ( Q ) . (1)Here, C ( Q ) is the set of grid cells with Q values in the top percentile, Q ( c i ) is the value of Q in cell c i ,and m ( c i ) is the total mass in cell c i . Using (cid:104) Q (cid:105) w instead of simply taking the maximum of the quantityacross the entire simulation domain allows us to track changes over the domain as a whole rather than at asingle localized point. This will therefore be a better reflection of the overall behavior of the flame ratherthan of a single localized fluctuation. Figure 3 shows (cid:104) T (cid:105) w and (cid:104) ˙ e nuc (cid:105) w as functions of time for the subsetof our runs that achieve a propagating flame. This figure is referenced throughout the subsequent sections. ATERAL F LAME D YNAMICS
II 5 T w [ K ] F500F500_2E8F500_3E8F1000 e nu c w [ e r g / g / s ] Figure 3.
Estimates of temperature (left panel) and nuclear energy generation rate (right panel) in the burning regionas functions of time. The quantities on the vertical axes are the mass-weighted averages defined in Equation 1.
Effect of Rotation Rate on Flame Structure
We run three models (
F250 , F500 , and
F1000 ) with the same initial model in terms of temperature butdiffering rotation rates. We saw in Eiden et al. (2020) that increasing the neutron star rotation rate reducesthe horizontal lengthscale of the flame. An estimate of this lengthscale is given by the Rossby radius ofdeformation, L R . The Rossby radius may be thought of as the scale over which the balance between theCoriolis force and horizontal pressure gradient becomes important, and is approximately given by L R ≈ √ gH Ω , (2)where g is the gravitational acceleration, H is the atmospheric scale height, and Ω is the neutron starrotation rate. In Figure 4 and Figure 5, we use ˙ e nuc measured at
50 ms and
100 ms to discern the horizontalextent of the flame at different rotation rates. Taking the edge at greatest radius of the bright teal/greenregion where the most significant energy generation is occurring as the leading edge of the flame in eachplot, we see that the horizontal extent of the flame (
F1000 ) appears to be reduced compared to thelower rotation
500 Hz run (
F500 ). From Equation 2, we can see that increasing the rotation rate from
500 Hz to should decrease L R by a factor of two, and that the greater confinement from the Coriolis forceshould reduce the horizontal extent of the flame by a similar factor. However, the Rossby radius is onlyan approximate measure of this horizontal lengthscale, and in our simulations we see that this scaling doesnot work so well for all rotation rates. The simulations seem to follow the theoretical scaling described inEquation 2 more closely at higher rotation rates ( and higher), based on the results of Eiden et al.(2020).The F500 and
F1000 runs both qualitatively resemble the flame structure in Eiden et al. (2020) — alaterally propagating flame that is lifted off of the bottom of the fuel layer — but they differ in their burningstructures. Figures 6 and 7 show time series of the mean molecular weight, ¯ A , for the F500 and
F1000 runs. Compared to those in Eiden et al. (2020), ashes behind the flame do not reach as high atomic weights.This is not surprising, since those previous runs artificially boosted the reaction rates. Comparing these H
ARPOLE ET AL . z ( c m ) × r (cm) z ( c m ) × e nu c ( e r gg s ) Figure 4.
Time series of the
500 Hz run
F500 showing the nuclear energy generation rate, ˙ e nuc . z ( c m ) × r (cm) z ( c m ) × e nu c ( e r gg s ) Figure 5.
Time series of the run
F1000 showing the nuclear energy generation rate, ˙ e nuc . z ( c m ) × r (cm) z ( c m ) × A b a r Figure 6.
Time series of the
500 Hz run
F500 showing the mean molecular weight, ¯ A . two new runs, the burning is much more evolved for the higher rotation rate, and the ash is actually ableto move ahead of the flame front (visible in the Figure 7
100 ms snapshot). We believe that this is becausethe increased rotation better confines the initial perturbation and subsequent expansion from the burning,increasing the temperature and density in the flame front such that the reaction rate increases, which allowsthe reactions to progress further. The ˙ e nuc plots in Figure 5 also support this interpretation, with the regionof the flame front nearest to the crust in the F1000 run reaching higher ˙ e nuc values than for the F500 runin Figure 4. In contrast to
F500 and
F1000 , the lowest rotation run —
F250 — failed to ignite. Thelack of ignition for
F250 also aligns with the reasoning given above, with the lower rotation in this casepotentially leading to insufficient confinement such that the temperature and density required for ignitionis not achieved. In this scenario, another source of confinement (e.g. magnetic fields, see Cavecchi et al.
ATERAL F LAME D YNAMICS
II 7 z ( c m ) × r (cm) z ( c m ) × A b a r Figure 7.
Time series of the run
F1000 showing the mean molecular weight, ¯ A . (2016)) would need to take over at lower rotation rates to allow a burst to occur, at least for the initial modelused here. Given that the size of our domain is ∼ L R for F250 (using Equation 2), it is also possiblethat we simply cannot confine the flame sufficiently with our current domain width. We see in Figure 11(discussed further in Section 3.2) that the
F500 flame took longer to achieve steady propagation than the
F1000 flame. It may therefore also be that we did not run our simulation for long enough to see the
F250 flame achieve the conditions required for ignition and steady propagation.Burning in the
F500 and
F1000 runs is concentrated in a dense region with circular motion. In Figure 8,which compares the horizontal r -velocity u , density ρ and the nuclear energy generation rate ˙ e nuc for the F250 , F500 , and
F1000 runs, most of the burning for each of the simulations occurs in a high densityregion ρ > × g cm − . The fluid in this dense, high energy generation region undergoes vorticalmotion, shown in the Figure 9 phase plots comparing u , the vertical z -velocity v and ˙ e nuc . This mostlikely corresponds to the leading edge of the flame where fresh fuel is being entrained. This feature is notdeveloped in the
250 Hz flame in Figure 9 (left panel); it could potentially develop at later times (past thepoint at which we terminated our simulation), or the burning could just fizzle out and the flame fail to igniteentirely.The mean molecular weight ¯ A within each of our simulations seems to grow along defined tracks con-fined to certain temperatures T , as shown in the Figure 10 phase plots. We believe that the tracks in the plotcorrespond to different burning trajectories in phase space resulting from different thermodynamic condi-tions. Comparing Figure 10 to Eiden et al. (2020), these tracks are much more neat and clearly defined. The“messiness” of the tracks may be dependent on how mixed the flame interior is. Since these new simulationsare un-boosted, they may be inherently less mixed than those in Eiden et al. (2020). F1000 aligns with thisinterpretation: its ¯ A tracks are somewhat disrupted compared to the slower rotation runs, possibly due tothe more vigorous mixing of the vortex at the flame front. Comparing the different runs, we also see thatas the rotation rate increases, so does the peak temperature. This makes sense if higher rotation leads to amore concentrated, intense vortex near the flame front. It also agrees with our earlier interpretation of theenhanced burning seen in Figure 7 for F1000 .3.2.
Effect of Rotation Rate on Flame Propagation
For the purpose of measuring the flame propagation speed and acceleration, we track the position ofeach of our flames as a function of time. We define the position in terms of a specific value of the energygeneration rate, ˙ e nuc , as we did in Eiden et al. (2020). To recapitulate: we first reduce the 2D ˙ e nuc data foreach simulation run to a set of 1D radial profiles by averaging over the vertical coordinate. After averaging,we take our reference ˙ e nuc value to be some fraction of the global ˙ e nuc maximum across all of these profiles. H ARPOLE ET AL . Figure 8.
Phase plot of the horizontal r -velocity u , density ρ and the nuclear energy generation rate ˙ e nuc for the 250,500 and runs ( F250 , F500 and
F1000 ) at t = 100 ms . The slightly lower ˙ e nuc values along the u = 0 axisare most likely a numerical artifact related to the density gradient setup and finite resolution of these simulations. Figure 9.
Phase plot of the horizontal r -velocity u , vertical z -velocity v and the nuclear energy generation rate ˙ e nuc for the 250, 500 and runs ( F250 , F500 and
F1000 ) at t = 100 ms . Since the flames in our simulations propagate in the positive horizontal direction, we then search the regionof each profile at greater radius than the local ˙ e nuc maximum for the point where the ˙ e nuc first drops belowthis reference value. This point gives us the location of our flame front.In Eiden et al. (2020), we used . of the global ˙ e nuc maximum for our reference value. For the hightemperature unboosted flames, however, we found that the ˙ e nuc profiles failed to reach that small a valueacross the domain at most times, which prevented us from obtaining reliable position measurements. Wetherefore use of the global ˙ e nuc maximum in this paper rather than . . This is still sufficiently smallthat our measurements are not overly sensitive to turbulence and other local fluid motions (the issue withsimply tracking the local maximum), but allows us to avoid the pitfall encountered by the . metric.Figure 11 gives the radial position of the flame front as a function of time for the F500 and
F1000 runs(blue and orange, respectively) to show the dependence on rotation rate. In Eiden et al. (2020), we applied alinear least-squares fit to the flame front position as a function of time to estimate the propagation velocity.As some of the flames in this set of simulation runs exhibit significant acceleration, for this study we instead
ATERAL F LAME D YNAMICS
II 9
Figure 10.
Phase plot of the mean molecular weight ¯ A , temperature T and the nuclear energy generation rate ˙ e nuc forthe 250, 500 and runs ( F250 , F500 and
F1000 ) at t = 100 ms . Table 2.
Flame speed and acceleration values measured for eachsimulation. v is the flame velocity in the + r direction at t =0 ms and a is the acceleration of the flame. The initial flamevelocities and accelerations are derived from a quadratic least-squares fit to each of the datasets for times t (cid:38)
20 ms . Usingthese fit parameters we calculate the velocity at t = 50 ms , v ,at which point the flames have reached steady propagation.run v (km s − ) a (km s − ) v (km s − ) F1000 . ± . − . ± .
13 3 . ± . F500 . ± .
013 5 . ± .
19 3 . ± . F500 2E8 . ± .
015 1 . ± .
25 3 . ± . F500 3E8 − . ± .
100 357 . ± .
15 12 . ± . fit the data with a quadratic curve of the form r ( t ) = 12 at + v t + r , (3)where the parameter a is the acceleration of the flame, v is the velocity at t = 0 ms , and r is the flamefront position at t = 0 ms . We do not include the data points with t (cid:46)
20 ms when performing the fit, sincethese correspond to the initial transient period before the flame has begun to propagate steadily. The valuesof a and v for the full suite of simulation runs are provided in Table 2. Note that v is only a parameter thatmay be used to calculate the flame speed at an arbitrary time. It is not an estimate of the true initial velocityof the flame, since the flame has not achieved ignition yet at t = 0 ms . We use the fit parameters to calculate0 H ARPOLE ET AL ..
20 ms when performing the fit, sincethese correspond to the initial transient period before the flame has begun to propagate steadily. The valuesof a and v for the full suite of simulation runs are provided in Table 2. Note that v is only a parameter thatmay be used to calculate the flame speed at an arbitrary time. It is not an estimate of the true initial velocityof the flame, since the flame has not achieved ignition yet at t = 0 ms . We use the fit parameters to calculate0 H ARPOLE ET AL .. r ( c m ) F500F1000
Figure 11.
Flame front position vs. time for the standard ( K )
500 Hz and runs (
F500 and
F1000 ).The dashed lines show quadratic least-squares fits to the data for t (cid:38)
20 ms . the flame speeds at t = 50 ms (when the flame has reached steady propagation), given in the fourth columnof Table 2.We found in Eiden et al. (2020) that the flame speed, s , obeys the following relation: s ∝∼ L R ( ˙ e nuc ) ∝ ( ˙ e nuc ) Ω , (4)where ˙ e nuc is the specific nuclear energy generation rate and L R is the Rossby length described by Equation2. This finding was consistent with the results of Cavecchi et al. (2013), who noted that the flame speeds intheir simulations scaled with the Rossby length. As seen in Figure 11, there is no obvious inverse scalingof the flame speed with rotation rate in this set of runs. We observed earlier, however, that nuclear reactionsprogress more quickly at higher rotation rate. This results in a higher ˙ e nuc — up to three to four times highernear the burning front after the flame ignites (see Figure 3) — which may counteract the reduction in flamespeed from the increased Coriolis confinement. We also observe that F500 accelerates faster than
F1000 ,which appears to experience a small deceleration at early times. This disparity may be a direct result of thedifference in Coriolis force. 3.3.
Effect of Temperature on Flame Structure
To explore the effect of different initial temperature configurations, we run four simulations fixed at arotation rate of
500 Hz with temperatures as shown in Table 1. For all the
500 Hz simulations (with theexception of the coolest run,
F500 ), we set T star = T hi , scaling δT accordingly to maintain a consistent valueof T hi + δT . If we let T star < T hi , the cooler neutron star surface might act as a heat sink and siphon away ATERAL F LAME D YNAMICS
II 11 z ( c m ) × F500 z ( c m ) × F500_2E8 r (cm) z ( c m ) × F500_3E8 A b a r Figure 12.
Comparison of ¯ A for 3 different
500 Hz models with neutron star temperatures T star of K , × K and × K ( F500 , F500 2e8 , and
F500 3e8 , respectively), and resulting envelope structures. Each flame is shownat
70 ms . energy that would otherwise go into heating the burning layer. By choosing T star = T hi , we can effectivelyexplore simulations with greater surface heating. There are several physically distinct mechanisms thatcould produce an increased temperature at the crust: crustal heating, some other form of shallow heatingor accretion-induced heating. In these simulations, we do not model the mechanism producing the heatingeffect, just the effect itself, so we do not distinguish between which of these mechanisms cause the heating.Figure 12 shows ¯ A for three
500 Hz simulations with different initial temperature structures ( T star ≤ × K ) at t = 70 ms . We do not plot F500 4E8 here because it fails to form a clear burning front.
F500 3E8 (Figure 12, bottom panel) — the hottest run to form a clear burning front — has a faster propagating flame(this will be discussed further in Section 3.4). It also reaches slightly higher ¯ A values than the two coolerruns. The Figure 18 ¯ A - T phase plots of F500 (left) and
F500 3E8 (middle) also reflect these ¯ A features,with F500 3E8 reaching slightly higher ¯ A values. F500 3E8 also reaches higher ˙ e nuc values at the low endof the temperature range ( < . × K ). There appear to be more causally connected regions across arange of ¯ A at low temperatures for F500 3E8 than for
F500 , suggesting that the higher ˙ e nuc for F500 3E8 generates burning in certain burning trajectories that are not present in the cooler
F500 run. Note that,although
F500 3E8 is the hottest run with a modified initial temperature configuration to form a distinctflame front, the highest rotation
F1000 flame actually reaches higher temperatures (see Figure 3, left panel)as well as higher ¯ A (see Figure 7).In contrast to the models with T star ≤ × K , F500 4E8 is so hot that the organized flame structureis lost. The fuel layer undergoes stable burning and essentially simmers in place, rather than forming aflame. This process has been observed in neutron stars with accretion rates that exceed the Eddington limit(Fujimoto et al. 1981; Bildsten 1998; Keek et al. 2009); in this model we have essentially recreated theconditions found in such stars with high accretion rates. This model burns so strongly that it is only run outto
40 ms . After an initial period when the burning moves across the domain, residual burning continues andeventually ignites the entire fuel layer at late times, as shown in Figures 13 and 14 for three snapshots takenin the last ten seconds of the simulation. In Figure 13, ˙ e nuc reaches values of − erg g − s − acrossthe domain and at heights up to ∼ . × cm . There is still significant burning occurring even higher,with ˙ e nuc reaching ∼ erg g − s − at heights up to ∼ . × cm . For comparison, the next hottest run( F500 3E8 ) only reaches maximum ˙ e nuc values on the order of erg g − s − (see Figure 3, right panel),2 H ARPOLE ET AL .even at the latest timesteps ( ∼
70 ms , when the flame is most developed). Significant ˙ e nuc values for allruns other than F500 4E8 are confined to the flame front, rather than spanning the entire domain.Figure 14 shows ¯ A for F500 4E8 . Again, burning extends across the domain and high into the atmosphereand fuel layer, lacking the characteristic flame structure shown in ¯ A plots for the lower temperature
500 Hz runs (see Figure 12). A distinct mass of material appears to have broken out of the atmosphere near theaxis of symmetry. The atmosphere edge is located at ∼ . × cm . A similar effect is also visible (to alesser extent) in the F500 3E8 plot in Figure 12, with a faint haze of material rising above the flame near theaxis of symmetry. However, this is likely to be a numerical artefact rather than a true physical effect, and aresult of the boundary conditions imposed at the axis of symmetry for these simulations.
F500 4E8 clearlyreaches much higher ¯ A values across the domain (especially at the latest snapshot, t = 40 ms ) compared toall the other runs described in this paper (see Figures 6, 7 and 12).Though F500 4E8 does not form a distinct burning front, it does achieve greatly enhanced burning, asshown in the Figure 18 ¯ A - T phase plot (right panel). Not only is F500 4E8 significantly hotter overallcompared to all other runs, but there is a large causally connected region across a wide range of ¯ A . Thisindicates that the hotter temperature has facilitated significant burning in burning trajectories that were notfavored at the lower temperatures. The burning trajectories are also very disrupted for F500 4E8 comparedto the cooler runs, suggesting that the hotter temperature leads to more vigorous mixing. Indeed, this appearsto be the case looking at Figure 14 and the right panel of Figure 17 (discussed further in section 3.4). Thedisrupted burning trajectories resemble those found in Figure 10 for the highest rotation
F1000 run, thoughthey are even more dramatically disrupted for
F500 4E8 . F500 4E8 is also noticeably hotter than
F1000 ,even though it is run for significantly less time ( t = 40 ms vs t = 100 ms ). Although the F500 4E8 run is clearly a special case in that it develops steady burning across the domain rather than a propagatingflame, enhanced burning in this hottest run aligns with results from the other simulations with differinginitial temperature structures. As T star is increased, the overall behavior and propagation of the flame issignificantly altered, implying that burning is very temperature-sensitive. We explore flame propagationfurther in Section 3.4. 3.4. Effect of Temperature on Flame Propagation
The method for measuring flame propagation described in Section 3.2 is applied in Figure 15 to the three
500 Hz runs with T star ≤ × K . Due to the lack of a clear burning front in F500 4E8 , we do notanalyze its propagation velocity and acceleration. As the initial T star is increased beyond ∼ × K , theflame becomes greatly accelerated. The initial flame velocities and accelerations derived from a quadraticleast-squares fit to each of the datasets, as well as the flame velocities at t = 50 ms calculated using the fitparameters, are provided in Table 2. Comparing the flame propagation at different initial temperatures, themost robust feature is the acceleration of the F500 3E8 run at t ∼
40 ms . The reason for the acceleration ofthe flame is not entirely clear. Whereas for the cooler runs, a state of steady flame propagation is achieved,for the
F500 3E8 run the flame speed continues to increase, suggesting that some instability driving theflame’s propagation persists to later times. It could be that energy released from burning begins to dominatethe flame’s propagation as it evolves, increasing the flame speed over time. Another possibility is that theincreased temperatures lead to enhanced turbulent mixing effects that pull in more fresh fuel for the flameto burn. Yet another possibility is that the higher initial T star leads to a greater average temperature in thefuel layer over time, making it easier for the flame to burn the fuel and propagate.In both the u - ρ phase plot in Figure 16 and the u - v phase plot in Figure 17, the horizontal u -velocityfor the F500 4e8 (right panel) burning region is slightly larger than that of the cooler runs. Additionally,
ATERAL F LAME D YNAMICS
II 13 z ( c m ) ×10
30 ms z ( c m ) ×10
35 ms z ( c m ) ×10
40 ms e nu c ( e r gg s ) Figure 13.
Time series of the
500 Hz run
F500 4E8 , with T star = 4 × K , showing ˙ e nuc . This model burns sostrongly that it is only run out to
40 ms ; the snapshots shown here are at T = 30 ms ,
35 ms , and
40 ms . z ( c m ) ×10
30 ms z ( c m ) ×10
35 ms z ( c m ) ×10
40 ms A b a r Figure 14.
Time series of the
500 Hz run
F500 4E8 , with T star = 4 × K , showing ¯ A . the cooler runs’ u -velocities are two to three orders of magnitude greater than the flame speeds listed inTable 2. These phase plot velocities therefore most likely correspond to vortical motion in the turbulentburning vortices rather than the propagation of a flame itself. The vertical z -velocity v in Figure 17 furthersuggests that F500 4e8 undergoes increased vortical motion, as we see that the velocity magnitude in theburning region is significantly larger than in the lower temperature runs. Also of interest in these plots isthat the coolest run shows high velocity material comprising both high and low ˙ e nuc with horizontal velocityasymmetry, both of which suggest the burning vortex in the coolest case is less well defined than at higher4 H ARPOLE ET AL ..
F500 4E8 , with T star = 4 × K , showing ¯ A . the cooler runs’ u -velocities are two to three orders of magnitude greater than the flame speeds listed inTable 2. These phase plot velocities therefore most likely correspond to vortical motion in the turbulentburning vortices rather than the propagation of a flame itself. The vertical z -velocity v in Figure 17 furthersuggests that F500 4e8 undergoes increased vortical motion, as we see that the velocity magnitude in theburning region is significantly larger than in the lower temperature runs. Also of interest in these plots isthat the coolest run shows high velocity material comprising both high and low ˙ e nuc with horizontal velocityasymmetry, both of which suggest the burning vortex in the coolest case is less well defined than at higher4 H ARPOLE ET AL .. r ( c m ) F500F500_2e8F500_3e8
Figure 15.
Flame front position vs. time for the three
500 Hz runs with T star of K , × K and × K ( F500 , F500 2e8 , and
F500 3e8 , respectively). The dashed lines show quadratic least-squares fits to the data for t (cid:38)
20 ms . Note that due to its rapid acceleration,
F500 3E8 (green) is only run out to ∼
70 ms to avoid surpassingthe domain boundary. temperatures. The coolest run thus does not seem to have developed the characteristic vortex structure atthe burning front (i.e. where ˙ e nuc is greatest) that can be clearly seen for the hotter runs at this time (
40 ms ).As can be seen in Figure 9 (which was plotted at
100 ms ), this does develop more at later times. Similar towhat we saw when comparing the runs with different rotation rates, it would therefore appear that it takeslonger for the flame to develop when T star is cooler. DISCUSSION AND CONCLUSIONSWe ran a number of simulations of laterally propagating flames in XRBs in order to explore the effects ofrotation and thermal structure. We found that increasing the rotation rate increased the energy generationrate within the flame and enhanced nuclear burning. Apart from the lowest rotation run (which failed to ig-nite), flame propagation was not noticeably impacted by rotation rate; by the time the different flame frontsreached steady propagation, they shared comparable velocities. These results are likely due to the rotation-dependent strength of the Coriolis force and its confinement of the flame balancing the enhanced nuclearburning. We explored several models with different crustal temperatures to determine what effects mecha-nisms such as higher accretion rates, crustal heating and shallow heating may have on flame propagation.We found that increasing the temperature of the crust significantly enhanced the flame propagation. This webelieve to be because a cooler crust allows heat to more efficiently be transferred away from the flame it-self, therefore reducing the flame’s temperature, slowing burning and consequently reducing its propagationspeed. At higher crustal temperatures, we saw that the inability for heat to be efficiently transported away
ATERAL F LAME D YNAMICS
II 15
Figure 16.
Phase plot of the horizontal r -velocity u , density ρ and the nuclear energy generation rate ˙ e nuc for the
500 Hz runs with T star of K , × K and × K ( F500 , F500 3e8 , and
F500 4e8 ) at t = 40 ms (the latesttime that F500 4E8 is run out to). The run with T star = 2 × K ( F500 2E8 ) closely resembles the cooler
F500 run,and is not shown here.
Figure 17.
Phase plot of the horizontal r -velocity u , vertical z -velocity v and the nuclear energy generation rate ˙ e nuc for the
500 Hz runs with T star of K , × K and × K ( F500 , F500 3e8 , and
F500 4e8 ) at t = 40 ms .The run with T star = 2 × K ( F500 2E8 ) closely resembles the cooler
F500 run, and is not shown here. from the flame front increased the flame temperature, driving unstable, accelerating flame propagation. Wesaw that if the crust temperature was too high, then instead of a flame the entire atmosphere would burnsteadily. This is reminiscent of what is seen for neutron stars with accretion rates exceeding the Eddingtonlimit.In future work, we would like to improve and expand our simulations in order to better understand theprocesses at play and to include more physics. This includes adding tracer particles to the simulations sowe can monitor the fluid motion and perform more detailed nucleosynthesis; extending our simulations to3D, which would hopefully alleviate some of the boundary effects we have observed in these simulationsbut will require significant computational resources; and exploring the resolution of our simulations moreso that we can ensure that we have resolved all of the necessary physical processes. We would also like tomodel H/He flames, as these are the sites of rp-process nucleosynthesis (Schatz et al. 2001). Initially wewill use the same reaction sequence explored in our previous convection calculations (Malone et al. 2014).6 H
ARPOLE ET AL ..
ARPOLE ET AL .. Figure 18.
Phase plot of the mean molecular weight ¯ A , temperature T and the nuclear energy generation rate ˙ e nuc forthe
500 Hz runs with T star of K , × K and × K ( F500 , F500 3e8 , and
F500 4e8 ) at t = 40 ms . Therun with T star = 2 × K ( F500 2E8 ) closely resembles the cooler
F500 run, and is not shown here.
Finally, we recently added an MHD solver to
Castro (Barrios Sazo 2020); this will allow us in the future toexplore the effects of magnetic fields on flame propagation in XRBs.ACKNOWLEDGMENTS
Castro is open-source and freely available at http://github.com/AMReX-Astro/Castro. The problemsetup used here is available in the git repo as flame wave . The work at Stony Brook was supportedby DOE/Office of Nuclear Physics grant DE-FG02-87ER40317. This material is based upon work sup-ported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific ComputingResearch and Office of Nuclear Physics, Scientific Discovery through Advanced Computing (SciDAC) pro-gram under Award Number DE-SC0017955. This research was supported by the Exascale ComputingProject (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and theNational Nuclear Security Administration. This material is also based upon work supported by the U.S.Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Departmentof Energy Computational Science Graduate Fellowship under Award Number DE-SC0021110. This workwas supported in part by the U.S. Department of Energy, Office of Science, Office of Workforce Develop-ment for Teachers and Scientists (WDTS) under the Science Undergraduate Laboratory Internship (SULI)program. MZ acknowledges support from the Simons Foundation. This research used resources of theNational Energy Research Scientific Computing Center, a DOE Office of Science User Facility supportedby 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 NationalLaboratory, which is supported by the Office of Science of the U.S. Department of Energy under ContractNo. DE-AC05-00OR22725, awarded through the DOE INCITE program. We thank NVIDIA Corporationfor the donation of a Titan X and Titan V GPU through their academic grant program. This research hasmade use of NASA’s Astrophysics Data System Bibliographic Services.
ATERAL F LAME D YNAMICS
II 17
Facilities:
NERSC, OLCF
Software:
Almgren, A., Sazo, M. B., Bell, J., et al. 2020, Journalof Open Source Software, 5, 2513,doi: 10.21105/joss.02513Almgren, A. S., Beckner, V. E., Bell, J. B., et al. 2010,ApJ, 715, 1221,doi: 10.1088/0004-637X/715/2/1221Altamirano, D., Watts, A., Kalamkar, M., et al. 2010,ATel, 2932, 1Barrios Sazo, M. G. 2020, PhD thesis, State Universityof New York at Stony BrookBernard-Champmartin, A., Braeunig, J.-P., &Ghidaglia, J.-M. 2012, Computers and Fluids, 7,doi: 10.1016/j.compfluid.2012.09.014Bhattacharyya, S., & Strohmayer, T. E. 2007, 666,L85, doi: 10.1086/521790Bildsten, L. 1998, in ASIC, Vol. 515, 419Bilous, A. V., & Watts, A. L. 2019, 245, 19,doi: 10.3847/1538-4365/ab2fe1Brown, P. N., Byrne, G. D., & Hindmarsh, A. C. 1989,SIAM J. Sci. Stat. Comput., 10, 1038Cavecchi, Y., Levin, Y., Watts, A. L., & Braithwaite, J.2016, MNRAS, 459, 1259,doi: 10.1093/mnras/stw728Cavecchi, Y., Watts, A. L., Braithwaite, J., & Levin, Y.2013, MNRAS, 434, 3526,doi: 10.1093/mnras/stt1273Cavecchi, Y., Watts, A. L., Levin, Y., & Braithwaite, J.2015, MNRAS, 448, 445,doi: 10.1093/mnras/stu2764Chakraborty, M., & Bhattacharyya, S. 2014, ApJ, 792,4, doi: 10.1088/0004-637X/792/1/4Cumming, A., Macbeth, J., Zand, J. J. M. i. T., & Page,D. 2006, The Astrophysical Journal, 646, 429,doi: 10.1086/504698Deibel, A., Cumming, A., Brown, E. F., & Page, D.2015, ApJL, 809, L31,doi: 10.1088/2041-8205/809/2/L31Eiden, K., Zingale, M., Harpole, A., et al. 2020, ApJ,894, 6, doi: 10.3847/1538-4357/ab80bc Fujimoto, M. Y., Hanawa, T., & Miyaji, S. 1981,Astrophysical Journal, 247, 267,doi: 10.1086/159034Galloway, D. K., in ’t Zand, J. J. M., Chenevez, J.,et al. 2020, arXiv e-prints, arXiv:2003.00685.https://arxiv.org/abs/2003.00685Hunter, J. D. 2007, Computing in Science and Engg.,9, 90, doi: 10.1109/MCSE.2007.55Johnston, Z., Heger, A., & Galloway, D. K. 2019,arXiv e-prints, arXiv:1909.07977.https://arxiv.org/abs/1909.07977Katz, M. P., Almgren, A., Sazo, M. B., et al. 2020, inProceedings of the International Conference forHigh Performance Computing, Networking, Storageand Analysis, SC ’20 (IEEE Press)Keek, L., & Heger, A. 2017, The AstrophysicalJournal, 842, 113, doi: 10.3847/1538-4357/aa7748Keek, L., Langer, N., et al. 2009, Astronomy &Astrophysics, 502, 871Kluyver, T., Ragan-Kelley, B., P´erez, F., et al. 2016, inPositioning and Power in Academic Publishing:Players, Agents and Agendas, 87–90,doi: 10.3233/978-1-61499-649-1-87Malone, C. M., Zingale, M., Nonaka, A., Almgren,A. S., & Bell, J. B. 2014, ApJ, 788, 115,doi: 10.1088/0004-637X/788/2/115Nethercote, N., & Seward, J. 2007, in Proceedings ofthe 28th ACM SIGPLAN Conference onProgramming Language Design andImplementation, PLDI ’07 (New York, NY, USA:ACM), 89–100, doi: 10.1145/1250734.1250746Oliphant, T. E. 2007, Computing in Science and Engg.,9, 10, doi: 10.1109/MCSE.2007.58Schatz, H., Aprahamian, A., Barnard, V., et al. 2001,Physical Review Letters, 86, 3471,doi: 10.1103/PhysRevLett.86.3471Spitkovsky, A., Levin, Y., & Ushomirsky, G. 2002,ApJ, 566, 1018, doi: 10.1086/338040Timmes, F. X. 2000, ApJ, 528, 913,doi: 10.1086/308203
ARPOLE ET AL ..