STELLAR AND GAS DYNAMICAL MODEL FOR TIDAL DISRUPTION EVENTS IN A QUIESCENT GALAXY

It is well known from observations that supermassive black holes (SMBHs) reside at the center of galactic nuclei (Kormendy & Richstone 1995; Kormendy & Ho 2001). If a star passes within the radius ${r}_{t}\sim {R}_{\star }{({M}_{\bullet }/{M}_{\star })}^{1/3}$ of the galaxy's central black hole (BH), the BH's tidal gravity exceeds the star's self-gravity and it is tidally disrupted (Hills 1975; Rees 1988). Stars in the galactic center move in the combined potential field of the SMBH and other stars in the galactic center. The star with orbital energy E is tidally captured if the orbital angular momentum is $J\leqslant {J}_{{lc}}\;=\;\sqrt{2{r}_{t}^{2}({\rm{\Phi }}({r}_{t})-E)}$, where J_lc is the angular momentum of the loss cone (Frank & Rees 1976) and ${\rm{\Phi }}(r)$ is the potential by the BH and other stars in the galactic center. The stellar interactions result in the diffusion of the stars into the loss cone. The stellar distribution function (DF) $f(E,J)$ follows the Fokker–Planck (FP) equation (Bahcall & Wolf 1976; Lightman & Shapiro 1977) and the rate of feeding of stars into the loss cone gives the theoretical tidal disruption event (TDE) rate ${\dot{N}}_{t}.$ The rate of TDE per galaxy depends on the stellar distribution in the galactic center, the SMBH mass M_•, and the structure of galactic nuclei that could be symmetric, axis-symmetric, or triaxial nuclei (Merritt 2013b). Cohn & Kulsrud (1978) obtained the numerical solution to the FP equation for spherical nuclei by means of a detailed boundary layer analysis and applied it to globular clusters. Wang & Merritt (2004) solved the steady state FP equation for the 51 galaxies with the Nuker profile by assuming a single mass star distribution and obtained the ${\dot{N}}_{t}\sim {10}^{-4}-{10}^{-5}\;{\mathrm{Yr}}^{-1}.$ They further predicted that ${\dot{N}}_{t}\propto {M}_{\bullet }^{-0.25}$ for the isothermal case (also see Merritt 2013a). Stone & Metzger (2014) employed a stellar mass function, $\xi (m)$, in their DF, applied it to a sample of 200 galaxies, and obtained ${\dot{N}}_{t}\propto {M}_{\bullet }^{-0.4}.$Magorrian & Tremaine (1999) solved the steady state FP equation for an axis-symmetric nuclei with stars on a centrophobic orbits and obtained ${\dot{N}}_{t}\sim {10}^{-4}{M}_{\bullet }^{-0.19}\;{\mathrm{Yr}}^{-1}.$ The non-spherical mass distribution in the galactic center provides an additional torque to the star's orbit, which results in additional diffusion of stars in the loss cone, and such orbits are called drain orbits. Vasiliev & Merritt (2013) solved the full FP equation numerically for an axis-symmetric nuclei with a slight deviation from the sphericity and found that ${\dot{N}}_{t}$ is two to three times higher than that for spherical geometry. The orbital evolution in the galactic nuclei is more complicated in the presence of triaxial potentials, which support two distinct families of tube orbits circulating about the long and short axes of the triaxial figure (Merritt & Vasiliev 2011). In addition, there is a new class of centrophilic orbits called the pyramids, and the defining feature of the pyramid orbit is that the minimum of J = 0 and a star on such an orbit will eventually find its way into the SBH even without the assistance of collisional relaxation (Merritt & Valluri 1999). Feeding rates due to collisional loss cone refilling are very large in such galaxies compared with the spherical galaxies (Merritt & Poon 2004). Due to the complexity of the orbits in the non-spherical galaxies, we assume that the galaxy is spherical in this paper and plan to study the more general ellipsoidal or triaxial case later.

Approximately two dozen TDE candidates have been observed and found with diverse mixture of optical (van Velzen et al. 2011; Cenko et al. 2012; Gezari et al. 2012; Arcavi et al. 2014; Chornock et al. 2014), UV (Gezari et al. 2008, 2009), and X-ray (Donley et al. 2002; Komossa 2002; Maksym et al. 2013) detection. The TDE rates are highly uncertain from observations due to low sample size, which are typically a value of ${10}^{-5}\;{\mathrm{Yr}}^{-1}$ per galaxy inferred from X-ray (Donley et al. 2002), UV (Gezari et al. 2008), and optically (van Velzen & Farrar 2014) selected events. Although uncertain, the values obtained are below the theoretical estimates.

The second theoretical aspect of TDE is the accretion dynamics of the disrupted debris and the resulting luminosity and spectral energy distribution. A star tidally captured by the BH is disrupted in a dynamical timescale. Rees (1988) calculated the energy of the disrupted debris for the initial parabolic orbit of the star. The debris is assumed to follow a Keplerian orbit around the BH and the mass in-fall rate is inferred to be $\dot{M}\propto {t}^{-5/3}.$ In general, the mass in-fall rate depends on the internal structure and properties of the star, and follows the ${t}^{-5/3}$ law only at the late stages of its evolution (Phinney 1989; Lodato et al. 2009; Guillochon & Ramirez-Ruiz 2013). The debris experience stream collision either due to incoming stream that intersects with the outflowing stream at the pericenter (Kochanek 1994) or due to relativistic precession at the pericenter (Hayasaki et al. 2013). These interactions result in circularization of the debris to form an accretion disk (Hayasaki et al. 2013; Bonnerot et al. 2015; Shiokawa et al. 2015). Strubbe & Quataert (2009) proposed a simplistic steady accretion model, with the fraction of mass outflow caused by the strong radiative pressure in the super Eddington phase constant. We follow the work of Strubbe & Quataert (2009) to obtain the light curve profile and duration of the flare detection $\delta {t}_{f}$ for given instrument parameters.

The energy of the debris E_d depends on the pericenter of the star orbit ${r}_{p}(E,J,{M}_{\bullet },m).$ Hence, the mass accretion rate $\dot{M}(E,J,{M}_{\bullet },m),$ the flux, and $\delta {t}_{f}$ depends on E and J. We build an accretion model based on the initial stellar system parameters and simulate light curve profiles as a function of E and J in the optical and X-ray bands. For a DF that depends on the E only, the stars are diffused into the loss cone through star–star interactions, which leads to the change in orbital angular momentum of the star (Lightman & Shapiro 1977); thus the DF of stars within the loss cone depends on both E and J and we calculate ${\dot{N}}_{t}(E,J,m,\gamma )$ and the detectable rates of TDE ${\dot{N}}_{D}$ for the various optical and X-ray surveys.

The crucial point is that J plays an important role in the stellar dynamical process through ${\dot{N}}_{t}(E,J,m,\gamma )$ and the accretion process through ${r}_{p}(E,J,{M}_{\bullet },m)$, which impacts the detectable TDE rates; this has not been taken into account in previous calculations.

The observed sample of candidate TDE is expanding rapidly, mainly at the optical frequencies due to the advent of the highly sensitive wide field surveys such as Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), observing in both the medium deep survey (MDS) mode and 3π survey mode (Kaiser et al. 2002). The study of TDE will be further revolutionized in the next decade by the Large Synoptic Survey Telescope (LSST) with high sensitivity in optical frequencies (LSST Science Collaboration et al. 2009) and the extended Roentgen Survey with an Imaging Telescope Array (eROSITA) in the X-ray band, which performs an all sky survey (ASS) twice a year (Merloni et al. 2012) and can detect hundreds or thousands of TDE per year (Gezari et al. 2008; Khabibullin et al. 2014; van Velzen & Farrar 2014). We calculate the detectable rate ${\dot{N}}_{D}$ for Pan-STARRS 3π, Pan-STARRS MDS, and LSST in the optical g band and eROSITA in soft X-ray band; instrument details are given in Table 2.

The Figure 1 shows the methodology adopted to calculate the detectable TDE with the initial parameters M_•, ${M}_{\star },{R}_{\star }$, E, J, and redshift z. In Section 2, we solve the steady state FP equation for a power law density profile $\rho (r)\propto {r}^{-\gamma }$ and stellar mass function $\xi (m).$ We obtain ${\dot{N}}_{t}(E,J,m,\gamma )$ for the typical parametric range of density profiles (γ = 0.6–1.4), energy, and angular momentum ($J\leqslant {J}_{{lc}}$), which we use later to calculate the detection rate, ${\dot{N}}_{D}.$ In Section 3, we calculate the energy E_d of the disrupted debris and the maximum radius ${R}_{l}(E,J)$ from the star center to the point where the debris is bound to the BH. We then simulate $\dot{M}(E,J,t).$ In Section 4, we compare the accretion t_a, ring formation t_r, viscous t_v, and radiation t_R timescales and discuss the conditions for the formation of an accretion disk. In Section 5, by equating the $\dot{M}$ and Eddington mass accretion rate ${\dot{M}}_{E},$ we obtain the critical BH mass ${M}_{c}(E,J)$, such that for ${M}_{\bullet }\lt {M}_{c},$ the accretion disk formed has a super Eddington phase. We then simulate the light curve profiles in the optical and X-ray as a function of E, J, and M_•, depending on whether the accretion disk is super Eddington or sub Eddington. The flux from the source at a redshift z is compared with the sensitivity of the mission instrument to obtain $\delta {t}_{f}.$ In Section 6, we calculate ${\dot{N}}_{D}$ for optical and X-ray missions by assuming the standard cosmological parameters and BH mass function to obtain the galaxy density. Table 1 shows a glossary of symbols we use in this paper.

Zoom In Zoom Out Reset image size

Table 1.  Glossary of Symbols

Common Parameters
M•	Black hole mass	J	Orbital angular momentum
M6	${M}_{\bullet }/({10}^{6}{M}_{\odot })$	Jlc	Loss cone angular momentum
M⋆	Stellar mass	Jc	Angular momentum of circular orbit
E	Orbital energy	ℓ	$J/{J}_{{lc}}$
$\xi (m)$	Stellar mass function	m	${M}_{\star }/{M}_{\odot }$
Stellar Dynamical Parameters
j	${J}^{2}/{J}_{c}^{2}$	${\mathcal{E}}$	$E/{\sigma }^{2}$
rt	Tidal radius	ρ	Galactic density
rh	Black hole influence radius	γ	Galaxy density power law
s	$r/{r}_{h}$	st	${r}_{t}/{r}_{h}$
Rs	Schwarzschild radius	${{\rm{\Phi }}}_{\star }$	Stellar potential
rb	Break radius of Nuker profile	${{\rm{\Phi }}}_{\bullet }$	Black hole potential
σ	stellar velocity dispersion	Φ	Total potential = ${{\rm{\Phi }}}_{\bullet }+{{\rm{\Phi }}}_{\star }$
q	Diffusion parameter	tMS	Main sequence lifetime
${{\mathcal{E}}}_{c}$	Critical energy for q = 1	Tr	Radial period of orbit
${\dot{N}}_{t}$	Theoretical TDE rate	f⋆	Probability of main sequence star capture
Accretion Dynamical Parameters
rp	Pericenter of the orbit	td	dynamical time
$\bar{e}$	$E/({{GM}}_{\bullet }/{r}_{t})$ = $({r}_{t}/{r}_{h}){\mathcal{E}}$	Γ	Adiabatic index
Ed	Energy of disrupted debris	k	Spin factor
Rl	Maximum radius from star center to bound debris	tm	Orbital period of inner-most debris
xl	${R}_{l}/{R}_{\star }$	ta	Accretion timescale
x	${\rm{\Delta }}R/{R}_{\star }$	${\rm{\Delta }}R$	Debris radius from star center
	${E}_{d}/{E}_{{dm}}$	Edm	Energy of inner-most bound debris
μ	$M/{M}_{\star }$	M	Debris mass with energy Ed
$\dot{M}$	Mass accretion rate	τ	$t/{t}_{m}$
${\dot{M}}_{E}$	Eddington mass accretion rate	tr	Ring formation timescale
fr	Fraction of star mass bound to black hole	tv	Viscous timescale
Mc	Critical black hole mass	tR	Radiation timescale
κ	Opacity of the medium	${{\mathcal{T}}}_{r}$	${t}_{r}/{t}_{a}$
rc	Circularization radius	${{\mathcal{T}}}_{v}$	${t}_{v}/{t}_{R}$
rL	Outflowing wind launch radius	Tph	Temperature of photosphere
rph	Radius of photosphere of outflowing wind	Te	Effective temperature of disk
Le	Luminosity emitted from the source	L	Luminosity
$\psi ({M}_{\bullet })$	Black hole mass function	LE	Eddington luminosity
$P({M}_{\bullet },z)$	Probability of detection	z	Redshift
ϒ	Detection efficiency of a detector	${\dot{N}}_{D}$	Detectable rate
Instrumental Parameters
fl	sensitivity of the detector	tcad	Cadence of instrument
tint	Integration time of detector	fs	Fraction of sky survey

Download table as:  ASCIITypeset image

The strength of tidal encounter is expressed in terms of parameter ${\eta }_{t},$ the ratio of surface gravity of the star ${{GM}}_{\star }/{R}_{\star }^{2}$ and tidal gravity due to the BH ${{GM}}_{\bullet }{R}_{\star }/{r}_{p}^{3}$ at the pericenter r_p, given by (Merritt 2013b)

$$\begin{eqnarray}&&{\eta }_{t}\;=\;{\left(\displaystyle \frac{{{GM}}_{\star }}{{R}_{\star }^{2}}\displaystyle \frac{{r}_{p}^{3}}{{{GM}}_{\bullet }{R}_{\star }}\right)}^{\frac{1}{2}}.\end{eqnarray} \tag{ 1 }$$

The tidal radius r_t is defined as the value of r_p that satisfies

$$\begin{eqnarray}&&{r}_{t}={\left({\eta }_{t}^{2}\displaystyle \frac{{M}_{\bullet }}{{M}_{\star }}\right)}^{\frac{1}{3}}{R}_{\star }\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&=\;2.25\times {10}^{-6}{\eta }_{t}^{\frac{2}{3}}{\left(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }}\right)}^{\frac{1}{3}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-\frac{1}{3}}\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)\;{\rm{pc}}.\end{eqnarray} \tag{ 2b }$$

The quantity ${\eta }_{t}$ is the form factor of order unity and we have taken ${\eta }_{t}\;=\;1.$ For a main sequence star with the mass-radius relation ${R}_{\star }\;=\;{R}_{\odot }{m}^{n}$, where $m\;=\;{M}_{\star }/{M}_{\odot },$ the tidal radius given by Equation (2b) reduces to

$$\begin{eqnarray}&&{r}_{t}({M}_{\bullet },\;m)\approx 2.25\times {10}^{-6}{\left(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }}\right)}^{\frac{1}{3}}{m}^{n-\frac{1}{3}}\;{\rm{pc}}\end{eqnarray} \tag{ 3 }$$

