Relativistic Tidal Disruption and Nuclear Ignition of White Dwarf Stars by Intermediate-mass Black Holes

Ignoring radiation, magnetic fields, and viscous contributions, none of which are exercised in this work, the stress energy tensor simplifies to

$$\begin{eqnarray}&&{T}^{\alpha \beta }=\rho {{hu}}^{\alpha }{u}^{\beta }+{{Pg}}^{\alpha \beta },\end{eqnarray} \tag{ 1 }$$

where ρ is the fluid mass density, $h=1+\epsilon /{c}^{2}+P/(\rho {c}^{2})$ is the specific enthalpy, c is the speed of light, ${u}^{\alpha }={u}^{0}{V}^{\alpha }$ is the contravariant 4-velocity, V^α is the transport velocity, g_αβ is the curvature metric, and P is the fluid pressure providing closure given an appropriate equation of state (see Section 2.2).

The four fluid equations (energy and three components of momentum) are derived from the conservation of stress energy: ${{\rm{\nabla }}}_{\mu }{T}_{\ \nu }^{\mu }={\partial }_{\mu }{T}_{\ \nu }^{\mu }+{{\rm{\Gamma }}}_{\alpha \mu }^{\mu }{T}_{\ \nu }^{\alpha }-{{\rm{\Gamma }}}_{\mu \nu }^{\alpha }{T}_{\ \alpha }^{\mu }={S}_{\nu }$, where ${{\rm{\Gamma }}}_{\mu \nu }^{\alpha }$ are the Christoffel symbols and S_ν represent arbitrary source terms. In addition to energy and momentum, we also require an equation for the conservation of mass ${{\rm{\nabla }}}_{\mu }(\rho {u}^{\mu })={\partial }_{0}(\sqrt{-g}{u}^{0}\rho )$ + ${\partial }_{i}(\sqrt{-g}{u}^{0}\rho {V}^{i})=0$.

Expanding out space and time coordinates, the four-divergence (${{\rm{\nabla }}}_{\mu }{T}_{\nu }^{\mu }={S}_{\nu }$) of the mixed index stress tensor is written as

$$\begin{eqnarray}&&{\partial }_{0}(\sqrt{-g}{T}_{\nu }^{0})+{\partial }_{i}(\sqrt{-g}{T}_{\nu }^{i})=\sqrt{-g}{T}_{\sigma }^{\mu }\ {{\rm{\Gamma }}}_{\mu \nu }^{\sigma }+\sqrt{-g}{S}_{\nu }.\end{eqnarray} \tag{ 2 }$$

Defining energy and momentum as ${ \mathcal E }=-\sqrt{-g}{T}_{0}^{0}$ and ${{ \mathcal S }}_{j}\,=\sqrt{-g}{T}_{j}^{0}$, the equations further take on a traditional transport formulation

$$\begin{eqnarray}&&{\partial }_{0}{ \mathcal E }+{\partial }_{i}(-\sqrt{-g}{T}_{0}^{i})=-\sqrt{-g}{T}_{\sigma }^{\mu }\ {{\rm{\Gamma }}}_{\mu 0}^{\sigma },\end{eqnarray} \tag{ 3 }$$

or

$$\begin{eqnarray}&&{\partial }_{0}{ \mathcal E }+{\partial }_{i}({ \mathcal E }{V}^{i})+{\partial }_{i}(\sqrt{-g}{{PV}}^{i})=-\sqrt{-g}{T}_{\sigma }^{\mu }\ {{\rm{\Gamma }}}_{\mu 0}^{\sigma },\end{eqnarray} \tag{ 4 }$$

for energy, and

$$\begin{eqnarray}&&{\partial }_{0}{{ \mathcal S }}_{j}+{\partial }_{i}(\sqrt{-g}{T}_{j}^{i})=\sqrt{-g}{T}_{\sigma }^{\mu }\ {{\rm{\Gamma }}}_{\mu j}^{\sigma },\end{eqnarray} \tag{ 5 }$$

or

$$\begin{eqnarray}&&{\partial }_{0}{{ \mathcal S }}_{j}+{\partial }_{i}({{ \mathcal S }}_{j}{V}^{i})+{\partial }_{i}(\sqrt{-g}P\,{g}_{j}^{0}\,{V}^{i})=\sqrt{-g}{T}_{\sigma }^{\mu }\ {{\rm{\Gamma }}}_{\mu j}^{\sigma }.\end{eqnarray} \tag{ 6 }$$

for momentum, and

$$\begin{eqnarray}&&{\partial }_{0}{ \mathcal D }+{\partial }_{i}({ \mathcal D }{V}^{i})=0,\end{eqnarray} \tag{ 7 }$$

for mass conservation, where ${ \mathcal D }=\sqrt{-g}{u}^{0}\rho =W\rho$ is the boost density.

Mesh motion is implemented by a straightforward replacement of generic advective terms:

$$\begin{eqnarray}&&{\partial }_{0}(\sqrt{-g}{T}_{\alpha }^{0})+{\partial }_{i}(\sqrt{-g}{T}_{\alpha }^{0}{V}^{i})\end{eqnarray} \tag{ 8 }$$

with

$$\begin{eqnarray}&&{\partial }_{0}(\sqrt{-g}{T}_{\alpha }^{0})+{\partial }_{i}(\sqrt{-g}{T}_{\alpha }^{0}({V}^{i}-{V}_{g}^{i}))+\sqrt{-g}{T}_{\alpha }^{0}{\partial }_{i}{V}_{g}^{i},\end{eqnarray} \tag{ 9 }$$

where ${V}_{g}^{i}$ is the grid velocity, and ${T}_{\alpha }^{0}$ is used here to represent any evolved field, including ${ \mathcal E }$, ${{ \mathcal S }}_{j}$, and ${ \mathcal D }$. The grid velocity can be set arbitrarily by the user, or computed internally by Lagrangian, softened Lagrangian, or potential-relaxing algorithms. At the end of each cycle, the zone objects (consisting of node positions, oriented face area vectors, cell volumes, etc.) are physically moved by projecting the zone-centered grid velocity to the nodes and moving them accordingly to the next time cycle, after which all local zone object elements (including the spacetime metric) are updated.

