Simulating realistic disc formation in tidal disruption events

TL;DR

This paper presents the first realistic simulation of disc formation in tidal disruption events, revealing how debris evolves into a thick, marginally-bound disc with spiral shocks and an extended envelope, illuminating early TDE radiation.

Contribution

It introduces a novel simulation approach using outflow initial conditions to model realistic stellar trajectories and black hole masses in TDEs.

Findings

Rapid formation of a thick, marginally-bound disc with spiral shocks.

Presence of an extended envelope surrounding the accretion flow.

Shock heating rate consistent with observed TDE luminosities.

Abstract

A star coming too close to a supermassive black hole gets disrupted by the tidal force of the compact object in a tidal disruption event, or TDE. Following this encounter, the debris evolves into an elongated stream, half of which coming back to pericenter. Relativistic apsidal precession then leads to a self-crossing shock that initiates the formation of an accretion disc. We perform the first simulation of this process considering a realistic stellar trajectory and black hole mass, which has so far eluded investigations for computational reasons. This numerical issue is alleviated by using as initial conditions the outflow launched by the self-crossing shock according the local simulation of Lu & Bonnerot (2019). We find that the gas leaving the intersection point experiences numerous secondary shocks that result in the rapid formation of a thick and marginally-bound disc. The mass distribution features two overdensities identified as spiral shocks that drive slow gas inflow along the mid-plane. Inward motion primarily takes place along the funnels of the newly-formed torus, from which a fraction of the matter can get accreted. Further out, the gas moves outward forming an extended envelope completely surrounding the accretion flow. Secondary shocks heat the debris at a rate of a few times $10^{44} \, \rm erg \, s^{-1}$ with a large fraction likely participating to the bolometric luminosity. These results pave the way towards a complete understanding of the early radiation from TDEs that progressively becomes accessible from observations.

Figures

Equations

$$R_{\mathrm{t}}=R_{\star }{\left(\frac{M_{\mathrm{h}}}{M_{\star }}\right)}^{1/3}\approx 15R_{\mathrm{g}}=5.5\times 10^{12}\mathrm{cm},$$

$$\Delta ϵ=\frac{G{M}_{\mathrm{h}}{R}_{\star }}{R_{\mathrm{t}}^2},$$

$$a_{\min }=\frac{G{M}_{\mathrm{h}}}{2\Delta ϵ}=\frac{R_{\star }}{2}{\left(\frac{M_{\mathrm{h}}}{M_{\star }}\right)}^{1/3}\approx 85R_{\mathrm{t}}.$$

$$t_{\min }=2^{-1/2}\pi {\left(\frac{G{M}_{\star }}{R_{\star }^3}\right)}^{-1/2}{\left(\frac{M_{\mathrm{h}}}{M_{\star }}\right)}^{1/2}\approx 40\mathrm{d},$$

$${\dot{M}}_{\mathrm{fb}}={\dot{M}}_{\mathrm{p}}{\left(1+\frac{t}{t_{\min }}\right)}^{-5/3},$$

$${\dot{M}}_{\mathrm{p}}=\frac{M_{\star }}{3t_{\min }}\approx 25{\dot{M}}_{\mathrm{Edd}}=1.5{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},$$

$${\bar{\rho }}_{\mathrm{sc}}(\vartheta )=-0.08436x^4+0.2434x^3-0.2017x^2+0.1103x+0.01051$$

$${\dot{M}}_{\mathrm{sc}}={\dot{M}}_{\mathrm{fb}},$$

$$\Phi =-\frac{G{M}_{\mathrm{h}}}{R}-\left(\frac{2R_{\mathrm{g}}}{R-2R_{\mathrm{g}}}\right)\left[\left(\frac{R-R_{\mathrm{g}}}{R-2R_{\mathrm{g}}}\right){\dot{R}}^2+\frac{R^2({\dot{\theta }}^2+{\sin }^2\theta {\dot{\varphi }}^2)}{2}\right],$$

$$ϵ=\frac{1}{2}\left[\frac{R^2{\dot{R}}^2}{(R-2R_{\mathrm{g}}{)}^2}+\frac{R^3({\dot{\theta }}^2+{\sin }^2\theta {\dot{\varphi }}^2)}{R-2R_{\mathrm{g}}}\right]-\frac{G{M}_{\mathrm{h}}}{R},$$

$$\ell =\frac{R^3({\dot{\theta }}^2+{\sin }^2\theta {\dot{\varphi }}^2{)}^{1/2}}{R-2R_{\mathrm{g}}},$$

$${\left[m(\Omega -{\Omega }_{\mathrm{p}})\right]}^2={\kappa }^2+c_{\mathrm{s}}^2{k}^2,$$

$$R\approx R_0{e}^{-A(\varphi -{\varphi }_0)},$$

$${\dot{M}}_{\mathrm{in}}^{\mathrm{st}}\approx \alpha ⟨\rho c_{\mathrm{s}}⟩R^2\omega ,$$

$$\begin{array}{rl}\rho & =\frac{{\dot{M}}_{\mathrm{out}}}{4\pi R^2{v}_{\mathrm{out}}}\\ & \approx 2\times 10^{-13}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}{\left(\frac{\mathrm{R}}{50{\mathrm{R}}_{\mathrm{t}}}\right)}^{-2}{\left(\frac{{\mathrm{v}}_{\mathrm{out}}}{0.01c}\right)}^{-1}\left(\frac{{\dot{\mathrm{M}}}_{\mathrm{out}}}{0.5{\dot{\mathrm{M}}}_{\mathrm{p}}}\right),\end{array}$$

$${\int }_{R_{\mathrm{tr}}}^{+\infty }\rho {\kappa }_{\mathrm{s}}\mathrm{d}R=c/v(R_{\mathrm{tr}}),$$

$$R_{\mathrm{tr}}\approx \frac{{\dot{M}}_{\mathrm{out}}{\kappa }_{\mathrm{s}}}{4\pi c}\frac{v}{v_{\mathrm{out}}}=23R_{\mathrm{t}}\left(\frac{{\dot{M}}_{\mathrm{out}}}{0.5{\dot{M}}_{\mathrm{p}}}\right)\left(\frac{v}{0.03c}\right){\left(\frac{v_{\mathrm{out}}}{0.01c}\right)}^{-1},$$

$$\frac{{\rho }_{\mathrm{s}}}{\rho }=\frac{4{\dot{M}}_{\mathrm{fb}}{v}_{\mathrm{out}}{R}^2}{{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{s}}{H}_{\mathrm{s}}^2}\approx 8000\left(\frac{0.5{\dot{M}}_{\mathrm{fb}}}{{\dot{M}}_{\mathrm{out}}}\right)\left(\frac{10v_{\mathrm{out}}}{v_{\mathrm{s}}}\right){\left(\frac{0.05R_{\mathrm{t}}}{H_{\mathrm{s}}}\right)}^2,$$

$${\int }_{R_{\mathrm{th}}}^{\infty }\rho ({\kappa }_{\mathrm{a}}{\kappa }_{\mathrm{s}}{)}^{1/2}\mathrm{d}R=1,$$

Full text

††pubyear: 2019††pagerange: Simulating realistic disc formation in tidal disruption events–Simulating realistic disc formation in tidal disruption events

1 Introduction

Tidal disruption events (TDEs) occur when a star is disrupted by a supermassive black hole, after which the debris fuels the compact object producing a luminous multi-wavelength flare lasting from months to years (Rees, 1988). These events represent unique tools to study black holes in the center of otherwise quiescent galaxies as well as their dust and gas environments. Furthermore, they hold promise for improving our understanding of multiple processes occurring around compact objects, such as super-Eddington accretion and the formation of relativistic jets. Because the victim star inherits its final trajectory from two-body interactions taking place in the surrounding nuclear stellar cluster, TDEs can also be used to study the various dynamical processes at play in this region.

The sample of TDE candidates is steadily growing as more of them are detected at optical, UV and soft X-ray wavelengths (Bade et al., 1996; Komossa & Greiner, 1999; van Velzen et al., 2011; Gezari et al., 2006; Holoien et al., 2016; Gezari et al., 2017). Thanks to the advent of high-cadence surveys, a larger fraction of these observed events have well-sampled lightcurves that include the early rise to peak (Gezari et al., 2012; Chornock et al., 2014; Arcavi et al., 2014; van Velzen et al., 2019; Blagorodnova et al., 2019; Holoien et al., 2019) emerging when the stellar debris start radiating at an observable level. A few of them also displays hard X-ray emission that is attributed to relativistic jets launched from the inner region of the accretion flow (Bloom et al., 2011; Cenko et al., 2012; Brown et al., 2015).

The star is first scattered onto a plunging trajectory very close to parabolic after a series of stellar encounters taking place at large distances (Frank & Rees, 1976). Its disruption by the black hole imparts a spread in orbital energy to the debris given by its depth inside the gravitational potential well (Evans & Kochanek, 1989; Lodato et al., 2009; Guillochon & Ramirez-Ruiz, 2013). The gas subsequently follows a wide range of ballistic and highly eccentric trajectories to evolve into an elongated stream that remains thin owing to the confining effect of its self-gravity (Kochanek, 1994; Guillochon et al., 2014; Coughlin et al., 2016). Roughly half this debris is bound and thus falls back towards the black hole at a rate entirely specified by its orbital energy distribution (Rees, 1988; Phinney, 1989). A wide consensus has been reached regarding the hydrodynamics at play during these early stages of evolution and it is not strongly affected by the inclusion of additional physical ingredients (Coughlin & Nixon, 2015; Kesden, 2012; Tejeda et al., 2017; Gafton & Rosswog, 2019; Bonnerot et al., 2017b; Steinberg et al., 2019).

