Tidal disruption discs formed and fed by stream–stream and stream–disc interactions in global GRHD simulations

Andalman, Zachary L; Liska, Matthew T P; Tchekhovskoy, Alexander; Coughlin, Eric R; Stone, Nicholas

doi:10.1093/mnras/stab3444

ABSTRACT

When a star passes close to a supermassive black hole (BH), the BH’s tidal forces rip it apart into a thin stream, leading to a tidal disruption event (TDE). In this work, we study the post-disruption phase of TDEs in general relativistic hydrodynamics (GRHD) using our GPU-accelerated code h-amr. We carry out the first grid-based simulation of a deep-penetration TDE (β = 7) with realistic system parameters: a black hole-to-star mass ratio of 10⁶, a parabolic stellar trajectory, and a non-zero BH spin. We also carry out a simulation of a tilted TDE whose stellar orbit is inclined relative to the BH midplane. We show that for our aligned TDE, an accretion disc forms due to the dissipation of orbital energy with ∼20 per cent of the infalling material reaching the BH. The dissipation is initially dominated by violent self-intersections and later by stream–disc interactions near the pericentre. The self-intersections completely disrupt the incoming stream, resulting in five distinct self-intersection events separated by approximately 12 h and a flaring in the accretion rate. We also find that the disc is eccentric with mean eccentricity e ≈ 0.88. For our tilted TDE, we find only partial self-intersections due to nodal precession near pericentre. Although these partial intersections eject gas out of the orbital plane, an accretion disc still forms with a similar accreted fraction of the material to the aligned case. These results have important implications for disc formation in realistic tidal disruptions. For instance, the periodicity in accretion rate induced by the complete stream disruption may explain the flaring events from Swift J1644+57.

accretion, accretion discs, black hole physics, MHD, software: simulations, transients: tidal disruption events

1 INTRODUCTION

In recent decades, several very bright flares in galactic nuclei have been observed and interpreted as tidal disruption events (TDEs), which occur when a star is scattered on to a nearly parabolic orbit around a supermassive black hole (BH) with a pericentre inside the tidal radius of the BH (Hills 1975; Frank & Rees 1976; Rees 1988). While these flares are typically discovered from quasi-thermal emission in the soft X-ray (Bade, Komossa & Dahlem 1996; Komossa & Bade 1999; Saxton et al. 2012), UV (Gezari et al. 2006, 2008), or optical (van Velzen et al. 2011; Gezari et al. 2012; Arcavi et al. 2014; Holoien et al. 2014) bands, they have been observed to emit radiation across the electromagnetic spectrum, from radio synchrotron (Zauderer et al. 2011; Alexander et al. 2017) to non-thermal hard X-rays and soft gamma rays (Bloom et al. 2011; Cenko et al. 2012; Brown et al. 2015).

Our current theoretical understanding of the tidal disruption process – the star’s first, terminal pericentre passage – is largely converged (Lacy, Townes & Hollenbach 1982; Carter & Luminet 1983; Guillochon & Ramirez-Ruiz 2013; Mainetti et al. 2017), at least for polytropic stars in Newtonian gravity. More recent simulations have explored how the immediate outcome of disruption depends on stellar spin (Golightly, Coughlin & Nixon 2019a; Kagaya, Yoshida & Tanikawa 2019), realistic models of the star’s internal structure (Golightly, Nixon & Coughlin 2019b; Ryu et al. 2020a), and general relativistic gravity (Gafton et al. 2015; Tejeda et al. 2017; Gafton & Rosswog 2019; Ryu et al. 2020b). However, we do not yet have a first principles understanding of how, or if, the stellar debris streams are able to form a nearly axisymmetric, or quasi-circular, accretion disc (for a recent review of this problem see Bonnerot & Stone 2020). Because the bound stellar debris has typical eccentricities 0.99 ≲ e ≲ 0.999 (Stone, Sari & Loeb 2013), an enormous excess of orbital energy must be dissipated for circularization to occur.

Early work conjectured that most of this energy dissipation arises from relativistic apsidal precession (Rees 1988): as the most tightly bound debris passes through pericentre, its apsidal angle measured in radians precesses by an order-unity amount, causing a large-angle collision with less tightly bound matter that has yet to return to pericentre. The shocks that thermalize bulk kinetic energy at the point of self-intersection offer a plausible mechanism for circularizing returning stellar debris (Hayasaki, Stone & Loeb 2013; Bonnerot et al. 2016). However, self-intersection shocks may be less efficient at circularizing the debris for inclined orbits around spinning Kerr BHs. In this regime, nodal precession, due to Lense–Thirring frame dragging may delay the onset of self-intersection by many orbits (Cannizzo, Lee & Goodman 1990; Kochanek 1994; Guillochon & Ramirez-Ruiz 2015; Hayasaki, Stone & Loeb 2016). Additionally, energy dissipation due to self-intersection shocks may be greatly limited for less relativistic orbital pericentres, with small-angle collisions occurring at self-intersection radii near the apocentre of the most tightly bound debris (Dai, McKinney & Miller 2015; Shiokawa et al. 2015).

An alternative dissipation site is at the stream pericentre itself, where the recompression of the returning debris generates ‘pancake’ or ‘nozzle’ shocks (Kochanek 1994). Recently, Bonnerot & Lu (2021) performed an in-depth study of the nozzle shock using a 2D simulation of a vertical slice of the stream and found that dissipation at the nozzle shock raises the entropy of the gas by two orders of magnitude. Newtonian hydrodynamic simulations by Ramirez-Ruiz & Rosswog (2009) have shown that this pericentre shock could feasibly circularize the tidal debris. However, these simulations were performed for a BH-to-star mass ratio of Q = 10³ and analytic estimates of Guillochon, Manukian & Ramirez-Ruiz (2014) suggest that pericentre recompression shocks become less efficient at more realistic mass ratios (e.g. Q = 10⁶). It is also possible that in many TDEs, efficient dissipation is lacking altogether, and the formation of an accretion disc is an inefficient process unfolding over many fallback times (Piran, Sa̧dowski & Tchekhovskoy 2015). We discuss further in Section 4.1.

TDE debris circularization and disc formation is a complex physical problem involving a large dynamic range, general relativistic orbital dynamics, the need for accurate treatment of hydrodynamic shocks, and possibly even magnetohydrodynamic (MHD) effects (Svirski, Piran & Krolik 2017). The many pieces of multiscale and non-linear physics involved in TDE disc formation mean that, for numerical reasons, almost all past simulations of this process employed major simplifying assumptions that cast doubt on the generality of their conclusions.

Ayal, Livio & Piran (2000) initiated the numerical study of TDE circularization using a post-Newtonian (PN) potential to simulate the lowest order level of apsidal precession in a finite-mass, smoothed particle hydrodynamics (SPH) framework, albeit with low (N ∼ 10³) particle number. More recently, global circularization simulations achieved much higher resolution by reducing the dynamic range of the problem in one of two ways. The first is to consider an unrealistically low mass ratio, typically Q ∼ 10³. In simulations of this type, general relativity is sometimes ignored completely (Guillochon et al. 2014), but when it is included, it has a minimal effect on the circularization process because the tidal radius around an intermediate-mass BH is not very relativistic (Ramirez-Ruiz & Rosswog 2009; Shiokawa et al. 2015).

The second option is to consider a realistic mass ratio (Q ∼ 10⁶) but an unrealistic pre-disruption stellar orbit. Tidally disrupted stars typically approach supermassive BHs on nearly parabolic orbits (Magorrian & Tremaine 1999), with initial eccentricities 1 − e₀ ∼ 10⁻⁵. For computational convenience, one may choose an unrealistic stellar eccentricity, e₀ ≲ 0.95, to reduce the debris stream apocentres. This approach was adopted by Hayasaki et al. (2013), who mimicked apsidal precession effects with a pseudo-Newtonian potential in a highly relativistic β = 5 TDE. They found rapid circularization due to orbital energy dissipation at stream self-intersections. Later simulations found that in less relativistic β = 1 TDEs, self-intersections are less efficient at energy dissipation compared to this initial work (Bonnerot et al. 2016).

The results of Hayasaki et al. (2013) were later confirmed and extended to different gas equations of state (Bonnerot et al. 2016; Hayasaki et al. 2016), as well as higher (but still sub-parabolic) eccentricities (Bonnerot et al. 2016; Sądowski et al. 2016). The low-e₀ limit of tidal disruption has also been used with PN potentials to include Lense–Thirring frame dragging, which was seen to substantially delay circularization provided debris streams remain thin (Hayasaki et al. 2016). More recently, Bonnerot & Lu (2020) have performed a TDE disc formation simulation with realistic astrophysical parameters using a different approximation: neglecting the returning debris streams entirely, and injecting mass, momentum, and energy (in the form of SPH particles) from the test-particle self-intersection point. The validity of this approach depends on the accuracy of the local injection scheme, and its independence from global gas evolution around the BH. We discuss this approach further and compare and contrast it to our results in Section 4.5.2.

In this paper, we use novel numerical techniques to capture the disc formation process in general relativistic hydrodynamics without sacrificing astrophysical realism in our choice of system parameters (e.g. Q, e₀). We use two-level adaptive mesh refinement (AMR) to resolve the relevant physics within our grid-based code. In Section 2, we outline our numerical scheme. In Section 3, we describe the general outcomes of our simulation, including the spatial properties of the nascent accretion flow. In Section 4, we more carefully analyse the specific physical mechanisms controlling the accretion and circularization process and provide a detailed comparison to the ZEro-BeRnoulli Accretion (ZEBRA) model of Coughlin & Begelman (2014). We conclude in Section 5.

2 NUMERICAL METHOD AND SETUP

We simulate the initial tidal disruption using the SPH code phantom (Price & Monaghan 2007) and we simulate the post-disruption evolution using our new GRMHD code h-amr (Liska et al. 2019), an approach analogous to those of Rosswog, Ramirez-Ruiz & Hix (2009) and Sądowski et al. (2016). With this hybrid method, we can account for the large range of spatial and temporal scales involved in the disruption process and debris stream formation while accurately capturing the essential shocks and general relativistic effects in the post-disruption evolution.

2.1 Initial disruption in phantom

The stellar disruption is initially followed by the smoothed-particle hydrodynamics code phantom (Price et al. 2018). The setup is identical to that described in Coughlin & Nixon (2015): a star of mass |$1\, {\rm M}_{\odot }$| is modelled as a γ = 5/3 polytrope, with the adiabatic index equal to the polytropic index, by placing ∼10⁷ particles on a close-packed sphere. The sphere is stretched to achieve roughly the correct polytropic density profile. The star is subsequently relaxed in isolation (i.e. without the external gravitational potential of the BH) for ten sound crossing times to smooth out numerical perturbations in the density profile. Self-gravity is included through the implementation of a tree algorithm alongside an opening angle criterion to calculate short-range forces (Gafton & Rosswog 2011). We also include the effects of shock heating in modifying the internal energy of the gas.

The relaxed polytrope is placed at a distance of |$5\, r_{\rm t}$| from the supermassive BH of mass |$10^{6}\, {\rm M}_{\odot }$| such that the centre of mass is on a parabolic orbit. To maintain hydrostatic balance initially, every particle comprising the star is given the velocity of the centre of mass. In its current version, phantom is a Newtonian code, and therefore has no direct means of implementing general relativistic effects. Instead, we mimic some of these effects with a pseudo-Newtonian ‘Einstein’ potential used by Nelson & Papaloizou (2000), Nealon, Price & Nixon (2015), Xiang-Gruess, Ivanov & Papaloizou (2016). The potential is given by

$$\begin{eqnarray*} \Phi _{\rm E} = -\frac{GM_{\rm BH}}{r}\left(1+\frac{3R_{\rm g}}{r}\right), \end{eqnarray*}$$

(1)

where R_g = GM_BH/c² is the gravitational radius and M_BH is the mass of the BH. This potential accurately reproduces the general relativistic apsidal precession angle at large radii relative to the gravitational radius, with deviations from the true precession angle becoming more pronounced as the radius r becomes comparable to R_g. However, for the large mass ratio considered here, the tidal approximation is upheld to a high degree of accuracy, meaning that the dominant effect of general relativity on the initial stellar encounter will be to rotate the entire star through the same precession angle.¹ Therefore, our usage of this potential, as opposed to a fully general relativistic treatment, is sufficient for the purpose of creating a realistic distribution of post-disruption debris. See Tejeda & Rosswog (2013) for a detailed evaluation of a similar potential used by Nowak & Wagoner (1991).

The initial, parabolic orbit of the star is established using the above potential (equation 1) to calculate the angular momentum necessary to achieve a pericentre distance of r_p = 7R_g. phantom uses an artificial viscosity prescription to mediate any strong shocks that may be present during the large compression suffered by the star and employs the standard switch proposed by Cullen & Dehnen (2010) (i.e. the artificial viscosity parameter is small when the star is far from pericentre and approaches values near unity as the star is compressed at pericentre). A non-linear term is also included to account for extremely strong shocks and prevent interparticle penetration (Price & Federrath 2010). The large number of particles (∼10⁷) was used to avoid the possibility of spurious numerical heating at pericentre caused by under-resolving the compression, predicted to be of the order of H_min/R_* ∼ β⁻³ ∼ 0.003 (Carter & Luminet 1983), though the compression could be smaller if shock heating halts the otherwise-adiabatic collapse (Bicknell & Gingold 1983), or due to 3D effects (Guillochon et al. 2009). Here β = r_t/r_p is the penetration factor, r_t is the tidal radius, and r_p is the pericentre radius.

Fig. 1 shows the density distribution of the disrupted stellar debris at 1.16 d after the disruption. At this time, we end the evolution of the TDE in phantom and use the distribution of debris as the initial conditions for our post-disruption simulation in h-amr. See Section 2.2 for a detailed description of how we map data from SPH to grid-based GRHD.

Figure 1.

Colour maps of the log of rest mass density and log κ (proportional to entropy, see equation 7), in the equatorial plane at the initial conditions of the post-disruption phase of the simulation in h-amr (1.16 d). The black contour on the right-hand panel outlines the area excluded by the entropy condition (κ < 10) which we use throughout our analysis to distinguish the material in the debris stream from the material in the accretion disc). The BH is located at the origin. The dotted lines indicate the x- and y-axes.

Open in new tab Download slide

Fig. 2 depicts the Bernoulli parameter of the tidally disrupted debris at this same time. The star approaches the BH on a parabolic orbit. When the BH tidally disrupts the star, the bound stellar debris falls back to the BH while the unbound debris continues on an outward trajectory. The star’s mass is split almost evenly between bound and unbound matter (Lacy et al. 1982; Evans & Kochanek 1989). We use the relativistic Bernoulli parameter to distinguish between bound and unbound material.

$$\begin{eqnarray*} b=-\frac{u_t(\rho c^2 + u_g \gamma)}{\rho c^2}-1 , \end{eqnarray*}$$

(2)

where b is the Bernoulli parameter, u is the 4-velocity, and ρ and u_g are the mass and internal energy densities in the fluid frame. At late times the pressure gradient within the disrupted debris becomes small and the Bernoulli parameter is approximately a conserved Lagrangian quantity. Even in a time-dependent system such as the one analyzed here, b > 0 corresponds to unbound material while b < 0 corresponds to bound material.

Figure 2.

A histogram of the Bernoulli parameter distribution at the initial conditions of the post-disruption phase of the simulation in h-amr (1.16 d). Each bin is weighted by solar masses and bin width. Unbound material and total unbound mass are shown in red; bound material and total bound mass are shown in blue. The mass-weighted average Bernoulli parameter is also shown. On average, the material in the initial conditions of the post-disruption phase is marginally bound (b = 0). The vertical lines represent the range of the Bernoulli parameter estimated from the frozen-in approximation (equation 3), and contain 98.1 per cent of the debris mass. The floor material is ignored using a density condition (ρ > 10⁻¹¹).

Open in new tab Download slide