At the beginning of each time cycle a series of coupled nonlinear equations are solved to extract primitive fields (mass density, internal energy, velocity) from evolved conserved fields (boost density, total energy, momentum), after which the equation of state is applied to compute the thermodynamic quantities: pressure, sound speed, and temperature. For Newtonian systems this procedure is straightforward, but relativity introduces a nonlinear interdependence of primitives, so their extraction from conserved fields requires special (Newton–Raphson) iterative treatment. We have implemented several procedures for doing this, solving one-, two-, or five-dimensional inversion schemes for magnetohydrodynamics, or a nine-dimensional fully implicit method (including coupling terms) when radiation fields are present (Fragile et al. 2012, 2014).

We adopt a two-component equation of state for this work:

$$\begin{eqnarray}&&P={P}_{\mathrm{rad}}+{P}_{\mathrm{gas}}={a}_{R}{T}^{4}/3+({\rm{\Gamma }}-1)Y\rho T,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&e-{e}_{{\rm{cold}}}={e}_{{\rm{rad}}}+{e}_{{\rm{gas}}}={a}_{R}{T}^{4}+Y\rho T,\end{eqnarray} \tag{ 11 }$$

where thermal radiation and ions both contribute to the energy density ($e\equiv \rho \epsilon$) and pressure (P). The temperature is calculated from the thermal component of the internal energy density, subtracting off the "cold" energy approximated by the initial barotropic equation of state that specifies the initial WD data in the absence of shocks: ${e}_{\mathrm{cold}}=K{\rho }^{{\rm{\Gamma }}}/({\rm{\Gamma }}-1)$. Additionally, Y is defined as

$$\begin{eqnarray}&&Y=\displaystyle \frac{{k}_{b}}{{m}_{n}}\left(\displaystyle \frac{1}{{\mu }_{W}}+{N}_{e}\right)\displaystyle \frac{1}{{\rm{\Gamma }}-1},\end{eqnarray} \tag{ 12 }$$

where k_b is Boltzmann's constant, m_n is the nucleon mass, ${\mu }_{W}=1/({\sum }_{i}{\rho }_{i}/(\rho \,{A}_{i}))$ is the mean molecular weight formed from isotopic densities, A_i is the atomic weight of isotope i, N_e is the number of electrons per nucleon, and Γ = 5/3 is the adiabatic index. The extraction of temperature from Equation (11) follows Haas et al. (2012) (see also Lee et al. 2005), but we have generalized (10) to account for radiation pressure. Haas et al. (2012) demonstrated (as we have here) this simplified multi-component treatment capable of reproducing qualitatively similar phase tracks as more complex degenerate equation-of-state models (RRH09) for the range and scope of stellar and interaction parameters considered in this work. This approximation might adequately capture hydrodynamic behavior, but it is unclear the extent to which it affects calculations of nuclear burn products, which can be highly sensitive to temperature. We anticipate in related future work to implement a more accurate Helmholtz model and assess sensitivities to choices of equation of state.

The nonlinear temperature dependence in Equation (11) requires Newton inversion each cycle to solve for the temperature given the fluid density and internal energy. Partial ionization is approximated by interpolating the number of electrons per nucleon N_e in Equation (12), between zero (no ionization) and 1/2 (full ionization), over the range of temperatures T = 5 × 10⁵ K to T = 5 × 10⁶ K. The mean molecular weight μ_W is updated each cycle from the isotopic fractions produced after solving the nuclear reaction network discussed in the following section.

Cosmos++ supports a number of nuclear reaction networks, including the 7- and 19-isotope α-chain and heavy-ion reaction models developed by Weaver et al. (1978), Timmes (1999), and Timmes et al. (2000). There is, of course, a trade-off between speed and accuracy as measured by nuclear energy production when choosing between these two networks. Timmes et al. (2000) report 30%–40% differences in energy production rates from the two models across different stellar environments and burn stages. We find that the more accurate 19-isotope model imposes a tolerable increase in computing time, so we adopt that network for this work. The 19-isotope network solves heavy-ion reactions using a standard α-chain composed of 13 nuclei (⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, and ⁵⁶Ni), plus additional isotopes (¹H, ³He, ¹⁴N, ⁵⁴Fe) to accommodate some types of hydrogen burning as well as neutrons and protons from photodisintegration.

To account for the advection of reactant species, the hydrodynamic equations described in Section 2.1 are supplemented with an additional set of mass continuity equations, one for each isotope j in the network, resulting in the following 19 partial differential equations that need to be solved each advance cycle:

$$\begin{eqnarray}&&{\partial }_{0}(W{\rho }_{j})+{\partial }_{i}(W{\rho }_{j}{V}^{i})=\displaystyle \sum _{k,l,m}{\rho }_{m}{\rho }_{l}{R}_{{lk}}(m)-{\rho }_{j}{\rho }_{k}{R}_{{kl}}(j).\end{eqnarray} \tag{ 13 }$$

Equation (13) is written in simplified form explicitly for, but not limited to, binary reactions, where R_lk(m) represents creation rates contributed by species m, and R_kl(j) are the corresponding destruction rates for species j. These rate coefficients are highly nonlinear in temperature and span many orders of magnitude. As a result, the differential equations for the reaction network (the ordinary differential equations consisting of the right side of Equation (13) after operator splitting off advection) are stiff and must be solved with implicit or semi-implicit integration methods. We have implemented a number of solvers and evaluated their performance in some of our calculations of tidal events, including fourth-, fifth-, and eighth-order explicit and semi-implicit adaptive time-step Runge–Kutta schemes, the variable-order Bader–Deufelhard method, and a fifth-order fully implicit method. The fully implicit solver proved faster and more robust for this class of problems, completing the 3D calculations faster (by a factor of several) and with fewer (by a factor of three) first-attempt network convergence failures than its closest competitor (Bader–Deufelhard).

