Tidal disruption events, main-sequence extreme-mass ratio inspirals and binary star disruptions in galactic nuclei

TL;DR

This paper presents an analytical model of stellar dynamics near supermassive black holes, estimating rates of tidal disruptions and extreme-mass ratio inspirals, and explores how binary disruptions influence these rates in galactic nuclei.

Contribution

It introduces a simplified analytical framework for stellar scattering processes, linking TDE and MS-EMRI rates and examining the impact of binary disruptions on these rates.

Findings

MS-EMRI rate is about 1% of TDE rate in the Milky Way.

Binary disruptions can increase MS-EMRI rates to nearly TDE rates if binary fraction is high.

Stellar collisions significantly affect phase space near SMBHs.

Abstract

The Galactic Center has been under intense scrutiny in the recent years thanks to the unprecedented missions aiming at measuring the gas and star dynamics near the supermassive black hole (SMBH) and at finding gravitational wave (GW) signatures of inspiralling stellar black holes. In the crowded environment of galactic nuclei, the two-body interactions alter the distribution of stars on long timescales, making them drift in energy and angular momentum. We present a simplified analytical treatment of the scattering processes in galactic stellar nuclei, assuming all stars have the same mass. We have discussed how the interplay between two-body relaxation and gravitational wave emission modifies the slope of the inner stellar cusp within the SMBH sphere of influence, and calculated the rates of tidal disruption events (TDEs) and main-sequence extreme-mass ratio inspirals (MS-EMRIs) of stars that are tidally disrupted by the SMBH. We find that typically the ratio of the TDE and MS-EMRI rates is the square of the ratio of the tidal and Schwarzschild radii. For our Galaxy, this implies that the rate of MS-EMRIs is just about a percent of the TDE rate. We then consider the role of stars injected on highly eccentric orbits in the vicinity of the SMBH due to Hills binary disruption mechanism, and show that the MS-EMRI rate can almost approach the TDE rate if the binary fraction at the SMBH influence radius is close to unity. Finally, we discuss that physical stellar collisions affect a large area of phase space.

Equations

$$T_{2B}(r)=\frac{P(r)}{\ln \Lambda }{\left(\frac{M}{m}\right)}^2\frac{1}{N(r)},$$

$$T_{2B}^J(r,r_p)=\frac{P(r)}{\ln \Lambda }{\left(\frac{M}{m}\right)}^2\frac{1}{N(r)}\left(\frac{r_p}{r}\right).$$

$$T_{2B}^E(r,r_p)=\frac{P(r)}{\ln \Lambda }{\left(\frac{M}{m}\right)}^2\frac{1}{N(r_p)}{\left(\frac{r_p}{r}\right)}^2={\left(\frac{r}{r_p}\right)}^{1/4}{T}_{2B}^J.$$

$$T_{GW}(r,r_p)=\frac{R_s}{c}\frac{M}{m}{\left(\frac{r_p}{R_s}\right)}^4{\left(\frac{r}{r_p}\right)}^{1/2},$$

$${\left(\frac{R_s}{r_p}\right)}^{5/2}\left(\frac{M}{m}\right)\frac{1}{\ln \Lambda N(r)}=1.$$

$$\frac{r_p}{R_s}=(\ln \Lambda {)}^{-2/5}{\left(\frac{r}{R_h}\right)}^{-1/2}.$$

$$\mathcal{N}=(\ln \Lambda {)}^{-7/5}\left(\frac{M}{m}\right)\left(\frac{R_s}{R_h}\right){\left(\frac{r}{R_h}\right)}^{-11/4}.$$

$$r=(\ln \Lambda {)}^{-4/15}{R}_s^{2/3}{R}_h^{1/3}.$$

$$\frac{N(r)}{T_{GW}(r)}=\frac{1}{P(R_h)}{\left(\frac{r}{R_s}\right)}^{-2}.$$

$$N(r)={\left(\frac{r}{R_s}\right)}^2{\left(\frac{R_h}{R_s}\right)}^{-3/2}\frac{M}{m}\propto r^2.$$

