Self-intersection of the Fallback Stream in Tidal Disruption Events

TL;DR

This paper presents a semi-analytical model explaining how stream self-intersection in tidal disruption events leads to observable optical and infrared features, unifying diverse behaviors and predicting new spectral characteristics.

Contribution

It introduces a new model for stream self-intersection in TDEs, linking it to observed optical/infrared properties and explaining the diversity in X-ray emissions.

Findings

CIO causes optical emission reprocessing and explains observed line widths.

Infrared spectrum of CIO-reprocessed emission follows a shallow power law.

Optical TDE rate is lower than total rate due to filtering by the reprocessing layer.

Abstract

We propose a semi-analytical model for the self-intersection of the fallback stream in tidal disruption events (TDEs). When the initial periapsis is less than about 15 gravitational radii, a large fraction of the shocked gas is unbound in the form of a collision-induced outflow (CIO). This is because large apsidal precession causes the stream to self-intersect near the local escape speed at radius much below the apocenter. The rest of the fallback gas is left in more tightly bound orbits and quickly joins the accretion flow. We propose that the CIO is responsible for reprocessing the hard emission from the accretion flow into the optical band. This picture naturally explains the large photospheric radius (or low blackbody temperature) and typical widths of the H and/or He emission lines seen in optical TDEs. We predict the CIO-reprocessed spectrum in the infrared to be L_\nu ~ \nu^{~0.5}, shallower than a blackbody. The partial sky coverage of the CIO also provides a unification of the diverse X-ray behaviors of optical TDEs. According to this picture, optical surveys filter out a large fraction of TDEs with low-mass blackholes due to lack of a reprocessing layer, and the volumetric rate of optical TDEs is nearly flat wrt. the blackhole mass for M < 10^7 solar masses. This filtering causes the optical TDE rate to be lower than the total rate by a factor of ~10 or more. When the CIO is decelerated by the ambient medium, radio emission at the level of that in ASASSN-14li may be produced, but the timescales and peak luminosities can be highly diverse. Finally, our method paves the way for global simulations of the disk formation process by injecting gas at the intersection point according to the prescribed velocity and density profiles.

Equations

$$\frac{r_{\mathrm{T}}}{r_{\mathrm{g}}}\equiv \frac{R_{\ast }}{r_{\mathrm{g}}}{\left(\frac{M}{M_{\ast }}\right)}^{1/3}=46.7M_6^{-2/3}{m}_{\ast }^{-1/3}{r}_{\ast }.$$

$${\mathcal{E}}_{\mathrm{T}}\equiv \frac{r_{\mathrm{g}}{R}_{\ast }}{r_{\mathrm{T}}^2}=2.13\times 10^{-4}{M}_6^{1/3}{m}_{\ast }^{2/3}{r}_{\ast }^{-1},$$

$$\ell =r_{\mathrm{p}}\sqrt{(1+\mathcal{E}{)}^2/{\mu }_{\mathrm{p}}-1},{\mu }_{\mathrm{p}}\equiv 1-2r_{\mathrm{g}}/r_{\mathrm{p}},$$

$$\begin{array}{rl}\dot{t}& =\frac{1+\mathcal{E}}{1-2r_{\mathrm{g}}/r},\dot{\varphi }=\ell r_{\mathrm{g}}/r^2,\\ {\dot{r}}^2& =(1+\mathcal{E}{)}^2-\left(1-\frac{2r_{\mathrm{g}}}{r}\right)\left(1+\frac{{\ell }^2}{r^2}\right),\\ \ddot{r}& =-\frac{r_{\mathrm{g}}^2}{r^2}+\frac{{\ell }^2{r}_{\mathrm{g}}}{r^3}\left(1-\frac{3r_{\mathrm{g}}}{r}\right).\end{array}$$

$$\begin{array}{rl}{\tilde{v}}_{\mathrm{r}}& =\frac{\mathrm{d}\tilde{r}}{\mathrm{d}\tilde{t}}={\mu }_{\mathrm{I}}^{-1}\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{\dot{r}(r_{\mathrm{I}})}{1+\mathcal{E}},\\ {\tilde{v}}_{\varphi }& =r_{\mathrm{I}}\frac{\mathrm{d}\tilde{\varphi }}{\mathrm{d}\tilde{t}}=\frac{r_{\mathrm{I}}{\mu }_{\mathrm{I}}^{1/2}}{1+\mathcal{E}}\dot{\varphi }(r_{\mathrm{I}}).\end{array}$$

$$\tan {\tilde{\theta }}_{\mathrm{I}}=\frac{{\tilde{v}}_{\mathrm{r}}}{{\tilde{v}}_{\varphi }}=\frac{\dot{r}(r_{\mathrm{I}})}{{\mu }_{\mathrm{I}}^{1/2}{r}_{\mathrm{I}}\dot{\varphi }(r_{\mathrm{I}})}.$$

$$f_{\mathrm{diss}}=\left(\frac{{\mu }_{\mathrm{I}}^{1/2}}{\sqrt{1-{\tilde{v}}_{\mathrm{r}}^2-{\tilde{v}}_{\varphi }^2}}-\frac{{\mu }_{\mathrm{I}}^{1/2}}{\sqrt{1-{\tilde{v}}_{\varphi }^2}}\right)(1-{\mu }_{\mathrm{I}}^{1/2}{)}^{-1},$$

$${\bar{v}}_{\mathrm{c}}=\frac{\tilde{\Gamma }{\tilde{v}}_r}{\sqrt{1+{\tilde{\Gamma }}^2{\tilde{v}}_r^2}}.$$

$$\bar{v}(\bar{\theta })={\bar{v}}_c,$$

$$\begin{array}{rl}{\bar{\rho }}_{\mathrm{n}}& =\sum _{i=0}^4{q}_i{x}^i,q_0=1.051\times 10^{-2},q_1=0.1103,\\ q_2& =-0.2017,q_3=0.2434,q_4=-0.08436,\end{array}$$

$$\begin{array}{rl}{\tilde{u}}_t& ={\tilde{\Gamma }}_{\varphi }{\bar{u}}_t+{\tilde{\Gamma }}_{\varphi }{\tilde{v}}_{\varphi }{\bar{u}}_x,{\tilde{u}}_y={\bar{u}}_y,\\ {\tilde{u}}_x& ={\tilde{\Gamma }}_{\varphi }{\tilde{v}}_{\varphi }{\bar{u}}_t+{\tilde{\Gamma }}_{\varphi }{\bar{u}}_x,{\tilde{u}}_z={\bar{u}}_z.\end{array}$$

$$\ell (\bar{\theta },\bar{\varphi })=r_{\mathrm{I}}\sqrt{{\tilde{u}}_x^2+{\tilde{u}}_y^2},1+\mathcal{E}(\bar{\theta },\bar{\varphi })={\mu }_{\mathrm{I}}^{1/2}{\tilde{u}}_t,$$

$$M_{\mathrm{cr}}\simeq (4.6\times 10^6{\mathrm{M}}_{\odot }){\eta }^{0.4}{\beta }^{-1.7}{m}_{\ast }^{-1/2}{r}_{\ast }^{3/2}.$$

$$r_{\mathrm{p},\mathrm{cr}}/r_{\mathrm{g}}\simeq 17{\eta }^{-4/15}{\beta }^{2/15},$$

$$M_{\max }\simeq (4.0\times 10^7{\mathrm{M}}_{\odot }){\beta }^{-3/2}{m}_{\ast }^{-1/2}{r}_{\ast }^{3/2}.$$

$$r_{\mathrm{p}}\approx r_{\mathrm{mb}}\approx 5^{1/3}{\xi }_{\ast }{r}_{\mathrm{T}},$$

$$f_{\mathrm{diss}}=1-\frac{1}{1+x},\mathrm{where}x=0.27{\eta }^{-1}{\beta }^3{M}_6^{5/3}{m}_{\ast }^{1/3}{r}_{\ast }^{-2}.$$

$$\Delta {\mathcal{E}}_{\mathrm{vis},\max }\sim \alpha v_{\mathrm{A}}^2(r_{\mathrm{T}})\sqrt{\frac{r_{\mathrm{a}}}{r_{\mathrm{T}}}}.$$

$${\dot{M}}_{\mathrm{fb},\max }\simeq (3{\mathrm{M}}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}){\eta }_{\max }^{3/2}{M}_6^{-1/2}{m}_{\ast }^2{r}_{\ast }^{-3/2},$$

$${\left(\frac{{\lambda }_{\mathrm{Edd}}{L}_{\mathrm{Edd}}}{4\text{uppi}{r}_{\mathrm{T}}^2{\sigma }_{\mathrm{SB}}}\right)}^{\frac{1}{4}}\simeq 2.6\times 10^5\mathrm{K}\frac{{\lambda }_{\mathrm{Edd}}^{1/4}{M}_6^{1/12}{m}_{\ast }^{1/6}}{r_{\ast }^{1/2}},$$

$$\begin{array}{rl}\rho (r)& \simeq \frac{{\dot{M}}_{\max }{f}_{\mathrm{unb}}}{\Omega r^{2}v}\simeq 1.5\times 10^{-12}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}\frac{K}{r_{14}^2},\\ K& \equiv \frac{2\text{uppi}}{\Omega }\frac{f_{\mathrm{unb}}}{0.5}\frac{{\eta }_{\max }^{3/2}{m}_{\ast }^2}{M_6^{1/2}{r}_{\ast }^{3/2}{v}_9},\end{array}$$

$$r_{\mathrm{tr}}\simeq 1.7\times 10^{14}K{v}_9\mathrm{cm},$$

$$U(r)=8.0\times 10^5\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}\frac{L_{\mathrm{in},44}{r}_{\mathrm{in},14}^{2/3}}{r_{14}^{8/3}{v}_9}\min \left[1,{\left(\frac{r_{\mathrm{tr}}}{r}\right)}^{1/3}\right].$$

$$T(r)\simeq 1.0\times 10^5\mathrm{K}\frac{L_{\mathrm{in},44}^{1/4}{r}_{\mathrm{in},14}^{1/6}}{r_{14}^{2/3}{v}_9^{1/4}}(\mathrm{for}r<r_{\mathrm{tr}}).$$

$$\eta (r)\equiv \frac{U(r)c/4\text{uppi}}{\int \mathrm{d}\nu \min \left[B_{\nu }(T_{\mathrm{BB}}),j_{\nu }(T_{\mathrm{BB}})c{t}_{\mathrm{diff}}\right]},$$

$$\eta (r_{\mathrm{th}})={\eta }_{\mathrm{crit}}.$$

$$\sigma =\max (\tau v/c,v/\sqrt{2}c),$$

$$j_{\nu }(r_{\mathrm{th},\nu })c{t}_{\mathrm{diff}}=B_{\nu }(T_{\mathrm{BB}}(r_{\mathrm{th},\nu })),$$

$$L_{\nu }\simeq 4\text{uppi}{r}_{\mathrm{th},\nu }^{3}4\text{uppi}{j}_{\nu }(r_{\mathrm{th},\nu })\simeq 4\text{uppi}{r}_{\mathrm{th},\nu }^2\frac{4\text{uppi}{B}_{\nu }(T_{\mathrm{BB}}(r_{\mathrm{th},\nu }))}{{\tau }_{\mathrm{s}}(r_{\mathrm{th},\nu })},$$

$${\kappa }_{\mathrm{low}\nu }\simeq (1.1\times 10^{17}\mathrm{c}{\mathrm{m}}^2{\mathrm{g}}^{-1})\rho T^{-3/2}(h\nu /\mathrm{eV}{)}^{-2},$$

Full text

1 Introduction

Tidal disruption events (TDEs) hold promise for probing the otherwise dormant supermassive blackholes (BHs) at the centers of most galaxies (Rees, 1988). The story starts with simple initial conditions: a star, of certain mass and radius, approaches the BH on a parabolic orbit of certain specific angular momentum. The star can be treated as a point mass until it reaches the tidal radius where the tidal forces exceed the star’s self-gravity. The hydrodynamical disruption phase, despite its complexity, is understood to at least order-unity level, thanks to decades of analytical and numerical studies (e.g., Lacy et al., 1982; Carter & Luminet, 1983; Rees, 1988; Evans & Kochanek, 1989; Laguna et al., 1993; Ayal et al., 2000; Lodato et al., 2009; Stone et al., 2013; Guillochon & Ramirez-Ruiz, 2013; Tejeda et al., 2017; Goicovic et al., 2019; Steinberg et al., 2019; Gafton & Rosswog, 2019). The result is that the post-disruption stellar debris acquires a spread of specific orbital energy, which is roughly given by the gradient of the BH’s gravitational potential across the star at the tidal radius. This means that roughly half of the stellar debris is unbound and the other half is left in highly eccentric bound orbits.

After the disruption phase, the star is tidally stretched into a very long thin stream and the evolution of the stream structure in the transverse and longitudinal directions are decoupled (Kochanek, 1994). Thus, the system enters the free-fall phase where each stream segment follows its own geodesic like a test particle (Coughlin et al., 2016). Then, after passing the apocenters of the highly eccentric orbits, the bound debris falls back towards the BH at a rate given by the distribution of specific energy (Evans & Kochanek, 1989; Phinney, 1989). Due to relativistic apsidal precession, the bound debris, after passing the pericenter, collides violently with the still in-falling stream (see Fig. 1). It has been shown that shocks at the self-intersection point is the main cause of orbital energy dissipation and the subsequent formation of an accretion disk (Rees, 1988; Kochanek, 1994; Hayasaki et al., 2013; Guillochon et al., 2014; Shiokawa et al., 2015; Bonnerot et al., 2016). However, the aftermath of the self-intersection is an extremely complex problem, which depends on the interplay among magnetohydrodynamics, radiation, and general relativity in 3D. No numerical simulations to date have been able to provide a deterministic model for TDEs with realistic star-to-BH mass ratio and high eccentricity (see Stone et al., 2018a, for a review). Many simulations consider either an intermediate-mass BH (e.g. Guillochon et al., 2014; Evans et al., 2015; Shiokawa et al., 2015; Sa̧dowski et al., 2016) or the disruption of a low-eccentricity (initially bound) star (e.g. Bonnerot et al., 2016; Hayasaki et al., 2016). It is unclear how to extrapolate the simulation results to realistic configurations and provide an answer to the following questions: How long does it take for the bound gas to form a circular accretion disk (if at all)? How much radiative energy is released from the system? What fraction of the radiation is emitted in the optical, UV or X-ray bands?