Fig. 2 shows that in the initial conditions, the majority of the material is marginally bound with Bernoulli parameter inside the range predicted by the frozen-in approximation of Stone et al. (2013). This approximation assumes that (i) the tidal forces outside of the tidal radius are negligible, so the star enters the tidal radius as an unperturbed sphere, and (ii) once the star crosses the tidal sphere, its fluid elements move ballistically, with a spread in orbital properties given by the potential gradient across the star. In reality, internal forces (e.g. self-gravity and hydrodynamics) are not totally negligible inside the tidal sphere, but previous simulations of deeply penetrating disruptions show that the frozen-in approximation reproduces the actual energy spread of the debris to within ≈20 per cent for γ = 5/3 polytropes (Steinberg et al. 2019). According to this impulsive disruption approximation, the spread of specific orbital energy Δb in Newtonian gravity is given by

$$\begin{eqnarray*} \Delta b = k \frac{GM_{\rm{BH}}R_*}{r_{\rm t}^2} , \end{eqnarray*}$$

(3)

where k is a constant of the order of unity related to stellar structure and rotation prior to disruption. If we let k = 1, we find that Δb = 2.12 × 10⁻⁴. Only 1.94 per cent of the mass in the initial conditions is outside the range predicted by the frozen-in approximation, verifying that the initial orbital energy distribution for the post-disruption phase is largely consistent with standard estimates.

A small fraction of the material has Bernoulli parameters well outside the range predicted by the frozen-in approximation. However, even though the most tightly bound debris (with specific energy |ε| > Δb) constitutes a small fraction of the total mass, it is the first matter to fall back, and therefore dominates the early stages of the circularization process studied here. Due to runtime limitations, these early stages are the primary focus of this paper. While these tails could be a byproduct of intense shock heating as the star is highly compressed near pericentre, we caution that they may also arise from numerical inaccuracies associated with the same highly compressed, and therefore difficult to resolve, configuration of gas.

Such broad-energy tails have been seen in high-β TDEs simulated with a range of codes and numerical algorithms. While a return time of 1.16 d for the most tightly bound debris might appear extreme, it is qualitatively consistent with these past simulations. For example, Guillochon & Ramirez-Ruiz (2013) find that the most tightly bound debris in Newtonian, grid-based, β = 4 simulations of n = 3 polytrope disruptions can return to pericentre after ≈3 d and that the time of first pericentre return decreases with increasing β. Steinberg et al. (2019) performed moving-mesh simulations of stellar disruptions and found that the extent of the high energy tail is also a function of β. For Newtonian disruptions of n = 3/2 polytropes, going from β = 5 to β = 7 moves the time of first mass return from ≈3 d to ≈1 d (Steinberg, private communication). Gafton & Rosswog (2019) used Newtonian and relativistic SPH simulations, with a code distinct from phantom, to disrupt a γ = 5/3 polytrope over a range of β and found that for large β the return time of the most bound debris was significantly earlier than the frozen-in prediction, with initial return times on the order of days.

These high-energy, low-mass debris tails have not been studied in detail, but their ubiquity across SPH, conventional grid-based, and moving mesh codes leads us to believe that they are likely physical. If, however, the high-energy tail of debris were primarily the result of numerical artefacts, then it would bias the earliest stages of mass return to (i) artificially early times and (ii) artificially low fallback rates relative to the time of first mass return.

However, our results depend solely on the relative values of mass fluxes rather than the absolute values, with the exception of the internal energy and density floors (Section 2.3). Therefore, even if the mass of the high-energy tail of debris in our initial conditions is an overestimate, our results can be straightforwardly rescaled to astrophysically realistic time and mass flux scales (i.e. our simulation would have started at a later time with similar values of relative mass flux). The qualitative features of the circularization process are therefore robust and should apply generically to systems with realistic physical parameters and β ≃ 7.

2.2 Mapping from phantom to h-amr

Before we begin our simulations in h-amr, we map SPH data from phantom to gridded data compatible with GRHD. First, smooth the particle properties on to a continuous domain. Secondly, we construct a grid on top of the smoothed particles.

For the first step, we use the SPH visualization tool splash and its built-in function ‘splash to grid,’ which is described in detail in Price (2007). The function smooths the density of each particle over a finite region according to a weighting function, or kernel, that is twice differentiable, maintains compact support, and decreases in magnitude from the location of the particle. This approach mirrors the standard procedure for SPH calculations. We use the default kernel in splash (and phantom): a cubic spline which vanishes at a distance of 2h from a given particle. The smoothing length h is spatially variables and set such that the mass inside the smoothing sphere is constant. We refer the reader to section 2.4 of Price (2012) for more details.

For the second step, the fluid variables of a given cell are determined by adding the contribution of every particle with a smoothing region that encompasses the cell itself. See fig. 7 of Price (2007) for an illustration. This method ensures that cells in high-density regions are sampled by a large number of particles and cells in low-density regions are sampled by relatively few particles. Because of their low density, sparsely sampled cells contribute minimally to the dynamics of the fluid and do not affect the bulk properties of the accretion flow simulated with h-amr.

The error in the first step is of the order of |$\mathcal {O}(h^2) = \mathcal {O}(\rho ^{-2/3})$|⁠. Therefore, in regions of extremely low density, the interpolation from SPH may introduce artefacts in the density and velocity fields. In particular, we find a region of enhanced density in the region around the black hole (Fig. 1). The density of this region is several orders of magnitude lower than the density of the most bound debris, so the trajectory of the returning debris is not significantly altered.

2.3 h-amr simulation parameters

As described in Section 1, our simulation parameters adhere to astrophysically realistic values (Q = 10⁶, e₀ ≈ 1). In phantom we insert a star on a parabolic trajectory with a pericentre distance of 7R_g and a penetration factor of β = r_t/r_p = 7. This high penetration encounter guarantees that self-gravity is negligible in the post-disruption evolution of the stream, although the influence of self-gravity on the stream structure may be revived at much later times than those simulated here due to the in-plane focusing of the debris; Coughlin et al. 2016; Steinberg et al. 2019). The circularization of the stellar debris likely occurs on a shorter time-scale for more relativistic encounters (Bonnerot & Stone 2020), decreasing the simulation duration required to study the circularization process.

h-amr uses a naturalized unit system where G = c = R_g = 1. The conversion factors from the simulation units to cgs units are given in Table 1. From now on, we will work in this naturalized unit system with the exception of time, which we will convert back to physical units of days since disruption. We will also restore G and c in equations to help keep track of units.

Although h-amr does not explicitly include viscosity, it is included implicitly through interactions at the cell level. Within the turbulent flows of our simulation, adjacent fluid elements are unlikely to move exactly parallel to one another and therefore will exchange momenta. Previous work has shown that the effective viscosity in early-time TDE accretion flows can be dominated by the Reynolds stress (Sądowski et al. 2016). If this is correct, we should not be significantly underestimating effective viscosity due to the absence of magnetic fields, though this question needs closer examination in future magnetohydrodynamic simulations.

We present two models, TDET0 and TDET30, corresponding to spin-orbit misalignment angles of zero and 30 deg, respectively (Table 2; see the 3D renderings in the Supporting Information). Because we begin the simulation in h-amr 1.16 d post disruption, the relationship between h-amr simulation time and time since disruption is affine linear.

$$\begin{eqnarray*} t_{\rm HAMR} = (t_{\rm days} - 1.16) \times 24 \times 3600 / 4.926 \end{eqnarray*}$$

(4)

Table 1.

Open in new tab

For each quantity, the number of cgs units per simulation unit is tabulated. Note that 0 R_g/c corresponds to 1.16 d after the disruption, so the relationship between simulation time and days since disruption is affine linear.

Quantity	cgs unit	h-amr unit	cgs unit / h-amr unit
Mass	g	R_gc²/G	2 × 10³⁹
Distance	cm	R_g	1.477 × 10¹¹
Time	s	R_g/c	4.926
Density	g cm⁻³	\|$c^2/R_{\rm g}^2 G$\|	6.207 × 10⁵

Quantity	cgs unit	h-amr unit	cgs unit / h-amr unit
Mass	g	R_gc²/G	2 × 10³⁹
Distance	cm	R_g	1.477 × 10¹¹
Time	s	R_g/c	4.926
Density	g cm⁻³	\|$c^2/R_{\rm g}^2 G$\|	6.207 × 10⁵

Table 1.

Open in new tab

For each quantity, the number of cgs units per simulation unit is tabulated. Note that 0 R_g/c corresponds to 1.16 d after the disruption, so the relationship between simulation time and days since disruption is affine linear.

Quantity	cgs unit	h-amr unit	cgs unit / h-amr unit
Mass	g	R_gc²/G	2 × 10³⁹
Distance	cm	R_g	1.477 × 10¹¹
Time	s	R_g/c	4.926
Density	g cm⁻³	\|$c^2/R_{\rm g}^2 G$\|	6.207 × 10⁵

Quantity	cgs unit	h-amr unit	cgs unit / h-amr unit
Mass	g	R_gc²/G	2 × 10³⁹
Distance	cm	R_g	1.477 × 10¹¹
Time	s	R_g/c	4.926
Density	g cm⁻³	\|$c^2/R_{\rm g}^2 G$\|	6.207 × 10⁵

Table 2.

Open in new tab

Simulation parameters for models TDET0 and TDET30, including black hole mass M_BH, stellar mass M_*, pericentre radius R_p, penetration factor β, and inclination angle of the stellar orbit i.

Model	M_BH (M_⊙)	M_* (M_⊙)	R_p (R_g)	β	i
TDET0	10⁶	1	7	7	0
TDET30	10⁶	1	7	7	30

Table 2.

Open in new tab

Simulation parameters for models TDET0 and TDET30, including black hole mass M_BH, stellar mass M_*, pericentre radius R_p, penetration factor β, and inclination angle of the stellar orbit i.

Model	M_BH (M_⊙)	M_* (M_⊙)	R_p (R_g)	β	i
TDET0	10⁶	1	7	7	0
TDET30	10⁶	1	7	7	30

In both models, we use a dimensionless BH spin of a = 0.9375 for the post-disruption evolution. Because the morphology and the fallback rate of our stream are not strongly dependent on spin (Tejeda et al. 2017), we rotate the initial data in h-amr about the y-axis by 30 deg for model TDET30 rather than repeating the simulation in phantom. We take this approach, rather than tilting the metric, to avoid the computational strain associated with a non-axisymmetric metric.

We run models TDET0 and TDET30 until 6.87 d (t = 10⁵R_g/c in h-amr) and 5.01 d (t = 6.7 × 10⁴R_g/c in h-amr) after the disruption, respectively. We evolve the models in the Kerr geometry using Kerr-Schild coordinates to avoid the coordinate singularity in the Boyer-Lindquist coordinates.

In this work, h-amr uses 2-level 3D adaptive mesh refinement (AMR) with a refinement criterion based on a threshold density. The total effective resolution is 2880 × 860 × 1200 in r × θ × ϕ. The cells are logarithmically spaced in the radial direction. The cell dimensions in R_g are given approximately as a function of position around the black hole:

$$\begin{eqnarray*} \Delta r \times \Delta \theta \times \Delta \phi \approx \frac{\ln (10^5)}{2880} r \times \frac{\pi }{860} r \times \frac{2\pi }{1200} r\sin \theta . \end{eqnarray*}$$

(5)

Our grid resolves the vertical extent of the stream (twice the stream scale height) with ∼12 cells near 500 R_g and ∼28 cells near pericentre, where the stream scale height is computed via the same method as the disc scale height discussed in Section 3.1. However, near pericentre, this calculation is artificially inflated by the stream structure’s complex ϕ-dependence. This resolution is sufficient to model the gas dynamics of the stream. However, we caution that we may underresolve pressure gradients within the stream.

In addition to AMR, h-amr uses local adaptive time stepping by setting the time-step in each cell to the smallest light crossing time in that cell. As a result, the time-steps decrease by a factor of Δϕ/Δθ near the pole due to the time-step limitation in φ. Therefore, we use the cylindrification method described by Tchekhovskoy, Narayan & McKinney (2011) to reduce the azimuthal extent of cells near the pole.

As we ran the simulation, we noticed that the stream disintegrates into the accretion disc after wrapping around the BH, a behaviour which we discuss more in Section 3.1. To verify that this stream disintegration was not a numerical artefact, we adjusted the refinement criterion for first-order refinement between 2.28 d (4 × 10⁴R_g/c) and 4.56 d (8 × 10⁴R_g/c). Near the BH, we decreased the cutoff density for first-order refinement to achieve the full effective resolution in a greater fraction of cells. Due to memory restrictions, we also had to increase the cutoff density at large radii, causing some parts of the outer stream to become unrefined. We found that the stream disintegration and other simulation properties were consistent across adjustments to the refinement criterion, suggesting that the physics converged for our resolution in h-amr.

In the post-disruption phase, we set floors for internal energy density and mass density at 2.27 × 10⁻¹² (3.75 × 10⁶ ergs cm⁻³ and 4.167 × 10⁻¹⁵ g cm⁻³ respectively). We assume an adiabatic index γ = 5/3 corresponding to the gas-pressure-dominated regime present in the star before it undergoes shocks. Although any accretion disc resulting from the TDE is expected to be radiation-pressure-dominated with an adiabatic index closer to γ = 4/3, h-amr currently does not allow for a variable adiabatic index so we chose to stay consistent with the initial simulation in phantom. Many previous works have also used an adiabatic index of γ = 5/3 (Guillochon & Ramirez-Ruiz 2015; Sądowski et al. 2016; Steinberg et al. 2019). In future simulations, we will implement a variable adiabatic index to more accurately model the thermodynamics.

3 RESULTS

3.1 Aligned disc formation and evolution

General relativistic apsidal precession near the BH causes the argument of pericentre to advance by approximately 31°, setting the incoming and outgoing streams on trajectories that intersect at R_SI = 142R_g. The precession angle and self-intersection radius are computed with the formalism in Appendix A4.

Shortly after the first pericentre passage, depicted in the first row of Fig. 3, the debris stream expands significantly and its density drops. The primary reason for this expansion is the different degrees of apsidal precession experienced by the inner- and outer-most edges of the stream. For instance, at t = 1.25 d (top right-hand panel of Fig. 3), the inner and outer edges of the stream at pericentre lie at r = 11R_g and r = 21R_g respectively, resulting in a difference of roughly 15° in their precession angle (Fig. 4). Bonnerot & Lu (2021) discuss this effect as being potentially responsible for a small amount of spreading of the debris following its first return to the initial point of disruption. Here, because the pericentre distance of the star is much more relativistic than the modest-β case that they considered, the differential precession is larger than they find and apparent by eye.

Figure 3.

Contour plots of the log of rest mass density in the equatorial plane in the TDET0 simulation during the debris’ initial pericentre passage (first row) and the first (second row) and third (third row) self-intersection events. In the initial pericentre passage, the stream falls back through near vacuum and matter from the stream begins to accumulate near the BH. In the self-intersection events, the stream undergoes apsidal precession and self-intersects close to the analytical self-intersection radius at 142R_g. As a result, the inner parts of the stream are completely disrupted. These violent events are a key dissipation mechanism in the early stages of the TDE evolution. Although powerful, we count only five such events. At late times in our simulation, dissipation occurs primarily through interactions with the newly formed disc (Fig. 5). For a more complete picture of the disc evolution, see the movies and 3D renderings linked in Section 6.

Open in new tab Download slide

Figure 4.

A contour plot of the log of rest mass density in the equatorial plane at 1.25 d (1725R_g/c; same as the top right-hand panel of Fig. 3) with equatorial geodesics of various pericentre radii. The geodesics diverge after pericentre, demonstrating the strong effect of differential precession in expanding the stream in the deeply penetrating (β = 7) and high-spin (a = 0.9375) encounter we consider.

Open in new tab Download slide

Over the course of the simulation, the self-intersection periodically becomes powerful enough to fully intercept the incoming debris stream. During these violent self-intersection events, the collision of the incoming and outgoing streams and the associated shock heating nearly destroys the stream interior to the self-intersection point and ejects material into a wide range of orbits. Rows two and three of Fig. 3 depict the time evolution of two such intersection/depletion cycles.

We find five violent self-intersection events in our simulation. They occur approximately 12 h apart (at 1.47, 2.01, 2.52, 2.92, and 3.68 d), and each lasts for roughly 2.74 h (2, 000R_g/c). The periodicity of the self-intersections is on the scale of the free-fall time from the self-intersection point,

$$\begin{eqnarray*} t_{\rm {ff}}=\frac{\pi }{2}\frac{R_{\rm SI}^{3/2}}{\sqrt{2GM_{\rm{BH}}}}\approx 3.64\ \rm {h}. \end{eqnarray*}$$