and $n\gt 1/3,\;$ ${r}_{t}({M}_{\bullet },m)$ increases with m. We take n = 0.8 for the entire range of stellar masses in our calculations (Kippenhahn & Weigert 1994). We take the lifetime of the main sequence star ${t}_{\mathrm{MS}}\propto {M}_{\star }^{-2.5}$ and the dynamical time of the star to fall into the BH ${t}_{\mathrm{dyn}}\;=\;\sqrt{{a}^{3}/{{GM}}_{\bullet }}$, where a is the semimajor axis of the star to the BH. For a star to be captured during its main sequence lifetime, ${t}_{\mathrm{dyn}}\lt {t}_{\mathrm{MS}}$, which gives

$$\begin{eqnarray}&&{m}_{l}\;=\;{\left({t}_{\odot }\sqrt{\displaystyle \frac{{{GM}}_{\bullet }}{{a}^{3}}}\right)}^{0.4},\end{eqnarray} \tag{ 4 }$$

where t_⊙ is the lifetime of the sun and is shown in Figure 2 for $a\;=\;{r}_{h}.$

Zoom In Zoom Out Reset image size

Stars in the galactic center move in the potential field of both SMBH and other stars in the galaxy. The DF is assumed to be a function of energy $E\;=\;{\rm{\Phi }}(r)-{v}^{2}/2$ only and is given by $f(E)\propto {E}^{p}$ for $r\leqslant {r}_{h}$ and ${\rm{\Phi }}(r)\;=\;{{GM}}_{\bullet }/r$, where ${r}_{h}\;=\;{{GM}}_{\bullet }/{\sigma }^{2}$ is the radius of influence and σ is the stellar velocity dispersion and is related to the M_• through the ${M}_{\bullet }-\sigma$ relation given by Ferrarese & Ford (2005)

$$\begin{eqnarray}&&{M}_{\bullet }\;=\;1.66\times {10}^{8}{M}_{\odot }{\left(\displaystyle \frac{\sigma }{200\;\mathrm{Km}\;{{\rm{s}}}^{-1}}\right)}^{4.86}.\end{eqnarray} \tag{ 5 }$$

Bahcall & Wolf (1976) introduced stellar scattering and diffusion and found that $p\;=\;3/4$ (Peebles 1972) gives a negatively divergent flux; they also obtained $p\;=\;1/4$ for the steady state distribution, which gives a constant energy flux. The stars are tidally captured if the angular momentum is $J\leqslant {J}_{{lc}}(E,{r}_{t})$, where ${J}_{{lc}}(E,{r}_{t})\;=\;\sqrt{2{r}_{t}^{2}({\rm{\Phi }}({r}_{t})-E)}$ is the loss cone angular momentum (Frank & Rees 1976). The maximum value of J is ${J}_{{lc}}(E,{r}_{t}).$ As ${J}_{{lc}}(E,{r}_{t})\geqslant 0,$ the maximum value of energy is ${E}_{{\boldsymbol{m}}}\;=\;{\rm{\Phi }}({r}_{t}).$

Because ${J}_{{lc}}(E,{r}_{t})$ depends on r_t, which varies with M_⋆, we consider a DF that depends on the stellar mass function $\xi (m)$ given by (Kroupa 2001)

$$\begin{eqnarray}\xi (m)\approx \left\{\begin{array}{ll}{{Hm}}^{-1.3} & 0.08\lt m\lt 0.5\\ {{Bm}}^{-2.3} & 0.5\lt m\lt 150\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray*}&&H\;=\;2B,B\;=\;\displaystyle \frac{1}{7.91-0.77\;{m}_{m}^{-1.3}}\;{\rm{and}}\;m\;=\;\displaystyle \frac{M}{{M}_{*}},\end{eqnarray*}$$

where m_m is the maximum mass of a main sequence star in the stellar distribution, taken to be 150.

Wang & Merritt (2004) and Stone & Metzger (2014) have taken the sample of galaxies with a Nuker law, which is basically a double power law profile with break radius r_b that separates inner and outer slopes. In Figure 3 we have shown the range of r_b and r_h for the sample of galaxies given in Wang & Merritt (2004). For most of the galaxies ${r}_{b}\gt {r}_{h},$ which is also true for the sample of galaxies given in Stone & Metzger (2014). Because the stellar dynamics in the galactic center is influenced by the BH for $r\leqslant {r}_{h},$ we consider a single power law density profile $\rho (r)\;=\;{\rho }_{0}{(r/{r}_{0})}^{-\gamma }$ for $r\leqslant {r}_{h}$, where γ is the inner slope of the Nuker law. We define r_h, where ${M}_{\star }({r}_{h})\;=\;2{M}_{\bullet }$ (Wang & Merritt 2004), such that

$$\begin{eqnarray}&&{\rho }_{0}{r}_{0}^{\gamma }\;=\;\displaystyle \frac{3-\gamma }{2\pi }{M}_{\bullet }{r}_{h}^{\gamma -3}.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

The potential due to the stellar distribution is obtained from the Poisson equation and given by

$$\begin{eqnarray}{{\rm{\Phi }}}_{\star }(r)=2{\sigma }^{2}\left\{\begin{array}{ll}\displaystyle \frac{1}{2-\gamma }\left[1-{\left(\displaystyle \frac{r}{{r}_{h}}\right)}^{2-\gamma }\right] & \gamma \ne 2\\ \mathrm{ln}(\displaystyle \frac{{r}_{h}}{r}) & \gamma =2.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

The total potential is given by ${\rm{\Phi }}(r)={{\rm{\Phi }}}_{\bullet }(r)+{{\rm{\Phi }}}_{\star }(r)$, where ${{\rm{\Phi }}}_{\bullet }(r)={{GM}}_{\bullet }/r$ is the potential due to the BH. We consider the DF as

$$\begin{eqnarray}&&F(E,m)=f(E)\xi (m).\end{eqnarray} \tag{ 9 }$$

The density of stars for $F(E,m)$ is given by

$$\begin{eqnarray}\rho (r) & = & \displaystyle \int {d}^{3}{{vM}}_{\star }F(E,m)\\ & = & \displaystyle \int {d}^{3}{vf}(E)\displaystyle \int {dm}\xi (m){M}_{\star },\end{eqnarray} \tag{ 10 }$$

where ${d}^{3}v=2\pi {v}_{t}{{dv}}_{t}{{dv}}_{r}$ with radial velocity v_r and tangential velocity v_t is given by

$$\begin{eqnarray}{v}_{t} & = & \displaystyle \frac{J}{r};\\ {v}_{r} & = & \sqrt{2({\rm{\Phi }}(r)-E)-\displaystyle \frac{{J}^{2}}{{r}^{2}}};\\ {d}^{3}v & = & \displaystyle \frac{2\pi }{{r}^{2}{v}_{r}}{JdEdJ}.\end{eqnarray} \tag{ 11 }$$

For a spherically isotropic galactic center, the function f(E) is obtained through the inverse transform of Equation (10), and is known as the Eddington formula, given by (Binney & Tremaine 2008)

$$\begin{eqnarray}&&f(E)=\displaystyle \frac{1}{\sqrt{8}{\pi }^{2}\left\langle{M}_{\star }\right\rangle}\displaystyle \frac{d}{{dE}}{\displaystyle \int }_{{E}_{\mathrm{min}}}^{E}\displaystyle \frac{d\rho }{d{\rm{\Phi }}}\displaystyle \frac{1}{\sqrt{E-{\rm{\Phi }}}}\;d{\rm{\Phi }},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&\left\langle{M}_{\star }\right\rangle={\displaystyle \int }_{0.08}^{150}{M}_{\star }\xi (m)\;{dm}\end{eqnarray} \tag{ 13 }$$

and E_min is the minimum of E taken to be −100. The number of stars in the cluster for a given $F(E,m)$ is

$$\begin{eqnarray}&&N=\int {d}^{3}r\int {d}^{3}v\int {dmF}(E,m),\end{eqnarray} \tag{ 14 }$$

where ${d}^{3}r=4\pi {r}^{2}{dr}$ for spherical galaxy. In terms of dimensionless variables ${\ell }=J/{J}_{{lc}}$ and ${\mathcal{E}}=E/{\sigma }^{2},$ the number of stars is given by

$$\begin{eqnarray}N({\mathcal{E}},{\ell })d{\mathcal{E}}d{\ell } & = & \displaystyle \int {dm}\;8{\pi }^{2}{J}_{{lc}}^{2}{\sigma }^{2}{\ell }{T}_{r}({\mathcal{E}},{\ell })\\ & & \times F({\mathcal{E}},m)d{\mathcal{E}}d{\ell }\end{eqnarray} \tag{ 15 }$$

where ${T}_{r}({\mathcal{E}},{\ell })=\oint \frac{{dr}}{{v}_{r}}$ is the radial period. For a spherical geometry, ${T}_{r}({\mathcal{E}},{\ell })$ is a function of ${\mathcal{E}}$ only, and the $N({\mathcal{E}},{\ell })$ increases with ${\ell }.$

From the given stellar density profile and potential ${\rm{\Phi }}(r),$ the DF f(E) is given by

$$\begin{eqnarray}&&f({\mathcal{E}})=\displaystyle \frac{{M}_{\bullet }}{\left\langle{M}_{\star }\right\rangle}\displaystyle \frac{1}{{r}_{h}^{3}{\sigma }^{3}}g({\mathcal{E}})\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&g({\mathcal{E}})=\displaystyle \frac{1}{4\sqrt{2}{\pi }^{3}}\\ &&\times \left\{\begin{array}{l}\gamma (3-\gamma )\frac{d}{d{\mathcal{E}}}{\int }_{{s}_{1}}^{{s}_{2}}\\ \;\times \frac{{s}^{\gamma -1}}{\sqrt{{\mathcal{E}}-s-\frac{2}{2-\gamma }(1-{s}^{\gamma -2})}}\;{ds} & \gamma \ne 2\\ {\int }_{{{\mathcal{E}}}_{\mathrm{min}}}^{{\mathcal{E}}}\frac{{\mathcal{L}}(2+{\mathcal{L}})}{{(1+{\mathcal{L}})}^{3}}\frac{1}{\sqrt{{\mathcal{E}}-{\rm{\Psi }}}}\;d{\rm{\Psi }} & \gamma =2\end{array}\right.\end{eqnarray} \tag{ 17 }$$

and ${\mathcal{L}}$ is the Lambert function given by ${\mathcal{L}}{e}^{{\mathcal{L}}}=(1/2){e}^{{\rm{\Psi }}/2}$, where ${\rm{\Psi }}={\rm{\Phi }}/{\sigma }^{2},\;$ $s={r}_{h}/r,\;$ ${\mathcal{E}}=E/{\sigma }^{2}$ (Wang & Merritt 2004) and s₁ and s₂ are obtained by solving

$$\begin{eqnarray}&&{s}_{1}+\displaystyle \frac{2}{2-\gamma }(1-{s}_{1}^{\gamma -2})={{\mathcal{E}}}_{\mathrm{min}}\end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray}&&{s}_{2}+\displaystyle \frac{2}{2-\gamma }(1-{s}_{2}^{\gamma -2})={\mathcal{E}}.\end{eqnarray} \tag{ 18b }$$

Figure 4 shows the plot of $g({\mathcal{E}})$ for various γ. For $\gamma =2,\;$ $g({\mathcal{E}})$ corresponds to Equation (17) of Wang & Merritt (2004). The BH potential dominates over the star potential for $r\ll {r}_{h}$ as shown in Figure 5 and thus $g({\mathcal{E}})\propto {{\mathcal{E}}}^{\gamma -3/2}$ for ${\mathcal{E}}\gg 1.$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. Loss Cone Theory

The loss cone is a geometrical region in phase space for which ${r}_{p}\leqslant {r}_{t}.$ We adopt the Cohn & Kulsrud (1978) formalism for computing the flux of stars into the loss cone in ${\mathcal{E}}$ and the $j={J}^{2}/{J}_{c}^{2}({\mathcal{E}})$ phase space, where ${J}_{c}({\mathcal{E}})$ is the angular momentum of circular orbit with energy ${\mathcal{E}}.$ We consider the diffusion in the  j space only and in the limit $j\to 0,$ the steady state FP equation reduces to (Merritt 2013b)

$$\begin{eqnarray}&&\displaystyle \frac{d{\mathcal{F}}(\chi ,y)}{d\chi }=\displaystyle \frac{d}{{dy}}\left(y\displaystyle \frac{d{\mathcal{F}}(\chi ,y)}{{dy}}\right)\end{eqnarray} \tag{ 19 }$$

with the boundary condition

$$\begin{eqnarray}&&{\mathcal{F}}(0,y)=0\forall y\lt {y}_{{lc}}\end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray}&&{\mathcal{F}}(0,y)={\mathcal{F}}(1,y)\forall y\geqslant {y}_{{lc}}\end{eqnarray} \tag{ 20b }$$

where

$$\begin{eqnarray}&&\chi =\displaystyle \frac{1}{\left\langle({\mathcal{E}})\right\rangle}{\displaystyle \int }_{{r}_{p}}^{r}\underset{j\to 0}{\mathrm{lim}}\displaystyle \frac{\left\langle{({\rm{\Delta }}j)}^{2}\right\rangle}{2j}\;\displaystyle \frac{{dr}}{{v}_{r}}\;{\rm{and}}\;y=\displaystyle \frac{j}{\left\langle({\mathcal{E}})\right\rangle}\end{eqnarray} \tag{ 21 }$$

and $\left\langle({\mathcal{E}})\right\rangle=\oint {\mathrm{lim}}_{j\to 0}\frac{\left\langle{({\rm{\Delta }}j)}^{2}\right\rangle}{2j}\;\frac{{dr}}{{v}_{r}}$ is the orbit averaged angular momentum diffusion coefficient and y_lc is y at $j={j}_{{lc}}.$  Merritt (2015b) has expressed ${\mathcal{F}}(\chi ,y)$ in terms of the distribution of the pericenters r_p and apocenters r_a, and calculated the capture rate ${\dot{N}}_{t}$ in terms of r_a and r_p, while we use appropriately scaled values of E and J² for calculating r_p and various other parameters required in the gas dynamical calculation in the following sections. Note that in the calculation we use the total potential inclusive of the stars and the BH.

The local diffusion coefficient in the limit $j\to 0$ is given by (Magorrian & Tremaine 1999)

$$\begin{eqnarray}&&\underset{j\to 0}{\mathrm{lim}}\displaystyle \frac{\left\langle{({\rm{\Delta }}j)}^{2}\right\rangle}{2j}=\displaystyle \frac{{r}^{2}}{{J}_{c}^{2}({\mathcal{E}})}\left\langle{({\rm{\Delta }}{v}_{\perp })}^{2}\right\rangle\end{eqnarray} \tag{ 22 }$$

where $\left\langle{({\rm{\Delta }}{v}_{\perp })}^{2}\right\rangle$ is given in Appendix L of Binney & Tremaine (2008). Thus, the orbit averaged diffusion coefficient is given by