The hope lies in the rapidly growing sample of TDE candidates discovered by recent UV-optical surveys, such as GALEX (Gezari et al., 2008, 2009), SDSS (van Velzen et al., 2011), Pan-STARRS (Gezari et al., 2012; Chornock et al., 2014; Blanchard et al., 2017), PTF (Arcavi et al., 2014; Blagorodnova et al., 2017; Hung et al., 2018), ASAS-SN (Holoien et al., 2014, 2016), and ZTF (van Velzen et al., 2018a), see the open TDE catalog http://tde.space. These events have highly diverse properties in terms of peak optical luminosities, lightcurve shapes, emission line profiles, and optical/X-ray flux ratios. Still, they provide a number of important clues for understanding the dynamics of UV-optical selected TDEs: (1) the photospheric radius of the (thermal) optical emission is typically \sim$$10^{14}–$10^{15}\rm\,cm$; (2) the typical widths of H$\alpha$ and/or He II emission lines in the optical band and CIV, NV, SiIV aborption lines in the UV band (e.g. Blagorodnova et al., 2018) are of order \sim$$10^{4}\rm\,km/s; (3) the rise/fade timescale is of order $\sim$months111We note a few exceptions such as iPTF16fnl (Blagorodnova et al., 2017) and ASASSN-15lh (Dong et al., 2016; Leloudas et al., 2016). We also note that current optical surveys are biased against detecting very fast ($\lesssim$week) and very long ($\gtrsim$year) transients, so the rise/fade timescales of detected events may not representative for the entire TDE family.; (4) the total energy radiated in the UV-optical band is typically \lesssim$$10^{51}\rm\,erg, which is much smaller than the energy budget of the system (\gtrsim$$10^{53}\rm\,erg even for disruption of low-mass stars).

The photospheric radius is much larger than the tidal radius (of order \sim$$10^{13}\rm\,cm), and the velocity inferred from line widths is much smaller than the Keplerian/escape velocity near the tidal radius. These properties are inconsistent with the wind-reprocessed emission from a circularized accretion disk near the tidal radius (Strubbe & Quataert, 2009; Miller, 2015). The low radiative efficiency in the optical band is known as the “missing energy” puzzle (Piran et al., 2015; Stone & Metzger, 2016; Lu & Kumar, 2018), whose solution depends on the source of the optical emission. Based on the arguments that the photospheric radius is of the same order as the semimajor axis of the most bound orbit and that the line width roughly agrees with the Keplerian velocity at the same radius, Piran et al. (2015) proposed that the optical emission is powered by the dissipation of orbital energy by stream self-intersection. An alternative phonomelogical model proposed by Metzger & Stone (2016) is that only a small fraction $f_{\rm in}\ll 1$ of the fall-back gas actually accretes onto the BH and the rest $(1-f_{\rm in})$ is blown away by the gravitational energy released from the accreting gas. In this model, if the energy efficiency of accreting gas is $\eta_{\rm acc}=0.1\eta_{\rm acc,-1}$, then the accretion fraction of order $f_{\rm in}\sim 10^{-2}\eta_{\rm acc,-1}^{-1}$. However, these models do not consider the detailed dynamics of the stream self-intersection and disk formation.

In this paper, we consider the stream kinematics in a semi-analytical way and explore the diverse consequences of the stream self-intersection. This approach is similar to Dai et al. (2015) who studied the location and gas velocity at the self-intersection point in a post-Newtonian way (only considering the lowest-order apsidal precession). However, we evolve the system in full general relativity before and after the self-intersection and study the properties of the shocked gas that are unbound, accreting, and plunging. More importantly, instead of assuming inelastic collision as in Dai et al. (2015), we use the realistic equation of state for radiation-dominated gas to model the intersection, motivated by the local simulation of colliding streams by Jiang et al. (2016). Thus, our approach provides a more comprehensive and self-consistent picture of the dynamics and multiwavelength emission from TDEs.

This paper is organized as follows. In §2, we calculate the location of the self-intersection point and the velocities of the two streams before the collision. In §3, we perform hydrodynamical simulation of the collision process. In §4, we consider the fate of the shocked gas after the self-intersection. Implications of TDE dynamics on the multiwavelength observations will be considered in §5. We discuss a number of issues in our modeling in §6. A summary is provided in §7. Unless otherwise specified, we use geometrical units where the gravitational constant and speed of light are $G=c=1$.

2 Self-Intersection of the Fallback Stream

We consider a star of mass $M_{*}=m_{*}\mathrm{M_{\odot}}$ and radius $R_{*}=r_{*}R_{\odot}$ interacting with a BH of mass $M=10^{6}M_{6}\mathrm{M_{\odot}}$. The gravitational radius of the BH is $r_{\rm g}\equiv M$. We take the pericenter of the star’s initial orbit to be $r_{\rm p}=r_{\rm T}/\beta$, where $\beta$ is a free impact parameter describing the depth of penetration and the $r_{\rm T}$ is the Newtonian Roche tidal radius defined as (Hills, 1975)

[TABLE]

The lower limit of the impact parameter $\beta_{\rm min}$ is of order unity, but to obtain its exact value corresponding to marginal disruption, one must take into account relativistic tidal forces and realistic stellar structure/rotation (these will be discussed later in §5.3). After the disruption, the stellar debris attains a spread of specific orbital energy for the stellar debris $\mathcal{E}\in(-\eta_{\rm max}\mathcal{E}_{\rm T},+\eta_{\rm max}\mathcal{E}_{\rm T})$, where we have defined the Newtonian tidal energy

[TABLE]

and $\eta_{\rm max}$ is a constant of order unity containing the uncertainties due to stellar structure/rotation and the detailed relativistic disruption process. The Newtonian orbital period of the leading edge ($\eta=\eta_{\rm max}$) is $P_{\rm min}\simeq(41\mathrm{\,d})\,\eta_{\rm max}^{-3/2}M_{6}^{1/2}m_{*}^{-1}r_{*}^{3/2}$.

The bound materials corresponding to $\mathcal{E}<0$ form an elongated thin stream which collides with itself due to apsidal precession (Fig. 1). Since the width of the stream is much smaller than the pericenter radius (e.g., Kochanek, 1994; Coughlin et al., 2016), a given stream segment, characterized by its specific energy $\mathcal{E}=\eta\mathcal{E}_{\rm T}$ and pericenter radius $r_{\rm p}=r_{\rm T}/\beta$ ($\eta\leq\eta_{\rm max}$ and $\beta\gtrsim 1$ are free parameters), moves along a certain geodesic until it collides with the still in-falling gas. Note that we define $\mathcal{E}$ and $r_{\rm p}$ based on Newtonian quantities $\mathcal{E}_{\rm T}$ and $r_{\rm T}$ only for convenience reason, our treatment of the orbital kinematics is fully general relativistic.

In this paper, we consider the simplest case of a non-spinning BH (the effects of BH spin will be discussed in §6). In spherical coordinates for the Schwarzschild spacetime, the initial position of the stream segment is $(t=0,\ r=r_{\rm p},\ \theta=0,\ \phi=0)$ and the proper time of the stream segment starts as $\tau=0$. We align the orbital plane with the equatorial plane of the coordinate system, so $\dot{\theta}\equiv\mathrm{d}\theta/\mathrm{d}\tau=0$. The specific angular momentum of is given by

[TABLE]

where $1+\mathcal{E}$ is the total energy including rest mass. Hereafter, the time derivative of any quantity $Q$ with respect to the stream segment’s proper time $\tau$ is denoted by $\dot{Q}$. Measuring the proper time in units of $r_{\rm g}$, we write the geodesic equations

[TABLE]

We use a Leapfrog method to integrate the above geodesic equations with timestep $\delta\tau=r_{\rm g}/30$. Since these two colliding flows have similar specific energies $\eta_{1}\approx\eta_{2}$, the radius for self-intersection $r_{\rm I}$ is approximately given by $\phi(r_{\rm I})=\uppi$.

For a stationary observer at the intersecting point $r_{\rm I}$, we define $\mu_{\rm I}\equiv 1-2r_{\rm g}/r_{\rm I}$, so the local differential length in the radial direction is $\mathrm{d}\widetilde{r}=\mu_{\rm I}^{-1/2}\mathrm{d}r$ and the local differential time is $\mathrm{d}\widetilde{t}=\mu_{\rm I}^{1/2}\mathrm{d}t=\mu_{\rm I}^{-1/2}(1+\mathcal{E})\mathrm{d}\tau$. In the following, we consider the stream-intersection process in the comoving frame of a local stationary observer at radius $r_{\rm I}$ (LSO frame hereafter), in which any quantity $Q$ is denoted with a tilde $\widetilde{Q}$. Then, the radial and transverse velocities of the colliding streams in the LSO frame are

[TABLE]

The intersecting half angle $\widetilde{\theta}_{\rm I}$ in the LSO frame is given by

[TABLE]

In Fig. 2, we compare the self-intersection radius, angle, and velocities from our general relativistic calculations with the corresponding lowest-order post-Newtonian results by Dai et al. (2015). We take $\beta=1.0$ and $\eta=1.0$ to be our fiducial parameters. We consider four different stellar masses of $m_{*}=0.2$, $0.5$, $1.0$, $1.5$ and the corresponding zero-age main-sequence stellar radii $r_{*}=0.23$, $0.46$, $0.89$, $1.63$ are taken from Tout et al. (1996) assuming solar metallicity, with errors of a few percent. As expected, we find that, for more massive BHs, the self-intersection occurs closer to the event horizon and the intersecting velocity is larger (the interaction is more violent).

If one assumes that two colliding flows have equal cross-sections and that the collision is completely inelastic, then the radial 4-velocity component $\dot{r}(r_{\rm I})$ gets dissipated and that the transverse 4-velocity component $r_{\rm I}\dot{\phi}(r_{\rm I})$ survives. In this case, we can quantify the efficiency of orbital energy dissipation by defining

[TABLE]

which describes the change in orbital energy divided by the gravitational binding energy at radius $r_{\rm I}$. This (maximum possible) dissipation efficiency is shown in the fourth panel of Fig. 2. In the low BH mass limit $M_{6}\ll 1$, the dissipation of orbital energy by shocks is extremely weak and we asymptotically have $f_{\rm diss}\propto M^{5/3}$ (marked as a black dashed line in the fourth panel), which agrees with the result of Bonnerot et al. (2017a). In those cases, if the circularization is still dominated by stream self-intersection, then the orbit stays highly eccentric for roughly $f_{\rm diss}^{-1}$ rounds and hence the circularization timescale is roughly $f_{\rm diss}^{-1}P_{\rm min}\propto M^{-7/6}$ (since $P_{\rm min}\propto M^{1/2}$). Other mechanisms, e.g. the magneto-rotational instability, may cause angular momentum exchange and drive circularization on a shorter timescale (Chan et al., 2018). As we discuss later in §5.1, TDEs by low-mass BHs typically generates long-lasting eccentric accretion disks which produce long-duration transients. On the other hand, for high-mass BHs $M_{6}\gtrsim 1$, stream intersection causes strong dissipation of orbital energy and hence the orbit may quickly circularize.

In the next section, we show that completely inelastic collision, as assumed by e.g. Dai et al. (2015), is a poor description of the stream dynamics, because the shocked gas is highly optically thick and hence evolves in a nearly adiabatic manner (Jiang et al., 2016).

3 Hydrodynamical Simulations of the Self-Intersecting Shocks

We numerically simulate the stream-stream collision in a special inertia frame described as follows. In the LSO frame, the 4-velocity of the outward-moving stream is $(\widetilde{u}_{\alpha})=(\widetilde{u}_{t},\widetilde{u}_{r},\widetilde{u}_{\theta},\widetilde{u}_{\phi})=\widetilde{\Gamma}(1,\widetilde{v}_{r},0,\widetilde{v}_{\phi})$, where the Lorentz factor is $\widetilde{\Gamma}\equiv(1-\widetilde{v}^{2})^{-1/2}$ and $\widetilde{v}^{2}=\widetilde{v}_{\rm r}^{2}+\widetilde{v}_{\phi}^{2}$. Our simulation box is centered at the self-intersecting point and is moving at velocity $\widetilde{v}_{\phi}$ with respect to the local stationary observer in the $\hat{\phi}$ direction. Thus, in the comoving frame of the simulation box (hereafter the SB frame), the two streams collide head-on with each of them moving at 4-velocity $\bar{u}_{\rm c}=\widetilde{\Gamma}\widetilde{v}_{r}$ and Lorentz factor $\bar{\Gamma}_{\rm c}=(1+\widetilde{\Gamma}^{2}\widetilde{v}_{\rm r}^{2})^{1/2}$, which means

[TABLE]

Hereafter, any quantity $Q$ in the SB frame is denoted with an overhead bar $\bar{Q}$. In all possible cases, the incoming velocities of the streams in the SB frame are sub-relativistic $\bar{v}_{\rm c}<0.3c$ (see Fig. 2). Due to extreme stretch and adiabatic cooling, the initial streams are dynamically cold with sound speed much less than the bulk velocity $\bar{v}_{\rm c}$. Another property of the initial streams is that the transverse size is much less than the orbital size \sim$$r_{\rm I} (Kochanek, 1994; Coughlin et al., 2016; Bonnerot et al., 2017a). These properties of the problem enable us to use a single non-relativistic hydrodynamic simulation in a flat spacetime to capture the structure of the shocked gas, which is self-similar within a region of size $\ll r_{\rm I}$.