(6)

These violent, periodic self-intersections events may produce natural, quasi-periodic variability in the inner disc accretion rate, possibly explaining the flaring events observed in TDEs such as SWJ1644+57 (Burrows et al. 2011; Zauderer et al. 2011), AT2018fyk (Wevers et al. 2019), and AT2019ehz (van Velzen et al. 2020). We discuss this hypothesis further in Section 4.3. The self-intersection events also create the initial accretion disc.

However, once the accretion disc becomes sufficiently dense and massive, no additional violent self-intersections occur because the outgoing stream completely disintegrates before the self-intersection point as in Fig. 5. At late times, the spread in angular momenta within the incoming stream grows due to angular momentum exchange with the accretion disc. As a result, the stream becomes thicker, enhancing the effect of differential precession discussed earlier in this section. The result is that the outgoing stream has a larger spread in trajectories and a lower density than the incoming stream. At even later times, the disc is sufficiently dense to absorb nearly all of the momentum from the weakened outgoing stream through shocks and hydrodynamic instabilities at the interface. In particular, the velocity difference between the outgoing stream and the disc leads to a Kelvin–Helmholtz instability (see Bonnerot & Stone 2020, section 2.2.4), which seeds turbulence in the outgoing stream causing it to break apart.

Figure 5.

Contour plots of the log of rest mass density density (left-hand panel) and κ (right-hand panel) in the equatorial plane at 5.7 d. At late times in the simulation, the debris stream undergoes shocks and instabilities near pericentre, causing it to disintegrate shortly after the pericentre passage. This is different from the stream’s early time evolution (Fig. 3), possibly because an inner accretion disc has formed with density comparable to the incoming stream. The black line depicts an equatorial geodesic (Appendix A2). Note that the self-intersection radius of the geodesic is much greater than the analytical self-intersection radius of 142R_g for a pericentre of 7R_g because the geodesic has a larger pericentre radius of ∼12R_g.

Open in new tab Download slide

At early times, especially before substantial disc formation, the pericentre radius undergoes fluctuations. From the movies linked in Section 6, we see that the outward movements of the pericentre coincide with the onset of a violent self-intersection. This suggests that the pericentre movement results from the azimuthal momentum added to the incoming stream by the outgoing stream during a self-intersection. After the self-intersections cease, the pericentre radius becomes more stable and takes on a value around 12R_g (corresponding to a self-intersection radius R_SI = 565R_g). This may be because the outgoing stream has a weaker, but more stable, contribution to the azimuthal momentum of the incoming stream at late times.

Instead of forming a standard, geometrically thin disc, the material surrounding the black hole is inflated into a geometrically thick structure that is both gas-pressure and centrifugally supported. We perform a more in-depth analysis of the force balance in the disc and the disc structure in relation to analytical models in Sections 4.4 and 4.5.

We use an entropy cutoff to distinguish between the matter in the stream and the matter in the disc. Throughout the remainder of this work, we use the quantity,

$$\begin{eqnarray*} \kappa =p\rho ^{-\gamma }=\frac{(\gamma -1)u_g}{\rho ^\gamma } \end{eqnarray*}$$

(7)

to track specific entropy, which is related to κ by

$$\begin{eqnarray*} \mathcal {S}=\frac{\rm{ln} \kappa }{\gamma - 1}. \end{eqnarray*}$$

(8)

The tidal compression of returning debris streams is approximately a reversible process, so entropy is nearly constant until the first shock. For the purposes of analysis, we define the stream as material with κ < 10, a definition we refer to as the entropy condition. Fig. 1 depicts an entropy profile of the stream at the initial conditions of the post-disruption phase in the equatorial slice. The black contour outlines the area covered by the entropy condition.

Fig. 6 shows the radial profiles of density, pressure, and φ-velocity within the disc. We compute the gas pressure p using the adiabatic equation of state,

$$\begin{eqnarray*} p=u_g(\gamma -1). \end{eqnarray*}$$

(9)

We compute the physical φ-velocity directly from the simulation as

$$\begin{eqnarray*} v_\varphi =\frac{u^{\phi }}{u^{t}} \sqrt{g_{\phi \phi }} , \end{eqnarray*}$$

(10)

where g is the metric tensor. For a given quantity |$\mathcal {Q}$|⁠, we compute the mass-weighted averages over two coordinates using

$$\begin{eqnarray*} \mathcal {Q}_{\rm{avg}}=\frac{\int {\mathcal {Q} \rho u^t {\rm d}A_{\mu \nu }}}{\int {\rho u^t {\rm d}A_{\mu \nu }}}, \end{eqnarray*}$$

(11)

where

$$\begin{eqnarray*} {\rm d}A_{\mu \nu }=\sqrt{- g}{\rm d}\mu {\rm d}\nu . \end{eqnarray*}$$

(12)

where g is the determinant of the metric tensor. Radial profiles have dA_μν ∝ dθdϕ and polar profiles have dA_μν ∝ drdϕ. In calculations involving radial averages, we restrict the region of integration radially to avoid capturing high-density material from the half of the star which escapes the black hole. The region of integration for a given calculation is described in more detail in the caption of the corresponding figure.

Figure 6.

Time averages of radial profiles of mass density, pressure, and φ-velocity, their power-law fits, and an inset plot of time-averaged rest mass density in the equatorial plane. Power-law fits (dashed lines) are calculated using a least-squares method (see Table B1 for more details). Time averages are over the simulation’s full duration. Mass density, pressure, and φ-velocity are averaged over spherical shells using equation (11). Pressure and φ-velocity are weighted by mass. The stream is ignored using the entropy condition. Mass, density, and pressure are multiplied by 5 × 10⁵ so that all three variables are roughly the same order of magnitude for ease of comparison. The vertical lines show the pericentre radius at 7R_g and the analytical self-intersection radius at 142R_g (Appendix A4). The analytical self-intersection radius is also shown on the inset plot. All three quantities follow power-law fits within the radii of the disc. The sub-unity coefficient on the φ-velocity indicates that the disc is sub-Keplerian. At radii less than the pericentre radius or greater than 400 R_g, there is minimal disc material, so the data at these radii does not reflect the large-scale properties of the disc. The origin of the flat density and pressure regions at large radii (r ≳ 1000R_g) is due to the choice of the floors. The origin of the dips in all quantities at small radii is due to the absence of disc material in the plunging region.

Open in new tab Download slide

Fig. 6 shows that the radial profiles (i.e. averaged over angles) of density and pressure closely follow power-law relationships between the inner and outer boundaries of the disc (10–200 R_g), hinting at a possible analytic description (see Section 4.5.1). The angle-averaged φ-velocity is fitted by v_φ ≃ 0.76r^−0.5, which implies a sub-Keplerian velocity distribution that may be due to thermal pressure support against gravity (see Section 4.4).

The internal energy density and mass density are floored at 2.27 × 10⁻¹² (see Section 2.3). These floors are responsible for the flat density and pressure regions at large radii in Fig. 6. While these floors would have a negligible effect on a TDE at peak fallback rate, they become significant for the early times and low fallback rates considered in our simulation. The floors may affect our results by providing external pressure support to the outer disc, artificially lowering its radial and vertical extent. We discuss this in Section 4.5.1.

Fig. 7 shows the polar profiles of density, pressure, and squared specific angular momentum within the disc. We calculate the pressure as above (equation 9) and the specific angular momentum as l = u_ϕ. The polar profiles of all three quantities are fit to a power α of sin ²θ near the equatorial plane. We analyse these relationships further and compare them to model predictions in Section 4.5.1.

Mass density, pressure, and angular momentum squared in our aligned TDET0 model plotted with respect to θ, their fits to a power law of sin 2θ, and an inset plot of rest mass density in the xz-plane at 5.7 d. Curve fits are estimated by eye and shown in dashed lines. α represents the exponent of a power law of sin 2θ. Angular momentum is normalized in radius with a factor of r−1/2. Density, pressure, and angular momentum are averaged over φ and 10 < r < 100 using equation (11). Pressure and angular momentum are weighted by mass. The stream is ignored using the entropy condition. Density, pressure, and angular momentum squared are multiplied by 0.158, 316, and 10−8 respectively so that all three quantities are roughly the same order for comparison purposes.

Figure 7.

Mass density, pressure, and angular momentum squared in our aligned TDET0 model plotted with respect to θ, their fits to a power law of sin ²θ, and an inset plot of rest mass density in the xz-plane at 5.7 d. Curve fits are estimated by eye and shown in dashed lines. α represents the exponent of a power law of sin ²θ. Angular momentum is normalized in radius with a factor of r^−1/2. Density, pressure, and angular momentum are averaged over φ and 10 < r < 100 using equation (11). Pressure and angular momentum are weighted by mass. The stream is ignored using the entropy condition. Density, pressure, and angular momentum squared are multiplied by 0.158, 316, and 10⁻⁸ respectively so that all three quantities are roughly the same order for comparison purposes.

Open in new tab Download slide

The vertical height of the accretion disc above the midplane is proportional to the distance from the centre of the BH, where the constant of proportionality is known as the scale height. We compute the scale height in two ways (Fig. 8).

$$\begin{eqnarray*} \frac{h_{\rm{sqrt}}}{r}=\sqrt{\frac{\int { (\theta -\pi /2)^2\rho u^t {\rm d}V}}{\int {\rho u^t {\rm d}V}}}, \end{eqnarray*}$$

(13)

$$\begin{eqnarray*} \frac{h_{\rm{abs}}}{r}=\frac{\int { |\theta -\pi /2| \rho u^t {\rm d}V}}{\int {\rho u^t {\rm d}V}}, \end{eqnarray*}$$

(14)

where

$$\begin{eqnarray*} {\rm d}V=\sqrt{-g}{\rm d}r{\rm d}\theta {\rm d}\phi \end{eqnarray*}$$

(15)

Figure 8.

The scale height of the disc plotted with respect to time in our aligned TDET0 model. Scale height is calculated by averaging the mass-weighted angle from the equatorial plane over |θ − π/2| < 0.3 and φ using two methods described by equations (13) and (14). The stream is ignored using the entropy condition. The generation of thermal energy due to stream–disc interactions and self-intersection shocks causes the gas in the disc to expand over time, increasing the scale height.

Open in new tab Download slide

Both methods show that h/r increases over time, indicating that the disc ‘puffs up’ from the midplane. Although it is possible that the vertical expansion of the disc is artificially slowed by the pressure and density floors, this effect should not significantly impact this general trend. The increase in the angular extent of the material is due to excess thermal energy generated in the disc by the dissipation of orbital energy. The scale height reaches a plateau around the time that the self-intersections stop (3.68 d), suggesting that the violent self-intersections play a crucial role in the early heating of the disc. We discuss the mechanisms of energy dissipation further in Section 4.1.

3.2 Tilted disc formation and evolution

The majority of TDE disc formation simulations use either Newtonian gravity or a general relativistic treatment (exact or approximate) of a non-spinning Schwarzschild BH. However, tidally disrupted stars approach the BH from a quasi-isotropic distribution of inclinations, highlighting the importance of more general disc formation simulations that account not just for BH spin, but also for spin-orbit misalignment. Various effects unique to tilted accretion discs, such as global precession, Bardeen–Petterson alignment, and disc tearing (Nixon & King 2012; Hawley & Krolkik 2019; Liska et al. 2021) may all manifest themselves in TDE accretion discs. Although previous studies have considered tilted TDEs analytically (Stone & Loeb 2012; Zanazzi & Lai 2019), only two numerical efforts have, to date, simulated the formation of an accretion flow following the disruption of a star on a misaligned orbit: the early work of Hayasaki et al. (2016), and the more recent simulations of Liptai et al. (2019). In both cases, the authors find that for adiabatic gas equations of state, it is challenging for nodal precession to cause significant delays in self-intersection. However, both of these simulations employed unrealistically eccentric stellar trajectories for computational convenience; the work presented in this section is the first numerical simulation of tilted TDEs with realistic astrophysical parameters.

The movie and 3D rendering (Section 6) of the tilted TDE shows that, unlike in the aligned TDE, the returning stream is never completely interrupted by self-intersections. The misalignment between the orbital plane of the stream and the rotational plane of the BH leads to strong nodal precession upon pericentre passage. The outgoing stream exits the BH in a separate plane from the incoming stream, so when the two streams collide, they are misaligned. Fig. 9 shows that this misalignment launches material from both streams out of their original planes.

A 3D contour of density (ρ = 10−8) visualized at 3.7 d in our tilted TDET30 model. We show the views from above and below the BH orbital plane in the left-hand and right-hand panels, respectively. The outgoing and incoming streams are misaligned at the self-intersection point due to the nodal precession at the pericentre passage. As a result, material is ejected out of the orbital plane of the star. The orbital plane of the BH is shown for reference. To see this effect in the full context of disc formation, see the 3D renderings linked in 6.

Figure 9.

A 3D contour of density (ρ = 10⁻⁸) visualized at 3.7 d in our tilted TDET30 model. We show the views from above and below the BH orbital plane in the left-hand and right-hand panels, respectively. The outgoing and incoming streams are misaligned at the self-intersection point due to the nodal precession at the pericentre passage. As a result, material is ejected out of the orbital plane of the star. The orbital plane of the BH is shown for reference. To see this effect in the full context of disc formation, see the 3D renderings linked in 6.

Open in new tab Download slide

Fig. 10 shows that the radial profiles of the tilted disc follow similar power-law relationships to those of the aligned disc, with the density and pressure falling off slightly faster in the aligned disc. As a result, the thermal pressure gradient forces in the tilted disc are larger than in the aligned disc, so the velocity distribution is more sub-Keplerian.

From a comparison of the polar density profiles of the aligned and tilted TDE (Figs 7 and 11), the accretion disc appears significantly thicker in the tilted simulation. However, there are two caveats to this interpretation of the data. First, the launching of stream material into different orbital planes as shown in Fig. 9 creates pockets of high density material at large polar angles which artificially increases our estimates of disc thickness. Secondly, the disc material may lie in multiple orbital planes because debris which falls back at early times has more time to undergo nodal precession than debris which falls back at late times. Any dependence of the inclination angle on r or ϕ would artificially increase our estimate of disc thickness. Additionally, the angular momenta of the gas may not have enough time to homogenize because the simulation was not run for multiple viscous times of the disc.

Figure 10.

Analogous to Fig. 6, but for the tilted TDE simulation, model TDET30. The data are tilted using our tilting algorithm (Appendix C) such that the star’s orbital plane coincides with θ′ = π/2. Coordinates in the tilted frame are denoted with a prime. The inset plot shows the rest mass density in the star’s orbital plane at 4.0 d. See Table B2 for more details about the power-law fits (dashed lines). Mass density and pressure are multiplied by 5 × 10⁵ so that all three variables are roughly the same order for comparison purposes.

Open in new tab Download slide

Figure 11.

Analogous to Fig. 7, but for the tilted TDE simulation, model TDET30. The data are tilted using our tilting algorithm (Appendix C) such that the star’s orbital plane coincides with θ′ = π/2. Coordinates in the tilted frame are denoted with a prime. Density and pressure are multiplied by 10^7.3 and 10^10.7 respectively so that all three quantities are roughly the same order for comparison purposes. The tilt angle varies over radius and time (Fig. 23), so the averages may not accurately reflect the thickness of the disc.

Open in new tab Download slide

4 DISCUSSION

4.1 Energy dissipation

The largest uncertainty in TDE evolution concerns the rate, location, and physical mechanisms that dissipate the orbital energy of dynamically cold debris streams. Past analytic models and numerical simulations examining early stages of a TDE generally focus on shock dissipation at different locations. The most important shock loci seen or proposed in past work are (i) compression shocks produced by the vertical collapse of the returning debris stream, located at the vertical caustics near pericentre (Guillochon et al. 2014; Shiokawa et al. 2015), (ii) self-intersection shocks produced at larger radii where an outgoing debris stream impacts an incoming one (Rees 1988; Hayasaki et al. 2013; Dai et al. 2015; Hayasaki et al. 2016), and (iii) ‘secondary shocks’ seen in the simulations of Bonnerot & Lu (2020), Bonnerot, Lu & Hopkins (2021) after the formation of an extended accretion flow. Each of these categories of shock have been seen to be the dominant energy dissipation mechanism at some times in some past numerical simulations of tidal disruption; however, the relative importance of these shocks is strongly affected by system parameters (e.g. SMBH mass, stellar eccentricity). The lack of published first principles TDE simulations with astrophysically realistic parameters renders the relative importance of these shocks unclear in typical TDEs.