When the first debris reaches pericenter, it experiences relativistic apsidal precession that causes its collision with the part of the stream still approaching the compact object. This self-crossing shock has been investigated through local simulations (Kim et al., 1999; Jiang et al., 2016; Lu & Bonnerot, 2019), which find that it is close to adiabatic and therefore causes the formation of an outflow from the intersection point. This initial source of dissipation has been proposed to initiate the formation of an accretion disc, during which the gas trajectories progressively become more circular. This process has been investigated semi-analytically (Kochanek, 1994; Dai et al., 2015; Bonnerot et al., 2017a) but simulations to date have only been able to follow it for encounters with an intermediate-mass black hole (Rosswog et al., 2009; Guillochon et al., 2014; Shiokawa et al., 2015) or involving a star with an initially bound trajectory (Hayasaki et al., 2013; Bonnerot et al., 2016a; Hayasaki et al., 2016; Sadowski et al., 2016). These conditions are adopted for entirely numerical reasons because they artificially decrease the length of the stream that alleviates the huge computational cost required to accurately resolve it in the realistic case of a parabolic encounter with a supermassive black hole. These numerical investigations led nevertheless to significant improvements in our understanding of how disc formation occurs. Most importantly, they put forward major issues with the early picture of a fast circularization into a compact ring located near the tidal radius (Rees, 1988). It was found that the accretion flow extends to much larger distances and could require a few to many orbital timescales to assemble (Shiokawa et al., 2015; Bonnerot et al., 2016a). However, it remains uncertain how to quantitatively extrapolate the results of these simplified studies to the physically-motivated situation.

A clear understanding of the disc formation at play in real TDEs is paramount since this process connects the well-understood phases of stellar disruption and early stream evolution to the properties of the nascent accretion flow, from which most of the radiation is expected to emerge. For instance, it can inform the various models for the origin of the optical emission from these events that relies on either self-crossing shocks happening at large radii (Piran et al., 2015) or reprocessing of radiation generated close to the tidal radius (Metzger & Stone, 2016). The resulting gas distribution also affects the diffusion of photons produced near the black hole and whose interactions with the surrounding matter determines the spectrum of these events, in particular the expected strengths of the various lines (Roth et al., 2016). Furthermore, the outcome of the disc formation process specifies the correct initial conditions to be used in simulations of gas accretion in this newly-formed structure where matter is thought to spiral inwards under the effect of magnetic stresses that can also be accompanied with the formation of relativistic jets (Dai et al., 2018; Curd & Narayan, 2018).

In this paper, we carry out the first realistic simulation of disc formation following the disruption of a star on a parabolic orbit by a supermassive black hole. This is achieved using the smoothed-particle-hydrodynamics (SPH) method by injecting gas inside the computational domain according to the properties obtained from the local study of Lu & Bonnerot (2019) for the outflow launched from the self-crossing shock. With this procedure, we are able to circumvent the high computational cost required to follow the entire evolution of the stream that has so far prohibited simulations of realistic disc formation in TDEs. Our simulation captures the complex interactions experienced by the shocked gas after leaving the intersection point that allowed us to follow the formation of the accretion flow as well as several other hydrodynamical phenomena.

This paper is organized as follows. In Section 2, we explain the initial conditions and numerical setup used for our simulation as well as the injection procedure used to model the outflow launched from the self-crossing shock. The results are presented in Section 3 focusing on the disc assembly through secondary shocks, gas inflow and outflow as well as heating and radiation from the system. Section 4 includes a resolution study and a discussion of several physical effects that are not included in our simulation. Finally, we summarize our main findings in Section 5.

2 SPH simulations

We simulate the evolution of the stellar debris after it has passed through the self-crossing shock. This is achieved numerically by continuously injecting gas from an outflowing box located at the intersection point according to the properties of the shocked gas obtained from the local simulation of Lu & Bonnerot (2019).

2.1 Initial conditions

As a representative example, we consider the disruption of a star with mass M_{\star}=0.5\,\mathrm{M}_{\hbox{\odot}} and radius R_{\star}=0.46\,\mathrm{R}_{\hbox{\odot}} on a parabolic orbit by a black hole of mass M_{\rm h}=2.5\times 10^{6}\,\mathrm{M}_{\hbox{\odot}}. The pericenter distance $R_{\rm p}$ is equal to the tidal radius

[TABLE]

where $R_{\rm g}=GM_{\rm h}/c^{2}$ is the gravitational radius of the black hole. The corresponding penetration factor, defined as the ratio of tidal radius to pericenter distance, is $\beta=1$. Contrary to all previous numerical study of disc formation, this choice of parameters is physically motivated for a realistic disruption of a main sequence star by a supermassive black hole. Because we consider this already a significant step forward in the understanding of this process, we defer a more exhaustive exploration of the parameter space to future studies.

For this set of parameters, the stellar debris gets imparted an energy spread (Stone et al., 2013)

[TABLE]

by the disruption, putting the most bound part on a orbit with semi-major axis

[TABLE]

According to Kepler’s third law, it corresponds to an orbital period

[TABLE]

which we will use as our reference time-scale throughout the paper. Note that the all gas elements keep instead an angular momentum close to that of the original star given by $\ell_{\star}=(2GM_{\rm h}R_{\rm t})^{1/2}$ for a parabolic orbit since the redistribution of this orbital parameter is negligible during the disruption process. As a result, the gas distribution evolves into an elongated stream, whose bound half comes back to the original stellar pericenter. This fallback of matter proceeds at a rate given by (Rees, 1988; Phinney, 1989)

[TABLE]

where $t=0$ corresponds to the peak of the fallback rate given by

[TABLE]

as obtained by setting to half a stellar mass the integral of equation (5) from 0 to infinity. The Eddington accretion rate is obtained from \dot{M}_{\rm Edd}=L_{\rm Edd}/(\eta c^{2})\approx 0.06\,\rm\mathrm{M}_{\hbox{\odot}}\,{yr}^{-1} for a radiative efficiency of $\eta=0.1$ where $L_{\rm Edd}=3.2\times 10^{44}\,\rm erg\,s^{-1}$ denotes the Eddington luminosity. We note that equation (5) assumes a flat energy distribution within the debris for simplicity keeping in mind that the associated decay law for the fallback rate varies when the density structure of the star is considered (Lodato et al., 2009; Guillochon & Ramirez-Ruiz, 2013).

When the first debris passes at pericenter, relativistic apsidal precession results in its intersection with the part of the stream still approaching the black hole. The resulting self-crossing shock causes the formation of an outflow from the intersection point. It is illustrated in the sketch of Fig. 1 along with the procedure used to model this process that we describe in the next section. The properties of the outflow are predicted from the local simulation carried out by Lu & Bonnerot (2019). For our choice of parameters, this study finds that the intersection radius is located at a distance $R_{\rm int}=381R_{\rm g}\approx 25R_{\rm t}$ from the black hole. Following the collision, the shocked gas expands spherically at a speed $v_{\rm e}=0.065\,c$ while retaining a net velocity $v_{\rm n}=0.0152\,c$ in the direction tangential to the black hole. In a frame co-moving at velocity $v_{\rm n}\mathbf{e}_{\rm x}$, the outflow has an axisymmetric density profile with a polar dependence given by

[TABLE]

where $x\equiv\min(\vartheta,\pi-\vartheta)$ and $\vartheta$ denotes the polar angle defined in Fig. 1. The normalization is chosen such that $2\pi\int_{0}^{\pi}\bar{\rho}_{\rm sc}(\vartheta)\sin\vartheta\mathop{}\!\mathrm{d}\vartheta=1$. Because the gas that passes through the self-crossing shock is immediately ejected in this quasi-spherical outflow, the corresponding total outflow rate is given by

[TABLE]

where the fallback rate on the right-hand side is specified by equation (5). Note that this prescription relies on a continuous ejection of gas from the self-crossing shock while, in reality, the mass outflow from the intersection point is intermittent. This is because the self-crossing shock is present only when provided with a stream component approaching the black hole and another going away from it after one more pericenter passage. While the former is always present, the latter gets exhausted and then replenished on the time $t_{\rm rep}$ needed for the gas element just upstream the intersection point at the moment when the self-crossing shock ends to get around the black hole and back to the intersection point. The outflow rate should then be $\dot{M}_{\rm sc}=2\dot{M}_{\rm fb}$ when the self-crossing shock is on and $\dot{M}_{\rm sc}=0$ when it is off. The replenishment timescale is roughly twice the free-fall time from the intersection point, that is $t_{\rm rep}=2(R_{\rm int}^{3}/GM_{\rm h})^{1/2}=2\,\rm d$ for our choice of parameters. The fact that $t_{\rm rep}/t_{\rm min}\approx 0.05\ll 1$ implies that the outflow rate can be well approximated by a continuous one, justifying our treatment and the use of equation (8).

2.2 Numerical setup

The simulation is carried out with the SPH code phantom (Price et al., 2018) by adding gas to the computational domain according to the properties described above for the debris leaving the self-crossing shock. This procedure is illustrated in the sketch of Fig. 1. It consists in injecting SPH particles (red dots) from a spherical box (grey sphere) of radius $R_{\rm b}=10R_{\rm g}\ll R_{\rm int}$ located at the intersection point. The particles have a total velocity (red solid arrow) given by $\mathbf{v}=v_{\rm n}\mathbf{e}_{\rm x}+v_{\rm e}\mathbf{e}_{\rm r}$ that is the vectorial sum of the net and expansion speeds (blue dashed arrows) directed along in the direction tangential to the black hole and radially from the injection point, respectively. Additionally, their specific thermal energy is set such that the outflow is highly supersonic.

