How empty is an empty loss cone?

TL;DR

This paper revisits the classical theory of empty loss cones in spherical systems, revealing that the occupancy within the loss cone is flatter than previously thought, affecting the predicted rates of tidal events.

Contribution

It demonstrates that the occupancy distribution in empty loss cones is flatter due to rare large scatterings, revising the exponential decay assumption of classical models.

Findings

Occupancy within the loss cone is approximately flat, not exponentially small.

Tidal events inside the loss cone are nearly as common as those at the edge, with probability decreasing as 1/β.

The results do not affect full loss cone phenomena like tidal disruption events of main sequence stars.

Abstract

We consider two body relaxation in a spherical system with a loss cone. Considering two-dimensional angular momentum space, we focus on "empty loss cone" systems, where the typical scattering during a dynamical time $j_{d}$ is smaller than the size of the loss cone $j_{\rm lc}$. As a result, the occupation number within the loss cone is significantly smaller than outside. Classical diffusive treatment of this regime predict exponentially small occupation number deep in the loss cone. We revisit this classical derivation of occupancy distribution of objects in the empty loss cone regime. We emphasize the role of the rare large scatterings and show that the occupancy does not decay exponentially within the loss cone, but it is rather flat, with a typical value $\sim [(j_d/j_{\rm lc})]^2\ln^{-2}(j_{\rm lc}/j_{\min})$ compared to the occupation in circular angular momentum (where $j_{\min}$ is the smallest possible scattering). Implication are that although the loss cone for tidal break of Giants or binaries is typically empty, tidal events which occurs significantly inside the loss cone ($\beta\gtrsim 2$), are almost as common as those with $\beta\cong 1$ where $\beta$ is the ratio between the tidal radius and the periastron. The probability for event with penetration factor $>\beta$ decreases only as $\beta^{-1}$ rather than exponentially. This effect has no influence on events characterized by full loss cone, such as tidal disruption event of $\sim 1m_\odot$ main sequence star.

Equations

$$\begin{array}[]{l}\frac{\partial f(\vec{j})}{\partial t}=-(t_{d}/2)^{-1}f(\vec{j})\Theta(j_{\text{lc}}-j)+\\ \\ \int d^{2}\vec{j_{\delta}}\ \left[f(\vec{j}-\vec{j_{\delta}})-f(\vec{j})\right]R(\vec{j_{\delta}}).\end{array}$$

$$R(j)=R_0|\overrightarrow{j_{\delta }}{|}^{-4}\text{ for }{j}_{\min }<|\overrightarrow{j_{\delta }}|<j_{\max },$$

$$\begin{array}[]{l}\frac{\partial f(\vec{j})}{\partial t}=-2t_{d}^{-1}f(\vec{j})\Theta(j_{\text{lc}}-j)+\\ \\ R_{0}\int d^{2}\vec{j_{\delta}}\ |\vec{j_{\delta}}|^{-4}\left[f(\vec{j}-\vec{j_{\delta}})-f(\vec{j})\right]\end{array}$$

$$\frac{\partial f}{\partial t}=\frac{1}{j}\frac{\partial }{\partial j}\left(j\cdot D_{\max }\frac{\partial f_{\text{null}}}{\partial j}\right)$$

$$D_{\max }\equiv \frac{1}{4}\int d^2\overrightarrow{j_{\delta }}|j_{\delta }{|}^{2}R(j_{\delta })=\frac{\pi }{2}{R}_0{\Lambda }_{\max }$$

