Exploring the origin of stars on bound and unbound orbits causing tidal disruption events

Tidal disruption events (TDEs) provide a clue to the properties of a central supermassive black hole (SMBH) and an accretion disk around it, and to the stellar density and velocity distributions in the nuclear star cluster surrounding the SMBH. Deviations of TDE light curves from the standard occurring at a parabolic encounter with the SMBH depends on whether the stellar orbit is hyperbolic or eccentric (Hayasaki et al. 2018) and the penetration factor ($\beta$, tidal disruption radius to orbital pericenter ratio). We study the orbital parameters of bound and unbound stars being tidally disrupted by comparison of direct $N$-body simulation data with an analytical model. Starting from the classical steady-state Fokker-Planck model of Cohn&Kulsrud (1978), we develop an analytical model of the number density distribution of those stars as a function of orbital eccentricity ($e$) and $\beta$. To do so fittings of the density and velocity distribution of the nuclear star cluster and of the energy distribution of tidally disrupted stars are required and obtained from $N$-body data. We confirm that most of the stars causing TDEs in a spherical nuclear star cluster originate from the full loss-cone region of phase space, derive analytical boundaries in eccentricity-$\beta$ space, and find them confirmed by $N$-body data. Since our limiting eccentricities are much smaller than critical eccentricities for full accretion or full escape of stellar debris, we conclude that those stars are only very marginally eccentric or hyperbolic, close to parabolic.