$$\begin{eqnarray}\left\langle({\mathcal{E}})\right\rangle & = & \displaystyle \frac{32\sqrt{2}}{3}\displaystyle \frac{{\pi }^{2}{G}^{2}\left\langle{M}_{f}^{2}\right\rangle\mathrm{log}{\rm{\Lambda }}}{{J}_{c}^{2}({\mathcal{E}})}\\ & & \times \displaystyle \frac{{M}_{\bullet }}{\left\langle{M}_{\star }\right\rangle}\displaystyle \frac{1}{{\sigma }^{2}}(2{h}_{1}({\mathcal{E}})+3{h}_{2}({\mathcal{E}})-{h}_{3}({\mathcal{E}}))\end{eqnarray} \tag{ 23 }$$

where M_f is the mass of field star with the maximum mass taken to be 150 M_⊙, ${\rm{\Lambda }}\approx {M}_{\bullet }/{M}_{\star }$, and

$$\begin{eqnarray}&&\left\langle{M}_{f}^{2}\right\rangle={M}_{\odot }^{2}{\displaystyle \int }_{0.08}^{150}{m}_{f}^{2}\xi ({m}_{f})\;{{dm}}_{f}\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}&&{h}_{1}({\mathcal{E}})={\displaystyle \int }_{0}^{s({\mathcal{E}})}{ds}^{\prime} \displaystyle \frac{{s}^{\prime 2}}{\sqrt{{\rm{\Psi }}({s}^{\prime })-{\mathcal{E}}}}{\displaystyle \int }_{-\infty }^{{\mathcal{E}}}d{{\mathcal{E}}}^{\prime }g({{\mathcal{E}}}^{\prime })\end{eqnarray} \tag{ 24b }$$

$$\begin{eqnarray}&&{h}_{2}({\mathcal{E}})={\displaystyle \int }_{0}^{s({\mathcal{E}})}{{ds}}^{\prime }\displaystyle \frac{{s}^{\prime 2}}{{\rm{\Psi }}({s}^{\prime })-{\mathcal{E}}}{\displaystyle \int }_{{\mathcal{E}}}^{{\rm{\Psi }}({s}^{\prime s})}d{{\mathcal{E}}}^{\prime }\sqrt{{\rm{\Psi }}({s}^{\prime })-{{\mathcal{E}}}^{\prime }}g({{\mathcal{E}}}^{\prime })\end{eqnarray} \tag{ 24c }$$

$$\begin{eqnarray}{h}_{3}({\mathcal{E}}) & = & {\displaystyle \int }_{0}^{s({\mathcal{E}})}{{ds}}^{\prime }\displaystyle \frac{{s}^{\prime 2}}{{({\rm{\Psi }}({s}^{\prime })-{\mathcal{E}})}^{2}}\\ & & \times {\displaystyle \int }_{{\mathcal{E}}}^{{\rm{\Psi }}({s}^{\prime })}d{{\mathcal{E}}}^{\prime }{({\rm{\Psi }}({s}^{\prime })-{{\mathcal{E}}}^{\prime })}^{\frac{3}{2}}g({{\mathcal{E}}}^{\prime })\end{eqnarray} \tag{ 24d }$$

where $s({\mathcal{E}})$ is obtained by solving ${\rm{\Psi }}(s)={\mathcal{E}},\;$ $s=r/{r}_{h}$ and

$$\begin{eqnarray}&&{\rm{\Psi }}(s)=\displaystyle \frac{{\rm{\Phi }}}{{\sigma }^{2}}=\displaystyle \frac{1}{s}+\displaystyle \frac{2}{2-\gamma }(1-{s}^{2-\gamma }).\end{eqnarray} \tag{ 25 }$$

Now ${J}_{c}^{2}({\mathcal{E}})$ is given by

$$\begin{eqnarray}&&{J}_{c}^{2}({\mathcal{E}})={\sigma }^{2}{r}_{h}^{2}[{s}_{c}({\mathcal{E}})+2\;{s}_{c}^{4-\gamma }({\mathcal{E}})]\end{eqnarray} \tag{ 26 }$$

where ${s}_{c}({\mathcal{E}})$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{2{s}_{c}}+\displaystyle \frac{2}{2-\gamma }(1-{s}_{c}^{2-\gamma })-{s}_{c}^{2-\gamma }={\mathcal{E}}.\end{eqnarray} \tag{ 27 }$$

The solution of Equation (19) with the boundary conditions (20b) is given by (Merritt 2013b)

$$\begin{eqnarray}&&{\mathcal{F}}(\chi ,\;y)=X({j}_{{lc}})\left(1-2\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{e}^{\frac{-{\alpha }_{n}^{2}q}{4}\chi }}{{\alpha }_{n}}\displaystyle \frac{{J}_{0}\left({\alpha }_{n}\sqrt{y/{y}_{{lc}}}\right)}{{J}_{1}({\alpha }_{n})}\right)\end{eqnarray} \tag{ 28 }$$

where ${\alpha }_{n}$ are the consecutive zeros of the Bessel function ${J}_{0}(\alpha ),\;$ $q=1/{y}_{{lc}}$, and $X({j}_{{lc}})$ is given by

$$\begin{eqnarray}X({j}_{{lc}}) & = & \displaystyle \frac{f({\mathcal{E}})}{1+{q}^{-1}\zeta (q)\mathrm{log}\left(\displaystyle \frac{1}{{j}_{{lc}}}\right)};\\ \zeta (q) & = & 1-4\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{e}^{\frac{-{\alpha }_{n}^{2}q}{4}}}{{\alpha }_{n}^{2}}\end{eqnarray} \tag{ 29 }$$

where $f({\mathcal{E}})$ is given by Equation (16). The Equation (28) gives the DF of stars in the loss cone in terms of energy ${\mathcal{E}}$ and angular momentum ${\ell }=J/{J}_{{lc}}=\sqrt{y/{y}_{{lc}}}.$ Figure 6 shows the plot of ${\mathcal{F}}(0,{\ell })$ for 1000 terms in the summation in Equation (28), which matches with the boundary condition given by Equation (20a) with an accuracy of $\sim {10}^{-3}.$ With an increase in the number of terms in summation, the order of accuracy of ${\mathcal{F}}(0,{\ell })$ to Equation (20b) increases; however, we use 1000 terms in summation to calculate ${\mathcal{F}}(\chi ,{\ell })$ because it is in close agreement with the boundary condition.

Zoom In Zoom Out Reset image size

The approximation to $\zeta (q)$ given by Cohn & Kulsrud (1978) is

$$\begin{eqnarray}\zeta (q)\approx {\zeta }_{{CK}}(q)=\left\{\begin{array}{ll}1 & q\gg 1\\ \displaystyle \frac{q}{0.186q+0.824\sqrt{q}} & q\ll 1.\end{array}\right.\end{eqnarray} \tag{ 30 }$$

We compared the $\zeta (q)$ obtained for a summation of 10,000 terms in Equation (29) with ${\zeta }_{{CK}}(q)$, and it does not fit very well for q close to unity, as shown in Figure 7. Thus we used a better approximation to $\zeta (q)$ given by

$$\begin{eqnarray}\zeta (q)\approx \left\{\begin{array}{ll}1 & q\geqslant 4\\ \displaystyle \frac{q}{\begin{array}{l}0.86{q}^{0.5}+0.384q-0.379{q}^{1.5}\\ \;\;+\;0.427{q}^{2}-0.095{q}^{2.5}\end{array}} & q\lt 4\end{array}\right.\end{eqnarray} \tag{ 31 }$$

which follows Equation (30) for $q\ll 1$ and gives a better fit to $\zeta (q)$ for q close to unity. The residual for Equation (30) is higher than the residual for Equation (31), as shown in Figure 7.

Zoom In Zoom Out Reset image size

The function $q({\mathcal{E}})$ given by

$$\begin{eqnarray}&&q({\mathcal{E}})=\displaystyle \frac{\left\langle({\mathcal{E}})\right\rangle}{{j}_{{lc}}}=\displaystyle \frac{\left\langle({\mathcal{E}})\right\rangle{J}_{c}^{2}({\mathcal{E}})}{{J}_{{lc}}^{2}({\mathcal{E}},{r}_{t})}=\displaystyle \frac{\left\langle({\mathcal{E}})\right\rangle{J}_{c}^{2}({\mathcal{E}})}{2{r}_{t}^{2}({\rm{\Phi }}({r}_{t})-{\mathcal{E}}{\sigma }^{2})}\end{eqnarray} \tag{ 32 }$$

can be interpreted as the ratio of the orbital period to the timescale for diffusional refilling of the loss cone. The regime $q({\mathcal{E}})\gt 1$ defines the pinhole or full loss cone in which stellar encounters replenish loss cone orbits much more rapidly than they are depleted, whereas $q({\mathcal{E}})\lt 1$ defines the diffusive or empty loss cone regime. The Figure 8 shows $q({\mathcal{E}})$ plotted as a function of ${\mathcal{E}}$ for $\gamma =1.0.$ The function $q({\mathcal{E}})$ decreases with ${\mathcal{E}}$, which implies that the high energy orbits have smaller diffusion angle. The smaller the diffusion angle, the higher the diffusion time, and thus the lower the feeding rate to the loss cone. The critical energy ${{\mathcal{E}}}_{c}$ defined by $q({{\mathcal{E}}}_{c})=1$ decreases with M_•, and m is shown in Figure 9 for $\gamma =1.$ The ${{\mathcal{E}}}_{c}$ is the energy from which the majority of the loss cone flux originates (Lightman & Shapiro 1977). With an increase in M_•, the relaxation time of the galaxy increases; thus the diffusion timescale increases (Frank & Rees 1976) and q decreases, which results in a decrease in ${{\mathcal{E}}}_{c}.$ As ${r}_{t}({M}_{\bullet },m)\propto {m}^{n-1/3},$ it increases with m for $n=0.8,$;  J_lc increases so q decreases, which results in a decrease of ${{\mathcal{E}}}_{c}.$ As γ increases, the diffusion timescale decreases due to an increase in the number of scatterers, and thus q and ${{\mathcal{E}}}_{c}$ increase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using Equation (28) and the mass function of stars in the galaxy $\xi (m)$ given by Equation (6), the loss cone feeding rate is given by (Merritt 2013b)

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\dot{N}}_{t}}{{dEdy}}=4{\pi }^{2}\int \;{dm}\xi (m){J}_{c}^{2}(E)\left\langle(E)\right\rangle{\mathcal{F}}(\chi =1,y).\end{eqnarray} \tag{ 33 }$$

The corresponding feeding rate in terms of E and j is given in Merritt (2015a).

The Jacobian of the transformation from $\{E,y\}$ space to dimensionless variables $\{{\mathcal{E}},{{\ell }}^{2}={(J/{J}_{{lc}}({\mathcal{E}},{r}_{t}))}^{2}$ = $j{({J}_{c}({\mathcal{E}})/{J}_{{lc}}({\mathcal{E}},{r}_{t}))}^{2}\}$ is given by

$$\begin{eqnarray}&&{dEdy}={\rm{Det}}\left[\displaystyle \frac{\partial (E,y)}{\partial ({\mathcal{E}},{{\ell }}^{2})}\right]d{\mathcal{E}}d{{\ell }}^{2},\\ &&{\rm{where}}\;\;\;\displaystyle \frac{\partial (E,y)}{\partial ({\mathcal{E}},{{\ell }}^{2})}=\displaystyle \frac{\partial (E,y)}{\partial (E,{J}^{2})}\cdot \displaystyle \frac{\partial (E,{J}^{2})}{\partial ({\mathcal{E}},{{\ell }}^{2})}\end{eqnarray} \tag{ 34 }$$

and the following result is obtained by calculating the product of the determinants of the two Jacobians in the above equation as

$$\begin{eqnarray}&&{dEdy}={\sigma }^{2}\displaystyle \frac{{J}_{{lc}}^{2}({\mathcal{E}},{r}_{t})}{\left\langle({\mathcal{E}})\right\rangle{J}_{c}^{2}({\mathcal{E}})}d{\mathcal{E}}d{{\ell }}^{2}.\end{eqnarray} \tag{ 35 }$$

Then, the feeding rate is given by

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\dot{N}}_{t}}{d{\mathcal{E}}d{{\ell }}^{2}{dm}}=4{\pi }^{2}{\sigma }^{2}\xi (m){J}_{{lc}}^{2}({\mathcal{E}}){\mathcal{F}}(\chi =1,{\ell })\end{eqnarray} \tag{ 36 }$$

Figure 10 shows the plot of ${\dot{N}}_{t}$ as a function of ${\mathcal{E}}$ for various ℓ for ${M}_{6}=1$ and $\gamma =0.8.$ The capture rates decreases with ${\mathcal{E}}$ because of the decrease in the diffusion coefficient $D({\mathcal{E}})$ and increases with ℓ due to the increase in $N({\mathcal{E}},{\ell })$ (see Equation (15)).

Zoom In Zoom Out Reset image size

Because ${J}_{{lc}}({r}_{t})=\sqrt{2{r}_{t}^{2}({\rm{\Phi }}({r}_{t})-E)}$ (see Section 2) this implies that $E\lt {\rm{\Phi }}({r}_{t})$ and the maximum value of $E={E}_{m}={\rm{\Phi }}({r}_{t}).$ Because ${r}_{t}({M}_{\bullet },m)/{r}_{h}\sim {10}^{-6}-{10}^{-5},$ and $r\ll {r}_{h},$ the potential is dominated by the BH potential, as shown in Figure 5, which implies that ${E}_{m}({M}_{\bullet },m)={\rm{\Phi }}({r}_{t}({M}_{\bullet },m))$ = ${{GM}}_{\bullet }/{r}_{t}({M}_{\bullet },m).$ The orbital motion of a star at the turning point of the orbit r_x, is given by

$$\begin{eqnarray}&&E={\rm{\Phi }}({r}_{x})-\displaystyle \frac{{J}^{2}}{2{r}_{x}^{2}}\end{eqnarray} \tag{ 37 }$$

where ${\rm{\Phi }}(r)={{\rm{\Phi }}}_{\bullet }(r)+{{\rm{\Phi }}}_{\star }(r),\;$ ${{\rm{\Phi }}}_{\bullet }(r)={{GM}}_{\bullet }/r$, and ${{\rm{\Phi }}}_{\star }(r)$ is given by Equation (8). In terms of dimensionless variables, ${\ell }=J/{J}_{{lc}}$ and $\bar{e}=E/{E}_{m}$, where ${E}_{m}={{GM}}_{\bullet }/{r}_{t}$, such that $\bar{e}=({r}_{t}/{r}_{h}){\mathcal{E}},$  and the Equation (37) using Equation (25) is given by

$$\begin{eqnarray}&&\displaystyle \frac{\bar{e}}{{s}_{t}}\;=\\ &&\displaystyle \frac{{s}_{x}-{{\ell }}^{2}{s}_{t}+\displaystyle \frac{2}{2-\gamma }\left[{s}_{x}^{2}(1-{s}_{x}^{2-\gamma })-{{\ell }}^{2}\;{s}_{t}^{2}(1-{s}_{t}^{2-\gamma })\right]}{{s}_{x}^{2}-{{\ell }}^{2}\;{s}_{t}^{2}}\end{eqnarray} \tag{ 38 }$$