To reproduce the outflow rate from the self-crossing shock prescribed by equation (8), a number of particles $\Delta N=\left[\dot{M}\Delta t/M_{\rm p}\right]$ is added to the simulation at each timestep $\Delta t$, where the brackets represent the nearest integer function.111Our simulation uses individual time-stepping, implying that SPH particles can have different timesteps. The timestep $\Delta t$ between two particle injections is the shortest among them. The mass of the SPH particles is set to M_{\rm p}=10^{-8}\,\mathrm{M}_{\hbox{\odot}} in the run we present in the paper. Our simulation ends at $t=t_{\rm min}$ when $37\%$ of the bound debris have passed through the self-crossing shock, which corresponds to a total number of $N=9.3\times 10^{6}$ injected particles. The particle mass sets the resolution of the simulation with a lower $M_{\rm p}$ resulting in a higher resolution since more particles are injected during a given amount of time. In Section 4.1, we provide a resolution study by considering SPH particles of different masses and evaluate the impact on our results.

The density profile of the outflow depends on the number of particles injected at each timestep along the direction specified by the polar and azimuthal angles $\vartheta$ and $\varphi$. As described in the previous section, the desired density profile is axisymmetric, which we obtain by choosing $\varphi$ from a uniformly random distribution. The polar dependence given by equation (7) is instead attained by selecting $\vartheta$ according to the probability distribution $\mathop{}\!\mathrm{d}N/\mathop{}\!\mathrm{d}\vartheta=\bar{\rho}_{\rm sc}(\vartheta)\sin\vartheta$ where $\mathop{}\!\mathrm{d}N$ is the number of particles injected with a polar angle between $\vartheta$ and $\vartheta+\mathop{}\!\mathrm{d}\vartheta$. This is realized in practice by calculating the polar angle from $\vartheta=F^{-1}(s)$, where $s$ is a uniformly random number such that $0\leq s\leq 1$ and the function $F^{-1}$ is computed numerically to be the inverse of $s=F(\vartheta)\equiv 2\pi\int_{0}^{\vartheta}\bar{\rho}_{\rm sc}(\vartheta^{\prime})\sin\vartheta^{\prime}\mathop{}\!\mathrm{d}\vartheta^{\prime}$ (Press et al., 1992).

After it has been injected inside the computational domain, the gas evolves according to the gravity of the black hole that we model with the gravitational potential given by (Tejeda & Rosswog, 2013)

[TABLE]

where $R$ is the distance from the black hole and dotted variables are derived with respect to time. Here, $\theta$ and $\phi$ denote the polar and azimuthal angles of a test particle in spherical geometry. Note that they are different from those $\vartheta$ and $\varphi$ defined earlier that denote the location on the injection box. This potential reproduces many features of the Schwarschild metric including relativistic apsidal precession. It has been been implemented in phantom and tested by Bonnerot et al. (2016a). We also ran the same simulation adopting a Keplerian potential and found no significant differences. A test particle evolving under the gravitational potential of equation (9) conserves the specific orbital energy and angular momentum given by

[TABLE]

[TABLE]

respectively (Tejeda & Rosswog, 2013). As expected, these two quantities converge to the Keplerian values in the limit $R\gg R_{\rm g}$.

We assume that the gas evolves adiabatically away from shocks where non-adiabatic heating occurs. Because pressure is provided by radiation, we adopt an adiabatic exponent of $\gamma=4/3$. To treat shocks, we rely on a standard implementation of artificial viscosity. In combination, we make use the switch designed by Cullen & Dehnen (2010) that accurately detects shocks and only switches on artificial viscosity at their location. The assumption of adiabaticity is due to the optical thickness of the gas that prevents it from radiatively cooling. We justify this assumption in Section 3.4 where we carry out post-processing analyses to evaluate the impact of photon diffusion. Self-gravity is insignificant for this problem because it is overwhelmed by pressure gradients due to heating by shocks and we therefore do not consider it. Accretion onto the compact object is modelled by an accretion radius located at $R_{\rm acc}=6R_{\rm g}$ that is the innermost stable circular orbit for a non-rotating black hole. When an SPH particle enters this radius, it is removed from the simulation.

2.3 Early dynamics of the injected gas

Before fully presenting our results, it is worth emphasizing the early dynamics of the gas under the gravity of the black hole after it leaves the self-crossing shock according to the outflow properties described in Section 2.1 for our choice of parameters. Before the self-crossing, the fluid elements within the stream share the same energy and angular momentum given by the energy spread $-\Delta\epsilon$ and the stellar angular momentum $\ell_{\star}$, respectively. The outflow launched from the intersection point leads to a redistribution of these orbital elements around these values.

We display in Fig. 2 the corresponding distributions of orbital energy (upper panel) and angular momentum projected along the $\mathbf{e}_{\rm z}$ direction orthogonal to the original stellar orbital plane (lower panel) in the outflowing gas normalized such that the area below the curves is one. They are computed both from the local simulation of Lu & Bonnerot (2019) (black solid line, also shown with the green histogram in their figure 8) and directly from our simulation (red dashed line) at $t/t_{\rm min}=5\times 10^{-3}$. As can be seen from the upper left snapshot of Fig. 3, the gas is still close to the injection point at this time and has not yet experienced significant interactions such that each fluid element has conserved their energy and angular momentum from the moment of injection. The fact that the distributions obtained from these two methods match to a good accuracy demonstrates that our injection procedure is correctly implemented.

The early motion of the gas is entirely specified by the distributions of Fig. 2. In particular, it can be seen that a significant fraction around $33\%$ of the shocked gas has a positive orbital energy $\epsilon>0$ and is therefore unbound. This unbound component corresponds to matter being accelerated in the positive $\mathbf{e}_{\rm x}$ direction of its net speed during the outflow that results in a total velocity large enough to escape the black hole (Lu & Bonnerot, 2019). Additionally, some of the gas has a negative angular momentum $\ell_{\rm z}<0$, meaning that it rotates in the direction opposite to that of the star. It corresponds to matter launched near the $-\mathbf{e}_{\rm x}$ direction opposite to its net speed. Although most of the gas leaving the self-crossing shock is still co-rotating, the majority moves along the $\mathbf{e}_{\rm x}$ direction such that it gets unbound and leaves the system. As a result, the bound gas has on average a negative angular momentum and orbits the black hole in the counter-rotating direction. These initial orbital properties are determinant for the subsequent gas evolution studied in the rest of the paper.

3 Results

We now present the results of the simulation to determine the fate of the gas injected from the intersection point by the self-crossing shock.222Movies of the simulation presented in this paper are available at http://www.tapir.caltech.edu/~bonnerot/realistic-disc.html. The simulation is carried out using the same conventions as displayed in the sketch of Fig. 1 with the star and stream of debris rotating in the clockwise direction. The origin of the coordinate system is set at the location of the black hole and $t=0$ corresponds to the moment when the first fluid element is injected into the computational domain.

3.1 Disc assembly

3.1.1 Secondary shocks

After being injected into the computational domain from the self-crossing point, the gas moves on a wide range of trajectories that are prone to crossing each others, leading to secondary shocks. The resulting gas evolution is displayed in Figs. 3 and 4 that show its density at different times in slices parallel to the orbital plane of the original star and orthogonal to it, respectively. The first interaction happens at $t/t_{\rm min}\approx 0.03$, involving bound matter with the lowest angular momentum promptly reaching pericenter with opposite directions of rotation. The trajectories of this gas converge on the side of the black hole opposite to the intersection point leading to almost head-on collisions. The shocked region is initially confined near the black hole but rapidly expands in all directions. By $t/t_{\rm min}\approx 0.3$, this process leads to the appearance of a clearly defined boundary between the cold infalling matter and the heated gas located at $y\approx 15R_{\rm t}$ that is approximately halfway between the black hole and the injection point. In addition, interactions take place at smaller distances that cause more dissipation owing to the larger relative velocities involved.

Furthermore, this gas distribution displays an overdense filament inside the orbital plane of the original star, most visible in Fig. 3 on the left-hand side of the snapshots corresponding to times $0.1\leq t/t_{\rm min}\leq 0.4$. It is due to weak shocks taking place at radii $R\gtrsim 20R_{\rm t}$ that involve bound gas getting away from the self-crossing point and debris injected earlier that approaches the compact object after having reached its apocenter. The gas making up this filament continuously approaches the black hole and perturbs the rotating structure forming at lower distances. Remarkably, the associated overdensity seems to extend deeply into the accretion disc, as is most evident at $t/t_{\rm min}=0.4$. We attribute this phenomenon to the formation of spiral shocks in the forming gaseous disc due to the external perturbation produced by the returning debris, as we further analyse in Section 3.2.2.

3.1.2 Circularization in a counter-rotating disc

The numerous secondary shocks forming as the gas reaches the black hole lead to kinetic energy dissipation. As a result, the gas progressively assumes more circular orbits to form a coherent rotating structure that settles around $t/t_{\rm min}\approx 0.3$. Remarkably, because most of the injected bound gas rotates in the direction opposite to that of the original star as explained in Section 2.3, the resulting disc is also counter-rotating. This can be seen from Fig. 5 that shows the gas density at $t/t_{\rm min}=0.5$ in the newly-formed disc overplotted with black arrows of the same size that denote the direction of the gas velocity field. While the star’s trajectory was in the clockwise sense, the disc forming from the debris turns in the opposite, counter-clockwise way.