Network failures occur when the solver does not converge to a user-specified fractional tolerance for each species abundance, which we set to 10⁻⁶. When the solver fails to meet this tolerance over a full hydrodynamic cycle, that particular zone is flagged and switched over to a subcycling method where the time step is systematically reduced until the solver converges. The equation of state is recalculated consistently at each subcycle, as is the internal energy production, properly updating the temperature and reaction rates before solving each new subcycled reaction. The local minimum time step needed for proper convergence is recorded and used to set the global time step for the next hydrodynamic cycle. The global time step is also constrained to limit relative energy production to 25% fractional change. These highly restrictive conditions are necessary for achieving stable and accurate solutions since nuclear reaction time steps can easily become much smaller than hydrodynamic timescales.

The specific energy released by nuclear reactions is calculated from

$$\begin{eqnarray}&&\delta {\epsilon }_{\mathrm{nuc}}={N}_{A}\displaystyle \sum _{j}{B}_{j}\delta {Y}_{j},\end{eqnarray} \tag{ 14 }$$

where N_A is Avogadro's number, B_j is the binding energy, and Y_j is the dimensionless molar abundance of isotope j. It is defined in terms of the mass fractions X_j = ρ_j/ρ, atomic weights A_j (proton + neutrons), total mass density ρ, and number densities n_j as ${Y}_{j}={X}_{j}/{A}_{j}={\rho }_{j}/(\rho {A}_{j})={n}_{j}(\rho {N}_{A})$. The corresponding change in internal energy (accumulated over all subcycled intervals) δe = ρδ_nuc is then used to update the relativistic momentum and total energy fields, after also updating the change in pressure δP from the equation of state:

$$\begin{eqnarray}&&{{ \mathcal E }}^{n+1}={{ \mathcal E }}^{n}-((\delta e+\delta P){{Wu}}_{0}+\sqrt{-g}\,\delta P),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{ \mathcal S }}_{j}^{n+1}={{ \mathcal S }}_{j}^{n}+(\delta e+\delta P){{Wu}}_{j}.\end{eqnarray} \tag{ 16 }$$

In this formulation, we have assumed that the WD material is sufficiently optically thick (to photons and neutrons) that released energy cannot be radiated away and is instead absorbed instantaneously in cells where it has been emitted.

The performance of the reaction solver and nuclear energy production model have been validated against the entire suite of problems and code tests proposed by Timmes et al. (2000), in addition to adiabatically coupled nuclear flows (e.g., standard Big Bang nucleosynthesis). The hydrodynamic coupling and subcycling procedure are identical to the methods we had previously developed and vetted for chemical reactive networks with radiative cooling, and we have subjected them to numerous tests, including self-similar reactive flows, molecular species formation in cosmological shocks, and hydrodynamic shock masking. Here we have only substituted nuclear rates for chemical ones.

All simulations use a 3D Cartesian grid centered on the WD, covering a spatial domain of 2R_WD along the z-axis and 4R_WD in the x–y orbital plane, where R_WD is the radius of the WD. For most of our calculations, the domain is resolved with three hierarchical grid layers with 96³ cells on the base level and two additional refinement levels (each level doubling the resolution) for an effective 384³ resolution covering the star. For a WD with radius 10⁹ cm, this results in an initial uniform cell resolution (Δx, Δy, Δz) of ∼(2.1, 2.1, 1.1) × 10⁷ cm, or in terms of scaled dimensions ∼(R_S,3/15, R_S,3/15, R_S,3/30), where R_S,3 is the Schwarzschild radius for a 10³ M_⊙ black hole. Adaptive mesh refinement (AMR) is used throughout the calculations, refining and de-refining every five hydrodynamic cycles. AMR is triggered by the local mass density, refining where the density exceeds ρ_max/20 and de-refining where the density drops below ρ_max/200, where ρ_max is the initial maximum cell density across the grid.

Since the TD events considered here have much in common with our earlier work on the disruption of the G2 cloud orbiting the center of our Galaxy (Anninos et al. 2012), we adopt a similar approach for implementing a grid velocity to track the trajectory of the WD. Grid velocity is used to move the base grid at the mean (density-weighted) Lagrangian velocity of the star, following its orbital motion in the x–y orbital plane while keeping the star centered on the grid. AMR is performed on the moving base mesh. The combination of grid velocity and AMR together allows us to use fewer refinement levels to achieve high resolution across the stars throughout their evolution. By comparison, HSBL12 built eight levels of refinement to achieve comparable resolution. We note, however, that we additionally improve on the vertical scale height resolution as described next.

Because the WD is compressed along the z-axis, perpendicular to the orbital plane, mesh motion along that direction requires special treatment. Tanikawa et al. (2017) argued the importance of capturing shocks within the scale height of the WD to resolve nuclear explosions. Estimating the scale height z_min of the WD in the z-direction as z_min/R_WD ∼ β^{−2/(Γ−1)} at periapsis, where β is the tidal strength parameter defined in the next section (Luminet & Carter 1986; Brassart & Luminet 2008; Tanikawa et al. 2017), we anticipate that a WD with an initial radius of 10⁹ cm will require a cell resolution of about 10⁶ cm in a tidal field with β = 10. We therefore augment the orbital mesh velocity by introducing a component converging toward the z = 0 plane with magnitude

$$\begin{eqnarray}&&| {V}_{g}^{z}(t)| ={\zeta }_{2}\,| {V}_{L}^{z}(t)| \,{\left|\displaystyle \frac{z(t)}{{L}_{z}(t)}\right|}^{{\zeta }_{N}},\end{eqnarray} \tag{ 17 }$$