In our simulation, we use an adiabatic ideal gas equation of state $P\propto\rho^{4/3}$. This is motivated by: (1) the high-density shocked gas is radiation pressure dominated, and (2) the shocked gas is highly optically thick before most of the heat is converted back to bulk motion via $PdV$ work (Jiang et al., 2016). The radiative efficiency of the shocked gas is estimated in Appendix B. Since we are concerned with the fate of the majority of the gas with fallback time $\lesssim 10P_{\rm min}$, the adiabatic assumption is a good one. For simplicity, we assume that the two colliding streams have the same cross-section and that there is no offset in the transverse direction. We will discuss the validity of these assumptions and consequences of relaxing them in §6.

We perform the simulation with the non-relativistic hydrodynamics module of PLUTO (Mignone et al., 2007), solving the mass and momentum conservation equations in 2D cylindrical coordinates ($R,\,z$). The $R$-axis corresponds to the $\hat{\phi}$ direction in the BH rest frame, and the $z$-axis is parallel to the $\hat{r}$ direction in the BH rest frame. The size of our simulation box is $0\leq R\leq 320$ and $-320\leq z\leq 320$. Two identical steady streams are injected in the form of top-hat jets moving in opposite directions at $z=-320$ and $z=320$ in the radius range $0\leq R\leq 1$. The other boundary conditions are as follows: $R=0$ axis-symmetric, $R=320$ outflow, $z=-320$ and $z=320$ outflow (except for the inner cylinder $R\leq 1$ where the streams are injected). The resolution222We also ran the same simulation with lower resolution $\delta R=\delta z=0.25$, and found the results to be similar. is $\delta R=\delta z=0.125$ ($N_{\rm R}=N_{\rm z}/2=2560$), which means the initial stream is resolved by 8 cells in the transverse direction. The initial streams have mass density 1, pressure $(4/3)^{-1}$, and velocities $\pm 3000$ (all in machine units, since the problem is scale-free in the non-relativistic limit). Since the adiabatic sound speed of the stream is unity, the Mach number is $3000\gg 1$ and hence the streams are effectively cold. The pressure of the ambient medium matches that of the streams. The mass density of the ambient medium is extremely small $10^{-8}$, so the shocked gas expands as if in vacuum. We run the simulation with time step $\delta t\simeq 6\times 10^{-6}$ for a sufficiently long time $t=0.8$ (or 7.5 times the domain crossing time) so that the structure of the shocked gas within a sphere of radius 320 has relaxed to a nearly stable configuration.

The large-scale structure of the system at $t=0.8$ is shown in Fig. 3. We see that the two streams collide at $z=0$ and the shocked gas expands in a roughly spherical way to radii much larger than the stream width (which equals to unity). In Fig. 4, we show the radial profiles of the mass flux at different polar angles $\bar{\theta}=30^{\rm o},\ 45^{\rm o},\ 60^{\rm o}$, and $90^{\rm o}$. The angular profiles of the velocity and mass flux at three different radii $r=75$, $150$ and $300$ are shown in Fig. 5. We can see that, at large distances from the shocks $r\gtrsim 100$, the velocity profile is very flat but the mass flux is heavily concentrated near the equatorial plane $60^{\rm o}\lesssim\bar{\theta}\leq 90^{\rm o}$. In the following, we simplify the velocity angular profile as isotropic

[TABLE]

and hence the density angular profile is the same as the mass flux profile. We use a fourth-order polynomial fit to the normalized density angular profile given by

[TABLE]

where $x\equiv\mathrm{min}(\bar{\theta},\,\uppi-\bar{\theta})$ and $2\uppi\int_{0}^{\uppi}\bar{\rho}_{\rm n}\sin\bar{\theta}\mathrm{d}\bar{\theta}=1$.

In the following, we Lorentz transform the velocity and mass flux angular profiles of the shocked gas from the SB frame back to the LSO frame. For a fluid element moving with speed $\bar{v}_{c}$ in an arbitrary ($\bar{\theta},\bar{\phi}$) direction ($\bar{\theta}$ being the polar angle and $\bar{\phi}$ being the azimuthal angle), we write its four-velocity in Cartesian components $(\bar{u}_{\alpha})=\bar{\gamma}(1,\bar{v}_{x},\bar{v}_{y},\bar{v}_{z})$, where $\bar{\gamma}=(1-\bar{v}^{2}_{c})^{-1/2}$, $\bar{v}_{x}=\bar{v}_{c}\sin\bar{\theta}\cos\bar{\phi}$, $\bar{v}_{y}=\bar{v}_{c}\sin\bar{\theta}\sin\bar{\phi}$, and $\bar{v}_{z}=\bar{v}_{c}\cos\bar{\theta}$. The simulation box is moving at velocity $\widetilde{v}_{\phi}$ and the corresponding Lorentz factor is $\widetilde{\Gamma}_{\phi}\equiv(1-\widetilde{v}_{\phi}^{2})^{-1/2}$, so the 4-velocity in the LSO frame is

[TABLE]

Then, the specific angular momentum and specific energy of this fluid element in the Schwarzschild spacetime are

[TABLE]

In the next section, we discuss what fraction of the shocked gas is unbound, plunging, or accreting.

4 Fate of the shocked gas after the self-intersection

When the shocked gas expands to a distance much greater than the stream width, the internal pressure becomes low enough that the motion of individual fluid elements is approximately ballistic. If the geodesic reaches infinity or inside the event horizon, we call the fluid element “unbound” or “plunging”, respectively. Those fluid elements with bound but non-plunging geodesics are denoted as “accreting.” The geodesic of a fluid element moving in the ($\bar{\theta},\bar{\phi}$) direction at $r_{\rm I}$ has specific angular momentum $\ell(\bar{\theta},\bar{\phi})$ and specific energy $1+\mathcal{E}(\bar{\theta},\bar{\phi})$, which are given by eq. (12). We note that the marginally bound parabolic orbit for the Schwarzschild spacetime has specific angular momentum $\ell_{\rm mb}=4r_{\rm g}$ and pericenter radius $r_{\rm mb}=4r_{\rm g}$ (Bardeen et al., 1972), so the stream self-intersection radius must always be greater than $4r_{\rm g}$ (see Fig. 2).

In Fig. 6 and 7, we show the Mollweide projection map of fate in terms of the polar angle $\bar{\theta}$ and azimuthal angle $\bar{\phi}$ in the simulation box frame. Here the polar angle $\bar{\theta}=0$ and $\uppi$ correspond to latitudes $90^{\rm o}$ N and $90^{\rm o}$ S, respectively. The azimuthal angle $\bar{\phi}=0$ and $\uppi$ correspond to longitudes $90^{\rm o}$ E and $90^{\rm o}$ W. The unbound (“unb”), accreting (“acc”), and plunging (“plg”) regions are shown in light grey, dark grey, and black, respectively. For the four cases with different BH masses $M_{6}=1.7,\ 2.5,\ 5,$ and $10$, we fix the star’s mass $m_{*}=0.5$, impact parameter $\beta=1.0$ ($r_{\rm p}=r_{\rm T}/\beta$) and orbital energy parameter $\eta=1.0$ ($\mathcal{E}=\eta\mathcal{E}_{\rm T}$). The cyan contours show the distribution of specific angular momentum (in units of $r_{\rm g}$) projected in the direction of the star’s initial orbital angular momentum. Note that, when determining the fate of a certain fluid element, we take into account its total specific energy and total specific angular momentum. Subsequently, the out-of-plane component of the angular moment will be further damped by shocks within the “accreting” gas. If cooling is efficient, then more gas is expected to plunge into the horizon. We also note that not all fluid elements marked as “plunging” will necessarily fall directly into the horizon. For instance, those moving in the ($\bar{\theta}\sim 0$, $\bar{\phi}\sim\uppi/2$) direction (near the north pole of Fig. 6 and 7) will most likely run into the “accreting” gas. The detailed dynamical evolution of the bound gas is studied in a separate work (Bonnerot & Lu, 2019).

In Fig. 8, we show the mass-weighted distributions of specific angular momentum projected along the star’s initial angular momentum and specific orbital energy, for four cases with $M_{6}=1.5$, 2.5, 5, and 10 (while keeping $\beta=1.0$, $\eta=1.0$, and $m_{*}=0.5$ fixed). Before the collision, the stream has specific angular momentum $\ell_{0}\approx\ell_{*}=r_{\rm p}\sqrt{1/\mu_{\rm p}-1}$ and orbital energy $-\mathcal{E}_{\rm T}$. The collision causes a spread of specific angular momentum by $\Delta\ell/\ell_{*}\sim$ a few, and the corresponding spread in the Keplerian circularization radius is about a factor of 10. In some cases (e.g. $M_{6}=1.5$ and 2.5), a large fraction of shocked gas is in counter-rotating orbits ($\ell<0$) and will subsequently collide with the forward-rotating gas ($\ell>0$) at a wide range of radii. The spread in specific orbital energy after the collision is very sensitive to the BH mass, with $\Delta\mathcal{E}/\mathcal{E}_{\rm T}\sim 1$ for $M_{6}=1.5$ but $\Delta\mathcal{E}/\mathcal{E}_{\rm T}\sim 100$ for $M_{6}=10$. For the $M_{6}=1.5$ case, there is no unbound gas. For BH mass $M_{6}\gtrsim 2$, a large fraction of the shocked gas is unbound ($\mathcal{E}>0$). For a highly eccentric Keplerian orbit, the eccentricity is given by $e=\sqrt{1-2|\mathcal{E}|\ell^{2}/r_{\rm g}^{2}}\approx 1-|\mathcal{E}|\ell^{2}/r_{\rm g}^{2}$. We can see that strong shocks due to self-intersection increase the product $|\mathcal{E}|\ell^{2}$ by one order of magnitude or more, and hence the accreting fraction of gas should quickly circularize.

In Fig. 9, we show the mass fractions of the unbound, accreting, and plunging gas as a function of BH mass, for four stellar masses $m_{*}=0.2$, $0.5$, $1.0$, and $1.5$ (while keeping $\beta=1.0$ and $\eta=1.0$ fixed). We find that, above a critical BH mass (to be quantified shortly), the unbound fraction quickly rises from $0\%$ to a maximum of $50\%$. At the same time, the accreting fraction drops from $100\%$ to $50\%$. For higher mass BHs, the plunging fraction333Note that the plunging fraction is non-zero even for very low BH masses, this is because the collision broadens the angular momentum distribution such that part of the shocked gas has almost zero angular momentum (see the upper panel of Fig. 8). The small bump (or dip) in the plunging (accreting) fraction for the $m_{*}=1.5$ case near BH masses $M_{6}\sim 0.3$ is because the velocity before the collision has comparable $\hat{r}$ and $\hat{\phi}$ components: $v_{r}\sim v_{\phi}$ (see the third panel of Fig. 2). We show the map of fate for $m_{*}=1.5$ and $M_{6}=0.3$ in Fig. 20 in the Appendix. quickly rises at the expense of the dropping accreting fraction, while the unbound fraction stays roughly unchanged at \sim$$50\%. There is a maximum mass $M_{\rm max}$ above which no accretion is possible, because the entire star plunges into the event horizon.

In this paper, we call the unbound fraction of the shocked gas the “collision-induced outflow” (CIO), which has important observational consequences (see the next section). The launching of CIO can be understood in the following Newtonian picture. If the intersection occurs at $r_{\rm p}\ll r_{\rm I}\lesssim r_{\rm a}/3$ ($r_{\rm a}$ being the apocenter radius), then the two streams typically collide at a large angle $2\theta_{\rm I}\gtrsim 90^{\rm o}$ near the local escape speed $|v_{\rm r}\hat{r}+v_{\rm\phi}\hat{\phi}|=\sqrt{v_{\rm r}^{2}+v_{\rm\phi}^{2}}\simeq v_{\rm esc}\simeq\sqrt{2r_{\rm g}/r_{\rm I}}$. The radial velocity component is dissipated by shocks, and then the shocked gas adiabatically expands at speed \sim$$v_{\rm r} in a roughly spherical manner in the SB frame moving at velocity $v_{\rm\phi}\hat{\phi}$. Going back to the LSO frame, we find the fastest moving shocked gas in the ($\bar{\theta}=\uppi/2$, $\bar{\phi}=0$) direction with velocity $(v_{\rm r}+v_{\rm\phi})\hat{\phi}$ and speed $v_{\rm r}+v_{\rm\phi}>\sqrt{v_{\rm r}^{2}+v_{\rm\phi}^{2}}\simeq v_{\rm esc}$. We see that CIO is a generic feature of gas streams colliding near the local escape speed of the intersection point444This feature was captured in the simulations of deeply penetrating TDEs by Evans et al. (2015); Sa̧dowski et al. (2016); Jiang et al. (2016). In many other works, the gas streams do not collide near the local escape speed of the collision point, either because the aspidal precession is so weak (due to small $M/M_{*}$ ratio) that the collision occurs near the apocenter (e.g. Shiokawa et al., 2015) or because the star is initially in a bound orbit with too low eccentricity (e.g. Bonnerot et al., 2016; Hayasaki et al., 2016). . In Fig. 10, we show the asymptotic kinetic energy and mass-weighted mean speed of the CIO for a number of cases, assuming the total amount of unbound mass to be $f_{\rm unb}m_{*}\mathrm{M_{\odot}}$.

We define the critical BH mass $M_{\rm cr}$ above which the mass fraction of unbound gas is more than $20\%$, i.e. $f_{\rm unb}(M_{\rm cr})=20\%$. After exploring an extensive grid of parameters (see Fig. 17 in the Appendix), we find

[TABLE]

This can be translated to a critical pericenter radius

[TABLE]

below which $f_{\rm unb}\gtrsim 20\%$. If we choose the critical unbound fraction to be $30\%$ (instead of $20\%$), the scalings in the above equations stay the same but the normalization changes to $M_{\rm cr}\simeq 5.3\times 10^{6}\,\mathrm{M_{\odot}}$ (and $r_{\rm p,cr}/r_{\rm g}\simeq 15$). The precise value of the critical unbound fraction is unimportant, because $f_{\rm unb}$ is very sensitive to the BH mass near $M\sim M_{\rm cr}$.