We analyse the energy dissipation of the system at both early and late times by tracking entropy and the ratio of the thermal to total energy flux along a streamline. The thermal and mass energy fluxes are given by

$$\begin{eqnarray*} \Phi _{\rm{thermal}}=-\sqrt{-g} (u_g + p) \eta u_t=-\sqrt{-g} u_g \gamma \eta u_t \end{eqnarray*}$$

(16)

$$\begin{eqnarray*} \Phi _{\rm{mass}}=-\sqrt{-g} \rho c^2 \eta u_t , \end{eqnarray*}$$

(17)

where |$\eta =\sqrt{v_{\rm i} v^{\rm i}}=\sqrt{g_{\rm ij}v^{\rm i}v^{\rm j}}$| is the magnitude of the 3-velocity (we are adopting the convention where Latin indices range from 1–3). We define the quantity

$$\begin{eqnarray*} \psi =\frac{\Phi _{\rm{thermal}}}{\Phi _{\rm{thermal}} + \Phi _{\rm{mass}}} \end{eqnarray*}$$

(18)

to characterize the generation of thermal energy in the disc (solid red lines in Figs 12 and 13). When thermal energy flux dominates, ψ approaches unity and when the mass energy flux dominates, ψ approaches zero. Shocks convert orbital energy into thermal energy, so ψ increases across the self-intersection shocks. We compare the entropy of the stream post-pericentre to the average entropy of the disc, with the latter defined as

$$\begin{eqnarray*} {\kappa} _{\rm{disc}}=\frac{\int {\kappa \rho u^t {\rm d}V}}{\int {\rho u^t {\rm d}V}}. \end{eqnarray*}$$

(19)

Here κ is defined by equation (7) and the region of integration is defined using a stricter version of the entropy condition (κ > 100) to ensure that none of the high-density, low-entropy stream material contributes to the average.

Figure 12.

At early times in our aligned TDET0 simulation, the stream heats up at the pericentre and the self-intersection. However, the majority of the entropy generation occurs at the self-intersection, suggesting that the pericentre heating is nearly adiabatic. To see this, we plot the relative amount of heating ψ (equation 18) and a proxy for entropy κ ∝ e^entropy against the distance along a streamline s_strm at an early time, 1.34 d. The streamline, depicted in the inset plot, is integrated from the velocity field using a second-order Runge-Kutta method. The location of the pericentre and self-intersections along the streamline are depicted by vertical dashed and dot-dashed lines, respectively. The thin grey line shows the distance between the streamline and the BH. The average entropy of the disc is shown by the blue horizontal dotted line. The thin and thick red lines show an approximation for the relative amount of heating ψ for a thin (h/r = 0.082) and thick (h/r = 1) disc, respectively (equation 21). Note that the returning gas will not necessarily follow the path of the streamline, so s_strm is not a perfect proxy for time. For instance, the entropy sometimes decreases mildly as s_strm increases.

Open in new tab Download slide

Figure 13.

Analogous to Fig. 12, but at a late time of 5.66 d. At late times in our aligned TDET0 simulation, the bulk of the heating and entropy generation occurs at the pericentre radius, suggesting that the pericentre is the most significant source of energy dissipation. The thin and thick red lines show an approximation for the relative amount of heating ψ for a thin (h/r = 0.132) and thick (h/r = 1) disc respectively (equation 21).

Open in new tab Download slide

Figs 12 and 13 show the early- and late-time dissipation profiles along a streamline, respectively. Although heating and entropy generation occur on similar levels at both times, the dissipation mechanisms are distinct. At early times, there is comparable heating as the stream passes through pericentre and as the outgoing stream reaches the intersection point. However, significant entropy generation occurs only at the self-intersection, suggesting that the heating at pericentre is nearly adiabatic while the heating at self-intersection is irreversible and shock-induced. The nearly adiabatic heating at pericentre implies that the nozzle shock is inefficient at dissipating orbital energy.

This finding differs from other numerical work which finds that the nozzle shock is associated with a large jump in entropy (Guillochon et al. 2014; Bonnerot & Lu 2021). The strong nozzle shock seen in these works requires a highly supersonic vertical collapse of the debris stream. However, we find that the vertical collapse of the debris stream at early times is subsonic, possibly due to the relatively short return time of the debris or additional heating produced in the deeply penetrating β = 7 encounter we consider. The average mass-weighted vertical Mach number in the stream at early times (distinguished by the entropy condition) is M_z ≈ 0.4. The vertical Mach number is estimated as M_z ≈ v^z/c_s, where v^z is the magnitude of the vertical velocity and c_s is the sound speed given by

$$\begin{eqnarray*} c_s^2 \approx \left(\frac{\partial p}{\partial \rho } \right)_{\mathcal {S}} = \frac{\gamma (\gamma - 1) u_g}{\rho } . \end{eqnarray*}$$

(20)

At late times, the bulk of the heating and the entropy generation occurs at the pericentre. After the pericentre passage, entropy increases to more than three quarters the entropy of the disc, with the remainder of the entropy generation occurring between the pericentre and self-intersection radii.

The greater value of ψ along the streamline in Fig. 13 compared to Fig. 12 indicates that more heating takes place along the streamline at late times, possibly due to the stream–disc interactions discussed in Section 4.1. However, this does not imply that more heating occurs at late times in the disc as a whole. In particular, we find that the average value of ψ within the disc is greater at time of Fig. 12 than at the time of Fig. 13, suggesting that more relative heating occurs within the disc at early times.

As discussed in Section 3.1, the accretion disc becomes thicker and more massive over time. Therefore, as time progresses, the disc absorbs a greater proportion of the momentum of the outgoing stream, reducing the impact of the self-intersections. When the rotational momentum of the disc surpasses the momentum of the outgoing stream, the outgoing stream disintegrates and the self-intersections stop all together.

In Fig. 14, we plot the rotational mass flux of the disc and the mass fallback rate. The mass fallback rate is a good approximation for the mass flow rate in the outgoing stream; particularly at early times when the incoming and outgoing stream mass flow rates are most similar. At early times, the ratio of the rotational mass flux to the mass fallback rate is less than unity (Fig. 14, top panel). This reflects that the accretion disc mass is small relative to that of the outgoing stream. The ratio of rotational mass flux to mass fallback rate only exceeds unity shortly after the first self-intersection event at 1.5 d. With each self-intersection, the disc grows more substantive until the disc completely intercepts the outgoing stream and self-intersections can no longer occur around 3.7 d in our simulation.

Figure 14.

The mass fallback rate and the rotational mass flux in the disc in our aligned TDE simulation. The mass fallback rate is the mass flux in the stream (distinguished by the entropy condition) through the tidal sphere. The rotational mass flux is estimated by azimuthal mass flux of the disc through the surface given by φ = 1.14π, 10 < r < 500, and π/2 − h/r < θ < π/2+h/r. The ratio of rotational mass flux to fallback rate is shown in the top panel. As the disc becomes denser and more massive over the course of the simulation, the rotational mass flux of the disc surpasses the mass fallback rate. Around 3.7 d, the disc completely intercepts the momentum of the outgoing stream and self-intersections can no longer occur.

Open in new tab Download slide

At early times, the incoming stream heats up before the pericentre passage at the self-intersection point. This is reflected in Fig. 12 by the dramatic increase in ψ in the incoming stream at the intersection point. At the late times, this intersection is too weak to appreciably heat the incoming stream. Instead, the heating before the pericentre passage is caused by the collision between the incoming stream and the dense inner accretion disc. This is reflected in Fig. 13 by the increase in ψ and entropy after the self-intersection point and before the pericentre passage. Together, these results suggest that the self-intersections play a larger role in the energy dissipation at early times in the TDE evolution.

To provide more context for the meaning of the quantity ψ, we compute an approximation for ψ in a thick (h/r = 1) and thin disc, where the thin disc approximation assumes a scale height equal to the scale height of the disc (0.069 at 1.4 d and 0.132 at 5.7 d).

$$\begin{eqnarray*} \psi \approx \frac{u_g \gamma }{\rho } \approx c_s^2 \approx \left(\frac{h}{r}\right)^2 v_k^2 \approx \left(\frac{h}{r}\right)^2 \frac{1}{r}, \end{eqnarray*}$$

(21)

where c_s is the sound speed in the disc. This approximation appears as a dotted line in Figs 12 and 13. Note that equation (21) is only a good approximation for small values of ψ. Contrary to the approximation, we find that ψ does not drop off with radius, especially at late times, due to the heating that occurs as the stream disintegrates into the disc.

The inclusion of a more realistic equation of state within the debris stream is expected to yield an even higher rate of dissipation due to hydrodynamical shocks (see Guillochon et al. 2014, section 3.3, paragraph 6). Additionally, the inclusion of magnetic fields is expected to yield extra dissipation through the action of magnetorotational instability (MRI; Balbus & Hawley 1991) when the inner and outer stream develop a strong shear in velocity near pericentre.

4.2 Circularization

In the standard TDE picture, the accretion disc circularizes efficiently as shocks dissipate orbital energy, resulting in a nearly axis-symmetric accretion disc with low eccentricity. However, not all TDE discs circularize completely (Piran et al. 2015). Cao et al. (2018) find that the optical emission lines of TDE ASASSN-14li are best modelled by an accretion disc with eccentricity e = 0.97. Furthermore, recent analytical work on TDEs has derived eccentric disc solutions which can produce radiation consistent with the X-ray and optical luminosities of many TDE candidates (Zanazzi & Ogilvie 2020).

In our simulation, the accretion disc tends towards circularization but never fully circularizes according to the criterion suggested by Bonnerot, Rossi & Lodato (2017): an average eccentricity lower than 1/3. Instead, our disc reaches an average eccentricity of 0.88 at late times, where eccentricity is given by

$$\begin{eqnarray*} e=\sqrt{1 + \frac{2 \varepsilon l^2}{G^2 M_{\rm{BH}}^2}}. \end{eqnarray*}$$

(22)

Here ε = −(u_t+1) is the total orbital energy, l = u_ϕ is the specific angular momentum, and the average is taken from r = 10R_g to r = 400R_g. It is important to note that this formula only provides an upper bound on the geometric eccentricity of particle trajectories because it assumes that a given fluid element is acted on only by a Newtonian gravitational force. For example, in the presence of internal pressure support, equation (22) gives values greater than zero for circular orbits. Despite this issue, test particle eccentricity is still a useful metric for the extent of circularization.

Due to the short duration of our simulation, we cannot confirm whether more complete circularization will occur at times after the end of our simulation. One factor that may inhibit circularization is the injection of high-eccentricity material into the disc by the returning debris stream. However, this effect will become negligible when the mass fallback declines to the point where energy and angular momentum input are negligible in analogy to the late-time behaviour of the disc mass in Cannizzo et al. (1990).

Fig. 15 shows a time-averaged streamline plot of the inner part of the disc coloured by eccentricity. The area immediately surrounding the stream contains high-eccentricity material, indicating that the stream continuously transfers its energy and angular momentum into the disc. Away from the stream, disc material orbits at more moderate eccentricities. Looking at snapshots throughout the duration of the simulation, we see a similar distribution of eccentricities in the disc (Fig. 16). Fig. 15 also shows that the eccentricity changes along each streamline, suggesting a continuous transport of energy and angular momentum within the disc itself.

Figure 15.

A time-averaged velocity streamline plot of the inner parts of the disc in the equatorial plane coloured by time-averaged eccentricity. Velocity and eccentricity are time-averaged over the simulation’s entire duration, averaged over |θ − π/2| < 0.5, and weighted by rest mass density. The stream is ignored using the entropy condition. The streamlines are integrated using a second-order Runge-Kutta method. Ellipses of various eccentricities are overlaid to provide context for the eccentricity data.

Open in new tab Download slide

Figure 16.

Snapshots of mass-weighted eccentricity in the equatorial plane at various times (see legend) averaged over θ = π/2 ± 0.5. The stream is ignored using the entropy condition. However, there still is an abundance of high-eccentricity material around the stream, suggesting that the stream constantly transfers its energy and angular momentum into the disc.

Open in new tab Download slide

Eccentricity is not evenly distributed across the different radii in the disc (Fig. 17). In particular, the inner parts of the disc are more circularized than the outer parts. This indicates that circularization is more efficient at smaller radii, possibly because the velocity shear between neighbouring radii is greater. Because eccentricity affects mean φ-velocity (Section 4.4), the uneven eccentricity distribution may contribute to the drop in φ-velocity near the outer edge of the disc in Fig. 6.

Figure 17.

Radial profile of eccentricity at various times (see legend; compare to Fig. 16). Eccentricity is mass-weighted and averaged over spherical shells. The stream is ignored using the entropy condition. Unbound material is ignored using the Bernoulli parameter. Only radii from 10R_g/c to 500R_g/c are shown. Eccentricity is unevenly distributed throughout the disc. In particular, we see greater circularization at smaller radii in the disc.

Open in new tab Download slide

In Fig. 17, there is a dip in the eccentricity at the self-intersection radius at 1.6 d. This may be a result of the second major self-intersection event, which occurs at 1.47 d, interrupting the incoming stream. This is visible in Fig. 16, where we can see that the effect of the stream on the disc eccentricity is much weaker than at any other time slice shown.

4.3 Accretion and outflow

Some TDE accretion models predict a period of super-Eddington accretion, the magnitude and duration of which depend on the fallback rate (Coughlin & Begelman 2014; Wu, Coughlin & Nixon 2018). This prediction is supported by observations; for instance, TDE Swift J1644+57 exhibits a super-Eddington luminosity (Burrows et al. 2011; Zauderer et al. 2011). We calculate the theoretical Eddington accretion rate for our simulation below. Assuming that the accreting material is mostly ionized hydrogen gas, the Eddington luminosity is

$$\begin{eqnarray*} L_{\rm{Edd}}=\frac{4 \pi G M_{\rm{BH}} c m_p}{\sigma _T}, \end{eqnarray*}$$

(23)

where m_p is the mass of the proton and σ_T is the Thomson cross-section. From the Eddington luminosity, the Eddington accretion rate is

$$\begin{eqnarray*} \dot{M}_{\rm{Edd}}=\frac{L_{\rm{Edd}}}{\epsilon c^2}, \end{eqnarray*}$$

(24)

where ϵ is the gravitational potential energy that is radiated as a fraction of the rest-mass energy. Combining the above two expressions yields

$$\begin{eqnarray*} \dot{M}_{\rm{Edd}}=\frac{4 \pi G M_{\rm{BH}} m_p}{\sigma _T \epsilon c}. \end{eqnarray*}$$

(25)

If we assume that ϵ = 0.1, then we find |$\dot{M}_{\rm{Edd}} \simeq 0.022\, {\rm M}_\odot \rm{yr}^{-1}$|⁠. Fig. 18 shows that the mass accretion rate at the BH reaches up to twice the Eddington limit, and the mass fallback rate reaches up to eight times the Eddington limit. This confirms that the TDE in our simulation exhibits the predicted period of super-Eddington accretion. We noted in Section 2.1 that the maximum fallback rates in our simulation are about an order of magnitude lower than the peak fallback rate due to the short duration of our simulation. Therefore, at its peak the accretion rate would be even more super-Eddington than seen in Fig. 18.

Figure 18.

Mass fallback rate, mass accretion rate at the event horizon, mass outflow rate, and the accretion efficiency plotted versus time in our aligned TDET0 simulation. The accretion efficiency settles to around 10 to 20 per cent at t ≳ 4 d. Positive mass fluxes are directed towards the BH. The accretion efficiency is calculated as a ratio of the mass flux at the event horizon to the mass fallback rate. Mass fallback rate is computed within the stream (distinguished by the entropy condition) through the tidal sphere. Mass outflow rate is computed as the unbound mass flux through the tidal radius. Bound matter is ignored using the Bernoulli parameter. All three mass fluxes increase roughly linearly with time, suggesting that the disc mass increases quadratically.

Open in new tab Download slide