INTRODUCTION
Most galaxies harbor supermassive black holes (SMBHs) with millions to billions of solar masses at their center.Tidal disruption events (TDEs) provide a good probe to identify dormant SMBHs in inactive galaxies.A star is tidally disrupted by an SMBH when the star approaches the SMBH closely enough that the black hole's tidal force exceeds the stellar self-gravity (Hills 1975).In classical TDE theory, a star on a parabolic orbit and is tidally disrupted by the SMBH at tidal disruption radius, r t = (M BH /m * ) 1/3 r * , where M BH , m * , and r * are the black hole mass, stellar mass and radius, respectively.Subsequently, half of the stellar debris falls back to the SMBH at a rate of t −5/3 so that the bolometric luminosity is proportional to t −5/3 if the mass fallback rate equals the mass accretion rate (Rees 1988;Evans & Kochanek 1989;Phinney 1989).However, recent observations have revealed that some observed TDEs show light curves, which deviate from the t −5/3 decay rate (Gezari et al. 2012;Holoien et al. 2014;Gezari et al. 2015;Miller et al. 2015;Holoien et al. 2016;van Velzen et al. 2019).Dozens of X-ray TDEs have light curves shallower than t −5/3 (Auchettl et al. 2017), while many optical/UV TDEs are well fit by t −5/3 (e.g.Hung et al. (2017)).
Some possible reasons for the deviation of the light curve from the t −5/3 law are discussed in current literature.The following are the three main reasons among them.First, the fallback debris would cause a self-crossing shock by a relativistic apsidal precession (Shiokawa et al. 2015;Piran et al. 2015;Ryu et al. 2020), outflowing a significant fraction of the debris Lu & Bonnerot (2020).Moreover, the secondary shock due to subsequently occurring collision forms an accretion disk.The bolometric luminosity at the photosphere clearly deviates from the t −5/3 decay (Bonnerot & Lu 2020).Second, even though the mass fallback rate follows the t −5/3 law, the radiative fluxes emitted from the accretion disk or disk wind modify the light curve variation.Lodato & Rossi (2011) have shown that the luminosity of the accretion disk observed in different bands may decay with diverse power-law indexes.For TDEs with well-evolved optically thick accretion disks, the observed X-ray light curves should decay following the form of a power-law multiplied by an exponential, which is caused by the Wien tail of the disk spectrum (Mummery & Balbus 2020).Moreover, when mass falls back to the SMBH at a super-Eddington rate, an outflow could be launched from the accretion disk.The luminosity of the outflow can decay with time shallower than t −5/3 (Strubbe & Quataert 2009).Final, the mass fallback rate can deviate from the standard t −5/3 rate due to external and internal properties of the tidally disrupted star, such as: its orbital eccentricity (Hayasaki et al. 2013(Hayasaki et al. , 2018;;Park & Hayasaki 2020); its orbital energy and angular momentum, the combination of these two quantities defines how deep its orbit reaches inside the tidal radius (we define a penetration factor β = r t /r p , where r p is the pericenter of the stellar orbit around the black hole); the detailed stellar internal structure and the possible survival of a stellar core during a partial tidal disruption (Guillochon & Ramirez-Ruiz 2013).
In this paper, we focus on the number density distribution of the stars, which cause TDEs, as a function of orbital eccentricity and the penetration factor.Since these orbital parameters leave imprints on the observable flux, the e and β could be obtained by fitting the light curves of the observed TDE with, e.g., MOSFiT, Guillochon et al. (2018).We can presume the number density by using this observed e and β in the end.Currently about a dozen of TDEs are known with β measurements (Mockler et al. 2019;Nicholl et al. 2019;Gomez et al. 2020), but they do not get the eccentricity independently, because the MOSFiT software used is based on hydrodynamic simulations of Guillochon & Ramirez-Ruiz (2013), who only simulated the e = 1 case.We note from this analysis that all the measured values of β are close to unity.Stone & Metzger (2016) suggested that the β value of the stars originating from the empty loss-cone regime should be very close to unity, while for the full loss-cone regime β could take larger values and the number density is proportional to β −2 (the definitions of empty and full loss-cone are given in Section 2.1).As we will discuss in Section 3.3, also in the full loss-cone regime many orbits may have β ∼ 1; larger β values occur in this case, but with a smaller probability.The eccentricity could provide further constraints on this issue, but as indicated before the eccentricity was set to exactly unity for all measurements so far.
In a preceding paper (Hayasaki et al. 2018) we have examined the distribution of tidally disrupted stars on the e-β plane by using N -body experiments of spherical nuclear star clusters.We found two interesting results: first, the eccentricities of the stars that causes TDEs usually take values between two critical eccentricities proposed by Park & Hayasaki (2020), which depend on β, and so there is some correlation between e and β.Second, the distributions of e and β vary with the mass ratio between stars and the central SMBH, as do the critical eccentricities.This is important for future studies using stars of different masses -here and in Hayasaki et al. (2018) we just limit ourselves to stars of equal mass.Light curve characteristics of tidal disruption flares depend on the eccentricity and its critical limits.This raises further interest in the distribution of orbital parameters.
In this paper, we analytically derive the number densities of bound and unbound stars that undergo TDEs in a spherical nuclear star cluster and test them by comparison with N -body simulations.We also examine the distribution of stars on the e-β plane by estimating the allowed eccentricity range for a given β.Predicting the relative frequency of TDEs with different eccentricity e and penetration factor β should help identifying realistic values of e and β for future hydrodynamic simulations of TDEs and also the interpretation of TDE observations (constraining the dynamic processes operating in the host cluster).We construct our analytical models in Section 2 and then describe the details of N -body simulations and compare the analytical number densities to the simulation results in Section 3. We discuss our results in Section 4 and draw our conclusions in Section 5.

ORBITAL PARAMETER DEPENDENCE OF THE STELLAR DISTRIBUTION
In this section we first briefly review the theory of the steady-state stellar distribution around a central black hole in a star cluster of Cohn & Kulsrud (1978) (hereafter CK78).It originates from a numerical solution of the orbit-averaged Fokker-Planck equation in energy-angular momentum space.The CK78 solution describes the distribution of stars inside the cusp surrounding the central SMBH, assuming the gravitational potential is dominated by the SMBH.It is well suited as a starting point to derive the stellar distribution in a phase space of orbital eccentricity e and penetration factor β (see Section 2.2).These quantities are of interest here, because they provide relevant input parameters for hydrodynamic simulations of TDEs, which in turn provide key information about the nature of TDE light curves.For our analysis of star cluster simulations with TDEs it is an advantage to use parameters closely related to the TDE and its observational characteristics, rather than the more conventional integrals of stellar motion.
In Section 2.3 we discuss the case of unbound stars experiencing TDEs.Earlier work by Magorrian & Tremaine (1999); Wang & Merritt (2004) and Stone & Metzger (2016) is based on a generalized treatment following CK78, using the Fokker-Planck equation also in the region of a galaxy unbound to the SMBH.For our purposes we choose a simpler but still useful approach in that regime.In what follows, we use subscript 'b' and 'u' to mark the quantities corresponding to bound and unbound cases, respectively.

Stellar distributions in energy-angular momentum space
The influence radius r h of an SMBH in the center of a star cluster is defined as the radius within which the enclosed stellar mass equals to the SMBH mass.Inside r h stars are considered as gravitationally bound to the black hole and a stellar density cusp forms.Following CK78 we characterize a stellar orbit in this region by both the specific orbital energy of a star E = v 2 /2 − GM BH /r and by its normalized squared angular momentum where v, G, J and J c are the velocity of the star, gravitation constant, specific angular momentum of the star and the corresponding circular angular momentum, respectively.In a spherically symmetric cusp with isotropic velocity dispersion, the stellar density distribution depends on the orbital energy only.1However, the assumption of isotropic velocity distribution breaks down because the SMBH removes stars with low angular momentum through TDEs.So, in classical loss-cone theory, the stellar number density n should also depend on R, and thus be a function of both E and R i.e. n = n(R, E) (Cohn & Kulsrud 1978).In phase space the loss-cone region is encompassed by R lc , where R lc is the square of the normalized loss-cone angular momentum (see equation A5).For the models described in this paper stars inside the loss-cone region can survive for no more than one orbital period, unless they find a way out before being disrupted (typically by being scattered out of the loss-cone by two-body relaxation).In reality partial tidal disruptions may occur (Zhong et al. 2022;MacLeod et al. 2013), stars may not be fully disrupted at the first passage near the tidal radius.In case of only full tidal disruptions considered in this paper the loss-cone runs out of stars quickly and n(R, E) vanishes to zero at R ≪ R lc .Simplified models based on moment equations of the Fokker-Planck equation (Amaro-Seoane & Spurzem 2001;Amaro-Seoane et al. 2004) assumed a sudden drop of n to zero at R lc , while the original work of CK78 shows the solution around and inside R lc follows a logarithmic profile and reaches zero at R 0 (defined by equation 6).Two-body relaxation encounters replenish the loss-cone by angular momentum diffusion.So, in steady state n(R, E) is determined by an equilibrium between disruption processes near the tidal radius and the replenishment process.By solution of the orbit averaged Fokker-Planck equation, taking both processes into account, CK78 found the following expression for the stellar density: where A(E) is an energy-dependent coefficient and R 0 is the square of the normalized angular momentum at the zero-boundary below which the number density goes to zero.The CK78 solution was limited to the Keplerian potential.Later works, such as Magorrian & Tremaine (1999) and Wang & Merritt (2004) that have taken the stellar potential into account also reported the logarithmic dependence on R.
Our focus is on the bound stars that cause TDEs (i.e.R 0 ≤ R ≤ R lc ), the cumulative number density, N TD,b (R, E), has the same ln R-dependence as equation 2 (because they originate from the stellar cusp described by the CK78 distribution) but the normalization coefficient is different, where the new coefficient is the number of bound stars that eventually enter the tidal radius with energy between E and E + dE.The quantity R and E in equation 3 shall take the values at the disruption, because these values are relevant to our theoretical models of the e and β distributions.The upper cutoff of R comes from the condition for tidal disruption: the separation between a star and the SMBH must be less than or equal to r t .This condition is translated to R ≤ R lc at the disruption according to the equations 1, A6 and A7.The value of N TD,b (E) results from the accumulation of TDEs with time, thus it is calculated as N TD,b (E) = F (E, t)dt, where F (E, t) is the flux of stars that enter the loss-cone at time t.In this work N TD,b (E) is obtained directly from the N -body simulation (for example using the orbital energy of the TDEs recorded in the simulation, see Figure 1), so we do not discuss F (E, t) here in detail-but see e.g.Section 2.2 in Merritt (2013).Note that the normalization of all number density distributions of tidally disrupted stars presented in this paper is to the total number of disrupted stars over the time of the simulation.From N TD,b (E) we obtain for the normalization coefficient By introducing where ∆R is the cumulative change of R over one orbital period of the star, Magorrian & Tremaine (1999) where is an approximation of the analytical solution derived by Merritt (2013), whose exact form is expressed in terms of Bessel series [also see equations (44,45) in Vasiliev & Merritt (2013)].Since R 0 is very close to R lc in the case of Q < 1, the number density almost goes to zero in the loss-cone region.Q < 1 is the empty loss-cone regime; while in the Q > 1 case R 0 ≪ R lc , we are in the full loss-cone regime.We will explain how to compute Q in Section 3.

Number density of bound stars
In this subsection we transform the number density of bound stars in the loss-cone according to equation 3 from the standard phase space variables R and E into new variables more suitable for our analysis of TDEs.We consider transformations into the following new pairs of variables: (R, E) into (e, E), (β, E) or (e, β).All Jacobian determinants are nonsingular (see Table 1), so all pairs can be used as new independent phase space variables.We focus in the following on (e, E) and (β, E)-after integration over E the resulting distributions in e and β can be compared with N -body simulations and also be used to analyze expected TDE characteristics (see Section 3).The details of the variable transformation are presented in Appendix A, which are derived in the Keplerian regime.For the marginally bound stars (E ≈ 0), the Keplerian assumption breaks down and the variable transformations given in the Appendix may become inaccurate.
Substituting the variable R in the expression of n TD,b (R, E) (equation 3) with 1 − e 2 (equation A4) and multiplying the corresponding Jacobian determinant provides the number density of the bound stars in the range of e ll ≤ e ≤ e ul : where are the boundaries of eccentricity obtained from the R = R lc limit and R = R 0 .Note that n TD,b (e, E) goes to zero outside of this eccentricity range.Substituting the variable R in the expression of N TD,b (R, E) with 2|E|/(|E t |β) (equation A8) and multiplying the corresponding Jacobian determinant, we obtain in the range of 1 ≤ β ≤ 1/g(Q), where we used equations (3) and (6) for the derivation.Note that the number density vanishes to zero at β = 1/g(Q), which corresponds to that R equals R 0 in the original number density (see equation 3).For Q ≫ 1, if ln β is negligible compared to − ln g(Q), equation 10 can be approximated as N TD,b (β, E) ≃ 2Q|E|A TD,b (E)/(|E t |β 2 ), verifying the β −2 dependence in the full loss-cone regime suggested by Stone & Metzger (2016).We also notice that for Q = 1 (g(Q) = 1/e; e: Euler's number), which corresponds to an energy value at the critical radius (Frank & Rees 1976;Amaro-Seoane et al. 2004), for fixed critical energy E = E crit we get the following results: (11)

Number density of unbound stars
After the pioneering work by Cohn & Kulsrud (1978); Shapiro & Marchant (1978) for bound stars there were more general papers, extending the domain of solution of the Fokker-Planck equation to the unbound stars in the galaxy.They focused on the tidal disruption event rate and studied the dependence of the event rate on the geometry of the host cluster (Magorrian & Tremaine 1999), the M BH − σ relation (Wang & Merritt 2004) and the stellar mass spectrum (Stone & Metzger 2016).All of them used the standard phase space variables of energy and angular momentum.The stellar number density in that case is written as n(E tot , R), with E tot = v 2 /2 − GM BH /r + Φ gal (r); the new term Φ gal (r) denotes the gravitational potential generated by all stars of the nuclear star cluster and the galaxy inside a radius r (here, for example, in the case of a spherical system).In the following we will argue that a direct variable transformation to our variables e and β as before is cumbersome and actually not really necessary.Let us first check the function R(β, E tot ), which defines the variable transformation from R and E tot to β and E tot .For β we need to find the pericenter distance r p as a function of E tot and R. From our definition of R in equation 1 we get from the expression for E tot above: r p is the smallest root of this equation in terms of r.Since that depends on the functional form of Φ gal (r) it is generally impossible to find an analytic solution; it has to be computed numerically for each galaxy.
In the case of eccentricity, another problem occurs-the definition of e for the two-body problem has no straightforward generalization for orbits in a more general star cluster or galactic potential.Typically, orbits in galactic potentials are not closed; generalized eccentricities may be defined using the angular momentum or pericenter and apocenter (r a ), e = (r a − r p )/(r a + r p ), but it is not always a conserved quantity except for in a spherical potential.For stars, which we are interested in, when they come close to the SMBH, the two-body eccentricity will be different from a value computed far out in the galaxy.Therefore, we look at the situation of a two-body problem only, for a hyperbolic encounter between a star and the SMBH.We compute e at a place near the tidal radius, and convert the orbital energy of the star to E = E tot − Φ gal (r t ), which is positive for an unbound star.Adopting the relation between e and β for the hyperbolic orbit (equation A9), we obtain There is no simple and universal relation between e, E tot and R, due to the complicated β(R, E tot ) term, resulting in a complex expression also for the transformation R(e, E tot ).Therefore, we do not use the number density of unbound stars in the form of Magorrian & Tremaine (1999), Wang & Merritt (2004) or Stone & Metzger (2016).Instead, we use a simpler approximation for the number density of unbound stars in terms of e and β near the tidal radius, which is also appropriate for the comparison with our N -body results (see Section 3).
Outside of the SMBH influence radius r h the loss-cone has negligible effect on the stellar distribution, because it is usually r t ≪ r h .Therefore the velocity distribution is close to a Gaussian along the principal axes of a velocity ellipsoid, also denoted as anisotropic Schwarzschild distribution; it allows for different velocity dispersions, e.g., in radial and tangential directions-for star clusters, see Amaro-Seoane et al. ( 2004), but see also Kazantzidis et al. (2004) for a counterexample in the galactic nuclei.In the following, we derive the N (β) based on a simple cross-sectional ansatz.The cross section for the stars that could pass within a distance r from the SMBH (with gravitational focusing) is Σ When the star's specific kinetic energy v 2 /2 at the orbital apocenter is much smaller than the gravitational potential GM BH /r at the pericenter of its orbit, as is the case discussed here, the cross section is approximately reduced to be Σ(r) ≈ 2πrGM BH /v 2 .Then the flow of stars that passes within r p can be estimated as f (< r p ) = nΣv ≈ 2πr p nGM BH /v, where n and v are the stellar number density and stellar velocity at the place from where these stars come.Thus we get the relation f (< r p ) ∝ (r p /r t )r t (see also equation 2 of Rees (1988), but note that we do not need to postulate isotropy here -it is sufficient to use the radial velocity only, because the tangential velocity is very small for loss-cone stars originating from a large distance to the black hole.Substituting β = r t /r p into the above relation, we find f hence the number density N (β) is proportional to β −2 and we write down the following expression: Then substituting the variable β in the expression of N TD,u (β, E) with E/[|E t |(e − 1)] (equation A9) and multiplying the corresponding Jacobian determinant, we have

COMPARISON WITH NUMERICAL EXPERIMENTS
Here we compare the distribution of tidally accreted stars in terms of e and β, which we have analytically derived in the preceding section, with N -body simulations.To get better statistical quality of the results we use in these comparisons the dependence on e and β only, rather than the joint 2D distribution in e, β (or a 3D distribution in E, e, β), for the reason of statistical noise, by integrating over all energies as follows: where we use the number densities of equations ( 8) and ( 10).For the number densities of unbound stars, we obtain N TD,u (e) and N TD,u (β) in the same way, but integrate from E = 0 to infinity.To evaluate the number densities from these equations we need to get an evaluation for Q.Equation 5shows that it requires the computation of the average angular momentum change per orbit ⟨∆J 2 ⟩.
To measure ⟨∆J 2 ⟩ from the simulation, one needs to record positions and velocities of all particles at very high frequency.We did not save such data from our models, but note that Vasiliev & Merritt (2013) have done such measurements and the results generally agree with the theoretical prediction, but have large scatters (see their Figure 8).Hence, we turn to the analytical solution of ⟨∆J 2 ⟩ to construct our theoretical model.⟨∆J 2 ⟩ is computed at the apocenter of a stellar orbit, because two-body relaxation affects the orbit most strongly at the apocenter passage.This assumption can be justified because the orbiting star passes its apocenter so slowly that it has more time to interact with the surrounding stars, and also the perturbing forces may exert a non-negligible torque to the passing star (Touma & Tremaine 1997;Zhong et al. 2015).From relaxation theory (Frank & Rees 1976;Merritt 2013), we get where J c (r a ) = √ GM r a is the specific angular momentum of a circular orbit at r a , t dyn = r a /σ(r a ) is the dynamical timescale, is the local relaxation time (Spitzer 1987), σ(r) is the velocity dispersion of the stars, ρ(r) is the radial density profile of the star cluster, m = M c /N is the mass of the star, and ln Λ = ln(0.11N) is the Coulomb logarithm (Giersz & Spurzem 1994).Note that equation 17 is correct only qualitatively, as there are some cases where the assumption used in equation 17 is invalid; e.g., in the ultrasteep stellar cusps (Fragione & Sari 2018;Stone et al. 2018), but such cusps are not presented in our simulations.
We introduce a dimensionless parameter, Q boost , of order unity to provide an approximate evaluation of Q. From equations ( 5) and ( 17), Q is then given by Here we have used R lc = J 2 lc /J c (r a ) 2 .One can approximate r a = a(1 + e) ≈ 2a for bound stars on highly eccentric orbits in the Keplerian potential, with the orbit's semimajor axis a.Since a = −GM BH /(2E), we conclude that for the given density profile and velocity dispersion, Q and g(Q) become a function of only E because of equation ( 7).However, the quantity −GM BH /(2E) diverges as E approaches 0. In practice, we compute the exact value of r a from the combined gravitational potential ϕ(r).
Note that for unbound tidally accreted stars, according to equations 14 and 15, the number density does not depend on Q.

Basic model
For comparison of N -body simulations with the analytical results, we use the data of our previously published study (Hayasaki et al. 2018).We choose from that paper two models; each has particle number N = 512K and r t = 10 −5 .They are the models with largest particle number and smallest tidal radius in that parameter study, we consider them as the ones closest to a real nuclear star cluster (though still not sufficient in terms of particle number).The two models differ only by their black hole mass: one has M BH = 0.01, while the other has M BH = 0.05 (referring to models 5 and 10 of Hayasaki et al. (2018), respectively), in units where the total cluster mass is unity.
Our spherical star cluster with N equal-mass stars and a star-accreting SMBH fixed at the center was initialized in the same way as in our previous papers (Hayasaki et al. 2018;Zhong et al. 2014)initially a Plummer model was used, which has a central flat core, which adjusts to the gravity of the central back hole during a few dynamical orbits, producing a cusp-like initial density distribution.More details about the time evolution and the profiles of density and velocity dispersion can be found in Zhong et al. (2014).We use dimensionless Hénon units, in which G = M c = 1 and the total energy of the system is E = −1/4 (Heggie 2014a,b).In the simulations, we take r t as a fixed accretion radius, in which all the stars are regarded as being tidally disrupted and removed from the simulations.More details can also be found in Hayasaki et al. (2018).Our N -body models adopt initially isotropic velocity distribution, therefore some stars are placed inside the loss-cone at the beginning.These stars shall cause a burst of TDEs.However, such a burst cannot last for a long time, because the loss-cone runs out of stars within one orbital period (at most, a few N -body time units), then the system enters the angular-momentum-diffusion-dominated phase.As a result, such initial surge of TDEs only accounts for less than a few percent of all the TDEs, hence, their impact on the validation of our theoretical models are negligible.After the initial adjustment, a central density cusp is established in the N -body simulation, even though the total simulation has been only about one-third of a half-mass relaxation time.Although the simulation times of these two models were less than one-third of the half-mass relaxation time, Preto et al. (2004) have shown that this is enough for the system to achieve the CK78 distribution.
Before approaching our final goal, to compare the number densities according to equations 16 with N -body simulations, we will first check the quantities Q = Q(E), A TD,b (E) and A TD,u (E) because they are required for the calculation of the analytical number densities.According to the definition, A TD,b (E) depends on N TD,b (E) and Q (equation 4) , while A TD,u (E) just equals to N TD,u (E).In order to evaluate these quantities, we measure N TD,b (E) and N TD,u (E) and approximate it by double-powerlaw function.Figure 1 illustrates the results obtained from the N -body model with M BH = 0.01, measured at the end of the simulation.In both panels, the red histograms shows N -body data of bound and unbound tidally disrupted stars, as a function of their energy.The stars are distributed between E min (< 0) and E max (> 0).Note that E min (roughly −2 in model unit) is much larger than E t (−500 in model unit).This is consistent with the loss-cone theory that the stars are originating far from the tidal disruption radius.N TD,b (E) and N TD,u (E) is used to compute the analytical expressions for A TD,b (E) and A TD,u (E).
To get an analytic expression for Q we use a similar method as before to approximate now the density profile and velocity dispersion profile of the star cluster using N -body data.To model the density profile, we use the following double-power-law function to fit the N -body data Then the 1D velocity dispersion is obtained via the Jeans equation (with G = 1), An example of the density and velocity dispersion profiles is depicted in Figure 2, which are measured from the simulation data when the density profile is stabilized.We also show the results of the double power law fitting on the density profile and the solution of Jeans equation for the velocity dispersion, which smooth the fluctuations in the data and are used to calculate Q.We observe that in the central part, well inside the influence radius (r ≤ 0.03), the simulated density profile show a steeper cusp (although very noisy due to low particle number) than the prediction of double-powerlaw fitting.This deviation only mildly affects our modeling, since the stellar mass in this cusp is less than 10 −3 and almost none of the disrupted stars are originated from this region.In our N -body model, the influence radius is defined as the radius within which the enclosed stellar mass equals the SMBH mass.The influence radius in the model with M BH = 0.01 (M BH = 0.05) roughly equals 0.1 (0.2).
From the fitted density profile (equation 20), we also compute the composited gravitational potential in the star cluster ϕ(r) and the apocenter r a for the (zero angular momentum) radial orbit with a given orbital energy E.An example of r a (E) is shown in the right panel of Figure 2.
Figure 3 depicts Q as a function of r a , which is evaluated by equation 19.Here, the density and velocity dispersion profile are modeled by equations 20 and 21 (see also Figure 2).For comparison purposes, we adopt two different values of Q boost = 1 and Q boost = 5.From the figure, we find that both the empty loss-cone regime (Q < 1) and full loss-cone regime (Q > 1) are present in our N -body data.The above criteria for the empty and full loss-cone regime are obtained by comparing the size of the loss-cone and the size of angular momentum diffusion (see equation 5).Magorrian & Tremaine (1999) have proposed another criterion, where the loss-cone regimes are separated by Q ≃ − ln R lc , which comes from the consideration of the loss-cone flux.This alternative criterion does not change our conclusion that both the empty and full loss-cone regimes exist in our N -body models.20), but replacing the density and radius quantities with number density and energy, respectively.

Distribution of tidally accreted stars in eccentricity and penetration factor
In the previous subsection, we have derived theoretical, analytical expressions for the distribution of tidally accreted stars (bound and unbound ones) in terms of eccentricity e and penetration factor β; in order to achieve that, we have used double power law functions for the stellar density, the Jeans equation for the stellar velocity dispersion, and double-power-law function for the energy distribution of tidally disrupted stars.Now we will check the final results of the previous subsection for N TD,b,u (e) and N TD,b,u (β) (see equation 16) directly against the N -body data of particles arriving at the tidal radius.Figure 4 shows the dependence of the number density of bound stars on the orbital eccentricity.The red histogram represents the simulated number densities.The uncertainties of the measurements are computed based on the Poisson error, the 1σ confidence level single-sided upper and lower limits are computed with equations ( 9) and ( 12) in Gehrels (1986).The black and cyan curves represent the theoretical number densities obtained with different Q boost values (note this specific color setting for Q boost is used in Figures 3, 4, 5 and 8).We find that Q boost can mildly affect the theoretical N TD,b (e).In both panels, the distributions are quite narrow near the parabolic case (e = 1).The simulated number densities are also in good agreement with the theoretical ones, except for some stronger fluctuations around e ∼ 0.996 and e ∼ 0.998.This is because the particle resolution of our N -body simulations is not sufficient there.The number density of the M BH = 0.01 case is wider for the orbital eccentricity than that of the M BH = 0.05 case.This trend can be interpreted as follows: from equation ( 9), the lowest eccentricity of the bound stars can be estimated to be e ll (M BH , E min ) = 1 − 4r t |E min |/GM BH .We find that E min = −2 in the M BH = 0.01 case, whereas E min = −4 (ignoring the isolated bins) in the M BH = 0.05 case.Substituting each quantity into the above equation, we find e ll (0.01, −2) = 0.9960 and e ll (0.05, −4) = 0.9984.These evaluations are consistent with the number density distributions shown in Figure 4.
Figure 5 shows the dependence of the number density of bound stars on the penetration factor β. The figure format of the two panels is the same as for Figure 4. We find Q boost has strong effect on the theoretical N TD,b (β).For bound stars, the maximum value of Q is achieved at r a = r h .Figure 3 shows that in the M BH = 0.01 case, Q(r h ) = (1.89,9.45) for Q boost = (1, 5), respectively.Since N TD,b (β, E) = 0 when β > 1/g(Q), the integrated N TD,b (β) shall vanish beyond β = 7.4 in the Q boost = 1 case, while for the other Q boost case, the vanishing point extends to much higher β, as shown in the left panel of Figure 5 [the vanishing behavior of the theoretical models in the right panel (M BH = 0.05) can be understood in the same way].We find that the analytical number densities obtained with Q boost = 5 are in good agreement with the simulated ones within the range of β ≲ 10.For β ≳ 10, the deviation between the analytical and simulated number densities gets larger, because of the poor numerical resolution of the N-body models.
Figures 6 and 7 compare the analytical number densities of N TD,u (e) and N TD,u (β) with those of the simulated number densities.Our theoretical predictions match well with the simulated number densities in both figures.As in the bound star case, the number densities of M BH = 0.01 case are Dependence of Q on the apocenter distance radius r a for bound stars.Hénon units are adopted for the plot.The solid ocher line and horizontal black dotted line delineate the − ln R lc curve ( see Magorrian & Tremaine 1999 for details) and the Q = 1 line, respectively, which provides the criterion for distinguishing the empty and full loss-cone regimes.The vertical dotted red line represents the influence radius of the SMBH.The solid black and light blue lines indicate the Q curves with Q boost = 1 and Q boost = 5, respectively, which are given by applying the fitted simulation results to equation ( 19). more widely distributed over the eccentricity than the M BH = 0.05 case.Substituting E max and |E t | = GM bh /r t into equation (A9) with β = 1, we obtain e max (M BH , E max ) = 1 + 2r t E max /(GM BH ).We find from the simulated value of E max that e max (0.01, 1.259) = 1.0025 and e max (0.05, 1.585) = 1.0006.These suggest that the number density is more widely distributed over the orbital eccentricity in the star cluster with the less massive black hole.

Distribution of stars on the eccentricity-penetration factor plane
At the time of tidal disruption, it is the eccentricity e and the penetration factor β that can be related to the observational characteristics (e.g., light curve).Also in our N -body simulations, we have a direct handle to determine these two quantities for any tidal disruption locally, without knowing anything about the large-scale distribution of stars and the gravitational potential.Therefore, we check here what we can deduce from our previous analytical results for the distribution of tidally disrupted stars on the e-β plane and compare again the expectations with the simulation data.For a given energy E = −GM BH /2a at the R = R 0 limit, we can define the minimum orbital eccentricity  of a bound star by using r t /β = a(1 − e) as where β = 1/g(Q) is obtained through equations ( 6) and (A8).Adopting β = 1, we find equation ( 22) corresponds to equation ( 9) at r t /a ≪ 1: e ll = 1 − 4r t |E|/(GM BH ) ≈ 1 − 2r t |E|/(GM BH ).The unbound stars all have E ≤ E max .We use again |E t | = GM bh /r t and equation (A9), to obtain For Q = 1, we obtain the orbital eccentricity on the boundary between the empty loss-cone and the full loss-cone regimes: where the corresponding semimajor axis, a lcb , can be obtained from Eq. 19 as where we have adopted r a = 2a, and r a is computed from the combined gravitational potential of stars and SMBH, see right panel of Figure 2 for the result.In this case, we find e lcb ≈ 1 for β ≥ 1. Figure 8 shows the distribution of bound and unbound stars that undergo TDEs on the e-β plane.The solid magenta curve denote e max .The e min (solid) and e lcb (dashed) lines are plotted with two different Q boost values that are indicated by the colors, respectively.The value of E max is obtained from our N -body simulations and a lcb ≈ r crit /2 is read out from Figure 3.We notice that Q boost = 1 gives poor estimate of e min : lots of data points are lying outside the boundary.On the other hand, the Q boost = 5 models show better agreement with the data.While the space between e lcb and e min corresponds to the empty loss-cone regime, the space between e lcb and e = 1 corresponds to the full loss-cone regime.The stars located above the e = 1 axis are supplied to the black hole from the Maxwellian distribution regime.We find that only a small fraction of the stars originate from the empty loss-cone regime and the corresponding values of β are distributed around unity, whereas most of the stars are originally supplied from the region between e lcb and e max over the much wider range of β.These results are consistent with loss-cone theory: due to the diffusion nature in the empty loss-cone regime, the loss-cone flux is much smaller than the full loss-cone regime, and the stars inside the empty loss-cone cannot penetrate the tidal radius too much.
We also find that there are some stars outside e min or e max .These outliers seem to violate the losscone theory.There are two possible reasons for the outliers behaving unexpectedly: first, the more energetic close two-body encounter occurs in the N -body simulations, leading to an enhancement of the angular momentum exchange, so that the stars have R < R 0 (Lin & Tremaine 1980).Second, the quantity Q depends on the density and the velocity dispersion of the star cluster, which can fluctuate with radius and evolve with time.The number of outliers is very small, so e min and e max are generally useful as limits.
The critical eccentricities are (Park & Hayasaki 2020): where k = 0 is the original case discussed in Hayasaki et al. (2018).They found by N -body experiments that stars on marginally eccentric and marginally hyperbolic orbits are the main source of TDEs in a spherical nuclear star cluster.
In our case, we have m ⋆ = 1/N ; for the simulations used in Figure 8 we find, for example, e crit,1 ≈ 0.885 and e crit,2 ≈ 1.12 for M bh = 0.01, β = 1, and k = 0 case and e crit,1 ≈ 0.933 and e crit,2 ≈ 1.067 for M bh = 0.05, β = 1, and k = 0 case.In any case, the maximum and minimum eccentricities predicted and found in the simulation, as shown in the figure, are much smaller than e crit,2 and much larger than e crit,1 , i.e. very close to the parabolic case.This means, in the terminology of Hayasaki et al. (2018) that all our TDEs are either marginally eccentric or marginally hyperbolic.This is not surprising, because our analysis is based on the same N -body data; however, we have derived in this subsection much narrower limits for e min and e max , which are based on the diffusion theory going back to CK78.

DISCUSSION
We have derived a semi-analytical model for the number distribution of eccentricity e and penetration factor β of tidally disrupted stars in a spherical nuclear star cluster around an SMBH.It has been compared to the results of our previously published (Hayasaki et al. 2018) direct N -body simulation with good agreement.To get our model, we use double-power-law functions to fit the stellar density of the N -body data and compute the velocity dispersion profile via the Jeans equation, as well as use double-power-law function to model the energy distribution of the bound and unbound tidally disrupted stars.Our method is based on the classical results of loss-cone diffusion for bound stars, using the Fokker-Planck equation (CK78).For unbound stars, we use a simple approximation based on the assumed Maxwellian character of the stellar distribution function.For an improved treatment of unbound stars, we need to consider the effects of the local self-gravity of stars and other external factors as the galactic potential.Our model is useful for discussing the scaling behavior of TDE statistics; current N -body simulations are still far away from realistic particle numbers and sizes of tidal radii (Hayasaki et al. 2018).Nevertheless, simulation data like the ones presented here have been used to extrapolate from our unphysically small particle numbers and large tidal disruption radii to real galactic nuclei -typical results for the TDE rate in N -body simulation models range around or little above 10 −6 per year per galaxy (Zhong et al. 2014;Panamarev et al. 2019;Li et al. 2023), which is the lower bound of Stone & Metzger (2016); the latter paper gives a range of up to 10 −4 per year per galaxy, in accord with recent work of Bortolas et al. (2023).When comparing such rates, one should bear in mind that the goal of our paper is not to give any accurate predictions of observed TDE rates.Cited papers include stars from a much wider range of origin, out into the bulge; our work has not yet properly accounted for a realistic mass spectrum and tidal disruption properties depending on stellar type and parameters (but this work is in progress).Another parameter not yet carefully checked in our models is the effect of the black hole mass (relative to the cluster mass and relative to the stellar particle mass).
Let us first qualitatively discuss how the number densities derived in our work would vary with the black hole mass relative to the star cluster.We assume a power-law density profile in the cusp ρ(r) = ρ h (r/r h ) −s , where ρ h = (3 − s)M BH /(4πr 3 h ) is the density at the influence radius according to the definition of r h .Note that the total stellar mass inside the influence radius is equal to the black hole mass.The velocity dispersion inside the influence radius follows σ 2 (r) ≈ GM BH /r.A critical radius, where the star consumption is balanced with its replenishment by two-body relaxation, is defined by the conditions θ lc = θ D (in the notation of Frank & Rees (1976) and Amaro-Seoane et al. ( 2004)) or Q = 1 (in our notation following CK78).Following Frank & Rees (1976) and Brockamp et al. (2011) we get As we have shown in the previous section, a Q boost factor is important for matching the theoretical model with the N -body results, so we also add the Q boost factor to equation 27 and in the following part we adopt Q boost = 5.Note that the left term yields for the traditional value s = 7/4 the result r crit ∝ (r t M BH ) (4/9) consistent with Baumgardt et al. (2004).The right-hand side delivers a different scaling, because we have used additionally results of the scaling procedure described in Zhong et al. (2014): BH , and we simply adopt 1/2 for the scaling of r h instead of 0.54, which is obtained from the M BH − σ relation, log(M BH /M ⊙ ) = 8.18 + 4.32 log[σ/(200km s −1 )] (Schulze & Gebhardt 2011).We get for the ratio of critical to influence radius using our scaling (ignoring the slowly varying logarithmic term).Therefore, for s < 4, the ratio of critical to influence radius increases with black hole mass, as was already observed for the standard case by Frank & Rees (1976).Figure 9 depicts r crit /r h as a function of M BH (see also equation 28) in the range of 10 3 M ⊙ ≤ M BH ≤ 10 8 M ⊙ , where we adopt s = 1.75 for the Bahcall-Wolf cusp (Bahcall & Wolf 1976) and s = 1 for the cusp obtained from the N -body simulations.
Assuming that the semimajor axis a min corresponding to E min is given by a fixed fraction of r crit as Since the ratio of r t to r crit is less than 3 × 10 −4 over the whole range of the black hole mass, e min is always larger: 0.994 for f = 0.05 and β = 1.Because of 1 − e min ∝ r t /r crit ∝ M (s−9)/[6(4−s)] BH , e min is closer to 1 for s < 4 as the black hole mass is larger.We also confirm that the black hole mass dependence of the pericenter radius is the same as that of the tidal disruption radius, i.e., r p,min = a min (1 − e min ) ∝ M Next, let us see how e min and e lcb (see equations 24 and 29) depend on the black hole mass on the e-β plane.The left panel of Figure 10 shows it for the M BH = 10 3 and 10 4 M ⊙ cases.The e lcb curve of M BH = 10 3 M ⊙ gets larger than M BH = 10 4 M ⊙ case (thick curve).In addition, as mentioned in Section 3.3 (see also Figure 8), most of the bound stars are distributed around the e lcb curve on the plane.These suggest that the star is distributed closer to e = 1 over the whole range of β as the black hole mass increases.The right panel of Fig. 10 depicts how the e min curve depends on the black hole mass on the e-β plane.It is clear from the panel that e min gets closer to 1 as the black hole mass is larger.This is consistent with 1 − e min ∝ M (s−9)/[6(4−s)] BH (s < 4) as we estimated in the previous paragraph.In summary, Figure 10 suggests that the stars are supplied into the black hole on extremely marginally eccentric to parabolic orbits for a spherical cluster with M BH > 10 7 M ⊙ black hole.
Our analysis of the e-β distributions could provide a good tool to probe the dynamical status of the stars causing TDEs in a star cluster.Hayasaki et al. (2013Hayasaki et al. ( , 2016) ) and Bonnerot et al. (2016) studied the accretion disk formation by performing hydrodynamic simulations, where the authors have adopted (β, e) = (5, 0.8) as initial values, although they also used other combinations of (e, β).This parameter set is clearly ruled out in our model.However, this does not mean that such a very tightly bound TDE cannot occur.The tightly bound stars are likely to supply to the loss-cone not by two-body encounters, but by other mechanisms: the tidal separation of stellar or compact binaries approaching the SMBH (Fragione & Sari 2018), accretion-disk-mediated TDEs (Kennedy et al. 2016), TDEs produced by a recoiling SMBH (Gualandris & Merritt 2008;Li et al. 2012) or by a merging SMBH binary (Hayasaki & Loeb 2016;Li et al. 2017).Finally, we note that star clusters possessing radially biased velocity distribution, especially at high binding energy, or having nonspherical gravitational potential may help to increase the rate of TDEs with e < 1 (not too close to 1) and extend the β distribution to β > 1 in the empty loss-cone regime; see also the detailed discussion at the end of Section 5.

CONCLUSION
Understanding TDEs and their light curves provides a clue to the properties of the central SMBH, the accretion disk around it, and to the stellar density and velocity distributions in the nuclear star cluster surrounding the SMBH.The link between TDEs and central star clusters in galactic nuclei has been a classical subject of seminal papers, such as Frank & Rees (1976); Bahcall & Wolf (1976) finding the classical density distribution near central SMBHs.Dokuchaev & Ozernoi (1977a,b) first noted that accretion and the tidal disruption of stars with low angular momentum cause the density profile to flatten out toward the black hole (the energy distribution function f (E) drops towards zero as they showed; today we would call this the empty loss-cone region-see also Ozernoi & Reinhardt (1978) for a summary of the topic at the time).CK78 put this on a more quantitative footing by using the technique of solving the orbit-averaged Fokker-Planck equation.Among the classical work in this field, also Rees's conjecture about the fate of tidal debris (Rees 1988) is most noteworthy; it has been expanded more recently by Hayasaki et al. (2018) and Park & Hayasaki (2020), looking for critical eccentricities that separate partial from full mass loss (hyperbolic) and partial from full mass accretion (eccentric).
Our study generalizes and expands this by computing approximate distributions of bound and unbound stars by predicting analytically (and comparing with N -body data) the number densities as a function of eccentricity e and penetration factor β, both of which are key parameters for the prediction of the observational appearance of TDEs.
By following the generalized model of CK78 for bound and unbound stars, we predict that tidally disrupted stars, for all penetration factors, occupy only a small range in eccentricities, much smaller than the critical eccentricities cited above.We estimate some minimum, maximum, and typical eccentricities for tidally disrupted stars.They are all very close to the parabolic case, either very marginally eccentric or very marginally hyperbolic (see Fig. 8).
Amaro-Seoane & Spurzem (2001) and Amaro-Seoane et al. ( 2004) proposed another model for the loss-cone in a spherical star cluster.It is interesting to note that they have also derived the density and velocity dispersion of bound and unbound loss-cone stars by using moment equations of the basic Fokker-Planck equation and very similar principles to CK78 and this paper.Due to the use of moment equations (the so-called gaseous model of star clusters; see, e.g., Giersz & Spurzem 1994) their analysis is completely based on density and velocity profiles rather than orbits with energy and angular momentum (or eccentricity and penetration factor).Their model takes into account an anisotropic velocity distribution also for the unbound stars.In the future, a more quantitative comparison of the two models could be done.Our theoretical predictions have all been tested against the data of our previously published direct N -body simulations (Hayasaki et al. 2018); comparison with our analytical model helps to understand the scaling behavior of the N -body simulations, since we can still not yet do realistic particle numbers for them.Our primary conclusions are summarized as follows: 1. Our results provide the number density of bound tidally disrupted stars as a function of orbital eccentricity e and penetration factor β; in practice, we used the cumulative numbers of disrupted stars N TD (e) and N TD (β) over some simulated time, since they can be directly compared with simulation results.To get them, we fit the stellar density with double power law profile and solved the corresponding velocity dispersion via the Jeans equation.We also use double-powerlaw functions to model the energy distribution of tidally disrupted stars (obtained from the N -body data).
2. From these results, we have analytically derived three characteristic orbital eccentricities: e min , e max , and e lcb in the loss-cone region, where e min and e max take the minimum and maximum values for a given β, respectively, whereas e lcb represents the orbital eccentricity which gives the boundary between the empty and full loss-cone regimes.These eccentricities are given by equations ( 22), (23), and (24), respectively.We have confirmed that the stars causing TDEs are distributed between e min and e max on the e − β plane by N-body experiments.Moreover, we find most of the bound stars are focused between e lcb and e = 1, i.e., in the full loss-cone regime, whereas the remaining bound stars are originating from the empty loss-cone regime.This result is consistent with the loss-cone theory.
3. We conclude from the limiting eccentricity values that they are very close to the parabolic case, and far away from the critical eccentricities for complete debris accretion or complete debris escape from the SMBH.We have shown that this conclusion holds also for larger more realistic black hole masses.
Our model of angular momentum diffusion at a given energy value E, as given in Eq, 19, uses a free parameter Q boost for fitting to our simulation results.Merritt (2013) suggested using the steadystate solution of a Fokker-Planck equation in angular momentum space, for every energy value, in order to obtain the flux across the loss-cone boundary, i.e. a value of Q boost in our terminology.This concept is based on the assumption that steady state in angular momentum space is achieved much faster than in energy space; it is used by the PhaseFlow 1D Fokker-Planck code (Vasiliev 2017).In our paper, we prefer not to follow such a two-timescale approach, rather keep all our fitting procedures in energy space and use the free factor Q boost .There may be several factors which could affect the clean separation of angular momentum and energy diffusion time scales.Most notable are rotation, axisymmetric gravitational potentials-see our earlier 2D Fokker-Planck model with full 2D representation of angular momentum diffusion near the loss-cone in Fiestas & Spurzem (2010); Fiestas et al. (2012), based on the 2D Fokker-Planck code used by Einsel & Spurzem (1999), but also strong anisotropy could have an impact here Szölgyén et al. (2019).Arguably the use of PhaseFlow will be a good method to quickly get the n(e) and n(β) distributions for TDEs in galaxy models, using energy distributions, and directly aimed at real systems (Pfister et al. 2019(Pfister et al. , 2020;;Bortolas et al. 2023).In any case, from the knowledge of n(e) and n(β) one could estimate the distribution of peak mass fallback rate in different galaxies, since the peak mass fallback rate depends on both β (Guillochon & Ramirez-Ruiz 2013) and e (Hayasaki et al. 2013;Park & Hayasaki 2020).
Some issues remain to be subject of further work; in our N -body simulations, we have adopted a fixed position and mass of the central black hole.For a star cluster with an intermediate mass black hole (IMBH), for example, the ratio of black hole mass to stellar mass will be smaller than in this paper, and the Brownian motion of the IMBH will not be suppressed.The Brownian motion of the black hole modifies the energy distribution of stars so that the number density and e-β distributions can be significantly affected, which is subject of our future work.
Our model is a high-resolution direct N -body model of a nuclear star cluster, following individual stellar orbits and TDEs.We have measured the distribution of orbital parameters of tidally disrupted stars and compared the results with a semi-analytical model.Our N -body models are not restricted to spherical symmetry, even though in this paper we do study only spherical nuclear star clusters.We do not intend to predict detailed TDE rates for specific galaxies, such as other models based on 1D Fokker-Planck theory (see, e.g., Stone & Metzger 2016;Pfister et al. 2019Pfister et al. , 2020;;Bortolas et al. 2023).Rather, we are interested in the analysis of the stellar distribution, relaxation and accretion processes only in the inner zone of a nuclear star cluster (stellar mass limited to ten times the black hole mass).Models based on 1D Fokker-Planck theory are computationally much faster and can extend much farther out, but rely on approximations such as spherical symmetry and a steady state in angular momentum diffusion.
Future work in the domain of our N -body simulation model is to include a stellar mass spectrum, stellar populations with different ages, direct stellar collisions and relativistic dynamics of stellar mass black holes in the nuclear star cluster.This goes along with a more realistic treatment of tidal disruptions (partial and full; Zhong et al. 2022;MacLeod et al. 2013) as well as direct plunges and relativistic or dissipative inspirals (see, e.g., our recent paper Li et al. 2023).Furthermore the assumption of spherical symmetry will be relaxed in favour of rotating, axisymmetric (Fiestas & Spurzem 2010;Zhong et al. 2015) and triaxial models (Norman & Silk 1983;Poon & Merritt 2002, 2004;Merritt & Poon 2004).All of these will disturb the steady-state picture underlying our current paper and have interesting consequences for TDEs and produce also gravitational-wave events instead of TDEs.The Jacobian determinants corresponding to the above variable transformations are summarized in Table 1.Note-The Jacobian determinants take a positive sign.
2 In the simulation, we record the position r and velocity v at the time when a star enters r t .The two-body eccentricity of a unbound star is e = 2EJ 2 /(GM BH ) 2 + 1, where J = |r × v| and E = |v| 2 /2 − GM BH /r.Then β is computed with equation A9.

Figure 1 .
Figure1.Energy dependence of the cumulative number densities for both the bound (left panel) and unbound (right panel) disrupted stars with M BH = 0.01 cases.In each panel, the red histogram shows the simulated number density (dN/dE), whereas the gray line represents the double-power-law curve, which is fitted to the simulated number density using equation taking the same functional form of equation (20), but replacing the density and radius quantities with number density and energy, respectively.

Figure 2 .
Figure2.Radial profiles of the density (left) and of the square of the stellar velocity dispersion (middle) for the model cluster with M BH = 0.01.In the left panel, the red dots and gray solid line represent the simulated and numerically fitted density profiles, respectively.In the middle panel, the red dots and gray solid line denote the square of stellar velocity dispersions evaluated by the simulation and by equation (21), i.e., the Jeans equation, respectively.The right panel depicts the apocenter distance radius r a as a function of the specific orbital energy E (= v 2 /2−GM BH /r) for a radial orbit with zero angular momentum.The horizontal blue dotted line indicates the influence radius r h , whereas the vertical black dotted line corresponds to the line of E = 0.

Figure 4 .
Figure 4. Dependence of the number density of the bound stars on the orbital eccentricity.The red histogram represents the simulated number density, while the solid black and light blue lines delineate the theoretical number densities with Q boost = 1 and with Q boost = 5, respectively.Note that the theoretical number density is computed by equation (16).The left panel corresponds to the M BH = 0.01 case, while the right panel corresponds to the M BH = 0.05 case.The error bars indicate the statistical uncertainty corresponding to the standard deviation.

Figure 5 .
Figure 5. Dependence of the number density of the bound stars on the penetration factor (β).The figure formats are the same as Figure 4 but for β.

Figure 6 .Figure 7 .
Figure 6.Same format as Figure 4 but for the unbound stars.Equation 15 is used to evaluate the theoretical number density quantitatively.

Figure 8 .
Figure 8. Distribution of the stars, which can cause TDEs, on the e-β plane (left: M BH = 0.01; right: M BH = 0.05).The solid magenta curve denote e max .The e min (solid) and e lcb (dashed) are plotted with two different Q boost values that are indicated by the colors, respectively (see also equations 22-24).
where f is a parameter determined by N -body simulations, we have E min ∝ −M 7/3, E min decreases as M BH decreases.Substituting E min = −GM BH /(2f r crit ) into equation (22), we obtain

Figure 9 .
Figure 9.Black hole mass dependence of r crit /r h .The two slopes for the density profiles: s = 1.75 (Bahcall-Wolf cusp) and s = 1 (N -body simulations) are adopted.
Figure10.Black hole mass dependence of two characteristic eccentricities: e min and e lcb in the e-β plane.We adopt the Bahcall-Wolf density cusp and Q boost = 5 for all the models.The left panel shows the β-dependence of e min and e lcb for the M BH = 10 3 M ⊙ and 10 4 M ⊙ cases.The dashed lines denote the eccentricity, e lcb , between the empty and full loss-cone regimes (see equation 24), whereas the solid lines denote the eccentricity, e min , of most tightly bound stars (see equation 29).The different color shows the different black hole mass.In the right panel, all the lines represent e min .We adopt different line styles for the different black hole masses for 10 3 M ⊙ ≤ M BH ≤ 10 8 M ⊙ .
S.Z. has been supported by the National Natural Science Foundation of China (NSFC 11603067) and acknowledges the support from Yunnan Astronomical Observatories, Chinese Academic of Sciences.K.H. has been supported by the Korea Astronomy and Space Science Institute (KASI) under the R&D program supervised by the Ministry of Science, ICT and Future Planning, and by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2016R1A5A1013277, 2017R1D1A1B03028580, and 2020R1A2C1007219 (K.H.)).K.H. has been supported by the National Supercomputing Center with supercomputing resources including technical support (KSC-2019-CRE-0082 (K.H.)) and in part by the National Science Foundation under Grant No. NSF PHY-1748958.K.H. has been also financially supported during the research year of Chungbuk National University in 2021.The authors acknowledge the Yukawa Institute for Theoretical Physics (YITP) at Kyoto University.Discussions during the YITP workshop YITP-T-19-07 on International Molecule-type Workshop "Tidal Disruption Events: General Relativistic Transients" were useful to complete this work.The authors also acknowledge support by the Chinese Academy of Sciences (CAS) through the Silk Road Project at NAOC.We are grateful for the support from the Sino-German Center (DFG/NSFC) under grant no.GZ1289.S.L., P.B. and R.S. acknowledge the Strategic Priority Research Program (Pilot B) Multi-wavelength gravitational wave universe of the Chinese Academy of Sciences (No. XDB23040100).S.L. and R.S. acknowledges Yunnan Academician Workstation of Wang Jingxiu (No. 202005AF150025).The work of PB was supported by the Volkswagen Foundation under the special stipend No. 9D154.PB thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023.

Table 1 .
Table of the variable transformation and the associated Jacobian determinant t |(e−1) 2 ]