The maximum Schwarzschild BH mass for tidal disruption to occur outside the event horizon can be estimated by requiring $r_{\rm p}=r_{\rm T}/\beta>r_{\rm mb}=4r_{\rm g}$ (since typically $|\mathcal{E}|\ll 1$, see eq. 2),

[TABLE]

Note that the maximum BH mass depends on the minimum impact parameter $\beta_{\rm min}$ at which the relativistic tidal forces exceeds the star’s self-gravity. In the limit $R_{*}\ll 4r_{\rm g}$, the local gravitational-field gradients can described by the relativistic tidal tensor in Fermi normal coordinates. For the Schwarzschild spacetime, the criterion for marginal tidal disruption can be written as

[TABLE]

where the $5^{1/3}$ factor comes from relativistic tidal stretching555It can be shown that, for the case where the star’s initial angular momentum is parallel to the spin of a Kerr BH, the $5^{1/3}$ factor stays the same for arbitrary spin. This is because the eigenvalues of the tidal tensor depends on the ratio $(\ell_{\rm mb}-a)/r_{\rm mb}\equiv 1$ (Kesden, 2012), where $\ell_{\rm mb}$ is the angular momentum of the marginally bound parabolic orbit and $-1<a<1$ is the dimensionless BH spin ($a<0$ corresponds to retrograde orbits). in the radial direction (Kesden, 2012) and the parameter $\xi_{*}\sim 1$ accounts for the internal structure of the star. We will discuss the choice of $\xi_{*}$ for different stellar masses in §5.3 on TDE demographics. The marginal disruption case occurs when $\ell\approx 4r_{\rm g}$ and $r_{\rm p}\approx 4r_{\rm g}$, which gives the minimum impact parameter $\beta_{\rm min}\approx 5^{-1/3}\xi_{*}^{-1}\simeq 0.6\,\xi_{*}^{-1}$ and the maximum mass for non-spinning BHs associated with TDEs $M_{\rm max,Sch}\approx(8.9\times 10^{7}\mathrm{\mathrm{M_{\odot}}})\,\xi_{*}^{3/2}m_{*}^{-1/2}r_{*}^{3/2}$. We can see that relativistic tidal forces are slightly better at disrupting stars than in the Newtonian approximation. In realistic situations, the precise $\beta_{\rm min}$ (and hence $M_{\rm max}$) will depend on the stellar structure, BH’s spin, star’s spin, and the misalignment between the star’s orbital and the BH’s spin angular momenta, etc. Fortunately, the precise value of $\beta_{\rm min}$ may not be important from the observational point of view, because those TDEs with BH mass close to $M_{\rm max}$ should be quite dim due to their low accreting fraction (most gas is either unbound or plunging, see Fig. 9).

5 Observations

In previous sections, we have described a semi-analytical model for the TDE dynamics, including the fluid properties at the stream self-intersection point and the fate of the shocked gas moving in different directions. In this section, we first discuss the circularization of the fallback stream and the formation of accretion disk in §5.1 and then we consider the observational implications of the unbound gas (when $f_{\rm unb}\neq 0$) in §5.2. TDE demographics will be discussed in §5.3.

5.1 Circularization of the fallback mass

For main-sequence stars disrupted by low-mass BHs ($M_{6}\ll 1$), the stream self-crossing occurs near the apocenter and the shocks only dissipate a small fraction of the orbital energy (see the fourth panel of Fig. 2, and Bonnerot et al., 2016; Chen & Shen, 2018). After exploring an extensive grid of parameters (see Fig. 18 in the Appendix), we find the dissipation efficiency in eq. (7) can be written in the following analytical form for $M_{6}\lesssim 1$

[TABLE]

Note that in the limit $x\ll 1$ ($f_{\rm diss}\approx x$), if the dissipation of orbital energy is only due to stream intersection, then the circularization timescale can be roughly estimated by $f_{\rm diss}^{-1}P_{\rm min}\simeq(152\mathrm{\,d})\,\eta^{-1/2}\beta^{3}M_{6}^{-7/6}m_{*}^{-4/3}r_{*}^{7/2}$. For an average star ($m_{*}\lesssim 1$) disrupted by low-mass BHs $M_{6}\ll 1$, this timescale is much longer than the durations of typical TDEs discovered in recent UV-optical surveys. We also see that tidal disruptions of red giants ($r_{*}\gg 1$) will likely have very long circularization timescales as well, unless they are in deeply penetrating orbits ($\beta\gg 1$).

If MHD turbulence develops rapidly, shear due to magnetic stresses may cause dissipation of orbital energy at a rate (per unit mass) $\dot{\mathcal{E}}_{\rm vis}\sim\alpha\Omega_{\ell}v_{\rm A}^{2}$ (Svirski et al., 2017), where $\alpha\sim 0.1$ is the viscous parameter (Shakura & Sunyaev, 1973), $\Omega_{\ell}(r)=\ell_{*}/r^{2}$ is the local orbital angular frequency, and $v_{\rm A}=\sqrt{B^{2}/4\uppi\rho^{2}}$ is the Alfvén speed. Due to conservation of flux along the stream, the magnetic field strength evolves with radius as $B\propto H(r)^{-2}$, where $H(r)$ is the stream thickness at radius $r$. Right after the disruption, the marginally bound part of the stream moves as $r\propto t^{2/3}$ and the stream length stretches666For nearly radial orbits, the Newtonian equation of motion is $r(t)=r_{0}+\int^{t}_{0}\mathrm{d}t\sqrt{2(r_{\rm g}/r+\mathcal{E})}$, where $\mathcal{E}$ is the binding energy and $r_{0}$ is the initial position. In the limit $\mathcal{E}\approx 0$ (marginally bound) and $r\gg r_{0}$, we have $r\propto t^{2/3}$. Consider two fluid elements with the same initial position $r_{0}$ but slightly different binding energy $\Delta\mathcal{E}$. After expanding for time $t$, they are separated by a distance $\Delta r=\Delta\mathcal{E}\int_{0}^{t}\mathrm{d}t[2(r_{\rm g}/r+\mathcal{E})]^{-1/2}$. In the limit $\mathcal{E}\approx 0$ and $r\gg r_{0}$, we have $\Delta r\propto tr^{1/2}\propto r^{2}$. as $r^{2}$, so we obtain the stream density evolution $\rho\propto r^{-2}H(r)^{-2}$ and hence the total viscous heating $\Delta\mathcal{E}_{\rm vis}\sim\dot{\mathcal{E}}_{\rm vis}t\propto r^{3/2}H(r)^{-2}$. The evolution of stream thickness $H(r)$ may be highly complex, depending on tidal forces, self-gravity, magnetic fields, and recombination of hydrogen (Kochanek, 1994; Guillochon et al., 2014; Coughlin et al., 2016). In the limiting case of equilibrium between self-gravity and gas pressure, we have $H\propto\sqrt{P/\rho^{2}}\propto\rho^{-1/6}$ and hence $H(r)\propto r^{1/2}$ (Coughlin et al., 2016), where we have taken a polytropic index of 5/3 which is appropriate before magnetic pressure overwhelms gas pressure or the recombination777For $\beta\sim 1$, recombination occurs at radius $r_{\rm rec}\sim\sqrt{T_{0}/10^{4}\mathrm{\,K}}\,r_{\rm T}\sim 30\,r_{\rm T}$, where $T_{0}\sim GM_{*}/kR_{*}\sim 10^{7}\rm\,K$ ($k$ being the Boltzmann constant) is the gas temperature right after tidal disruption. of hydrogen. In this regime, $\Delta\mathcal{E}_{\rm vis}\propto r^{1/2}$ and most dissipation occurs near the largest radii at which the scaling $H(r)\propto r^{1/2}$ holds. In the other limit where tidal forces dominate over self-gravity, we have $H(r)\propto r$ and $\Delta\mathcal{E}_{\rm vis}\propto r^{-1/2}$, which means that most dissipation occurs near the smallest radii at which the scaling $H(r)\propto r$ holds.

For $\beta\sim\eta\sim 1$, we provide a conservative estimate of the viscous dissipation by assuming the scaling $\Delta\mathcal{E}_{\rm vis}\propto r^{1/2}$ up to the apocenter radius of the most tightly bound orbit $r_{\rm a}$ and obtain

[TABLE]

Since $\alpha\sqrt{r_{\rm a}/r_{\rm T}}\sim\alpha(M/M_{*})^{1/6}\sim 1$, we obtain $\Delta\mathcal{E}_{\rm vis,max}\sim v_{\rm A}^{2}({r_{\rm T}})$, where $v_{\rm A}$ is the Alfvén speed near radius $r_{\rm T}$. The magnetic field may be amplified in the tidal disruption process due to forced differential rotation (Bonnerot et al., 2017b) and the total magnetic energy may be written as $f_{\rm B}GM_{*}^{2}/R_{*}$, where $GM_{*}^{2}/R_{*}$ is the work done by tidal forces and $f_{\rm B}\ll 1$ is the conversion efficiency. Then, we obtain $\Delta\mathcal{E}_{\rm vis,max}\lesssim f_{\rm B}GM_{*}/R_{*}\sim 2\times 10^{-6}f_{\rm B}\ll\mathcal{E}_{\rm T}$, where $\mathcal{E}_{\rm T}$ is the typical orbital energy of the stream (eq. 2). Therefore, the dissipation of orbital energy due to viscous shear is highly inefficient over the orbital timescale.

We also note that dissipation by the nozzle shock operating near the pericenter may also be inefficient, because the ratio between the velocity components perpendicular and inside the star’s orbital plane is of order $H/r\ll 1$. However, this picture may be changed by strong apsidal precession (which causes oblique compression) if the pericenter is close to the horizon $r_{\rm p}\lesssim 10r_{\rm g}$.

Therefore, we conclude that TDEs by low-mass BHs ($M_{6}\ll 1$) have circularization timescale $t_{\rm cir}\gg P_{\rm min}$, which is much longer than the typical duration of the current sample of TDEs discovered by UV-optical surveys. This has important implications on TDE demographics, which will be discussed in §5.3.

In the following, we focus on TDEs by relatively high-mass BHs $M_{6}\gtrsim 1$ where the dissipation of orbital kinetic energy is dominated by stream self-intersection because $f_{\rm diss}\sim 1$. As shown in Fig. 8, the distributions of specific angular momentum and orbital energy are broadened by the collision. The eccentricity $e\simeq 1-|\mathcal{E}|\ell^{2}/r_{\rm g}^{2}$ drops to the level of $1-e\sim 0.1$ due to the increase of the product $|\mathcal{E}|\ell^{2}$ by typically one order of magnitude. Subsequently, the bound gas (and a small fraction of the unbound gas, see §5.2) will collide violently at a wide range of radii between \sim$$r_{\rm p} and \sim$$r_{\rm I} over a timescale $\Omega_{\rm K}^{-1}(r_{\rm I})=c^{-1}\sqrt{r_{\rm I}^{3}/r_{\rm g}}<P_{\rm min}$, where $\Omega_{\rm K}(r_{\rm I})$ is the Keplerian angular frequency for a circular orbit at $r_{\rm I}$. Thus, orbital circularization due to exchange of energy and angular momentum occurs rapidly after the initial stream self-intersection. The detailed dynamics is highly complicated due to the interplay among gas, radiation (providing cooling), and magnetic fields (providing viscosity). This is explored in a separate work (Bonnerot & Lu, 2019).

The most interesting situation is when a significant fraction of the shocked gas becomes unbound in the form of CIO, which occurs for BH masses in the range $M_{\rm cr}<M<M_{\rm max}$. The CIO carries away (positive) energy of $\mathcal{E}\gtrsim\mathcal{E}_{\rm T}$ (see Fig. 8 and 10), and the rest of the shocked gas is left in more tightly bound orbits. The (positive) angular momentum carried away by the CIO is a factor of a few greater than that before the collision, so the remaining bound gas typically have negative angular momentum and hence rotates in the opposite direction of the initial star. Due to subsequent shocks, the counter-rotating gas will quickly settle into circular orbits not at radius $2r_{\rm p}$ but with a radial spread of at least one order of magnitude (even without viscosity).

5.2 Collision-induced Outflow (CIO)

For BH masses in the range $M_{\rm cr}<M<M_{\rm max}$, we find a large fraction of gas is launched in the form of a wide-angle CIO. In the following, we first study the morphology of the CIO (§5.2.1), and then discuss the observational implications of the CIO, including reprocessing the extreme-UV (EUV) or soft X-ray disk emission into the optical band (§5.2.2) and radio emission from the shock driven into the ambient medium (§5.2.3).

5.2.1 Morphology of the CIO

We discretize the unbound cone uniformly into $\sim$200 beams and then integrate the geodesics of each beam over longer timescales. We ignore the internal pressure of the CIO based on strong adiabatic cooling during the expansion. The CIO morphologies are shown in Figs. 11 and 12. We see that, within a distance of \sim$$10^{14}\,cm, the unbound gas expands into complex morphology and covers a large fraction of the sky viewed from the BH. A small portion of the unbound gas will collide with each other (all along the negative z axis due to BH’s gravitational focusing) and further dissipate their kinetic energy via shocks. For the low BH mass (but $M>M_{\rm cr}$ or $f_{\rm unb}\gtrsim 20\%$) cases, the self-intersection point is far from the event horizon and the gas in the unbound cone is ejected mildly above the local escape speed, so their trajectories are strongly affected by the BH’s gravity (see Fig. 11). For the high BH mass cases, the violent shocks at the intersection point launches unbound gas well above the local escape speed, so their trajectories are almost a straight line (see Fig. 12).

In the Newtonian picture (appropriate at distances $\gg r_{\rm g}$), the peak mass fallback rate of the stream can be estimated

[TABLE]