As can be seen from Fig. 3, this gaseous disc extends to a distance of $R\approx 15\,R_{\rm t}=0.2\,a_{\rm min}$. This property results from the spread in angular momentum imparted by the self-crossing shock as well as secondary shocks. It confirms that the gas does not settle into a circular ring at $R=2R_{\rm t}$. This latter assumption is commonly made in the literature but is only expected if angular momentum exchange is negligible from the stellar disruption onwards, which we find is not the case.333If all the fluid elements were able to circularize while keeping a common angular momentum equal to that of the star $\ell_{\star}=(2GM_{\rm h}R_{\rm t})^{1/2}$, the resulting disc would be a ring located at the circularization radius $R=\ell_{\star}^{2}/(GM_{\rm h})=2R_{\rm t}$. Because angular momentum gets efficiently redistributed among the debris through multiple shocks, we find that the trajectories of the gas within the disc span instead a large range of radii. Furthermore, Fig. 4 shows that the newly-formed disc is very thick with an aspect ratio $H/R\approx 1$ as a consequence of the thermal energy injection imparted by the circularizing shocks. It also displays two clearly defined funnels where the gas density is lower than along the disc mid-plane by up to two orders of magnitude.

Gas circularization can be evaluated more quantitatively from the evolution of the eccentricity distribution shown in Fig. 6 at different times and considering the debris located at a distance $R\leq 20R_{\rm t}$ from the black hole and within a polar angle $\pi/8$ from the equatorial plane, that is such that $|\theta-\pi/2|\leq\pi/8$. This region contains an increasing mass of gas but, for clarity, all the distributions are normalized such that the area below them are equal to one. Although the gas trajectories are sub-Keplerian because of the large pressure gradients pointing outward, we define an effective eccentricity that approximately describes the shape of the orbits by replacing the azimuthal component of the velocity $v_{\phi}$ by the larger value given by $v^{\prime}_{\phi}=(v^{2}_{\phi}+R|\nabla P/\rho|)^{1/2}$, where the pressure gradient is estimated by $R|\nabla P/\rho|\approx 8P/3\rho$ using the fact that the pressure scales as $P\propto\rho^{4/3}\propto R^{-8/3}$ since the gas evolves adiabatically and the density decreases roughly as $R^{-2}$ in most of the domain considered (see Section 3.3). From this increased azimuthal component, the eccentricity is then computed according to Keplerian dynamics through $e=\sqrt{1+2\epsilon\ell^{2}/(GM_{\rm h})^{2}}$, which is valid since most of the matter is located at distances from the black hole much larger than $R_{\rm g}$. At $t/t_{\rm min}=0.01$, none of the gas has circularized yet since it has not experienced any interaction. The peak around $e\approx 1$ is present because only the most eccentric gas can reach small radii within this time. We emphasize that it is not representative of the whole distribution of injected bound gas whose eccentricities are closer to $e\approx 0.5$. The distribution widens around $t/t_{\rm min}=0.03$ as more gas can make it in the region considered and the secondary shocks cause some of the orbits to become close to circular. As time goes on, the eccentricity distribution shifts to lower values while remaining significantly broad. It settles down around $t/t_{\rm min}=0.3$ that corresponds to the completion of disc formation according to the gas evolution presented above. The final distribution peaks around $e=0.2$ with a large spread around this value.

The gas retains non-negligible eccentricities as a result of highly irregular motion even after the debris has settled into a coherent rotating structure. This point is illustrated by the overplotted blue line in Fig. 5 that represents a typical trajectory for the gas inside the disc obtained by following a single SPH particle of our simulation. After joining the disc from the injection point, this gas element experiences large and violent variations in its orbital elements such that its trajectory covers a wide range of locations both within the mid-plane and away from it. Note that the line has sharp corners when the fluid element reaches small radii. This feature is artificial and results from the fact that this trajectory is recorded for time intervals larger than the orbital timescale of the particle. This sporadic motion of the gas inside the disc is also accompanied by variations of its angular momentum vector whose direction randomly changes that can result in a misalignment of up to a few degrees with the $\mathbf{e}_{\rm z}$ direction orthogonal to the stellar orbital plane.

3.1.3 Overall disc dynamics

The global dynamics of the disc is analysed more in details in Fig. 7 that shows with solid black lines the azimuthally-averaged profiles of specific orbital energy (upper panel) and angular momentum projected along the original stellar orbital plane (lower panel) obtained from equation (10) and (11) at $t/t_{\rm min}=0.5$ and considering the gas located within a polar angle $\pi/8$ from the equatorial plane. Note that the gas angular momentum is opposite to that of the star with $\ell_{\rm z}/\ell_{\star}<0$ because the disc is counter-rotating as explained above. The red lines denote the values of the orbital energy and angular momentum assuming circular motion in a Keplerian potential given by $\epsilon_{\rm c}=-GM_{\rm h}/(2R)$ and $\ell_{\rm c}=(GM_{\rm h}R)^{1/2}$, respectively. They are larger than that of the gas at all distances that is a consequence of both its eccentric trajectories and the fact that the disc is supported by pressure gradients in addition to the centrifugal force. For $R\lesssim 5R_{\rm t}$, the gas angular momentum nevertheless approaches its circular value with $\ell_{\rm z}\approx\ell_{\rm c}$, implying that centrifugal support dominates in this inner region. At larger distances, this profile is nearly flat with an inversion around $R\approx 20R_{\rm t}$ due to the contribution from the low angular momentum gas directly coming from the injection point on almost radial trajectories.

The blue dotted line in the upper panel of Fig. 7 shows the Bernoulli number of the gas given by $\mathcal{B}=\epsilon+4P/\rho$ for the adiabatic exponent $\gamma=4/3$. The second term takes into account the energy increase resulting from heating by secondary shocks. This number is close to zero throughout the disc, slowly declining outwards. The decrease at $R\approx R_{\rm t}$ is caused by the pressure drop associated to the removal of the gas crossing the accretion radius. At small radii $R\lesssim 5R_{\rm t}$, it becomes positive due to the large internal energy injection happening in this region. The fact that $\mathcal{B}\approx 0\gg\epsilon_{\rm c}$ for most of the gas implies that it can reach large distances or even get unbound from the system if its heat is converted into kinetic energy. This situation is prone to strong gas outflow that will be further examined in Section 3.3. The marginally-bound nature of this gaseous structure is also consistent with analytical models of TDE discs, such as the one developed by Coughlin & Begelman (2014).

3.2 Inflow and accretion

3.2.1 Gas inflow

The gas present in the inner regions moves on average towards the black hole. We compute the rate at which this matter flows inward from $\dot{M}_{\rm in}=(\Delta M_{\rm in}-\Delta M_{\rm out})/\Delta t$ where $\Delta M_{\rm in}$ and $\Delta M_{\rm out}$ denote the masses of gas that cross a sphere of given radius during a time-step $\Delta t$ going inward and outward, respectively. The radial profile of this total net inflow rate time-averaged on the interval $0.4\leq t/t_{\rm min}\leq 0.5$ is displayed in Fig. 8 (black line) that also shows the same calculation considering only the gas originating from the funnels (red dashed line) and the disc (blue dotted line). As before, gas is defined to belong to the disc if it is located within a polar angle $\pi/8$ from the equatorial plane while the rest of the matter is part of the funnels. The total inflow rate increases only slowly with radius with a value $\dot{M}_{\rm in}\lesssim 0.1\,\dot{M}_{\rm p}\approx 2.5\,\dot{M}_{\rm Edd}$. However, its origin changes drastically with distance from the black hole. At $R\gtrsim 10R_{\rm t}$, most of the inward motion occurs close to the mid-plane while the gas inside the funnels moves either slower inward or outward. This inflow is mostly caused by matter joining the outer edges of the disc either coming directly from the injection point or from larger distances after being previously expelled by secondary shocks. This hierarchy inverses closer in where inflow along the funnels dominates that inside the disc by about an order of magnitude, the latter decreasing to $\sim 10^{-2}\dot{M}_{\rm p}\approx 0.25\,\dot{M}_{\rm Edd}$. We highlight that this difference cannot be accounted for by the different solid angles covered by the two regions since that of the funnels is larger by a factor of only $(1-\sin(\pi/8))/\sin(\pi/8)\approx 1.6$. The larger inflow rate along the funnels in the inner region is therefore caused by a higher mass flux. This polar angle dependence of the inflowing mass flux has also been found in the simulation of disc formation by Sadowski et al. (2016).

The main characteristics of the gas radial motion can be seen directly from Fig. 9 that shows the gas density in a slice orthogonal to the equatorial plane at $t/t_{\rm min}=0.5$ with black arrows of different lengths indicating both direction and magnitude of the velocity field. While there is hardly any inward motion visible inside the disc, the matter contained in the funnels is clearly falling towards the black hole that is possible due to lower pressure and centrifugal support in these regions. The blue line indicates the trajectory of a given SPH particle that is representative of the inflowing gas. After being injected from the intersection point, this fluid element joins the disc and revolves several times around the black hole. However, it eventually gets ejected from the mid-plane to the upper funnel, along which it quickly falls inward all the way inside the accretion radius. The orange line shows the trajectory of another gas element that enters the lower funnel and gets accreted in one single infall after being injected from the self-crossing point with only a mild change in its orbital elements. We find that this second possibility is, however, less likely with most of the mass inside the funnels being injected at least a few dynamical timescales earlier.

3.2.2 Spiral shocks