where z(t) is the time-dependent coordinate, L_z(t) is the changing vertical length of the grid along the positive z-axis (distance from the orbital plane to the outer boundary), ζ₂ ≥ 0 is a constant multiplier, and $| {V}_{L}^{z}(t)|$ is the average magnitude of the downward-directed z-component of the star's internal Lagrangian velocity. The factor $| z(t)/{L}_{z}(t){| }^{{\zeta }_{N}}$ smoothly scales back the vertical grid motion to zero at the midplane, preventing mesh singularities from forming. Large values of the exponent ζ_N > 1 tend to collapse the outer grid edges faster than the midplane cells. Fractional values 0 < ζ_N < 1 tend to collapse the midplane zones at a faster rate. We choose ζ_N = 0.85 as a reasonable compromise. The multiplier ζ₂ is initially set to 5 but deactivated (reset to zero) when the vertical scale height resolution improves tenfold to ∼10⁶ cm, or equivalently R_S,3/300. This prescription allows the innermost and outermost cells to move inward toward the orbital plane at different velocities, collapsing the inner zones by a factor of 10, while cutting the distance to the outer z-boundary planes by about two-thirds, depending on specific case parameters.

When the WD passes periapsis, grid motion in the x–y plane is modified to keep both the star and black hole on the grid for the remaining evolution to capture flow circulation. This is accomplished by continuing to track the orbital velocity of the star, while also allowing the grid to stretch according to

$$\begin{eqnarray}&&{V}_{g}^{x}=\min (0,{v}_{L}^{x}(t))\displaystyle \frac{{x}_{\max }(t)-x(t)}{{x}_{\max }(t)-{x}_{\min }(t)}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{V}_{g}^{y}=\min (0,{v}_{L}^{y}(t))\displaystyle \frac{{y}_{\max }(t)-y(t)}{{y}_{\max }(t)-{y}_{\min }(t)},\end{eqnarray} \tag{ 19 }$$

where x_max/min(t) and y_max/min(t) are the evolving maximum/minimum x and y grid positions. This choice of grid velocity holds the upper right edge of the grid fixed to retain the black hole on the grid, while allowing the lower left edge to stretch with the tidal debris and approximately track the stellar centroid as it flows past and around the black hole.

The boundary conditions on the outer edges of the grid are set to outflow conditions. All ghost-zone quantities are equated to the values of their adjacent internal-zone neighbors, except that the velocity component normal to the boundary is set to zero if it points onto the grid. Note that this does not necessarily prevent material from coming onto the grid because the parabolic reconstruction of the fluxes may still produce negative or positive values.

Several deformation modes are active during typical disruption encounters. In the x–y orbital plane of the binary system, the star will experience tidal stretching along the direction connecting the black hole to the WD, and tidal compression perpendicular to this direction. The star will also undergo tidal compression in the direction perpendicular to the orbital plane. Compression perpendicular to the orbital plane is significantly greater than that in the orbital plane, and is the primary mechanism that could, in principle, ignite nuclear reactions in the WD. These deformation modes are demonstrated in Figures 3 and 4, which image the logarithm of mass density and temperature in the x–y (orbital) and y–z planes for run B3M2R06. Stretching and compression in the orbital plane are both evident in the series of images in Figure 3, conspiring to simultaneously thin and elongate the WD in orthogonal directions whose orientations change as the WD approaches the BH. Compression acting perpendicular to the orbital plane is on display in Figure 4, where the WD collapses to a fraction of its original size. Both figures demonstrate the effectiveness of the prescription for mesh velocity described in Section 2.4 to track the orbital trajectory and at the same time collapse the grid in response to tidal compression along the vertical direction. Figures 3 and 4 are representative of all encounter scenarios we have studied in that nuclear burning is not triggered and thus has little effect on stellar disruption until the approach and passage of periapsis. Ignition is triggered mainly along the dense filamentary-like nozzle structure that forms where the stellar material is pinched to maximal compression as it passes periapsis, and is clearly evident in the third plates of Figure 3 projecting radially outward from the black hole. Apart from the creation of burn products, the important hydrodynamical effect of nuclear ignition is the deposition of released energy that contributes to the overall disruption of the WD.

Zoom In Zoom Out Reset image size

We attribute the small-scale structures present in the third panel of Figure 3 to transitory hydrodynamic instabilities. Material compressing onto the orbital plane before reaching the nozzle is stable, but as Figure 11 shows, this material bounces after passing through the nozzle and accelerates outward from the orbital plane, reversing the direction of its vertical velocity and reshocking as it hits infalling material. In addition, the infalling material is not perfectly aligned with the dense thin shell of material forming along the nozzle, which creates tangentially aligned streamlines within the dense layer that can potentially induce perturbations (though we cannot conclusively say that this actually happens due to the lack of sufficient resolution in the orbital plane). In all, this creates an interpenetrating, multi-shock morphology inside the disrupted star that gives rise to interspersed Rayleigh–Taylor unstable regions and strong shear flows that also trigger Kelvin–Helmholtz instabilities. The result is a complex flow pattern downstream of the nozzle characterized by debris separation into filamentary structures. We have verified that these features are robust (present) regardless of numerical treatments, including activation of AMR, mesh motion, thermonuclear reactions, or equation of state (they are present even with a simple ideal gas law with no ionization temperature thresholds). These features, however, are strongly sensitive to grid resolution and are less pronounced when the central plane is inadequately resolved. In fact, when we increase resolution, the structures become more coherent, numerous, and smaller in scale, and we expect this trend to continue at resolutions beyond our current limits. Upon close inspection these features are evident also in the work of other authors (see for example Kawana et al. 2018), but they are more pronounced in our calculations due to increased spatial resolution and perhaps our use of less dissipative, high-order shock-capturing numerical methods. Although we can say with some degree of confidence that these structures are to be expected, clearly more careful studies with much higher resolution will be needed to better understand their cause and effects.

