Simulating the Collapse of a Thick Accretion Disk due to a Type I X-Ray Burst from a Neutron Star

Type I X-ray bursts result from unstable thermonuclear burning of accreted matter on the surface of a neutron star (e.g., Strohmayer & Bildsten 2006; Galloway & Keek 2017). The sudden release of energy is thermalized by the neutron star atmosphere and radiated away as X-rays over a timescale of several seconds to minutes (e.g., Belian et al. 1976; Lewin et al. 1993; Galloway et al. 2008). During its peak, the luminosity of the burst significantly outshines the surrounding accretion disk, and may even exceed the local Eddington luminosity of the star (e.g., Hansen & van Horn 1975; Lewin et al. 1984; Tawara et al. 1984; Galloway et al. 2008). Following the burst, material resumes accumulating on the stellar surface, leading to subsequent X-ray bursts with a recurrence time of hours to days, depending on the metallicity and accretion rate of the infalling gas (e.g., Fujimoto et al. 1981; Bildsten 1998). As the neutron star surface is the origin of the observed X-rays, there is significant interest in using the spectra and luminosities of X-ray bursts to measure the physical size of the star, a measurement necessary to constrain the equation of state of nuclear matter (e.g., Güver et al. 2012; Kajava et al. 2014; Nättilä et al. 2017).