where ${s}_{x}={r}_{x}/{r}_{h},\;$ ${s}_{t}={r}_{t}/{r}_{h}.$ Since $\bar{e}$ is a monotonically decreasing function of s_x, and both the pericenter and apocenter of the orbit should lie below r_h, the minimum value of $\bar{e}$ is at s_x = 1 for ${r}_{x}={r}_{h};$ taking $\bar{e}({s}_{x}=1)={\bar{e}}_{h},$ Equation (38) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{\bar{e}}_{h}}{{s}_{t}}=\displaystyle \frac{1-{{\ell }}^{2}{s}_{t}-\displaystyle \frac{2}{2-\gamma }{{\ell }}^{2}\;{s}_{t}^{2}(1-{s}_{t}^{2-\gamma })}{1-{{\ell }}^{2}\;{s}_{t}^{2}}.\end{eqnarray} \tag{ 39 }$$

Now ${s}_{t}={r}_{t}/{r}_{h}\sim {10}^{-5}$–${10}^{-6},$ ${\bar{e}}_{h}\simeq {r}_{t}/{r}_{h}.$ Because ${J}_{{lc}}({r}_{x})=\sqrt{2{r}_{x}^{2}({\rm{\Phi }}({r}_{x})-E)}$, this implies that $E\lt {\rm{\Phi }}({r}_{x})$ and the maximum value of $E={\rm{\Phi }}({r}_{x})\simeq {\rm{\Phi }}({r}_{t})$, which corresponds to $\bar{e}\simeq 1$ and thus $\bar{e}\lesssim 1.$ The total potential is dominated by the BH near r_t, as shown in Figure 5. After ignoring the second term in the numerator in the right-hand side of Equation (38), which is a factor ${s}_{t}={r}_{t}/{r}_{h}\simeq {10}^{-5}-{10}^{-6}\ll 1$ lower than the first, Equation (38) reduces to

$$\begin{eqnarray}&&\bar{e}=\displaystyle \frac{{x}_{x}-{{\ell }}^{2}}{{x}_{x}^{2}-{{\ell }}^{2}},\end{eqnarray} \tag{ 40 }$$

where ${x}_{x}={s}_{x}/{s}_{t}.$ If x_p is lower of the two roots of x_x, it is given by

$$\begin{eqnarray}&&{x}_{p}=\displaystyle \frac{1}{2\bar{e}}\left(1-\sqrt{1-4\bar{e}(1-\bar{e}){{\ell }}^{2}}\right).\end{eqnarray} \tag{ 41 }$$

For a star to be tidally disrupted, ${x}_{p}\lt 1$, which results in

$$\begin{eqnarray}&&(1-\bar{e})(1-{{\ell }}^{2})\gt 0\end{eqnarray} \tag{ 42 }$$

and ${x}_{p}\gt 0$ results in

$$\begin{eqnarray}&&{{\ell }}^{2}\bar{e}(1-\bar{e})\gt 0.\end{eqnarray} \tag{ 43 }$$

The Equation (43) restricts the range of $\bar{e}$ to $\bar{e}\lt 1$ and thus Equation (42) implies that ${\ell }\lt 1.$ Thus, the applicable ranges are ${\bar{e}}_{h}\lt \bar{e}\lt 1$ and $0\lt {\ell }\lt 1.$ We derived the turning points s_x by solving Equation (38) subject to the constraint ${r}_{p}\lt {r}_{t}$ and ${r}_{t}\lt {r}_{a}\lt {r}_{h}$, and verified the range for $\bar{e}$ and ℓ derived above.

While we ignored relativistic effects in the analysis above, we plan to include them in the future.

The lifetime of a main sequence star is ${t}_{\mathrm{MS}}={t}_{\odot }{m}^{-2.5}$ where ${t}_{\odot }={10}^{10}\;{\rm{years}}$ is the lifetime of a solar type star and the radial period is given by

$$\begin{eqnarray}&&{T}_{r}={\displaystyle \int }_{{r}_{p}}^{{r}_{a}}\displaystyle \frac{1}{{v}_{r}}{dr}.\end{eqnarray} \tag{ 44 }$$

In terms of $s=r/{r}_{h}$ and using Equation (25), T_r is given by

$$\begin{eqnarray}&&{T}_{r}(\bar{e},{M}_{\bullet },m)=\displaystyle \frac{{r}_{h}}{\sqrt{2}\sigma }{\displaystyle \int }_{{s}_{p}}^{{s}_{a}}\\ &&\times \displaystyle \frac{s\sqrt{{s}_{t}}}{\sqrt{\begin{array}{l}{s}_{t}\left(s-{{\ell }}^{2}{s}_{t}+\displaystyle \frac{2}{2-\gamma }\left[{s}^{2}(1-{s}^{2-\gamma })\right.\right.\\ \left.\left.-{{\ell }}^{2}\;{s}_{t}^{2}(1-{s}_{t}^{2-\gamma })\right]\right)-\bar{e}({s}^{2}-{{\ell }}^{2}\;{s}_{t}^{2})\end{array}}}{ds},\end{eqnarray} \tag{ 45 }$$

where s_a and s_p are the dimensionless apocenter and pericenter obtained by solving Equation (38). We find numerically that the radial period T_r is approximated by