In general, and as expected, we find that more massive WDs compress to greater densities and temperatures, as do WDs on trajectories approaching closer to the black hole. The density–temperature phase tracks calculated for the densest 10% of stellar material and plotted in Figure 5 demonstrate this for a few cases involving encounters with a 10³ M_⊙ black hole. Notice the positive correlation of increasing compression and heating with decreasing perihelion radius. Also shown in this figure is the curve separating dynamical and helium burn reaction timescales. If the dynamical timescale ${\tau }_{D}=1/\sqrt{G{\bar{\rho }}_{\mathrm{WD}}}$ is shorter than the helium combustion timescale ${\tau }_{B}\approx 9\times {10}^{-4}{({T}_{9})}^{3}{({\rho }_{6})}^{-2}\exp (4.4/{T}_{9})$ s (Khokhlov & Ergma 1986), where ${\bar{\rho }}_{\mathrm{WD}}$ is the mean WD density, T₉ = T/10⁹ K, and ρ₆ = ρ/10⁶ g cm⁻³, the star can expand rapidly enough to quench burning by cooling the WD interior and reducing its density. Most of the cases we have studied track into the regime with a fast combustion timescale, achieving density compression ratios of more than an order of magnitude and temperatures exceeding those needed to burn carbon–oxygen chains. The maximum densities and (separately) temperatures observed in our parameter studies are summarized in Table 2. The resolution case studies of B3M6R20, B3M6R12, and B3M6R06 collectively suggest that the peak densities and temperatures converge to about 20%. This might be adequate if the only consideration were hydrodynamic behavior or global nuclear energy production (see Section 4.2). However, nuclear reaction rates are highly sensitive to temperature and we find somewhat greater uncertainty in the isotopic distributions, as discussed in more detail below.

Zoom In Zoom Out Reset image size

Mass accretion rates are plotted as a function of time in Figure 6 for a sample of cases. There are two distinct signature regions found in all scenarios except those with the most distant interaction (largest R_P): a first burst that peaks for the closest encounters between 10⁶ and 10⁷ M_⊙ yr⁻¹ and lasts for 0.5–1 s centered at periapsis (t ∼ 0). This initial burst is followed by a quasistatic phase with magnitudes in the range 10²–10⁴ M_⊙ yr⁻¹ and a time dependence that, at least for the closest encounters, scales roughly as $\dot{M}\propto {10}^{-0.9t}$ as the flow begins to circulate around the black hole one to two seconds after periapsis. The range of accretion rates found in the quasistatic phase is considerably greater than the rate (∼50 M_⊙ yr⁻¹) associated with background gas accretion represented by the horizontal line in Figure 6. The background accretion rate is highly sensitive to the background density set in the calculations, but the quasistatic rate is not and appears to be affected only by grid resolution, which dictates coverage of the BH surface boundary, and numerical dissipation. Comparing against purely hydrodynamic calculations, we find that nuclear burning enhances tracks of mass accretion rates, total accreted mass, and peak density and temperature by modest amounts, roughly 5%–10%.

Zoom In Zoom Out Reset image size

These accretion rates are many orders of magnitude greater than the rates expected from Eddington-limited accretion, which, assuming a standard 10% radiative efficiency, implies a rate

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Edd}}=\displaystyle \frac{{L}_{\mathrm{Edd}}}{0.1{c}^{2}}=4\times {10}^{-5}\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{10}^{3}\,{M}_{\odot }}\right)\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 24 }$$

for hydrogen-poor gas with electron scattering opacity ${\kappa }_{\mathrm{es}}\,=0.2\,(1+{X}_{H})\to 0.2$ cm² g⁻¹, where ${L}_{\mathrm{Edd}}=4\pi {{cGM}}_{\mathrm{BH}}/{\kappa }_{\mathrm{es}}$. Assuming a late-time fall-back scaling of t^−5/3 (Rees 1988) applied at a return time of ≈10² s, we can reasonably expect these disruption events to emit X-ray and γ-ray transients and shine at Eddington luminosity for a year or more.

One of the objectives of this work is to extend the results of RRH09 to the relativistic regime, and provide independent validation concerning the nucleosynthetics of TDEs using alternative high-resolution grid-based numerical techniques and expanded nuclear reaction networks. RRH09 used smoothed particle hydrodynamics methods coupled with a seven-isotope α-network to explore a multi-dimensional parameter space of TDEs covering a range of black hole masses, WD masses, and tidal strengths in search of conditions under which explosive nuclear burning can occur. They found a generically universal dependence on tidal strength and to a lesser degree on the WD mass. In particular they found that ignition occurs when β ≥ 3 (β ≥ 5) for WDs of 1.2 (0.2 and 0.6) solar masses. Their results also suggest a nominal dependence on the central potential model, generally finding stronger reactions (greater iron-group production) with a pseudopotential that accounts for some near-field relativistic effects compared to purely Newtonian treatments. Hence our interest in exploring fully general relativistic interactions.

Additionally, recent studies by Tanikawa et al. (2017) have cast some doubt on the reliability (or limitation) of multi-dimensional numerical calculations due to the difficulty of allocating sufficient resources in the central ignition region. They argue that coarse resolution does not adequately resolve shock heating and can produce spurious heating and artificially trigger ignition. At tidal strengths β ∼ 5, Tanikawa et al. (2017) suggest that a spatial resolution of about 10⁶ cm is required to resolve the scale height sufficiently in order to produce reliable or convergent nuclear reactants. We have adapted our hybrid AMR/moving mesh technique to achieve slightly better than the recommended resolution (≲10⁶ cm) in the central orbital plane, improving on previous work by about an order of magnitude.