Fig. 19 shows that bound and unbound material accrete on to the BH at increasingly higher rates as the simulation evolves. However, outside of the innermost radii in the disc, more material flows away from the BH than towards the BH. This outward mass flux drives the radial expansion of the disc. The figure also shows that the radial mass fluxes of bound and unbound material are similar throughout the simulation, indicating that a significant fraction of the initially bound material gets unbound. We can loosely estimate this fraction by comparing the outflow rate to the mass flux fallback in Fig. 18. Due to the short duration of our simulation, we compute outflow rates as the mass flux of unbound disc material through the tidal sphere. In a longer simulation, we could compute a more precise outflow rate by computing the mass flux at larger radii and at later times. Fig. 18 shows that the outflow rate is approximately 1/3 of the fallback rate, suggesting that this same fraction of infalling material becomes unbound. However, it is possible that the outflows may be artificially suppressed by the density and pressure floors discussed in Section 3.1, so this may be an underestimate for the fraction of material which becomes unbound in a physical TDE.

$Radial profiles of the mass flux in the disc in units of the Eddington accretion rate. We distinguish between the bound and unbound material at three different times, spread evenly across the duration of the simulation. Positive mass fluxes are directed towards the BH. Bound and unbound material is distinguished using the Bernoulli parameter. The mass flux of bound material is shown in solid lines, and the mass flux of unbound material is shown in dotted lines. The stream is ignored using the entropy condition. The horizontal dotted lines show $\dot{M}=0$. Bound and unbound material accrete on to the BH at increasingly higher rates as the simulation progresses.$

Figure 19.

Radial profiles of the mass flux in the disc in units of the Eddington accretion rate. We distinguish between the bound and unbound material at three different times, spread evenly across the duration of the simulation. Positive mass fluxes are directed towards the BH. Bound and unbound material is distinguished using the Bernoulli parameter. The mass flux of bound material is shown in solid lines, and the mass flux of unbound material is shown in dotted lines. The stream is ignored using the entropy condition. The horizontal dotted lines show |$\dot{M}=0$|⁠. Bound and unbound material accrete on to the BH at increasingly higher rates as the simulation progresses.

Open in new tab Download slide

The rate at which mass is added to the disc is given approximately by subtracting the mass flux outflow and mass accretion rate from the mass flux fallback. Because all three mass fluxes increase linearly with time (Fig. 18), the mass of the disc increases quadratically over the course of our simulation. We verify this by fitting the disc mass to a quadratic time-dependence using the least-squares method. We find that the disc mass in units of grams is given approximately by |$M_{\rm disc}/{\rm M}_\odot = 4.85 \times 10^{-8} t_{\rm days}^2$| with a coefficient of determination R² = 0.994.

We also quantify the accretion efficiency (Fig. 18) using the ratio of the mass fallback rate to the mass accretion rate at the event horizon, where the mass flux is computed as

$$\begin{eqnarray*} \dot{M}=-\int {\rho u^r {\rm d}A_{\theta \phi }}. \end{eqnarray*}$$

(26)

The accretion efficiency exhibits some periodicity due to the periodic self-intersection of the stream (Section 4.3). However, after the initial spike at two days into the disruption, it settles into a range of 10 to 20 per cent, suggesting that our TDE is reasonably efficient at getting the gas from the debris stream to the black hole.

If all accreting material were fully circularized, the accretion efficiency would place a lower limit on the extent of circularization. Then, if no circularization occurred, then the mass accretion rate at the event horizon would vanish, and if the disc circularized completely, then part of the material that falls back would eventually accrete (the rest would fly out as an outflow), giving the lower limit to extent of circularization.

However, we find that a significant fraction of the stellar debris accretes with moderate or high eccentricities, as shown in the left-hand panel of Fig. 20. Here, we plot a time-averaged histogram of the eccentricity of material inside the innermost stable circular orbit, R_ISCO = 2.04R_g (for prograde orbits in the equatorial plane and our black hole spin value of a = 0.9375; Bardeen, Press & Teukolsky 1972).

$The left-hand panel shows a time-averaged histogram of the eccentricity of material inside RISCO, weighted by rest-mass. Each bin is normalized by total accreting mass and bin width. The mean eccentricity, 〈e〉 ≈ 0.73, is shown by the vertical dashed line. The right-hand panel shows a time averaged colour map of the eccentricity of material inside RISCO, weighted by rest mass. The overlaid grayscale contours outline regions containing 25 and 50 per cent of the total mass flux. Floor material is ignored using an angular momentum condition (|uϕ| > 10−5). The left-hand panel indicates that a significant fraction of material accretes with moderate eccentricities, allowing for a high accretion efficiency without complete circularization. However, from the right-hand panel, we see that highly eccentric material preferentially accretes at high latitudes, except for the material at azimuthal angles just beyond that of the pericentre (see the text for discussion). Most of the low-latitude (disc) accretion occurs at moderate eccentricity.$

Figure 20.

The left-hand panel shows a time-averaged histogram of the eccentricity of material inside R_ISCO, weighted by rest-mass. Each bin is normalized by total accreting mass and bin width. The mean eccentricity, 〈e〉 ≈ 0.73, is shown by the vertical dashed line. The right-hand panel shows a time averaged colour map of the eccentricity of material inside R_ISCO, weighted by rest mass. The overlaid grayscale contours outline regions containing 25 and 50 per cent of the total mass flux. Floor material is ignored using an angular momentum condition (|u_ϕ| > 10⁻⁵). The left-hand panel indicates that a significant fraction of material accretes with moderate eccentricities, allowing for a high accretion efficiency without complete circularization. However, from the right-hand panel, we see that highly eccentric material preferentially accretes at high latitudes, except for the material at azimuthal angles just beyond that of the pericentre (see the text for discussion). Most of the low-latitude (disc) accretion occurs at moderate eccentricity.

Open in new tab Download slide

Eccentric accretion has been found in several earlier works (Shiokawa et al. 2015; Sądowski et al. 2016; Bonnerot & Lu 2020) and may explain the low luminosity and temperature observed in TDE candidates. A more extreme model was proposed by Svirski et al. (2017) in which magnetic stresses transport the angular momentum away from the black hole, driving eccentric accretion. However, this model does not apply to our simulation, which does not include magnetic fields.

There are two ways that material can accrete with moderate or high eccentricities. Gas in the disc may lose angular momentum through, e.g. turbulent viscosity, shocks, or spiral waves, until it plunges into the black hole. This process is analogous to accretion in a quasi-circular accretion disc. Alternatively, gas at high latitudes may get torqued and free fall directly into the black hole, an effect which has been observed in previous works (see fig. 12 of Sądowski et al. 2016). We call the first type of accretion eccentric disc accretion and the second type of accretion ballistic polar accretion.

We can differentiate between these two types of accretion by the latitude at which matter enters the region r < R_ISCO. We find that highly eccentric material (e > 0.7) preferentially accretes at high latitudes, suggesting that the dominant method of highly eccentric accretion is ballistic polar accretion (Fig. 20). Material accreting at low latitudes generally has more moderate eccentricities (0.4 < e < 0.7), and is thus likely driven by eccentric disc accretion. However, there is a patch of high eccentricity material at low latitudes in the region 2π/3 < ϕ < π, just past the pericentre at an azimuthal angle of ϕ = π/2. Accretion in this region may be due to coherent chunks of the incoming stream that undergo turbulent exchange angular momentum with the disc and accrete directly on to the black hole (see also Fig. 16). The innermost (25 per cent) mass flux contour in Fig. 20 indicates that the majority of accretion occurs in this region.

As we describe in Section 3.1, the debris stream in our aligned TDE simulation collides with itself in five violent self-intersection events that occur approximately 12 h apart and last for roughly 2000R_g/c, or 2.74 h. We apply these results to TDE Swift J1644+57, which exhibits quasi-periodic flaring during the first few days of its initial evolution. Other authors have proposed that this flaring is due to a precessing jet (Stone & Loeb 2012; Tchekhovskoy et al. 2013). However, our simulations show that even without precession, the flaring due to violent self-intersections can explain both the number of flares and their time-scale. Swift J1644+57 is a |$10^5\!-\!10^6\, {\rm M}_{\odot }$| BH, so 1 d corresponds to a time-scale of 5000–50000 R_g/c, similar to the time-scale of the self-intersections in our simulation.

It is unlikely that this flaring is a direct consequence of self-intersection events because the material at the self-intersection point is too optically thick to produce X-rays without adiabatic cooling (Jiang, Guillochon & Loeb 2016). Instead, we propose that the periodicity of the self-intersections leads to a periodicity in the accretion that feeds the jets, an effect that we see in our aligned TDE simulation.

As we discuss in Section 4.2, we can normalize the mass accretion rate at the event horizon by the mass fallback rate for the aligned and tilted simulations. We see quasi-periodic behaviour only in the aligned case (Fig. 18) where violent, periodic self-intersections occur. This behaviour does not perfectly correlate with the major self-intersection events in the simulation, which may be due to the similar time-scale of the self-intersections (∼12 h apart lasting ∼3 each) and the fallback time from the self-intersection point (∼4 h). However, the large fluctuations in accretion rate stop after the last major self-intersection event at 3.7 d.

Sądowski et al. (2016) found a marginally bound torus after self-intersection with some unbound material at high polar angles. They also found periodic behaviour due to the interactions of the outgoing stream with the incoming stream. However, this interaction was not as violent as in our simulation, which may be due to the differences in the orbital properties of the initial star.

4.4 Force balance

In Section 3.1, we show that the φ-velocities in the accretion disc are 76 per cent of the expected φ-velocities in a circularized Keplerian accretion disc. There are two possible explanations for our findings. First, the non-zero eccentricity of the disc decreases the average φ-velocity of the disc relative to a completely circularized disc. Secondly, the disc is internally supported by non-gravitational forces.

As we discuss in Section 4.2, the average eccentricity of the disc at late times is e = 0.88. To determine the effect of this eccentricity on the velocity distribution, we set up artificial velocity fields with constant eccentricities of 0 and 0.88 using a method described in Appendix A3. For each velocity field, we compute the radial profile of φ-velocity with equation (11), where the density weight is determined as a function of radius by the power-law relationship depicted in Fig. 6. On average, we find that the φ-velocities in the eccentric disc are 93.8 per cent of the φ-velocities in the circularized disc. Therefore, the φ-velocities in our accretion disc are at most |$76{{\ \rm per\ cent}} / 0.938 = 81{{\ \rm per\ cent}}$| of the φ-velocities in a Keplerian accretion disc at the same eccentricity. Therefore, the disc must also be externally supported by non-gravitational forces.

The only non-gravitational forces in our simulation are thermal pressure-gradient forces. Because pressure drops off as the distance from the BH increases (Fig. 6), these forces are directed away from the BH in the equatorial plane, reducing the centripetal force on the accretion disc. To analyse the force balance in the disc, we compute the ratio of the pressure-gradient force density, ∇p, to the centripetal force density required to maintain Keplerian orbits (Fig. 21), where the Keplerian centripetal force density is

$$\begin{eqnarray*} f_{\rm{kep}}=\rho \frac{v_K^2}{r}=\rho \frac{GM_{{\rm{BH}}}}{r^2}. \end{eqnarray*}$$

(27)

As the ratio of force densities increases, the matter in the disc forms stable orbits at increasingly sub-Keplerian velocities. We find that the gradient force is a substantial fraction of the Keplerian centripetal force at all radii in the disc allowing the disc material into maintain sub-Keplerian velocities. In the inner parts of the disc, the pressure-gradient force ranges from 25 to 40 per cent of the Keplerian centripetal force, which fully accounts for the sub-Keplerian φ-velocity distribution.

$A time-averaged radial profile of the pressure gradient force density normalized by the Keplerian force density. Within the disc, we see that pressure gradient forces are a significant fraction of the Keplerian force, accounting for the sub-Keplerian φ-velocity distribution that we find in Fig. 6. Pressure is averaged over spherical shells and mass-weighted. Time averages are over the simulation’s full duration. The stream is ignored using the entropy condition. We apply a smoothing function to the data to improve readability.$

Figure 21.

A time-averaged radial profile of the pressure gradient force density normalized by the Keplerian force density. Within the disc, we see that pressure gradient forces are a significant fraction of the Keplerian force, accounting for the sub-Keplerian φ-velocity distribution that we find in Fig. 6. Pressure is averaged over spherical shells and mass-weighted. Time averages are over the simulation’s full duration. The stream is ignored using the entropy condition. We apply a smoothing function to the data to improve readability.

Open in new tab Download slide

4.5 Comparison to disc models

4.5.1 Comparison to ZEBRA model

We compare the properties of the post-intersection accretion flow in our simulation with those predicted by the ZEro-BeRnoulli Accretion (ZEBRA) model proposed by Coughlin & Begelman (2014). The key assumptions of the ZEBRA model are that

the Bernoulli parameter b is zero everywhere,
the potential is Newtonian,
the magnetic energy density is not sufficient to destabilize the disc with respect to the Høiland criteria.

Coughlin & Begelman (2014) show that assumption (i) ensures that the disc is gyrentropic; that is, surfaces of constant entropy, angular momentum, and Bernoulli parameter coincide. From assumption (ii), Newton’s law, and the Bernoulli equation, Blandford & Begelman (2004) derive self-similar solutions for gyrentropic discs with an arbitrary Bernoulli parameter. The ZEBRA solutions are a special case when the Bernoulli parameter is zero everywhere.

Assumption (iii) trivially holds in our simulations due to the absence of magnetic fields. We show that assumption (i) holds in Fig. 2, which depicts a histogram of the Bernoulli parameter weighted by mass in the initial conditions. We calculate that the average mass-weighted Bernoulli parameter in the initial conditions is 3.6 × 10 − 5. The Bernoulli parameter in the disc is larger than this initial parameter because only the most bound debris reaches the black hole over the duration of our simulation. However, even at late times, the Bernoulli parameter is smaller than the gravitational binding energy within the disc. At 5.7 d into the simulation, the average mass-weighted Bernoulli parameter is approximately 1.5 × 10⁻³, or 38 per cent of the binding energy in the disc at a characteristic radius r₀ = 259R_g determined by the mass-weighted mean radial coordinate of the disc.

Assumption (ii) becomes less accurate as the radial coordinate approaches the gravitational radius. However, we find that the power-law relationships predicted by the ZEBRA model extend nearly to the inner boundary of the disc as shown in Fig. 6. The ZEBRA model self-similar solutions are

$$\begin{eqnarray*} \rho (r, \theta)=\rho _0\left(\frac{r}{r_0} \right)^{ - q}(\sin ^2\theta)^{\alpha }, \end{eqnarray*}$$

(28)

$$\begin{eqnarray*} p(r, \theta)=\beta \frac{GM_{\rm{BH}}\rho _0}{r}\left(\frac{r}{r_0}\right)^{ - q}(\sin ^2\theta)^{\alpha }, \end{eqnarray*}$$

(29)

$$\begin{eqnarray*} l^2(r,\theta)=aGM_{\rm{BH}}r\sin ^2\theta . \end{eqnarray*}$$

(30)

These solutions describe the accretion flow density, pressure, and squared specific angular momentum, respectively, with the additional definitions

$$\begin{eqnarray*} \alpha =\frac{1 - q(\gamma - 1)}{\gamma - 1}, \end{eqnarray*}$$

(31)

$$\begin{eqnarray*} \beta =\frac{\gamma - 1}{1+ \gamma - q(\gamma - 1)}, \end{eqnarray*}$$

(32)

$$\begin{eqnarray*} a=2\frac{1 - q(\gamma - 1)}{1+ \gamma - q(\gamma - 1)}, \end{eqnarray*}$$

(33)

where r₀ is some characteristic radius in the disc and ρ₀ is the density at that radius in the midplane. Of particular importance to our analysis are the following relationships:

ρ ∝ r^−q,
p ∝ r^{−q − 1},
1/2 < q < 3,
l ∝ sin ²θ,
ρ ∝ (sin ²θ)^α,
p ∝ (sin ²θ)^α,
Equation (31).