$$\begin{eqnarray}&&{T}_{r}(\bar{e},{M}_{\bullet },m)\simeq \displaystyle \frac{\pi }{2\sqrt{2}}\displaystyle \frac{{r}_{h}}{\sigma }\\ \times \ \left\{\begin{array}{ll}0.57{e}^{[0.27(1.47-\frac{\bar{e}}{{s}_{t}})]} & \bar{e}\lt 1.47{s}_{t}\\ {\left(\displaystyle \frac{\bar{e}}{{s}_{t}}\right)}^{-\frac{3}{2}} & \bar{e}\geqslant 1.47{s}_{t}.\end{array}\right.\end{eqnarray} \tag{ 46 }$$

Because the BH potential dominates at high energy, the corresponding orbits are Keplerian and the radial period $\propto {\bar{e}}^{(-3/2)}.$ Using the ${M}_{\bullet }-\sigma$ relation given by Equation (5) and ${r}_{h}={{GM}}_{\bullet }/{\sigma }^{2},$ the radial period in Keplerian regime is ${T}_{r}\propto {M}_{\bullet }^{0.38}{{\mathcal{E}}}^{-3/2}$, where ${\mathcal{E}}=\bar{e}/{s}_{t}.$ The stars on the loss cone orbits are captured in the radial period timescale and the number of stars in the loss cone orbit is ${N}_{{lc}}({\mathcal{E}})d{\mathcal{E}}$ (Merritt 2013b). Thus the capture rate is approximately given by ${\dot{N}}_{t}=\int ({N}_{{lc}}({\mathcal{E}})/{T}_{r}({\mathcal{E}}))d{\mathcal{E}}\propto {M}_{\bullet }^{-0.38}$ in the regime $\bar{e}\gt 1.47{s}_{t}$, which is consistent with the average best-fit slope of $-0.3$ over the entire range of $\bar{e}$ that was found numerically.

A tidally captured star is on main sequence if its main sequence lifetime is ${t}_{\mathrm{MS}}\gt {T}_{r}$, where T_r is the radial period of the orbit. Considering all possible radial phases of the star in an orbit, the probability that a star of mass m is tidally captured as a main sequence is given by

$$\begin{eqnarray}&&{f}_{\star }(\bar{e},{M}_{\bullet },m)={\rm{Min}}\left[1,\displaystyle \frac{{t}_{\mathrm{MS}}(m)}{{T}_{r}(\bar{e},{M}_{\bullet },m)}\right].\end{eqnarray} \tag{ 47 }$$

Using Equations (28), (36), and $d{\mathcal{E}}d{{\ell }}^{2}={s}_{t}^{-1}d\bar{e}d{{\ell }}^{2},$ the capture rate is given by

$$\begin{eqnarray}\displaystyle \frac{{d}^{2}{\dot{N}}_{t}}{d\bar{e}d{{\ell }}^{2}{dm}} & = & 4{\pi }^{2}\;{s}_{t}^{-1}{\sigma }^{2}\xi (m){f}_{\star }(\bar{e},{M}_{\bullet },m)\\ & & \times {J}_{{lc}}^{2}(\bar{e}){\mathcal{F}}(\chi =1,{\ell }).\end{eqnarray} \tag{ 48 }$$

Figure 5 of Freitag & Benz (2002) gives the maximum mass of BH as function of mass m of the star that is disrupted. We observe that stars with mass $m\gt 0.8$ are tidally disrupted in the entire range of BH mass ${10}^{6}{M}_{\odot }\lt {M}_{\bullet }\lt {10}^{8}{M}_{\odot }$ and for $m\lt 0.8,$ a substantial fraction is tidally captured without disruption. Thus we take the effective star mass range to be $0.8\lt m\lt 150.$ The net capture rate is given by

$$\begin{eqnarray}{\dot{N}}_{t}(\gamma ,\;{M}_{\bullet }) & = & 4{\pi }^{2}{\displaystyle \int }_{0.8}^{150}{dm}{\displaystyle \int }_{{\bar{e}}_{h}}^{1}d\bar{e}{\displaystyle \int }_{0}^{1}d{{\ell }}^{2}{\sigma }^{2}\;{s}_{t}^{-1}\\ & & \times ({M}_{\bullet },\;m)\xi (m){f}_{\star }(\bar{e},{M}_{\bullet },m){J}_{{lc}}^{2}(\bar{e})\\ & & \times {\mathcal{F}}(\chi =1,{\ell }).\end{eqnarray} \tag{ 49 }$$

We solved Equation (28) to obtain ${\mathcal{F}}(1,\;{\ell })$ and used it in Equation (48) to calculate the capture rate. The integration of Equation (48) over the energy range ${\bar{e}}_{h}\lt \bar{e}\lt 1$ and angular momentum range $0\lt {\ell }\lt 1$ results in a capture rate per unit mass $d{\dot{N}}_{t}/{dm}$, which is a decreasing function of m as $\xi (m)$ decreases with m and is shown in Figure 11(a) for various ${M}_{\bullet }={10}^{6}{M}_{\odot }{M}_{6}$ and $\gamma =1.0.$ Similarly, integrating Equation (48) over energy ${\bar{e}}_{h}\lt \bar{e}\lt 1$ and test star mass $0.8\lt m\lt 150$ results in $d{\dot{N}}_{t}/d{\ell }$, which is an increasing function ℓ as shown in Figure 11(b).

Zoom In Zoom Out Reset image size

The net ${\dot{N}}_{t}$ obtained using Equations (28) and (49) increases with γ and decreases with M₆ as shown in Figure 12. For ${M}_{6}\geqslant 10,$ ${\dot{N}}_{t}\propto {M}_{6}^{\beta }$ where $\beta =0.3\pm 0.01$, as shown in Figure 12(b) for various γ. The increase with γ is nonlinear as shown in Figure 13 and for ${\dot{N}}_{t}\propto {\gamma }^{p},$ the best-fit value of $p\sim 2.1.$ The galaxies with larger γ posses higher rates because their denser central stellar populations naturally feature shorter relaxation times and faster rates of energy and angular momentum diffusion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The classical description of a TDE was outlined by Rees (1988). In this picture, a star on parabolic orbit is tidally captured and disrupted at the pericenter and the distribution of mass of disrupted debris with respect to specific binding energy ${dM}/{{dE}}_{d}$ is roughly flat, where E_d is the energy of the disrupted debris. For stars on a parabolic orbit, Lodato et al. (2009) found that ${dM}/{{dE}}_{d}$ depends on the properties of the star and adiabatic index Γ. Using Equation (41), the pericenter is given by

$$\begin{eqnarray}{r}_{p}(\bar{e},{\ell },{M}_{\bullet },m) & = & \displaystyle \frac{{r}_{t}}{2\bar{e}}\left(1-\sqrt{1-4\bar{e}(1-\bar{e}){{\ell }}^{2}}\right)\\ & = & {r}_{t}\left(\displaystyle \frac{2{{\ell }}^{2}(1-\bar{e})}{1+\sqrt{1-4{{\ell }}^{2}\bar{e}(1-\bar{e})}}\right).\end{eqnarray} \tag{ 50 }$$

Equivalently

$$\begin{eqnarray}{r}_{p}(E,J,{M}_{\bullet }) & = & \displaystyle \frac{{{GM}}_{\bullet }}{2E}\left[1-\sqrt{1-\displaystyle \frac{2{{EJ}}^{2}}{{G}^{2}\;{M}_{\bullet }^{2}}}\right]\\ & = & \displaystyle \frac{{J}^{2}}{{{GM}}_{\bullet }}{\left[1+\sqrt{1-\displaystyle \frac{2{{EJ}}^{2}}{{G}^{2}\;{M}_{\bullet }^{2}}}\right]}^{-1}.\end{eqnarray} \tag{ 51 }$$

The stars on the loss cone orbits are captured within the dynamical time ${t}_{d}={({r}_{p}^{3}/{{GM}}_{\bullet })}^{0.5}$ with tidal acceleration ${a}_{t}={{GM}}_{\bullet }{\rm{\Delta }}R/{r}_{p}^{3}$ where ${\rm{\Delta }}R$ is the debris distance from the star center at the moment of breakup. Then, the energy of the disrupted debris is given by

$$\begin{eqnarray}&&{E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,{\rm{\Delta }}R)\approx \bar{e}{E}_{{\boldsymbol{m}}}({M}_{\bullet },m)-\displaystyle \frac{2{{kGM}}_{\bullet }{\rm{\Delta }}R}{{r}_{p}^{2}(\bar{e},{\ell },{M}_{\bullet },m)}\end{eqnarray} \tag{ 52 }$$

where ${\rm{\Delta }}R\in \{-{R}_{\star },{R}_{\star }\},$ the negative sign corresponds to the region toward the BH and k is the spin-up factor due to tidal torque by the SMBH on a star, given by (Alexander & Kumar 2001)

$$\begin{eqnarray*}k=\left\{\begin{array}{ll}1 & \text{non spin up (no change in angular velocity)}\\ 3 & \text{spin up to break up angular velocity.}\end{array}\right.\end{eqnarray*}$$

The maximum distance from the center of a star to the point where the debris is bound to the BH at the moment of disruption is obtained by setting E_d = 0 in Equation (52) and is given by

$$\begin{eqnarray}&&{R}_{l}(\bar{e},{\ell },{M}_{\bullet },m)=\displaystyle \frac{{r}_{p}^{2}(\bar{e},{\ell },{M}_{\bullet },m)\bar{e}}{2{{kr}}_{t}({M}_{\bullet },m)}.\end{eqnarray} \tag{ 53 }$$

Figure 14 shows the contour plot of ${x}_{l}\equiv {x}_{l}(\bar{e},{\ell },{M}_{\bullet },m)={R}_{l}(\bar{e},{\ell },{M}_{\bullet },m)/{R}_{\star }$ for ${M}_{6}=1$ and m = 1. The value of ${r}_{p}(\bar{e},{\ell },{M}_{\bullet },m)$ is less than Schwarzschild radius ${R}_{{\rm{s}}}({M}_{\bullet })$ for ${\ell }\leqslant 0.2$, whereas ${x}_{l}(\bar{e},{\ell },{M}_{\bullet },m)$ increases with $\bar{e}$ and the increase with ℓ is significant for high energy orbits. With the increase in the value of ${x}_{l}(\bar{e},{\ell },{M}_{\bullet },m),$ the mass of the debris bound to the BH increases and for ${x}_{l}(\bar{e},{\ell },{M}_{\bullet },m)=1,$ the entire debris is bound to the BH.

Zoom In Zoom Out Reset image size

The time taken for the most tightly bound debris to return its pericenter after disruption is given by

$$\begin{eqnarray}&&{t}_{m}(\bar{e},{\ell },{M}_{\bullet },m)=\displaystyle \frac{2\pi {{GM}}_{\bullet }}{{[2{E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,-{R}_{\star })]}^{3/2}}.\end{eqnarray} \tag{ 54 }$$

As the bound material falls back to its pericenter, it loses its energy and angular momentum, thus accreting into the SMBH and giving rise to the flare (Phinney 1989). The in-fall mass accretion rate at time t after disruption for the debris following Keplerian orbits is given by

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=\displaystyle \frac{{dM}}{{{dE}}_{d}}\displaystyle \frac{{{dE}}_{d}}{{da}}\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{1}{3}{(2\pi {{GM}}_{\bullet })}^{\displaystyle \frac{2}{3}}\displaystyle \frac{{dM}}{{{dE}}_{d}}{t}^{-\frac{5}{3}}\end{eqnarray} \tag{ 55 }$$

where a is the semimajor axis of the debris with orbital energy ${E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,{\rm{\Delta }}R).$ The term ${dM}/{{dE}}_{d}$ is the energy distribution within the bound matter and depends on the internal structure of the star (Phinney 1989; Lodato et al. 2009). We now write it in terms of dimensionless quantities using Equations (52) and (53), and modify the dimensionless quantities given in Lodato et al. (2009) by including the dependence on ${r}_{p}(\bar{e},{\ell },{M}_{\bullet },m)$ through ${R}_{l}(\bar{e},{\ell },{M}_{\bullet },m)$ and express

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{{x}_{l}-x}{{x}_{l}+1},x={x}_{l}-{\tau }^{-2/3}(1+{x}_{l})\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\mu }{d\varepsilon }=(1+{x}_{l})\displaystyle \frac{d\mu }{{dx}},\displaystyle \frac{d\mu }{d\tau }=\displaystyle \frac{2}{3}\displaystyle \frac{d\mu }{d\varepsilon }{\tau }^{\displaystyle \frac{-5}{3}}\end{eqnarray} \tag{ 57 }$$

where $\varepsilon (\bar{e},{\ell },{M}_{\bullet },m)$ = ${E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,{\rm{\Delta }}R)/$ ${E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,-{R}_{\star }),\;$ $x={\rm{\Delta }}R/{R}_{\star },\;$ $\mu =M/{M}_{\star }$ and $\tau (\bar{e},{\ell },{M}_{\bullet },m)$ $\;\equiv \;t/{t}_{m}(\bar{e},{\ell },{M}_{\bullet },m)$, where ${E}_{{dm}}={E}_{d}(\bar{e},{\ell },{M}_{\bullet },m,-{R}_{\star })$ is the energy of inner-most tightly bound debris. The term $d\mu$

$$\begin{eqnarray}&&\displaystyle \frac{d\mu }{{dx}}=\displaystyle \frac{3}{2}b{\displaystyle \int }_{x}^{1}{\theta }^{u}(x^{\prime} )x^{\prime} \;{dx}^{\prime} \end{eqnarray} \tag{ 58 }$$

where b is the ratio of central density ${\rho }_{c}$ to mean density $\bar{{\rho }_{\star }}=3{M}_{\star }/4\pi {R}_{\star }^{3}$ and θ is the solution of Lane–Emden equation for the given polytrope u related to the density by $\rho ={\rho }_{c}{\theta }^{u}.$ The total mass accretion rate is given by $\dot{M}(\bar{e},{\ell },{M}_{\bullet },m,t)=({M}_{\star }/{t}_{m})(d\mu /d\tau )$, which depends on the orbital parameters through x_l and t_m. We simulated the mass fallback rate for u = 1.5, which corresponds to ${\rm{\Gamma }}=5/3.$ Figure 15 shows the plot of $d\mu /d\tau$ for various values of x_l. With an increase in $\bar{e},$ the x_l increases and the orbital period of the debris decreases, which implies that the mass in-fall rate increases. Thus, the peak accretion rate increases with x_l. Next, we examine the conditions for formation of an accretion disk.

Zoom In Zoom Out Reset image size

The debris of the disrupted star follows a Keplerian orbit around the BH. This debris experiences stream–stream collision either due to an incoming stream that intersects with the outflowing stream at the pericenter (Kochanek 1994) or due to relativistic precession at the pericenter (Hayasaki et al. 2013). This stream–stream collision results in a shock breakout that circularizes the debris and forms an accretion disk (Ramirez-Ruiz & Rosswog 2009).

Even though the circularization timescale (time required for the debris to circularize to form an accretion disk) is not accurately known, it is roughly given by ${t}_{c}\approx {n}_{\mathrm{orb}}{t}_{m}$, where n_orb is the minimum number of orbits required for circularization (Ulmer 1999). As the debris falls toward the pericenter, it is accreted with an accretion rate given by Equation (55). The formation of an accretion disk depends on t_c and the accretion timescale t_a, which we define as the time required to consume 99% (at the 3σ level) of bound debris. If this timescale is less than t_c, the matter is accreted before the disk is formed. We approximated $d\mu /{dx}$ for convenience with a Gaussian function because it depends on the solution of the Lane–Emden equation, which is symmetric about the center of the star and is given by $d\mu /{dx}\simeq 1.192{e}^{-4.321{x}^{2}}.$ The total mass consumed in dimensionless time ${\tau }_{a}={t}_{a}/{t}_{m}$ is given by

$$\begin{eqnarray}&&{\rm{\Delta }}\mu ={\displaystyle \int }_{1}^{{\tau }_{a}}\displaystyle \frac{d\mu }{d\tau }\;d\tau \end{eqnarray} \tag{ 59 }$$

where $d\mu /d\tau$ is given by Equation (57). If f_r is the fraction of debris bound to the BH, then in time ${\tau }_{a},$ the mass accreted by the BH is ∼0.99f_r. Then, the accretion timescale t_a is given by

$$\begin{eqnarray}&&{t}_{a}(\bar{e},{\ell },{M}_{\bullet },m)\\ &&=\;{t}_{m}{\left(\displaystyle \frac{1+{x}_{l}}{{x}_{l}+\displaystyle \frac{1}{2.0787}{{\rm{Erf}}}^{-1}[0.997-1.962{f}_{r}]}\right)}^{\frac{3}{2}}\end{eqnarray} \tag{ 60 }$$

where

$$\begin{eqnarray}&&{f}_{r}\equiv {f}_{r}(\bar{e},{\ell },{M}_{\bullet },m)={\displaystyle \int }_{-1}^{{x}_{l}}\displaystyle \frac{d\mu }{{dx}}\;{dx}.\end{eqnarray} \tag{ 61 }$$

The orbital period of debris will vary due to energy. The energy gradient between the bound debris will fill out a ring and the initial spatial distance between the bound debris will determine the ring formation timescale t_r. Let the dispersion in the energy around the initial energy E be ${\rm{\Delta }}E.$ A ring is formed in the timescale ${t}_{r}=2\pi /{\rm{\Delta }}{\rm{\Omega }}$, where ${\rm{\Delta }}{\rm{\Omega }}$ is the dispersion in the orbital frequency (Hadrava et al. 2001). The orbital frequency ${\rm{\Omega }}={(2E)}^{3/2}/({{GM}}_{\bullet })$ and dispersion ${\rm{\Delta }}{\rm{\Omega }}\propto {E}^{1/2}{\rm{\Delta }}E$ and ${\rm{\Delta }}E$ are given by

$$\begin{eqnarray}{\rm{\Delta }}E(\bar{e},{\ell },{M}_{\bullet },m) & = & \displaystyle \frac{2{{kGM}}_{\bullet }}{{r}_{p}^{2}(\bar{e},{\ell },{M}_{\bullet },m)}\\ & & \times ({\rm{Min}}[{R}_{l}(\bar{e},{\ell },{M}_{\bullet },m),{R}_{\star }]+{R}_{\star }).\end{eqnarray} \tag{ 62 }$$

Then, t_r is given by

$$\begin{eqnarray}{t}_{r}(\bar{e},{\ell },{M}_{\bullet },m) & = & \displaystyle \frac{\pi }{3\sqrt{2}k}\displaystyle \frac{{r}_{p}^{2}(\bar{e},{\ell },{M}_{\bullet },m)}{{\bar{e}}^{1/2}{E}_{{\boldsymbol{m}}}^{1/2}({M}_{\bullet },m)}\\ & & \times {({R}_{\star }+{\rm{Min}}[{R}_{l}(\bar{e},{\ell },{M}_{\bullet },m),{R}_{\star }])}^{-1}.\end{eqnarray} \tag{ 63 }$$

The ratio ${{\mathcal{T}}}_{r}(\bar{e},{\ell },{M}_{\bullet },m)={t}_{r}(\bar{e},{\ell },{M}_{\bullet },m)/{t}_{a}(\bar{e},{\ell },{M}_{\bullet },m)$ is given by

$$\begin{eqnarray}&&{{\mathcal{T}}}_{r}(\bar{e},{\ell },{M}_{\bullet },m)=\displaystyle \frac{1}{3\sqrt{2}}\displaystyle \frac{{x}_{l}}{1+\mathrm{Min}[1,{x}_{l}]}\\ &&\qquad \times {\left(\displaystyle \frac{{x}_{l}+\displaystyle \frac{1}{2.0787}{{\rm{Erf}}}^{-1}[0.997-1.962{f}_{r}]}{{x}_{l}}\right)}^{\frac{3}{2}}\end{eqnarray} \tag{ 64 }$$

and an accretion disk is formed if ${{\mathcal{T}}}_{r}(\bar{e},{\ell },{M}_{\bullet },m)\lt 1.$ In Figure 16, the top panel (a) shows the contour plot of ${{\mathcal{T}}}_{r}(\bar{e},{\ell },{M}_{\bullet },m)$ for M₆ = 1, m = 1, and ${r}_{p}(\bar{e},{\ell },{M}_{6},m)\gt {R}_{{\rm{s}}},\;$ ${{\mathcal{T}}}_{r}\lt 1.$ The bottom panel (b) shows Max[${{\mathcal{T}}}_{r}({10}^{-6}\leqslant \bar{e}\leqslant 1,0\leqslant {\ell }\leqslant 1,1\leqslant {M}_{6}\leqslant 100)$] < 1, which implies that the bound debris will form an accretion disk.

Zoom In Zoom Out Reset image size

The radiation timescale of the disk is given by

$$\begin{eqnarray}&&{t}_{R}(\bar{e},{\ell },{M}_{\bullet },m)=\displaystyle \frac{{f}_{r}{M}_{\star }{c}^{2}}{\eta \dot{M}(\bar{e},{\ell },{M}_{\bullet },m)}{c}^{2}\end{eqnarray} \tag{ 65 }$$

where c is the light speed, $\dot{M}$ is the accretion rate, and η is the radiative efficiency of the disk. The viscous timescale t_v of the disk formed is given by

$$\begin{eqnarray}&&{t}_{v}={\displaystyle \int }_{{r}_{\mathrm{in}}}^{{r}_{c}}\displaystyle \frac{1}{{V}_{r}}\;{dr}\end{eqnarray} \tag{ 66 }$$

where V_r is the radial inflow velocity of matter in the disk, r_in is the inner radius of disk, and ${r}_{c}(\bar{e},{\ell },{M}_{\bullet },m)=2{r}_{p}(1-{r}_{p}\bar{e}/{r}_{t})$, where ${r}_{c}(\bar{e},{\ell },{M}_{\bullet },m)$ is the circularization radius. The radial inflow velocity is given by (Shakura & Sunyaev 1973; Strubbe & Quataert 2009)

$$\begin{eqnarray}&&{V}_{r}=-\displaystyle \frac{3}{2}\displaystyle \frac{\nu }{r}\displaystyle \frac{1}{f}\end{eqnarray} \tag{ 67 }$$

where $f=1-\sqrt{\frac{{r}_{\mathrm{in}}}{r}}$ and ν is the viscosity of the medium given by $\nu =\alpha {c}_{s}H.$ The parameter α is taken to be 0.1 and ${c}_{s}=H\sqrt{{{GM}}_{\bullet }/{r}^{3}}$, where H is the disk scale height. Strubbe & Quataert (2009) calculated the disk scale height for a slim disk (see Section 5.1) to be

$$\begin{eqnarray}&&\displaystyle \frac{H}{r}(\bar{e},{\ell },{M}_{\bullet },m,r)=\displaystyle \frac{3}{4}f\left(\displaystyle \frac{\dot{M}}{\eta {\dot{M}}_{E}}\right){\left(\displaystyle \frac{r}{{R}_{s}}\right)}^{-1}\\ &&\times {\left\{\displaystyle \frac{1}{2}+{\left\{\displaystyle \frac{1}{4}+\displaystyle \frac{3f}{2}{\left(\displaystyle \frac{\dot{M}}{\eta \dot{{M}_{E}}}\right)}^{2}{\left(\displaystyle \frac{r}{{R}_{s}}\right)}^{-2}\right\}}^{\displaystyle \frac{1}{2}}\right\}}^{-1}\end{eqnarray} \tag{ 68 }$$

where ${\dot{M}}_{E}$ is the Eddington mass accretion rate and $\dot{M}$ is taken to be the time averaged accretion rate. Using Equations (66) and (67), the viscous timescale is given by

$$\begin{eqnarray}{t}_{v}(\bar{e},{\ell },{M}_{\bullet },m) & = & \displaystyle \frac{2}{3\alpha }\displaystyle \frac{1}{\sqrt{{{GM}}_{\bullet }}}{\displaystyle \int }_{{r}_{\mathrm{in}}}^{{r}_{c}(\bar{e},{\ell },{M}_{\bullet },m)}\\ & & \times {\left(\displaystyle \frac{H(\bar{e},{\ell },{M}_{\bullet },m,r)}{r}\right)}^{-2}\sqrt{r}f\;{dr}.\end{eqnarray} \tag{ 69 }$$

The ratio ${{\mathcal{T}}}_{v}(\bar{e},{\ell },{M}_{\bullet },m)={t}_{v}(\bar{e},{\ell },{M}_{\bullet },m)/{t}_{R}(\bar{e},{\ell },{M}_{\bullet },m)$ is given by

$$\begin{eqnarray}{{\mathcal{T}}}_{v}(\bar{e},{\ell },{M}_{\bullet },m) & = & \displaystyle \frac{4}{9\alpha }\displaystyle \frac{1}{\sqrt{{{GM}}_{\bullet }}}\displaystyle \frac{{\eta }^{3}}{{f}_{r}{M}_{\star }}\displaystyle \frac{{\dot{M}}_{E}^{2}}{\dot{M}}{\displaystyle \int }_{{r}_{\mathrm{in}}}^{{r}_{c}}{\left(\displaystyle \frac{r}{{R}_{{\rm{s}}}}\right)}^{2}\displaystyle \frac{\sqrt{r}}{f}\\ & & \times {\left(\displaystyle \frac{1}{2}+{\left\{\displaystyle \frac{1}{4}+\displaystyle \frac{3f}{2}{\left(\displaystyle \frac{\dot{M}}{\eta {\dot{M}}_{E}}\right)}^{2}{\left(\displaystyle \frac{r}{{R}_{s}}\right)}^{-2}\right\}}^{\displaystyle \frac{1}{2}}\right)}^{2}\;{dr}\end{eqnarray} \tag{ 70 }$$

where radiative efficiency is typically $\eta =0.1$ and an accretion disk formed is a slim disk if ${{\mathcal{T}}}_{v}\lt 1.$ Figure 17 shows the contour plot of ${{\mathcal{T}}}_{v}(\bar{e},{\ell },{M}_{\bullet },m)$ for m = 1 and ℓ = 1 and 0.6. We conclude that the accretion disk formed is a slim disk for ${M}_{6}\leqslant 31.6.$ For higher mass SMBHs, a thin disk forms from the disrupted debris of a star on low energy orbit and ${\ell }\sim 1$, and a thick disk for a star on high energy orbit.

Zoom In Zoom Out Reset image size

The Eddington mass accretion rate is given by

$$\begin{eqnarray}&&{\dot{M}}_{E}=\displaystyle \frac{4\pi {{GM}}_{\bullet }}{\eta \kappa c}\end{eqnarray} \tag{ 71 }$$

where κ is the opacity of the medium taken to be Thompson opacity and η is the radiative efficiency. For a given tidally disrupted star, the accretion disk formed has a super Eddington phase if ${M}_{c}(\bar{e},{\ell },m)\gt {M}_{\bullet }$ where ${M}_{c}(\bar{e},{\ell },m)$ is the critical BH mass. We have numerically equated the peak accretion rate ${\dot{M}}_{p}$ and ${\dot{M}}_{E}$, and obtained ${M}_{c}(\bar{e},{\ell },m)$ as shown in Figure 18 for m = 1. The ${M}_{c}(\bar{e},{\ell },m)$ decreases with ℓ and increases with $\bar{e}.$ For a given ${\ell },$ an increase in $\bar{e}$ increases E_d and thus the orbital period of the disrupted debris decreases, which results in an increase in the $\dot{M}$ and hence the peak accretion rate ${\dot{M}}_{p}.$ For a given $\bar{e},$ the pericenter r_p increases with ℓ, which results in a decrease in E_d and thus ${\dot{M}}_{p}$ decreases. For a given $\bar{e}$ and ${\ell },\;$ ${r}_{p}(\bar{e},{\ell },m)$ increases with m, which results in a decrease in E_d. The decrease in E_d implies an increase in fallback time and a decrease in peak accretion rate ${\dot{M}}_{p}$, which results in a decrease in ${M}_{c}(\bar{e},{\ell },m).$

Zoom In Zoom Out Reset image size

5.1. Super Eddington Phase

For ${M}_{\bullet }\leqslant {M}_{c}(\bar{e},{\ell },m),$ the radiation produced by viscous stress in the rotating disk is trapped by electron scattering and the disk is radiatively inefficient. The time for the photon to diffuse out of the gas is longer than both the inflow time in the disk and the dynamical time for the outflow. Thus the disk is radiation pressure dominated and the opacity is given by the electron scattering. We assume that the opacity due to the Thompson scattering and in our flux calculation adopt the work of Strubbe & Quataert (2009). The strong radiative pressure induces an outflowing wind. At the launch radius ${r}_{L}(\bar{e},{\ell },{M}_{\bullet },m)={r}_{c}(\bar{e},{\ell },{M}_{\bullet },m)=2{r}_{p}(1-{r}_{p}\bar{e}/{r}_{t}),$ the internal energy is converted into the kinetic energy of the outflows, and the material leaves the disk with the temperature at the launch radius as determined by ${{aT}}_{L}^{4}\approx (1/2)\rho ({r}_{L}){v}_{w}^{2}$ (Lodato & Rossi 2011), where a is the radiation constant. The outflow geometry is assumed to be spherical. The photons are trapped up to the radius where $\rho \kappa r\sim 1$ and the radius of photosphere ${r}_{\mathrm{ph}}\sim {(\rho \kappa )}^{-1}$ is given by

$$\begin{eqnarray}&&{r}_{\mathrm{ph}}(\bar{e},{\ell },{M}_{\bullet },m,t)=\displaystyle \frac{{f}_{\mathrm{out}}\dot{M}\kappa }{4\pi {v}_{w}},\end{eqnarray} \tag{ 72 }$$

where ${M}_{\bullet }={M}_{6}{10}^{6}{M}_{\odot },\;$ ${f}_{\mathrm{out}}=({\dot{M}}_{\mathrm{out}}/\dot{M})$, and ${v}_{w}={f}_{v}\sqrt{{{GM}}_{\bullet }/{r}_{L}}$ is the velocity of outflowing wind where f_v is taken to be unity. The outflowing wind is assumed to expand adiabatically so that the density is $\rho (r)\propto {T}^{3}(r)$, where T is the temperature of the outflowing wind. Using this scaling relation, the temperature at the photosphere is ${T}_{\mathrm{ph}}={T}_{L}{(\rho ({r}_{\mathrm{ph}})/\rho ({r}_{L}))}^{1/3}$, and using Equation (72), T_ph is given by

$$\begin{eqnarray}&&{T}_{\mathrm{ph}}(\bar{e},{\ell },{M}_{\bullet },m,t)={(4\pi )}^{\frac{5}{12}}{\left(\displaystyle \frac{1}{2a}\right)}^{\frac{1}{4}}\\ &&\times {\kappa }^{-\frac{2}{3}}{f}_{\mathrm{out}}^{-\frac{5}{12}}{f}_{v}^{\frac{11}{12}}{\dot{M}}^{-\frac{5}{12}}{r}_{L}^{-\frac{7}{24}}{({{GM}}_{\bullet })}^{\frac{11}{24}}.\end{eqnarray} \tag{ 73 }$$

Dotan & Shaviv (2011) have calculated the fraction of outflowing material from the super Eddington slim disk with $\dot{M}/{\dot{M}}_{E}$ = 1, 5, 10, and 20, respectively. We approximated their result by the following relation (Lodato & Rossi 2011)

$$\begin{eqnarray}&&{f}_{\mathrm{out}}(\bar{e},{\ell },{M}_{\bullet },m,t)=\displaystyle \frac{2}{\pi }\mathrm{arctan}\left[\displaystyle \frac{1}{4.5}\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{E}}-1\right)\right].\end{eqnarray} \tag{ 74 }$$