where we have assumed a flat mass distribution over specific energy between $-\eta_{\rm max}\mathcal{E}_{\rm T}$ and $\eta_{\rm max}\mathcal{E}_{\rm T}$ after tidal disruption (e.g. Evans & Kochanek, 1989; Guillochon & Ramirez-Ruiz, 2013). This peak fallback rate lasts for a duration roughly given by the period of the most bound orbit $P_{\rm min}\simeq(41\mathrm{\,d})\,\eta_{\rm max}^{-3/2}M_{6}^{1/2}m_{*}^{-1}r_{*}^{3/2}$. During this time, the time-averaged mass feeding rate to the self-intersecting point is $\dot{M}_{\rm max}\simeq\dot{M}_{\rm fb,max}$ (from both colliding streams). We note that this feeding rate is not constant but modulated by twice the free-fall timescale $2t_{\rm ff}(r_{\rm I})=2\sqrt{r_{\rm I}^{3}/r_{\rm g}}/c=(3.7\mathrm{\,d})(r_{\rm I}/10^{3}r_{\rm g})^{3/2}M_{6}$ in the Newtonian picture, because each segment of length \sim$$2r_{\rm I} will collide with the next segment of the same length. From Figs. 2 and 9, we see that for BH masses $M\gtrsim M_{\rm cr}$, the intersection radius is much below the apocenter radius $r_{\rm I}\ll r_{\rm a}$, so the modulation timescale is much less than the orbital period $2t_{\rm ff}(r_{\rm I})\ll P_{\rm min}$. This discrete mass injection may modulate the optical lightcurve during the early rise segment but not near the peak, because the CIO has highly complex structure with a broad velocity distribution such that the optical flux near the peak is contributed by multiple shells (via photon diffusion, see §5.2.2). On the other hand, if the inner accretion disk is not blocked by the large CIO column for some viewing angles, then the X-ray lightcurve may be strongly affected by the variable mass feeding rate to the accretion disk, provided that the viscous timescale is comparable or shorter than $2t_{\rm ff}(r_{\rm I})$. We also note that hydrodynamic interaction between the fallback stream and the accretion flow may modify the stream’s trajectory and cause the modulation to be non-periodic.

In the next subsection, we show that the CIO generates the optical emission from TDEs. We take an order-of-magnitude approach by assuming that the mass outflowing rate from the self-intersecting point to be steady and the unbound gas expands in a roughly spherical manner at radius $\gg r_{\rm I}$. We ignore the hydrodynamical effects of the wind/radiation from the inner accretion disk. Modeling the full radiative hydrodynamics is left for future works.

5.2.2 Optical emission from TDEs

Optically bright TDEs came as a surprise because the radiation from the inner disk has characteristic temperature

[TABLE]

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant and $\lambda_{\rm Edd}$ is the Eddington factor given by the disk bolometric luminosity over the Eddington luminosity $L_{\rm Edd}\simeq 1.5\times 10^{44}\,M_{6}\rm\,erg\,s^{-1}$ (for solar metallicity). TDEs selected as UV-optical transients have photospheric radii \sim$$10^{14}–$10^{15}\,\mathrm{cm}\gg r_{\rm T}$ and color temperatures $\sim$a few$\times 10^{4}\,\mathrm{K}$ much less than that given by eq. (20). In the following, we show that the CIO naturally provides the long-sought “reprocessing layer” which absorbs the higher frequency radiation from the inner disk and re-emits at lower frequencies (Loeb & Ulmer, 1997; Guillochon et al., 2014; Metzger & Stone, 2016).

We study the temperature structure of the CIO by assuming a steady-state spherically symmetric structure heated from the bottom at radius $r_{\rm in}\sim\,$a few$\times r_{\rm I}$. We assume the received heating power to be $L_{\rm in}$, which could be in the form of hard emission or wind from the accretion flow888The evolution of the EUV and soft X-ray luminosity from the inner accretion disk and its wind kinetic power on timescale $\lesssim 1$yr is still uncertain due to our limited understanding of multi-dimensional super-Eddington accretion flow, analytically (Begelman, 1979; Narayan & Yi, 1994; Blandford & Begelman, 2004) or numerically (Sa̧dowski et al., 2014; McKinney et al., 2014; Jiang et al., 2017). We remain agnostic about the heating source’s nature and make the (highly simplied) assumption that the velocity and density profiles of the CIO are not strongly modified by the energy injection. This assumption breaks down when the energy injection significantly accelerates the CIO, which should be studied in future works..

When the CIO reaches distances $\gg r_{\rm I}$, for a crude estimate, we assume the density and velocity distributions of the outflowing gas to be roughly uniform within a cone of solid angle $\Omega$. Then, the density profile is given by

[TABLE]

where $f_{\rm unb}$ is the unbound fraction (Fig. 9), $v=10^{9}v_{9}\rm\,cm\,s^{-1}$ is the mass-weighted mean velocity (see the third panel of Fig. 10), and $K$ is the dimensionless “wind constant” which depends on many parameters.

The photon-trapping radius $r_{\rm tr}$, where photon diffusion time equals to the dynamical expansion time, is given by the scattering optical depth $\tau_{\rm s}\simeq\kappa_{\rm s}\rho(r_{\rm tr})r_{\rm tr}\simeq c/v$, i.e.

[TABLE]

where $\kappa_{\rm s}=0.34\rm\,cm^{2}\,g^{-1}$ is the Thomson scattering opacity for solar metallicity. Here we have assumed that the Rosseland-mean opacity roughly equals to the scattering opacity. The scattering photospheric radius $r_{\rm scat}=5.1\times 10^{15}\rm\,cm$ is typically not important in determining the optical appearance of a TDE.

In the radius range $r_{\rm in}<r<r_{\rm tr}$, photons are advected by the expanding wind and the radiation energy density evolves as $U(r)\propto\rho^{4/3}\propto r^{-8/3}$ (Strubbe & Quataert, 2009). Above the radius $r_{\rm tr}$, photons rapidly diffuse away from the local fluid. Since the diffusive flux is given by $F_{\rm diff}\simeq Uc/\tau_{\rm s}\propto r^{-2}$, we have $U(r)\propto r^{-3}$ for $r>r_{\rm tr}$. The normalization for the above scalings for radiation energy density is given by $L_{\rm in}=4\uppi r_{\rm in}^{2}U(r_{\rm in})v$, which means

[TABLE]

We assume that the radiation is well thermalized near $r_{\rm in}$, so the radiation spectrum is nearly a blackbody up to $r_{\rm tr}$ and the radiation temperature profile is

[TABLE]

At larger radii $r>r_{\rm tr}$, the temperature profile depends on whether the majority of photons get thermalized due to a combination of scattering and absorption. In the following, we describe a semi-analytical way of capturing the effect of frequency-dependent thermalization.

At each radius, we define a blackbody temperature $T_{\rm BB}\equiv(U/a)^{1/4}$, which is the temperature the radiation field would have if LTE is achieved. Since the emissivity and absorption opacity are strongly frequency-dependent (due to bound-free edges and lines), it is difficult to achieve an equilibrium between emission and absorption at all frequencies. Instead, we define a rough LTE criterion (see Nakar & Sari, 2010) which is applicable at $r>r_{\rm tr}$,

[TABLE]

where $B_{\nu}(T_{\rm BB})$ is the Planck function at temperature $T_{\rm BB}$, $j_{\nu}(T_{\rm BB})$ is obtained from Cloudy999Version 17.01 of the code last described Ferland et al. (2017). by assuming the gas is under a thermal radiation bath of temperature $T_{\rm BB}$, and the diffusion time is given by $t_{\rm diff}=\tau_{\rm s}r/c$. Then we dopt a critical value $\eta_{\rm crit}=5/4$ such that the radiation is considered to be in LTE at radii where $\eta(r)<\eta_{\rm crit}$ and non-LTE otherwise. This critical value means that equilibrium between emission and absorption is achieved at about 80% of the frequencies near the peak of the overall spectrum. Thus, the frequency-averaged thermalization radius $r_{\rm th}$ is given by

[TABLE]

Then, the two characteristic radii $r_{\rm tr}$ and $r_{\rm th}$ determine the radial profile of the radiation temperature $T(r)$, which has three power-law segments: $T\propto r^{-2/3}$ ($r_{\rm in}<r<r_{\rm tr}$), $T\propto r^{-3/4}$ ($r_{\rm tr}<r<r_{\rm th}$, assuming $r_{\rm th}<r_{\rm scat}$), and $T=\rm const$ ($r>\mathrm{max}(r_{\rm tr},r_{\rm th})$). Note that in the case where $r_{\rm tr}>r_{\rm th}$, the middle segment does not exist. The mean photon energy the observer sees is given by $2.7k_{\rm B}T[\mathrm{max}(r_{\rm tr},r_{\rm th})]$. With the radiation temperature $T(r)$, energy density $U(r)$, and density $\rho(r)$ at each radius (for a logarithmic radial grid), we use Cloudy to compute the degree of ionization for each chemical species and their energy-level population, under Solar abundance.

We make use of the volumetric emissivity $j_{\nu}(r)$ (for a logarithmic frequency grid) output from Cloudy. At radius $r>r_{\rm tr}$, the energy of photon are still significantly modified by electron scattering. This is because the local intensity distribution is anisotropic with an outwards diffusive flux. This intensity anisotropy means that, at a given radius, an electron scatters more outward-going photons than inward-going ones, and hence photons overall exert a force on this electron. Since the electron is moving outwards at velocity $v$, this force due to photon scattering is doing work to accelerate the electron (of course, this electron is dynamically coupled with a proton such that the actual acceleration is small). The net effect of the photon-electron momentum transfer is that, photons lose a fraction \sim$$v/(c\tau) of energy over each scattering (see Appendix C). Since it takes $\tau^{2}$ scatterings for each photon to escape, the total amount of energy loss is \sim$$\tau v/c. We are interested in the region at $r>r_{\rm tr}$ or $\tau<c/v$, so photons lose energy by less than a factor of 2 (and hence overall adiabatic cooling is not important) but this energy shift is important for the transport of line photons by effectively broadening the lines (Pinto & Eastman, 2000; Roth et al., 2016). We take a broadening factor of

[TABLE]

and perform a Gaussian kernel smoothing over the Cloudy output of $j_{\nu}(r)$ at each radius.

Now we have all the ingredients to calculate the specific luminosity of the escaping photons from the wind. For each frequency $\nu$, the thermalization radius $r_{\rm th,\nu}$ is given by the equilibrium between emission and absorption, i.e.

[TABLE]

which is equivalent to the effective absorption optical depth $\tau_{*,\nu}\simeq\sqrt{\tau_{\rm a,\nu}(\tau_{\rm a,\nu}+\tau_{\rm s})}\simeq 1$ (Rybicki & Lightman, 1979). Then, the specific luminosity is roughly given by

[TABLE]

where $\tau_{\rm s}(r_{\rm th,\nu})=\kappa_{\rm s}\rho(r_{\rm th,\nu})r_{\rm th,\nu}$. As shown in Fig. 13, our model can reproduce the optical and UV spectral-energy distributions (SEDs) of typical TDE candidates such as ASASSN-14li and PS1-10jh. One robust prediction of our wind reprocessing model is that the SED in the NIR band is softer than that in the optical-UV, typically $L_{\nu}\propto\nu^{\sim 0.5}$. This can be explained as follows.

The absorption opacity in the NIR continuum is dominated by free-free transitions (ignoring the Gaunt factor, Rybicki & Lightman, 1979)

[TABLE]

where the density $\rho$ and temperature $T$ are in units of $\rm g/cm^{3}$ and $\rm K$, respectively. In the limit $\kappa_{\rm low\,\nu}\ll\kappa_{\rm s}$, the effective opacity is given by $\kappa_{\rm*,\nu}\simeq\sqrt{\kappa_{\rm low\,\nu}\kappa_{\rm s}}$ (Rybicki & Lightman, 1979), so the frequency-dependent thermalization radius is given by $\kappa_{\rm*,\nu}\rho r_{\rm th,\nu}\simeq 1$, i.e.,

[TABLE]

where $T=10^{4.5}T_{4.5}\rm\,K$ is the electron temperature at the thermalization radius (the final results depends very weakly on $T$). The above equation agrees reasonably well with Fig. 13.

At frequencies with $r_{\rm th,\nu}<r_{\rm tr}$, thermalization occurs below the trapping radius, and the escaping specific luminosity is given by $L_{\nu}=4\uppi r_{\rm tr}^{2}(4\uppi)B_{\nu}(T(r_{\rm tr}))/\tau_{\rm s}(r_{\rm tr})$, which has a blackbody shape at temperature $T(r_{\rm tr})$. However, at frequencies with $r_{\rm th,\nu}>r_{\rm tr}$, thermalization occurs above the trapping radius, and eq. (29) gives

[TABLE]

which applies for at low frequencies (such that $r_{\rm th,\nu}>r_{\rm tr}$)

[TABLE]

This behavior $L_{\nu}\propto\nu^{0.5}$ should be observable in the NIR (see Figs. 4 and 5 of Roth et al., 2016). This effect is analogous to the radio/infrared free-free absorption in the wind of Wolf-Rayet stars (Wright & Barlow, 1975; Crowther, 2007). The weak dependence on the electron temperature $T^{-1/8}$ means that eq. (32) can be used to measure the “wind parameter” $K\propto\dot{M}/v$ for individual TDEs, similar to measuring the mass-loss rate from Wolf-Rayet stars101010The CIO is likely clumpy (due to e.g. episodic mass ejection), so a further correction for the volume filling factor $f_{\rm V}<1$ is needed (Osterbrock & Flather, 1959)..

5.2.3 Other pieces of information — lines and X-rays

The observed H$\alpha$ and HeII emission lines have complex and sometimes double-peaked or boxy structures (e.g. Arcavi et al., 2014; Holoien et al., 2016, 2018a; Blagorodnova et al., 2018). They have been modeled with the reprocessed emission from an elliptical disk (Liu et al., 2017; Holoien et al., 2018a). However, these elliptical disks may be highly unstable on timescales $\sim$months because each annulus undergoes apsidal precession at a different rate. In our picture, the emission line profiles are mainly controlled by the bulk motion of the line formation region of the CIO (at a few times the trapping radius $r_{\rm tr}$), which can either be blue- or red-shifted depending on the observer’s line of sight. We also note that, if the line formation region has large scattering optical depth, the line profile may be further modified by Comptonization (Roth & Kasen, 2018). The N III and O III emission lines in some TDEs, e.g. AT2018dyb (Leloudas et al., 2019), are probably due to Bowen fluorescence, which requires a large flux of (unseen) EUV photons.