Between each X-ray burst cycle the neutron star system is characterized as an accretion-powered low-mass X-ray binary, and appears in one of a number of different X-ray spectral and timing states, which are likely related to the geometry of the accretion flow onto the neutron star (e.g., Hasinger & van der Klis 1989; Shirey et al. 1999). The accretion disk and its accompanying corona are therefore additional X-ray sources that are present during a burst event, and must be taken into account when X-ray burst data is utilized to measure neutron star radii. In a traditional X-ray burst analysis, the X-ray spectrum produced by the disk and corona before the burst (also known as the "persistent" spectrum) is used as the background spectrum to be subtracted from the burst spectrum (e.g., Galloway et al. 2008). The implicit assumption behind this procedure is that the disk and corona are not affected by the burst and their contribution to the spectrum does not change. However, recent analyses of both a large number of bursts in the Rossi X-Ray Timing Explorer (RXTE) archive (Worpel et al. 2013, 2015) and two long individual bursts (in't Zand et al. 2013; Keek et al. 2014) show strong evidence that the normalization of the persistent spectrum may be elevated during the burst by factors of ∼1–10, potentially indicating an increase in the accretion rate, possibly due to Poynting–Robertson (PR) drag (Blumenthal 1974; Walker & Meszaros 1989; Walker 1992). In addition, careful examination of a small number of bursts has shown that the high-energy (>30 keV) persistent emission appears to drop during the burst (Maccarone & Coppi 2003; Chen et al. 2013, 2018; Ji et al. 2013; Kajava et al. 2017), perhaps due to increased Compton cooling in the corona. Taken together, both of these results suggest that the accretion disk and corona may be significantly affected by the strong photon field during an X-ray burst (Degenaar et al. 2018).

The physical response of an accretion disk and corona to the impact of an X-ray burst is not currently well understood. As discussed by Ballantyne & Everett (2005), a wide range of processes and responses are possible, including outflows (driven by radiation pressure), inflows (driven by PR drag), and heating and cooling effects. This situation must be improved if X-ray bursts are to be exploited as effective probes of neutron star physics. Toward this end, this Letter presents the first numerical simulation of an accretion disk subject to the sudden, intense radiation field of an X-ray burst. The results clearly show the impact of the burst on the disk and should provide important guidance for correctly including the effects of the disk when modeling burst spectra. The next section describes the numerical setup, with the results described in Section 3. A short summary and discussion are presented in Section 4.

The simulations presented in this Letter use the Cosmos++ computational astrophysics code (Anninos et al. 2005; Fragile et al. 2012, 2014). The numerical treatments of the hydrodynamics, radiation, and viscosity are nearly identical to those described by Fragile et al. (2018), so we forgo providing extensive details here. The simulations are two-dimensional, axisymmetric, and cover a radial range from r_min = 10.7 km to r_max = 1531 km, with exponential spacing. We use "outflow" boundaries at both the inner and outer radial boundaries, which effectively means that we are assuming any boundary layer effects happen inside r_min. The full range of θ is covered with a latitude coordinate, x₂, related to θ by θ = x₂ + 0.25 sin (2x₂), which concentrates resolution toward the midplane. A layer of static mesh refinement further enhances the resolution in the region of interest. The base resolution is 96², with the additional layer of refinement covering the region r_min ≤ r ≤ 535 km and 0.1π ≤ x₂ ≤ 0.9π (275 ≤θ ≤ 1525) for an effective peak resolution of 192² over the disk.

Both the X-ray burst properties (Galloway & Keek 2017) and the strength of the burst-disk interaction (Ji et al. 2014) are observed to depend on whether the system is in the low/hard or high/soft state when the burst is ignited. These two states are thought to correspond to different accretion regimes and disk geometries. In the low/hard state, the accretion rate is relatively low and the thin accretion disk is possibly truncated at small radii into a geometrically thick, hot flow that acts as a large corona (e.g., Done et al. 2007). In contrast, the high/soft state is related to more rapid accretion through an untruncated, geometrically thin disk extending down to the stellar surface. In this Letter, a geometrically thick disk is assumed for the simulation, which would correspond to the system being in the low/hard state. The results for a thin disk in the high/soft state will be presented in subsequent work.

In constructing the thick disk, we have assumed the spacetime around our 1.45 M_⊙ neutron star can be adequately described by the Schwarzschild metric (in Kerr-Schild polar coordinates), appropriate for non- or slowly rotating neutron stars (e.g., Hartle 1967). We then construct the thick disk as an initially stationary hydrodynamic torus, of the type first described by Kozlowski et al. (1978). We choose a constant specific angular momentum profile, with the inner edge of the torus initially at r_in = 47 km and its center at r_cen = 86 km. This fully specifies the torus solution, other than an arbitrary density normalization, which we adjust to give us the desired mass accretion rate. Accretion is driven by a Shakura & Sunyaev (1973) type α-viscosity with the viscosity coefficient calculated following standard disk theory as μ = νρ = αρc_sH, where α = 0.025, ${c}_{{\rm{s}}}=\sqrt{{P}_{\mathrm{tot}}/\rho }$ is the thermal sound speed, P_tot = P_gas + P_rad is the total pressure, $H={c}_{{\rm{s}}}/{V}^{\phi }$ is the disk height, and V^ϕ is the azimuthal coordinate velocity of the gas. This viscosity rapidly (on the dynamical timescale) redistributes the angular momentum of the initial torus into a nearly Keplerian profile.

Although this torus setup is fairly standard in accretion simulation literature, it does present two related disadvantages for our current purposes. (1) Because the torus is isolated, mass is not replenished as it leaves the simulation domain, either as a result of accretion through the inner boundary (i.e., onto the neutron star) or expulsion beyond the outer boundary (i.e., as part of a wind). (2) This also means that the mass accretion rate through the inner boundary is never truly steady, instead rising sharply at early times and then gradually decaying as the simulation proceeds. For our current setup, the mass accretion rate peaks around $\dot{m}=\dot{M}{c}^{2}/{L}_{\mathrm{Edd}}={10}^{-3}$, where ${L}_{\mathrm{Edd}}=1.7\,\times {10}^{38}$ erg s⁻¹, but decays to below 10⁻⁵ by the end of the simulation. This leads to a net loss of mass in the domain of >97% over the first 2.25 s of the simulation. This dramatic loss of mass is one reason why these simulations were not run for the entire duration of the burst.

X-ray bursts are typically not observed at such low accretion rates (Galloway & Keek 2017), but numerically simulated disks at more realistic mass accretion rates tend to spontaneously become thin (see, e.g., model B in Ohsuga & Mineshige 2011), in contrast with what appears to be the case in nature. This is an issue that has yet to be resolved. As our interest in this Letter is to see what happens to a disk in the low/hard state during an X-ray burst, we make the choice to use an accretion rate that we know will produce a geometrically thick (radiatively inefficient) disk. We suspect that our results would still hold at higher $\dot{m}$ (provided the disk is still in the low/hard state), but this will need to be checked.

Most Type I X-ray bursts exhibit a fast rise and power law or exponential decay (Lewin et al. 1993), so we model the burst luminosity, L(t), using the so-called Norris model from the Gamma-Ray Burst community (Norris et al. 2005):

$$\begin{eqnarray}&&L(t)={L}_{0}{e}^{2{({\tau }_{1}/{\tau }_{2})}^{1/2}}{e}^{\tfrac{-{\tau }_{1}}{t-{t}_{s}}-\tfrac{t-{t}_{s}}{{\tau }_{2}}},\end{eqnarray} \tag{ 1 }$$

where L₀ is the peak burst luminosity, t_s is the burst start time, and τ₁ and τ₂ characterize the burst rise and decay. In this Letter we choose ${L}_{0}={10}^{38}$ erg s⁻¹, t_s = −0.4 s, τ₁ = 6 s, and τ₂ = 1 s, which produces a burst that peaks at a simulation time of t = 2.05 s, lasts for about 10 s, and has a total energy output of 3 × 10³⁸ erg. The simulation duration (2.25 s) and burst peak (2.05 s) times are much longer than the dynamical and thermal timescales in the inner parts of the disk (t_dyn = 2π/Ω ∼ 0.001 s and t_th = 2π/(αΩ) ∼ 0.05 s, respectively), and comparable to the viscous timescale (${t}_{\mathrm{vis}}={r}^{2}/(\alpha {c}_{s}H)\sim 10$ s).

The burst luminosity is introduced into the simulation domain as an outward pointing radiation flux through the inner radial boundary, ${R}_{t}^{r}(t,{r}_{\min })=-L(t)/(4\pi {r}_{\min }^{2})$. Any contribution to this luminosity coming from accretion onto the neutron star surface is neglected, as the accretion luminosity from this disk is quite low (<10³³ erg s⁻¹). We also ignore the angular velocity of the emitting material, consistent with our use of the Schwarzschild metric. Prograde (in the same direction as the orbital motion of the disk) rotation of the neutron star would cause PR drag to be less effective. For the fastest rotating neutron stars, with spin periods of a few ms, the observable effects of PR drag are expected to virtually disappear (Walker 1992). On the other hand, while the M₁ closure scheme that we use for the radiation correctly captures PR drag in the limit of a distant point source (e.g., Wielgus et al. 2016), for a finite-size source such as a neutron star, because M₁ does not properly handle light coming from extended sources, it underestimates the drag (see Figure 1 of Sok Oh et al. 2011). It remains for more accurate radiation treatments to fully assess the quantitative impact of PR drag during Type I X-ray bursts.

As we will show, the dominant effect of the burst radiation comes from Compton cooling of the hot, thick disk. The strength of the Compton cooling goes as $4\rho {\kappa }^{{\rm{s}}}({T}_{\mathrm{gas}}-{T}_{\mathrm{rad}})/{m}_{e}$, where ρ is the gas density, κ^s = 0.34 cm² g⁻¹ is the electron-scattering opacity, m_e is the mass of the electron, and T_gas and T_rad are the gas and radiation temperatures, respectively. Typical values are ${T}_{\mathrm{gas}}\sim {10}^{10}$ K (for disk plasma prior to the burst) and ${T}_{\mathrm{rad}}\approx 4\times {10}^{6}\,{\rm{K}}=0.3\,\mathrm{keV}$ (near the inner radial boundary at the peak of the burst).

To establish the importance of the X-ray burst on the evolution of our accretion disk, two otherwise identical simulations were run, one with the burst activated, the other without. Figure 1 shows the dramatic contrast in density structure and temperature between the two simulated accretion disks when averaged over the final 0.5 s (close to the peak of the burst). The disk without the burst (left) is hot (T ≳ 10¹⁰ K)⁴ and geometrically thick (H/R ≈ 0.3), as expected based on the initialization of the simulation. The case with the burst (right) is almost diametrically opposite: cold (T ≲ 10⁷ K) and geometrically thin (H/R ≲ 0.01).⁵ Clearly, the burst radiation had a conspicuous effect on the structure of the accretion flow.

Zoom In Zoom Out Reset image size

In this low accretion rate, thick-disk limit, we expect the dominant physical effect of the burst to come from inverse Compton cooling of the disk plasma by the burst photons. Degenaar et al. (2018) showed that this should begin to happen once the burst luminosity becomes roughly equal to the accretion luminosity, which is almost immediately in our simulations. The added cooling can be enough to dramatically lower the temperature and remove much of the thermal pressure support in the disk, resulting in its vertical collapse, consistent with what is seen in Figure 1. As further proof, we ran another simulation that included a burst, but with the Compton cooling term turned off in our radiation package. The results of that simulation, with a disk temperature and scale height within a factor of two of the no-burst simulation, confirm the principal role of Compton cooling.

Another possible impact of the burst on the disk could come through radiation pressure driving material away. However, we find that at no time is more than 4% of the total mass in our simulation domain caught up in an unbound (u₀ < −1) outflow (v^r > 0), i.e., a wind. Further, that 4% is a transient phase early in the simulation (t < 0.1 s), well before the burst peak; at later times, the mass caught up in an unbound outflow becomes completely negligible. The mass in bound, but outflowing, material is much larger (≳50%), but this is to be expected as the initial torus is viscously spreading (about half of the material moving outward and half moving inward). Thus, we conclude that the differences that we see in Figure 1 cannot be attributed to radiation pressure driving, which is consistent with the sub-Eddington nature of our burst.

A third predicted effect of Type I X-ray burst photons on the surrounding accretion disk is to drive additional accretion through PR drag. The red dashed line in Figure 2 shows that the X-ray burst does enhance the mass accretion rate through the inner simulation boundary by factors of ≈3–4, but it is obvious that the enhancement to $\dot{m}$ does not happen on the same timescale as the burst (which peaks at 2.05 s). Before addressing this discrepancy, let us first look for direct evidence of PR drag in the simulations. Figure 3 shows a spacetime diagram of the ratio between the density-weighted orbital angular frequency of gas in the disk for the simulations with and without a burst. In the case without a burst, the orbital angular frequency maintains a steady Keplerian profile throughout the simulation (not shown). However, this ratio plot shows that the simulation with a burst has orbital angular frequencies that begin to decay (orbit slower) on exactly the timescale of the burst, as would be expected if the disk material were being acted upon by an additional torque.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, this PR drag cannot be the only cause for the enhanced accretion seen in our burst + Compton simulation (Figure 2; red dashed curve). First, the enhanced accretion starts almost immediately and peaks well before Figure 3 shows significant deviations in the angular momentum profile (t ≳ 0.75 s). Instead, this early accretion appears to be tied to the Compton-driven cooling and collapse of the disk, which happens on the local thermal timescale. As evidence of this link, note that the simulation with a burst but without Comptonization (Figure 2; green dotted curve) does not show the same early enhancement of accretion. In both cases, though, the accretion rate peaks before the burst does. This has to do with the fact that the inner part of the disk cannot refill (by accreting material from larger radii) as fast as it is being drained by these additional processes. The mass accretion rate onto a neutron star during a burst will, thus, be sensitive to both the strength of the PR drag and the time-varying surface density profile of the disk.

Using general relativistic, viscous, radiation hydrodynamic simulations, we have confirmed that the emission from a Type I X-ray burst can affect a surrounding hot, thick accretion flow in at least two important ways: (1) the burst can Compton cool it by orders of magnitude, resulting in a temporary transition to a much cooler, thinner disk; and (2) the combination of the collapse of the disk and PR drag can amplify the mass accretion rate by a factor of a few during the burst. This is the first time that these effects have been demonstrated in such simulations.

This study focused exclusively on Type I X-ray bursts occurring during the low/hard state. We find that the relative thickness of the disk in this state presents a conveniently large target for the burst photons. The rapid (t ≲ a few t_th) cooling in the simulations is in good agreement with observations of a prompt drop in the high energy (>30 keV) flux during real X-ray bursts in the hard state (e.g., Chen et al. 2018). The three order of magnitude drop in temperature, however, is significantly larger than what is inferred from observations. Part of this discrepancy can be explained by the fact that we did not include the effects of pair production, which would have acted as a thermostat, regulating the maximum temperature of the plasma. Additionally, the unrealistically low disk mass accretion rate that we considered probably contributed to the excessively high temperatures.

Another highlight of our simulations was the confirmation of PR drag as an important process acting on an accretion disk during a burst (Figure 3). However, its impact was relatively modest, as it only enhanced the mass accretion rate onto the neutron star by factors of ≈3–4, less than the observed increases in the persistent flux of ≳4 (Worpel et al. 2015). However, we have to be careful here. Although we have chosen the most ideal scenario for PR drag, in that the applied radiation field is from a non-rotating neutron star, the M₁ closure scheme employed underestimates PR drag from extended sources. Therefore, future simulations will be needed to better quantify the effect of PR drag.

In summary, this Letter presents the first simulation of an accretion disk impacted by a powerful X-ray burst. The parameters were chosen so that the effects of the burst luminosity and PR drag would be maximized. This was sufficient to demonstrate that a Type I X-ray burst going off in the hard state may lead to important structural and thermodynamic changes in an accretion disk. These changes need to be understood before, for example, burst spectroscopy could be used to measure neutron star radii.

P.C.F. and J.W.L.W. acknowledge support from SC NASA EPSCoR RGP 2017 and National Science Foundation grant AST-1616185. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. D.R.B. and T.J.M. thank the International Space Science Institute, where part of this work was carried out, for their hospitality.