The disc displays two spirals of higher density that form around $t/t_{\rm min}=0.3$ and remain present until $t/t_{\rm min}\approx 0.6$. They can be seen by looking at the density distribution displayed in Fig. 10 for $t/t_{\rm min}=0.4$ (see also the snapshots in Fig. 3). In particular, the inset shows a zoom-in on a part of the most evident left spiral with gas density rescaled to make it more visible with the density only increasing by a factor of $\sim 2$ from blue to orange. The external part of this feature coincides with the location where the filament of bound matter joins the disc outer edge, as discussed in Section 3.1.1. Although clearly present, the spiral on the opposite side is weaker and more intermittent. This density pattern remains at the same location around the black hole with no rotation. We identify these overdensities as spiral shocks driven by the external perturbation caused by the stream of injected bound gas that continuously strikes the outer edge of the disc on the left-hand side. This is supported by the fact that the gas dissipates part of its kinetic energy when crossing the density features, proving that they are shocks.

Spiral shocks can be localized using the dispersion relation of the density wave from which they develop that is given by

[TABLE]

equation (6.55) of Binney & Tremaine (2008) without considering self-gravity where $\Omega$ and $\Omega_{\rm p}$ are the gas and pattern angular frequencies while $\kappa^{2}=R^{-3}\mathop{}\!\mathrm{d}(R^{4}\Omega^{2})/\mathop{}\!\mathrm{d}R$ represents the squared epicyclic frequency. Here, $k$ and $m$ denote the radial and azimuthal wave numbers of the perturbation, respectively. We adopt $\kappa\approx\Omega$ that corresponds to a gas angular momentum profile close to Keplerian. Even though it is less accurate in the outer parts of the disc where the orbits are significantly sub-Keplerian (lower panel of Fig. 7), this assumption only weakly affects our conclusion on the location of the spiral shocks. Additionally, $\Omega_{\rm p}\approx 0$ because the perturbation does not rotate around the black hole. Equation (12) then leads to $k\approx\Omega(m^{2}-1)^{1/2}/c_{\rm s}$ such that the location of the density wave is given by $\mathop{}\!\mathrm{d}\phi/\mathop{}\!\mathrm{d}R=-k/m\approx-\mathcal{M}(m^{2}-1)^{1/2}/(Rm)$ where $\mathcal{M}=R\Omega/c_{\rm s}$ denotes the Mach number and $\phi$ is the azimuthal angle measured from the positive $\mathbf{e_{\rm x}}$ direction and increasing in the counter-clockwise way for the point of view adopted in Fig. 10. The solution describes a trailing spiral of equation

[TABLE]

where $A=m/(\mathcal{M}(m^{2}-1)^{1/2})$ while $R_{0}$ and $\phi_{0}$ are the inner radius and azimuthal angle of the feature. For a given value of $m$, the shape of the resulting spiral is determined by the Mach number. Increasing it decreases $A$ that makes the spiral more tightly wound.

Following Ju et al. (2016), we fit the overdensities visible in the disc using equation (13) assuming $m=2$ since two spirals are visible. The result is shown in Fig. 10 where the two dashed red lines correspond to the fitted spirals that use $R_{0}=R_{\rm t}$ and $\phi_{0}=0$ while the Mach number is set to $\mathcal{M}=1.3$ that is similar to that measured from the simulation, which is expected since $\mathcal{M}\approx H/R\approx 1$ in accretion discs (Frank et al., 2002). The quality of this fit is satisfying considering the violent gas motion and large variations of its properties within the forming disc. It therefore favours our interpretation that the density pattern seen in the simulation is due to spiral shocks.

3.2.3 Viscosity

The gas passing through the spiral shocks experiences a change in velocity that can be seen directly from the inset of Fig. 10 whose black arrows display both the magnitude and direction of the velocity field. To evaluate this effect more quantitatively, we compute the perturbations in the radial and azimuthal components of the velocity compared to the background flow from $\delta v_{\rm R}=v_{\rm R}-\langle v_{\rm R}\rangle$ and $\delta v_{\rm\phi}=v_{\rm\phi}-\langle v_{\rm\phi}\rangle$, where the brackets represent time and azimuthal averages. Fig. 11 shows the evolution of the radial (black solid line) and azimuthal (red solid line) velocity components time-averaged in the interval $0.4\leq t/t_{\rm min}\leq 0.5$ as a function of azimuthal angle for the gas contained in a small annulus centred at radius $R=10R_{\rm t}$ and within a polar angle $\pi/8$ from the equatorial plane. The horizontal dashed lines of the same color denote the azimuthal averages of these quantities, namely $\langle v_{\rm R}\rangle$ and $\langle v_{\phi}\rangle$. The velocity perturbations $\delta v_{\rm R}$ and $\delta v_{\phi}$ are simply the difference between these two curves as indicated by the black and red arrows. The vertical grey dotted lines represent the azimuthal angles of the two spirals obtained from the fitted red dashed curves of Fig. 10. Both radial and azimuthal velocities decrease around these locations with the associated velocity perturbations taking negative values, that is $\delta v_{\rm R}<0$ and $\delta v_{\rm\phi}<0$. The most prominent spiral feature on the left causes the largest and most extended velocity variations while the weaker overdensity on the opposite side leads to smaller and more localized perturbations.

It has long been known that spiral shocks can lead to outward angular momentum transport in astrophysical discs. The reason is that the components of the induced velocity perturbations are correlated such that $\langle\delta v_{\rm R}\delta v_{\rm\phi}\rangle>0$, as can be clearly seen from Fig. 11. To evaluate the amount of associated inward motion, we compute the viscous parameter $\alpha=T_{\rm R\phi}/\langle P\rangle$, where $\langle P\rangle$ denotes the average pressure while $T_{\rm R\phi}=\langle\rho\delta v_{\rm R}\delta v_{\rm\phi}\rangle$ is the Reynolds stress tensor. We find that $\alpha\approx 0.02$ within $R\lesssim 15R_{\rm t}$ that corresponds to the location of the disc, thus confirming that the velocity perturbations are correlated at all distances in this region. The corresponding stress results in an inflow of gas inside the disc at a rate

[TABLE]

evaluating the mass flux from $\langle\rho v_{\rm r}\rangle\approx\alpha\langle\rho c_{\rm s}\rangle$ where the radial velocity is obtained from $v_{\rm r}\approx\nu/R\approx\alpha c_{\rm s}$ that evaluates the kinematic viscosity from the prescription $\nu\approx\alpha c_{\rm s}H$ (Frank et al., 2002) and uses the fact that $H/R\approx 1$ in the disc. Here, $\omega$ denotes the solid angle of the disc given by $\omega=4\pi\sin(\pi/8)$ since this region contains by definition the gas located within a polar angle $\pi/8$ from the mid-plane. Using the radial dependence of the parameter $\alpha$ derived above, we calculate the inflow rate induced by the stress as a function of radius. The resulting profile time-averaged in the interval $0.4\leq t/t_{\rm min}\leq 0.5$ is shown with a blue dash-dotted line in Fig. 8 to compare it with that obtained directly from the gas motion inside the disc (blue dashed line). At $R\lesssim 10R_{\rm t}$, the two rates roughly agree with $\dot{M}^{\rm st}_{\rm in}\approx\dot{M}_{\rm in}$, demonstrating that gas inflow is induced by the passage of the debris through the spiral shocks where their angular momentum is transported outward. As anticipated in Section 3.2.1, the inward motion at larger radii is instead dominated by matter falling from large distances and penetrating through the disc with a residual radial velocity.

Because magnetic fields are not included in our simulation, we cannot capture the development of the magneto-rotational instability MRI, Balbus & Hawley 1991. The gas orbital period in the disc is of order $t/t_{\rm min}\approx 0.01$ within $R\leq 5R_{\rm t}$, which implies that MRI rapidly reaches saturation there. The resulting angular momentum transport is likely to significantly enhance the viscosity compared to that computed above and only due to spiral shocks. This point is discussed further in Section 4.5. A definitive evaluation of the amount of viscosity inside the disc requires magneto-hydrodynamics simulations of the process, which we defer to future investigations.

3.2.4 Accretion

A fraction of the inflowing gas can make it close enough to be accreted onto the black hole. In our simulation, a fluid element is accreted if it crosses the accretion radius located at $R_{\rm acc}=6R_{\rm g}$. The evolution of the corresponding accretion rate is shown with a black solid line in Fig. 12 that also displays the inflow rates measured at a distance of $R=3R_{\rm t}\approx 45R_{\rm g}$ through the funnels (red dashed line) and the disc (blue dotted line). Accretion onto the black hole is dominated by gas moving along the funnels while inward motion inside the disc provides only a small contribution, as anticipated in Section 3.2.1. At early times, the mass accretion rate exhibits a small peak at $\dot{M}_{\rm acc}\approx 0.06\,\dot{M}_{\rm p}=1.5\dot{M}_{\rm Edd}$ due to the infall of debris that has recently been injected and can rapidly reach the inner regions. This is possible since the disc acting as a centrifugal barrier has not yet fully formed, as can be seen in Fig. 4 at $t/t_{\rm min}=0.1$. Matter accumulates in the inner region as evidenced by the green dash-dotted line showing the total inflow rate at $R=20R_{\rm t}$, which is larger than that measured at lower radii. This accumulation leads to an increase of the accretion rate to $\dot{M}_{\rm acc}\approx 0.2\,\dot{M}_{\rm p}=5\dot{M}_{\rm Edd}$ that largely compensates the outer inflow of matter from $t/t_{\rm min}\approx 0.5$. As a result the mass of non-accreted gas enclosed inside $R=20R_{\rm t}$ settles to $M_{\rm enc}\approx 0.02M_{\star}$ that remains roughly constant until the end of the simulation at $t/t_{\rm min}=1$.