Thus, the luminosity from the outflowing wind is given by ${L}_{\nu }^{\mathrm{out}}(\bar{e},{\ell },{M}_{\bullet },m,t)=4\pi {r}_{\mathrm{ph}}^{2}{B}_{\nu }({T}_{\mathrm{ph}})$ and ${B}_{\nu }({T}_{\mathrm{ph}})$ is the intensity obtained assuming the outflowing wind as a black body.

In the super Eddington phase, the time for the photon to diffuse out of the disk is longer than the viscous time, so that the disk that is formed is thick and advective, whereas in the case when $\dot{M}\leqslant {\dot{M}}_{E},$ the disk is thin and cools by radiative diffusion. Strubbe & Quataert (2009) considered a slim disk model by introducing an additional advection term in the energy conservation equation, where the effective temperature profile of the disk as a function of radius is given by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{SB}}{T}_{e}^{4}(\bar{e},{\ell },{M}_{\bullet },m,r,t)=8.54\times {10}^{17}\\ &&\times \displaystyle \frac{{M}_{6}^{-1}{\left(\displaystyle \frac{\eta }{0.1}\right)}^{-1}{\left(\displaystyle \frac{r}{{R}_{{\rm{s}}}}\right)}^{-3}\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{E}}\right)f}{\displaystyle \frac{1}{2}+{\left\{\displaystyle \frac{1}{4}+\displaystyle \frac{3f}{2}{\left(\displaystyle \frac{\dot{M}}{\eta {\dot{M}}_{E}}\right)}^{2}{\left(\displaystyle \frac{r}{{R}_{s}}\right)}^{-2}\right\}}^{\displaystyle \frac{1}{2}}}\;{\rm{W}}\;{{\rm{m}}}^{-2}\end{eqnarray} \tag{ 75 }$$

where $f=1-\sqrt{\frac{{r}_{\mathrm{in}}}{r}},\;$ ${\sigma }_{\mathrm{SB}}$ is Stefan–Boltzmann constant, r_in is the inner radius of the disk, and R_s is the Schwarzschild radius BH.

5.2. Sub-Eddington Phase

The disk is sub-Eddington for the BH mass ${M}_{\bullet }\gt {M}_{c}(\bar{e},{\ell },m)$ and $\dot{M}\lt {\dot{M}}_{E}.$ We then consider the disk as the radiative thin disk whose the temperature profile is given by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{SB}}{T}_{e}^{4}(\bar{e},{\ell },{M}_{\bullet },m,r,t)=8.54\times {10}^{17}{M}_{6}^{-1}\\ &&\times {\left(\displaystyle \frac{\eta }{0.1}\right)}^{-1}{\left(\displaystyle \frac{r}{{R}_{{\rm{s}}}}\right)}^{-3}\left(\displaystyle \frac{\dot{M}}{{\dot{M}}_{E}}\right)f\;{\rm{W}}\;{{\rm{m}}}^{-2}.\end{eqnarray} \tag{ 76 }$$

We see that Equation (75) is the modified temperature profile of the thin disk and follows the thin disk for $\dot{M}\lt {\dot{M}}_{E}.$ Since the sub-Eddington phase exists only for a certain duration, we assume that Equation (75) is the temperature profile for the entire duration as it approaches the thin disk model for $\dot{M}\lt {\dot{M}}_{E}.$ Assuming the disk to be a black body, the intensity of the disk is given by ${B}_{\nu }({T}_{e}(\bar{e},{\ell },{M}_{\bullet },m,r,t))$ and thus the disk luminosity is given by

$$\begin{eqnarray}{L}_{\nu }^{\mathrm{Disk}}(\bar{e},{\ell },{M}_{\bullet },m,t) & = & {\displaystyle \int }_{{r}_{\mathrm{in}}}^{{r}_{c}(\bar{e},{\ell },{M}_{\bullet },m)}\\ & & \times {B}_{\nu }({T}_{e}(\bar{e},{\ell },{M}_{\bullet },m,r,t))2\pi r\;{dr}.\end{eqnarray} \tag{ 77 }$$

The total luminosity can be written as

$$\begin{eqnarray*}&&{L}_{\nu }(\bar{e},{\ell },{M}_{\bullet },m,t)\\ &&=\;\left\{\begin{array}{ll}\begin{array}{l}{L}_{\nu }^{\mathrm{Disk}}(\bar{e},{\ell },{M}_{\bullet },m,t)\\ +\;{L}_{\nu }^{\mathrm{out}}(\bar{e},{\ell },{M}_{\bullet },m,t)\end{array} & {M}_{\bullet }\lt {M}_{c}(\bar{e},{\ell },m)\\ {L}_{\nu }^{\mathrm{Disk}}(\bar{e},{\ell },{M}_{\bullet },m,t) & {M}_{\bullet }\geqslant {M}_{c}(\bar{e},{\ell },m).\end{array}\right.\end{eqnarray*}$$

If ${\nu }_{l}$ and ${\nu }_{h}$ are the minimum and maximum frequency of the spectral band, then the luminosity of the emitted radiation in the given spectral band in the rest frame of the galaxy is given by

$$\begin{eqnarray}&&{L}_{e}(\bar{e},{\ell },{M}_{\bullet },m,z,t)={\displaystyle \int }_{{\nu }_{l}(1+z)}^{{\nu }_{h}(1+z)}{L}_{\nu }(\bar{e},{\ell },{M}_{\bullet },m,t)\;d\nu .\end{eqnarray} \tag{ 78 }$$

The observed flux ${f}_{\mathrm{obs}}(\bar{e},{\ell },{M}_{\bullet },m,z,t)$ = ${L}_{e}(\bar{e},{\ell },{M}_{\bullet },m,z,t)/(4\pi {{\rm{d}}}_{L}^{2}(z))$, where z is the redshift and d_L is the luminosity distance, and the radiation is observed only if

$$\begin{eqnarray}&&{f}_{l}\lt {f}_{\mathrm{obs}}(\bar{e},{\ell },{M}_{\bullet },m,z,t)=\displaystyle \frac{{L}_{e}(\bar{e},{\ell },{M}_{\bullet },m,z,t)}{4\pi {d}_{L}^{2}(z)}\end{eqnarray} \tag{ 79 }$$

where f_l is the sensitivity of the detector. The Equation (79) is utilized to generate a digital signal A(t) such that

$$\begin{eqnarray}A(t)=\left\{\begin{array}{ll}1 & \mathrm{if}\;\mathrm{eqn}\;(79)\;\mathrm{holds}\;\mathrm{true}\\ 0 & \mathrm{if}\;\mathrm{eqn}\;(79)\;\mathrm{does}\;\mathrm{not}\;\mathrm{hold}\;\mathrm{true}.\end{array}\right.\end{eqnarray} \tag{ 80 }$$

The width of the digital signal gives the duration of the flare detection used in the event rate calculation (see Section 6.3).

As an example, the observed flux in optical g band is shown in Figure 19 for ${M}_{6}=1$, z = 0.1, and ℓ = 0.6 (blue), 0.8 (red), 1.0 (orange). For ${M}_{6}=1,$ both the outflowing wind and disk contribute to the observed flux and the flux from the wind dominates in the initial time. We can observe a dip in flux due to the outflowing wind whose ${r}_{\mathrm{ph}}\propto \dot{M}$ and ${T}_{\mathrm{ph}}\propto {\dot{M}}^{-5/12}$ and the occurrence time of dip is nearly at the time of peak accretion rate ${\dot{M}}_{p}$ (see Figure 15). With an increase in $\bar{e},$ the ${\dot{M}}_{p}$ increases, which results in an increase in r_ph and a decrease in T_ph, and thus a decrease in the intensity of radiation $B({T}_{\mathrm{ph}}).$ The flux is due to outflowing wind $\propto {r}_{\mathrm{ph}}^{2}B({T}_{\mathrm{ph}})$ and decreases with ${\dot{M}}_{p}$ if the decline in $B({T}_{\mathrm{ph}})$ is higher than the rise in r_ph. Ulmer (1999) predicted that the minimum value of n_orb for the disrupted debris to get circularized is ∼2–3, and thus we have utilized the f_obs starting from the time $\tau =3$ to generate a digitized signal and calculate the duration of the detection.

Zoom In Zoom Out Reset image size

For any transient survey, the net detectable TDE rate depends on the number density of non-active galaxies, the theoretical capture rate per galaxy (see Section 2), the luminosity distance of galaxies, the sensitivity of the detector, and the duration of flare detection. In this section, we will carry out the detailed calculation of each quantity separately and then combine them in Section 6.3.

6.1. Number Density of Non Active Galaxies

The number density of quasars is a function of redshift and luminosity, where the quasars emitting radiation of low intensity ($L\lt {10}^{40}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$) are non-active galaxies as compared to the normal quasars ($L\approx {10}^{45}-{10}^{46}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$). According to the Soltan (1982) argument, if quasars were powered by accretion onto an SMBH, then such SMBH must exist in our local universe as "dead" quasars or non-active galaxies. The number density of galaxies (quasars) can be obtained using the quasar luminosity function (QLF; Hopkins et al. 2007),

$$\begin{eqnarray}&&\displaystyle \frac{d\psi }{d\mathrm{log}L}=\displaystyle \frac{{\psi }_{*}}{{(\displaystyle \frac{L}{{L}_{*}})}^{{\gamma }_{1}}+{(\displaystyle \frac{L}{{L}_{*}})}^{{\gamma }_{2}}},\end{eqnarray} \tag{ 81 }$$

where ${\psi }_{*},{L}_{*},{\gamma }_{1},{\gamma }_{2}$ are a function of redshift z. As TDEs are the main observational signatures of quiescent galaxies, we need to determine the number density of quiescent galaxies at any redshift z. Chen et al. (2007) used the QLF to obtain the duty cycle $\delta (z)={10}^{-3}{\left(z/0.1\right)}^{2.5}$, where $\delta (z)$ is defined as the ratio of the number of active galaxies to the total number of galaxies. Thus, the BH mass function of quiescent galaxies is given by

$$\begin{eqnarray}&&\displaystyle \frac{d\psi }{{{dM}}_{\bullet }}=(1-\delta (z))\displaystyle \frac{d\psi }{d\mathrm{log}L}\left|\displaystyle \frac{d\mathrm{log}L}{{{dM}}_{\bullet }}\right|\end{eqnarray} \tag{ 82 }$$

where L is the luminosity of the quasars, which is taken to be $L=\eta {L}_{E}=\eta 4\pi {{GM}}_{\bullet }c/\kappa$, where L_E is the Eddington luminosity, η is 0.1, and κ is the opacity due to Thompson scattering. This gives the number density of quiescent galaxies as a function of BH mass M_• and redshift z, as

$$\begin{eqnarray}&&\displaystyle \frac{d\psi }{{{dM}}_{\bullet }}({M}_{\bullet },z)=(1-\delta (z))\displaystyle \frac{1}{{M}_{\bullet }}\displaystyle \frac{d\psi }{d\mathrm{log}L}\end{eqnarray} \tag{ 83 }$$

6.2. Luminosity Distance

We assume ΛCDM cosmology with ${{\rm{\Omega }}}_{m}=0.315,{{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.685$, ${H}_{o}=67.3\;\mathrm{Km}\;{{\rm{s}}}^{-1}\;{\mathrm{Mpc}}^{-1}$ (Planck Collaboration et al. 2013). The luminosity distance as a function of redshift z is given as

$$\begin{eqnarray}&&{{\rm{d}}}_{L}(z)=(1+z)\displaystyle \frac{c}{{H}_{o}}{\displaystyle \int }_{0}^{z}\displaystyle \frac{1}{{({(1+z^{\prime} )}^{3}{{\rm{\Omega }}}_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}})}^{0.5}}\;{dz}^{\prime} .\end{eqnarray} \tag{ 84 }$$

Consider now a small volume of the universe at redshift z with radial width ${dz}$ covering an opening angle ω on the observer's sky (Khabibullin et al. 2014). The comoving volume of the slice is

$$\begin{eqnarray}&&{{dV}}_{c}=\omega {d}_{H}^{3}\displaystyle \frac{{I}^{2}(z)}{W(z)}{dz}\end{eqnarray} \tag{ 85 }$$

where $\omega =4\pi {f}_{s},{d}_{H}=c/{H}_{o},W(z)={({(1+z)}^{3}{{\rm{\Omega }}}_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}})}^{0.5},$