$$N(r)=\frac{M}{m}\{\begin{array}{ll}{\left(\frac{r}{R_s}\right)}^2{\left(\frac{R_h}{R_s}\right)}^{-3/2}& \mathrm{for}r<R_s^{2/3}{R}_h^{1/3}\\ & \\ {\left(\frac{r}{R_h}\right)}^{5/4}& \mathrm{for}r>R_s^{2/3}{R}_h^{1/3}\end{array}$$

$$N(r,r_p)=\frac{M}{m}\{\begin{array}{ll}{\left(\frac{r_p}{R_s}\right)}^{3/2}{\left(\frac{r}{R_s}\right)}^{1/2}{\left(\frac{R_h}{R_s}\right)}^{-3/2}& \mathrm{for}{r}_p<R_s(r/R_h{)}^{-1/2}\\ & \\ {\left(\frac{r}{R_h}\right)}^{1/4}\left(\frac{r_p}{R_h}\right)& \mathrm{for}{r}_p>R_s(r/R_h{)}^{-1/2}\end{array}$$

$$\frac{r_p}{R_s}={\left(\frac{m}{M}\right)}^{2/3}\left(\frac{R_h}{R_s}\right){\left(\frac{r}{R_s}\right)}^{-1/3}.$$

$${\mathcal{R}}_{TDEs}=\frac{N(R_h)}{\ln (J_c/J_{LC})T_{2B}(R_h)}\approx \frac{1}{P(R_h)},$$

$$r_0=R_h{\left(\frac{R_s}{R_T}\right)}^2,$$

$${\mathcal{R}}_{MS-EMRIs}\approx \frac{N(r_0)}{T_{2B}(r_0)}=\frac{1}{P(R_h)}{\left(\frac{R_s}{R_T}\right)}^2.$$

$$\frac{{\mathcal{R}}_{MS-EMRIs}}{{\mathcal{R}}_{TDEs}}\approx {\left(\frac{R_s}{R_T}\right)}^2.$$

$$T_{\mathrm{coll}}^{-1}=\frac{N(r)}{r^3}{\left(\frac{GM}{r}\right)}^{1/2}{R}_{\ast }^2\left[1+{\left(\frac{Gm}{R_{\ast }}\frac{r}{GM}\right)}^2\right].$$

$$T_{\mathrm{coll}}^{-1}=\frac{N(r_p)}{r_p^{2}r}{\left(\frac{GM}{r}\right)}^{1/2}{R}_{\ast }^2.$$

$$r_p={\left[R_{\ast }\left(\frac{M}{m}\right)\ln {\Lambda }^{-1/2}\right]}^{-8}{r}^9.$$

$$r_p={\left(\frac{M}{m}\right)}^{-8/11}{\left(\frac{R_h}{R_s}\right)}^{5/11}{\left(\frac{R_{\ast }}{R_s}\right)}^{8/11}{\left(\frac{r}{R_s}\right)}^{4/11}{R}_s$$

$$r_{p,inj}$$

$$r_{inj}$$

$$\frac{r_{p,inj}}{R_s}<(\ln \Lambda {)}^{-2/5}{\left(\frac{r_{inj}}{R_h}\right)}^{-1/2}.$$

$$\alpha <{\alpha }_c=(\ln \Lambda {)}^{-4/15}\frac{R_s}{R_T}{\left(\frac{R_h}{R_s}\right)}^{1/3}{\left(\frac{m}{M}\right)}^{1/9}$$

$$\frac{r_{p,c}}{R_s}=(\ln \Lambda {)}^{-4/15}{\left(\frac{R_h}{R_s}\right)}^{1/3}{\left(\frac{m}{M}\right)}^{1/9},$$

$$\frac{r_c}{R_s}=(\ln \Lambda {)}^{-4/15}{\left(\frac{R_h}{R_s}\right)}^{1/3}{\left(\frac{M}{m}\right)}^{2/9}.$$

$$N(r)=\frac{M}{m}\times \{\begin{array}{ll}{\left(\frac{\eta R_h}{r_{min}}\right)}^{1/2}{\left(\frac{r}{R_h}\right)}^{5/4}& \mathrm{for}r<r_{min}\\ & \\ {\eta }^{1/2}{\left(\frac{r}{R_h}\right)}^{3/4}& \mathrm{for}{r}_{min}<r<\eta R_h\\ & \\ {\left(\frac{r}{R_h}\right)}^{5/4}& \mathrm{for}\eta R_h<r\end{array}$$

$$N(r)={\left(\frac{M}{m}\right)}^{5/9}{\left(\frac{R_s}{R_h}\right)}^{1/6}{\eta }_B\times \{\begin{array}{ll}{\left(\frac{r}{r_{p,c}}\right)}^4& r<r_{p,c}\\ {\left(\frac{r}{r_{p,c}}\right)}^{1/2}& r_{p,c}<r\end{array}$$

$$R_h=\frac{GM}{{\sigma }^2}\approx 2\mathrm{pc},$$

$$r_{in}=R_{\ast }{\left(\frac{M}{m}\right)}^{1/3}\approx 1\mathrm{AU}.$$

Full text

††thanks: E-mail: sari@phys.huji.ac.il; giacomo.fragione@mail.huji.ac.il

1. Introduction

Most of the galaxies over the whole Hubble sequence harbour super massive black holes (SMBHs), with masses in the range $10^{6}\ \mathrm{M}_{\odot}\lesssim M\lesssim 10^{10}\ \mathrm{M}_{\odot}$, in their innermost regions (Kormendy & Ho, 2013). Dense and complex structures of stars, stellar remnants and gas surround SMBHs (see Alexander (2017) for a recent comprehensive review). The Milky Way’s Galactic Center (GC) is the only nucleus close enough to resolve scales smaller than $\sim$ pc (Gillessen et al., 2017; Baumgardt et al., 2018; Gallego-Cano et al., 2018). The recent big advance in dedicated instruments, e.g. GRAVITY111https://www.eso.org/sci/facilities/paranal/instruments/gravity.html (Eisenhauer et al., 2011; Gravity Collaboration et al., 2018a, b), allows to observe with an unprecedented precision the GC, which marks an unique opportunity of improving the understanding of our Galaxy and galactic nuclei in general.

In the GC, stars and compact remnants move in the smooth near-Keplerian potential of the SMBH, which dominates the dynamics within the radius of influence $R_{h}$, beyond which the potential of the SMBH is overcome by the galactic field (Merritt, 2013). On timescales much longer than the orbital period, the microscopic fluctuations of the potential make stars energy and angular momentum diffuse, as a result of continuous non-coherent scatterings with other stars. Stars are subject to a net residual specific force $\propto\sqrt{N}$ ($N$ is the number of stars), and their energy and angular momentum diffuse on the typical two-body timescale $T_{\rm 2B}\gtrsim 10^{9}$-$10^{10}$ yr (Bar-Or & Alexander, 2014, 2016). As a consequence, a population of equal-mass stars rearrange their orbits and relax into a cuspy density profile $n\propto r^{-7/4}$ ($r$ is the radial distance with respect to the SMBH; Bahcall & Wolf, 1976, BW). Mass spectrum, and the subsequent dynamical friction, affects the typical slope, with more massive objects that relax in a steeper profile (Bahcall & Wolf, 1977; Hopman & Alexander, 2006a). Massive perturbers may play a role as well (Perets et al., 2007). Only galactic nuclei harbouring SMBHs less massive than $\sim 10^{7}\ \mathrm{M}_{\odot}$ have typical evolutionary timescales small enough to make the effects of the uncorrelated stellar interactions important within a Hubble time.

On timescales smaller than $T_{\rm 2B}$, but longer than the stars orbital period, the residual torque $\propto\sqrt{N}$ resulting from fluctuations of the average potential becomes relevant in shaping the stars orbits (Rauch & Tremaine, 1996). Both the direction and magnitude of the angular momentum (hence eccentricity) diffuse on a resonant relaxation timescale, $T_{\rm RR}\gtrsim 10^{7}$-$10^{9}$ yr. On even smaller timescales, only the transverse component of the residual torque has a non-negligible effect, thus shaping the inclinations of the orbital planes of the stars on a vector resonant relaxation timescale, $T_{\rm VRR}\gtrsim 10^{5}$-$10^{7}$ yr (Kocsis & Tremaine, 2011, 2015).

In addition to star-star interactions, the gravitational wave (GW) emission radiation can dissipate energy. When energy is efficiently dissipated in the innermost regions of the cusp, stars and stellar black holes gradually inspiral and illuminate the GW sky as an EMRI (extreme-mass ratio inspiral; Hopman & Alexander, 2006b; Aharon & Perets, 2016). In the last orbits, relativistic precession decouples the GW inspiral from the residual torques of the background stars, thus resulting in a slow inspiral that roughly conserves the pericenter of the orbit (Peters, 1964), while shrinking its semimajor axis by losing energy at each pericenter passage, until the object is swallowed by the SMBH. This object can either be a stellar black hole (BH-EMRI; Hopman & Alexander, 2005) or a main-sequence star (MS-EMRI; Linial & Sari, 2017). Characterizing BH-EMRIs and MS-EMRI is of extremely interest for the upcoming Laser Interferometer Space Antenna (LISA)222https://lisa.nasa.gov/, which is expected to probe SMBH demographics in galactic nuclei and cosmological parameters even at large redshifts (Amaro-Seoane et al., 2017).

While in the GW regime stars are disrupted after gradually losing energy, stars can be swallowed by directly plunging onto the SMBH. Unlike the GW mergers, the plunge disruption scenario requires the star to remain on a plunging orbit only long enough to pass through its periapsis once, where the star is disrupted and shines as a tidal disruption event (TDE) (Stone et al., 2013). The typical radius within which a star is disrupted is the tidal disruption radius $R_{T}\approx R_{*}(M/m)^{1/3}$, where $R_{*}$ and $m$ are the radius and mass of the star, respectively. Stars are driven onto plunging loss-cone orbits by the continues two-body scattering events within the SMBH sphere of influence, which randomizes the stars energy and angular momentum. However, the angular momentum evolves much faster than energy for eccentric orbits, thus rendering the angular momentum diffusion as the driving mechanism for producing TDEs (Alexander, 2017).

Both observations and theory suggest that stars can be transported to the innermost part of the cusp of stars by binary star disruptions (Hills, 1988; Brown et al., 2014; Brown, 2015). In the inner regions of galactic nuclei, binary stars undergo three-body exchange interactions with the SMBH, where one of the stars is ejected as hypervelocity star (HVS) with velocities of hundreds km s*-1*, while the former companion remains bound to the SMBH (Yu & Tremaine, 2003; Sari et al., 2010; Kobayashi et al., 2012; Rossi et al., 2014). Recently, Koposov et al. (2019) reported the observation of a HVS with a velocity of $1755\pm 50\,\rm km\,s^{-1}$ in the Galactic frame. Triple and quadruples stars may undergo the same fate as well (Fragione & Gualandris, 2018; Fragione, 2018). Other mechanisms have been proposed to explain the observed population of HVSs in our Galaxy (Yu & Tremaine, 2003; Capuzzo-Dolcetta & Fragione, 2015; Fragione & Capuzzo-Dolcetta, 2016; Fragione et al., 2017), but binary disruptions remain the favored scenario. However, recent analyses show a growing evidence of high-velocity objects possibly not originated in the GC (Boubert et al., 2018; Marchetti et al., 2018; de la Fuente Marcos & de la Fuente Marcos, 2019). Little attention has been devoted to the long term effects of binary disruptions in the ecology of galactic nuclei. Coleman Miller et al. (2005) showed that the tidal breakup of black hole binaries can produce events observable with low eccentricity in the LISA band. Bromley et al. (2012) examined the fate of the stars that remain bound to the SMBH after the binary disruption, and found that most of them undergo TDEs, thus fuelling the growth of the SMBH. Recently, Fragione & Sari (2018) have shown that continuous injection of stars enabled by binary disruptions may make the cusp slope steeper, whose extent and importance depend on the injection rate and survival fraction of the injected stars.

In this paper, we provide a simplified analytical treatment of the scattering processes in galactic stellar nuclei, assuming all stars have the same mass. This problems contains two dimensionless numbers. First, the mass ratio $M/m$, which is also the number $N_{h}$ of stars within the SMBH influence radius. Second, the ratio between $R_{h}$ and the SMBH Schwarzchild radius $R_{s}=2GM/c^{2}$. By coincidence, these dimensionless numbers are equal in our Galaxy, $M/m\sim R_{h}/R_{s}\sim 4\times 10^{6}$. More massive galaxies tend to have $N_{h}>R_{h}/R_{s}$, and vice versa. We also discuss the role of collisions, and focus on the effect of binary disruptions on the rate of TDEs and MS-EMRIs.

The paper is organized as follows. In Section 2, we discuss the typical two-body relaxation and GW timescales, and how they shape the cusp of stars surrounding the SMBH. In Section 3, we explain the role of binary disruptions, which is then investigated by means of numerical simulations in Section 4. Finally, in Section 5, we discuss the implications of our findings and draw our conclusions.

2. Two-body relaxation and gravitational waves

Consider a stellar cusp, with $N(r)$ stars, each of mass $m$ and with semimajor axis $r$ around a SMBH of mass $M$. This is related to the number density of stars by $N(r)\approx n(r)r^{3}$. The two-body relaxation time at radius $r$ is given by (Bahcall & Wolf, 1976)

[TABLE]

where $P(r)$ is the star’s orbital period and $\ln\Lambda\sim 10$ is the Coulomb logarithm. This is the timescale by which a star at a roughly circular orbit changes it angular momentum and energy by a factor of order unity.

For stars on very eccentric orbits, the orbital angular momentum $J$ is much smaller than the circular angular momentum at the same semimajor axis $J_{c}$, and the timescale to significantly change their angular momentum is smaller by $(J/J_{c})^{2}$. If we denote their periapsis distance by $r_{p}\approx(J/J_{c})^{2}r$, then we can define the eccentric orbit relaxation time as (Binney & Tremaine, 1987)

[TABLE]

For very eccentric orbits with $r_{p}\ll r$, this is the typical timescale to change the periapsis distance of the star’s orbit. In the previous expression, we assumed that angular momentum changes are dominated by interactions at $r$. This holds as long as $N(r)r$ is an increasing function of $r$, so as long as the density profile is shallower than $n(r)\propto r^{-4}$, which we will assume is the case.

Very eccentric stars also relax their energy faster than stars on circular orbits of the same semimajor axis. This is because the kicks they get from the inner denser cluster are more significant to their energy change. The cross section for some energy change $\Delta E$ is given by $(Gm/\Delta E)^{2}$, independent of the velocities of the stars. As a consequence, the inner cusp will dominate the scatterings as long as $n(r)r$ is a decreasing function of $r$, which we will assume to be true. The typical timescale to change the semimajor axis for an eccentric star is therefore given by (Binney & Tremaine, 1987)

[TABLE]

Energy relaxation for highly eccentric orbits is therefore slightly less efficient than that of angular momentum, if the BW density profile ($n\propto r^{-7/4}$) is assumed. Density profiles steeper than $\rho(r)\propto r^{-2}$ will result in energy relaxation being faster than that of angular momentum.

As discussed, GWs circularize the orbit of highly eccentric orbits, while keeping the periapsis distance roughly fixed. This happens on a timescale (Peters, 1964; Hopman & Alexander, 2006a)

[TABLE]

where $R_{s}=2GM/c^{2}$ is the SMBH Schwarzchild radius.

Equating the timescales of GWs (Eq. 4) to that of periapsis evolution by two-body encounters (Eq. 2)

[TABLE]

In the case of a BW cusp profile, $N(r)=(M/m)(r/R_{h})^{5/4}$, hence

[TABLE]

Therefore, a star at the radius of influence with pariapse close to the Schwarzschild radius has equal time to shrink its semimajor axis by GWs and to change its periapsis by scatterings. In the $r_{p}$-$r$ plane, stars which are above the line given by Eq. 6, evolve primarily due to two-body scatterings, while stars below it mostly shrink in semimajor axis by emission of GWs. This is given by line $T_{GW}=T_{2B}$ in Fig. 1. It is remarkable that this line passes through the point $(r=R_{h},r_{p}=R_{s})$ regardless of the two dimensionless parameters in the problem, $M/m$ and $R_{h}/R_{s}$ (neglecting the factor $(\ln\Lambda)^{-2/5}$).

We note that the above calculation assumes that it takes more than one orbit to change the stars’ semimajor axis or periapsis significantly. On the line given by Eq. (6), the number of orbits in such a timescale is (assuming an BW profile)

[TABLE]

For our Galaxy, the first product is of order unity (by coincidence), so it takes about a single orbit for a star at the radius of influence and $r_{p}=R_{s}$ to change it periapsis or apoapsis significantly. However, anywhere within the radius of influence, $\mathcal{N}\gg 1$, thus our previous assumption is satisfied for our Galaxy. Using the $M\propto\sigma^{4}$ relation (Merritt & Ferrarese, 2001), this is even more readily satisfied for any larger galaxy, but only marginally at the radius of influence of smaller galaxies. Note that $\mathcal{N}>1$ is also the condition that a loss cone defined by the line (6) is empty. Usually, the loss cone is defined by the tidal radius or the Schwarzschild radius. Yet, it makes more sense to treat Eq. 6 as defining the loss cone, as for smaller angular momenta gravitational waves deplete the stars and provide an effective zero boundary condition for the diffusion the occurs for larger angular momenta by two body interactions.

As mentioned, there are two important dimensionless parameters in this problem, $\left(M/m\right)$ and $\left(R_{h}/R_{s}\right)$. The relation between these parameters is important. The above discussion shows that $\left(M/m\right)\gg\left(R_{h}/R_{s}\right)$ results in an empty loss cone at the radius of influence, and with a loss cone defined by the Schwarzschild radius. For our Galaxy, by coincidence, both are $\left(M/m\right)\approx\left(R_{h}/R_{s}\right)\approx 4\times 10^{6}$, then it is a borderline situation. However, as we discuss later, the tidal radius in our Galaxy for main sequence stars is larger than the event horizon, and the loss cone associated with the tidal radius is therefore empty at the radius of influence.

The line of equal timescale cuts that of a circular orbit, i.e. $r=r_{p}$, for

[TABLE]

For our Galaxy, this is $\sim 100$ times the Schwarzschild radius. Below this distance, the BW cusp is modified and GWs shrinks the orbit. Interestingly, a constant flux of stars through circular orbits, would imply $N(r)\propto r^{4}$, i.e. a decreasing density towards the center. However, the flux of stars on circular orbits decaying by GWs is not constant as it is supplemented by the flux of starts evolving first by two-body encounters and then circularizing by GWs (along the curving dark cones of Fig. 1). Therefore, instead, we demand that the flux of stars through circular orbits of size $r<R_{s}^{2/3}R_{h}^{1/3}$ 333We are neglecting the factor $(\ln\Lambda)^{-4/15}\sim 0.6$. would equal the flux of stars supplied by two-body encounter into $r_{p}=r$, i.e. from initial circular orbits of size $r_{i}=R_{h}(r/R_{s})^{-2}$. The flux of these is

[TABLE]

This implies

[TABLE]

The overall cusp is therefore given by

[TABLE]

For our Galaxy, this implies a single star below a distance of $\sim 50R_{s}$. Yet, this point is within the regime of GWs, which act on each star separately. The low number of stars, even when it falls below unity, does not invalidate our solution and represents the expectation value for the number of stars.

It is easy to generalize this equation to calculate also the number of eccentric orbits, i.e., instead of $N(r)$, we could estimate $N(r,r_{p})$ - the number of stars with semimajor axis smaller than $r$ and periapsis smaller than $r_{p}$. This is given by

[TABLE]

Now, $N(r,r_{p})=1$ is given by

[TABLE]

Accidentally, for the parameters of our Galaxy, there is a single object with semimajor axis equal to the radius of influence, but periapsis as small as the Schwarzschild radius. This is another consequence of the coincidence $\left(M/m\right)\approx\left(R_{h}/R_{s}\right)$.

2.1. Tidal disruption events and main-sequence extreme-mass ratio inspirals

For small enough galaxies, tidal disruption of main sequence stars occurs before the star hits the Schwarzschild radius (Stone et al., 2013). If this encounter happens on an extremely eccentric orbit, this result in a violent TDE of the star, over its dynamical time. If, instead, the star approaches on a slowly shrinking circular orbit, the star transfers mass to the SMBH on the GW timescale (Dai & Blandford, 2013; Linial & Sari, 2017). These events are know as main sequence MS-EMRIs. Linial & Sari (2017) showed that for main-sequence stars mass transfer may result is an expanding orbit causing an inverted ”Chirp” signal which (spelled backwards) we call a ”Prich”, where the GW frequency decreases with time. These events may also appear as a sequence of TDE-like flares, if two (or more) consecutive MS-EMRIs collide (Metzger & Stone, 2017). Figure 1 illustrates the different origins of TDE and MS-EMRI events.

In this picture, the rate of TDEs is dominated by the supply from roughly circular orbits at the radius of influence $R_{h}$

[TABLE]

simply the inverse period at the radius of influence. In the previous equation, $J_{c}$ and $J_{LC}=\sqrt{2GMR_{T}}$ are the circular and loss-cone angular momentum, respectively. On the other hand, the rate of MS-EMRIs is given by the supply from roughly circular orbits of the largest radius that does not result in a TDE. The periapsis of these orbits, once their evolution is dominated by GWs, is just above the tidal radius $R_{T}$. Their initial radius is therefore

[TABLE]

and the rate of these events is given by

[TABLE]

The ratio between the rate of MS-EMRIs and that of TDEs is simply

[TABLE]

For our Galaxy, and Solar-like stars, $R_{T}\sim 10R_{s}$, so the MS-EMRI rate is about a percent of the TDE rate.

2.2. Collisions

Close enough to the SMBH, the orbital velocities become larger than the typical escape velocity from the surface of a star ($v_{esc}\sim 600\,\rm km\,s^{-1}$ for $1\,\rm M_{\odot}$ star), and collisions become more likely than scattering events. Whether it results in a merger or a destruction, depends on the mass ratio of the colliding stars, the ratio between their relative velocity to their surface escape speed as well as the impact parameter of the collision (grazing or head on) (Benz & Hills, 1987; Trac et al., 2007; Gaburov et al., 2010). This lead Alexander (2017) to define the collision radius where the velocity dispersion equals the escape velocity. He estimated the collision rate, the inverse collision time, as

[TABLE]

where the term in parenthesis is the gravitational focusing term, increasing the effective cross section for collisions beyond the physical size of the star. However, this term is significantly larger than unity only if the velocity dispersion of the stars is smaller than their escape speed and in that case the collision time is longer than the relaxation time. We can therefore ignore gravitational focusing when collisions are important.

However, as we have shown, the interesting objects observationally, those that lead to tidal disruption events or those that start evolve quickly due to gravitational waves, are extremely eccentric. The collisional time is not determined only by their semimajor axis, but also depends on their periapsis distance. We therefore have to generalize the result of Alexander (2017) to highly eccentric orbits. Collisions around pericenter would dominate as long as the optical depth, $N(r)/r^{2}$, decreases with radius. Such is the case for a Bahcall-Wolf cusp or any density profile with $\alpha>1$. We therefore use the eccentric collisional time as

[TABLE]

Equating $T_{\rm coll}$ to $T_{2B}^{J}$ (Eq. 2), we find

[TABLE]

Note that the quantity in square brackets is the collisional radius result of Alexander (2017).

For stars on small enough periapsis and semimajor axis GW rather than two body scattering dominate the evolution. Equating the collision time to $T_{GW}$ (Eq. 4) we obtain:

[TABLE]

We show the phase space where collisions are important in Fig. 1. We do not take into account collisions in our calculations, but it is clear that they would affect a significant portion of the phase space. A realistic description of galactic nuclei must take collisions into account. We leave this to future work.

3. The effect of binaries

Binary disruption results in one star ejected, perhaps as a HVS, while the other remains tightly bound to the SMBH Hills (1988). Fragione & Sari (2018) showed that binary breakup is an important source of stars close to the SMBH, and may modify the cusp structure if binaries are sufficiently abundant.

These stars are injected with periapsis distance ($r_{p,inj}$) comparable to that of the binary tidal radius, and a semimajor axis ($r_{inj}$) larger than that by a factor of $\sim(M/m)^{1/3}\approx 100$, almost independently of the original semimajor axis (as long as it is larger than that)

[TABLE]

In the previous equations, $\alpha$ is the initial binary semimajor axis $a_{b}$ in terms of the stellar radius. The stars injected into the tightest orbits are those who broke from an almost contact binary, i.e. $\alpha$ is of the order of a few.

The criterion for these injected stars to be mostly dominated by GW evolution is (see Eq. 6)

[TABLE]

We now define a critical binary separation parameter $\alpha_{c}$ above which the captured member of the binary is mostly affected by two-body scatterings

[TABLE]

The binaries with orbital separation parameter of $\alpha_{c}$ dissolve on a periapsis distance of

[TABLE]

while their semimajor axis is larger by a factor of $(M/m)^{1/3}$

[TABLE]

For the parameters of the Milky Way and Solar-mass main sequence stars, $R_{T}/R_{s}\approx 10$, hence $\alpha_{c}\approx 3$.

If the fraction of disrupted binaries is large enough, the injected stars modify the cusp profile. This increases the two-body encounter rate, and shortens the two-body relaxation time, somewhat reducing the value of $\alpha_{c}$. Fragione & Sari (2018) showed that the cusp obtains a shape of

[TABLE]

where $\eta$ is the fraction of binaries at the influence radius $R_{h}$, and $r_{min}$ is the minimum injected radius, which depends on the minimum separation of binaries. The outer section of the cusp, $\eta R_{h}<r<R_{h}$, is unaffected by the binaries. We also note that the extent of the modified cusp also depends on the fraction of the injected stars that survive without being dissolved either as TDEs or due to GW emission. Fragione & Sari (2018) estimated that $\sim 70$-$80$% of the injected stars diffuse to TDE orbits.

However, once we include the influence of GWs, initially tight binaries that got dissolved can be circularized by two-body-scatterings into orbits smaller than $r_{p,c}$. We can derive this further population by requiring that the flux of stars $N(r)/T_{GW}$ equals the rate of injected from binary disruptions. This adds an additional population of stars on very tight orbits in the following amount:

[TABLE]

In the above equation, $\eta_{B}=\eta/\ln(a_{\max}/a_{\min})$ is the fraction of tight binaries per logarithmic unit of separation. Assuming the fraction of binaries is $\eta=0.1$ and a log-uniform binary separation distribution (see Sect. 4.1) in the range $a_{\rm min}=0.01$ AU-$a_{\rm max}=1$ AU for solar-mass stars, then $\eta_{B}\sim 0.04$. For the Milky Way, the coefficient in Eq. 29 for this extra-component is $400\eta_{B}\approx 17$. As a consequence, we expect a population of $\sim 20$ Solar-mass stars in the GC, orbiting at $r_{p,c}\approx 30R_{S}\approx 3{\rm AU}$ at orbital periods below a day.

4. Numerical simulations

Direct $N$-body simulations would represent the ideal tool to investigate the evolution and distribution of stars in the proximity of a SMBH. Unfortunately, prohibitive computational times limit $N$-body simulations to a small number of stars, of the order of $\sim 50$k stars, thus far from the $\sim 10^{6}$ stars that reside in the SMBH sphere of influence (Baumgardt et al., 2004a, b; Baumgardt, 2017).

In what follows, we describe the computational method we adopt to perform long-term evolution of the stellar cusp around the Milky Way’s SMBH. Compared to Fragione & Sari (2018), we upgraded the scheme by including angular momentum evolution and GW energy loss. To summarize, our method has four main features

• •

It follows the 2-D evolution of the energy and angular momentum (or semimajor axis and eccentricity, respectively) of each star, therefore conserving the number of stars;

• •

The rate of scattering changes with energy and angular momentum, and with the number of particles that have similar energy and angular momentum;

• •

It takes into account relativistic GW effects;

• •

It takes into account possible source terms, e.g. from breaking of binaries.

We focus our attention on the region inside the sphere of influence, where the stellar dynamics is dominates by the SMBH field

[TABLE]

where $\sigma$ is the velocity dispersion external to the radius of influence. In all our calculations, we assume $M=4\times 10^{6}\ \mathrm{M}_{\odot}$ (Gillessen et al., 2017) and a single-mass population of stars of mass $m=1\ \mathrm{M}_{\odot}$. The region of interest spans a wide range of energy, angular momentum and distances with respect to the SMBH. The innermost radius we take into account is the tidal disruption radius of $1\ \mathrm{M}_{\odot}$ star (Stone et al., 2013)

[TABLE]

The semimajor axis and eccentricity of the stars are continuous variables, but we divide the $(r;r_{p})$ space into bins, thus considering a 2D grid. Every time step we count how many stars $N(i,j)$ ($i$ is the index for the semimajor axis, $r_{i}$, and $j$ is the index for the pericenter, $r_{p,j}$) are in a given square of the grid, and how many stars are in a given semimajor axis bin, regardless of the pericenter

[TABLE]

To determine the timestep, we compute the angular momentum two-body timescale for each square of the grid (see Eq. 2)

[TABLE]

where

[TABLE]

Here, we set $f_{t}=0.1$, a reduction factor for the timescale, and $\gamma=1.5$, a factor that takes into account the bin size (Fragione & Sari, 2018). In analogy to Fragione & Sari (2018), we choose as overall timestep the minimum of the two-body timescales

[TABLE]

such that the total time at the k-th step is $T_{k}=T_{k-1}+\Delta T$. Updating at every time step all the stars in each square of the 2D grid is time-consuming. Hence, we check for each square of the 2D grid if

[TABLE]

where $L(i,j)$ is the moment when we made the last update (with respect to the overall time) of a given square $(i,j)$ of the grid. If Eq. 36 is satisfied, we update the energy and angular momentum of the stars in the square $(i,j)$. For each of these stars, we find the different bins explored during its orbital motion and identify $W1$ and $W2$ such that the orbit of the star is entirely within the bins $(r_{W1},r_{W2})$. Then, we compute the average impact parameter $B(W)$ in each of the bins it went through as

[TABLE]

where $P_{*}$ is the orbital period of the given star. The typical shifts in energy and angular momentum are computed as

[TABLE]

where $v(W)=\sqrt{GM/r_{W}}$ is the velocity of the star in the bin. Finally, we update the energy and the angular momentum of each star in each bin

[TABLE]

where $0\leq\chi<2\pi$ and $0\leq\Phi<2\pi$ are drawn randomly from a uniform distribution.

After updating energy and angular momentum, we check the following conditions

• •

if $J_{new}\leq J_{LC}=\sqrt{2GMR_{T}}$, the star is considered a TDE and removed

• •

if $r_{new}=GM/E_{new}\geq R_{h}$, the star is considered escaped from the cusp and removed

In these cases, a new star is generated randomly in the last bin and with eccentricity drawn from a thermal distribution (Fragione & Sari, 2018).

In the case GW effects are taken into account444We take into account GWs for all the stars with $r_{p}<0.1R_{h}$., we update the energy and angular momentum of the stars as (Hopman & Alexander, 2006a, b)

[TABLE]

where

[TABLE]

and

[TABLE]

If binary injection is taken into account, we parametrize the rate of disrupted binaries with the dimensionless parameter $\eta$. Following Fragione & Sari (2018), we generate in each timestep

[TABLE]

injected stars as a consequence of the binary tidal disruption, that are added to the pre-existing population. We sample injected star semimajor according to $f(a)\propto 1/a$, and we set their eccentricities to $1-(m/M)^{1/3}\approx 0.99$ (Eq. 23; see also Brown, 2015; Alexander, 2017).

4.1. Cusp slope, tidal disruption events and main-sequence extreme-mass ratio inspirals

In our simulations, we consider $0\leq\eta\leq 1$ ($\eta=0$ means no dissolved binaries are injected), and study how the cusp slope, the TDE rate and the MS-EMRI rate change as a function of the binary fraction $\eta$. We consider $10^{6}$ stars and start with a stable Bahcall & Wolf (1976) cusp, i.e star semimajor axis follow a distribution $N(r)\propto r^{5/4}$. For initial momenta, we sample uniformly $0\leq(J/J_{c})^{2}\leq 1$ ($J_{c}=\sqrt{GMr}$ is the circular angular momentum of a given semimajor axis $r$), i.e. thermal distribution of eccentricities. When dissolved binaries are injected, we assume that the binary semimajor axis are distributed according to (Duquennoy & Mayor, 1991)

[TABLE]

in the interval $(a_{min}$,$a_{max})$ at the SMBH radius of influence, and set the semimajor axis of the captured stars according to the same distribution scaled by a factor $\sim(M/m)^{2/3}$ (Fragione & Sari, 2018)

[TABLE]

in the range $r_{min}=(M/m)^{2/3}a_{min}$ and $r_{max}=(M/m)^{2/3}a_{max}$, with eccentricity $e\sim 0.99$. For Solar-mass stars, we chose the minimum as $a_{min}=0.01$ AU, while for the maximum we set $a_{max}=0.1$ AU, since binaries with larger semimajor axis are typically disrupted by the background stars rather than by the SMBH (Hopman, 2009).

Figure 3 illustrates the profile $N(r)/N_{h}$ of the simulations along with the theoretical curves $\propto r^{5/4}$ and $\propto r^{3/4}$ from Eq. 28. For $\eta=0$, our results recover the standard Bahcall & Wolf (1976) solution, where energy conservation dictates $\alpha=5/4$. We also report the results of the simulations if stars are injected as a consequence of tidal binary disruption when $\eta=1$. As discussed, while the number density has the same slope of the Bahcall-Wolf solution for $r<r_{inj}$ (with a higher normalization), the cusp develops a steeper profile for $r_{inj}<r<\eta R$. The extent of the steeper profile depends both on the binary fraction near the SMBH influence radius, $\eta$, and on the fraction of injected stars that are not disrupted by the SMBH, $\omega<1$. As previously discussed (see Eq. 29), a further $\ln(a_{\max}/a_{\min})$ factor has to be taken into account. As a consequence, the effective fraction to consider in Eq. 28 would be

[TABLE]

In Fragione & Sari (2018), we discussed the fraction of the stars injected in the vicinity of the SMBH from dissolved binaries that would circularize even though the high eccentricity of its orbit. We considered a Brownian process governed by a continuous diffusion equation, and found that the ratio of fluxes upwards to circular orbits compared to that downward to tidal disruption orbits is

[TABLE]

where $J_{0}$ is the angular momentum of the injected star. For the parameters of our Galaxy, $\ln(J_{0}/J_{LC})\sim\sqrt{\alpha}\sim 2$ and $J_{c}/J_{0}\sim\sqrt{r_{inj}/r_{p}}\sim(M/m)^{1/6}\sim 10$. The ratio is therefore $\omega\sim 0.3$, and $\eta_{eff}\sim 0.1\eta$. As a consequence, $\sim 70$% of the injected stars diffuse to TDE orbits, enhancing their rates. From Eq. 51, only $\sim 10\%$ of the injected stars are not disrupted and the cusp results steeper from $r_{min}\sim 250$ AU up to $R=0.1\eta R_{h}$. In Fig. 3, we show that the results of our simulation nicely follow Eq. 28, with $\eta_{eff}=0.1\eta$.

In the previous Sections, we discussed the orbits that origin TDEs and MS-EMRIs. The stars that undergo TDEs come from the outskirt of the SMBH influence radius on very eccentric orbits ($1-e\sim 10^{-5}$). In this case, the dominant effect is the continuous two-body interactions with other stars, which scatters stars onto plunging orbits, thus disrupting them within one orbital period. The stars that end their lives as MS-EMRIs inspiral by losing gradually their energy due to GW emission within the region $T_{GW}<T_{2B}$, while keeping the pericenter roughly constant. For $\eta=0$, we find that the rates are $\Gamma_{TDE}\sim 3.3\times 10^{-6}$ yr*-1* for TDEs and $\Gamma_{MS-EMRI}\sim 4\times 10^{-8}$ yr*-1* for MS-EMRIs, in nice agreement with the analytical predictions of Eq. 17. The ratio is $\Gamma_{MS-EMRI}/\Gamma_{TDE}\sim 1.2$ %.

Injection of dissolved binaries may modify the above picture. As shown in Fig. 2, disrupted binaries would leave injected stars both in the two-body driven region and in the GW driven region ($T_{GW}>T_{2B}$ and $T_{GW}<T_{2B}$, respectively), thus possibly enhancing both rates. Moreover, the injection of binaries may modify the distribution of stars (Fragione & Sari, 2018), thus shifting the $T_{GW}=T_{2B}$ boundary line. In general, we expect different TDE and MS-EMRI rates from the case with no injection of stars ($\eta=0$). Figure 4 illustrates the TDE and MS-EMRI rates as a function of the binary fraction $\eta$. While even for small $\eta$’s the MS-EMRI rate becomes larger than the case of no binaries, the TDE rate is significantly larger only when $\eta\gtrsim 0.5$, and becomes $\sim 2$ times larger than the case of no binary disruptions. While with no injection of binaries $\Gamma_{MS-EMRI}$ is just a few percent of $\Gamma_{TDE}$, the MS-EMRI rate becomes $\sim 10$ times larger for $\eta\sim 0.1$ and $\sim 1/5$ of the TDE rate if the fraction of binaries is significant ($\eta\gtrsim 0.5$).

5. Discussions and Conclusions

The recent big advance in instruments dedicated to the observation of our GC SMBH and the future GW LISA mission offer the unprecedented opportunity to test theoretical models of the densest environment in the Universe (Amaro-Seoane et al., 2017; Gravity Collaboration et al., 2018a, b). In galactic nuclei harbouring an SMBH, the short-term dynamics is dominated by the deep potential well of the SMBH, while two-body interactions shape the distribution of stars on longer timescales, modifying their energy and angular momentum. Stars revolving around the SMBH near the edge of its influence sphere may be scattered into extremely-eccentric orbits, which eventually end up with the tidal disruption of the star itself if it crosses the tidal disruption sphere. Moreover, stars sufficiently close to the SMBH also suffer from energy loss due to GW radiation emission. If the energy loss is rapid enough, these stars will gradually inspiral onto the SMBH on a Peters timescale, becoming luminous in GWs in the mHz frequency band (Alexander, 2017).

In this paper, we have provided a simplified analytical treatment of the scattering processes in galactic stellar nuclei, assuming all stars have the same mass. We have discussed how the interplay between two-body relaxation and gravitational wave emission modifies the slope of the inner cusp within the SMBH sphere of influence. We have calculated the TDE rate of stars that are tidally disrupted by the SMBH and the rate of stars that end up their lives as MS-EMRI, disrupted by the SMBH due to the slow emission of GW radiation. We have found that typically the rate in the latter case is just a few percent of the former, thus implying a few MS-EMRI events for hundred TDEs.

We have also discussed the role of binary disruptions in the ecology of the galactic nuclei dynamics. In Fragione & Sari (2018), we showed that stars injected on highly-eccentric orbits in the vicinity of the SMBH due to Hills binary disruption may modify the shape of the density cusp, if the injection rate is large enough. The high eccentricity of these injected stars makes them of high relevance in the context of TDE and MS-EMRI events, since they can either be disrupted by the SMBH on plunging orbits or on gradually-circularized orbits due to GW emission. We have shown that the MS-EMRI rate can almost approach the TDE rate in case the binary fraction at the SMBH influence radius is close to unity.

We note that in our model we do not take into account the coherent torques between slowly precessing orbits, i.e. the resonant relaxation process (Kocsis & Tremaine, 2011, 2015). Such mechanisms are expected to be relevant for $r\sim 1000$ AU, until relativistic precession decouples the GW inspiral from the residual torques of the background stars. This may be the case of the S-stars, whose orbits may have evolved due to the resonant relaxation process (Perets et al., 2009; Perets & Gualandris, 2010; Gillessen et al., 2017). Nevertheless, Bar-Or & Alexander (2016) showed that the inclusion of the resonant relaxation in the dynamical processes has little impact on the estimated rates of TDEs and EMRIs.

We have ignored the role of physical collisions. Close enough to the SMBH, the orbital velocities become larger than the typical escape velocity from the surface of a star ($v_{esc}\sim 600\,\rm km\,s^{-1}$ for $1\,\rm M_{\odot}$ star), and collisions become more likely than scattering events (Alexander, 2017). We generalize this analysis and take into account highly eccentric orbits which are the sources of dynamical events like TDEs and MS-EMRIs. We find the phase space where collisions are important, and plot that in Fig. 1. The likely outcome of a stellar collision depends, apart from the relative mass ratio of the colliding stars, on the ratio between their relative velocity to their surface escape speed (Benz & Hills, 1987; Trac et al., 2007; Gaburov et al., 2010). Due to the large velocity dispersion in galactic nuclei, the physical collision of two equal-mass $1\,\rm M_{\odot}$ stars will likely eject half of the total, thus creating a blue straggler, in the outermost regions of the cusp. Closer to the SMBH, the collisions could lead to stellar destruction rather than to mergers since the kinetic energy in the colliding star exceeds the binding energy. Collisions have been shown to possibly play some role in the depletion of the red giant population (Dale et al., 2009; Bartko et al., 2010; Schödel et al., 2018), but they may have also some importance in shaping the distribution of solar-mass stars Sills et al. (2005); Dale & Davies (2006), whose impact on the galactic nuclei economy deserves future attention.

Finally, we have not considered a mass function both for the stars in the cusp and the injected stars, but only a single-mass population of $1\ \mathrm{M}_{\odot}$ stars. As a consequence, we also do not account for mass segregation as due to dynamical friction. Objects of different masses have generally different slopes of the cusp density. The more massive the object the steeper the cusp profile due to more efficient segregation towards the center, with stellar black holes leading with the largest value $\rho_{BH}(r)\propto r^{-\alpha_{BH}}$, where $\alpha_{BH}\approx 1.5$-$2$ Bahcall & Wolf (1976, 1977); Alexander & Hopman (2009); Aharon & Perets (2016); Baumgardt et al. (2018). Massive stellar remnants would also reduce the two-body relaxation time at small radii (Vasiliev, 2019), thus changing the relative ratio of MS-EMRIs and TDEs (Eq. 17) in favour of MS-EMRIs. We leave a detailed calculation to a future study. This also might effect the formation of binary stars and black holes that can merge as a consequence of the Kozai-Lidov effect (Antonini & Perets, 2012; Stephan et al., 2016; Fragione et al., 2018b; Grishin et al., 2018; Hoang et al., 2018; Fragione & Antonini, 2019; Fragione et al., 2019; Stephan et al., 2019). All the picture is even more complicated by the possible presence of intermediate mass black holes and other remnants brought by inspiralling star clusters (Fragione et al., 2018a, c).

6. Acknowledgements

We thank Scott Tremaine and Brian Metzger for useful comments and discussions. RS is supported by an iCore and an ISF grant. GF is supported by the Foreign Postdoctoral Fellowship Program of the Israel Academy of Sciences and Humanities. GF also acknowledges support from an Arskin postdoctoral fellowship at the Hebrew University of Jerusalem.