The partial sky coverage of the CIO provides a unification of the diverse X-ray properties of optically selected TDEs. When the line of sight to the inner accretion disk is not blocked by the CIO, the observer should see optical emission as well as the EUV or soft X-ray emission from the inner accretion disk or its wind (Strubbe & Quataert, 2009; Dai et al., 2018; Curd & Narayan, 2019). When the line of sight is only blocked by the region of the CIO with modest optical depth, the observer may see blueshifted absorption lines from high ionization species (e.g. Brown et al., 2018; Blagorodnova et al., 2018). When the line of sight is blocked by the highly optically thick region of CIO, the observer only sees optical emission initially. Then, as the CIO’s mass outflowing rate drops with time, the trapping radius shrinks and hence the EUV and soft X-ray photons from the inner disk suffer less adiabatic loss. As a result, the soft X-ray flux (on the Wien tail) should gradually rise and the spectrum hardens on timescales of $\sim 1\,$yr (Margutti et al., 2017; Gezari et al., 2017; Holoien et al., 2018b).

5.2.4 Radio emission from non-jetted TDEs

In this subsection, we discuss the radio emission from the adiabatic shock driven by the CIO into the circum-nuclear medium (CNM). As shown in Fig. 10, the CIO has kinetic energies from E_{\rm k}\sim$$10^{50}\rm\,erg up to \sim$$10^{52}\rm\,erg and mean speed between v_{0}\sim$$0.01c and \sim$$0.1c. In the following, we simplify the complex CIO structure as a thin shell covering a solid angle $\Omega$ within which the density and velocity distributions are uniform. We assume that the ambient medium has a power-law density profile in the radial direction $n=n_{\rm pc}r_{\mathrm{pc}}^{-k}$ ($k<3$), where $r_{\rm pc}=r/\mathrm{pc}$. We also ignore sideway expansion of the shocked region since $\Omega\sim 2\uppi$, so the system is one dimensional.

When the CIO reaches a radius $r$, the total number of shocked electrons from the CNM is given by

[TABLE]

where $N_{\rm pc}\equiv n_{\rm pc}\times(1\mathrm{\,pc})^{3}$ is a reference number of electrons. We ignore the acceleration of particles by the reverse shock (driven into the ejecta) because it is much weaker than the forward shock (driven into the CNM). The deceleration radius $r_{\rm dec}$ is given by $E_{\rm k}=(1/2)N(r_{\rm dec})m_{\rm p}v_{0}^{2}$ ($m_{\rm p}$ being proton mass), which means

[TABLE]

We smoothly connect the free-expansion phase with the Sedov-Taylor phase by using the following velocity profile

[TABLE]

and hence the shock reaches radius $r$ at time

[TABLE]

The electron number density in the shocked region is $4n(r)$ and the mean energy per proton is $(1/2)m_{\rm p}v(r)^{2}$, so the thermal energy density is $2n(r)m_{\rm p}v(r)^{2}$. We assume that a fraction $\epsilon_{\rm B}\ll 1$ of the thermal energy is shared by magnetic fields, so the magnetic field strength is

[TABLE]

We assume that electrons share a fraction $\epsilon_{\rm e}\ll 1$ of the thermal energy and that they are accelerated to a power-law momentum distribution with index $p$. We expect particle acceleration from non-relativistic shocks to give $2<p<3$, both theoretically (Bell, 1978; Blandford & Eichler, 1987; Malkov & Drury, 2001; Park et al., 2015; Caprioli & Spitkovsky, 2014) and observationally (Chevalier, 1998; Green, 2014; Zanardo et al., 2014). For fast shocks where the mean energy per electron $\epsilon_{\rm e}m_{\rm p}v^{2}/2\gg m_{\rm e}c^{2}$ ($m_{\rm e}$ being electron mass), the majority of the particle number and kinetic energy are both concentrated near a relativistic minimum momentum $\gg m_{\rm e}c$. For slow shocks where $\epsilon_{\rm e}m_{\rm p}v^{2}/2\lesssim m_{\rm e}c^{2}$, most particles have non-relativistic momenta but the majority of kinetic energy is in mildly relativistic particles with Lorentz factor $\gamma\sim 2$. We are interested in the number density of ultra-relativistic electrons. These two regimes above can be smoothly connected by assuming a power-law Lorentz factor distribution $\mathrm{d}N/\mathrm{d}\gamma\propto\gamma^{-p}$ above the minimum Lorentz factor (Granot et al., 2006; Sironi & Giannios, 2013)

[TABLE]

Then the normalization is given by the total energy of these relativistic electrons being $\epsilon_{\rm e}N(r)m_{\rm p}v(r)^{2}/2$, i.e.,

[TABLE]

An electron of Lorentz factor $\gamma\gg 1$ has characteristic synchrotron frequency

[TABLE]

where $e$ is the electron charge. Since the peak specific power is $P_{\nu,\rm max}\simeq e^{3}B/m_{\rm e}c^{2}$, the specific luminosity at frequency $\nu$ in the optically thin regime is given by

[TABLE]

The synchrotron self-absorption frequency $\nu_{\rm a}$ and the corresponding Lorentz factor $\gamma_{\rm a}$ are defined where the optical depth $\alpha_{\nu_{\rm a}}\Delta\ell_{\rm r}\sim 1$ ($\Delta\ell_{\rm r}$ being the radial thickness of the emitting region). The volumetric emissivity at $\nu_{\rm a}$ is given by $j_{\nu_{\rm a}}=\alpha_{\nu_{\rm a}}2kT\nu_{\rm a}^{2}/c^{2}$ (in the Rayleigh-Jeans limit $h\nu_{\rm a}\ll kT$), where $kT\simeq\gamma_{\rm a}m_{\rm e}c^{2}$ is the temperature of electrons responsible for absorption. Assuming $\gamma_{\rm a}>\gamma_{\rm m}$ (which will later be shown to be true for non-relativistic shocks), we can write the specific luminosity as $4\uppi j_{\nu_{\rm a}}\Omega r^{2}\Delta\ell_{\rm r}$, and hence

[TABLE]

Combining eqs. (42) and (43), we see that the Lorentz factor $\gamma_{\rm a}$ corresponding to the self-absorption frequency is determined by

[TABLE]

which gives

[TABLE]

where we have defined a reference column density $\mathcal{N}_{\rm pc}=n_{\rm pc}\times 1\,$pc. If $\nu(\gamma_{\rm m})\ll\nu_{\rm a}$ and synchrotron/inverse-Compton cooling are negligible, the synchrotron spectrum when the shock is at radius $r$ is given by (Granot & Sari, 2002)

[TABLE]

In Fig. 14, we show the radio emission from CIO colliding with the CNM for a number of cases. We denote the average velocity $v_{0}$ in units of $c$ and the kinetic energy (with unit $\rm erg$) in log-scale. The three cases with $(v_{0},E_{\rm k})=(0.01,50),\ (0.03,51),\ (0.1,52)$ are motivated by the mean velocities and kinetic energies in Fig. 10. The case with $(v_{0},E_{\rm k})=(0.2,52)$ is for comparison with that with $(v_{0},E_{\rm k})=(0.1,52)$, and we see that an outflow with higher velocity generates brighter and earlier-peaked radio emission. The cases with $(v_{0},E_{\rm k})=(0.1,49)$ and (0.1, 50) are motivated by the fact that the CIO velocity profile is non-uniform with a fraction of the mass moving faster than the mean velocity. We find that the faster portion of the ejecta generates bright radio emission at early time. For each combination of $(v_{0},E_{\rm k})$, we take two different CNM density normalizations $n_{\rm pc}=100$ and $3\rm\,cm^{-3}$. As expected, we find that, for higher CNM densities, the radio emission is brighter and peaks earlier.

We also show the data from several TDEs for comparison but do not intend to search for the best-fit parameters for individual cases. The upper limits for iPTF16fnl (Blagorodnova et al., 2017) reported at $15\,$GHz have been scaled by a factor of $(5/15)^{(1-p)/2}$ (assuming $\nu_{\rm a}>5\,$GHz). The upper limits for the X-ray selected TDEs reported at $3\,$GHz by Bower et al. (2013) are not scaled.

Even though we keep the following parameters fixed $k$=1.5, $p$=2.4, $\epsilon_{\rm e}$=0.1, $\epsilon_{\rm B}$=0.01, and $\Omega$=2$\uppi$, the radio luminosity and duration are extremely diverse. Generally, we expect TDEs with CIO to have some radio emission at the level of ASASSN-14li lasting for years up to centuries. We also note that radio emission from the jetted TDE Swift J1644+57 (Bloom et al., 2011; Burrows et al., 2011; Zauderer et al., 2011) is much brighter (and peaks earlier) than that from the CIO, because this source was powered by a relativistic jet pointing towards the observer. For off-axis jetted TDEs, the radio emission due to the CIO may be mistaken as a signature of jets (a possible way of distinguishing between them is to resolve the motion of the radio emitting region by long-baseline interferometry).

Another possible source of wide-angle outflow is the wind expected from super-Eddington accretion in TDEs with BH masses $M\lesssim 10^{7}\mathrm{M_{\odot}}$ (Strubbe & Quataert, 2009; Sa̧dowski et al., 2014; Jiang et al., 2017). In fact, the super-Eddington wind may be more powerful than the CIO, because the energy efficiency of the CIO is only \sim$$r_{\rm g}/r_{\rm I}\sim 0.001 to 0.01. Thus, we expect TDEs with strong super-Eddington wind to generate radio emission comparable to or even brighter than that in our $(v_{0},E_{\rm k})=(0.2,52)$ case. Late-time radio observations can potentially test whether super-Eddington accretion flows generate jets or winds. We note that the unbound tidal debris typicall has very small solid angle (Guillochon et al., 2014), so its radio emission (and reprocessing of the high-energy photons from the disk) is much weaker than that of the CIO. It is less likely that the radio emission from ASASSN-14li is caused by the unbound tidal debris (Krolik et al., 2016), unless the star was in a very deeply penetrating $\beta\gg 1$ orbit (Yalinewich et al., 2019).

Finally, we note that the CIO may interact with a pre-existing accretion disk (or the broad line region), if the BH was active before the TDE. If the accretion disk gas is sufficiently dense, the shocks become radiative and bright optical emission like in PS16dtm (Blanchard et al., 2017) may be generated.

5.3 TDE demographics

TDE demographics, in terms of the total TDE rate as a function of BH mass and properties of the disrupted star, has been considered by Stone & Metzger (2016) and Kochanek (2016). In this section, we focus on the rate of optically bright TDEs only, based on the picture that the CIO reprocesses the disk emission from the EUV into the optical band.

We differentiate the TDE rate with three parameters, stellar mass $m_{*}=M_{*}/\mathrm{M_{\odot}}$, impact parameter $\beta$, and BH mass $M$, in the following way

[TABLE]

where $\mathcal{R}$ is the normalization rate per BH in units of $\mathrm{yr}^{-1}$, the normalized stellar mass function satisfies $\int(\mathrm{d}n_{*}/\mathrm{d}m_{*})\mathrm{d}m_{*}=1$, the probability distribution of the impact parameter has also been normalized $\int(\mathrm{d}P/\mathrm{d}\beta)\mathrm{d}\beta=1$, and the BH mass function (BHMF) $\mathrm{d}n_{\rm BH}/\mathrm{d}\mathrm{log}M$ has unit $\rm[Mpc^{-3}\,dex^{-1}]$. The factor $m_{*}^{-1/12}r_{*}^{1/4}$ is because stars with a larger tidal radius are slightly preferred roughly by a factor of $r_{\rm T}^{1/4}$ (MacLeod et al., 2012).

The power-law dependence on the BH mass $M^{\alpha}$ depends, in a non-trivial way, on the stellar density and velocity profiles near the center of individual galaxies. The index $\alpha$ is empirically derived by combining the surface brightness profiles of a sample of galaxies with BH masses inferred from galaxy scaling relations (e.g. Magorrian & Tremaine, 1999; Wang & Merritt, 2004; Stone & Metzger, 2016). There is a core-cusp bimodal distribution of central surface brightness profiles of early-type galaxies used for TDE rate calculations (Lauer et al., 2007). The most recent work by Stone & Metzger (2016) gives $\alpha\simeq-0.25$ for samples of only111111The TDE rates for cusp galaxies are typically $\sim$10 times higher than that for core galaxies of the same BH mass. cusp or core galaxies. For comparison, we also show the results for $\alpha\simeq-0.5$ which do not affect our conclusions qualitatively. We caution that the above studies typically assume a spherically symmetric and time-independent galactic potential, nearly isotropic stellar velocity distribution (except for the loss cone), and the refilling of the loss cone by two-body relaxation only. Other factors, such as massive perturbers, aspherical potential, binary BHs, resonant relaxation, may strongly affect the estimated TDE rate (e.g. Vasiliev & Merritt, 2013; Merritt, 2013). Therefore, we leave the normalization factor $\mathcal{R}$ as a free parameter, which roughly means the (per-BH) rate of TDEs for M-dwarf stars disrupted by \sim$$10^{6}\mathrm{M_{\odot}} BHs.

In loss-cone dynamics, the probability distribution for the impact parameter $\mathrm{d}P/\mathrm{d}\beta$ has two regimes. In the “pinhole” regime (far from the BH), the change in stars’ angular momentum per orbit $\Delta\ell$ is much larger than the size of the loss-cone $\ell_{\rm lc}\simeq\sqrt{2r_{\rm g}r_{\rm T}}$, so $\mathrm{d}P/\mathrm{d}\beta$ simply depends on the “area” of the loss cone per unit change in $\beta$, i.e. $\mathrm{d}P/\mathrm{d}\beta\propto\beta^{-2}$. In the “diffusive” regime (near the BH), $\Delta\ell\ll\ell_{\rm lc}$ and hence stars are always disrupted near the boundary of the loss-cone with minimum penetration depth, i.e. $\mathrm{d}P/\mathrm{d}\beta$ is nearly a $\delta$-function. The fraction of TDEs in the pinhole regime $f_{\rm pin}$ depends on the detailed stellar density profile near the BH and has large uncertainty at each BH mass. Following Kochanek (2016), we take