$$\begin{eqnarray}&&I(z)={\displaystyle \int }_{0}^{z}\displaystyle \frac{1}{{({(1+z^{\prime} )}^{3}{{\rm{\Omega }}}_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}})}^{0.5}}\;{dz}^{\prime} \end{eqnarray} \tag{ 86 }$$

and f_s is the fraction of sky observed.

6.3. Probability of Flare Detection

We generate the spectrum in the form of digital signal using Equation (80), and the width of the digital signal provides the duration of flare detection $\delta {t}_{f}(\bar{e},{\ell },{M}_{\bullet },m,z).$ If ${t}_{\mathrm{cad}}\;{\rm{and}}\;{t}_{\mathrm{int}}$ are the cadence and integration time of the detector, then the probability of detection of an event is given by

$$\begin{eqnarray}&&P(\bar{e},{\ell },{M}_{\bullet },m,z)={\rm{Min}}\left[1,\displaystyle \frac{\delta {t}_{f}(\bar{e},{\ell },{M}_{\bullet },m,z)}{{t}_{\mathrm{cad}}+{t}_{\mathrm{int}}}\right]\end{eqnarray} \tag{ 87 }$$

Using Equations (28), (48), (83), (85), and (87), the net detectable event rate by the detector is given by

$$\begin{eqnarray}\displaystyle \frac{{d}^{5}{\dot{N}}_{D}(\gamma ,\bar{e},{\ell },{M}_{\bullet },m,z)}{{{dM}}_{\bullet }{dmd}\bar{e}d{\ell }{dz}} & = & \omega {d}_{H}^{3}\left(\displaystyle \frac{d{\rm{\Psi }}}{{{dM}}_{\bullet }}\right)\displaystyle \frac{{d}^{2}{\dot{N}}_{t}}{d\bar{e}d{\ell }{dm}}\\ & & \times (\gamma ,\bar{e},{\ell },{M}_{\bullet },m)\displaystyle \frac{{I}^{2}(z)}{W(z)}\\ & & \times P(\bar{e},{\ell },{M}_{\bullet },m,z).\end{eqnarray} \tag{ 88 }$$

With the given initial parameters ${M}_{\bullet },\;m,\;\bar{e},\;{\ell }$, and z, we generate the light curves using Equation (78) in the optical g and soft X-ray bands. The generated spectrum is compared with the sensitivity of the detector f_l to generate a digital signal using Equation (80), and the width of the digital signal gives the duration of flare detection. The range of initial parameters in the calculation are taken to be ${M}_{6}={M}_{\bullet }/{10}^{6}{M}_{\odot }$ = $1-100,\;m=0.8-150,\;\bar{e}$ = ${\bar{e}}_{h}-1,\;{\ell }=0-1$, and $z=0-{z}_{s}$, where ${z}_{s}({M}_{\bullet },\;m,\;\bar{e},\;{\ell })$ is the detection limit of the survey. Then, using Equation (88), we calculated the net detectable rate by integrating in steps over redshift $z,\;{\ell },\;\bar{e},\;m$, and finally over M_•, such that

$$\begin{eqnarray}{\dot{N}}_{D} & = & {\displaystyle \int }_{1}^{100}{{dM}}_{6}{\displaystyle \int }_{0.8}^{150}{dm}{\displaystyle \int }_{{\bar{e}}_{h}}^{1}d\bar{e}{\displaystyle \int }_{0}^{1}d{\ell }{\displaystyle \int }_{0}^{{z}_{s}}\\ & & \times {dz}\displaystyle \frac{{d}^{5}{\dot{N}}_{D}(\gamma ,\bar{e},{\ell },{M}_{6},m,z)}{{{dM}}_{6}{dmd}\bar{e}d{\ell }{dz}}.\end{eqnarray} \tag{ 89 }$$

The detectable rate per M₆ integrated over z, $\bar{e},$ℓ, and m for various γ is shown in Figure 20 for LSST and Pan-STARRS 3π detectors parameters. As the ${\dot{N}}_{t}$ and BH mass functions decrease with M₆, the $d{\dot{N}}_{D}/{{dM}}_{6}$ decreases with M₆. The net ${\dot{N}}_{D}$ integrated over $\bar{e},\;$ ${\ell }$, m, and M₆ is plotted as a function of γ in Figure 21 for various missions. The resulting ${\dot{N}}_{D}\propto {\gamma }^{{\rm{\Delta }}}$ where Δ is given in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Mission Instrument Parameters and Predicted Rate of the Surveys

Survey	Band	fs	Sensitivity/Flux	Cadence Time	Integration Time	${\dot{N}}_{D}$ (${\mathrm{Yr}}^{-1}$)a	${\gamma }_{s}$ b	${\dot{N}}_{\mathrm{obs}}$ (${\mathrm{Yr}}^{-1}$)c	ϒd	Δe
				(s)	(s)	$\gamma =0.7\;\pm \;0.1$
LSST	Optical	0.5	24.5 AB mag (g band)	2.6 × 105	10	5003 ± 1421	0.63	4131	0.91	1.97
Pan-STARRS (MDS)	Optical	0.0012	24.8 AB mag (g band)	3.46 × 105	30	12.3 ± 3.5	0.77	15	0.92	1.98
Pan-STARRS $3\pi$	Optical	0.75	24 AB mag (g band)	6.05 × 105	30	6337 ± 1800	0.48	3106	0.85	1.94
eROSITA	X-ray	1	$2.4\times {10}^{-14}(\;\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{cm}}^{-2})$	1.58 × 107	$1.6\times {10}^{3}$	679.5 ± 195	...	...	0.7	2.06