We compare the ZEBRA model predictions to the radial and polar profiles of our simulated disc (Figs 6 and 7). (a) and (b) imply that density and pressure depend on r as a power law and (e) and (f) imply that density and pressure depend on θ as a power law of sin ²θ. These predicted dependencies provide a reasonable fit for our data within the boundaries of the disc (Tables B1 and B2). Our fitted power-law exponents for the radial profiles density and pressure differ by 1.22, which nearly matches the difference of 1.0 predicted by (a) and (b). In the model of Coughlin & Begelman (2014), the power-law index q can be constrained by the mass inflow rate and prescribed disc physics (e.g. the fact that angular momentum is efficiently transported in the disc), but in general it is expected to be on the order of ∼1 − 2 (see fig. 8 of Coughlin & Begelman 2014 and fig. 4 of Wu et al. 2018), which is exactly what we see in our simulation.

Table B1.

Open in new tab

Curve-fitting results for the radial profiles of mass density, pressure, and φ-velocity from 20–250 R_g shown in Fig. 6, including the power-law parameters (ax^b) and their relative standard deviation errors. The exponent for φ-velocity is fixed at -0.5.

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	4.13E-7	4.71 per cent	−1.10	1.26 per cent
Pressure	2.54E-7	7.37 per cent	−2.32	1.00 per cent
φ-velocity	0.759	1.64 per cent

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	4.13E-7	4.71 per cent	−1.10	1.26 per cent
Pressure	2.54E-7	7.37 per cent	−2.32	1.00 per cent
φ-velocity	0.759	1.64 per cent

Table B1.

Open in new tab

Curve-fitting results for the radial profiles of mass density, pressure, and φ-velocity from 20–250 R_g shown in Fig. 6, including the power-law parameters (ax^b) and their relative standard deviation errors. The exponent for φ-velocity is fixed at -0.5.

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	4.13E-7	4.71 per cent	−1.10	1.26 per cent
Pressure	2.54E-7	7.37 per cent	−2.32	1.00 per cent
φ-velocity	0.759	1.64 per cent

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	4.13E-7	4.71 per cent	−1.10	1.26 per cent
Pressure	2.54E-7	7.37 per cent	−2.32	1.00 per cent
φ-velocity	0.759	1.64 per cent

Table B2.

Open in new tab

Curve-fitting results for the tilted TDE radial profiles of mass density, pressure, and φ-velocity from 20–250 R_g shown in Fig. 10, including the power-law parameters (ax^b) and their relative standard deviation errors. The exponent for φ-velocity is fixed at -0.5.

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	5.19E-7	1.73 per cent	−1.34	0.397 per cent
Pressure	2.91E-7	6.75 per cent	−2.62	0.840 per cent
φ-velocity	0.593	2.17 per cent

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	5.19E-7	1.73 per cent	−1.34	0.397 per cent
Pressure	2.91E-7	6.75 per cent	−2.62	0.840 per cent
φ-velocity	0.593	2.17 per cent

Table B2.

Open in new tab

Curve-fitting results for the tilted TDE radial profiles of mass density, pressure, and φ-velocity from 20–250 R_g shown in Fig. 10, including the power-law parameters (ax^b) and their relative standard deviation errors. The exponent for φ-velocity is fixed at -0.5.

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	5.19E-7	1.73 per cent	−1.34	0.397 per cent
Pressure	2.91E-7	6.75 per cent	−2.62	0.840 per cent
φ-velocity	0.593	2.17 per cent

Variable	a	σ_a/μ_a	b	σ_b/μ_b
ρ	5.19E-7	1.73 per cent	−1.34	0.397 per cent
Pressure	2.91E-7	6.75 per cent	−2.62	0.840 per cent
φ-velocity	0.593	2.17 per cent

However, the alpha parameter does not match its predicted value from the ZEBRA model. The power-law exponent for density indicates that q ∼ 1. Therefore, α ∼ 0.5 by equation (31). Instead, we find values for α of unity and 12 from pressure and density, respectively. In addition, l² is proportional to (sin ²θ)^2.2 rather than sin ²θ. These discrepancies indicate that the disc must be thinner than predicted by the ZEBRA model.

The ZEBRA model predicts that the specific angular momentum of the disc must be at least 76 per cent of the Keplerian value with our assumption of a polytropic index of 5/3 (Coughlin & Begelman 2014). Coincidentally, we find that the φ-velocities in the disc are 76 per cent of the Keplerian values for circular orbits. As we discuss in Section 4.4, the non-zero eccentricity of our disc automatically decreases the φ-velocities in the disc relative to the Keplerian velocity. Adjusting for this effect, the φ-velocities in the disc are 81 per cent of the Keplerian values.

As we mention in Section 3.1, the internal energy density and mass density floors at 2.27 × 10⁻¹² may artificially decrease the radial and vertical extent of the disc by providing external pressure support. Therefore, our results at radii within the disc boundaries (≲ 500R_g) are more reliable than at larger distances. Without the floors, it is possible that the power-law curves for density and pressure in Figs 6 and 10 would continue past 500R_g. This additional pressure confinement may also be responsible for the flattening of the disc as compared to the ZEBRA model. Because the floors are non-rotating, the external pressure could decrease the angular momentum of material at the edges of the disc, possibly leading to artificially efficient accretion.

4.5.2 Bonnerot and Lu model

Recently, Bonnerot & Lu (2020) performed a TDE simulation with a realistic stellar trajectory and mass ratio. They found that self-intersections launch outflows. These outflows undergo extensive ‘secondary shocks’ that ultimately result in the formation of an accretion disc. This contrasts with our results, in which the formation and circularization of the accretion disc results primarily from stream–disc interactions near pericentre (Section 4.1). Even when violent self-intersections do not occur, as in model TDET30, an accretion disc still forms.

Bonnerot & Lu (2020) overcame the numerical challenges of simulating a TDE with a realistic stellar trajectory and mass ratio by using a non-spinning BH and incorporating the local simulation of Lu & Bonnerot (2019) into their initial conditions to describe the outflows produced by self-intersection shocks. However, that local simulation includes assumptions that maximize the impact of the self-intersection shocks.

Lu & Bonnerot (2020) perform their simulation in a special inertial frame, which they refer to as the simulation box (SB) frame, in which the incoming and outgoing streams collide head-on. The frame is related to the lab frame by a boost to the comoving frame of a local stationary observer at the self-intersection radius followed by a boost to a frame where the φ-velocity of the outgoing stream vanishes. This technique relies on three assumptions.

First, it requires that the incoming and outgoing streams have equal (φ-velocities, aspect ratios) or precisely opposite (radial velocities) properties at the point of self-intersection. However, significant pericentre dissipation, Lense–Thirring frame dragging, or hydrodynamic instabilities at the boundary of the stream and the disc could cause the outgoing stream to have a lower density and velocity and a different trajectory relative to the incoming stream, decreasing the violence of the self-intersections.

Secondly, the incoming and outgoing streams are only completely parallel in the SB frame at the intersection point. As the radial distance from the self-intersection point increases, the head-on gas trajectories used in Lu & Bonnerot (2019) diverge from the physical trajectories. Therefore, the approach is only accurate in the case that there are minimal interactions between pre- and post-intersection material and, e.g. the accretion disc. By ignoring these interactions, the head-on approach increases the relative importance of the self-intersections in the overall TDE evolution.

Thirdly, the Lu & Bonnerot (2019) simulation uses 2D cylindrical coordinates, implicitly assuming axisymmetry of the colliding streams. This 2D approach also cannot fully capture 3D fluid instabilities and turbulence inherent in the violent interaction.

Many of the novel effects observed by Bonnerot & Lu (2020) are tied to the strong outflows sourced at the self-intersection point: for instance, the formation of a retrograde accretion disc with respect to the star’s initial orbital angular momentum is due to the preferential loss of prograde debris in the self-intersection outflows. In our simulations, stream–disc interactions near pericentre rapidly become the primary locus of energy dissipation and efficiently suppress the return of coherent outgoing streams to the self-intersection site. This behaviour is difficult to reconcile with local mass injection schemes near the self-intersection radius.

One source of the discrepancy between our work and that of Bonnerot & Lu (2020) may be the high-β encounter analyzed here, which causes the outgoing stream to significantly expand after the pericentre passage due to differential precession (Section 3.1). For a more standard TDE with β ≈ 1, the approach of Bonnerot & Lu (2020) may become more accurate. It is also possible that increasing the resolution across the vertical extent of the stream in our simulation would result in stronger nozzle shocks. We leave it to future work to apply our methods to less deeply penetrating TDEs.

4.6 Analysis of the tilted TDE simulation, model TDET30

Fig. 22 shows that the accretion efficiency of the tilted TDE is only slightly less than the aligned TDE, hovering from 10 to 15 per cent. Just like the aligned scenario, all three mass fluxes increase roughly linearly in time, implying a quadratically increasing disc mass.

Figure 22.

Analogous to Fig. 18, but for the tilted TDE simulation, model TDET30. The accretion efficiency is similar to that of the aligned run.

Open in new tab Download slide

In Fig. 23, we compute the tilt and precession angle of the disc using the method of Fragile et al. (2007), Fragile & Anninos (2005), Nelson & Papaloizou (2000) which we describe below for convenience. The tilt angle |$\mathcal {T}$| is given by

$$\begin{eqnarray*} \mathcal {T} (r) = \arccos \bigg [ \frac{\boldsymbol {J}_{\rm{BH}} \boldsymbol\cdot \boldsymbol {J}_{\rm{disc}}(r)}{|\boldsymbol {J}_{\rm{BH}}||\boldsymbol {J}_{\rm{disc}}|} \bigg ] , \end{eqnarray*}$$

(34)

where |$\boldsymbol {J}_{\rm{BH}} = a M_{\rm BH} \hat{\boldsymbol{z}}$| is the angular momentum vector of the BH and |$\boldsymbol {J}_{\rm{disc}}(r)$| is the angular momentum vector of the disc in an asymptotically flat space. |$\boldsymbol {J}_{\rm{disc}}(r)$| is given component-wise by

$$\begin{eqnarray*} (J_{\rm{disc}})_\rho =\frac{\epsilon _{\mu \nu \sigma \rho } L^{\mu \nu } S^\rho }{2 \sqrt{-S^\alpha S_\alpha }}, \end{eqnarray*}$$

(35)

where

$$\begin{eqnarray*} L^{\mu \nu } = \int {\left(x^\mu T^{\nu 0} - x^\nu T^{\mu 0}\right) {\rm d}V} \end{eqnarray*}$$

(36)

$$\begin{eqnarray*} S^\rho =\int {T^{\rho 0} {\rm d} V} \end{eqnarray*}$$

(37)

T is the stress–energy tensor, and ϵ is the 4D Levi–Civita symbol. The unit vector |$\hat{y}$| points along the axis about which the initial conditions are initially tilted and |$\hat{z}$| points along the angular momentum axis of the BH. The precession angle is computed similarly using the definition

$$\begin{eqnarray*} \mathcal {P}(r) = \arccos \bigg [ \frac{\boldsymbol {J}_{\rm{BH}} \boldsymbol\times \boldsymbol {J}_{\rm{disc}}(r)}{|\boldsymbol {J}_{\rm{BH}} \boldsymbol\times \boldsymbol {J}_{\rm{disc}}|} \boldsymbol\cdot \hat{\boldsymbol{y}} \bigg ] \end{eqnarray*}$$

(38)

$The tilt angle $\mathcal {T}$ and precession angle $\mathcal {P}$ of the tilted disc over time in our titled TDE simulation TDET30. Panels (a) and (b) depict the tilt and precession angles calculated from the net angular momentum of the disc for material r < 500Rg. Panels (c) and (d) depict the tilt and precession angles calculated from the net angular momentum of the disc material at the radius indicated on the y-axis. Both tilt and precession angles are continuous within the range of radii in the disc outside of stream’s pericentre (10Rg ≲ r ≲ 200Rg at the beginning of the simulation and 10Rg ≲ r ≲ 400Rg at the end of the simulation), indicating that there is no disc-tearing. Below 10Rg, tilt and precession angles fluctuate rapidly due to stream–disc interactions at pericentre. We do not observe significant precession of the disc in the duration of the simulation. It is possible that disc precession at small radii is inhibited by the constant injection of stream material in the orbital plane of the star.$

Figure 23.

The tilt angle |$\mathcal {T}$| and precession angle |$\mathcal {P}$| of the tilted disc over time in our titled TDE simulation TDET30. Panels (a) and (b) depict the tilt and precession angles calculated from the net angular momentum of the disc for material r < 500R_g. Panels (c) and (d) depict the tilt and precession angles calculated from the net angular momentum of the disc material at the radius indicated on the y-axis. Both tilt and precession angles are continuous within the range of radii in the disc outside of stream’s pericentre (10R_g ≲ r ≲ 200R_g at the beginning of the simulation and 10R_g ≲ r ≲ 400R_g at the end of the simulation), indicating that there is no disc-tearing. Below 10R_g, tilt and precession angles fluctuate rapidly due to stream–disc interactions at pericentre. We do not observe significant precession of the disc in the duration of the simulation. It is possible that disc precession at small radii is inhibited by the constant injection of stream material in the orbital plane of the star.

Open in new tab Download slide

Figure A1.

The ratio of φ-velocity to circular Keplerian velocity for orbits of various eccentricities. We compute φ-velocity by setting up artificial velocity fields in the equatorial plane with a constant eccentricity. Then, we average φ-velocity over radius using a mass weight determined by the power law of best fit for the mass density radial profile in the aligned TDE simulation.

Open in new tab Download slide

Due to Lense–Thirring precession, we would expect the precession angle of the inner disc to increase monotonically with time. The precession rate of the line of nodes of a particle orbiting with eccentricity e and semimajor axis a_orb around a BH with dimension-less spin parameter a is given by Merritt (2013) as

$$\begin{eqnarray*} \frac{{\rm d} \Omega _{\rm LT}}{{\rm d} t} = \frac{2G^2M_{\rm BH}^2 a}{c^3 a_{\rm orb}^2 (1 - e)^{3/2}} . \end{eqnarray*}$$

(39)

Plugging in the dimension-less spin parameter a = 0.9375 and the mean eccentricity at each radius (Fig. 17), we find that the precession period exceeds the simulation duration of 5 d for radii r > 85R_g. The disc in our simulation radially extends to 400R_g at late times, so we should only expect significant precession within the inner regions of the disc. However, even at radii r < 85R_g, Fig. 23 shows that the precession angle remains consistent throughout the duration of our simulation. One explanation is that precession is inhibited by the angular momentum supplied by the debris stream which is in the initial orbital plane of the star.

At the times well before peak fallback time considered in our simulation, the fallback rate of the stellar debris increases linearly with time (Fig. 22). Therefore, the debris stream accounts for a significant proportion of the total angular momentum budget. By summing the components of Newtonian angular momentum in the region r < 500R_g, we estimate that the net angular momentum of the stream in this region is approximately 7 per cent of the net angular momentum of the disc at late times.

Previous GRMHD simulations of tilted accretion discs have shown disc tearing, where the accretion disc occupies separate planes over different ranges in radii (Liska et al. 2021). disc tearing occurs when the torque exerted on the disc by differing rates of Lense–Thirring precession at different radii surpasses the viscous forces holding the disc together. In particular, the inner part of the disc may become aligned with the equatorial plane of the BH, a phenomenon known as Bardeen–Peterson alignment. The continuity of the tilt and precession angles over the range of radii in the disc (Fig. 23) suggests that the disc remains intact. Similarly to precession, disc tearing and Bardeen–Peterson alignment may be inhibited by the contribution of angular momentum from the stream.

At times well after the end of our simulation, the mass accretion rate will drop below the Eddington limit. Therefore, the disc may cool and begin precessesing and/or tearing. Since this simulation was not run for multiple viscous times of the accretion disc, the presence of precession in tilted TDEs remains an open question.

5 CONCLUSIONS

In this work, we simulate a tidal disruption of a Sun-like star by a supermassive BH for a realistically large mass ratio (Q = 10⁶) and for a realistic stellar orbit (e₀ ≈ 1). We simulate the initial disruption in post-Newtonian SPH and migrate to full GRHD as the debris stream approaches the BH. We also present the first simulation of a tilted TDE in GRHD (Section 6). Our use of realistic parameters poses a number of challenges. A high mass ratio leads to a thin stellar debris stream that is difficult to resolve. We accommodate this difficulty using two levels of AMR. A parabolic stellar trajectory necessitates a large range of temporal and spatial scales. As the initial eccentricity of the star increases, the fallback time of the stellar debris and the apocentre of the debris stream orbit grow. The unprecedented efficiency of h-amr due to GPU-acceleration and AMR allows us to cover the necessary range of scales to simulate the earliest stages of accretion disc formation.