One might expect that most of the gas accreted in ballistic infall through the funnels retains a significant fraction of its orbital energy when disappearing inside the black hole. For example, the gas parcel following the blue trajectory in Fig. 9 appears to get accreted in a single ballistic infall starting around $R_{\rm ff}\approx 10R_{\rm t}=150R_{\rm g}$ conserving its orbital energy $\epsilon_{\rm ff}$ during the process. If this were true, this fluid element would dissipate only a small fraction $\epsilon_{\rm ff}/\epsilon_{\rm c}\approx R_{\rm acc}/R_{\rm ff}=0.04$ of the maximum possible energy $\epsilon_{\rm c}$ of a circular orbit at $R_{\rm acc}$. However, we do not find that this is the case in our simulation. Instead, the large majority of the mass is accreted with an orbital energy similar to $\epsilon_{\rm c}$ after it has dissipated its original excess energy by shocks taking place in the inner region of the disc. This behaviour is different from that found by Sadowski et al. (2016) who put forward this low dissipation efficiency as a possible origin for the low radiative efficiency of the disc. Instead, our simulation shows that accretion is often associated to an almost optimal thermal energy injection into the surrounding gas. Note that it is nevertheless possible that the heated gas itself enters the accretion radius before it is able to significantly cool, thus again reducing the radiative efficiency of accretion, as further examined in Section 3.4.

3.3 Outflow

As explained in Section 2, the self-crossing shock drives a massive outflow. Out of this injected mass, about $33\%$ is initially unbound from the black hole and gets ejected near the $\mathbf{e_{\rm x}}$ direction. An additional source of outflow relates to secondary shocks that can launch part of the initially bound debris to large distances with a fraction being unbound. In fact, unbinding most of the matter is energetically possible if a small fraction of mass circularizes into a Keplerian orbit near $R_{\rm t}$ and transfers its dissipated orbital energy to the rest of the gas. We may therefore expect that this supplementary heating significantly contributes to the outflowing gas motion.

We compute the outflow rate through a sphere located at a distance of $R=50R_{\rm t}$ from the black hole from $\dot{M}_{\rm out}=\Delta M_{\rm out}/\Delta t$ where $\Delta M_{\rm out}$ is the mass crossing this radius going outward within a time $\Delta t$. Fig. 13 displays its evolution for the gas reaching this distance inside the right (blue solid line) and left (red solid line) hemispheres at $x>0$ and $x<0$, respectively. Both outflow rates decline following the injection rate shown with the grey dot-dashed line and analytically computed from equation (8). This is because most of the injected mass reaches distances of order $a_{\rm min}\approx 85R_{\rm t}>50R_{\rm t}$ as it goes away from the self-crossing point in a quasi-spherical outflow. The contribution from the right hemisphere is larger by a factor of a few than that from the left one because the gas has a net velocity along the $\mathbf{e}_{\rm x}$ direction.

Only a fraction of this outward-moving gas is unbound from the black hole. This can be seen from Fig. 13 that also shows the outflowing rates through the right (blue dashed line) and left (red dashed line) hemispheres considering only the gas with positive Bernoulli number $\mathcal{B}>0$. Even though we later refer to this gas as being unbound, we highlight that it is not necessarily the case if this matter experiences interactions at larger radii or radiates a significant fraction of its internal energy (see Section 3.4). Both of these processes are likely to result in $\mathcal{B}<0$ after the gas has crossed the sphere such that it in fact remains bound to the black hole. The component crossing the right hemisphere is dominated by gas originally unbound by the self-crossing shock that is confined to the region with $x>0$. It remains close to the initial value of $33\%$ of the injection rate at all times. The unbound outflow rate through the left hemisphere cannot be due to gas coming directly from the self-crossing shock since this matter does not reach $x<0$. It therefore has necessarily gained its energy from interactions happening near the black hole. The spike at $t/t_{\rm min}\lesssim 0.1$ is produced by gas involved in the first secondary shocks that gains sufficient energy to escape near the $-\mathbf{e}_{\rm y}$ direction as can be seen from Fig. 4 at $t/t_{\rm min}=0.03$. Later on, the large majority of the outflow through this hemisphere is bound due to the confinement by later-arriving debris. This unbound outflow rate increases again from $t/t_{\rm min}\approx 0.5$ due to the enhanced dissipation occurring inside the disc after it has settled. Nevertheless, this unbound component remains at all times about an order of magnitude lower than that coming directly from the self-crossing point.

The energy content of the outward-moving matter can be further analysed by looking at Fig. 14 that shows the distribution of Bernoulli number for the gas contained in a slice parallel to the equatorial plane at $t/t_{\rm min}=0.7$. The black arrows of the same size represent the direction of the velocity field while the orange dashed circle is located at $R=50R_{\rm t}$ where the outflow rates of Fig. 13 are evaluated. Most of the gas with $\mathcal{B}>0$ (red) is produced by the self-crossing shock and able to cross this sphere at $x>0$. Unbound gas is also present near the black hole due to heating by secondary shocks as expected from the upper panel of Fig. 7. Some of it forms pockets that are able to move outward to larger radii with $x<0$. Nevertheless, most of the matter in this region remains bound with $\mathcal{B}<0$ (blue) after being injected.

The outflowing matter reaches large distances from the black hole that produces a large-scale envelope around it. Fig. 15 shows density profiles at $t/t_{\rm min}=0.5$ for different radial directions obtained by considering the gas inside cones of opening angle $\pi/8$ and centred along the unit vectors $\mathbf{e}_{\rm x}$ (black line), $-\mathbf{e}_{\rm x}$ (red line), $(\mathbf{e}_{\rm z}-\mathbf{e}_{\rm x})/\sqrt{2}$ (blue line), $(\mathbf{e}_{\rm x}+\mathbf{e}_{\rm z})/\sqrt{2}$ (green line) and $\mathbf{e}_{\rm z}$ (orange line). For $R\lesssim 20R_{\rm t}$, the density clearly decreases for lower polar angles, being the highest along the $\pm\mathbf{e}_{\rm x}$ ($\theta=\pi/2$) directions and the lowest along the $\mathbf{e}_{\rm z}$ ($\theta=0$) direction. This is due to the formation of the funnels that evacuates the polar regions during disc formation. Around $R\approx 30R_{\rm t}$, the density is similar along all directions that reflects the quasi-spherical nature of the outflow. In this region, the densities roughly agree with the grey dashed line that displays the analytical profile corresponding to a steady-state spherical wind given by

[TABLE]

where the outflow rate is set to $\dot{M}_{\rm out}=0.5\dot{M}_{\rm p}$ and the gas outward velocity to $v_{\rm out}=0.01c$ as estimated from the simulation. The density profiles keep following this analytical prescription at larger radii for the $\mathbf{e}_{\rm x}$, $(\mathbf{e}_{\rm x}+\mathbf{e}_{\rm z})/\sqrt{2}$ and $\mathbf{e}_{\rm z}$ directions corresponding to $x\geq 0$ since the gas in this region is either weakly bound or unbound and can therefore reach large distances. However, the density declines faster for $R\gtrsim 40R_{\rm t}$ for the $-\mathbf{e}_{\rm x}$ and $(\mathbf{e}_{\rm z}-\mathbf{e}_{\rm x})/\sqrt{2}$ directions that correspond to $x<0$. This is because the gas in this region is the most bound and has therefore the lowest apocenter distances that prevent the gas from moving to large radii. Overall, the properties of this outflow are consistent with the description adopted by Metzger & Stone (2016) who assume that most of the matter gets launched in a wind after falling back near the black hole. However, we emphasize that the outflow is not steady-state but possesses instead complex energy and density distributions as can be seen from Fig. 14.

3.4 Heating and radiation

A large amount of heat is injected by secondary shocks. Some of this thermal energy excess can be radiated since photons are transported outward by a combination of advection and diffusion until some of them eventually leak out entirely from the system. While it incorporates heating by shocks, our simulation does not capture photon diffusion and assumes instead that radiation is advected with the flow at every location by adopting an adiabatic equation of state with pressure dominated by radiation. Nevertheless, we perform simple post-processing calculations both to evaluate the validity of this assumption and approximately estimate the radiative luminosity that emerges during the disc formation process from the shock-heating obtained from the simulated gas evolution.

The heating rate $\dot{E}_{\rm sh}$ resulting from all the shocks experienced by the gas can be computed directly from our simulation using equation (42) of Price et al. (2018) thanks to the shock detector implemented in our code. Its evolution is shown with a black solid line in Fig. 16 along with that of the initial self-crossing shock (red dashed line) computed analytically from $\dot{M}_{\rm sc}v^{2}_{\rm e}/2\approx 2\times 10^{44}\,\rm erg\,s^{-1}(1+t/t_{\rm min})^{-5/3}$. The heating rate due to secondary shocks exhibits a small bump to $\dot{E}_{\rm sh}\approx 2\times 10^{44}\,\rm erg\,s^{-1}$ visible at $t/t_{\rm min}\approx 0.1$ related to the early inflow of gas near the black hole also visible in the net inflow rate along the funnels at $R=3R_{\rm t}$ shown with a dashed red line in Fig. 12. As the disc builds up, more mass gets funnelled in the inner regions and shock-heated that leads to an increase of the associated heating rate to its maximal value of $\dot{E}_{\rm sh}\approx 4\times 10^{44}\,\rm erg\,s^{-1}\approx L_{\rm Edd}$ at $t/t_{\rm min}\approx 0.7$. From that moment, it remains roughly constant until the end of the simulation as the disc density settles while the rate of injected mass drops. The heating rate due to the initial self-crossing shock drops following the decline in the fallback rate to become about an order of magnitude lower than that produced by secondary shocks. This is because most shocks in our simulation take place closer to the black hole than the intersection point where relative velocities are larger leading to more heating.