[TABLE]

which is very similar to the fitting result by Stone & Metzger (2016) in the range of BH masses of interest. Then the probability distribution of $\beta$ is given by

[TABLE]

According to eq. (16), the minimum impact parameter is $\beta_{\rm min}\simeq 0.6\,\xi_{*}^{-1}$, which includes relativistic tidal forces for the Schwarzschild spacetime (Kesden, 2012) and $\xi_{*}$ depends on the star’s internal structure. We note that $\xi_{*}$ is not well measured in general relativity even for polytropic stars. Hydrodynamic simulations of disruptions with polytropic or realistic stellar structures in the Newtonian limit ($r_{\rm g}\ll r_{\rm p}$, Guillochon & Ramirez-Ruiz, 2013; Mainetti et al., 2017; Goicovic et al., 2019) show that the star loses about half of the mass when $\xi_{*}\simeq 0.5$ (for polytropic index 4/3) or $\xi_{*}\simeq 1.0$ (for polytropic index 5/3). The former is appropriate for radiative stars with $m_{*}>1.2$ and the latter is good for convective stars with $m_{*}<0.3$ (see a similar treatment by Phinney, 1989). For stars in between $0.3<m_{*}<1.2$, we take a linear interpolation in log$\,m_{*}$ space. Thus,

[TABLE]

We also note that the maximum impact parameter is taken to be infinity, because a star’s orbit can have arbitrarily low angular momentum. The effect of stars being swallowed by the event horizon will be taken into account later when integrating over the BHMF.

We take the Kroupa initial mass function (IMF, Kroupa, 2001) truncated at $m_{\rm*,max}$ (related to the age of the stellar population)

[TABLE]

The two constants $a_{1}$ and $a_{2}$ are given by the continuity at $m_{*}=0.5$ and normalization $\int(\mathrm{d}n_{*}/\mathrm{d}m_{*})\mathrm{d}m_{*}=1$. We ignore compact stellar remnants since they are fewer in number and are typically swallowed as a whole for $M\gtrsim 2\times 10^{5}\mathrm{M_{\odot}}$. We also ignore red giants, because they have long fallback time $P_{\rm min}\sim 11\mathrm{\,yr}\,M_{7}^{1/2}m_{*}^{-1}(r_{*}/10)^{3/2}$ (and even longer circularization time) and do not have an optically thick layer of gas to reprocess the hard disk emission into the optical band (see §5.2.2). The rate of TDEs contributed by binary stars is lower than that from single stars by a factor of \sim$$f_{\rm bi}f_{\rm pin}(R_{*}/a)^{3/4}\ll 1, where $f_{\rm bi}$ is the binary fraction near the galactic center and $a$ is the semimajor axis of the binary orbit. Tidal breakup of the binary has a larger Roche radius $r_{\rm T,b}\simeq(a/R_{*})r_{\rm T}$ and hence occurs at a higher rate than that for single stars by a factor of \sim$$(a/R_{*})^{1/4} (MacLeod et al., 2012). However, stellar disruption is only possible at high impact parameter $\beta\gtrsim a/R_{*}$ in pinhole regime, which means the disruption rate is a factor of $f_{\rm pin}R_{*}/a$ smaller than the tidal breakup rate.

The Kroupa IMF extends down to $m_{*}=0.08$ and then becomes shallower $\mathrm{d}n_{*}/\mathrm{d}m_{*}\propto m_{*}^{0.3}$ for lower mass brown dwarfs. However, TDEs of such low-mass objects likely do not generate much optical emission, the reason being as follows. The mass of the reprocessing CIO can be estimated by $M_{\rm CIO}\sim\rho A\Delta t$ (Lu & Kumar, 2018), where $\rho\simeq(\kappa r_{\rm ph})^{-1}$ is the gas density, $A\simeq L_{\rm opt}/\sigma_{\rm SB}T^{4}\simeq\Omega r_{\rm ph}^{2}$ is the surface area of the optical photosphere, $\kappa$ is the effective absorption opacity, $r_{\rm ph}$ is the photospheric radius, $L_{\rm opt}$ and $T$ are the optical luminosity and blackbody temperature, and $\Delta t$ is the peak duration. Since half of the star’s mass is in unbound tidal debris and only half of the bound mass may be ejected as CIO, we obtain a lower limit for the star’s mass $M_{*}\gtrsim 4M_{\rm CIO}$. Putting in conservative numbers, we obtain

[TABLE]

Fast transients with $\Delta t\lesssim 10\,$d and $L_{\rm opt}\lesssim 10^{43}\rm\,erg\,s^{-1}$ are increasingly likely to have been missed by current surveys. In the following, we take the conservative minimum stellar mass of $m_{\rm*,min}=0.1$. Larger $m_{\rm*,min}$ will lead to lower rates of optically bright TDEs.

We plug eqs. (49), (51) and a given BH mass function into eq. (47) and calculate the integrated volumetric TDE rate

[TABLE]

where the minimum BH mass for CIO launching $M_{\rm cr}$ is given by eq. (13) and the BH mass above which the entire star gets swallowed is given by eq. (15).

The BHMF for $M\lesssim 10^{6.5}\mathrm{M_{\odot}}$ is highly uncertain even in the local Universe. Evolutionary models are constructed by inferring BH growth by the “observed” bolometric luminosity function of active galactic nuclei (AGN). Various treatments of bolometric corrections, radiative efficiency of the accretion disks, and AGN duty cycles may give different results. In this paper, we take two different BHMFs for the local Universe by Merloni & Heinz (2008, MH08) and Shankar et al. (2009, SWM09), as shown in Fig. 19 in the Appendix. The main difference between the two lies in the low-mass end: the MH08 mass function is nearly flat while the SWM09 mass function rapidly diverges121212We see that TDE demographics provide a valuable, direct probe of the BH mass function on the low-mass end. as $\mathrm{d}n_{\rm BH}/\mathrm{d}\mathrm{log}M\propto M^{-0.6}$. Fig. 15 and Fig. 16 shows the TDE demographics for these two BHMFs, respectively. The BHMF can also be directly calculated by applying correlations between BH mass, bulge luminosity and stellar velocity distribution for galaxies in the local Universe, as done131313The two methods of obtaining the BHMF are not independent. Typically, the radiative efficiency of AGN is calibrated by the total BH mass density in the local Universe inferred from galaxy scaling relations (Soltan, 1982; Marconi et al., 2004). by Marconi et al. (2004). We also tried using their BHMF and found that it gives similar results as the MH08 mass function, as shown in Fig. 21 in the Appendix.

On the low BH-mass side, the predicted rate of optically bright TDEs is nearly flat with respect to the BH mass. This is because those TDEs with $M<M_{\rm cr}$ have been filtered out due to insufficient amount of CIO being launched. Our results roughly agree with the rate given by van Velzen (2018), which was based on the “V/Vmax” method and the BH masses are inferred from galaxy scaling relations with updated stellar velocity dispersion by Wevers et al. (2017). We also show the total TDE rate without requiring $M>M_{\rm cr}$ (red curves), which rises more rapidly towards the low-mass end. This is because TDEs favor smaller BHs by the factor $M^{\alpha}$ (we have taken $\alpha=-0.25$ or $-0.5$) and the BHMF model of Shankar et al. (2009) diverges towards the low-mass end (the MH08 model has a shallower behavior). Unfortunately, the current small-number statistics are not able to discriminate between the two scenarios (shown in blue and red curves) at a significant confidence level.

Thus, our picture predicts that the majority of TDEs by BHs with $M\lesssim 10^{6}\mathrm{M_{\odot}}$ are not optically bright and will hence be missed by current optical transient surveys. The rate of optically bright TDEs is a factor of $\sim$10 or more141414 In Figs. 15 and 16, if we take away the $m_{*}\in(0.1,0.2)$ bin due to insufficient mass for the reprocessing layer, the rate of optically bright TDEs will be lower by a factor of $\sim$2 (but the overall shape of the rate as a function of BH mass stays nearly the same). We also note that there could be a large population of TDEs hidden from optical view by dust extinction (Wang et al., 2018). lower than the total TDE rate. Some of these missing TDEs should be observable by wide field-of-view soft X-ray surveys like eROSITA (Cappelluti et al., 2011) and Einstein Probe (Yuan et al., 2015). Our picture can be tested by comparing the detection rates of TDEs in the X-ray and optical bands, although one should keep in mind that low BH-mass TDEs may have a long circularization timescale due to weak apsidal precession (see §5.1).

On the high BH-mass end, the optically bright TDE rate is strongly suppressed due to stars being swallowed by the event horizon, which has been used as a supportive evidence for the existence of BH event horizon (Lu et al., 2017) and that the observed candidates are actually TDEs (van Velzen, 2018). We note that the grey data point near $\mathrm{log}M\sim 8.3$ only contains the TDE candidate ASASSN-15lh, whose nature is still being debated (Dong et al., 2016; Leloudas et al., 2016; Krühler et al., 2018). In our picture, it can be explained by disruption of a relatively massive star $m_{*}\sim 5$ by a non-spinning BH. Disruption of a Sun-like star by a rapidly spinning BH is also possible, because $\beta_{\rm min}$ can be smaller than $0.6\,\xi_{*}^{-1}$ for a prograde orbit (Kesden, 2012; Leloudas et al., 2016).

Finally, we note that rare post-starburst galaxies are over-represented in the current sample of TDE host galaxies by a factor of $\sim$20 to $\sim$100 (Arcavi et al., 2014; French et al., 2016; Law-Smith et al., 2017; Graur et al., 2018), which may be due to higher stellar density concentration near the galactic centers (e.g. Law-Smith et al., 2017; Stone et al., 2018b). Our method also applies to the group of post-starburst galaxies (with a higher rate normalization constant $\mathcal{R}$), as long as their BHMF is similar to that of the entire galaxy population. An important difference is that the age of the stellar population near the centers of post-starburst galaxies may be significantly younger than that for the other normal galaxies, ranging from 100 Myr to 1 Gyr. This will affect the TDE rate on the high BH-mass end. Another potential difference is that the pinhole fraction $f_{\rm pin}$ may be lower for more cuspy (steeper) stellar density distribution in post-starburst galaxies (Stone et al., 2018b).

6 Discussion

In this section, we discuss a number of issues that require further thoughts in future works.

(1) The stream self-intersection may be delayed due to Lense-Thirring (LT) precession, if the BH’s spin is misaligned with the angular momentum of the stellar orbit (e.g. Kochanek, 1994; Dai et al., 2013; Guillochon & Ramirez-Ruiz, 2015; Hayasaki et al., 2016). For highly eccentric orbits, the angle by which the orbital angular momentum vector precesses over one period is $(\Delta\omega)_{\rm LT}\approx 4\uppi a(r_{\rm g}/2r_{\rm p})^{3/2}\sin i$ (to leading post-Newtonian order), where $a$ is the dimensionless spin of the BH and $i$ is the inclination angle ($i=0$ for spin-orbit alignment). For a given orbit, we express the maximum ratio of the stream width over the distance to the BH as $(H/r)_{\rm max}=f_{\rm H}\beta R_{*}/r_{\rm T}$, where $f_{\rm H}$ describes possible broadening of the stream due to apsidal/LT precession151515Without LT precession, the ratio between the velocity perpendicular to the orbital plane $v_{\perp}$ and the velocity within the orbital plane $v_{\parallel}$ is of order $v_{\perp}/v_{\parallel}\sim(M_{*}/M)^{1/3}\ll 1$. However, for strong LT precession $(\Delta\omega)_{\rm LT}\gtrsim v_{\perp}/v_{\parallel}$, a fraction of the $v_{\parallel}$ component is aligned with the direction of vertical compression, so the stream width after the bounce may be broader than in the case without LT precession. Strong apsidal precession can also cause the tidal compression in the orbital plane to be oblique and hence part of the orbital velocity may be dissipated near the pericenter., hydrogen recombination, and magnetic fields. Then, intersection may be avoided for a particular orbit when $(\Delta\omega)_{\rm LT}\gtrsim(H/r)_{\rm max}$, i.e.

[TABLE]

We can see that TDEs by slowly spinning $a\ll 1$ low-mass $M_{6}\lesssim 1$ BHs are expected to have prompt intersection between the first and second orbits (as shown in Fig. 1). For rapidly spinning high-mass BHs, intersection may be avoided promptly (if $f_{\rm H}\sim 1$) but will eventually occur with a delay. From the point of view of an observer who defines $t=0$ as the moment of stream intersection, the delay itself is not important, since the mass flux of the stream stays unchanged. On the other hand, as long as the intersection occurs between two adjacent orbits (the n-th and the n+1-th, as found by Guillochon & Ramirez-Ruiz, 2015), the intersection radius ($r_{\rm I,LT}$) under LT precession is roughly the same as that without LT precession ($r_{\rm I}$). This is because for most TDEs the apsidal precession angle ($3\uppi r_{\rm g}/r_{\rm p}$) is much larger than the LT precession angle. Thus, the intersection radius $r_{\rm I}$, angle $\widetilde{\theta}_{\rm I}$ and velocity $\widetilde{v}_{\rm I}$ calculated in the Schwarzschild spacetime are similar to those for spinning BHs. Therefore, our model for the hydrodynamical collision process, including redistribution of angular momentum/specific energy and the possibility of launching the CIO, should be largely applicable.

(2) TDE demographics on the high BH-mass end depends on the spin distribution. In the case where the star’s initial angular momentum is parallel to the BH’s spin angular momentum, the pericenter radius of the marginally bound parabolic orbit is $r_{\rm mb}=(1+\sqrt{1-a})^{2}r_{\rm g}$ (Bardeen et al., 1972), where $-1<a<1$ is the spin parameter of the BH ($a<0$ for retrograde orbits). The marginal disruption case corresponds to $r_{\rm p}\approx r_{\rm mb}\approx 5^{1/3}\xi_{*}r_{\rm T}$, which gives the maximum mass for Kerr BHs hosting TDEs