We find that the TDE naturally and efficiently forms an accretion disc, although the high-eccentricity material constantly supplied by the stream inhibits circularization. The accretion efficiency fluctuates between 10 and 20 per cent over the duration of our simulation (Fig. 18). We also find that a significant fraction of material accretes at moderate eccentricities (0.4 < e < 0.7), with highly eccentric material (e > 0.7) preferentially accreted at high latitudes (Fig. 20).

During the post-disruption phase of our aligned TDE simulation, the debris stream undergoes a series of violent self-intersection events in which the incoming and outgoing streams collide. We propose that these self-intersections are the phenomena responsible for the early-time flaring of TDE Swift J1644+57 and other TDEs. The self-intersections account for both the number of flares and their time-scale.

At early times (i.e. during the first 3 d of the simulation), self-intersections play a crucial role in orbital energy dissipation. At late times, the newly formed accretion disc completely intercepts the outgoing stream, causing the violent self-intersection events to cease. As a result, stream–disc interactions near pericentre are the dominant dissipation mechanism. These interactions raise the entropy of the debris stream three quarters of the way to the final entropy of the accretion disc (Fig. 13). Consequently, thermal energy flux dominates over mass energy flux in the debris stream post-pericentre compared to before the pericentre passage.

We find that the newborn disc exhibits super-Eddington accretion. The radial and polar dependencies of density and pressure within the disc closely reflect the self-similar solutions proposed by Coughlin & Begelman (2014) in the ZEBRA model. Non-zero eccentricity has a small effect on the deviation from Keplerian velocities in the disc, and the more prominent effects are from thermal pressure gradient forces (Fig. 21). The thermal energy generated by accretion heats up the inner part of the disc. The temperature distribution through the disc creates a thermal pressure gradient force that supports the disc against gravity, leading to a sub-Keplerian velocity distribution.

For a TDE with a 30 deg BH-spin–stellar-orbit misalignment angle, we find that nodal precession causes the incoming and outgoing streams to intersect off-centre (Fig. 9). This ejects gas from the outgoing stream on to orbits with larger tilt angles and results in less violent self-intersections. However, an accretion disc still forms with a similar accreted fraction of the material to the aligned case.

The largest drawback of our simulation is its short duration, of about one week, which allows only the small fraction of the stellar debris with specific energy well outside the frozen-in approximation to accrete on to the BH (Section 2.1). However, these early stages of TDE disc formation are crucial because they capture the initial disc formation and the emergence of the orbital energy dissipation mechanism. Our results suggest that disc formation in TDEs may be a runaway process; once sufficient mass has partially circularized at small radii, stream–disc interactions become the dominant dissipation mechanism, further growing the (initially eccentric) accretion disc. Clearly, the late-time evolution of TDE discs requires further study, particularly as the fallback rate approaches its peak and deviates from linear growth.

As discussed in Section 2.3, it is possible that we under resolve pressure gradients within the stream. This will be corrected in future simulations by including additional levels of AMR. In future work, we plan to study how TDE accretion disc formation is affected by the magnetic field of the disrupted star. This will allow us to model the MRI of the disc. We also plan to incorporate a variable polytropic index based on the equation of state used by Shiokawa et al. (2015) to more accurately model the thermodynamics of the disc. Finally, we will either use the pseudo-Newtonian potential from Tejeda & Rosswog (2013) in our initial SPH simulation or we will model the entire disruption in full general relativity.