Notes. The parameters of the survey are taken from (1) LSST (Strubbe & Quataert 2009 and http://www.lsst.org/lsst/overview/), (2) Pan-STARRS (MDS; Medium Deep survey; van Velzen et al. 2011), (3) Pan-STARRS 3π (Strubbe & Quataert 2009 and http://pan-starrs.ifa.hawaii.edu/public/), (4) eROSITA (SRG; Khabibullin et al. 2014 and http://www.mpe.mpg.de/eROSITA)

^aOur predicted values along with the error estimates for an assumed range of ${\rm{\Delta }}\gamma =0.1$ around a typically observed median of $\gamma =0.7.$ ^b ${\dot{N}}_{D}({\gamma }_{s})={\dot{N}}_{D}$ estimated by van Velzen et al. (2011). ^cResults from van Velzen et al. (2011). ^dDetection efficiency of the detector given by Equation (91). ^eDetectable rate ${\dot{N}}_{D}\propto {\gamma }^{{\rm{\Delta }}}$.

Download table as:  ASCIITypeset image

Using Equations (28), (36), (83), and (85), the occurrence rate of TDE is given by

$$\begin{eqnarray}{\dot{N}}_{o} & = & {\displaystyle \int }_{1}^{100}{{dM}}_{6}{\displaystyle \int }_{0.8}^{150}{dm}{\displaystyle \int }_{{\bar{e}}_{h}}^{1}d\bar{e}{\displaystyle \int }_{0}^{1}d{\ell }{\displaystyle \int }_{0}^{{z}_{s}}{dz}\omega {d}_{H}^{3}\\ & & \times \left(\displaystyle \frac{d{\rm{\Psi }}}{{{dM}}_{6}}\right)\displaystyle \frac{{d}^{3}{\dot{N}}_{t}}{d\bar{e}d{\ell }{dm}}(\gamma ,\bar{e},{\ell },{M}_{6},m)\displaystyle \frac{{I}^{2}(z)}{W(z)}\end{eqnarray} \tag{ 90 }$$

where integration limits are same as those taken for Equation (89). We define the detection efficiency of TDE for a detector to be

$$\begin{eqnarray}&&{\rm{\Upsilon }}=\displaystyle \frac{{\dot{N}}_{D}}{{\dot{N}}_{o}}.\end{eqnarray} \tag{ 91 }$$

The ${\dot{N}}_{D}$ calculated for LSST, Pan-STARRS 3π, Pan-STARRS MDS, and eROSITA mission along with their detection efficiency are given in Table 2. van Velzen et al. (2011) estimated the event rates on the basis of observational studies in the optical bands for the Sloan Digital Sky survey (SDSS) and scaled the result of SDSS to the other missions using the relation $\dot{N}\propto {f}_{s}{f}_{l}^{-3/2}$ (Gezari et al. 2009). This relation is valid only if we assume all other parameters, such as cadence and integration time of the detector, to be same and Table 2 shows the estimated rates by van Velzen et al. (2011) ${\dot{N}}_{\mathrm{obs}}$ (column with c notation) and our predicted rates ${\dot{N}}_{D}$ (column with a notation). The value of γ from the observed density profile is in the range ∼ 0.5–1.2 (Wang & Merritt 2004; Stone & Metzger 2014). Our results are in reasonable agreement with their results. Khabibullin et al. (2014) calculated the number of events N detectable at any moment in the X-ray band for the eROSITA mission assuming a constant theoretical rate for ${M}_{6}\sim 1-10$ and light curve profile to follow the ${t}^{-5/3}$ law, whereas we have followed a more rigorous calculation to obtain ${\dot{N}}_{D}.$ Our prediction for eROSITA does not include the limitations in Khabibullin et al. (2014), but are more precise and in agreement with their rough estimate.

The values of γ for which our predictions of ${\dot{N}}_{D}$ match with the scaled-up values in van Velzen et al. (2011) are shown as ${\gamma }_{s}$ in Table 2. The only free parameter in our estimate is γ and this is likely to vary from source to source. Not knowing the expected distribution of γ as a function of redshift, we calculated the error in our estimation of ${\dot{N}}_{D}$ by taking a fiduciary range in the observed median of $\gamma =0.7\pm 0.1$, as is shown in Table 2.

The star's initial orbital parameters E and J have significant effects on both stellar and accretion dynamics. We have seen that the effect of J, which has not been included previously, plays a crucial role in constructing the shape of light curve profiles.

We employed a single power density model because most of the galaxies given in Wang & Merritt (2004) and Stone & Metzger (2014) have a break radius ${r}_{b}\gt {r}_{h}$ and we calculated the ${\dot{N}}_{t}\sim {10}^{-5}-{10}^{-4}$ Yr⁻¹, which shows a nonlinear dependence with M_• as shown by Wang & Merritt (2004) for single stellar mass distribution with a Nuker profile. Using Equations (28), (49), (83), and (85) the galaxy average ${\dot{N}}_{t}$ is given by

$$\begin{eqnarray}&&\left\langle{\dot{N}}_{t}\right\rangle(\gamma )=\displaystyle \frac{{\displaystyle \int }_{{10}^{6}{M}_{\odot }}^{{10}^{8}{M}_{\odot }}{{dM}}_{\bullet }{\displaystyle \int }_{0}^{1}{dz}\left(\displaystyle \frac{d\psi }{{{dM}}_{\bullet }}\right)\displaystyle \frac{{I}^{2}(z)}{W(z)}{\dot{N}}_{t}(\gamma ,{M}_{\bullet })}{{\displaystyle \int }_{{10}^{6}{M}_{\odot }}^{{10}^{8}{M}_{\odot }}{{dM}}_{\bullet }{\displaystyle \int }_{0}^{1}{dz}\left(\displaystyle \frac{d\psi }{{{dM}}_{\bullet }}\right)\displaystyle \frac{{I}^{2}(z)}{W(z)}}\end{eqnarray} \tag{ 92 }$$

and is shown in Figure 22. For $\gamma =0.7,\;$ $\left\langle{\dot{N}}_{t}\right\rangle\sim 6.8\times {10}^{-5}$ yr⁻¹, which is close to the observational inferred value $\sim {10}^{-5}$ yr⁻¹ and $\left\langle{\dot{N}}_{t}\right\rangle\approx {\gamma }^{2}$ for $\gamma \leqslant 1.2.$ The sample of galaxies taken by Stone & Metzger (2014) have γ varying over all ranges up to $\gamma =1.2$, which implies that their $\langle {\dot{N}}_{t}\rangle$ is γ independent, whereas we have calculated $\langle {\dot{N}}_{t}\rangle$ assuming that all galaxies have the same γ. The discrepancy in theory and observation is smaller in our model for $\gamma \leqslant 1$ compared with Stone & Metzger (2014), who have predicted $\langle {\dot{N}}_{t}\rangle \sim {\rm{few}}\times {10}^{-4}$ yr⁻¹ by taking into account the Schechter BH mass function and a Nuker profile. The γ averaged $\left\langle{\dot{N}}_{t}\right\rangle$ integrated over the range $0.6\leqslant \gamma \leqslant 1.2$ is $\sim 1.3\times {10}^{-4}$ yr⁻¹.

Zoom In Zoom Out Reset image size

Our calculation of ${\dot{N}}_{t}$ is based on a simple model of two body relaxations in a spherical potential. The alternative relaxation mechanism such as resonant relaxation (Rauch & Tremaine 1996), which dominates for high energy orbits and anomalous relaxation for highly eccentric orbits $e\to 1$ (Hamers et al. 2014) and non-spherical potential can moderately change the values of ${\dot{N}}_{t}.\;$Merritt (2015b), taking into account the resonant relaxation theory in a potential dominated by the BH, has shown that the enhancement of angular momentum diffusion at large binding energies results in a depletion of DF at those energies, and results in a density deficit core, but the effect of resonant relaxation on ${\dot{N}}_{t}$ is moderate for a power law model, such as the one used here; however, resonant relaxation for more general DF can be important. The assumption of isotropic velocity distribution overestimates the TDE rates if the true velocity distribution is anisotropic and tangentially biased, and underestimates the rates if the aniostropy is radially biased. This is because a tangentially biased distribution will have longer angular momentum relaxation time compared with a radially biased distribution. Recent N-body simulations by Zhong et al. (2014) have indicated that the presence of a loss cone will bias the orbits toward tangential anisotropy at a smaller distance and radial anisotropy at a larger distance, and that the tangential anisotropy is minor near ${{\mathcal{E}}}_{c}.$ From both observational and theoretical perspectives, it is unclear whether the galactic nuclei are sufficiently anisotropic (and overwhelmingly in the tangential direction) so as to reduce TDE rates to the observationally inferred values.

Rees (1988) and others have considered the stellar orbit to be nearly parabolic. We included the angular momentum J in the calculation and studied the effect of J on accretion dynamics. We modified the dimensionless quantities given in Lodato et al. (2009) and for the low eccentric orbits, which results in an increase in peak accretion rate. Strubbe & Quataert (2009) and Lodato & Rossi (2011) calculated the spectral profile for a parabolic orbit that does not have any dip in their luminosity profile. The inclusion of J induces a dip in the light curve profile, which gets deeper with increased energy. We can also see that our results in the optical band match the result of Lodato & Rossi (2011) for $\bar{e}\ll 1.$

In general, the accretion of matter into the BH is non-steady because the mass at the outer radii are higher than the mass at inner radii. Montesinos Armijo & de Freitas Pacheco (2011) evaluated the surface density and temperature profile assuming the accretion disk to be thin and the accretion rate $\propto {t}^{-5/3}.$ A model for a non-steady accretion mechanism that includes both the super- and sub-Eddington phase is required to better understand the evolution and emission from the disk. The α viscosity prescription used by Strubbe & Quataert (2009) is not applicable in the super-Eddington phase due to low efficiency and high opacity of the disk, so a general viscosity prescription, such as $\nu \propto {{\rm{\Sigma }}}^{d}(r){r}^{e}$ where d and e are constants, can be used to evaluate the accretion disk, where ${\rm{\Sigma }}(r)$ is surface density profile (Mangalam 2001; Shen & Matzner 2014).

We have built simple analytical expressions to evaluate the condition for the formation of an accretion disk. We did not include the stream interactions and relativistic effects in the calculations, which can speed up the rate of formation of the disk through the angular momentum exchange. The hydrodynamical simulations by Ramirez-Ruiz & Rosswog (2009) have shown that the debris interactions result in the formation of an accretion disk with mass accretion rate showing deviation from Lodato et al. (2009) in early time and following ${t}^{-5/3}$ in the late stage. Very recently, Bonnerot et al. (2015) performed hydrodynamical simulations for a star on a highly elliptical orbit with the resulting debris undergoing apsidal precession; they found that the higher the eccentricity and/or the deeper the encounter, the faster the circularization. For an efficient cooling, the debris forms a thin and narrow ring of gas. For an inefficient cooling, they settle in a thick and extended torus, mostly centrifugally supported against gravity. The general relativistic hydrodynamical simulations by Shiokawa et al. (2015) have shown that the accretion rate still rises sharply and then decays as a power law; however, its maximum is 10% smaller than the previous expectation, and timescale of the peak accretion is longer than previously predicted values. This is due to the mass accumulation at higher radius because of angular momentum exchange at large radii. The overall conclusion is that the resulting debris will form an accretion disk. The thickness of the disk formed and the circularization timescale as a function of stellar parameters and M_• still need to be evaluated. We think that while it is important to include relativistic effects and stream collisions, it is unlikely that it will change the conclusion, and will only slightly change the conditions derived for the disk formation.

Strubbe & Quataert (2009) predicted ${\dot{N}}_{D}$ assuming a constant capture rate, stellar orbits to be parabolic, and the flare duration to be the duration of Eddington phase, which is obtained assuming ${dM}/{{dE}}_{d}$ as a constant. Gezari et al. (2008) and van Velzen et al. (2011) used the observed detectable rate for the GALEX mission in the Near Ultra Violet and SDSS in optical band, respectively, and scaled it to the other missions assuming survey parameters such as cadence and integration time to be same. van Velzen et al. (2011) observationally estimated higher rates compared with the estimation by Gezari et al. (2008) due to low sample size. We performed a detailed calculation, taking in account both stellar and accretion dynamics, and predicted the detectable rates that are in agreement with the prediction by van Velzen et al. (2011). We have not included the filter transmission in generating the spectrum. As the filter transmission varies over the wavelength in the given spectral band, and is less than unity, the simulated flux gets reduced, which results in the reduction in the $\delta {t}_{f}$ and hence the detectable rate ${\dot{N}}_{D}.$

India's space mission ASTROSAT, which was launched recently, has a payload Sky Scanning Monitor (SSM) to follow up the transient universe in the X-ray band by scanning nearly half the sky in about 6 hours for a continued same stellar pointing of the spacecraft (http://astrosat.iucaa.in/?q=node/13). The sensitivity of the instrument is $\sim 7.2\times {10}^{-10}\;\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{cm}}^{-2}$ with the integration time of the detector to be 10 minutes. With these parameters, the detectable rate for ASTROSAT is expected to be less than ∼1 yr⁻¹. For the optical surveys in the g band, namely LSST and PAN-STARRS, the TDE may not be resolved and the corresponding rates predicted could be an overestimate by a factor of a few.

We studied in detail, the model of the TDE, taking into account the stellar dynamical and gas dynamical inputs. The overall system parameters include BH mass M_•, specific orbital energy E and angular momentum J, star mass M_⋆ and radius R_⋆, and the pericenter of the star orbit ${r}_{p}(E,J,{M}_{\bullet }).$ We solved the steady state FP equation using the standard loss cone theory for the galactic density profile $\rho (r)\propto {r}^{-\gamma }$ and stellar mass function $\xi (m)$, where $m={M}_{\star }/{M}_{\odot }$, and obtained the feeding rate of stars to the BH $\dot{N}(E,J,m,\gamma )$ that is an increasing function of J and γ, but a decreasing function of E and m. Because the stars evolve along their orbits toward the BH, we compared the lifetime of the main sequence star to the radial period of its orbit and calculated the probability f_⋆ for a star to be captured as a main sequence given by Equation (47). Using this we model the in-fall rate of the disrupted debris, $\dot{M}(E,J,m,t)$, and discuss the conditions for the formation of an accretion disk considering accretion, viscous, ring formation, and radiation timescales. We find that the accretion disk is almost always formed for the fiduciary range of the physical parameters. By equating the peak of $\dot{M}(E,J,m,t)$ to the Eddington rate, we derive the critical black mass ${M}_{c}(E,J,m).$ We have simulated the light curve profiles in the relevant optical g band and soft X-rays for both super- and sub-Eddington accretion disks as a function of $\dot{M}(E,J,m,t)$, taking typical stellar system parameters. Specifically, we found the following key results:

• 1.  
  In Section 2.1, we approximated the radial period of an orbit with Equation (46) and the ${M}_{\bullet }-\sigma$ relation. The radial period of an orbit in Keplerian potential is ${T}_{r}\propto {M}_{\bullet }^{-0.38}{{\mathcal{E}}}^{-3/2}.$ The capture rate ${\dot{N}}_{t}=\int ({N}_{{lc}}({\mathcal{E}})/{T}_{r})$ $d{\mathcal{E}}\propto {M}_{\bullet }^{-0.38}.$
• 2.  
  The applicable ranges of dimensionless energy $\bar{e}$ and angular momentum ℓ are given by {${\bar{e}}_{h}\lt \bar{e}\lt 1,$ $0\lt {\ell }\lt 1$} where ${\bar{e}}_{h}={s}_{t}={r}_{t}/{r}_{h}.$
• 3.  
  We solved the steady state FP equation in Section 2.1 and obtained the capture rate using Equation (49). We found that the capture rate ${\dot{N}}_{t}$ does not show a power law dependence with M_• and it increases with γ. Even though the increase in ${\dot{N}}_{t}$ with γ is non linear, an approximate fit gives ${\dot{N}}_{t}\propto {\gamma }^{p}$, where $p\sim 2.1$ (see Figure 13). For ${M}_{6}\gt 10,\;$ ${\dot{N}}_{t}\propto {M}_{6}^{-\beta }$ and $\beta \sim 0.3\pm 0.01$ (see Figure 12).
• 4.  
  In Section 3, we show that the fractional radius from the star center ${x}_{l}(\bar{e},{\ell },{M}_{\bullet },m)$ to the point where the debris is bound to the BH increases with $\bar{e}$ and ℓ (see Figure 14). The increase with ℓ is significant only for high energy orbits. The peak accretion rate increases with x_l. The decline to later ${t}^{-5/3}$ law is steeper if the energy of the initial orbit is higher, as shown in Figure 15.
• 5.  
  In Section 5, by equating ${\dot{M}}_{p}$ and ${\dot{M}}_{E},$ the critical BH mass ${M}_{c}(\bar{e},{\ell },m)$ is found to increase with $\bar{e}$ and decrease with ℓ and m (see Figure 18). With the decrease in ${\ell },$ the r_p decreases and thus $\dot{M}$ increases, which results in an increase in ${\dot{M}}_{p}.$ For higher $\bar{e}$ and lower ${\ell },\;$ ${M}_{c}(\bar{e},{\ell },m)$ can exceed the BH mass limit for TDE to occur (i.e., ∼ 3 $\times {10}^{8}{M}_{\odot }.$).
• 6.  
  In Section 4, we found that Max[${{\mathcal{T}}}_{r}$ $({10}^{-6}\leqslant \bar{e}\leqslant 1,0\leqslant {\ell }\leqslant 1,1\leqslant {M}_{6}\leqslant 100)$] < 1, which implies that the debris will form an accretion disk (see Figure 16). The ratio ${{\mathcal{T}}}_{v}\lt 1$ for ${M}_{6}\leqslant 31.6$ implies that the accretion disk formed is a slim disk. The ${{\mathcal{T}}}_{v}$ increases with M₆ and decreases with $\bar{e}$ as shown in Figure 17. The higher mass SMBHs form a thin disk from the disrupted debris of a star on low energy orbit and ${\ell }\sim 1$, and a thick disk for a star on a high energy orbit.
• 7.  
  In Section 5.2, we derive the observed flux as a function of $\bar{e}$ and ${\ell }.$ Figure 19 shows the observed fluxes f_obs in the optical g band and the peak observed flux increases with a decrease in ${\ell }.$ The decline of the light curve profile to the later stage gets steeper with increasing ${\ell }.$
• 8.  
  In Section 6, the net detectable rate ${\dot{N}}_{D}$ is calculated for the various missions observing in optical and X-ray bands. Using standard cosmological parameters and mission instrument details, we predict the detectable tidal disruption rates for $\gamma =0.7$ for LSST to be ∼5003 yr⁻¹; Pan-STARRS in the optical g band performing in either the ASS mode or the deep imaging survey mode were predicted to be ∼ 6337 yr⁻¹ for operation in 3 π mode and ∼12.3 yr⁻¹ in the MDS mode, which are in reasonable agreement with scaled-up values based on SDSS detection. Our prediction for eROSITA in the soft X-ray band is about ∼679.5 yr⁻¹, which is consistent with Khabibullin et al. (2014). The values of γ for which our predictions of ${\dot{N}}_{D}$ match with the scaled-up values in van Velzen et al. (2011) are shown as ${\gamma }_{s}$ in Table 2. We also estimated the error in ${\dot{N}}_{D}$ for an error in fit to γ, which is taken to be 0.1 and is also shown in Table 2.
• 9.  
  Our results are in reasonable agreement with the scaled-up values from the SDSS observations (van Velzen et al. 2011), as given in Table 2, along with the detection efficiency of the detector ϒ. The ϒ is lowest for the eROSITA mission due to the high cadence of half year and is highest for Pan-STARRS MDS due to very high sensitivity. The ${\dot{N}}_{D}\propto {\gamma }^{{\rm{\Delta }}}$, where ${\rm{\Delta }}\sim 1.95$ in optical band is shown in Figure 21.

We can use the TDE rate as a proxy to estimate the BH mass distribution as a function of redshift and relate it to the occupation fraction of SMBHs in galaxies (Stone & Metzger 2014). The TDE can influence jet activity from BH systems in two ways. The stellar debris can quench the jet (Gopal-Krishna et al 2008) for BH masses below $3\times {10}^{8}{M}_{\odot }$ or the Poynting flux from a spinning hole can be boosted by the ingestion of a star. There have been several observations (e.g., Burrows et al. 2011) indicating radio transients after a TDE and these have been studied by simulations (Tchekhovskoy et al. 2014). In the future, we plan to include more details of the underlying physics, such as resonant relaxation, non-spherical systems, relativistic gas dynamics, and viscosity prescription, to improve our predictions.

We thank the anonymous referee for insightful and helpful comments that improved the paper significantly. We also thank Andrea Merloni and M.C. Radevi for useful discussions.