We now want to estimate the fraction of this generated heat that participates to the emerging luminosity. Close to the black hole where densities are the largest, numerous scatterings by free electrons prevent photons from moving with respect to the gas. At larger radii, the decreased density allows radiation to decouple from the matter and diffuse away. The transition occurs at the trapping radius where the diffusion timescale $t_{\rm diff}\approx\tau R/c$ equals the dynamical timescale $t_{\rm dyn}\approx R/v$. Here, $\tau$ denotes the optical depth while $v$ is the gas velocity. This equality yields $\tau\approx c/v$ such that the trapping radius $R_{\rm tr}$ is given by

[TABLE]

where $\kappa_{\rm s}=0.34\,\rm cm^{2}\,g^{-1}$ denotes the electron scattering opacity and the velocity on the right-hand side is evaluated at the trapping radius. We compute this critical distance by performing the integral of equation (16) from outside in until it reaches the local ratio of speed of light to gas velocity. This calculation is carried out considering different radial paths corresponding to given values of $\theta$ and $\phi$ that allows us to specify the trapping surface, from which photons are able to diffuse from the gas.

The trapping surface is shown in Fig. 17 at different times with lines overplotted on the gas distribution of radiation temperature evaluated at $t/t_{\rm min}=0.5$ along slices parallel (upper panel) and orthogonal (lower panel) to the equatorial plane and computed from $T_{\rm rad}=(e/a)^{1/4}$ where $e$ is the internal energy density while $a=7.6\times 10^{-15}\,\rm erg\,cm^{-3}\,K^{-4}$ denotes the radiation constant. The red solid line corresponds to $t/t_{\rm min}=0.1$ when the trapping surface is growing primarily along the $-\mathbf{e}_{\rm y}$ direction as freshly-injected gas reaches this location, as can be seen directly from the corresponding density snapshots of Figs. 3 and 4. The maximal extent of this region is reached around $t/t_{\rm min}=0.5$ (blue dashed line), from which it starts receding inward as the density drops due to the decline in the injection rate. As a result, the trapping surface is the smallest at $t/t_{\rm min}=1$ (orange dot-dashed line). The trapping radius is only weakly dependent on the radial path considered because the gas envelope surrounding the black hole has an approximately spherical density structure as explained in Section 3.3. In particular, we note that the trapping surface is not affected by the presence of the funnels at low radii. The average value of the trapping radius can be recovered by plugging the analytical density profile of equation (15) into equation (16) that yields

[TABLE]

where the total velocity of $v=0.03c$ is evaluated based on the simulation near the trapping radius at $t/t_{\rm min}=0.5$. This value is indicated by a vertical grey dotted line on the density profiles of Fig. 15. At all times, the trapping radius encloses the region where most of the interactions between debris take place, which implies that our assumption of adiabaticity is valid for most of the hydrodynamics described in the above sections. One exception concerns the unbound outflow that we define in Section 3.3 based on the conservation of the Bernoulli number that assumes adiabaticity. Taking into account energy loss through radiation outside the trapping radius would reduce the unbound fraction compared to what we estimated. Similarly, the radiation temperature is computed accurately from the inner regions where $T_{\rm rad}\approx 10^{6}\rm K$ to the trapping surface where it has decreased to $T_{\rm rad}\approx 10^{5}\rm K$ due to the adiabatic losses in the expansion driven by radiation pressure. However, our assumption of advective transport breaks down outside the trapping radius where photons decouple from the gas. This implies that our calculation of the radiation temperature is not self-consistent in these outer parts as discussed further in Section 4.4.

The trapping radius is only marginally larger than the outer edge of the disc where most of the dissipation takes place. We therefore do not expect significant adiabatic losses as radiation is advected to this surface before it diffuses outward. To evaluate this effect, we show in Fig. 16 with blue crosses the evolution of the emerging bolometric luminosity obtained by integrating on the trapping surface the diffusion flux given by $F_{\rm dif}=c|\nabla e|/(3\rho\kappa_{\rm s})\approx\xi ec/\tau$ equation (I.24) of Shapiro & Teukolsky (1983) computing the energy density directly from our simulation.444Radiation can be absorbed outside the trapping radius before it reaches the observer. However, these photons are rapidly re-emitted without energy loss and the bolometric luminosity is therefore unaffected by this process. In the second term, we fix $\xi=8/9$ as obtained from the approximate scaling $e\propto\rho^{4/3}\propto R^{-8/3}$. The bolometric luminosity being close to the heating rate from secondary shocks (solid black line) at all times confirms our expectation of small energy loss due to adiabatic expansion. Additionally, it suggests that advection and subsequent accretion of hot gas onto the black hole does not lead to a significant reduction of the energy output.

The calculations presented above to estimate how photon diffusion occurs in the forming disc are approximate and we do not pretend to have solved radiation transport in TDEs. To achieve this goal, the disc formation process has to be studied by radiation-hydrodynamics simulations, which we defer to the future.

4 Discussion

4.1 Resolution study

To evaluate the influence of numerical resolution on our results, we repeat the same simulation with both a larger and lower number of particles. This is done by modifying their mass to M_{\rm p}=3\times 10^{-9}\,\mathrm{M}_{\hbox{\odot}} and M_{\rm p}=3\times 10^{-8}\,\mathrm{M}_{\hbox{\odot}} that corresponds to a higher and lower resolution compared to that M_{\rm p}=10^{-8}\,\mathrm{M}_{\hbox{\odot}} adopted in the rest of the paper. We find that the dynamics of disc formation is quantitatively similar for the two larger resolutions with a well-defined torus forming that has the same mass $M_{\rm enc}=0.02M_{\star}$ enclosed within $R\leq 20R_{\rm t}$ at $t/t_{\rm min}=0.5$. The gas evolution is however significantly different in the lower resolution simulation, for which the gas clusters near the black hole but does not form a coherent structure. As a result, the enclosed mass is smaller than in the two other runs with $M_{\rm enc}=0.008M_{\star}$ at the same time. This influence of resolution is similar to that obtained in the early SPH simulation by Ayal et al. (2000), who find that a disc does not form when the matter falls back near pericenter with a very low number of particles. We conclude that the resolution adopted in the simulation presented in the paper is high enough that the results are not affected by increasing it further.

4.2 Validity of the injection procedure

Our treatment of the self-crossing shock assumes that it is continuously joined by two streams that remain identical at all times. This approach assumes that the returning debris is unaffected by the presence of the matter that has already accumulated near the black hole. The validity of this assumption can be evaluated from the density ratio between these two components. As can be seen from Fig. 15, the density profile of the gas in our simulation is accurately described by equation (15). The density of the stream can be estimated in a similar way from $\rho_{\rm s}=\dot{M}_{\rm fb}/(\pi H^{2}_{\rm s}v_{\rm s})$ where $H_{\rm s}$ and $v_{\rm s}$ are its width and velocity assuming cylindrical symmetry. The density ratio is then

[TABLE]

evaluating it at $R\approx 5R_{\rm t}$. The outflow rate is $\dot{M}_{\rm out}\approx 0.5\dot{M}_{\rm fb}$ at all times as evidenced by Fig. 13. The stream moves at a speed given by $v_{\rm s}\approx(2GM_{\rm h}/5R_{\rm t})^{1/2}\approx 0.1c$ that is about an order of magnitude larger than that of the outflowing gas, of $v_{\rm out}\approx 0.01c$. Furthermore, the stream is able to maintain a thin profile until it reaches the vicinity of the black hole (Coughlin et al., 2016) such that $H_{\rm s}\approx 10R_{\star}\approx 0.05R_{\rm t}$. The three ratios on the right-hand side of equation (18) are therefore close to unity such that the stream is significantly denser than the mass already present around the black hole when it arrives. This calculation shows that it is legitimate to neglect in first approximation the dynamical impact of the accumulated mass on the freshly-returning stream. It is nevertheless possible that this debris can be affected by hydrodynamical instabilities taking place at the boundary between the stream and the gas present in the disc or surrounding envelope. These interactions are similar to those investigated by Bonnerot et al. (2016b) during the early evolution of the stream. When the stream becomes very tenuous at late times while mass keeps accumulating near the compact object, the ratio between the two densities is however likely to decrease, thus enhancing the impact of the surrounding matter on the trajectories of the returning debris.

4.3 Extrapolation to other parameters

Our choice of parameters corresponds to an intersection radius $R_{\rm int}\approx 25R_{\rm t}$ that is much lower than the apocenter of the most bound debris of $2a_{\rm min}\approx 170R_{\rm t}$. This is due to the large apsidal relativistic precession experienced by the stream when reaching its pericenter at a distance of $R_{\rm p}\approx 15R_{\rm g}$. As a result, the initial self-crossing shock launches the debris on a wide range of trajectories with the bound ones promptly reaching small radii. As demonstrated by Lu & Bonnerot (2019), this description likely applies to more relativistic disruptions such as those involving penetration factors $\beta>1$ or black holes of larger masses than the one M_{\rm h}=2.5\times 10^{6}\,\mathrm{M}_{\hbox{\odot}} we consider. We therefore argue that our simulation provides at least the correct qualitative gas evolution at play during disc formation for this region of parameter space. Similar simulations have been carried out by Bonnerot et al. (2016a) and Sadowski et al. (2016) who consider the disruption of bound stars assuming a pericenter well within the tidal radius such that $R_{\rm p}\approx 10R_{\rm g}$, only slightly smaller than what we use. Our results are in several respects similar to theirs, most importantly the fact that the disc forms on a timescale of a few tenths of $t_{\rm min}$.