A number of diagnostics can be used to quantify the activation or efficiency of nuclear ignition, including production of calcium- and iron-group elements and the total amount of energy released from thermonuclear reactions. We discuss element production in subsequent paragraphs, but focus first on energy release. The total energy produced by nuclear reactions is presented in Table 2 for all the cases we considered, and in Figure 7 from a select sample. This quantity is not strongly sensitive to spatial resolution, and is thus a fairly robust measure of thermonuclear ignition. The trend for energy production is generally monotonically increasing with decreasing R_P (or increasing β) as the WD experiences greater compression, so long as significant chunks are not captured by the BH to limit the amount of nuclear fuel available for production. This effect is observed when comparing B3M2R09 and B3M2R06. The latter produces slightly less nuclear energy and iron elements despite achieving higher peak density and temperature.

Zoom In Zoom Out Reset image size

We find that all of our calculations (β ≥ 2.6) ignite, but only those cases with β ≳ 4.5 release nuclear energies comparable to the star's binding energies (1.3 × 10⁴⁹ and 1.3 × 10⁵⁰ erg for the 0.2 M_⊙ and 0.6 M_⊙ WDs, respectively). These results are generally consistent with Rosswog et al. (2009), as are in select (or common) cases the magnitudes of released nuclear energy. For example, comparing case B3M2R20 to run 6 from RRH09 (with identical tidal strengths β = 5), we find that e_nuc agrees to within a factor of two. However, as we point out below, this level of agreement is not found in all diagnostics.

An exploding star that releases more energy than its binding energy may not necessarily produce copious amounts of iron-group elements. Hence we look separately at the mass of iron-group products as a complementary measure of ignition efficiency. For the 19-isotope network used in our calculations, iron-group is defined as the sum of iron and nickel isotopes. The peak abundance of iron summed over the entire grid in each simulation is shown in Table 2, and plotted as a function of time in Figure 8 for a couple representative cases along with the other isotopic species. We additionally introduce silicon-group (sum of Si, S, and Ar) and calcium-group (sum of Ca, Ti, and Cr) elements to downsample and simplify the number of species curves in Figure 8. A comparison of the particular case B3M2R09 (left panel) can be made to Figure 8 in RRH09 by further combining our silicon- and calcium-group curves. Although not identical, run 2 in RRH09 and B3M2R09 set the same WD and BH masses and comparable β (11 versus 12). We observe a similar mass deficit for helium elements, an iron abundance that agrees to within a factor of a few, but very different distributions of all the other minor isotopes. These differences might be explained by the adoption (or neglect) of trace elements in the initial data or simplifications necessitated by the different α-chain networks (neither the 7- nor 19-isotope network can be expected to accurately predict precise distributions, because each makes simplifications to the reality that would otherwise require modeling hundreds or perhaps thousands of isotopes and reaction paths). The right panel in Figure 8 plots the species distributions for a moderately disrupted massive (0.6 M_⊙) WD, case B3M6R12. This is an interesting example where the nuclear flows are choked off before iron production can be completed. In this case, both calcium- and silicon-group elements dominate nickel and iron.

Zoom In Zoom Out Reset image size

Figure 9 plots a time sequence of three isosurface contours showing the distribution of helium, Ca-group, and Fe-group elements from the same case (B3M2R06) and times as Figures 3 and 4. Ca-group elements are first formed when stellar matter is compressed to the proper density and temperature at and around periapsis, and preferentially along the inner surface closest to the black hole. Iron-group elements are then subsequently formed in their (downstream) wake. Comparing the third time images in Figures 3 and 9 shows that nucleosynthesis is initiated where the star is maximally compressed along the dense radial filament cut across the star in the orbital plane as it passes periapsis. Nuclear products are then subsequently transported in a circular fashion and mixed alongside unignited stellar and background material around the black hole.

Figures 10 and 11 present a different orthogonal view of the B3M6R06 event that highlights the complexity of reactive flows and demonstrates the need for high grid resolution. The left and right panels in Figure 10 are taken at the same time (t = 0) and represent the same quantities, but are calculated at different grid resolutions: the left panel at L2, the right at L3. In both images, red represents the mass fraction of combined carbon and oxygen (the initial composition for this event), and blue represents the combined mass fractions of nuclear products (Si+Ca+Fe groups). Figure 11 shows the mass density (in color) and the velocity flow (aligned with the arrows) corresponding to the high-resolution panel. In this orientation, the black hole is behind the (zoomed in) stellar matter, which is moving from right to left around the black hole, as it pinches to maximum compression then bounces off the orbital plane. The dense radial filament observed in Figure 3 (though for a different case, B3M2R06) is orthogonal to the image plane and projects toward the black hole at the point of maximum compression located approximately in the middle of these images (vertically and horizontally). The vertical scale of the image is zoomed in to roughly 6% of the initial star domain. The physical length scale between tick marks is the same on both vertical and horizontal axes. These images complement Figure 9, providing a more complete picture of the disruption process, and begin to demonstrate the complex structure of the WD core when it ignites. As noted earlier, nuclear products first form at and along the pinched feature where density and temperature are highest, responding in time to the evolving shock structure off the orbital plane and producing interpenetrating isotopic distributions, then transport downstream as they mix with other material.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The total mass of iron-group elements calculated by RRH09 for 0.2 M_⊙ (0.6 M_⊙) WDs ranged anywhere from negligible up to 0.034 M_⊙ (0.0003 M_⊙) for the BH masses we have considered, peaking at tidal strength β = 12 with a 0.2 M_⊙ WD, and at β = 5 with a 0.6 M_⊙ WD. Interestingly we find significantly greater iron-group production in our calculations: 0.089 M_⊙ (0.37 M_⊙) for 0.2 M_⊙ (0.6 M_⊙) WDs at β = 11 (8.9). We thus observe three times greater iron production from 0.2 M_⊙ WDs, and at least an order of magnitude more from 0.6 M_⊙ WDs even after scaling down the tidal strength to the values quoted in RRH09.