[TABLE]

We can see that the $M_{\rm max,Kerr}$ is strongly affected by the BH spin only when $a\gtrsim 0.5$. The maximum BH mass is also affected by the age of the stellar population ($m_{\rm*,max}$, stellar interior structure $\xi_{*}$, and the number of evolved subgiants). It may be difficult to extract the information on the BH spin distribution from TDE rate on the high BH-mass end. We also note that the critical mass $M_{\rm cr}$ (above which significant amount of CIO is launched) is mainly affected by the spin-independent de Sitter term of the apsidal precession, so the TDE demographics on the low BH-mass end should be insensitive to the BH spin distribution.

(3) We have assumed that the two colliding streams have the same cross-section and that there is no offset in the transverse direction. This is reasonable provided that (i) all processes occurring when the fallback stream passes near the pericenter $r\sim r_{\rm p}$ before the collision are largely reversible161616This means that, if we denote the two colliding ends as $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ (in chronological order) and reverse the velocity at $\mathcal{C}_{2}$, the stream will evolve back to the conditions at $\mathcal{C}_{1}$ (except for the velocity being in the opposite direction). and that (ii) the angular stream widths $H/r$ are larger than than the amount of LT precession per orbit. However, there could be many irreversible processes occurring near the pericenter, including: (i) apsidal and LT precession causing the tidal compression to be oblique (instead of perpendicular to the orbital velocity); (ii) mass loss along with the bounce following the tidal compression; (ii) viscosity causing exchange of angular momentum between adjacent shear layers. The general relativistic evolution of the fallback stream over multiple orbits is still an open question, mainly because the extremely large aspect ratio makes it a challenging task for numerical simulations. If these irreversible processes are indeed important and eq. (54) is satisfied, our current model needs two additional parameters: the ratio of the cross-sections of the two colliding streams and the fractional offset in the transverse direction. The hydrodynamics of the stream-stream collision and the subsequent expansion of the shocked gas may be largely modified. This is out of the scope of the current work and should be studied in the future.

(4) The energy radiated in the UV-optical band is typically \lesssim$$10^{51}\,erg, which is much smaller than the energy budget of the system, even assuming radiatively inefficient accretion (Piran et al., 2015; Lu & Kumar, 2018). In our picture, this “missing energy” puzzle may be explained in two possible scenarios. The first is that the disk bolometric emission is capped near the Eddington level for an extended amount of time171717The late-time (5-10 yrs) UV-optical emission from a number of TDEs reported by van Velzen et al. (2018b) supports this scenario, but the it is also possible that the late-time excess is due to dust scattering echo (which has been seen in many supernovae). ($\gg P_{\rm min}$) but the CIO reprocesses the disk emission in to the UV-optical band only for a timescale of order $P_{\rm min}$ (then the reprocessed emission moves into the EUV and soft X-ray as the trapping radius shrinks). Since the observed peak UV-optical luminosity is also near the Eddington limit (e.g. Wevers et al., 2017), the efficiency for reprocessing, defined as the observed UV-optical luminosity divided by the intrinsic disk luminosity, is required to be of order unity in this case. The second scenario is that the disk bolometric emission significantly exceeds the Eddington limit (as in the simulations by Jiang et al., 2017), but the reprocessing efficiency is much less than unity. The reason for a low reprocessing efficiency could be that, if $r_{\rm tr}\gg r_{\rm in}$, photons are trapped in the expanding CIO and hence their energy is adiabatically lost in the form of $PdV$ work. However, detailed radiation-hydrodynamic simulations are needed to distinguish between these two scenarios.

7 Summary

We have described a semi-analytical model for the dynamics of TDEs, including the properties of the fallback stream before the self-intersection and the fate of the shocked gas after the intersection. We circumvent the computational challenge faced by previous TDE simulation works by assuming that the post-disruption bound stream follows the geodesics in the Schwarzschild spacetime until the self-intersection. Then we numerically simulate the (non-relativistic) hydrodynamical collision process in a local box at the intersection point. Since the cross-sections of the two colliding streams are much smaller than the size of the orbit and the streams are pressureless (or cold) before the collision, the collision process and the expanding structure of the shocked gas are self-similar. This allows us to explore a wide range of TDE parameter space in terms of the stellar mass, BH mass, and impact parameter. Our method provides a way for global simulations of the disk formation process by injecting gas at the intersection point according to the velocity and density profiles (eqs. 9 and 10) shown in this paper.

The most important observational implication is that a large fraction of the fallback gas can be launched in the form of a collision-induced outflow (CIO) when the BH mass is above a critical value $M_{\rm cr}$ (eq. 13). We propose that the CIO is responsible for reprocessing the accretion disk emission from the EUV or soft X-ray to the optical band. This picture can naturally explain the large photospheric radius of \sim$$10^{14}–$10^{15}\,$cm (or low blackbody temperature of a few$\times 10^{4}\,$K), and the typical widths of the H and/or He emission lines. We predict the CIO-reprocessed spectrum in the infrared to be $L_{\nu}\propto\nu^{\sim 0.5}$, shallower than a blackbody. A blackbody fit to the optical SED, as commonly done in the literature, may underestimate the true color temperature. Our picture is different from that of Piran et al. (2015) in that the radiation energy ultimately comes from the accretion flow rather than the stream collision (which is shown to be nearly adiabatic, Jiang et al., 2016). Our model is also different from that of Metzger & Stone (2016) in that we identify the physical origin of the “reprocessing layer” and that this layer is aspherical.

The partial sky coverage of the CIO provides a natural unification of the diverse X-ray behaviors of the optically selected TDEs. Depending on the observer’s line of sight, an optically bright TDE may show strong X-ray emission (when the inner disk is not veiled) or weak/no X-ray emission (when the inner disk is veiled), which agrees with the large range of X-ray to optical peak flux ratios: \sim$$10^{-4} for iPTF16fnl (Blagorodnova et al., 2017), \sim$$10^{-2} for AT2018zr (van Velzen et al., 2018a), and $\sim$1 for ASASSN-14li (Holoien et al., 2016). As the CIO’s mass outflowing rate drops with time, the X-ray fluxes for veiled TDEs may gradually rise with time, as observed in ASASSN-15oi and -15lh (Margutti et al., 2017; Gezari et al., 2017; Holoien et al., 2018b). Our picture is different from those of Dai et al. (2018) and Curd & Narayan (2019) which describe that the X-ray to optical flux ratio is controlled by the observer’s viewing angle with respect to rotational axis of the accretion disk (instead of the CIO’s outflowing direction).

In cases where the CIO is launched (BH mass $M>M_{\rm cr}$), the rest of the fallback gas is left in more tightly bound orbits with higher (sometimes negative) specific angular momentum than the original star, and hence the circularization process is expected to occur on timescale of order \sim$$P_{\rm min} after the onset of intersection. If this is confirmed by future simulations, then it explains the rise/fade timescale ($\sim$months) of optically bright TDEs. We note that the circularization radius of the accreting gas may be different from $2r_{\rm p}$ (as generally assumed in the literature, e.g. Rees, 1988; Strubbe & Quataert, 2009; Shen & Matzner, 2014). Another unexpected prediction is that, in some cases, the accretion disk rotates in the opposite direction as that of the initial star.

The total kinetic energy of the CIO spans a wide range from \sim$$10^{50} up to a few$\times 10^{52}\,$erg (in rare cases). The mass-weighted mean speed varies from \sim$$0.01c to \sim$$0.1c. The shocks driven into the ambient medium by this outflow can produce radio emission with highly diverse timescales and peak luminosities, depending on the density profile of the ambient medium, CIO’s velocity and energy, and microphysics of particle acceleration/magnetic field amplification by the shocks. The radio emission from ASASSN-14li and a few other TDE candidates may be from the afterglow of the CIO (instead of the unbound tidal debris, which typically has a much narrower solid angle).

We also find that the volumetric rate of optically bright TDEs is nearly flat with respect to the BH mass in the range $M\lesssim 10^{7}\mathrm{M_{\odot}}$. This is because TDEs with $M<M_{\rm cr}$ have been filtered out due to lack of significant amount of CIO. Our results roughly agree with the BH mass function of optically selected TDEs obtained by van Velzen (2018). This filtering leads to an optical TDE rate that is a factor of $\sim$10 or more lower than the total TDE rate (without requiring $M<M_{\rm cr}$). For TDEs by BHs with $M<M_{\rm cr}$, the stream self-intersection becomes less and less efficient at dissipating the orbital energy and other mechanisms such as MHD turbulence may be responsible for driving the formation of a circular disk. The circularization timescale of these TDEs may be much longer than $P_{\rm min}$. Some of them should be observable by wide field-of-view X-ray surveys like eROSITA (Cappelluti et al., 2011) and Einstein Probe (Yuan et al., 2015). Our model can be tested by comparing the rates of TDEs in the X-ray and optical bands.

8 acknowledgments

We are grateful for the discussions with with Eliot Quataert and Tony Piro on the radiative transfer in a reprocessing wind, with Brian Metzger and Dan Kasen on the absorption opacity, and with Sterl Phinney on the effects of stellar interior structure. We also thank Phil Hopkins, Bing Zhang, Shri Kulkarni, Brad Cenko, Sjoert van Velzen, and Jim Fuller for useful conversations. We are indebted to Brian Metzger and Pawan Kumar for reading an earlier version of the manuscript and providing valuable comments. We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. This research benefited from interactions at the ZTF Theory Network Meeting, funded by the Gordon and Betty Moore Foundation through Grant GBMF5076 and by the National Science Foundation under Grant No. NSF PHY-1748958. WL was supported by the David and Ellen Lee Fellowship at Caltech.

Appendix A Supplemental Figures

In the Appendix, we provide a number of figures to support the main content. Their descriptions are in the captions.

Appendix B Low Radiative Efficiency of Self-intersection Shocks

We justify the usage of adiabatic equation of state for our hydrodynamic simulation in §3. At the collision point, the Thomson scattering optical depth of the stream in the transverse direction is given by $\tau_{\rm st}\simeq 2\kappa_{\rm s}\rho H$, where $H$ is the radius of the stream (assuming cylindrical shape), $\rho$ is the gas density, and $\kappa_{\rm s}=0.34\rm\,cm^{2}\,s^{-1}$ is the scattering opacity. The mass flowing rate of the stream is $\dot{M}_{\rm fb}\simeq\uppi H^{2}\rho v$, and hence

[TABLE]

where we have used conservative values for the flow velocity $v$, transverse radius $H$, and mass flowing rate (see eq. 19). After the collision, the pressure of the shocked gas is dominated by radiation, and adiabatic expansion converts internal energy back into (roughly spherical) bulk motion over a radius of r$$\,\sim\,a few$\times H$. Beyond this radius, the gas expands with nearly constant velocity $v$ and the density drops with radius as $\rho\propto r^{-2}$ (see Fig. 4). Radiation is advected by the expanding shocked gas until the photon trapping radius $r_{\rm tr}\simeq\kappa_{\rm s}\dot{M}/(4\uppi c)$. The radiative efficiency, i.e. the ratio between the emergent luminosity and the total kinetic power, is roughly given by

[TABLE]

We see that the stream self-intersection is radiatively inefficient.

Appendix C Photon Down-scattering in Diffusion Region

Consider photons diffusing through a scattering slab in the $\hat{\boldsymbol{x}}$ direction. For a general form of angular dependence for the intensity

[TABLE]

where $I=\int\mathrm{d}\nu I_{\nu}$, $P_{n}(\mu)$ are the Legendre polynomials (only the first three terms are important here) $P_{0}=1$, $P_{1}=\mu$, $P_{2}=(3\mu^{2}-1)/2$, and the orthonormality gives $\int P_{m}P_{n}\mathrm{d}\mu=2\delta_{mn}/(2n+1)$ ($\delta_{mn}$ being the Kronecker delta). The energy density, flux, and pressure of the radiation field are given by the different moments of intensity $U=(2\uppi/c)\int I(\mu)\mathrm{d}\mu=4\uppi I_{0}/c$, $F=2\uppi\int I(\mu)\mu\mathrm{d}\mu=4\uppi I_{1}/3$, and $P=(2\uppi/c)\int I(\mu)\mu^{2}\mathrm{d}\mu=(4\uppi/3c)(2I_{2}/5+I_{0})$. Now take an electron moving at velocity $\beta c$ and Lorentz factor $\gamma$ in the $\hat{\boldsymbol{x}}$ direction. We Lorentz transform the radiation field from the lab frame to the comoving frame of the electron, where quantities are denoted with a prime (′). In the comoving frame, the electron gains momentum $\Delta p^{\prime}$ over time $\Delta t^{\prime}$ due to scattering

[TABLE]

where $\sigma_{\rm T}$ is the Thomson cross-section. Making use of the orthonormality relations, the integral above can be expressed in terms of $I_{0}$, $I_{1}$, and $I_{2}$ and hence $U$, $F$, and $P$. Going back to the lab frame, the electron gains energy at a rate given by

[TABLE]

which is at the expense of radiation energy. Therefore, the fractional energy loss of a photon under each scattering (over a timescale $\lambda_{\rm mfp}/c$) is given by

[TABLE]

where the mean free path is $\lambda_{\rm mfp}=(n_{\rm e}\sigma_{\rm T})^{-1}$. In the limit of isotropic radiation field $F=0$ and $P/U=1/3$, the result $\delta\nu/\nu=4\gamma^{2}\beta^{2}/3$ agrees with Rybicki & Lightman (1979, eq. 7.16a). In the diffusion region of the (non-relativistic) CIO, $r_{\rm tr}<r<r_{\rm scat}$, the diffusive flux is given by $F\simeq Uc/\tau_{\rm s}\gg\beta$ ($\tau_{\rm s}$ being the scattering optical depth), so the fractional energy shift per scattering is

[TABLE]

Since a typical photon undergoes \sim$$\tau^{2} scatterings before escaping, the cumulative fractional energy shift is $\tau^{2}\delta\nu/\nu\sim\beta\tau$. This justifies eq. (27).