Decreasing the black hole mass to M_{\rm h}\lesssim 10^{6}\,\mathrm{M}_{\hbox{\odot}} while keeping $\beta=1$ reduces the apsidal precession angle such that the intersection radius rapidly increases to $R_{\rm int}\gtrsim a_{\rm min}$. As described by Lu & Bonnerot (2019), the initial self-crossing shock becomes weaker in this situation such that the debris that passes through it moves on a narrower range of trajectories. This type of mildly relativistic encounters has been investigated by several authors considering either initially bound stars (Bonnerot et al., 2016a; Hayasaki et al., 2016) or a black hole of low mass with M_{\rm h}=500\,\mathrm{M}_{\hbox{\odot}} (Shiokawa et al., 2015). These works find that about ten $t_{\rm min}$ are required for the debris to reach a significant level of circularization. It is legitimate to think that this property carries over to standard TDEs when the self-crossing shock is weak. However, extrapolating the detailed hydrodynamics found in these simulations is uncertain due to the simplifications they make. Most importantly, decreasing the stellar eccentricity results in all the debris stream being bound that artificially suppresses the tail of later-arriving matter. Considering an intermediate-mass black hole modifies the location of the self-crossing shock due to the dependence on the black hole mass of both the apsidal precession angle and the semi-major axis $a_{\rm min}$. It is also likely to increase the dynamical importance of dissipation occurring during the stream compression at pericenter (Guillochon et al., 2014). For these reasons, the detailed dynamics at play for this region of parameter space remains largely unclear.

4.4 Properties of the emerging radiation

Our simulation demonstrates that the heating rate by secondary shocks quickly increases to becomes larger than that of the initial self-crossing shock by almost an order of magnitude (Fig. 16). Furthermore, we argue based on post-processing calculations that the bolometric luminosity is likely similar to this shock heating rate due to the proximity of the trapping radius to the newly-formed disc. This favours the interpretation of Metzger & Stone (2016) for the origin of the optical emission from TDEs that relies on dissipation near the tidal radius and then radiative transport by advection and diffusion. Piran et al. (2015) argue instead that this emission is caused by the initial self-crossing shock when it occurs on larger scales similar to $a_{\rm min}$. However, our simulation demonstrates that this initial source of heating is weaker in addition to being more susceptible to adiabatic losses (Lu & Bonnerot, 2019). For this reason, we argue that this option is less likely to account for observations.

We find that the radiation temperature is $T_{\rm rad}\approx 10^{5}\rm K$ at the trapping surface (Fig. 17). As photons are transported further out by diffusion, they can be absorbed and re-emitted that reduces their temperature until the thermalization radius $R_{\rm th}$ given by equation (1.98) of Rybicki & Lightman (1979)

[TABLE]

where $\kappa_{\rm a}$ denotes the absorption opacity that is much lower than that $\kappa_{\rm s}$ due to scattering. In order to reproduce the optical emission detected from TDEs, this radius has to be a factor of a few larger than the trapping radius figure 13 of Lu & Bonnerot (2019). A more accurate determination of the spectrum could be achieved using the density distribution obtained from our simulation by carrying out post-processing radiative transfer calculations similar to those by Roth et al. (2016).

4.5 Influence of magnetic fields

Our simulation does not consider magnetic fields that the debris inherits from the star upon disruption. As mentioned in Section 3.2.3, including them would likely enhance angular momentum transport due to the development of the MRI compared to that produced by hydrodynamical effects alone. Furthermore, the presence of magnetic fields inside the newly-formed disc could lead to additional outflows through the Blandford & Payne (1982) mechanism. However, it has been found that a weak or zero net poloidal magnetic field reduces both the angular momentum transport resulting from the MRI and the strength of magnetically-driven winds (Hawley et al., 1995; Salvesen et al., 2016; Jiang et al., 2017; Zhu & Stone, 2018).

This situation is likely to apply to TDEs because the magnetic field lines get aligned with the direction of stream stretching prior to disc formation that suppresses the vertical field (Guillochon & McCourt, 2017; Bonnerot et al., 2017b). For this reason, it is possible that the dynamical influence of magnetic fields is reduced compared to other accreting systems that do not lack this poloidal component. This argument is supported by the simulation of Sadowski et al. (2016) who finds that the dynamics of the newly-formed disc is unaffected by the inclusion of magnetic fields even after the MRI has reached saturation. The possibility that the early accretion has an hydrodynamical origin has also been proposed by Nealon et al. (2018).

4.6 Impact of black hole spin

Our initial conditions assume that the self-crossing shock happens after the tip of the stream has passed pericenter and moves outward to intersect with matter still falling back. However, this initial collision can be prevented if the black hole has a non-zero spin due to Lense-Thirring precession that modifies the orbital plane of the debris during pericenter passage (Dai et al., 2013; Guillochon & Ramirez-Ruiz, 2015; Hayasaki et al., 2016). In this case, the stream may experience a large number of revolutions before it intersects with itself. However, Guillochon & Ramirez-Ruiz (2015) show that this collision is more likely between neighbouring windings that differ by only one pericenter passage. In this situation, the properties of the delayed self-crossing shock may be similar to the non-rotating case that we are simulating. One important difference relates to the presence of multiple streams that did not collide and could occupy a significant fraction of the volume directly surrounding the black hole. Interactions involving this gas could affect the dynamics of the forming disc, making it significantly more complex.

5 Summary

We carried out the first simulation of disc formation in TDEs considering realistic initial conditions valid for a star disrupted on a parabolic trajectory by a supermassive black hole. This is achieved by following the gas evolution only after it has passed through the self-crossing shock that we model by injecting SPH particles inside the computational domain according to the properties derived from the numerical study of the initial collision by Lu & Bonnerot (2019). This technique allows us to circumvent the high computational cost required to accurately capture the previous passage of the stream at pericenter. The results can be summarized as follows.

Due to the large range of trajectories followed by the gas leaving the self-crossing shock, multiple secondary shocks take place that lead to the formation of an accretion disc at $t/t_{\rm min}\approx 0.3$. This structure rotates in the direction opposite to that of the original star and retains significant eccentricities broadly distributed around $e\approx 0.2$ even after the overall gas evolution has settled. It is very thick with an aspect ratio $H/R\approx 1$ and extends to distances $R\approx 15R_{\rm t}$. 2. 2.

Pressure support is important in most of the disc such that its angular momentum profile is significantly sub-Keplerian with an almost flat radial dependence. Most of the gas in this region has a Bernoulli number $\mathcal{B}\approx 0$ that is much larger than its orbital energy due to the excess heat injected by the circularizing shocks that makes the disc prone to outflows. 3. 3.

Two overdensities develop inside the newly-formed disc that we identify as spiral shocks excited by the arrival of bound mass from the injection point that continuously strikes its outer edge. The velocity changes experienced by the debris moving through these features cause outward angular momentum transport that can account for the inward motion of gas along the mid-plane at a slow rate of $\sim 10^{-2}\dot{M}_{\rm p}=0.25\dot{M}_{\rm Edd}$. 4. 4.

Gas inflow is dominated on small scales by matter falling along the funnels of the torus at a rate of $\dot{M}_{\rm in}\approx 0.1\dot{M}_{\rm p}=2.5\dot{M}_{\rm Edd}$ that can result in accretion onto the black hole. The majority of this gas stayed in the disc for several orbital periods after it has been injected from the self-crossing shock. Moreover, accretion does not happen on a single near-radial infall but after most of the gas orbital energy has been deposited in the inner region of the disc. 5. 5.

At larger distances $R\gtrsim 50R_{\rm t}$, the gas moves outward at a rate a factor of a few lower than the fallback rate. It is dominated by matter directly coming from the self-crossing shock that contains a significant fraction of unbound gas with a lower contribution from debris gaining energy from secondary shocks. The resulting density profile can be approximately described by a steady-state spherical wind coming from the black hole vicinity. 6. 6.

Secondary shocks cause gas heating at a rate peaking at $\dot{E}_{\rm sh}\approx 4\times 10^{44}\,\rm erg\,s^{-1}\approx L_{\rm Edd}$ that is larger than that due to the self-crossing shock by almost an order of magnitude. Post-processing calculations suggest that a large fraction contributes to the bolometric luminosity due to the proximity of the trapping radius to the outer edge of the disc that reduces adiabatic losses.

This study and the numerical method it uses open several promising lines of research that we intend to pursue in the future. This includes the exploration of the hydrodynamics of disc formation for a wider range of parameter space, most importantly less relativistic encounters than we assume here. Additionally, we envisage simulations of this process treating in a self-consistent way radiative diffusion and angular momentum transport including magnetic fields. Major advances are foreseen that will improve our understanding of the early emission from TDEs now accessible by observational facilities.

Acknowledgments

We thank Eliot Quataert, Sterl Phinney, Phil Hopkins, Nathan Roth, Luc Dessart, Chris White and Greg Salvesen for useful discussions. We acknowledge the use of splash (Price, 2007) for producing most of the figures in this paper. This research benefited from interactions at the ZTF Theory Network Meeting, partly funded by the National Science Foundation under Grant No. NSF PHY-1748958. The research of CB was funded by the Gordon and Betty Moore Foundation through Grant GBMF5076. WL was supported by the David and Ellen Lee Fellowship at Caltech.