SUPPORTING INFORMATION

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their detailed feedback which has greatly enhanced the clarity of this work. In particular, we thank the referee for the suggestion to analyse the eccentricity of accreting material in more detail. This research was made possible by National Science Foundation (NSF) Petascale Computing Resource Allocation awards no. 1615281 and OAC-1811605 as part of the Blue Waters sustained-petascale computing project, which is supported by the NSF (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Centre for Supercomputing Applications. Additional computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program under award PHY129. This research used resources of the Oak Ridge Leadership Computing Facility, which is a Department of Energy Office of Science User Facility supported under Contract DE-AC05-00OR22725.

Eric Coughlin ackowledges support from the National Aeronautics and Space Administration (NASA) through the Hubble Fellowship Program, grant #HST-HF2-51433.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555, and NSF grant AST-2006684. Alexander Tchekhovskoy is supported by the National Science Foundation grants AST-1815304 and AST-1911080. Nicholas Stone received financial support from NASA, through both Einstein Postdoctoral Fellowship Award Number PF5-160145 and the NASA Astrophysics Theory Research Program (Grant NNX17AK43G; Principal Investigator B. Metzger). He also received support from the Israel Science Foundation (Individual Research Grant 2565/19). Matthew Liska was supported by the John Harvard Distinguished Science Fellowship and the Institute for Theory and Computation Fellowship.

Our simulation, post-processing, and analysis relied on tools from the scipy and numpy Python libraries (Jones et al. 2001; Oliphant 2006). Our 2D and 3D visualizations relied on the matplotlib Python Library and the visit visualization software, respectively (Hunter 2007; Childs et al. 2012).

DATA AVAILABILITY

Movies and 3D renderings are available at MNRAS online and on YouTube. Simulation data are available upon request from Alexander Tchekhovskoy at [email protected].

Footnotes

1

The tidal tensor implied by the Einstein potential produces a moderately different tidal shear than the exact general relativistic value. However, the general relativistic corrections to the Newtonian tidal shear (either exact or from equation 1) are only |$\mathcal {O}(5{{\ \rm per\ cent}})$| at the tidal radius where the star is disrupted, and ballistic motion sets in.

REFERENCES

Alexander

K. D.

,

Wieringa

M. H.

,

Berger

E.

,

Saxton

R. D.

,

Komossa

S.

,

2017

,

ApJ

,

837

,

153

10.3847/1538-4357/aa6192

Crossref

Search ADS

Arcavi

I.

et al. ,

2014

,

ApJ

,

793

,

38

10.1088/0004-637X/793/1/38

Crossref

Search ADS

Ayal

S.

,

Livio

M.

,

Piran

T.

,

2000

,

ApJ

,

545

,

772

10.1086/317835

Crossref

Search ADS

Bade

N.

,

Komossa

S.

,

Dahlem

M.

,

1996

,

A&A

,

309

,

L35

Balbus

S. A.

,

Hawley

J. F.

,

1991

,

ApJ

,

376

,

214

10.1086/170270

Crossref

Search ADS

Bardeen

J. M.

,

Press

W. H.

,

Teukolsky

S. A.

,

1972

,

ApJ

,

178

,

347

10.1086/151796

Crossref

Search ADS

Bicknell

G. V.

,

Gingold

R. A.

,

1983

,

ApJ

,

273

,

749

10.1086/161410

Crossref

Search ADS

Blandford

R. D.

,

Begelman

M. C.

,

2004

,

MNRAS

,

349

,

68

10.1111/j.1365-2966.2004.07425.x

Crossref

Search ADS

Bloom

J. S.

et al. ,

2011

,

Science

,

333

,

203

10.1126/science.1207150

Crossref

Search ADS

PubMed

Bonnerot

C.

,

Lu

W.

,

2020

,

MNRAS

,

495

,

1374

10.1093/mnras/staa1246

Crossref

Search ADS

Bonnerot

C.

,

Lu

W.

,

2021

,

preprint (arXiv:2106.01376)

Bonnerot

C.

,

Stone

N.

,

2020

,

Formation of an Accretion Flow

.

Springer

,

Berlin

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bonnerot

C.

,

Rossi

E. M.

,

Lodato

G.

,

Price

D. J.

,

2016

,

MNRAS

,

455

,

2253

10.1093/mnras/stv2411

Crossref

Search ADS

Bonnerot

C.

,

Rossi

E. M.

,

Lodato

G.

,

2017

,

MNRAS

,

464

,

2816

10.1093/mnras/stw2547

Crossref

Search ADS

Bonnerot

C.

,

Lu

W.

,

Hopkins

P. F.

,

2021

,

MNRAS

,

504

,

4885

10.1093/mnras/stab398

Crossref

Search ADS

Brown

G. C.

,

Levan

A. J.

,

Stanway

E. R.

,

Tanvir

N. R.

,

Cenko

S. B.

,

Berger

E.

,

Chornock

R.

,

Cucchiaria

A.

,

2015

,

MNRAS

,

452

,

4297

10.1093/mnras/stv1520

Crossref

Search ADS

Burrows

D. N.

et al. ,

2011

,

Nature

,

476

,

421

Cannizzo

J. K.

,

Lee

H. M.

,

Goodman

J.

,

1990

,

ApJ

,

351

,

38

10.1086/168442

Crossref

Search ADS

Cao

R.

,

Liu

F. K.

,

Zhou

Z. Q.

,

Komossa

S.

,

Ho

L. C.

,

2018

,

MNRAS

,

480

,

2929

10.1093/mnras/sty1997

Crossref

Search ADS

Carter

B.

,

Luminet

J.-P.

,

1983

,

A&A

,

121

,

97

Cenko

S. B.

et al. ,

2012

,

ApJ

,

753

,

77

10.1088/0004-637X/753/1/77

Crossref

Search ADS

Childs

H.

et al. ,

2012

, in

Bethel

E. W.

,

Childs

H.

,

Hansen

C.

, eds,

High Performance Visualization–Enabling Extreme-Scale Scientific Insight

.

Taylor & Francis

,

Oxfordshire

, p.

357

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Coughlin

E. R.

,

Begelman

M. C.

,

2014

,

ApJ

,

781

,

82

10.1088/0004-637X/781/2/82

Crossref

Search ADS

Coughlin

E. R.

,

Nixon

C.

,

2015

,

ApJ

,

808

,

L11

10.1088/2041-8205/808/1/L11

Crossref

Search ADS

Coughlin

E. R.

,

Nixon

C.

,

Begelman

M. C.

,

Armitage

P. J.

,

Price

D. J.

,

2016

,

MNRAS

,

455

,

3612

10.1093/mnras/stv2511

Crossref

Search ADS

Cullen

L.

,

Dehnen

W.

,

2010

,

MNRAS

,

408

,

669

10.1111/j.1365-2966.2010.17158.x

Crossref

Search ADS

Dai

L.

,

McKinney

J. C.

,

Miller

M. C.

,

2015

,

ApJ

,

812

,

L39

10.1088/2041-8205/812/2/L39

Crossref

Search ADS

Evans

C. R.

,

Kochanek

C. S.

,

1989

,

ApJ

,

346

,

L13

10.1086/185567

Crossref

Search ADS

Fragile

P. C.

,

Anninos

P.

,

2005

,

ApJ

,

623

,

347

10.1086/428433

Crossref

Search ADS

Fragile

P. C.

,

Blaes

O. M.

,

Anninos

P.

,

Salmonson

J. D.

,

2007

,

ApJ

,

668

,

417

10.1086/521092

Crossref

Search ADS

Frank

J.

,

Rees

M. J.

,

1976

,

MNRAS

,

176

,

633

10.1093/mnras/176.3.633

Crossref

Search ADS

Gafton

E.

,

Rosswog

S.

,

2011

,

MNRAS

,

418

,

770

10.1111/j.1365-2966.2011.19528.x

Crossref

Search ADS

Gafton

E.

,

Rosswog

S.

,

2019

,

MNRAS

,

487

,

4790

10.1093/mnras/stz1530

Crossref

Search ADS

Gafton

E.

,

Tejeda

E.

,

Guillochon

J.

,

Korobkin

O.

,

Rosswog

S.

,

2015

,

MNRAS

,

449

,

771

10.1093/mnras/stv350

Crossref

Search ADS

Gezari

S.

et al. ,

2006

,

ApJ

,

653

,

L25

10.1086/509918

Crossref

Search ADS

Gezari

S.

et al. ,

2008

,

ApJ

,

676

,

944

10.1086/529008

Crossref

Search ADS

Gezari

S.

et al. ,

2012

,

Nature

,

485

,

217

Golightly

E. C. A.

,

Coughlin

E. R.

,

Nixon

C. J.

,

2019a

,

ApJ

,

872

,

163

10.3847/1538-4357/aafd2f

Crossref

Search ADS

Golightly

E. C. A.

,

Nixon

C. J.

,

Coughlin

E. R.

,

2019b

,

ApJ

,

882

,

L26

10.3847/2041-8213/ab380d

Crossref

Search ADS

Guillochon

J.

,

Ramirez-Ruiz

E.

,

2013

,

ApJ

,

767

,

25

10.1088/0004-637X/767/1/25

Crossref

Search ADS

Guillochon

J.

,

Ramirez-Ruiz

E.

,

2015

,

ApJ

,

809

,

166

10.1088/0004-637X/809/2/166

Crossref

Search ADS

Guillochon

J.

,

Ramirez-Ruiz

E.

,

Rosswog

S.

,

Kasen

D.

,

2009

,

ApJ

,

705

,

844

10.1088/0004-637X/705/1/844

Crossref

Search ADS

Guillochon

J.

,

Manukian

H.

,

Ramirez-Ruiz

E.

,

2014

,

ApJ

,

783

,

23

10.1088/0004-637X/783/1/23

Crossref

Search ADS

Hawley

J. F.

,

Krolkik

J. H.

,

2019

,

ApJ

,

878

,

149

10.3847/1538-4357/ab1f6e

Crossref

Search ADS

Hayasaki

K.

,

Stone

N.

,

Loeb

A.

,

2013

,

MNRAS

,

434

,

909

10.1093/mnras/stt871

Crossref

Search ADS

Hayasaki

K.

,

Stone

N.

,

Loeb

A.

,

2016

,

MNRAS

,

461

,

3760

10.1093/mnras/stw1387

Crossref

Search ADS

Hills

J. G.

,

1975

,

Nature

,

254

,

295

10.1038/254295a0

Crossref

Search ADS

Holoien

T. W.-S.

et al. ,

2014

,

MNRAS

,

445

,

3263

10.1093/mnras/stu1922

Crossref

Search ADS

Hunter

J. D.

,

2007

,

Comput. Sci. Eng.

,

9

,

90

10.1109/MCSE.2007.55

Crossref

Search ADS

Jiang

Y.-F.

,

Guillochon

J.

,

Loeb

A.

,

2016

,

ApJ

,

830

,

125

10.3847/0004-637X/830/2/125

Crossref

Search ADS

Jones

E.

,

Oliphant

T.

,

Peterson

P.

et al. ,

2001

,

SciPy: Open Source Scientific Tools for Python

.

Available at: http://www.scipy.org/

Kagaya

K.

,

Yoshida

S.

,

Tanikawa

A.

,

2019

,

preprint (arXiv:1901.05644)

Kochanek

C. S.

,

1994

,

ApJ

,

422

,

508

10.1086/173745

Crossref

Search ADS

Komossa

S.

,

Bade

N.

,

1999

,

A&A

,

343

,

775

Lacy

J. H.

,

Townes

C. H.

,

Hollenbach

D. J.

,

1982

,

ApJ

,

262

,

120

10.1086/160402

Crossref

Search ADS

Liptai

D.

,

Price

D. J.

,

Mandel

I.

,

Lodato

G.

,

2019

,

preprint (arXiv:1910.10154)

Liska

M.

et al. ,

2019

,

preprint (arXiv:1912.10192)

Liska

M.

,

Hesp

C.

,

Tchekhovskoy

A.

,

Ingram

A.

,

van der Klis

M.

,

Markoff

S. B.

,

Van Moer

M.

,

2021

,

MNRAS

,

507

,

983

10.1093/mnras/staa099

Crossref

Search ADS

Lu

W.

,

Bonnerot

C.

,

2020

,

MNRAS

,

492

,

686

10.1093/mnras/stz3405

Crossref

Search ADS

Magorrian

J.

,

Tremaine

S.

,

1999

,

MNRAS

,

309

,

447

10.1046/j.1365-8711.1999.02853.x

Crossref

Search ADS

Mainetti

D.

,

Lupi

A.

,

Campana

S.

,

Colpi

M.

,

Coughlin

E. R.

,

Guillochon

J.

,

Ramirez-Ruiz

E.

,

2017

,

A&A

,

600

,

A124

10.1051/0004-6361/201630092

Crossref

Search ADS

Merritt

D.

,

2013

,

Dynamics and Evolution of Galactic Nuclei

.

Princeton Univ. Press

,

Princeton

Merritt

D.

,

Alexander

T.

,

Mikkola

S.

,

Will

C. M.

,

2010

,

Phys. Rev. D

,

81

,

062002

10.1103/PhysRevD.81.062002

Crossref

Search ADS

Nealon

R.

,

Price

D. J.

,

Nixon

C. J.

,

2015

,

MNRAS

,

448

,

1526

10.1093/mnras/stv014

Crossref

Search ADS

Nelson

R. P.

,

Papaloizou

J. C. B.

,

2000

,

MNRAS

,

315

,

570

10.1046/j.1365-8711.2000.03478.x

Crossref

Search ADS

Nixon

C. J.

,

King

A. R.

,

2012

,

MNRAS

,

421

,

1201

10.1111/j.1365-2966.2011.20377.x

Crossref

Search ADS

Nowak

M. A.

,

Wagoner

R. V.

,

1991

,

ApJ

,

378

,

656

10.1086/170465

Crossref

Search ADS

Oliphant

T. E.

,

2006

,

A guide to NumPy

, vol.

1

.

Trelgol Publishing

,

USA

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Piran

T.

,

Sa̧dowski

A.

,

Tchekhovskoy

A.

,

2015

,

MNRAS

,

453

,

157

10.1093/mnras/stv1547

Crossref

Search ADS

Price

D. J.

,

2007

,

Publ. Astron. Soc. Austr.

,

24

,

159

10.1071/AS07022

Crossref

Search ADS

Price

D. J.

,

2012

,

J. Comput. Phys.

,

231

,

759

10.1016/j.jcp.2010.12.011

Crossref

Search ADS

Price

D. J.

,

Federrath

C.

,

2010

,

MNRAS

,

406

,

1659

10.1111/j.1365-2966.2010.16810.x

Crossref

Search ADS

Price

D. J.

,

Monaghan

J. J.

,

2007

,

MNRAS

,

374

,

1347

10.1111/j.1365-2966.2006.11241.x

Crossref

Search ADS

Price

D. J.

et al. ,

2018

,

Publ. Astron. Soc. Austr.

,

35

,

e031

10.1017/pasa.2018.25

Crossref

Search ADS

Ramirez-Ruiz

E.

,

Rosswog

S.

,

2009

,

ApJ

,

697

,

L77

10.1088/0004-637X/697/2/L77

Crossref

Search ADS

Rees

M. J.

,

1988

,

Nature

,

333

,

523

10.1038/333523a0

Crossref

Search ADS

Rosswog

S.

,

Ramirez-Ruiz

E.

,

Hix

W. R.

,

2009

,

ApJ

,

695

,

404

10.1088/0004-637X/695/1/404

Crossref

Search ADS

Ryu

T.

,

Krolik

J.

,

Piran

T.

,

Noble

S. C.

,

2020a

,

ApJ

,

904

,

98

10.3847/1538-4357/abb3cf

Crossref

Search ADS

Ryu

T.

,

Krolik

J.

,

Piran

T.

,

Noble

S. C.

,

2020b

,

ApJ

,

904

,

101

10.3847/1538-4357/abb3cc

Crossref

Search ADS

Saxton

R. D.

,

Read

A. M.

,

Esquej

P.

,

Komossa

S.

,

Dougherty

S.

,

Rodriguez-Pascual

P.

,

Barrado

D.

,

2012

,

A&A

,

541

,

A106

10.1051/0004-6361/201118367

Crossref

Search ADS

Shiokawa

H.

,

Krolik

J. H.

,

Cheng

R. M.

,

Piran

T.

,

Noble

S. C.

,

2015

,

ApJ

,

804

,

85

10.1088/0004-637X/804/2/85

Crossref

Search ADS

Steinberg

E.

,

Coughlin

E. R.

,

Stone

N. C.

,

Metzger

B. D.

,

2019

,

MNRAS

,

485

,

L146

10.1093/mnrasl/slz048

Crossref

Search ADS

Stone

N.

,

Loeb

A.

,

2012

,

Phys. Rev. Lett.

,

108

,

061302

10.1103/PhysRevLett.108.061302

Crossref

Search ADS

Stone

N.

,

Sari

R.

,

Loeb

A.

,

2013

,

MNRAS

,

435

,

1809

10.1093/mnras/stt1270

Crossref

Search ADS

Svirski

G.

,

Piran

T.

,

Krolik

J.

,

2017

,

MNRAS

,

467

,

1426

10.1093/mnras/stx117

Crossref

Search ADS

Sądowski

A.

,

Tejeda

E.

,

Gafton

E.

,

Rosswog

S.

,

Abarca

D.

,

2016

,

MNRAS

,

458

,

4250

10.1093/mnras/stw589

Crossref

Search ADS

Tchekhovskoy

A.

,

Narayan

R.

,

McKinney

J. C.

,

2011

,

MNRAS

,

418

,

L79

10.1111/j.1745-3933.2011.01147.x

Crossref

Search ADS

Tchekhovskoy

A.

,

Metzger

B. D.

,

Giannios

D.

,

Kelley

L. Z.

,

2014

,

MNRAS

,

437

,

2744

10.1093/mnras/stt2085

Crossref

Search ADS

Tejeda

E.

,

Rosswog

S.

,

2013

,

MNRAS

,

433

,

1930

10.1093/mnras/stt853

Crossref

Search ADS

Tejeda

E.

,

Gafton

E.

,

Rosswog

S.

,

Miller

J. C.

,

2017

,

MNRAS

,

469

,

4483

10.1093/mnras/stx1089

Crossref

Search ADS

van Velzen

S.

et al. ,

2011

,

ApJ

,

741

,

73

10.1088/0004-637X/741/2/73

Crossref

Search ADS

van Velzen

S.

et al. ,

2021

,

ApJ

,

908

,

4

10.3847/1538-4357/abc258

Crossref

Search ADS

Wevers

T.

et al. ,

2019

,

MNRAS

,

488

,

4816

10.1093/mnras/stz1976

Crossref

Search ADS

Wevers

T.

,

van Velzen

S.

,

Jonker

P. G.

,

Stone

N. C.

,

Hung

T.

,

Onori

F.

,

Gezari

S.

,

Blagorodnova

N.

,

2017

,

MNRAS

,

471

,

1694

10.1093/mnras/stx1703

Crossref

Search ADS

Wu

S.

,

Coughlin

E. R.

,

Nixon

C.

,

2018

,

MNRAS

,

478

,

3016

10.1093/mnras/sty971

Crossref

Search ADS

Xiang-Gruess

M.

,

Ivanov

P. B.

,

Papaloizou

J. C. B.

,

2016

,

MNRAS

,

463

,

2242

10.1093/mnras/stw2130

Crossref

Search ADS

Zanazzi

J. J.

,

Lai

D.

,

2019

,

MNRAS

,

487

,

4965

10.1093/mnras/stz1610

Crossref

Search ADS

Zanazzi

J. J.

,

Ogilvie

G. I.

,

2020

,

MNRAS

,

499

,

5562

10.1093/mnras/staa3127

Crossref

Search ADS

Zauderer

B. A.

et al. ,

2011

,

Nature

,

476

,

425

APPENDIX A: DERIVATIONS

In this appendix, we provide several derivations used in our analysis which did not fit into the main body of the text.

A1 Bernoulli parameter

Throughout our analysis, we used the Bernoulli parameter to distinguish between bound and unbound material. The Bernoulli parameter is the ratio of total energy flux to mass energy flux. The total energy flux in spatial coordinate xⁱ is given by

$$\begin{eqnarray*} \Phi _{\rm{total}}^i = -T_t^i = -(\rho c^2 + p + u_g) u_t u^i = -(\rho c^2 + \gamma u_g) u_t u^i, \nonumber\\ \end{eqnarray*}$$

(A1)

where T is the stress–energy tensor and g is the metric tensor. The total mass energy flux in spatial coordinate xⁱ is given by

$$\begin{eqnarray*} \Phi _{\rm{mass}}^i = \rho c^2 u^i \end{eqnarray*}$$

(A2)

Therefore, the relativistic Bernoulli parameter is given by

$$\begin{eqnarray*} b^{\prime } = \frac{\Phi _{\rm{total}}^i}{\Phi _{\rm{mass}}^i} = -\frac{(\rho c^2 + \gamma u_g) u_t u^i}{\rho c^2 u^i} = -\frac{u_t (\rho c^2 + \gamma u_g)}{\rho c^2} \end{eqnarray*}$$

(A3)

Thus defined, the relativistic Bernoulli parameter is counted off from unity, which corresponds to the rest-mass contribution: b′ > 1 corresponds to hydrodynamically unbound and b′ < 1 bound material. To follow the more familiar non-relativistic convention, we subtract the rest-mass contribution so that positive and negative values correspond to unbound and bound material, respectively:

$$\begin{eqnarray*} b = b^{\prime } - 1 = -\frac{u_t (\rho c^2 + \gamma u_g)}{\rho c^2} - 1 \end{eqnarray*}$$

(A4)

A2 Equatorial geodesics

In Fig. 5, we depict a geodesic in the equatorial plane. In the Kerr geometry, equatorial geodesics remain in the equatorial plane, so we set θ = π/2 and u^θ = 0. We then solve for u^t, u^r, and u^ϕ using the following equations.

$$\begin{eqnarray*} E= -g_{t\mu }u^\mu \end{eqnarray*}$$

(A5)

$$\begin{eqnarray*} L= g_{\phi \mu }u^\mu \end{eqnarray*}$$

(A6)

$$\begin{eqnarray*} g_{\mu \nu } u^\mu u^\nu = \kappa , \end{eqnarray*}$$

(A7)

where g is the metric tensor, E is energy, L is angular momentum, and κ = −1 for time-like geodesics. For the geodesic in Fig. 5, E and L are taken from their simulation values at the Cartesian point (−500, −200, 0). At each point along the geodesic, we compute u^r and u^ϕ and integrate the resulting differential equations. We linearly interpolate the covariant metric to the points along the geodesic.

A3 Artificial velocity fields

In Section 4.4, we create artificial velocity fields to control for the effect of the non-zero eccentricity of the disc in our analysis of its velocity distribution. We set up artificial velocity fields of constant eccentricity e and aligned pericentres under a Newtonian regime. For a given point in the midplane, we calculate the semimajor axis of the orbit a and the eccentric anomaly E from the distance from the BH r and the true anomaly ν = θ − π/2.

$$\begin{eqnarray*} a= \frac{r(1+e\cos {\nu })}{1-e^2} \end{eqnarray*}$$

(A8)

$$\begin{eqnarray*} E= \arctan {\frac{\sqrt{1-e^2}\sin {\nu }}{e+\cos {\nu }}} . \end{eqnarray*}$$

(A9)

Then, we compute the Cartesian state vectors.

$$\begin{eqnarray*} \mathbf {x}= r {\begin{pmatrix}\cos \nu \\ \sin \nu \end{pmatrix}} \end{eqnarray*}$$

(A10)

$$\begin{eqnarray*} \dot{\mathbf {x}}= \sqrt{\frac{a}{r}} {\begin{pmatrix}-\sin E\\ \sqrt{1-e^2} \cos E \end{pmatrix}} . \end{eqnarray*}$$

(A11)

Finally, we compute |$\dot{r}$| and |$\dot{\varphi }$|⁠.

$$\begin{eqnarray*} \dot{r} = \frac{x_1\dot{x}_1+x_2\dot{x}_2}{r} \end{eqnarray*}$$

(A12)

$$\begin{eqnarray*} \dot{\varphi } = \frac{x_1\dot{x}_2-\dot{x}_1 x_2}{r^2} . \end{eqnarray*}$$

(A13)

A4 Analytical self-intersection radius

The analytical self-intersection radius is given by Wevers et al. (2017). Consider the orbit of a mass-less test particle in the equatorial plane around the BH. Averaged over one orbit, general relativistic apsidal precession causes the argument of pericentre to advance by approximately an amount

$$\begin{eqnarray*} \delta \omega = A_S - 2 A_J , \end{eqnarray*}$$

(A14)

where A_S and A_J are the contributions to the apsidal precession of BH mass and spin-induced frame dragging, respectively, and the precession due to the BH’s quadrupole moment is ignored. To the lowest post-Newtonian order, A_S and A_J are given by Merritt et al. (2010) as

$$\begin{eqnarray*} A_S = \frac{6 \pi }{c^2}\frac{GM_{\rm{BH}}}{r_{\rm{p}}(1+e)} \end{eqnarray*}$$

(A15)

$$\begin{eqnarray*} A_J = \frac{4 \pi a}{c^3}\left(\frac{GM_{\rm{BH}}}{r_{\rm{p}}(1+e)}\right)^{3/2} , \end{eqnarray*}$$

(A16)

where e is the orbital eccentricity. From δω, we find the self-intersection radius with equation (A17).

$$\begin{eqnarray*} R_{\rm{SI}}=\frac{R_{\rm{p}} (1+e)}{1+e \cos (\pi + \delta \omega / 2)} . \end{eqnarray*}$$

(A17)

For a marginally bound stellar orbit e ≈ 1 with a pericentre radius r_p = 7R_g, we find a self-intersection radius of 142R_g. For a similar orbit with a pericentre radius r_p = 12R_g as seen in our simulation at late times, we find an analytical self-intersection radius of 565R_g.

APPENDIX B: CURVE FITTING DATA

In this appendix, we provide the curve fitting data for Figs 6 and 10 in Tables B1 and B2, respectively.

APPENDIX C: TILTING ALGORITHM

In this appendix, we describe the tilting algorithm used in our analysis of model TDET30. For our analysis of the tilted TDE, we untilt the data so that the orbital plane of the star lies in the equatorial plane. For each point on our original spherical grid, we convert to Cartesian coordinates and multiply by the rotation matrix |$\rm {R}_y(\pi /6)$|⁠. Then, we use a third-order spline method to interpolate our data to each point on the rotated grid.

In Fig. C1, we test our tilting algorithm by tilting and untilting one time slice. While the edges of the stream lose some of their definition, the overall structure of the system remains intact.

Figure C1.

Contour plots of the log of rest mass density in the equatorial plane at 5.7 d. The left-hand panel shows the unaltered data and the right-hand panel shows the data after two applications of the tilting algorithm at angles of π/6 and −π/6.

Open in new tab Download slide

Author notes

John Harvard Distinguished Science and ITC Fellow

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
December 2021	5
January 2022	45
February 2022	17
March 2022	7
April 2022	21
May 2022	11
June 2022	21
July 2022	10
August 2022	13
September 2022	8
October 2022	13
November 2022	9
December 2022	29
January 2023	19
February 2023	27
March 2023	27
April 2023	31
May 2023	37
June 2023	10
July 2023	15
August 2023	29
September 2023	16
October 2023	32
November 2023	25
December 2023	34
January 2024	31
February 2024	31
March 2024	30
April 2024	14
May 2024	53
June 2024	55
July 2024	49
August 2024	37
September 2024	51
October 2024	55
November 2024	31
December 2024	47
January 2025	14
February 2025	26
March 2025	39
April 2025	43
May 2025	33

Article Contents

Tidal disruption discs formed and fed by stream–stream and stream–disc interactions in global GRHD simulations

ABSTRACT

1 INTRODUCTION

2 NUMERICAL METHOD AND SETUP

2.1 Initial disruption in phantom

2.2 Mapping from phantom to h-amr

2.3 h-amr simulation parameters

3 RESULTS

3.1 Aligned disc formation and evolution

3.2 Tilted disc formation and evolution

4 DISCUSSION

4.1 Energy dissipation

4.2 Circularization

4.3 Accretion and outflow

4.4 Force balance

4.5 Comparison to disc models

4.5.1 Comparison to ZEBRA model

4.5.2 Bonnerot and Lu model

4.6 Analysis of the tilted TDE simulation, model TDET30

5 CONCLUSIONS

SUPPORTING INFORMATION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

REFERENCES

APPENDIX A: DERIVATIONS

A1 Bernoulli parameter

A2 Equatorial geodesics

A3 Artificial velocity fields

A4 Analytical self-intersection radius

APPENDIX B: CURVE FITTING DATA

APPENDIX C: TILTING ALGORITHM

Author notes

Supplementary data

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only