We have additionally considered interactions with large periapsis radii and small to moderate β in order to explore the possibility of producing calcium-rich debris and evaluate whether TDEs are viable sources of observed calcium-rich gap transients (Sell et al. 2015), despite the fact that many of these trajectories are further from the black hole and do not require full relativistic treatment. As Table 2 shows, we have systematically varied parameter sets to approach the low-density regime that is most likely to produce conditions capable of prematurely terminating nuclear reactions and give rise to dominant calcium-group elements (Ca, Ti, Cr) or IMEs rather than iron or nickel. Holcomb et al. (2013) calculated the smallest critical spatial scale ξ_crit that is required for helium detonation at these densities, finding ξ_crit ≈ 10⁹ (10⁷) cm at ρ ≈ 10⁵ (10⁶) g cm⁻³, well above our grid resolution. They argue that incomplete burn occurs if the size of the fuel reservoir is comparable to ξ_crit, which it is at densities approaching 10⁵ g cm⁻³. We find, as shown in Table 2, a direct correlation of systematically weaker reactions as this density is approached, but a regime clearly exists where Ca-group elements dominate over iron and nickel production before nuclear ignition eventually fails at large R_P. Table 2 suggests that conditions that give rise to Ca-rich debris are not particularly efficient at nucleosynthesis. So although calcium-rich outflows can indeed be produced by TD with IMBHs at low interaction strengths β ≲ 5, the resulting outflows, at least for the encounter parameters we have considered, will consist of ≲0.03 M_⊙ in the form of IMEs, converting roughly 14% of available fuel mass.

In addition to distant encounters with 10³ M_⊙ BHs, calcium-rich transients may also arise from close encounters with more massive 10⁴ M_⊙ BHs. The range of tidal strengths accommodated by more massive BHs is severely restricted and consequently less energetic from a nucleosynthesis perspective. Stars cannot be initialized with periapsis radius closer than about four Schwarzschild radii, or tidal strengths greater than β ∼ 5, without being promptly captured by the black hole (see Figure 2). This is a regime of relatively weak interaction where we observe little nuclear activity, so we limit the number of case studies involving 10⁴ M_⊙ BHs to just two calculations. These two cases, however, prove to be excellent candidates for calcium-rich gap transients precisely because of their unique property of sampling both weak tidal strengths and ultraclose encounters. In fact, case B4M2R04 is an example of a moderate strength interaction (β = 5.4) at the closest perihelion radius we have considered, and produces the greatest mass fraction of calcium-group elements in a calcium-rich environment, as demonstrated in Figure 12. Figure 13 shows that these transient candidates develop similar layered compositions to their iron-rich counterparts, except the dense inner core is composed of calcium-group elements, surrounded by elements of smaller atomic mass (e.g., silicon).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 suggests that density–temperature phase tracks are fairly well converged at our highest resolution. We have therefore extracted time histories of density and temperature for the densest 10% of matter in the WD for a few cases, and run a much larger nuclear reaction network on those point-dimensional data to get a better picture of the full range of element synthesis in the dense WD core and to compare it against the results for the smaller inline network. The cases we have chosen for post-processing in this way include B3M2R06 and B3M6R06, the strongest reacting 0.2 M_⊙ and 0.6 M_⊙ WDs, respectively. Because we do not evolve the entire star, but only a sample of its dense core, we simplify the initial isotopic composition and set it to either pure helium or a 50% carbon/oxygen mixture for the 0.2 M_⊙ and 0.6 M_⊙ WDs, respectively. The data are evolved in time using the Torch code (Timmes 1999) to solve standalone 19- and 640-isotope reaction networks with time-dependent density and temperature track profiles as inputs, solving the network system to convergence over each discrete time interval before adjusting density and temperature for the next cycle. Of course, these track histories do not necessarily correspond to contiguous fluid lines, but they are nonetheless useful for probing reactive transitions in the dense core.

The two cases (B3M6R06 and B3M2R06) are similar in that the 640-isotope network converts more than 90% of nuclear fuel to iron-group elements, a percentage that is respectively 1.5 and 3 times greater than the global fractions shown in Table 2. Global fractions are understandably lower since they sample all stellar material at very different densities and temperatures beyond the core where iron elements predominantly reside. This is supported by independently post-processing the dense core data using the same 19-isotope network adopted in the simulations. The 19- and 640-isotope models predict roughly the same iron-group fraction, differing only by a few per cent and giving greater confidence in the ability of the smaller network to provide reasonable predictions. Although the end product from the post-processing comes out almost entirely as ⁵⁶Ni, the 3D simulations predict that 84% (B3M6R06) and 71% (B3M2R06) of iron-group elements are in the form of ⁵⁶Ni, most likely attributed to the treatment of photodisintegration. At the high temperatures achieved in these ultraclose encounters, photodisintegration is very efficient. Burn products that form while the WD is undergoing initial compression are immediately disassembled when temperatures approach (6–7) × 10⁹ K. Iron-group elements do not form in substantial quantity until stellar material decompresses and cools to where α-capture can proceed uninhibited. Although this effect is observed in some of the species plots from the 3D simulations, it is not as dramatic since those data represent globally integrated mass fractions and thus sample a range of temperatures at any single time. It is in these photoprocesses where the small and large networks differ the most: the 19-species model tends to overestimate photodisintegration and produce more stable iron isotopes at the expense of nickel. Hence the nickel fractions found in the 3D simulations and quoted above are likely a lower bound to what a larger network might predict.