$$f_{\text{null}}(j)=\left\{\begin{array}[]{c c}0&j<j_{0}\\ f_{\text{cir}}(\delta\Lambda^{\max}_{\text{lc}})^{-1}\ln(j/j_{0})&j>j_{0}\end{array}\right.\ .$$

$$\begin{array}[]{l l}\delta\Lambda^{\max}_{\text{lc}}\equiv&\Lambda_{\max}-\Lambda_{\text{lc}}=\\ &\ln(j_{\max}/j_{0})\simeq\ln(j_{\max}/j_{\text{lc}})\end{array}$$

$$\begin{array}[]{c}0=\frac{1}{j}\frac{\partial}{\partial j}\left(j\cdot D_{\max}\frac{\partial f_{\text{exp}}}{\partial j}\right)-2t_{d}^{-1}f_{\text{exp}}(j)\Theta(j_{\text{lc}}-j)\\ \end{array}$$

$$f_{\text{exp}}(j)=\left\{\begin{array}[]{c c}f_{\text{cir}}\cdot\frac{\frac{j_{2}}{j_{\text{lc}}}\exp\left\{\frac{j-j_{\text{lc}}}{j_{2}}\right\}}{\left[\frac{j_{2}}{j_{\text{lc}}}+\delta\Lambda^{\max}_{\text{lc}}\right]}&j<j_{\text{lc}}\\ \\ f_{\text{cir}}\cdot\frac{\left[\frac{j_{2}}{j_{\text{lc}}}+\ln\left(\frac{j}{j_{\text{lc}}}\right)\right]}{\left[\frac{j_{2}}{j_{\text{lc}}}+\delta\Lambda^{\max}_{\text{lc}}\right]}&j>j_{\text{lc}}\end{array}\right.$$

$$j_2\equiv (t_d{D}_{\max }/2{)}^{1/2}.$$

$$0=\frac{1}{j}\frac{\partial }{\partial j}\left(jD(j)\frac{\partial f_{\text{diff}}}{\partial j}\right)-2t_d^{-1}{f}_{\text{diff}}(j)\Theta (j_{\text{lc}}-j)$$

$$D(j)\equiv \frac{\pi }{2}{R}_0\ln \left(\frac{\max \{j-j_{\text{lc}},j_d\}}{j_{\min }}\right)$$

$$j_d\equiv {\left[\frac{\pi }{4}{t}_d{R}_0\ln \left(\frac{j_{\text{lc}}}{j_{\min }}\right)\right]}^{1/2}$$

$$\begin{array}[]{l}f_{\text{diff}}(j)=\\ \\ \left\{\begin{array}[]{c c}f_{\text{cir}}\cdot\frac{\frac{j_{d}}{j_{\text{lc}}}\Lambda_{\text{lc}}^{-1}\exp\left\{\frac{j-j_{\text{lc}}}{j_{d}}\right\}}{\frac{j_{d}}{j_{\text{lc}}}\Lambda_{\text{lc}}^{-1}+\ln\left[1+\Lambda_{\text{lc}}^{-1}\delta\Lambda^{\max}_{\text{lc}}\right]}&j<j_{\text{lc}}\\ \\ f_{\text{cir}}\cdot\frac{\frac{j_{d}}{j_{\text{lc}}}\Lambda_{\text{lc}}^{-1}+\ln\left[1+\Lambda_{\text{lc}}^{-1}\ln\left(\frac{j}{j_{\text{lc}}}\right)\right]}{\frac{j_{d}}{j_{\text{lc}}}\Lambda_{\text{lc}}^{-1}+\ln\left[1+\Lambda_{\text{lc}}^{-1}\delta\Lambda^{\max}_{\text{lc}}\right]}&j>j_{\text{lc}}\end{array}\right.\end{array}$$

$$\frac{{\frac{\partial f_{\text{diff}}}{\partial j}|}_{j_{\text{cir}}}}{{\frac{\partial f_{\text{exp}}}{\partial j}|}_{j_{\text{cir}}}}=\frac{\delta {\Lambda }_{\text{lc}}^{\max }}{{\Lambda }_{\max }\ln \left[1+\left(\frac{\delta {\Lambda }_{\text{lc}}^{\max }}{{\Lambda }_{\text{lc}}}\right)\right]}$$

$$\begin{array}[]{c}f_{\text{$1$step}}(j)=\frac{t_{d}}{2}\int_{j_{d}}^{j_{\max}}d^{2}\vec{j_{\delta}}\ R(\vec{j_{\delta}})f_{\text{diff}}(\vec{j}-\vec{j_{\delta}})=\\ \\ \frac{t_{d}R_{0}}{2}\int_{j_{d}}^{j_{\max}}d^{2}\vec{j_{\delta}}\ |\vec{j_{\delta}}|^{-4}f_{\text{diff}}(\vec{j}-\vec{j_{\delta}})\end{array}$$

$$\begin{array}[]{l l}f_{\text{$1$step}}^{deep}&\cong f_{\text{cir}}\left(\frac{j_{d}}{\Lambda_{\text{lc}}j_{\text{lc}}}\right)^{2}\frac{\left[1-\frac{1}{2\Lambda_{\text{lc}}}+\frac{2j_{d}}{j_{\text{lc}}}\right]}{\frac{j_{d}}{j_{\text{lc}}}\Lambda_{\text{lc}}^{-1}+\ln[1+\Lambda_{\text{lc}}^{-1}\delta\Lambda^{\max}_{\text{lc}}]}\\ &\sim\left(\frac{j_{d}}{\Lambda_{\text{lc}}j_{\text{lc}}}\right)^{2}f_{\text{cir}}\end{array}$$

$$\frac{\pi {\left(\frac{j_{\text{lc}}}{2}\right)}^2\cdot f_{\text{1step}}^{deep}}{2\pi j_{\text{lc}}{j}_d\cdot f_{\text{diff}}(j_{\text{lc}})}\sim \frac{1}{8{\Lambda }_{\text{lc}}}$$

$$\frac{\pi {\left(\frac{j_{\text{lc}}}{2}\right)}^2\cdot f_{\text{diff}}(j=1/2)}{2\pi j_{\text{lc}}{j}_d\cdot f_{\text{diff}}(j_{\text{lc}})}\sim \frac{j_{\text{lc}}}{8j_d}\exp \left\{-\frac{j_{\text{lc}}}{2j_d}\right\}$$

$$t_{\min }\equiv {\left[{\int }_{j_{\min }}^{j_{\max }}{d}^2\overrightarrow{j_{\delta }}R(j_{\delta })\right]}^{-1}$$

$$j_{\min }\sim \left(\frac{m_0}{M_{\text{bh}}}\right)j_{\text{cir}}.$$

$$j_d\sim {\Lambda }_{\text{lc}}^{1/2}{\left(\frac{m_0}{M_{\text{bh}}}\right)}^{1/2}{j}_{\text{cir}}.$$

$$j_{\text{lc}}\simeq {\left(\frac{2r_t}{r_{\text{bh}}}\right)}^{1/2}{j}_{\text{cir}}.$$

$$r_t\simeq {\left(\frac{M_{\text{bh}}}{m_{\ast }}\right)}^{1/3}{R}_{\ast }.$$

$$\frac{j_d}{j_{\text{lc}}}\simeq \frac{{\Lambda }_{\text{lc}}}{2\cdot 0.34}\frac{G{m}_0{m}_{\ast }^{1/3}}{M_{\text{bh}}^{1/3}{R}_{\ast }{\sigma }^2}.$$

$$\frac{j_d}{j_{\text{lc}}}\simeq 0.41\cdot {\left(\frac{{\Lambda }_{\text{lc}}}{10}\right)}^{1/2}{\left(\frac{m_{\ast }}{m_{\odot }}\right)}^{-0.233}{\left(\frac{M_{\text{bh}}}{10^8{m}_{\odot }}\right)}^{-0.4}.$$

$$\frac{j_d}{j_{\text{lc}}}\simeq 0.25\cdot {\left(\frac{{\Lambda }_{\text{lc}}}{10}\right)}^{1/2}{\left(\frac{m_{\ast }}{m_{\odot }}\right)}^{-0.233}{\left(\frac{M_{\text{bh}}}{10^8{m}_{\odot }}\right)}^{-0.4}{\left(\frac{a}{2R_{\ast }}\right)}^{-1/2}$$

$${\Lambda }_{\text{lc}}=\frac{1}{2}\ln \left[{\Lambda }_{\text{lc}}\left(\frac{M_{\text{bh}}}{m_0}\right)\right]+\ln \left(\frac{j_{\text{lc}}}{j_d}\right)$$

$${\Lambda }_{\text{cir}}=\ln \left(\frac{M_{\text{bh}}}{m_0}\right)<2{\Lambda }_{\text{lc}}$$

$$\beta \equiv \frac{\text{tidal radius}}{\text{peri-astron distance}}={\left(\frac{j_{\text{lc}}}{j}\right)}^2.$$

$${\gamma }_{net}^{\text{small}}(j_{\text{scat}},j_{\text{targ}})=2\pi j_{\text{targ}}{D}_{j_{\text{scat}}}{\frac{\partial f}{\partial j}|}_{j_{\text{targ}}}\sim R_{0}f(j_{\text{targ}})$$

Full text

We revisit this classical derivation of occupancy distribution of objects in the empty loss cone regime. We emphasize the role of the rare large scatterings and show that the occupancy does not decay exponentially within the loss cone, but it is rather flat, with a typical value $\sim[(j_{d}/j_{\text{lc}})]^{2}\ln^{-2}(j_{\text{lc}}/j_{\min})$ compared to the occupation in circular angular momentum (where $j_{\min}$ is the smallest possible scattering).

Implication are that although the loss cone for tidal break of Giants or binaries is typically empty, tidal events which occurs significantly inside the loss cone ($\beta\gtrsim 2$), are almost as common as those with $\beta\cong 1$ where $\beta$ is the ratio between the tidal radius and the periastron. The probability for event with penetration factor $>\beta$ decreases only as $\beta^{-1}$ rather than exponentially. This effect has no influence on events characterized by full loss cone, such as tidal disruption event of $\sim 1m_{\odot}$ main sequence star.

††pubyear: 2016††pagerange: How empty is an empty loss cone?–LABEL:lastpage

1 Introduction

Many physical phenomena, especially in the dynamics around the galactic center involve destructive processes. Perhaps the simplest problem is that of tidal disruption of a star when it is close enough to the galactic center massive black hole. Orbits which penetrate that distance are called loss cone orbits and are characterized by low specific angular momentum. If orbital elements of objects around the galactic center black hole were fixed, all loss cone orbits have already been destroyed and the loss cone would be completely empty.

Many previous works have discussed the loss cone filling rate (e.g. Frank & Rees, 1976; Lightman & Shapiro, 1977; Cohn & Kulsrud, 1978; Magorrian & Tremaine, 1999; Yu & Tremaine, 2002, 2003; Milosavljević & Merritt, 2003), and see Merritt (2013) for a comprehensive discussion. Most of these works assume that loss cone filling is well described by the Fokker-Planck equation, i.e by diffusion in the angular momentum - energy phase space. The rational is that if the typical change in energy or angular momentum during a dynamical time are small compared to the size of the loss cone, diffusion is an accurate description. However, it is well know that hard scatterings contributes equally to the diffusion coefficient. Such scatterings are not captured correctly by diffusion or Fokker-Plank equation (Bar-Or et al., 2013).

To address that, we start in section §2 with a Boltzmann equation, which correctly describes small as well as large scatterings. We go over a series of possible approximate solution in §3. Finally in §4 we compare those approximation with numerical experiments which simulate the full Boltzmann equation.

2 The Boltzmann equation with loss cone

We consider the evolution due to two body scattering of the otherwise conserved values of the stellar orbits, i.e the angular momentum $\vec{j}$ and the energy $E$. The loss cone is represented by $||\vec{j}||<j_{\text{lc}}(E)$, where $j_{\text{lc}}(E)$ is the angular momentum related to some critical periastron.

Since scatterings driven by gravitational interactions makes stars enter the loss cone mainly due to angular momentum loss (Lightman & Shapiro, 1977; Merritt, 2013), we treat the whole dynamics as two dimensional $\vec{j}$ phase space (two tangential components of the angular momentum), and ignore the energy. Moreover, since the system have polar symmetry, the distribution function only depends on the magnitude of $\vec{j}$ ($0\leq j\lesssim j_{\text{cir}}$).

We define the occupancy distribution $f(t,\vec{j})$ to be proportional to the number of stars at time $t$, per unit angular momentum square around angular momentum $\vec{j}$, i.e $f(t,\vec{j})d^{2}\vec{j}$ is proportional to the number of stars with angular momentum $\vec{j}$. Therefore $f$ has units of $1/j^{2}$. We use the normalization $f(j_{\text{cir}})\equiv f_{\text{cir}}$ where $f_{\text{cir}}$ is approximately the total number of objects in the system over $4\pi j_{\text{cir}}^{2}$.

We account for two processes: scattering and destruction, and focus on orbits with low angular momentum compared to that of a circular orbit. This is an excellent assumption for the galactic center where loss cone orbits have angular momentum over a thousand time smaller than those of circular orbits. We assume that the rate of scattering by an amount $R(\vec{j_{\delta}})$, is independent of the current angular momentum of the object . This is justified since for the almost radial orbits most of the scattering occurs far from pericenter and the exact position of the pericenter is thus unimportant. The destruction process, eliminates once every dynamical time objects of low enough angular momentum $j<j_{\text{lc}}$. Exact description of such process requires following the orbital phase of each star. To avoid that, and have the distribution function depend on angular momentum only, we use the following approximate Boltzmann equation:

[TABLE]

where $t_{d}$ is the orbital time, or dynamical time, of the particles, $\Theta(j_{\text{lc}}-j)$ is the Heaviside function which gets $1$ for $j<j_{\text{lc}}$ and [math] for $j>j_{\text{lc}}$ and $R(\vec{j_{\delta}})$ is the differential rate of scatterings of size and direction $\vec{j_{\delta}}$, i.e. number of scattering per unit time per $d^{2}\vec{j_{\delta}}$.

The sink term $-(t_{d}/2)^{-1}f(\vec{j})\Theta(j_{\text{lc}}-j)$, which takes the loss process into account, and especially the factor of $1/2$, warrants some discussion. A star may get scattered into the loss cone, $j<j_{\text{lc}}$, at any phase of its orbit. It will be lost when its orbital phase brings it close enough to the central object. To mimic that without following the phases of individual stars, we have to make an assumption on the phase distribution for particles in the loss cone.

If scattering in and out of the loss cone is faster than the orbital time, the loss cone will be uniformly populated in orbital phase. In such a case, the rate of destruction of particles is equal to the number of particles in the loss cone divided by the orbital time. Therefore, the sink term in this limit is $-t_{d}^{-1}f(\vec{j})\Theta(j_{\text{lc}}-j)$.

However, we are interested in the opposite limit, where scattering deep into the loss cone are rare. In this “empy loss cone" limit, a particle that was scattered deep into the loss cone will remain there until it gets destroyed If we assume that it entered the loss cone with an arbitrary phase, then a star survives within the loss cone for half its orbital period on average. The sink term is therefore given by $-(t_{d}/2)^{-1}f(\vec{j})\Theta(j_{\text{lc}}-j)$, which is what we used in our equation 1.

Our choice of the sink term, mimics well the behavior deep inside the LC and allows us to find the occupancy distribution deep inside the LC. Close to the LC’s edges on the other hand, where the typical star may scattered inward and outward of the LC during one dynamical time, this sink term differs from the real physical behavior by numerical factor of order $\sim 2$. We therefore expects the calculated occupancy distribution to be a good approximation inside the LC, and to err by some on its edges.

We find the scattering’s rate term $R(\vec{j_{\delta}})$. We consider small LC, i.e. $j_{\text{lc}}$ much smaller than the circular angular momentum $j_{\text{cir}}$ . This is the relevant case for HVSs or tidal disruption events. In this case, the stars which are about to get scattered into the loss cone, are moving almost radially with respect to the central BH during most of their period. Their scatterers on the other hand, are the more numerous stars moving on much more circular orbits.

In this case, of gravitational scatterings of stars with low angular momentums, the velocity’s change is in the direction perpendicular to the BH’s direction. Therefore $j_{\delta}\equiv a\Delta v_{\perp}\simeq Gma/bv$. where here $\Delta v_{\perp}$ is the velocity change in the perpendicular direction, $m$ is the mass of an individual scatterer, $a$ is the distance from the central BH and $v$ the relative velocity of the encountering objects. The rate of such encounters is proportional to $b^{2}$, so that the rate of encounters resulting in angular momentum change of order $j_{\delta}$ is proportional to $j_{\delta}^{-2}$.

Moreover, since stars only interact with stars in their close surrounding, we may assume that in any phase of the orbit, the subject star interact with an isotropic background. As a result, the probability that an individual subject star experience a scattering due to an impact parameter $b$, does not depend on the direction. Thus, $R(\vec{j_{\delta}})$ which is the rate per unit $d^{2}\vec{j_{\delta}}$ is isotropic and proportional to $||\vec{j_{\delta}}||^{-4}$.

We therefore define

[TABLE]

where $R_{0}$ is some function of the masses of the particles and the velocity dispersion and $j_{\min},\ j_{\max}$ are the minimal and the maximal possible scattering: $j_{\min}$ corresponds to an impact parameters of order the distance between stars, and $j_{\max}\approx j_{\text{cir}}$ corresponds to an impact parameter that results in a velocity change of order unity

Substitute the rate term (equation 2) into the Boltzmann equation (equation 1) our Boltzmann equation becomes

[TABLE]

It is the process described by this equation that we simulate numerically in §4.

3 Solution to the Boltzmann equation

The system we describe has no steady state solution, since it contains a sink (the loss cone) without any specified supply term. We may force a steady state solution to the Boltzmann equation by normalizing $f(j_{\text{cir}})\equiv f_{\text{cir}}$ independent of time. Physically this could be thought as a result of supply of new objects to the outskirts of the loss cone, e.g. from orbits of higher energy which we do not consider, or an approximation of low angular momenta $j\ll j_{\text{cir}}$ whose intrinsic evolution should have been fast.

In this section, we list approximate solutions to equation 3 with different levels of accuracy. We denote the solution sunder these approximate solutions by $f_{\text{null}}(j)$, $f_{\text{exp}}(j)$, $f_{\text{diff}}(j)$ and f_{\text{1step}}(j) while the colors: blue (dashed), green (dashed), orange and red are used to mark those four solutions respectively in all figures.

3.1 Diffusion with null loss cone - $f_{\text{null}}(j)$

The diffusion approximation assumes that scatterings can be considered small, i.e., the distribution function $f(\vec{j})$ does not change significantly over the size of the scattering. In this case the integrand on the right hand side of equation 3 may be expanded using Taylor approximation and equation 3 for $j>j_{\text{lc}}$ may be rewritten as

[TABLE]

where:

[TABLE]

Where $\Lambda_{\max}\equiv\ln(j_{\max}/j_{\min})$ 111 We use the notaion $\Lambda_{i}$ to define $\ln(j_{i}/j_{\min})$ at other cases as well: $\Lambda_{d}\equiv\ln(j_{d}/j_{\min})$, $\Lambda_{\text{lc}}\equiv\ln(j_{\text{lc}}/j_{\min})$, $\Lambda_{\text{cir}}=\ln(j_{\text{cir}}/j_{\min})$. and $\Lambda_{\max}=\ln(j_{\max}/j_{\min})$. We also use the notation $\delta\Lambda^{\max}_{\text{lc}}$ to define $\ln(j_{\max}/j_{\text{lc}})\equiv\Lambda_{\max}-\Lambda_{\text{lc}}$ . The expression for the diffusion coefficient $D_{\max}$ shows that scattering of all scale contribute equally to the diffusion, leading to the familiar coulomb logarithm $\Lambda_{\max}=\ln(j_{\max}/j_{\min})$.

We find the steady state solution of equation 4. In systems where the typical scattering over the dynamical time is smaller than the size of the loss cone, the loss cone is relatively empty. In the current approximation we therefore take the boundary condition $f_{\text{null}}(j_{0})=0$, where $j_{\text{lc}}-j_{0}\ll j_{\text{lc}}$. (Lightman & Shapiro, 1977; Magorrian & Tremaine, 1999; Merritt, 2013).

This approximation gives us a steady state solution of the form

[TABLE]

where

[TABLE]

The solution of the form $f(j)\propto\ln(j/j_{0})$ was found by Lightman & Shapiro (1977) (see also Merritt (2013)). The blue dashed line in figures 1,2 and 3 marks $f_{\text{null}}(j)$. Here, we normalized $f$ to $f_{\text{cir}}$ at some large angular momentum $j_{\max}$. Typically $j_{\max}\sim j_{\text{cir}}$, if more than the relaxation time has elapsed and $j_{\max}\sim(Dt)^{1/2}$ at earlier times.

This solution describes a constant flux of objects towards the loss cone, given by $\pi^{2}R_{0}f_{\text{cir}}\Lambda_{\max}/\delta\Lambda^{\max}_{\text{lc}}$. However, it was derived using the zero boundary conditions. By using this boundary condition, we give up the ability to inspect the distribution of particles within the loss cone. The question of how empty is an empty loss cone, how deep particles penetrate into it (what is exactly $j_{\text{lc}}-j_{0}$?) and what is the distribution of $j$ of particles that suffer distruction can not be answered in this framework.

3.2 Diffusion with Delayed Distractio - $f_{\text{exp}}$

The first correction, which was already used by Cohn & Kulsrud (1978), allows to account for particles inside the loss cone. We drop the null boundary condition $f(j_{\text{lc}})=0$, and instead take a sink term inside the loss cone just as in the original Boltzmann equation (eq. 3). The resulting diffusion equation for steady state is

[TABLE]

This could be easily solved both inside and outside the loss cone. The two solutions are then matched to ensure continuity of $f$ and $\partial f/\partial j$ at $j_{\text{lc}}$. The solution inside the loss cone is given by a Bessel function, but since it decays so strongly within the loss cone, we ignore the cylindrical geometry and approximate the Bessel function as an exponential. We obtain

[TABLE]

with

[TABLE]

Note that the parameter $j_{2}$ is alike Lightman & Shapiro (1977)’s $j_{2}$ which they define as $(t_{d}/t_{\text{relax}})^{1/2}j_{\text{cir}}$, since $t_{\text{relax}}=j_{\text{cir}}^{2}/D_{\max}$. Moreover, $f_{\text{exp}}(j)$ outside the LC, may be written as $\ln(j/j_{0})$, where $j_{0}\equiv j_{\text{lc}}\exp(-j_{2}/j_{\text{lc}})\cong j_{\text{lc}}-j_{2}$, i.e $j_{2}$ is the typical penetration of particles into the LC. The green dashed line in figures 1,2 and 3 marks $f_{\text{exp}}(j)$.

We compare our occupancy $f_{\text{exp}}$ close to the loss cone’s edge with those of Cohn & Kulsrud (1978) in the empty LC regime $j_{2}\ll j_{\text{lc}}$. The functional form of the solution is similar. However, their penetration of particles into the loss cone ($j_{\text{lc}}-j_{0}$ where $j_{0}$ is set according to the $\propto\ln(j/j_{0})$ behavior of the profile outside the LC) is larger by a numerical factor of $1.165$ than ours222Milosavljevic solution of Cohn & Kulsrud (1978)’s Fokker-Planck equation as presented in Merritt (2013), results in a factor of $(\pi/2)^{1/2}\cong 1.25$ . The reason for the difference is that at depth of $\lesssim j_{2}$ inside the loss cone, our assumption that the average time that a particle spends inside the loss cone is $t_{d}/2$ (used in equation 3 and 9, and therefore in the definition of $j_{2}$ above) is not accurate (see discussion below equation 1). A star that resides very close to the boundary of the loss cone, has a chance to be scattered in and out of it within a dynamical time. Their phases are therefore more uniform and the factor of $2$ we have used should be closer to unity.

3.3 $f_{\text{diff}}(j)$ - Variable diffusion coeficient

The diffusion approximation is valid as long as the steps size is smaller than the scale over which the function $f$ changes. Our solution in §3.2 shows that far outside the loss cone, the scale over which the function $f$ changes is roughly the distance from the loss cone, while very close to the loss cone edge, it changes on a scale $j_{2}$.

Higher orders of the Fokker-Planck equation can partially account for this effect (Bar-Or et al., 2013). However, where these effects are strong, infinite numbers of terms would need to be considered. We therefore take another approach. In appendix A, we argue that outside the loss cone, large scatterings bigger than this scale, turn out to be less important.

Therefore, taking into account only the relevant small scatterings, reduces the coulomb logarithm that appears in the diffusion coefficient. This leads as to an effective diffusion equation with a variable diffusion coefficient:

[TABLE]

with

[TABLE]

where $j_{d}$ is the scale on which $f(j)$ is changing close to the LC surrounding. Notice that due to the logarithmic dependency of $D(j)$, $f(j)$ close to the LC is changing on scale shorter than $j_{2}=(t_{d}D_{\max}/2)^{1/2}$ and therefore $j_{d}<j_{2}$. We define $j_{d}\equiv(t_{d}D(2j_{\text{lc}})/2)^{1/2}$ or

[TABLE]

In this description, we have not focused on the best estimate for the diffusion coefficient inside the loss cone, since as we shall see in the next section, the loss cone is mostly populated by non diffusive processes.

Solving this equation inside and outside the loss cone under continuity conditions, we find the steady state solution to this equation to be approximately

[TABLE]

Our solution 14 is approximate. It is an exact solution to equation 11 with diffusion coefficient of the form $D(j)=(\pi R_{0}/2)\cdot\ln(j/j_{\min})$, rather than the diffusion coefficient presented in 12. As a result, the solution is inaccuraate very close to the LC edge where $0<|j-j_{\text{lc}}|\lesssim(j_{\text{lc}}j_{\min})^{1/2}$. It is, however, quite accurate elsewhere. The orange line in figures 1,2 and 3 marks $f_{\text{diff}}(j)$.

Note that very far from the loss cone, where $\ln(j/j_{\text{lc}})\gg\Lambda_{\text{lc}}$, $f_{\text{diff}}(j)$ profile satisfy $f(j)\propto\ln[\ln(j/j_{\text{lc}})]$ rather than $f(j)\propto\ln(j/j_{\text{lc}})$. This flattening of the profile compared to the expectation from constant diffusion coefficient ($f_{\text{exp}}(j)$), comes from the logarithmic reduction of the diffusion coefficient at low $j$’s: $D\sim R_{0}\ln(j_{\text{lc}}/j_{\min})$ rather than $D\sim R_{0}\ln(j_{\max}/j_{\min})$. This reduction, dictate lower flux also at large $j$’s where $D(j)\sim D(j_{\text{cir}})\sim D_{\max}$.

For similar normalization of $f(j_{\text{cir}})\equiv f_{\text{cir}}$, we find the ratio between the flux according to $f_{\text{diff}}$ and the flux according to $f_{\text{exp}}$ by comparing $\partial f_{\text{diff}}/\partial j|_{j_{\text{cir}}}$ and $\partial f_{\text{exp}}/\partial j|_{j_{\text{cir}}}$. We obtain:

[TABLE]

This reduction can never be more significant than order unity. If $\delta\Lambda^{\max}_{\text{lc}}\ll\Lambda_{\text{lc}}$, the reduction is simply the $\Lambda_{\text{lc}}/\Lambda_{\max}$ equivalent to the diffusion coefficient reduction in the LC surroundings. If on the other hand $\delta\Lambda^{\max}_{\text{lc}}\gg\Lambda_{\text{lc}}$, the reduction is $\sim\ln(\Lambda_{\max}/\Lambda_{\text{lc}})$ - which for realistic parameter can not be larger than $\sim 3$.

Specific examples are given in §5.

3.4 f_{\text{1step}} - The effect of large scatterings

The delayed destruction by the loss cone, predicts exponentially small amount of particles deep inside the loss cone. The typical decay length of the exponent, is $j_{d}$, the typical diffusion during half dynamical time (see equation 9), which is necessarily $\ll j_{\text{lc}}$ in the case of empty LC. On the other hand, scatterings of size $j\sim j_{\text{lc}}\gg j_{d}$ are rare only by factor $\sim(j_{\text{lc}}/j_{d})^{2}$ compared to $j_{d}$ scatterings. Via those scatterings, particle may enter anywhere in the loss cone with roughly constant probability. Therefore, we expect that inside the loss cone, at depth larger than $\sim j_{d}$, the profile is dominated by particles that entered into the loss cone due to one large scattering rather than by diffusion on time scales $\sim t_{d}$.

We define f_{\text{1step}}(j) as the steady state solution $f_{\text{diff}}(j)$ evolved by a single step taken from the probability distribution of steps over time $t_{d}/2$ (the average time particles survive inside the loss cone). Since the fraction of particles that preform scattering of size and direction $\vec{j_{\delta}}$ during time $t_{d}/2$, is $R(\vec{j_{\delta}})t_{d}/2$ (true only for $j\gg j_{d}$), we find that deep inside the loss cone, at $j$’s satisfies $j_{\text{lc}}-j\gg j_{d}$, the profile is:

[TABLE]

This equation may be solved numerically everywhere inside the loss cone and this solution noted by the red dashed line in figures 1, 2 and 3.

Specifically, in the central part of the loss cone ($j\ll j_{\text{lc}}$), one may estimate f_{\text{1step}}(j), by substituting $j=0$ in the integrand and solve the integrals. Than, in first order in $\Lambda_{\text{lc}}^{-1}$ and $j_{d}/j_{\text{lc}}\ll 1$ the expression for f_{\text{1step}} may be approximated as:

[TABLE]

This value is marked by the horizontal dashed-doted black line in figures 1 and 3.

One major consequence of the last result is that the number of objects in the deep LC (say $j<j_{\text{lc}}/2$) is lower only by a factor of order unity than the number of objects on the LC’s edge:

[TABLE]

instead of the more significant factor

[TABLE]

predicted by diffusive model.

3.5 Analytical solutions - Summary

Although we could not find an exact solution to equation 3, we provide approximate solutions in different regimes.

Outside the loss cone and at depth of $\lesssim j_{d}\equiv[(\pi/4)t_{d}R_{0}\ln(j_{\text{lc}}/j_{\min})]^{1/2}$ into the loss cone the profile is dominated by diffusive behavior with delayed destruction. Therefore the solution is approximately $f_{\text{diff}}(j)$ (equation 14). Deeper than $\sim j_{d}$ inside the loss cone, only large scatterings contribute to the population and therefore the solution is approximately f_{\text{1step}}(j) (equation 16).

Since this solution takes into account large scatterings which are not considered in a Fokker-Planck equation, it is qualitatively different than previous solutions. In particular:

• •

The profile far outside the LC, in angular momentums $j\gg j_{\text{lc}}^{2}/j_{\min}$ (if such angular momentums are yet $\lesssim j_{\text{cir}}$), is not $\propto\ln(j/j_{\text{lc}})$, but flatter, with asymptotic behavior $\propto\ln[\ln(j/j_{\text{lc}})]$ (equation 14). As a result, the LC filling rate is reduced by logarithmic factor (equation 15)

• •

The profile deep inside the LC is not exponentially small with a scale of $j_{d}$, but rather polynomially $\sim[j_{d}/j_{\text{lc}}\ln(j_{\text{lc}}/j_{\min})]^{2}f_{\text{cir}}$.

Special attention to large scatterings was already discussed in Bar-Or et al. (2013). Bar-Or et al. (2013) focussed on energy rather than angular momentum relaxation. They used higher order terms of the Master equation (or Boltzmann equation) than those leading to the familiar Fokker-Planck equation.

4 Numerical Experiments

We compare our analytical results to numerical Monte-Carlo simulation results in order to demonstrate that $f_{\text{diff}}(j)$ and f_{\text{1step}}(j) are more accurate than the previous solutions. In this simulation, we let a large number of particles $N_{p}$ evolve over time from a flat distribution to an effective steady state distribution. Each particle starts from some random location in $j$’s space and earn random scatterings over long time. We run this simulation over large virtual time $t_{\text{sim}}$, so that a quasi steady state is achieved in some finite region $j<j_{\text{sim}}\equiv[t_{\text{sim}}D(2j_{\text{lc}})]^{1/2}$ while $j_{\text{lc}}\ll j_{\text{sim}}$).

4.1 Simplified simulation description

We employ the following simple algorithm for finding the steady state solution of equation 3 in low angular momentums ($j<j_{\text{sim}}$). In practice, since this algorithm is slow, we employ an accelerated version of it as described in appendix §B. We presented here the simplified (non accelerated) algorithm for simplicity.

• •

A single point particle has initial angular momentum $j<j_{\max}$ and a random direction. The probability for the particle to start in a ring $\{j,j+dj\}$ is proportional to $j$.

• •

Every step each particle preform a single jump. The direction of the jump is randomly picked with equal probability to each direction. The size of the jump $j$ varies between $j_{\min}$ and $j_{\max}$ with probability $\propto j^{-3}$ (corresponding to $R(\vec{j_{\delta}})\propto j_{\delta}^{-4}$ - equation 2). We change the position of the particle according to the chosen size and direction of the jump.

After each jump, the time is advanced by a constant

[TABLE]

which correspond to the time it takes to a particle to make a step larger than $j_{\min}$ according to equation 2.

• •

Whenever a particle has penetrated to the loss cone, i.e. it has $j<j_{lc}$, it has a probability of $1-e^{-2t_{\min}/t_{d}}$ to be expelled. This mimics the loss cone term which appears in the Boltzmann equation (equation 3). If the particle has been expelled, it does not enter into the statistics of final locations.

Notice that as has been mentioned in the discussion following equation 1, such a sink term mimics the true physics deep inside the LC (at depth larger than $\sim j_{d}$). Thus, since a star that have entered those regions, typically does not scattered out of the LC, but rather remains there until its destruction.

Since this simulates the angular momentum evolution over time $t_{\text{sim}}$, it effectively gives the steady state profile up to angular momentums $\gtrsim j_{\text{sim}}\equiv[t_{\text{sim}}D(2j_{\text{lc}})]^{1/2}$.

4.2 Simulation parameters

The parameters of the simulation all in units of $j_{\text{lc}}$ and $t_{d}/2$ are as follows:

• •

$j_{d}$ - the typical scattering at the LCs surrounding, which determents how empty the LC is. The value of $j_{d}$, set the rates of all scattering via its relation to $R_{0}$ (equation 13).

• •

$j_{\min}$ - the minimal possible scattering.

• •

$j_{\max}$ - the maximal possible initial angular momentum. We use $j_{\max}$ also as the maximal possible scattering. Ideally $j_{\max}\rightarrow\infty$

• •

$t_{\text{sim}}\text{ or }j_{\text{sim}}$ - $t_{\text{sim}}$ is the virtual time we let the simulation evolve. We define $j_{\text{sim}}\equiv[t_{\text{sim}}D(2j_{\text{lc}})]^{1/2}$ and demand $j_{\text{lc}}\ll j_{\text{sim}}\ll j_{\max}$.

• •

$N_{p}$ - Number of particles. This number should be as large as possible in order to ensure large number of particles deep in the LC.

We have preformed three simulations for the following parameters:

[TABLE]

Simulations A and B represent the same physics (same $j_{\min},j_{d}$). The main difference between those two, is at $t_{\text{sim}}$ and correspondingly $j_{\text{sim}}$ and at $N_{p}$. Simulation A is designed to give the details of the profile at the loss cone surrounding while simulation B gives an insight on relatively large distances (the deviations from $f_{\text{null}}$ and $f_{\text{exp}}$ in high $j$’s). Simulation C demonstrates the case of a mildly empty loss cone.

4.3 results

The results of the three simulations, and the theoretical models for the relevant parameters, are all presented in figures 1, 2 and 3. In general, the suggested model, which is represented by the $f_{\text{diff}}(j)$ outside the loss cone and close to its edge and f_{\text{1step}}(j) at depths $\gtrsim j_{d}$ into it (equations 14 and 16), agrees with the data.

Figure 1, shows the occupancy distribution in the LC’s vicinity. Our main result is that inside the LC, at $j_{\text{lc}}-j\gg j_{d}$, all the diffusive profiles ($f_{\text{null}},f_{\text{exp}},f_{\text{diff}}$) fail to describe the actual profile of simulation A. As has been claimed in §3.4 the profile in this region is dominated by objects who have entered the LC via one large scattering and indeed f_{\text{1step}}(j) describes the profile in those regions well.

As expected, the profile in the LC’s edges surroundings is best described by $f_{\text{diff}}(j)$. $f_{\text{null}}(j)$ under estimates the profile at distance of $\sim j_{d}$ from the LC’s edge, while $f_{\text{exp}}(j)$ overestimate it. Figure 1 also shows a flattening of the profile in the high $j^{\prime}s$ ($j\gg j_{\text{lc}}$) with respect to $f_{\text{null}}$ or $f_{\text{exp}}$. To better capture this feature, we have preformed another simulation (Simulation B) which has the same physical parameters ($j_{d},j_{\min}$) but a much larger $j_{\text{sim}}$ (on the price of a much smaller number of particles). Figure 2 presents the results of this simulation. Notice the profile deviate significantly from any constant diffusion coefficient model ($f_{\text{null}}$ or $f_{\text{exp}}$). $j$’s dependent diffusion coefficient on the other hand, led us to $f_{\text{diff}}(j)$ and that profile seems to agrees well with the actual profile.

Figure 3 present the results of simulation C. It shows the applicability of our approximations in a more marginal empty LC’s case, were $j_{d}/j_{\text{lc}}=0.23$.

5 Realistic examples

In sections 3.3 and 3.4 we have described two effects. The first, is the flux reduction as a consequence of the lower diffusion coefficient (equation 15) and the second is the presence of objects deep inside the LC, even if $j_{d}\ll j_{\text{lc}}$. The significance of the described effects, depends upon the values of $j_{\min},j_{d}$ and $j_{\text{lc}}$, we therefore discuss the values of these parameters in realistic systems.

5.1 Parameters in an equal mass realistic systems

Consider a massive black hole of mass $M_{\text{bh}}$, surrounded by spherical bulge. We assume that the massive black hole’s sphere of influence contains $\sim M_{\text{bh}}/m_{0}$ objects of mass $m_{0}$.

While experiencing a scattering of impact parameter $b$ with an object of mass $m_{0}$, the typical object changes its angular momentum by $\delta j\sim j_{\text{cir}}\cdot Gm_{0}/bv^{2}$. Since scattering at the sphere of influence, $r\sim r_{\text{bh}}$ dominate the scatterings, we find that the minimal possible scattering is:

[TABLE]

The number of scatterings of size $j_{\min}$ during one dynamical time is roughly the number of the stars ($M_{\text{bh}}/m_{0}$). Therefore, the contribution from scatterings of size about $j_{\min}$ to the diffusion coefficient is $\sim(M_{\text{bh}}/m_{0})j_{\min}^{2}t_{d}^{-1}$. Since other step sizes have equal contribution we find from equation (13)

[TABLE]

The loss cone size, $j_{\text{lc}}$, depends upon the ratio between the tidal radius, $r_{t}$, and the radius of the sphere of influence, $r_{\text{bh}}$:

[TABLE]

where for tidal disruption event (TDE) of star of mass $m_{*}$ and radius $R_{*}$

[TABLE]

Comparing $j_{\text{lc}}$ with $j_{d}$ taking numerical factor exact for isothermal sphere 333 We have used $j_{d}=(\Lambda_{\text{lc}}/\Lambda_{\max})^{1/2}j_{2}$ where $j_{2}=(t_{d}/t_{\text{relax}})^{1/2}j_{\text{cir}}$ and $t_{\text{relax}}$ is the Chandrasekhar relaxation time $t_{\text{relax}}=0.34\sigma^{3}/G^{2}m_{0}\rho$ (Lightman & Shapiro, 1977). Than, assuming an isothermal sphere, we have substituted $r_{\text{bh}}=GM_{\text{bh}}/\sigma^{2}$ and $\rho=M_{\text{bh}}/2\pi r_{\text{bh}}^{3}$ (Merritt, 2004). , we find:

[TABLE]

Loss cone in this case turn out to be only marginally empty considering a typical main sequence star. For a main sequence star $R_{*}\simeq(m_{*}/m_{\odot})^{0.8}R_{\odot}$ and an $M-\sigma$ relation of the form $M_{BH}=M_{0}\left(\sigma/200\text{km s}^{-1}\right)^{\beta}$ where $M_{0}=1.32\cdot 10^{8}m_{\odot},\ \beta=4.24$ (Gültekin et al., 2009) we find

[TABLE]

In principle this may lead to empty LC at large $M_{\text{bh}}$, however, since the maximal $M_{\text{bh}}$ for which a TDE occurs outside the horizon is $\simeq 10^{8}m_{\odot}\cdot(m/m_{\odot})^{0.7}$, empty LC’s scenario is only marginally relevant only for relatively massive stars and MBHs. However for Giants, empty LC scenario become relevant. Since the radius of a giant of mass $m_{\odot}$, is $\gtrsim 100R_{\odot}$, $r_{t}$ also grows by factor of $\sim 100$ and $j_{\text{lc}}$ grows by factor of $10$. This makes $j_{d}/j_{\text{lc}}\sim 0.04$ (LC emptier than presented in figure 1) 444

Notice however that although $j_{d}/j_{\text{lc}}$ in the simulation is approximately similar to that of the realistic system, $j_{\min}/j_{\text{cir}}$ in the simulation is much larger than in the reality (reduction of $j_{\min}$ and enhancement of the maximal $j$ for which the system is relaxed is numerically expansive). As a result, $\delta\Lambda^{\max}_{\text{lc}}$ and $\Lambda_{\text{lc}}$ in a real system are larger by a factor of $\simeq 2$ than in the simulation.

For binary’s breakup which leads to hypervelocity stars (HVSs) production, $r_{t}\simeq(M_{\text{bh}}/2m_{*})^{1/3}a$ where $a$ is the binary’s separation and $2m_{*}$ is the binary’s mass. In this case:

[TABLE]

and therefore the LC may be empty also for low mass MBHs such as Sagittarius A* or for low mass stars. For instance, a binary of B stars ($4m_{\odot}$ each, see Brown et al. (2009)) and separation of $15R_{*}$, which interact with an MBH of mass $4.3\cdot 10^{6}m_{\odot}$, has a mildly empty LC with $j_{d}/j_{\text{lc}}\simeq 0.23$ (same as in 3); and a binary of two solar type stars with separation of $15R_{\odot}$, interacting with Andromeda’s MBH ($M_{\text{bh}}\simeq 1.2\cdot 10^{8}m_{\odot}$) has $j_{d}/j_{\text{lc}}\simeq 8.4\cdot 10^{-2}$ (smaller than $j_{d}/j_{\text{lc}}$ of figure 1).

5.2 LC filling rate reduction in a realistic system

Notice that the relation between $j_{\min}$ and $j_{d}$ as given by equations 21 and 22, restrict the LC’s influx rate reduction presented in equation 15. According to equation 15, the flux reduction is $\sim\Lambda_{\text{lc}}/\Lambda_{\text{cir}}$ in the regime we later find to be relevant of $\Lambda_{\text{lc}}\gtrsim\delta\Lambda^{\max}_{\text{lc}}$. Since our model only assumes empty LC we require $j_{d}\ll j_{\text{lc}}$ and obtain:

[TABLE]

therefore, in a system in which all objects have the same mass, the influx reduction (for $\ln(j_{\text{lc}}/j_{\min})\ll\ln(M_{\text{bh}}/m_{0})$) is about a factor of $2$.

6 summery and applications

We discussed the effect of large scatterings on the occupancy distribution of stars in the presence of an empty loss cone. We expect the occupancy distribution to be of the form of $f_{\text{diff}}$ (equation 14) outside the LC and at depth $\lesssim j_{d}$ into it, and f_{\text{1step}} (equation 16) at larger depths inside the LC. Those expectations were deduced from an analytical model and verified using Monte Carlo simulations.

Our solution, contain two major differences compared to the standard Fokker-Planck approach (Lightman & Shapiro, 1977; Cohn & Kulsrud, 1978; Magorrian & Tremaine, 1999; Merritt, 2013).

First, at large angular momentums $j\gtrsim j_{\text{lc}}^{2}/j_{\min}$ (if those exist) the profile is not $\propto\ln(j/j_{\text{lc}})$ as expected from a Fokker-Planck approximation, but flatter, with an asymptotic behavior $\propto\ln[\ln(j/j_{\text{lc}})]$ (equation 14). This flatter profile make the flux into low angular momentums a little lower. More accurately, if $j_{\text{cir}}<j_{\text{lc}}^{2}/j_{\min}$ the Coulomb logarithm $\ln(j_{\max}/j_{\min})$ should be replaced with $\ln(j_{\text{lc}}/j_{\min})$ while calculating the flux, and if $j_{\text{cir}}>j_{\text{lc}}^{2}/j_{\min}$ the rate is even slower than that (equation 15). This effect reduces the theoretical loss cone filling rate of the Milky way by a factor of $\approx 2$ for the case of HVS where only two-body relaxation in spherical symmetry is take into account.

Second, empty loss cone is much less empty at depth $\gg j_{d}$ than expected from any kind of diffusive model. According to diffusive approximations, the profile inside the loss cone decays exponentially over a scale of $\sim j_{d}$ (similar to $j_{2}$ of Cohn & Kulsrud (1978)). See our equations (9) or (14) or the classical work of Cohn & Kulsrud (1978). However, according to our new understanding, large scattering uniformly populate the loss cone. The profile deep inside the loss cone does not decay exponentially but is rather constant at a value of $\sim f(j_{\text{cir}})(j_{d}/j_{\text{lc}})^{2}\ln^{-2}(j_{\text{lc}/j_{\min}})$.

This second effect may have some observational significance. Let us discuss the nature of tidal events of giant stars in which the LC for tidal disruption event is empty . The hardness of tidal disruption events are characterized by the parameter:

[TABLE]

Previous analysis based on the diffusion approximation predicts the typical tidal event to have $\beta-1\sim j_{d}/j_{\text{lc}}$, were large $\beta$’s tidal events are exponentially rare with probability $\sim\exp\{-(\beta-1)j_{\text{lc}}/j_{d}\}$. Instead, we find that large scatterings populate the distribution deep in the LC, making it lower only by factor of $(j_{d}/j_{\text{lc}})\ln^{-1}(j_{\text{lc}}/j_{\min})$ compared to the occupancy distribution on the LC’s edge. Since the typical width of the edge is $\sim j_{d}$, the total number of particles on the edge is only $\sim\ln(j_{\text{lc}}/j_{\min})$ higher than the number of particles at high $\beta$’s. As a result, the typical fraction of high $\beta$’s tidal events out of all tidal events is $\sim\ln^{-1}(j_{\text{lc}}/j_{\min})$. Moreover, since the profile deep inside the LC is roughly flat, the probability of high $\beta$ events does not decays exponentially but scales as the surface within $j<j_{\text{lc}}\beta^{-1/2}$, i.e as $\beta^{-1}$.

This prediction may be verified observationally if deeply penetrating tidal events look different from shallow ones.

Acknowledgements

This research was partially supported by iCore, ISF and ISA grant

Appendix A Which scatterings are important ?

In section 3.3, we argue that scatterings $\gtrsim j-j_{\text{lc}}$ are less dominant outside the LC, and may be neglected. This assumption lead us to $j$’s dependency of the diffusion coefficient (equation 12). In this section we explain this assumption

We start with estimating the total contribution of scatterings of size $j_{\text{scat}}$ to the net influx toward angular momentums lower than $j_{\text{targ}}$ in the case $j_{\text{scat}}\ll j_{\text{targ}}-j_{\text{lc}}$ as well as in the case where $j_{\text{scat}}\gg j_{\text{targ}}-j_{\text{lc}}$. For simplicity, we assume $j_{\text{targ}}\gg j_{\text{lc}}$ ($j_{\text{targ}}\sim j_{\text{targ}}-j_{\text{lc}}$). However, same principles hold also for the case of $j_{\text{targ}}-j_{\text{lc}}\lesssim j_{\text{lc}}$.

For small scatterings, $j_{\text{scat}}\ll j_{\text{targ}}$, the diffusion approximation holds and therefore

[TABLE]

where $D_{j_{\text{targ}}}$ is the contribution to the diffusion coefficient of scatterings of size $\sim j_{\text{targ}}$ which satisfies $D_{j_{\text{targ}}}\sim R_{0}$.

For large scatterings, $j_{\text{scat}}\gg j_{\text{targ}}$, the diffusive approximation does not holds and the net influx have to be estimated from comparison of the influx and the outflux. The influx into $j<j_{\text{targ}}$ by scatterings of size $\sim j_{\text{scat}}$, is approximately the product of the rate of $j_{\text{scat}}$ scatterings ($\sim R_{0}j_{\text{scat}}^{-2}$), the number of particle at $j_{\text{scat}}$ surroundings ($j_{\text{scat}}^{2}f(j_{\text{scat}})$) and the probability that scattering of size $j_{\text{scat}}$ scatters an object into a sling-shot of size $j_{\text{scat}}\ll j_{\text{targ}}$ ($(j_{\text{targ}}/j_{\text{scat}})^{2}$). Therefore the influx into $j<j_{\text{targ}}$ by scatterings of size $\sim j_{\text{scat}}$ is

[TABLE]

Similar estimation of the outflow from angular momentums $<j_{\text{targ}}$ (without the probability term which in this case is unity) gives:

[TABLE]

Since $j_{\text{targ}}<j_{\text{scat}}$, $f(j_{\text{targ}})$ is necessarily $\leq f(j_{\text{scat}})$. Therefore $\gamma_{in}^{\text{large}}(j_{\text{scat}},j_{\text{targ}})$ is an upper limit for the net influx:

[TABLE]

In order to understand which scattering dominate the profile. we would like to compare for each $j_{\text{targ}}$, the total contribution of scatterings smaller than $j_{\text{targ}}$ ($\Gamma_{net}^{\text{small}}(j_{\text{targ}})$), with that of scatterings larger than $j_{\text{targ}}$ ($\Gamma_{net}^{\text{large}}(j_{\text{scat}})$):

[TABLE]

At this stage we are ready to formulate properly the question of "which scattering dominate the profile?". There are two possibilities:

• •

If $f(j)$ grows slower than $\propto j^{2}$ in $j_{\text{targ}}$’s neighborhood, than $\max\{\gamma_{net}^{\text{large}}\}=\gamma_{net}^{\text{large}}(j_{\text{targ}}\sim j_{\text{scat}})$ In this case all scatterings at sizes $<j_{\text{targ}}$ has equal contribution, while larger scatterings has smaller contribution.

• •

On the other hand, if $f(j)$ grows faster than $\propto j^{2}$ in $j_{\text{targ}}$’s neighborhood, than $\max\{\gamma_{net}^{\text{large}}\}>\gamma_{net}^{\text{large}}(j_{\text{targ}}\sim j_{\text{scat}})\sim\Gamma_{net}^{\text{small}}(j_{\text{targ}})$. In this case there is some scale$>j_{\text{targ}}$, which dominate the scatterings everywhere.

Understanding which scattering dominate the profile is therefore not simple, since it requires some assumptions about $f(j)$. We therefore discuss this question under the assumption that small scatterings dominate the profile, and afterwards under the assumption that large scatterings are dominant. We show that the assumption of small scattering dominance is self consistent and give arise to a steady state profile, while assumption of large scattering dominance, is not.

If small scattering ($j_{\text{scat}}\ll j_{\text{targ}}$) dominate the profile, it is set by diffusion with diffusion coefficient of the form 12. In this case the discussion of §3.3 is relevant and $f(j)$ has a steady state form of $f_{\text{diff}}(j)$ (equation 14). Since $f_{\text{diff}}(j)$ grows slower than $\propto j^{2}$, small scatterings are dominant and the argument is self-consistent.

On the other hand, if large scattering dominate the profile, the profile has to grow faster than $\propto j^{2}$ at $j_{\text{targ}}$’s neighborhood. Lets assume for instance that it grows as $f(j)=f(j_{\text{targ}})(j/j_{\text{targ}})^{N}$ ($N>2$) up to some angular momentum $j_{\gamma}\gg j_{\text{targ}}$, for which $j^{-2}f(j)$ is maximal. Under this assumption, the influx rate into $j<j_{\text{targ}}$ is: $R_{0}f(j_{\text{targ}})(j_{\gamma}/j_{\text{targ}})^{N-2}$. This dependency of the influx rate upon $j_{\text{targ}}$ prevent zero divergence and as a result, the profile can not be in steady state.

Steady-state profile therefore grantee dominance of the small scatterings outside the LC: In some particular $j_{\text{targ}}$, all scatterings of size $\lesssim j_{\text{targ}}$ have the same contribution to the diffusion coefficient while scatterings $\gtrsim j_{\text{targ}}$ can be ignored.

Appendix B Accelerating the numerical scheme

As a working example, we consider a system with $j_{d}\sim 10^{-1}$, $j_{\min}\sim 10^{-2}$ (while $j_{\text{lc}}\equiv 1$), evolved enough so that $j_{\text{sim}}\sim 10$. Since we are interested in the occupancy distribution inside the loss cone, we require that the total number of particles in its inner half (with surface of $\pi 0.5^{2}\sim 1$) would be at list $10$ particles for sufficient statistical accuracy. This means a particle density of $\sim 10$ deep in the LC.

The largest angular momentum in the simulation, $j_{\max}$ has to be a few times larger than $j_{\text{sim}}$ and we take $j_{\max}=30$. The particle density at $j\sim j_{\max}$ is higher by factor of $\sim(j_{\text{lc}}/j_{d})^{2}\ln^{2}(j_{\text{lc}}/j_{\min})\sim 2.5\cdot 10^{3}$ compared to the density deep in the LC which was assumed to be $\sim 10$ and therefore, the total number of particles is $10^{8}$.

Relaxation has to take place over a scale of $j_{\text{sim}}\sim 10$ meaning that each particle experience $\sim(j_{\text{sim}}/j_{\min})^{2}\ln^{-1}(j_{\text{sim}}/j_{\min})\sim 10^{5}$ steps which, leads to a total number of $\sim 10^{13}$ computations. We therefore introduce two major improvements to accelerate the numerical scheme.

B.1 Focusing our efforts on the loss cone surrounding

Most of the scatterings are small and occurs far from the LC. In these regions, calculating the frequent scattering of size $j_{\min}$ is not necessary. We therefore replace those frequent scatterings with a few larger ones, which are still smaller than the typical angular momentum over which $f$ is changing. This is done as follows. We set the minimal size of the scatterings at angular momentum $j$ to be

[TABLE]

and replaced $R(j_{\delta})$, the rate of each scatterings to:

[TABLE]

with

[TABLE]

Those definitions reduced significantly the number of calculations on one hand, and on the other hand mimic well the real behavior: the total contribution of the small scatterings to the diffusion does not changed as well as the rate of large scatterings. Moreover, the scatterings’ rate close to the LC has not been changed at all. In the working example above, this reduces the number of computations by three orders of magnitude.

In order to implement $R_{\text{eff}}(j_{\delta})$ in the simulation, The time length attributed to each step has to be recalculated (it is no longer $t_{\min}$ since steps of size larger than $j_{s}$ are more rare than once every $t_{\min}$). We therefore attribute to each time step a period of $t_{s}(j_{s})$ where $t_{s}$ is set so that the rate of scatterings of size $>ej_{s}$ would not be influenced by the presented adaptation and we find $t_{s}$ as a function of $j_{s}$ to be:

[TABLE]

Notice that for $j_{s}=j_{\min}$ we get $t_{s}=t_{\min}$ as expected.

B.2 Using relaxed profile to account for more particles

After time $t_{\text{sim}}$, the profile is already relaxed up to angular momentums $\sim j_{\text{sim}}$. However, if after achieving this steady state we let the system evolve for a period of another $t_{\text{add}}\equiv(j_{\text{lc}}/j_{\text{sim}})^{2}t_{\text{sim}}$, all particles would change their location by typically angular momentum $\sim j_{\text{lc}}$. This means that there is no correlation between the two situations on scale of $\lesssim j_{\text{lc}}$ and we can use this result as an independent profile.

Repeating this process over and over again, we may increase our statistics on the same physical condition without creating dependence between the different repeats. This may go on and on as long as $N_{\text{rep}}\cdot t_{\text{add}}\lesssim t_{sim}$ where $N_{\text{rep}}$ is the number of repeats.