Figures 14 and 15 compare production factors as a function of atomic mass from the 3D simulations (red symbols) against the 640-isotope Torch model applied to the dense core tracks (black symbols). All data points represent the final mass fractions in each run, accounting for radioactive decay collapsing all species with a common mass number to their stable form. These are then normalized to the mass fraction in the Sun of the stable species. Connecting lines indicate isotopes of a given element (listed in the figures); a star indicates the most abundant isotope of a given element. Circles represent isotopes that were made as themselves, while triangles indicate that the stable species was made as a radioactive progenitor. As reactions proceed along the proton-rich side of the α-chain, they produce stable nuclei up to ⁴⁰Ca, after which unstable nuclei are formed, beginning with ⁴⁴Ti and produced in our calculations at very high relative (to solar) abundance as indicated in both figures. Notice that although calcium represents a small fraction of elements produced in the 3D simulations, it is comparable to iron and nickel when normalized by solar abundance. Notice also that the extracted phase track results (represented by black symbols) come out almost entirely as iron-group elements and cannot properly predict the creation of IMEs forming outside the WD core since they sample only the densest parts of the WD. This point was made in a previous paragraph, but is emphasized nicely in both Figures 14 and 15. In general the post-processed results suggest conditions similar to Type Ia supernovae (SNe) and to incomplete silicon burning and nuclear statistical equilibrium in Type II SNe: the dominant species produced are ⁵⁶Fe and ⁵⁷Fe made as radioactive nickel, both of which have been observed through γ-line emission in SNe Ia and SNe II. Yet heavier species of cobalt and nickel (made as copper and zinc) are co-produced, but the interesting species ⁶⁴Zn (made as ⁶⁴Ge) is not made in amounts that would suggest tidally disrupted WDs as a possible production site.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Though many details are lost with the 19-isotope model adopted in this work, this exercise demonstrates, among other things, that it is a fairly good approximation to the overall solution and accurately predicts the dominant end products. But it is also clear that significant nuclear reactions take place outside the dense WD core where IMEs are most likely to reside, and that a more accurate post-processing treatment will require fluid tracers to resolve nucleosynthesis on a global scale across the entire WD.

TD events are potential sources of strong gravitational waves, particularly during the early phase when the star is sufficiently compact and the interaction highly asymmetric. After passing periapsis, however, the resulting tidal forces disrupt the star into an extended debris tail and the wave signal diminishes significantly. Gravitational wave emissions from TDEs will therefore have burst-like behavior with estimated amplitude (RRH09)

$$\begin{eqnarray}\begin{array}{rcl}h & \approx & \displaystyle \frac{{{GM}}_{{\rm{WD}}}{R}_{{\rm{S}}}}{{c}^{2}{R}_{{\rm{P}}}D}\approx 1\times {10}^{-22}\beta {\left(\displaystyle \frac{D}{10{\rm{Mpc}}}\right)}^{-1}\\ & & \times \,{\left(\displaystyle \frac{{M}_{{\rm{WD}}}}{0.6{M}_{\odot }}\right)}^{4/3}{\left(\displaystyle \frac{{R}_{{\rm{WD}}}}{{10}^{9}{\rm{cm}}}\right)}^{-1}{\left(\displaystyle \frac{{M}_{{\rm{BH}}}}{{10}^{3}{M}_{\odot }}\right)}^{2/3},\end{array}\end{eqnarray} \tag{ 25 }$$

and frequency

$$\begin{eqnarray}&&f\approx {\left(\displaystyle \frac{{{GM}}_{{\rm{BH}}}}{{R}_{{\rm{P}}}^{3}}\right)}^{1/2}\approx 0.1\,{\beta }^{3/2}{\left(\displaystyle \frac{{M}_{{\rm{WD}}}}{0.6{M}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{{R}_{{\rm{WD}}}}{{10}^{9}{\rm{cm}}}\right)}^{-3/2}\,{\rm{Hz}},\end{eqnarray} \tag{ 26 }$$

placing the encounter scenarios we have considered within the anticipated frequency range of LISA (Amaro-Seoane et al. 2017).

More precisely, we calculate the retarded wave amplitudes of both polarizations, h₊ and h_×, for an observer along the z-axis at distance D from the source using the quadrupole approximation

$$\begin{eqnarray}&&{{Dh}}_{+}=\displaystyle \frac{G}{{c}^{4}}({\ddot{I}}_{{xx}}-{\ddot{I}}_{{yy}}),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{Dh}}_{\times }=\displaystyle \frac{2G}{{c}^{4}}{\ddot{I}}_{{xy}},\end{eqnarray} \tag{ 28 }$$

where ${\ddot{I}}_{{ij}}$ is the second time derivative of the reduced quadrupole moment of the mass distribution. In practice, the second time derivatives are eliminated through substitution of the hydrodynamic evolution equations and integration by parts, to produce an expression easily calculated from evolved fields (Camarda et al. 2009). We plot in Figure 16 a sample of computed waveforms that showcase their diversity with different combinations of WD and BH mass and interaction. Figure 17 plots the strain amplitude (assuming a source distance of 10 Mpc) and frequency for most of our simulations. The strain amplitude is computed as $\sqrt{{h}_{+}^{2}+{h}_{\times }^{2}}$, while the frequency is defined simply as the inverse of twice the time interval between the first minimum and first maximum in the h_× signal. The solid lines in Figure 17 represent the estimated dependences predicted by Equations (25) and (26). We find that the β¹ and β^3/2 analytic scalings for the amplitude and frequency, respectively, are, to within small normalization factors, well represented by the data at sufficiently small β, but deviate significantly at large tidal strengths (β > 10). Although our parameter space does not sample black hole or WD masses sufficiently to verify their analytic scalings (we considered only two values for each), we do note that both strain and frequency scalings are generally consistent with predictions, increasing with both M_WD and M_BH. We also note that the addition of nuclear burning does not significantly affect wave signals. We find, for instance, that wave amplitudes with and without active nucleosynthesis match to within a fraction of 1%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Gravitational wave attributes for all events we considered are packed tightly between strain amplitudes of 0.5 × 10⁻²² to 3.5 × 10⁻²² at a source distance of 10 Mpc and between frequencies of 0.1–0.7 Hz, just below and outside the frequency range of LIGO (Abbott et al. 2016). Although these signals fall within LISA's frequency band, the amplitudes are sufficiently small that they will not likely be observed except at source distances within 10–100 kpc, which, if current estimates of the WD–IMBH disruption rate (500 yr⁻¹ Gpc⁻³) are correct, would be extremely rare events.