Relaxation near supermassive black holes driven by nuclear spiral arms: anisotropic hypervelocity stars, S-stars and tidal disruption events

TL;DR

This paper investigates how nuclear spiral arms in galactic centers influence star and binary star dynamics near supermassive black holes, leading to increased tidal disruption events, hypervelocity stars, and possibly affecting black hole mergers.

Contribution

It introduces a model for diffusion caused by nuclear spiral arms, quantifies their impact on star disruption rates, and explains the anisotropic distribution of hypervelocity stars.

Findings

TDE rates increase by up to 20% due to spiral arms.

Binary disruption rates are enhanced by up to a factor of 100.

The scenario explains the observed anisotropic distribution of HVSs.

Abstract

Nuclear spiral arms are small-scale transient spiral structures found in the centers of galaxies. Similarly to their galactic-scale counterparts, nuclear spiral arms can perturb the orbits of stars. In the case of the Galactic Center (GC), these perturbations can affect the orbits of stars and binaries in a region extending to several hundred parsecs around the supermassive black hole (MBH), causing diffusion in orbital energy and angular momentum. This diffusion process can drive stars and binaries to close approaches with the MBH, disrupting single stars in tidal disruption events (TDEs), or disrupting binaries, leaving a star tightly bound to the MBH, and an unbound star escaping the galaxy, i.e., a hypervelocity star (HVS). Here, we consider diffusion by nuclear spiral arms in galactic nuclei, specifying to the Milky Way GC. We determine nuclear spiral arm-driven diffusion rates using test-particle integrations, and compute disruption rates. Our TDE rates are up to 20% higher compared to relaxation by single stars. For binaries, the enhancement is up to a factor of ~100, and our rates are comparable to the observed numbers of HVSs and S-stars. Our scenario is complementary to relaxation driven by massive perturbers. In addition, our rates depend on the inclination of the binary with respect to the Galactic plane. Therefore, our scenario provides a novel potential source for the observed anisotropic distribution of HVSs. Nuclear spiral arms may also be important for accelerating the coalescence of binary MBHs, and for supplying nuclear star clusters with stars and gas.

Tables

Table 1. Table 1 Overview of the initial conditions for the Monte-Carlo integrations in Section  3 . The values for i 𝑖 i , r init subscript 𝑟 init r_{\mathrm{init}} and ϕ ˙ init subscript ˙ italic-ϕ init \dot{\phi}_{\mathrm{init}} have linear, logarithmic, and linear spacings, respectively. In the bottom row, x 𝑥 x is a random number between 0 and 1.

Symbol	Description	Value(s)
Grid parameters
$i$	Initial orbital inclination w.r.t. the Galactic plane	$0^{\circ }$-${180}^{\circ }$ (20 values; linear)
$N_s$	Number of spiral arm events	[1, 10]
${\tilde{\Sigma }}_{s,0,i}$	Spiral arm gas surface density	[100,500,1000] $M_{\odot }{\mathrm{pc}}^{-2}$
Fixed parameters
$r_{\mathrm{init}}$	Initial distance to the center	10-1000 pc (20 values; log.)
${\dot{\varphi }}_{\mathrm{init}}$	Initial $\varphi$ time derivative	$\left(0.1-1\right)\frac{v_{c,\mathrm{init}}}{r_{\mathrm{init}}}$ (5 values; linear)
$t_{\mathrm{int}}$	Integration time	$100\mathrm{Myr}$
$m$	Number of spiral arms (per event)	2
${\Omega }_s$	Spiral arm pattern speed	$50\mathrm{km}{s}^{-1}{\mathrm{kpc}}^{-1}$
${\sigma }_{s,i}$	Spiral arm event temporal Gaussian width	$50\mathrm{Myr}$ ($N_s=1$)
		$10\mathrm{Myr}$ ($N_s=10$)
$t_{s,i}$	Spiral arm event times	$50\mathrm{Myr}$ ($N_s=1$)
		$10-90\mathrm{Myr}$ ($N_s=10$)
$R_0$	Spiral arm phase parameter	$x$ pc

Table 2. Table 2 An overview of the binary models adopted for the rate calculations in Section  4 .

Model	Description	$M_{\mathrm{bin}}$	$M_{1,\mathrm{low}}$	$M_{1,\mathrm{up}}$	$f_{\mathrm{PDMF}}$	$f_{\mathrm{bin}}$	$t_{\star }$	Period	$a_{\mathrm{bin},\min }$	$a_{\mathrm{bin},\max }$
		($M_{\odot }$)	($M_{\odot }$)	($M_{\odot }$)			($\mathrm{Myr}$)	distribution	($\mathrm{AU}$)	($\mathrm{AU}$)
M1	Solar type	2	1	2	$9.41\times {10}^{-1}$	0.4	$1.1\times {10}^4$	lognormal	$3.35\times {10}^{-2}$	50
M2	HVS	8	2.5	4.5	$2.15\times {10}^{-2}$	0.6	179.2	MDS1	$5.32\times {10}^{-2}$	0.984
M3	S-star	16	5	15	$1.79\times {10}^{-3}$	0.8	37.2	MDS1	$6.70\times {10}^{-2}$	6
M4	O star	30	15	120	$3.66\times {10}^{-6}$	1.0	12.8	MDS1	$8.27\times {10}^{-2}$	50

Table 3. Table 3 Top part: disruption rates Γ Γ \Gamma of single stars (‘Single’) and binary stars (‘Binary’) due to relaxation driven by stars (first data row), and driven by nuclear spiral arms (in addition to stars; all other data rows). In the latter case, each row corresponds to a different inclination i 𝑖 i (expressed in degrees), number of spiral arm events N s subscript 𝑁 s N_{\mathrm{s}} , and (constant) spiral arm surface density Σ ~ s , 0 , i subscript ~ Σ s 0 𝑖 \tilde{\Sigma}_{\mathrm{s},0,i} (in units of M ⊙ ​ pc − 2 subscript M direct-product superscript pc 2 \mathrm{M}_{\odot}\,\mathrm{pc^{-2}} ). For binary stars, four models, ‘M1’ through ‘M4’, are assumed (cf. Table  2 ). Details of the rate calculations are given in Section  4.1 . The largest enhancement factor of the rates with respect to relaxation by single stars is given in the row ‘Maximum enhancement’. Only a small selection of the parameters is shown here; the complete list of the disruption rates is given in Table  4 . Bottom part: the expected number N ⋆ subscript 𝑁 ⋆ N_{\star} of HVSs and S-stars (models ‘M2’ and ‘M3’). The bottom row shows the observed number (for the S-stars; e.g., Genzel et al. 2010 ) and the estimated number (for HVSs; e.g., Brown 2015 ). The complete list of numbers is given in Table  5 .

			$\Gamma /\left({10}^{-6}{\mathrm{yr}}^{-1}\right)$
			Single	Binary
				M1	M2	M3	M4
By stars			24.5	6.8	0.130	0.031	0.002
By stars and spirals
$i$	$N_s$	${\tilde{\Sigma }}_{s,0,i}$
(deg)		$(M_{\odot }{\mathrm{pc}}^{-2}$)
0.0	1	500.0	29.5	504.6	0.518	0.227	0.027
9.5	10	1000.0	30.4	673.8	0.636	0.298	0.038
90.0	1	500.0	26.3	251.7	0.269	0.100	0.011
180.0	1	500.0	28.2	470.6	0.384	0.160	0.021
Maximum enhancement			1.2	99.1	4.9	9.6	24.6
					$N_{\star }=\Gamma t_{\star }$
By stars					23	1
By stars and spirals
0.0	1	500.0			93	8
9.5	10	1000.0			114	11
90.0	1	500.0			48	4
180.0	1	500.0			69	6
Observed/estimated number					$\sim 300$	$\sim 20$

Table 4. Table 4 Disruption rates Γ Γ \Gamma of single stars (‘Single’) and binary stars (‘Binary’) due to relaxation driven by stars (first data row), and driven by nuclear spiral arms (in addition to stars; all other data rows). In the latter case, each row corresponds to a different inclination i 𝑖 i (expressed in degrees), number of spiral arm events N s subscript 𝑁 s N_{\mathrm{s}} , and (constant) spiral arm surface density Σ ~ s , 0 , i subscript ~ Σ s 0 𝑖 \tilde{\Sigma}_{\mathrm{s},0,i} (in units of M ⊙ ​ pc − 2 subscript M direct-product superscript pc 2 \mathrm{M}_{\odot}\,\mathrm{pc^{-2}} ). The left (right) table applies to prograde (retrograde) orbits. For binary stars, four models, ‘M1’ through ‘M4’, are assumed (cf. Table  2 ). Details of the rate calculations are given in Section  4.1 . Bottom rows: the largest enhancement factor of the rates with respect to relaxation by single stars.

Prograde
			$\Gamma /\left({10}^{-6}{\mathrm{yr}}^{-1}\right)$
			Single	Binary
				M1	M2	M3	M4
By stars			24.5	6.8	0.130	0.031	0.002
By stars and spirals
$i$	$N_s$	${\tilde{\Sigma }}_{s,0,i}$
(deg)		$(M_{\odot }{\mathrm{pc}}^{-2}$)
0.0	1	100.0	25.5	259.3	0.237	0.091	0.011
0.0	1	500.0	29.5	504.6	0.518	0.227	0.027
0.0	1	1000.0	30.3	459.8	0.604	0.264	0.029
0.0	10	100.0	26.1	358.6	0.295	0.123	0.016
0.0	10	500.0	29.4	384.8	0.498	0.211	0.023
0.0	10	1000.0	29.8	356.5	0.511	0.214	0.023
9.5	1	100.0	26.0	181.4	0.228	0.077	0.008
9.5	1	500.0	28.0	457.4	0.377	0.159	0.020
9.5	1	1000.0	29.5	601.1	0.516	0.231	0.030
9.5	10	100.0	27.6	359.4	0.323	0.126	0.015
9.5	10	500.0	29.7	636.2	0.540	0.246	0.033
9.5	10	1000.0	30.4	673.8	0.636	0.298	0.038
19.0	1	100.0	26.0	164.8	0.223	0.074	0.007
19.0	1	500.0	28.2	446.9	0.387	0.161	0.020
19.0	1	1000.0	29.5	583.0	0.517	0.231	0.030
19.0	10	100.0	27.7	311.5	0.329	0.125	0.014
19.0	10	500.0	29.7	568.2	0.550	0.244	0.030
19.0	10	1000.0	30.1	610.2	0.592	0.271	0.034
28.5	1	100.0	25.8	146.9	0.211	0.069	0.007
28.5	1	500.0	28.1	405.4	0.382	0.157	0.019
28.5	1	1000.0	29.5	519.6	0.500	0.215	0.026
28.5	10	100.0	27.1	260.9	0.293	0.108	0.012
28.5	10	500.0	29.3	494.0	0.492	0.210	0.025
28.5	10	1000.0	30.0	535.2	0.554	0.244	0.029
38.0	1	100.0	25.5	127.3	0.196	0.062	0.006
38.0	1	500.0	27.7	361.1	0.347	0.139	0.016
38.0	1	1000.0	29.1	474.0	0.461	0.197	0.024
38.0	10	100.0	26.8	218.4	0.276	0.099	0.010
38.0	10	500.0	29.1	455.6	0.468	0.198	0.023
38.0	10	1000.0	29.4	486.6	0.497	0.216	0.026
47.5	1	100.0	25.3	110.8	0.183	0.057	0.005
47.5	1	500.0	27.5	315.2	0.325	0.126	0.014
47.5	1	1000.0	28.7	401.8	0.421	0.174	0.020
47.5	10	100.0	26.9	212.0	0.273	0.096	0.010
47.5	10	500.0	28.9	422.1	0.444	0.184	0.021
47.5	10	1000.0	29.3	450.0	0.487	0.208	0.024
57.0	1	100.0	25.2	100.8	0.176	0.054	0.005
57.0	1	500.0	27.1	291.1	0.313	0.121	0.013
57.0	1	1000.0	28.4	404.1	0.402	0.167	0.019
57.0	10	100.0	26.0	184.3	0.225	0.078	0.008
57.0	10	500.0	28.8	385.0	0.439	0.178	0.019
57.0	10	1000.0	29.7	449.3	0.506	0.212	0.024
66.5	1	100.0	25.0	95.6	0.170	0.051	0.005
66.5	1	500.0	26.9	259.9	0.296	0.112	0.012
66.5	1	1000.0	28.2	383.4	0.384	0.157	0.018
66.5	10	100.0	25.8	183.0	0.228	0.079	0.008
66.5	10	500.0	29.2	415.3	0.470	0.194	0.021
66.5	10	1000.0	29.5	401.1	0.482	0.199	0.022
76.0	1	100.0	24.9	85.1	0.163	0.048	0.004
76.0	1	500.0	26.6	245.6	0.277	0.102	0.011
76.0	1	1000.0	28.1	362.8	0.370	0.151	0.017
76.0	10	100.0	25.3	151.6	0.205	0.071	0.007
76.0	10	500.0	28.6	372.0	0.401	0.160	0.018
76.0	10	1000.0	29.6	453.8	0.488	0.205	0.023
85.5	1	100.0	24.9	79.4	0.162	0.047	0.004
85.5	1	500.0	26.3	245.1	0.268	0.098	0.011
85.5	1	1000.0	28.2	381.8	0.383	0.156	0.018
85.5	10	100.0	25.5	155.3	0.223	0.077	0.007
85.5	10	500.0	27.8	380.1	0.389	0.162	0.018
85.5	10	1000.0	29.5	438.6	0.492	0.205	0.023
Max. enhancement			1.2	99.1	4.9	9.6	24.6

Table 5. Table 5 Similar to Table  4 , here showing the expected number of HVSs (model ‘M2’) and S-stars (model ‘M3’).

Prograde
			$N_{\star }=\Gamma t_{\star }$
			Binary
			M2	M3
By stars			23	1
By stars and spirals
$i$	$N_s$	${\tilde{\Sigma }}_{s,0,i}$
(deg)		$(M_{\odot }{\mathrm{pc}}^{-2}$)
0.0	1	100.0	43	3
0.0	1	500.0	93	8
0.0	1	1000.0	108	10
0.0	10	100.0	53	5
0.0	10	500.0	89	8
0.0	10	1000.0	92	8
9.5	1	100.0	41	3
9.5	1	500.0	68	6
9.5	1	1000.0	93	9
9.5	10	100.0	58	5
9.5	10	500.0	97	9
9.5	10	1000.0	114	11
19.0	1	100.0	40	3
19.0	1	500.0	69	6
19.0	1	1000.0	93	9
19.0	10	100.0	59	5
19.0	10	500.0	99	9
19.0	10	1000.0	106	10
28.5	1	100.0	38	3
28.5	1	500.0	69	6
28.5	1	1000.0	90	8
28.5	10	100.0	52	4
28.5	10	500.0	88	8
28.5	10	1000.0	99	9
38.0	1	100.0	35	2
38.0	1	500.0	62	5
38.0	1	1000.0	83	7
38.0	10	100.0	49	4
38.0	10	500.0	84	7
38.0	10	1000.0	89	8
47.5	1	100.0	33	2
47.5	1	500.0	58	5
47.5	1	1000.0	76	6
47.5	10	100.0	49	4
47.5	10	500.0	80	7
47.5	10	1000.0	87	8
57.0	1	100.0	32	2
57.0	1	500.0	56	4
57.0	1	1000.0	72	6
57.0	10	100.0	40	3
57.0	10	500.0	79	7
57.0	10	1000.0	91	8
66.5	1	100.0	30	2
66.5	1	500.0	53	4
66.5	1	1000.0	69	6
66.5	10	100.0	41	3
66.5	10	500.0	84	7
66.5	10	1000.0	86	7
76.0	1	100.0	29	2
76.0	1	500.0	50	4
76.0	1	1000.0	66	6
76.0	10	100.0	37	3
76.0	10	500.0	72	6
76.0	10	1000.0	88	8
85.5	1	100.0	29	2
85.5	1	500.0	48	4
85.5	1	1000.0	69	6
85.5	10	100.0	40	3
85.5	10	500.0	70	6
85.5	10	1000.0	88	8

Equations

$${\Phi }_{\bullet }(R,\varphi ,z)=-\frac{G{M}_{\bullet }}{\sqrt{R^2+z^2+{ϵ}_{\mathrm{soft}}^2}},$$

$${\Phi }_{\mathrm{bulge}}(R,\varphi ,z)={\Phi }_0{\left(1+\frac{a}{a_0}\right)}^{2+\beta }$$

$$a={\left(\frac{R^2{\cos }^2(\varphi )}{{\gamma }_x^2}+\frac{R^2{\sin }^2(\varphi )}{{\gamma }_y^2}+\frac{z^2}{{\gamma }_z^2}\right)}^{1/2}.$$

$$v_{\mathrm{c}}^2=G{M}_{\bullet }/R.$$

$$v_{\mathrm{c}}^2$$

$$\approx {\Phi }_0(2+\alpha )(R/a_0{)}^{2+\beta },$$

$${\Sigma }_{\mathrm{s}}(R,\varphi ,t)={\tilde{\Sigma }}_{\mathrm{s}}(R,t)\exp \left[\mathrm{i}m\varphi +\mathrm{i}f(R,t)\right]$$

$$m\varphi +f(R,t)=\mathrm{constant}$$

$$k(R,t)\equiv \frac{\partial f(R,t)}{\partial R},$$

$$\cot (\alpha )\equiv \left|R\frac{\partial \varphi }{\partial R}\right|=\left|\frac{kR}{m}\right|,$$

$${\Sigma }_{\mathrm{s}}(R,\varphi ,t)\approx {\Sigma }_{\mathrm{a}}(t)\exp \left[\mathrm{i}k(R_0,t)(R-R_0)\right],$$

$${\Sigma }_{\mathrm{a}}(t)={\tilde{\Sigma }}_{\mathrm{s}}(R_0,t)\exp \left[\mathrm{i}m{\varphi }_0+\mathrm{i}f(R_0,t)\right].$$

$${\Phi }_{\mathrm{plane}\mathrm{wave}}$$

$$\times \exp \left[-|k(R_0,t)z|\right].$$

$${\Phi }_{\mathrm{s}}(R,\varphi ,z,t)$$

$$\times \exp \left[-|k(R,t)z|\right].$$

$$|k(R,t)|=\frac{m\cot (\alpha )}{R}.$$

$$f(R,t)=m\cot (\alpha )\ln (R)+f_0(t),$$

$$f(R,t)=m\cot (\alpha )\ln (R/R_0)-m{\Omega }_{\mathrm{s}}t.$$

$${\Phi }_{\mathrm{s}}(R,\varphi ,z,t)\approx -2\pi G{\tilde{\Sigma }}_{\mathrm{s}}(R,t)\frac{R\tan (\alpha )}{m}$$

$$\times \exp \left[\mathrm{i}m\left\{\varphi +\cot (\alpha )\ln (R/R_0)-{\Omega }_{\mathrm{s}}t\right\}\right]$$

$$\times \exp \left[-|z|\frac{m\cot (\alpha )}{R}\right].$$

$${\tilde{\Sigma }}_{\mathrm{s}}(R,t)=\sum _{i=0}^{N_{\mathrm{s}}}{\tilde{\Sigma }}_{\mathrm{s},0,i}(R)\exp \left[-\frac{(t-t_{\mathrm{s},i}{)}^2}{2{\sigma }_{\mathrm{s},i}^2}\right].$$

$$\ddot{R}$$

$$\ddot{\varphi }$$

$$\ddot{z}$$

$$J_{\mathrm{c}}=r_{\mathrm{c}}{v}_{\mathrm{c}},$$

$$v_{\mathrm{c}}^2/r_{\mathrm{c}}=|\nabla \Phi {|}_{r=r_{\mathrm{c}}},$$

$$2[\psi (r_{\mathrm{c}})-\mathcal{E}]-J_{\mathrm{c}}^2/r_{\mathrm{c}}^2=2[\psi (r_{\mathrm{c}})-\mathcal{E}]-r_{\mathrm{c}}|\nabla \Phi {|}_{r=r_{\mathrm{c}}}=0.$$

$$\cos i=\frac{\mathbf{J}\cdot \hat{\mathbf{z}}}{|\mathbf{J}\cdot \hat{\mathbf{z}}|}=\mathrm{sign}\left(\dot{\varphi }\right)\frac{R}{r}.$$

$$\Sigma \sim 1.36\Delta n_e{m}_{\mathrm{p}}0.25R\approx 42{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2}\left(\frac{R}{100\mathrm{pc}}\right).$$

$$\Sigma$$

$$\sim 1.3\times 10^3{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2}{\left(\frac{R}{100\mathrm{pc}}\right)}^{-0.9},$$

$${\mathcal{D}}_{\mathcal{E}}$$

Full text

1. Introduction

The Galactic Center (GC) is well known for harboring Sgr A*, a supermassive black hole (MBH) of $M_{\bullet}=4\times 10^{6}\,\mathrm{M}_{\odot}$ (Eisenhauer et al., 2005; Ghez et al., 2008; Gillessen et al., 2009; Boehle et al., 2016). The GC region lies in the central region of the Galactic Bulge ($\sim 0.3-3$ kpc from the GC), which consists mainly of a population of old ($\sim 10\,\mathrm{Gyr}$) stars. The inner 300 pc are known as the Nuclear Bulge (e.g., Mezger et al. 1996), which, in addition to old stars, also contains young stars in dense environments (Philipp et al., 1999; Mezger et al., 1999; Figer et al., 2004). In particular, there are two massive clusters, the young ($\sim 5$ Myr old) Quintuplet cluster ($\sim 10^{4}\,\mathrm{M}_{\odot}$) at $\sim 30$ pc from the center (e.g., Figer et al. 1999a, b) and the even younger Arches cluster ($\sim 10^{4}\,\mathrm{M}_{\odot}$) at $\sim 30$ pc from the center (e.g., Figer et al. 1999a). Furthermore, the central parsec contains, in addition to old stars, a population of several hundred young O and B stars (e.g., Schödel et al. 2014; see Alexander 2005; Genzel et al. 2010 for reviews).

These young stellar populations are indicative of recent star formation several million years ago, and this is consistent with the presence of a large supply of molecular gas in the inner $\sim 200\,\mathrm{pc}$ (the central molecular zone; Serabyn & Morris 1996). On scales of several hundred pc, gas can be perturbed by torques from nuclear bars (e.g., Shlosman et al. 1989; Englmaier & Shlosman 2004), driving inflows to the central regions (Englmaier & Shlosman, 2000). Simulations show that these inflows can result in transient features in the gas density, closely resembling spiral arms (Englmaier & Shlosman, 2000; Maciejewski et al., 2002; Maciejewski, 2004b; Ann & Thakur, 2005; Namekata et al., 2009; Thakur et al., 2009; Kim et al., 2012; Li et al., 2015; Ridley et al., 2017). These nuclear spiral arms are indeed observed in 50 to 80 per cent of both active and quiescent galaxies (Laine et al., 1999; Regan & Mulchaey, 1999; Martini & Pogge, 1999; Pogge & Martini, 2002; Martini et al., 2003a, b). In particular, for NGC 1097 (Fathi et al., 2006; Davies et al., 2009; van de Ven & Fathi, 2010) and NGC 6951 (Storchi-Bergmann et al., 2007), streaming motions along the nuclear spiral arms have been mapped.

Nuclear spiral arms can be considered as small-scale variants of galactic-scale spiral arms. Whereas galactic-scale spiral arms consist of both stars and gas, nuclear spiral arms are observed to consist of gas only (e.g., Fathi et al. 2006; Davies et al. 2009; van de Ven & Fathi 2010; Storchi-Bergmann et al. 2007). Galactic-scale spiral arms can perturb stars orbiting in the Galaxy (Goldreich & Tremaine 1979; Binney & Tremaine 2008, S6.2.6). If the spiral-structure features are transient, then these perturbations can drive heating of stars in the Galactic disk, i.e., increasing the velocity dispersion in the radial and transverse directions, while not affecting velocities in the vertical direction (Barbanis & Woltjer, 1967). This has been invoked (Sellwood & Carlberg, 1984; De Simone et al., 2004) to explain the observed age-velocity-dispersion relation in the Galactic disk and the ratio of the vertical-to-radial velocity dispersion, although a combination with heating by giant molecular clouds (GMCs; Spitzer & Schwarzschild 1951, 1953) is likely required to best match the observations (Carlberg, 1987; Jenkins & Binney, 1990; Jenkins, 1992).

Similarly, nuclear-scale spiral arms can perturb the orbits of stars in certain regions in the central few hundred pc around the MBH. In particular, nuclear spiral structures can drive stars onto orbits bringing them very close to the MBH, triggering strong interactions.

These interactions include the disruption of a single star by the tidal force of the MBH resulting in a luminous outburst known as a tidal disruption event (TDE; Hills 1975), or the tidal disruption of a stellar binary (Hills, 1988). In the latter case, one of the stars in the binary remains bound to the MBH in a close and eccentric orbit, whereas the other star is ejected with high speeds of $\sim 1000\,\mathrm{km\,s^{-1}}$. This scenario can explain both the S-stars in the GC, orbiting the MBH in close and eccentric orbits (Genzel et al., 2003; Eisenhauer et al., 2005; Ghez et al., 2008; Gillessen et al., 2009), and hypervelocity stars (HVSs), stars observed to be moving with unusually high velocities and no longer bound to the Galaxy (e.g., Brown et al. 2005; see Brown 2015 for a review). HVSs can also be produced by the disruption of clusters of stars such as globular clusters (Arca-Sedda et al., 2016; Capuzzo-Dolcetta & Fragione, 2015) and young clusters like the Arches cluster (Fragione et al., 2017).

A pertinent issue in the disruptions of single stars and binaries by the MBH is how these objects are (continuously) injected to loss cone orbits of the MBH, i.e., to orbits with small pericenter distances (see, e.g., Merritt 2013). Such perturbations can come from other stars; the resulting disruption rates for single stars in earlier studies were found to be $\sim 10^{-5}\,\mathrm{yr^{-1}}$ (Syer & Ulmer, 1999; Magorrian & Tremaine, 1999), whereas rates of $\sim 10^{-4}\,\mathrm{yr^{-1}}$ are found in more recent calibrated studies (Wang & Merritt, 2004; Stone & Metzger, 2016). For stellar binaries, the rates arising from relaxation by stars are too low to explain the observed number of S-stars and HVSs (Yu & Tremaine, 2003). However, relaxation can also be driven by massive structures in the GC such as GMCs and massive stellar clusters, which can enhance the rate to be consistent with the observed numbers (Perets et al., 2007; Perets, 2009).

Given the likely presence of nuclear spiral arms in the GC and in galactic centers in general, this motivates an investigation of the possibility that nuclear spiral structures within the regions of the central few hundred pc around the central MBHs of galaxies can drive stars and binaries onto loss-cone orbits. This provides a channel, in addition to GMCs and massive stellar clusters, to supply the MBH with stars and binaries on close orbits, and hence produce TDEs, S-star-like stars and HVSs.

In addition to enhancing event rates, binary disruptions driven by nuclear spiral arms may result in different signatures of HVSs that can help pinpoint their origin. Plausibly, nuclear spiral arm structures are highly flattened and lie within the plane of the Galaxy. This suggests that relaxation driven by nuclear spiral arms is a strong function of the inclination of the orbit of the (barycenter of the) stellar binary with respect to the Galactic disk, and the latter should be reflected in the inclination distribution of both the captured and ejected stars. For the captured stars, relaxation by the surrounding stellar cluster rapidly randomizes the orbital orientation and modifies the eccentricity distribution to be less eccentric; this is consistent with the S-stars (Perets et al., 2009; Antonini & Merritt, 2013; Antonini, 2014; Hamers et al., 2014). The ejected stars, on the other hand, do not experience orientation randomization, and the inclination dependence of relaxation of binaries by nuclear spiral arms should be reflected in the HVS inclination distribution. The latter is indeed not completely consistent with a random orientation with respect to the GC (Brown et al., 2014).

In this paper, we investigate relaxation driven by nuclear spiral arm structures, and consider the implications for TDEs, S-star-like stars and HVSs. In Sections  2 and  3, we carry out test-particle integrations in an assumed potential of the inner region of the GC and transient nuclear spiral arms, and investigate diffusion driven by the latter. These results are used in Section 4 to compute the resulting disruption rates for both single stars and binaries. We discuss our results in Section 5, and conclude in Section 6.

2. Orbits around the GC in the presence of nuclear spiral arms

Our methodology consists of two main parts. First, we integrate orbits in the GC environment considering also the contribution to the potential due to nuclear spiral arms (Section 2), and we extract relaxation time-scales and diffusion coefficients from these integrations (Section 3). Subsequently, we use our results to compute the disruption rates by the MBH for both single and binary stars (Section 4).

2.1. Potential

We model the potential of the inner several hundred pc of the GC with the following components. Throughout Section 2, we adopt a cylindrical coordinate system $(R,\phi,z)$ with the GC as the origin; the distance to the GC is given by $r\equiv\sqrt{R^{2}+z^{2}}$.

2.1.1 MBH

For the MBH, we adopt the potential

[TABLE]

where $G$ is the gravitational constant, and with an MBH mass $M_{\bullet}=4\times 10^{6}\,\mathrm{M}_{\odot}$ adopted from the MBH in the Milky Way GC (Eisenhauer et al., 2005; Ghez et al., 2008; Gillessen et al., 2009; Boehle et al., 2016). For numerical reasons, the MBH potential is softened with a (conservative) softening length of $\epsilon_{\mathrm{soft}}=10^{-3}\,\mathrm{pc}$. Our aim is not to accurately model the dynamics in the innermost region around the MBH (i.e., collisional relaxation by stars near the sphere of influence), but our interest is in larger regions extending to several hundred pc around the MBH. Therefore, such softening is justified. Note that the potential equation (1), equivalent to a Plummer potential (Plummer, 1911), gives rise to steady precession of the orbit around the MBH, and does not lead to diffusion of the orbital angular momentum nor energy.

2.1.2 Bulge

The bulge is modeled with a potential similar to that of Kim & Morris (2001), i.e.,

[TABLE]

with

[TABLE]

Here, we adopt $\beta=-1.8$, $a_{0}=1\,\mathrm{pc}$ and $\Phi_{0}=3\times 10^{4}\,\mathrm{pc^{2}\,Myr^{-2}}$ from Kim & Morris (2001).

In contrast to Kim & Morris 2001, who assumed a triaxial and rotating bulge potential with $\gamma_{y}=0.9$ and $\gamma_{z}=0.8$ in most of their models, we assume a spherically symmetric potential, i.e., we set $\gamma_{x}=\gamma_{y}=\gamma_{z}=1$. In the latter case, orbits experience precession only, and the energy and angular momentum are conserved. Any deviations from energy and angular-momentum conservation can then be ascribed to perturbations from nuclear spiral arms. Triaxial potentials, although interesting on their own, are beyond the scope of this paper.

In Fig. 1, we show the circular-speed curve $v_{\mathrm{c}}(R)$ from the potential $\Phi_{\bullet}+\Phi_{\mathrm{bulge}}$ as a function of $R$ (with $\phi=0$ and $z=0$, and setting $\gamma_{x}=\gamma_{y}=\gamma_{z}=1$). The inner few pc are dominated by the MBH (dotted lines), in which case

[TABLE]

The outer regions are dominated by the bulge, for which, in our model,

[TABLE]

where the last line applies if $R\gg a_{0}$, i.e., in that case, $v_{c}\propto R^{1+\beta/2}=R^{0.1}$. The circular speed in our model is $\sim 100\,\mathrm{km\,s^{-1}}$ at $\sim$ 10 pc from the center, and rises to $\sim 150\,\mathrm{km\,s^{-1}}$ at $\sim$ 1 kpc. This is consistent with, e.g., the potential Model I in S2.7 of Binney & Tremaine (2008). At larger radii, our model clearly does not give a realistic representation of the Galaxy potential: the circular-speed curve of the latter flattens around $220\,\mathrm{km\,s^{-1}}$ (see, e.g., fig. 2.20 of Binney & Tremaine 2008), whereas our model gives an unbounded increase of $v_{\mathrm{c}}$. Here, we are only interested in the central few hundred pc around the MBH; therefore, this inconsistency should not affect our conclusions.

2.1.3 Nuclear spiral arms

General potential

Let $\Sigma_{\mathrm{s}}(R,\phi,t)$ be the surface density associated with spiral structure at the coordinate $(R,\phi)$ at time $t$. We assume that the spiral structure is located within an infinitely thin disk $z=0$. Following Section 6.2.2b of Binney & Tremaine (2008), we write $\Sigma_{\mathrm{s}}$ in a form separating the variations in density along a spiral arm, described by $\tilde{\Sigma}_{\mathrm{s}}(R,t)$, and the variations experienced while passing between arms, i.e.,

[TABLE]

(see also, e.g., Maciejewski 2004a). Here, the surface density has been written in complex form ($\mathrm{i}\equiv\sqrt{-1}$). The number of spiral arms is given by $m>0$, and the shape function $f(R,t)$ describes the spiral shape. Along a spiral arm,

[TABLE]

(e.g., Section 6.1.3c of Binney & Tremaine 2008), and, therefore, the radial surface-density dependence along the spiral arms is described by $\tilde{\Sigma}_{\mathrm{s}}(R,t)$. The spiral wave number $k$ is defined from the shape function as

[TABLE]

and is related to the winding angle $\alpha$ according to

[TABLE]

where we used equation (7) after the last equality.

To find the potential at any point $(R,\phi,z)$ and time $t$ due to the surface density $\Sigma_{\mathrm{s}}(R,\phi,t)$, we follow Section 6.2.2b of Binney & Tremaine (2008) and expand $\Sigma_{\mathrm{s}}(R,\phi,t)$ around a local point $(R_{0},\phi_{0})$, neglecting any variations with respect to $\phi$ and variations in $\tilde{\Sigma}_{\mathrm{s}}(R,t)$, i.e., $f(R,t)\approx f(R_{0},t)+k(R_{0},t)(R-R_{0})$ and $\tilde{\Sigma}_{\mathrm{s}}(R,t)\approx\tilde{\Sigma}_{\mathrm{s}}(R_{0},t)$. This expansion is appropriate for tightly-wound spiral arms ($|k|R\gg 1$) because the short wavelengths, $\lambda=2\pi/|k|\ll R$, imply that only the surface density close to a point $(R_{0},\phi_{0})$ contributes to the potential, whereas the contribution from the more distant surface density cancels out. Therefore,

[TABLE]

with

[TABLE]

Equation (10) represents a plane wave in a razor-thin homogenous disk with the wave vector $\bf{k}$ oriented radially. The potential due to this distribution is (Binney & Tremaine, 2008, S5.6.1)

[TABLE]

Here, $G$ is the gravitational constant. Setting $R_{0}=R$ and $\phi_{0}=\phi$, this gives the approximate potential due to surface density equation (6),

[TABLE]

The error made in equation (2.1.3) is $\mathcal{O}(|kR|^{-1})$ (Binney & Tremaine, 2008, S6.2.2b); for tightly-wound spirals, $|kR|^{-1}\ll 1$. Also, in equation (2.1.3) we assumed that $\tilde{\Sigma}_{\mathrm{s}}(R,t)$ is a slowly varying function of $R$.

Logarithmic spirals

Nuclear spiral arms are observed to be approximately logarithmic (e.g., van de Ven & Fathi 2010). For a logarithmic spiral, $R\propto\exp(\phi)$. From equation (7), this implies that the shape function $f\propto\ln(R)$. Furthermore, from equation (8), $k\propto 1/R$ implying that $\alpha$ is constant (cf. equation 9). Therefore, the wave number can be written as

[TABLE]

Integrating equation (8), the shape function is given in terms of the constant $\alpha$ and $m$ by

[TABLE]

where $f_{0}(t)$ is a function of time only. We define the radius $R_{0}$ such that the logarithmic term in equation (15) vanishes (i.e., the value of $R_{0}$ only determines the overall phase), and introduce a pattern speed $\Omega_{\mathrm{s}}$ at which the spiral pattern rotates. Equation (15) can then be written as

[TABLE]

Substituting the expressions for $k$ and $f$ into equation (2.1.3) we find, for logarithmic spirals,

[TABLE]

The $z$-scale height of equation (2.1.3) is $H_{z}=R\tan(\alpha)/m$, i.e., equation (2.1.3) implies a constant $H_{z}/R=\tan(\alpha)/m$. For logarithmic spirals, the fractional error in the potential is $|kR|^{-1}=\tan(\alpha)/m$. For example, for $m=2$ spirals with $\alpha=10^{\circ}$, $\tan(\alpha)/m\approx 0.088$.

For illustration, we show in Fig. 2 isosurface density curves for logarithmic spirals by combining equations (6) and (16). We assume that $\tilde{\Sigma}_{\mathrm{s}}$ is a constant, and plot curves at different times $t$ for which $\Sigma_{\mathrm{s}}=\tilde{\Sigma}_{\mathrm{s}}$, setting $m=2$, $\alpha=10^{\circ}$, $R_{0}=10\,\mathrm{pc}$ and $\Omega_{\mathrm{s}}=50\,\mathrm{km\,s^{-1}\,kpc^{-1}}$. Note that $R_{0}$ only sets the phase of the spiral arms; the arms have no intrinsic scale.

Transient spirals

The spiral surface density as described above is assumed to last for a certain amount of time, i.e., the nuclear spiral arms are transient and occur in episodes, referred to as nuclear spiral arm events. The total surface density due to nuclear spirals at a given time $t$ is given by summing over all events, i.e.,

[TABLE]

Here, $i$ sums over the spiral arm events, and $\tilde{\Sigma}_{\mathrm{s,0},i}(R)$, $t_{\mathrm{s},i}$, and $\sigma_{\mathrm{s},i}$ are the radial surface density profile, central time and approximate duration, respectively, of the $i^{\mathrm{th}}$ spiral arm event.

In Fig. 3, we show how gravitationally relevant nuclear spiral arms are compared to the MBH and other stars, by plotting the associated mass, $\pi R^{2}\tilde{\Sigma}_{\mathrm{s},0,i}$, as a function of $R$. Here, we compute the stellar mass from the stellar density implied by the bulge potential (Section 2.1.2), assuming a (single) stellar mass of $M_{\star}=1\,\mathrm{M}_{\odot}$. We include three values of $\tilde{\Sigma}_{\mathrm{s},0,i}$ which we adopt in the numerical integrations of Section 3 below (see Section 3.1 for a motivation of these values). The mass fraction in nuclear spirals arms compared to stars is $\sim 10^{-3}$ to $\sim 10^{-2}$ at $R=10\,\mathrm{pc}$ depending on the surface density; at $R=10^{2}\,\mathrm{pc}$, the mass fractions are $\sim 10^{-2}$ to $10^{-1}$. This already shows that nuclear spiral arms are important at relatively large distances from the MBH ($R\gtrsim 10\,\mathrm{pc}$), which we also find in the integrations in Section 3.

2.2. Orbit integrations

We integrate the equations of motion in cylindrical coordinates for the position vector $\boldsymbol{r}=R\hat{\boldsymbol{R}}+z\hat{\boldsymbol{z}}$, $\boldsymbol{\ddot{r}}=-\boldsymbol{\nabla}\Phi$, in the total potential $\Phi=\Phi_{\bullet}+\Phi_{\mathrm{bulge}}+\Phi_{\mathrm{s}}$, i.e.,

[TABLE]

where dots denote derivatives with respect to time. Here, when modeling the motion of a binary in the GC, we neglect the quadrupole moment of the stellar binary, i.e., we treat the binary as a point mass, and $\boldsymbol{r}$ should be interpreted as the binary barycenter position.

Equations (19) can be cast in a system of ordinary differential equations (ODEs), which we solve numerically using odeint from the Python Scipy library. The latter is an interface to the LSODA Fortran routine, which uses variable time-steps to integrate the system of ODEs. The routine automatically and dynamically detects between stiff and nonstiff methods; for stiff cases it uses the backward differentiation formula method (with a dense or banded Jacobian), and for nonstiff cases it uses the Adams method. Error control within the solver is determined by the input parameters rtol and atol, such that the error in each ODE variable $y_{i}$ is less than or equal to $\textsc{rtol}\times\mathrm{abs}(y_{i})+\textsc{atol}$. We set the relative and absolute tolerance parameters to $\textsc{rtol}=10^{-9}$ and $\textsc{atol}=10^{-9}$, respectively.

In the case without nuclear spiral arms and a spherical (non-rotating) bulge, orbits experience precession only, and the energy and angular momentum are conserved. In Fig. 4, we show relative energy errors (solid lines) and angular momentum errors (dashed lines) as a function of time for an orbit initially at $(R,\phi,z)=(10\,\mathrm{pc},0,0)$. The integration time is 100 Myr, which corresponds to $\approx 455$ orbits (the orbital period is $\approx 0.22\,\mathrm{Myr}$). Different line widths correspond to different initial values of $J/J_{\mathrm{c}}$, where $J$ is the orbital angular momentum with respect to the origin, and $J_{\mathrm{c}}$ is the angular momentum of a circular orbit (see also Section 3). As expected, the lower $J/J_{\mathrm{c}}$ (the more eccentric the orbit), the larger the relative energy and angular-momentum errors. Note that the highly eccentric orbits experience rapid precession (with the argument of periapsis changing of the of order of $\sim 90^{\circ}$ per orbit).

After 100 Myr, which is also the time-span used in the integrations in Section 3, the energy errors remain $\lesssim 3\times 10^{-6}$, and the angular-momentum errors remain $\lesssim 10^{-4}$. We are interested in the statistical properties for an ensemble of orbits. A typical change in $J/J_{\mathrm{c}}$ in the simulations of Section 3 is $\sim 0.1$, which is $\gg 10^{-4}$. Therefore, we believe these errors are sufficiently small for our purposes.

2.3. Example evolution

We give a number of examples of orbit integrations to illustrate the ability of transient nuclear spiral structures to diffuse the orbits of objects (single stars and binaries) in regions around the central few hundred pc of the GC. In Figs  5, 6 and 7, we plot the energy (normalized to the initial energy), angular momentum (normalized to $J_{\mathrm{c}}$, cf. Section 3) and orbital inclination, respectively, for four orbits at different initial radii and initial orbital inclinations. The orbital inclination is measured with respect to the $z$ axis; the nuclear spiral arms are assumed to be confined to the plane $z=0$ (cf. Section 2.1.3). In the four panels, we take two values of the initial distance to the center, $r_{\mathrm{init}}=5\,\mathrm{pc}$ and $r_{\mathrm{init}}=30\,\mathrm{pc}$, and two values of the initial inclination, $i_{\mathrm{init}}=0^{\circ}$ and $i_{\mathrm{init}}=72^{\circ}$.

For the nuclear spiral arm parameters, we assume a constant spiral gas surface density, i.e., $\tilde{\Sigma}_{\mathrm{s,0},i}$ is independent of radius, and take $\tilde{\Sigma}_{\mathrm{s,0},i}=500\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$. The number of spiral arm events is $N_{\mathrm{s}}=10$, with $t_{\mathrm{s},i}$ at fixed equal intervals between 0 and 100 Myr, and $\sigma_{\mathrm{s},i}=10\,\mathrm{Myr}$ (i.e., 10 spiral arm events, each lasting about 10 Myr). For each spiral arm event, the spiral pattern speed is set to $50\,\mathrm{km\,s^{-1}\,pc^{-1}}$, $m=2$ arms are assumed, and the winding angle is $\alpha=10^{\circ}$. We refer to Section 3.1 for motivation of these parameter choices.

Close to the MBH (cf. the top panels in Figs  5, 6), the potential is dominated by the MBH, and spiral structure leads to only small variations in the energy and angular momentum. Further away from the MBH (cf. the bottom panels in Figs  5 and 6), nuclear spiral arms can have a large effect on both energy and angular momentum, with the largest variations occurring in $J/J_{\mathrm{c}}$. There is also a strong dependence on the initial inclination: for $r_{\mathrm{init}}=30\,\mathrm{pc}$, $J/J_{\mathrm{c}}$ changes drastically over time if the orbit lies in the same plane as the spiral arms ($i_{\mathrm{init}}=0^{\circ}$), whereas the changes are much smaller if $i_{\mathrm{init}}=72^{\circ}$. In the latter case, $E$ and $J$ change substantially only when the orbit crosses the $z=0$ plane and spiral arms happen to be active near this moment.

Fig. 7 shows how the inclination changes over time. If the orbit is initially coplanar, it remains coplanar. In the highly inclined case, the inclination is affected by the nuclear spiral arms, with $i$ changing ‘impulsively’ at each passage of the Galactic plane, increasing and decreasing in a periodic fashion. In our example, the inclination changes remain modest, with overall changes of less than $\approx 10^{\circ}$.

3. Diffusion by nuclear spiral arms

Adopting the potential model described in Section 2, we carry out a large number (1,200,000) of orbital integrations representing single stars and binaries in the GC (Section 3.1). From these integrations, we determine diffusion coefficients (Section 3.3), which are then used to compute disruption rates in Section 4.

Throughout, we use the orbital angular momentum $J\equiv|\boldsymbol{r}\times\dot{\boldsymbol{r}}|$ and the (negative) energy $\mathcal{E}\equiv-\frac{1}{2}\dot{\boldsymbol{r}}^{2}+\psi(r)$ to describe the relaxation process. Here, $\psi(r)\equiv-\Phi(r)$ is the negative of the potential; note that $\mathcal{E}$ is defined with a minus sign, contrary to convention. The orbital angular momentum is normalized to the orbital angular momentum of a circular orbit, $J_{\mathrm{c}}$. The latter is determined from the potential via

[TABLE]

where $v_{\mathrm{c}}$ is the speed of a circular orbit at radius $r_{\mathrm{c}}$, i.e.,

[TABLE]

and $r_{\mathrm{c}}$ is determined implicitly by setting the radial speed $v_{\mathrm{r}}$ to zero (cf. Section 4.1), i.e.,

[TABLE]

3.1. Initial conditions

An overview of the initial conditions is given in Table 1. We determine energy and angular-momentum relaxation time-scales and diffusion coefficients on a grid of three parameters: the inclination $i$, the number of spiral arm events $N_{\mathrm{s}}$, and the nuclear spiral arm gas surface density $\tilde{\Sigma}_{\mathrm{s},0,i}$. For the latter, we assume, for simplicity, a surface density that is independent of radius $R$. Also, for simplicity, the initial orbits are realized with $\dot{R}_{\mathrm{init}}=\dot{z}_{\mathrm{init}}=0$. In this case, the orbital inclination $i$ is given by

[TABLE]

Also, note that the nodal orientation of the orbit with respect to the plane $z=0$ is described by $\phi$, which is sampled randomly (see below).

For each grid parameter combination $(i,N_{\mathrm{s}},\tilde{\Sigma}_{\mathrm{s},0,i})$, we take 20 values of the initial distance to the center $r_{\mathrm{init}}$ ranging between 10 and 1000 pc with logarithmic spacing, and 5 values of $\dot{\phi}_{\mathrm{init}}$, the initial $\phi$ time derivative. The latter values are taken to be $(0.1-1)\,v_{\mathrm{c,init}}/r_{\mathrm{init}}$, with five values sampled linearly; these values correspond to different values of the initial angular momentum. For each resulting parameter combination, i.e., $(i,N_{\mathrm{s}},\tilde{\Sigma}_{\mathrm{s},0,i},r_{\mathrm{init}},\dot{\phi}_{\mathrm{init}})$, we carry out $N_{\mathrm{MC}}=100$ integrations for a duration of $t_{\mathrm{int}}=100\,\mathrm{Myr}$, each with a different initial phase angle $\phi$ sampled randomly in the range $[0,2\pi)$. This approach results in a total of 1,200,000 integrations.

The nuclear spiral arm parameters are currently not well constrained. For the active galaxy Arp 102B, Fathi et al. (2011) found a two-armed spiral structure (i.e., $m=2$) in the inner kpc. Measurements of the velocity field in NGC 1097 (van de Ven & Fathi, 2010) indicate a logarithmic spiral with two arms and a pitch angle of $\alpha=52^{\circ}\pm 4^{\circ}$. The spiral overdensity in electrons was found to be $\Delta n_{e}\sim 50\,\mathrm{cm^{-3}}$ at radii of several hundred pc (cf. Fig. 4 of van de Ven & Fathi 2010), and the radius-to-height ratio at similar radii was inferred to be $h/R\sim 0.25$. The electron overdensity corresponds to a gas overdensity of $\Delta\rho\sim 1.36\,m_{\mathrm{p}}\Delta n_{e}$, where $m_{\mathrm{p}}$ is the proton mass, and the factor 1.36 takes into account the presence of helium (van de Ven & Fathi, 2010). This gives a gas surface overdensity of $\Sigma\sim\Delta\rho\,h$, where $h$ is the height associated with the nuclear spiral arm gas, i.e.,

[TABLE]

An estimated theoretical upper limit on $\Sigma$ can be obtained by requiring stability with respect to gravitational collapse, i.e., by setting the Toomre $Q$ parameter to unity (Toomre, 1964). Approximately, this implies $Q=c_{\mathrm{s}}\kappa/(\pi G\Sigma)=1$, where $c_{\mathrm{s}}$ is the sound speed, here taken to be a constant $c_{\mathrm{s}}=10\,\mathrm{km\,s^{-1}}$, and $\kappa$ is the epicyclic frequency, here computed as a function of radius $R$ using Eq. (3.79a) from Binney & Tremaine 2008 with the potential $\Phi(R)=\Phi_{\bullet}(R)+\Phi_{\mathrm{bulge}}(R)$ and setting $z=0$, $\epsilon_{\mathrm{soft}}=0$, and $\gamma_{x}=\gamma_{y}=\gamma_{z}=1$ (cf. Sections 2.1.1 and 2.1.2). The surface density associated with $Q=1$ is then given by

[TABLE]

where the last line applies if $R\gg 1\,\mathrm{pc}$, and we substituted the values from Section 2.1.

In the simulations of Englmaier & Shlosman (2000), the pitch angle was found to be a function of radius, with $\alpha\sim 10^{\circ}-20^{\circ}$ between $\sim 200$ and 400 pc from the center. In other simulations by Kim et al. (2012), typical spiral arm surface densities were $\sim 100\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$, and the spiral arm structures lasted on the order of 100 Myr. In the recent simulations of Ridley et al. (2017), two-armed spirals were found in their GC models, with surface densities of $\sim 100\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$.

Here, we take a constant $m=2$ and a pitch angle of $\alpha=10^{\circ}$. Note that the approximation of tightly-wound spirals, made in the potential (cf. Section 2.1.3), breaks down for large $\alpha$. For the surface density $\tilde{\Sigma}_{\mathrm{s},0,i}$, we take three values, 100, 500 and 1000 $\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$. We consider the lower value as typical, whereas the upper value is close to the maximum value for gravitational stability (cf. equation 3.1), and should be interpreted as an extreme case. The number of spiral arm events is taken to be either $N_{\mathrm{s}}=1$ or 10. For $N_{\mathrm{s}}=1$, we assume there is one spiral arm event at $t=50\,\mathrm{Myr}$, with a characteristic width of $\sigma_{\mathrm{s},i}=50\,\mathrm{Myr}$. For $N_{\mathrm{s}}=10$, the spiral event times are set between 10 and 90 Myr with linear spacing, and $\sigma_{\mathrm{s},i}=10\,\mathrm{Myr}$. The spiral arm pattern speed $\Omega_{\mathrm{s}}$ is set to $\Omega_{\mathrm{s}}=50\,\mathrm{km\,s^{-1}\,kpc^{-1}}$.

3.2. Data reduction

From the integrations, we extract time series of the energy $\mathcal{E}$ and orbital angular momentum $J$. The latter is normalized to the angular momentum of a circular orbit, $J_{\mathrm{c}}$ (cf. equation 20). For each realization $i$, we record the energy and angular-momentum changes, $\Delta\mathcal{E}_{i}$ and $\Delta[J/J_{\mathrm{c}}(\mathcal{E})]_{i}$, respectively, at the end of the integration time with respect to the initial time, i.e., for a time interval of $\Delta t=t_{\mathrm{int}}=100\,\mathrm{Myr}$. Subsequently, the energy and angular-momentum diffusion coefficients are computed from the rms changes according to

[TABLE]

The underlying assumption in equations (26) is that diffusion in both energy and angular momentum occurs in a random-walk process, i.e., $\mathcal{E}$ and $J/J_{\mathrm{c}}$ grow proportionally to $\Delta t^{1/2}$.

In determining the diffusion coefficients, we exclude integrations in which the test particle becomes essentially unbound from the GC. We consider this to be the case if the final distance from the center $r_{\mathrm{final}}$ (i.e., after 100 Myr) satisfies $r_{\mathrm{final}}>1000\,\mathrm{pc}$. This can occur if the spiral arm perturbations are strong, particularly for the highest surface density considered, $\tilde{\Sigma}_{\mathrm{s},0,i}=1000\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$. In the case of an absence of sufficient data points among the 100 Monte-Carlo integrations ($N\leq 40$), the diffusion coefficients are set to zero, i.e., no relaxation.

3.3. Results

Although the diffusion coefficients are used in the rate calculations in Section 4, a (perhaps) physically more intuitive quantity is the relaxation time-scale, $t_{\mathrm{r},i}\equiv 1/\mathcal{D}_{i}$, where $i$ is either $\mathcal{E}$ or $J$. In this section, results are shown in terms of the relaxation time-scales. In Fig. 8, the angular-momentum relaxation time-scales $t_{\mathrm{r},J}$ are shown as a function of $r_{\mathrm{init}}$. In each panel, the different symbols correspond to different inclinations, restricting to prograde inclinations (retrograde inclinations are included in Fig. 9). Each panel corresponds to a different combination of $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. The time-scales are averaged over the initial $\dot{\phi}$, or, equivalently, the initial $J/J_{\mathrm{c}}$. Similarly, in Fig. 10, the energy relaxation time-scales are plotted, in this case for a single value of $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. Typically, the energy relaxation time-scales are roughly an order-of-magnitude longer than the angular-momentum time-scales. Therefore, in the following, we will consider angular-momentum relaxation only.

In addition, we show in Fig. 11 the fractions of unbound orbits in the simulations as a function of the initial distance to the GC. Here, we define an orbit to be unbound if $r>10^{3}\,\mathrm{pc}$ and the orbital speed is larger than the local escape speed at the end of the integration. The fractions are computed over the $N_{\mathrm{MC}}=100$ realizations, and averaged over the initial $J/J_{\mathrm{c}}$.

Generally, the following trends are revealed in Figures 8 and 11.

Relaxation by spirals is inefficient close to the center ($r\lesssim 10\,\mathrm{pc}$), where the MBH and bulge potentials dominate. At larger distances($r\sim 300\mathrm{pc}$), $t_{\mathrm{r},J}$ can be as short as $\sim 500\,\mathrm{Myr}$. 2. 2.

Relaxation by spirals increases in efficiency for larger $\tilde{\Sigma}_{\mathrm{s},0,i}$ (larger overall spiral gas surface density). In addition, relaxation is more effective for a larger number of spiral arm events ($N_{\mathrm{s}}=10$ vs. 1) with a shorter typical duration (10 vs. 50 Myr). However, as the strength of the spiral perturbations increases, the probability for becoming unbound increases as well. In fact, for all values of $\tilde{\Sigma}_{\mathrm{s},0,i}$ considered here, particles become unbound for $r_{\mathrm{init}}>r_{\mathrm{init,crit}}$, where the critical initial radius depends on $i$, $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. For this reason, larger spiral perturbations do not necessarily lead to higher disruption rates (for the latter, we require the orbit to remain bound). This effect will be quantified in Section 4. 3. 3.

Spirals are very efficient at strongly perturbing orbits at zero inclination $(i=0^{\circ}$), in which case many test particles become unbound at typically several tens of pc. At larger inclinations, perturbations are weaker and particles remain bound at larger radii.

4. Disruption rates

Having computed diffusion coefficients associated with nuclear spiral arms in Section 3, we here proceed to compute the disruption rates by the MBH. First, for reference, we briefly review some aspects of loss-cone theory (see also, e.g., Merritt 2013).

Below, we assume that objects (i.e., single stars or binaries) are disrupted when they pass within the tidal disruption or loss cone radius,

[TABLE]

For single stars, $M$ is the stellar mass and $R$ the stellar radius, whereas for binaries, $M=M_{\mathrm{bin}}=M_{1}+M_{2}$ is the binary mass and $R$ is the binary semimajor axis $a_{\mathrm{bin}}$. Any object initially on an orbit with a pericenter distance less than $r_{\mathrm{lc}}$ is rapidly depleted (within the orbital time $P$); therefore, the loss cone needs to be refilled in order to drive steady disruption. For reference, $r_{\mathrm{lc}}$ is plotted as a function of $a_{\mathrm{bin}}$ with the black lines in Fig. 12 for several values of $M_{\mathrm{bin}}$.

Close to the MBH, the orbital time-scale is short but relaxation is typically inefficient (this applies both to relaxation by stars, and by nuclear spiral arms). This implies that the loss-cone orbits are depleted, and refilling is driven by the slow process of relaxation (empty loss-cone regime). Far away from the MBH, relaxation (by stars) occurs more efficiently, but the orbital time-scale is longer as well. Rather than being determined by relaxation, loss cone refilling is limited by the orbital time-scale (full loss-cone regime). Therefore, the disruption rate in the full loss-cone regime cannot be enhanced by relaxation processes such as nuclear spiral arms. Any enhancement must originate from the empty loss-cone regime, where refilling is limited by relaxation. Typically, the flux of objects into the loss cone peaks near the transition radius between the empty and full loss cones.

In the case of the GC and assuming that relaxation is driven by stars, the loss-cone flux for single stars peaks around the sphere of influence of the MBH, i.e., at a few pc. As shown in Section 3, this is inside the regime where nuclear spiral arms are effective. Therefore, for single stars, nuclear spiral arms can only increase the relaxation rate in the full loss-cone regime, and therefore do not lead to significantly higher disruption rates. For binaries, $r_{\mathrm{lc}}$ can be much larger, depending on the binary semimajor axis. This shifts the transition from the empty to full-loss-cone regimes to much larger radii. Consequently, if there is enhanced relaxation by nuclear spiral arms (compared to relaxation by single stars) within the empty loss cone regime, then the disruption rates can be substantially increased. The aim of Section 4 is to calculate by how much the rates are increased, using the numerical integrations of Section 3.

4.1. Methodology

The disruption rates are computed using standard loss-cone theory, and assuming spherical symmetry. The latter is not self-consistent with the integrations in Section 3, in which the assumed geometry was nonspherical. Rather than modifying the standard loss-cone theory to non-spherical geometries, we here choose an approximate approach in which the rates are computed assuming spherical symmetry, taking different values of the relaxation rate by nuclear spiral arms depending on the (initial) inclination. Despite this inconsistency in the geometry, we expect that our approach still captures the most important aspects of the dependence of the disruption rates on the inclination. A self-consistent, anisotropoic loss-cone calculation is left for future work.

Let $\mathcal{E}\equiv-\frac{1}{2}v^{2}+\psi(r)$ be the (negative) energy, where $\psi(r)$ is the negative of the potential. In our disruption rate calculations, we include in the potential the contributions from the MBH and the spherical bulge potential, i.e.,

[TABLE]

with $\epsilon_{\mathrm{soft}}=0$ and $\gamma_{x}=\gamma_{y}=\gamma_{z}=1$ (cf. Section 2.1). The density $\rho(r)$ of single stars is computed from $\Phi_{\mathrm{bulge}}$ using the Poisson equation, and we assume the same mass for all background stars, $M_{\star}=1\,\mathrm{M}_{\odot}$, i.e., the number density $n(r)=\rho(r)/M_{\star}$.

We adopt the Cohn-Kulsrud boundary layer formalism (Cohn & Kulsrud, 1978) which describes the angular-momentum flux into the loss cone, i.e., $F_{\mathrm{lc}}(\mathcal{E})$, the number of stars lost per unit time and energy $\mathcal{E}$. The formalism is based on matching the non-averaged and orbit-averaged solutions to the Fokker-Planck equation in angular-momentum space. The resulting steady-state angular-momentum flux into the loss cone is

[TABLE]

Here, $f(\mathcal{E})$ is the distribution function in energy space of the scattered population (single stars or binaries). It is computed from the number density $n(r)$ using Eddington’s formula (Eddington, 1916), i.e.,

[TABLE]

The (radial) orbital period is given by

[TABLE]

where $v_{\mathrm{r}}(r,\mathcal{E},\mathcal{R})=\{2[\psi(r)-\mathcal{E}]-\mathcal{R}\,J_{\mathrm{c}}^{2}(\mathcal{E})/r^{2}\}^{1/2}$ is the radial orbital speed as a function of the angular-momentum variable $\mathcal{R}\equiv(J/J_{\mathrm{c}})^{2}$. In equation (29), the dependence of the distribution function of the scattered population in angular-momentum space (i.e., $\mathcal{R}$), is assumed to be logarithmic with $\mathcal{R}$ as $\mathcal{R}\rightarrow 0$.

The quantity $\bar{\mu}(\mathcal{E})$ is the orbit-averaged angular-momentum diffusion coefficient that describes the diffusion rate in angular-momentum space. For diffusion by stars, it can be computed explicitly from the distribution function $f(\mathcal{E})$ and the potential $\psi(r)$; for completeness, the expression for $\bar{\mu}_{\mathrm{stars}}$ is given in Appendix A.

In the case of relaxation by nuclear spiral arms, the associated $\bar{\mu}(\mathcal{E})$ is computed from the numerical integrations of Section 3, which yielded $\mathcal{D}_{J}(r)$, the angular-momentum diffusion coefficient due to spiral arms as a function of distance to the GC. As mentioned in Section 3.3, diffusion in energy space is much less efficient than diffusion in angular momentum, therefore, we neglect the former. We orbit average $\mathcal{D}_{J}(r)$ over radius according to (e.g., Merritt 2013, 5.5.2)

[TABLE]

where $v(r,\mathcal{E})=\sqrt{2[\psi(r)-\mathcal{E}]}$ is the orbital speed at radius $r$ for an orbit with energy $\mathcal{E}$, and $p(\mathcal{E})$ is the phase-space volume element per unit energy111The factor 4 in equations (32) and (33) is adopted from Merritt (2013); it is omitted by some authors, e.g., Spitzer (1987).,

[TABLE]

Subsequently, $\bar{\mu}(\mathcal{E})$ is computed from

[TABLE]

i.e., we include both relaxation by stars and spirals; any enhancement of the disruption rate by stars is due to nuclear spiral arms.

For illustration, we show in Fig. 13 the non-averaged $\mathcal{D}_{J}(r)$ as a function of $r$ (black dotted lines), and the averaged $\overline{\mathcal{D}}_{J}(\mathcal{E})$ as a function of $r_{\mathrm{c}}=r_{\mathrm{c}}(\mathcal{E})$ (solid red lines), for three different inclinations, and fixed values of $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. We recall that $r_{\mathrm{c}}(\mathcal{E})$ is the radius of a circular orbit (cf. equation 22). The orbit-averaging procedure tends to smooth out sharp kinks in the non-averaged diffusion coefficients.

The factor in equation (29) depending on $\mathcal{R}_{0}(\mathcal{E})$ takes into account the empty and full loss-cone regimes. It is given by

[TABLE]

where $\mathcal{R}_{\mathrm{lc}}(\mathcal{E})$ is the value of $\mathcal{R}$ corresponding to the loss cone; it can be found by setting $v_{\mathrm{r}}(r,\mathcal{E},\mathcal{R}_{\mathrm{lc}})=0$ (the orbit just grazing the loss cone), i.e.,

[TABLE]

The function $q(\mathcal{E})$, defined as

[TABLE]

characterizes the full and empty loss-cone regimes, and the function $\xi$ of $q$ is given by

[TABLE]

Here, $\alpha_{m}$ is the $m^{\mathrm{th}}$ zero of the Bessel function of the first kind, $J_{0}(\alpha)$. For $q\ll 1$, $\xi(q)\approx(2/\sqrt{\pi})\sqrt{q}$, whereas for $q\gg 1$, $\xi(q)\approx 1$.

For illustration, we briefly discuss the two limits of $q$ in equation (29). In the empty loss-cone regime ($q\ll 1$, or, more precisely, $q\ll-\ln\mathcal{R}_{\mathrm{lc}}$), $F_{\mathrm{lc}}\approx 4\pi^{2}P_{\mathrm{r}}(\mathcal{E})J_{\mathrm{c}}^{2}(\mathcal{E})\bar{\mu}(\mathcal{E})f(\mathcal{E})(-\ln\mathcal{R}_{\mathrm{lc}})^{-1}$. Using $p(\mathcal{E})\approx J_{\mathrm{c}}^{2}(\mathcal{E})P_{\mathrm{r}}(\mathcal{E})$ and $N(\mathcal{E})=4\pi^{2}p(\mathcal{E})f(\mathcal{E})$ (Merritt, 2013, 3.2.1), this can be written as

[TABLE]

where $N(\mathcal{E})$ is the number of stars per unit energy, and $\overline{t}_{\mathrm{r}}(\mathcal{E})\equiv 1/\bar{\mu}(\mathcal{E})$ is the relaxation time-scale.

In the full loss-cone regime ($q\gg 1$, or, more precisely, $q\gg-\ln\mathcal{R}_{\mathrm{lc}}$), $F_{\mathrm{lc}}\approx 4\pi^{2}J_{\mathrm{lc}}^{2}(\mathcal{E})f(\mathcal{E})$, which can be written as

[TABLE]

Equations (39) and (40) are the well-known relations for the empty and full loss-cone fluxes (e.g., equations 6.8 and 6.9 of Alexander 2005).

From the loss-cone flux computed using the Cohn-Kulsrud formalism (equation 29), we compute the total disruption rate by integrating over all energies and binary semimajor axes $a_{\mathrm{bin}}$ (in the case of binaries), i.e.,

[TABLE]

Here, the function $g_{\mathrm{bin}}(\mathcal{E},a_{\mathrm{bin}})$ takes into account the binary occurrence rate, the binary semimajor axis distribution, and the process of binary evaporation (with a similar approach as Perets et al. 2007). For single stars, formally $g_{\mathrm{bin}}(\mathcal{E},a_{\mathrm{bin}})=\delta(a_{\mathrm{bin}}-r_{\mathrm{lc}})$, where $\delta$ is the delta function, and we set $r_{\mathrm{lc}}=1\,\mathrm{AU}$. For binaries,

[TABLE]

where $f_{\mathrm{PDMF}}$ is the fraction of stars with respect to the assumed background stellar distribution $n(r)$ for the mass range of a binary model given the present-day mass function (PDMF; see Section 4.2 below), $f_{\mathrm{bin}}$ is the fraction of binaries with respect to $n(r)$, $\mathrm{d}N/\mathrm{d}a_{\mathrm{bin}}$ is the initial normalized binary semimajor axis distribution, $t_{\mathrm{H}}=10\,\mathrm{Gyr}$ is an (approximate) Hubble time, $t_{\star}$ is the stellar lifetime (depending on the binary component masses), and $\overline{t}_{\mathrm{evap}}(\mathcal{E},a_{\mathrm{bin}})$ is an orbit-averaged evaporation time-scale. The latter time-scale captures the effect of perturbations of binaries by the stellar background distribution – soft binaries tend to become softer due to this process, eventually unbinding them (Heggie, 1975; Hut & Bahcall, 1983; Hut, 1983). The local evaporation time-scale is computed from (Binney & Tremaine, 2008, S7.5.7)

[TABLE]

Here, the velocity dispersion $\sigma(r)$ is computed from the isotropic Jeans equation,

[TABLE]

where the enclosed mass is

[TABLE]

The binary Coulomb factor $\Lambda_{\mathrm{bin}}$ is taken to be

[TABLE]

For hard binaries ($\Lambda_{\mathrm{bin}}<1$), we assume $t_{\mathrm{evap}}=10^{8}\,\mathrm{Gyr}$, i.e., essentially no evaporation; we neglect the process of binary hardening. Subsequently, the local evaporation time-scale equation (43) is orbit averaged (cf. equation 32, with the angular-momentum diffusion coefficient $\mathcal{D}_{J}$ replaced by the evaporation time-scale), yielding $\overline{t}_{\mathrm{evap}}(\mathcal{E},a_{\mathrm{bin}})$. The limits of the integration over $a_{\mathrm{bin}}$ are discussed below, in Section 4.2.

In practice, the calculations are carried out on a grid in energy; the energies in this grid are computed according to $\mathcal{E}=\psi(r)$ from a logarithmic grid in radius with $5\times 10^{-2}\,\mathrm{pc}<r<10^{3}\,\mathrm{pc}$. Integrals over radius are carried out over the range $10^{-8}\,\mathrm{pc}<r<10^{10}\,\mathrm{pc}$.

4.2. Binary star models

We adopt four models representing different types of binary stars. An overview is given in Table 2. The parameters included in each model are the PDMF fraction $f_{\mathrm{PDMF}}$, the total binary mass $M_{\mathrm{bin}}$, the semimajor axis/orbital period distribution, the fraction $f_{\mathrm{bin}}$ of binaries with respect to the assumed background stellar distribution $n(r)$, and the stellar main-sequence lifetime $t_{\star}$. For simplicity, we assume equal-mass-ratio binaries, and we use a ‘typical’ mass, $M_{\mathrm{bin}}$, to compute the evaporation time-scales (cf. equation 43) and the limits on the integration over binary separations (see below), whereas we assume a range of primary star masses when computing $f_{\mathrm{PDMF}}$.

The background distribution $n(r)$ is assumed to consist of $1\,\mathrm{M}_{\odot}$ stars, but the binaries considered here are more massive, and, therefore, less common. We take the relative frequency of the stars of the binaries into account with the fraction $f_{\mathrm{PDMF}}$. This fraction is computed assuming a multicomponent power law Miller-Scalo initial mass function (IMF; Miller & Scalo 1979) which is converted to the PDMF assuming a constant star-formation history (Figer et al. 2004, see also Genzel et al. 2010) over $t_{\mathrm{H}}=10\,\mathrm{Gyr}$, and using the main-sequence lifetime $t_{\star}$. For a given $M_{\mathrm{bin}}$, we compute $t_{\star}$, the main-sequence life-time of a star with mass $M_{\mathrm{bin}}/2$ and metallicity $z=0.02$, using SSE (Hurley et al., 2000) within the AMUSE framework (Portegies Zwart et al., 2013; Pelupessy et al., 2013).

The relative frequency of stars with respect to the $1\,\mathrm{M}_{\odot}$ population is then given by

[TABLE]

where $M_{\mathrm{1,low}}$ and $M_{\mathrm{1,up}}$ are the adopted lower and upper ranges on the mass of the primary star of the binary (depending on the binary model), $M_{\star}=1\,\mathrm{M}_{\odot}$, $M_{\star,\mathrm{up}}=125\,\mathrm{M}_{\odot}$, and $h(M)=t_{\star}(M)/t_{\mathrm{H}}$ if $t_{\star}(M)<t_{\mathrm{H}}$, and $h(M)=1$ otherwise. The masses $M_{\mathrm{1,low}}$ and $M_{\mathrm{1,up}}$ and the fractions $f_{\mathrm{PDMF}}$ are included in Table 2.

The binary models include Solar-type binaries, binaries capable of producing HVSs and S-stars, and O-star binaries. For Solar-type stars, we adopt a lognormal period distribution with mean $\mu_{\log_{10}(P/\mathrm{d})}=5$ and standard deviation $\sigma_{\log_{10}(P/\mathrm{d})}=2.3$ (Duquennoy & Mayor 1991; Raghavan et al. 2010; Moe & Di Stefano 2016, hereafter MDS16). For the more massive binaries, we adopt the analytic fits from MDS16 (cf. equations 20-23 from that paper) for the distribution of $\log_{10}(P/\mathrm{d})$ as a function of primary mass $M_{1}=M_{\mathrm{bin}}/2$. The binary fractions $f_{\mathrm{bin}}$ are adopted from visual inspection of Fig. 37 of MS16.

The following limits for $a_{\mathrm{bin}}$ are adopted in equation (41). The lower limit $a_{\mathrm{bin,min}}$ is computed assuming an orbital period of $\log_{10}(P/\mathrm{d})=0.2$, the lowest value considered in Moe & Di Stefano (2016). With the exception of our HVS and S-star models, the upper limit $a_{\mathrm{bin,max}}$ is taken to be $a_{\mathrm{bin,max}}=50\,\mathrm{AU}$, approximately the largest semimajor axis for which, after tidal disruption by the MBH, the bound star orbits the MBH within several pc. The semimajor axis of the bound orbit is approximately given by (Hills, 1988)

[TABLE]

Here, we assume $M_{2}=M_{\mathrm{bin}}/2$. For reference, we plot equation (48) for our assumed binary models with the blue lines in Fig. 12. For $a_{\mathrm{bin}}=50\,\mathrm{AU}$, $a_{\mathrm{bound}}$ is on the order of 1 pc (with some dependence on $M_{\mathrm{bin}}$).

For our HVS model (cf. model ‘M2’ in Table 2), we impose a smaller $a_{\mathrm{bin,max}}$ because we are only interested in stars that can escape the Galaxy potential. The escape speed is approximately given by (Hills, 1988)

[TABLE]

Assuming an escape speed from the Galactic Bulge of $900\,\mathrm{km\,s^{-1}}$ (Haardt et al., 2006), equation (49) for our model ‘M2’ implies a maximum binary separation of $a_{\mathrm{bin,max}}\approx 0.984\,\mathrm{AU}$. We impose this maximum separation for the HVS model in equation (41).

For our S-star model (cf. model ‘M3’ in Table 2), we impose $a_{\mathrm{bin,max}}=6\,\mathrm{AU}$, corresponding to a maximum semimajor axis of the bound orbit around the MBH of $\approx 0.044\,\mathrm{pc}$ (cf. equation 48 and Fig. 12), approximately consistent with the S-star orbits (Genzel et al., 2003; Eisenhauer et al., 2005; Ghez et al., 2008; Gillessen et al., 2009).

4.3. Results

4.3.1 Differential loss rates

To illustrate the enhancement of the loss rates due to relaxation by nuclear spiral structure in the empty loss-cone regime, we show in Figs. 14 and 15, for two combinations of the grid parameters $(i,N_{\mathrm{s}},\tilde{\Sigma}_{\mathrm{s},0,i})$, the loss-cone flux $F_{\mathrm{lc}}$. For the case of the disruption of binary stars, we here do not yet integrate over the binary semimajor axis distribution although we do take into account the PDMF and the binary fraction (cf. $f_{\mathrm{PDMF}}$ and $f_{\mathrm{bin}}$ in equation 42). Rather than plotting $F_{\mathrm{lc}}(\mathcal{E})$ as a function of energy $\mathcal{E}$, we plot the loss-cone flux as a function of distance to the GC. We define $F_{\mathrm{lc}}(r)$ as the differential flux per unit distance, i.e., $F_{\mathrm{lc}}(\mathcal{E})\,\mathrm{d}\mathcal{E}=F_{\mathrm{lc}}(r)\,\mathrm{d}r$. Also, we set $r=r_{\mathrm{c}}(\mathcal{E})$, the radius of a circular orbit, and multiply $F_{\mathrm{lc}}(r_{\mathrm{c}})$ by $r_{\mathrm{c}}$ to get the loss rate per unit $\ln r_{\mathrm{c}}$ (i.e., $r_{\mathrm{c}}\,F_{\mathrm{lc}}(r_{\mathrm{c}})$ has dimensions of one over time).

The top panel in Fig. 14 applies to the disruption of single stars ($r_{\mathrm{lc}}=1\,\mathrm{AU}$). The solid black (red) lines correspond to relaxation by stars (nuclear spiral arms). For reference, the associated relaxation time-scales $\overline{t}_{\mathrm{r}}(\mathcal{E})=1/\overline{\mathcal{D}}_{J}(\mathcal{E})$ are shown with the black and red dotted lines; in this case, the right-hand axes apply. For relaxation by single stars, the differential rate peaks at several pc, near the radius of influence. In the case of relaxation by nuclear spiral arms, the peak occurs at larger radii, where the relaxation rate is enhanced compared to relaxation by single stars. The relaxation time-scale by spirals is shorter at radii beyond the radius of the peak; however, the latter regime lies within the full loss-cone regime for single stars. Therefore, the rates are not enhanced, and the integrated rate due to spirals only, $\Gamma_{\mathrm{by\,spirals}}\approx 2\times 10^{-6}\,\mathrm{yr^{-1}}$, is actually lower in this example (by an order of magnitude) compared to the integrated rate due to stellar relaxation, $\Gamma_{\mathrm{by\,stars}}\approx 2.5\times 10^{-5}\,\mathrm{yr^{-1}}$ (these rates are indicated in the top panel in Fig. 14; note that the rates given below in Section 4.3.2 always include relaxation by stars).

The situation is different for binaries. In the bottom panel of Fig. 14, we consider the disruption of binaries (binary model ‘M3, i.e., representing S-stars) with a semimajor axis of $\approx 3.93\,\mathrm{AU}$, for which the loss-cone radius is $r_{\mathrm{lc}}\approx 1.2\times 10^{-3}\,\mathrm{pc}$ (cf. equation 27 or Fig. 12). For relaxation by single stars, the integrated rate over energies (indicated in the panels; note that these rates are not integrated over $a_{\mathrm{bin}}$) would be larger compared to the case of disruptions of single stars if the PDMF were ignored. In Fig. 14, the PDMF is taken into account, however; consequently, the integrated binary disruption rate is lower compared to the single star disruption rate. Nuclear spirals do significantly enhance the binary disruption rate compared to single stars. Comparing relaxation by nuclear spirals to relaxation by stars, the binary disruption rate is higher by a factor of $\approx 4.4$.

For the disruption of binaries and radii $r_{\mathrm{c}}\gtrsim 30\,\mathrm{pc}$, spirals sufficiently enhance the relaxation rate to drive relaxation into the full loss-cone regime. This is clear from the solid red line in the bottom panel of Fig. 14, which approaches the black dashed line (representing the full lone-cone regime, cf. equation 40) for $r_{\mathrm{c}}\gtrsim 30\,\mathrm{pc}$. In contrast, the loss-cone flux for relaxation by stars (black solid line in the bottom panel) does not reach the full loss-cone value for the radial range shown.

A similar figure is shown in Fig. 15, now for a larger, nonzero inclination. The spiral relaxation time-scales at intermediate radii ($\sim 30\,\mathrm{pc}$) are longer compared to the coplanar case, but at large radii ($\gtrsim 50\,\mathrm{pc}$), the loss cone is still full. The integrated rates are lower compared to the case of $i=0^{\circ}$ in Fig. 14.

4.3.2 Integrated loss rates

For each grid parameter $(i,N_{\mathrm{s}},\tilde{\Sigma}_{\mathrm{s},0,i})$, the loss-cone flux, equation (29), is integrated over energy, and, for binaries, over the binary semimajor axis taking into account the PDMF, the semimajor axis distribution, the finite stellar lifetime, and binary evaporation (cf. equation 41). An abbreviated list of the rates is given in the top part of Table 3; a complete list is given in Table 4 in Appendix B. We include results taking into account relaxation by stars only (first data row), and taking into account relaxation driven by nuclear spiral arms in addition to stars (all other data rows). For the disruption of binary stars, the four assumed models are included in the right four columns.

The highest binary disruption rate obtained is $\Gamma\approx 6\times 10^{-4}\,\mathrm{yr^{-1}}$ for binary model ‘M1’ with $i=9.5^{\circ}$, $N_{\mathrm{s}}=10$, and $\tilde{\Sigma}_{\mathrm{s},0,i}=1000\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$. This rate is a factor $\approx 100$ higher compared to relaxation by stars only, the largest enhancement compared to relaxation by single stars. For the other binary disruption models, nuclear spirals can enhance the rates by factors of $\sim$ 5, 10 and 25 for ‘M2’, ‘M3’ and ‘M4’, respectively (cf. the bottom rows in Table 3).

Although less significantly, nuclear spiral arms also enhance the disruption rate of single stars. Our highest TDE rate is $\approx 3.0\times 10^{-5}\,\mathrm{yr^{-1}}$, which is a factor of 1.2 higher compared to relaxation by stars only.

In Fig. 16, the disruption rates by spirals and stars are plotted as a function of $\cos\,i$, assuming binary model ‘M2’ (representing HVSs). Different colors and symbols correspond to different values of $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. There is a clear dependence of $\Gamma$ on $i$, with the overall normalization depending on $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$. For large $N_{\mathrm{s}}$ and/or $\tilde{\Sigma}_{\mathrm{s},0,i}$, the rates for coplanar binaries ($i=0^{\circ}$) are lower compared to the rates for slightly inclined binaries ($i\approx 10^{\circ}$). For more inclined binaries, the rates decline again. This dependence on $i$ be understood by noting that for $i=0^{\circ}$, many binaries become unbound due to spiral perturbations (cf. Section 3.3). For larger $i$, the unbinding effect becomes less dominant, whereas the diffusion rate also decreases, resulting in a maximum of the disruption rate. If $N_{\mathrm{s}}$ and/or $\tilde{\Sigma}_{\mathrm{s},0,i}$ are small, the overall strength of the spiral perturbations is smaller. Therefore, orbits can remain bound if $i=0^{\circ}$.

Apart the above intricacy of the $i$-dependence for different $N_{\mathrm{s}}$ and $\tilde{\Sigma}_{\mathrm{s},0,i}$, there is generally a clear trend of decreasing $\Gamma$ with less inclined orbits. Nevertheless, even for high inclinations, the rates are still significantly higher (by factors of a few) compared to relaxation by stars only (in which $\Gamma\approx 1.3\times 10^{-7}\,\mathrm{yr^{-1}}$, cf. Table 3 and the horizontal black dashed line in Fig. 16). The disruption rates are approximately symmetric around $\cos\,i=0$.

4.3.3 Expected number of HVSs and S-stars

In the bottom part of Table 3, we give the expected number of HVSs and S-stars computed from $N_{\star}=\Gamma\,t_{\star}$, where the lifetime $t_{\star}$ is taken from Table 2. A complete list is given in Table 5 in Appendix B. The number of HVSs in our models ranges between $\sim 50$ to $\sim 100$; the estimated number of HVSs within 100 kpc of the Galaxy is $\sim 300$ (Brown, 2015), consistent within a factor of a few. Our predicted S-star numbers are a few, up to $\sim 10$, consistent within a factor of a few with the observed number of $\sim 20$ (e.g., Genzel et al. 2010).

5. Discussion

5.1. Overall rates

In our calculations, the disruption rate of single stars driven by relaxation by stars is $\approx 2.5\times 10^{-5}\,\mathrm{yr^{-1}}$, which is consistent with previous studies in which rates of $\sim 10^{-5}\,\mathrm{yr^{-1}}$ were found (Syer & Ulmer, 1999; Magorrian & Tremaine, 1999). In our models with nuclear spiral arms included as well, the disruption rate of single stars is at most $\approx 3.0\times 10^{-5}\,\mathrm{yr^{-1}}$, which is $20\%$ higher compared to relaxation by single stars only. Although the enhancement by nuclear spiral arms is not as significant as for binaries, it complements other models for enhanced TDE rates such as triaxial nuclei (Vasiliev & Merritt, 2013).

For our model ‘M1’ (Solar-type stars), the binary disruption rate driven by stars is $\approx 0.7\times 10^{-5}\,\mathrm{yr^{-1}}$, which is lower compared to the rate of Yu & Tremaine 2003, who found, for $1\,\mathrm{M}_{\odot}$ stars, a rate of $10^{-5}(\eta/0.1)\,\mathrm{yr}^{-1}=4\times 10^{-5}\,\mathrm{yr^{-1}}$, where $\eta$ is the binary fraction, which we took to be 0.4 (cf. Table 2). For more massive binaries ($M_{\mathrm{bin}}=4\,\mathrm{M}_{\odot}$ and higher), we found rates of at most $10^{-7}\,\mathrm{yr^{-1}}$. Zhang et al. (2013) found rates of $\sim 10^{-4}$ to $10^{-5}\,\mathrm{yr^{-1}}$, for various binary injection models and stellar masses between 3 and 15 $\mathrm{M}_{\odot}$. The much higher rates found by the latter authors could be explained by noting that Zhang et al. (2013) assumed an IMF and not the PDMF. In addition, Zhang et al. (2013) also considered top-heavy IMFs.

For our HVS and S-star model binaries, our disruption rates for relaxation by nuclear spiral arms are on the order of $10^{-6}\,\mathrm{yr}$, implying numbers of HVSs and S-stars which are roughly consistent with the observed numbers. Assuming relaxation by stars only, the expected numbers are inconsistent with the observed numbers. This result is similar to relaxation by massive perturbers (Perets et al., 2007). Nuclear spiral arms can be considered as a complementary channel to supply binaries to the loss cone of an MBH, with a roughly equal contribution compared to relaxation driven by massive perturbers.

5.2. Implications for HVSs, and caveats

We assumed that the nuclear spiral arm structures are confined to the plane of the Galaxy. Consequently, the disruption rates, in particular for binaries able to produce HVSs (binary model ‘M2’), are dependent on the inclination $i$ of the (original) binary with respect to the Galactic plane; this dependence was shown in Fig. 16. The difference in disruption rates between a nearly coplanar and nearly perpendicular orbit is modest, i.e., a factor of $\sim 2$. Even for nearly perpendicular orbits, the disruption rates are still higher compared to relaxation by stars.

An anisotropic distribution of HVSs was also considered by Lu et al. (2010), who suggested that HVSs originate from the inner clockwise-rotating stellar disk within half a parsec of the MBH in the GC. This implies an angular distribution which is distinct from ours: our inclination distribution is not linked to the clockwise disk but to the Galactic plane, suggesting that HVSs should be less common far above the this plane, compared to within it.

A caveat in our results is that our quoted inclination is the initial orbital inclination of the binary with respect to the Galactic plane. In addition to stochastically changing the orbital angular momentum, transient nuclear spiral arms can also affect the orbital inclination (see, e.g., Fig. 7). To simplify our analysis, a changing inclination was not taken into account.

In Fig. 17, we show, for two values of the parameter $N_{\mathrm{s}}$, the change of the orbital inclination over the integration averaged over the $N_{\mathrm{MC}}=100$ realizations, as a function of the initial inclination $i_{0}$, for initially prograde orbits. The standard deviations of the change of the inclination are shown with the error bars. The data are averaged over the initial orbital angular momentum, $J/J_{\mathrm{c}}$, and results for different initial radii are shown.

Depending on the parameters, $\Delta i$ can be large, even approaching $90^{\circ}$, although the standard deviation can be substantial as well. There is a strong dependence on the initial radius and inclination, with larger changes for larger radii, and smaller initial inclinations. These effects could modify the true inclination dependence of the disruption rates, since objects on initially nearly coplanar can be transferred to highly inclined orbits, for which the disruption rates are lower. Effectively, this could smear out the inclination dependence shown in Fig. 16.

Another caveat, mentioned in Section 4.1, is that spherical symmetry was assumed in the rate computations, which is clearly not self-consistent with the underlying nuclear spiral potential (the latter was derived assuming a razor-thin surface density distribution). Also, the binary was treated as a point mass, i.e., the quadrupole moment of the binary was neglected. These caveats should be addressed in future, more detailed work. Another aspect that merits future investigation is the effect of triaxial potentials.

5.3. Implications for S-stars

Binary model ‘M3’ represents the origin of S-stars (B-type stars with masses of $\sim 8\,\mathrm{M}_{\odot}$). Their disruption rates by nuclear spiral arms are similar to the HVS model. The dependence of the rates on the inclination is the same compared to HVS, apart from an overall different normalization factor. This may appear to be inconsistent with the S-stars, which are observed to be highly isotropic in their orbital distribution (e.g., Eisenhauer et al. 2005; Ghez et al. 2008). However, the time-scale for the S-stars orbits to randomize their orientation, i.e., the ‘vector’ or ‘2d’ resonant relaxation time-scale (e.g., Merritt 2013, eq. 5.237) is

[TABLE]

which is short compared to the typical S-star lifetime of $\sim 40\,\mathrm{Myr}$ (e.g., Table 2). Here, we assumed a semimajor axis of $a=10\,\mathrm{mpc}$, a background stellar mass of $M_{\star}=1\,\mathrm{M}_{\odot}$, and $N=10^{3}$. Therefore, our result of the inclination dependence of the disruption rate is not inconsistent with the observed isotropy of the S-star orbits.

5.4. Binary MBHs

Relaxation by nuclear spiral arms can also be important for binary MBH systems. Stars approaching a binary MBH at a close distance statistically extract energy from the binary MBH orbit, gradually shrinking the latter, whereas the stars are ejected through the slingshot effect (Saslaw et al. 1974). In the standard model, this process halts when the binary MBH reaches a separation of order 1 parsec, giving rise to the ‘final parsec problem’ (see,.e.g, Merritt & Milosavljević 2005 for a review). Among the proposed mechanisms to solve this problem (see, e.g., Vasiliev 2014 for an overview), is enhanced relaxation by massive perturbers (Perets & Alexander, 2008). Similarly, nuclear spiral arms can drive stars into the ‘loss cone’ of the binary MBH, thereby accelerating its coalescence.

Here, we only briefly consider the effects of nuclear spiral arms, and do not consider the evolution of the binary MBH in response to stellar encounters. Also, for simplicity, we assume the same background population of stars as before (cf. Section 4.1), whereas in a more realistic situation, the stellar distribution is dynamic (in particular, it is depleted near the MBH). Note, however, that in our case the enhancement by nuclear spiral arms is at large radii of several hundred parsec, at which stars could be supplied from the Galactic Bulge. Here, we adopt a loss-cone radius of $r_{\mathrm{lc}}=0.5\,\mathrm{pc}$ (approximately the separation of a ‘stalled’ binary MBH).

In Fig. 18, we show the differential rates at which stars are supplied to the binary MBH, assuming relaxation driven by stars (solid black line), and relaxation driven by nuclear spiral arms (solid red line). The integrated rates are $\approx 1.3\times 10^{-5}$ and $\approx 2.7\times 10^{-2}\,\mathrm{yr^{-1}}$ for relaxation by stars and spirals, respectively. This suggests that nuclear spiral arms are extremely efficient at supplying stars to the binary MBH, although it remains to be seen how this rate is affected when binary MBH evolution and the depletion of stars are taken into account.

5.5. Other implications

In principle, the perturbations from nuclear spiral arms act not only on stars and stellar binaries, but also on any object or system in similar regions as the nuclear spiral arms. In particular, the orbits of stellar clusters like the Quintuplet cluster could be excited to high eccentricities, leading to the breakup of clusters and building up nuclear star clusters (Tremaine et al., 1975). Similarly, gas clouds could be excited to near-radial orbits, supplying nuclear star clusters with gas for star formation (e.g., Milosavljević 2004; McLaughlin et al. 2006). In the GC in particular, such processes may have contributed to the formation of the innermost O stars, which are concentrated in a disk-like structure (e.g., Genzel et al. 2000; Levin & Beloborodov 2003; Tanner et al. 2006; Paumard et al. 2006; Lu et al. 2009; Bartko et al. 2009).

Also, planets and planetesimals could be driven to loss-cone orbits, giving rise to their tidal breakup, and relatively low-intensity flares (Zubovas et al. 2012; in the case of rocky planetesimals, the injected energy arises from friction of small tidally disrupted fragments with the surrounding ambient gas). In the case of the GC, relaxation of planetesimals by stars in this model is, however, already sufficient to account for the observed flaring rate of about once per day, irrespective whether the planetesimals formed in a large-scale disk around the MBH, or in disks around stars (Hamers & Portegies Zwart, 2015).

6. Conclusions

We have considered diffusion of single and binary stars into the loss cone of an MBH, driven by perturbations from transient nuclear spiral arm structures which are observed in the centers of galaxies. Disruptions of single stars are believed to be observable as TDEs. Disruptions of binary stars are thought to result in a highly eccentric star on a tight orbit around the MBH, and a star with high velocity escaping from the galaxy (Hills, 1988). The latter mechanism is currently favored to explain the origin of the S-stars in the GC, and HVSs.

Focussing on the GC, we adopted a simple model for the potential associated with nuclear spiral arms, and we carried out test-particle integrations within this potential to model the dynamics of single and binary stars in regions of the central few hundred pc around the GC. From these integrations, we extracted energy and angular-momentum diffusion coefficients and relaxation time-scales due to transient nuclear spiral arms, for different combinations of plausible parameters. Using the angular-momentum diffusion coefficients, we computed disruption rates of single stars and binaries, giving estimates of the formation rates of S-stars and HVSs due nuclear spiral arms in the GC. Our conclusions are as follows.

1. Nuclear spiral arms are ineffective at driving orbital diffusion in the regions of the central $\lesssim 10\,\mathrm{pc}$ around the MBH, where the potential is dominated by the MBH and the stellar bulge. At larger radii, angular-momentum relaxation time-scales can be as short as $\sim 500\,\mathrm{Myr}$ depending on the inclination, the spiral arm surface density $\tilde{\Sigma}_{\mathrm{s},0,i}$, and the number of spiral arm events $N_{\mathrm{s}}$. Relaxation in energy is typically an order-of-magnitude less efficient compared to relaxation in angular momentum.

2. The relaxation rate due to nuclear spiral arms depends on the inclination of the orbit of the binary with respect to the plane of the Galaxy (we assumed that the nuclear spiral arm structures are confined within the plane of the Galaxy). Generally, a higher inclination implies less efficient relaxation. However, if the inclination is low and spiral perturbations are strong, then the binary can become unbound from the GC, in which case the binary can no longer be disrupted by the MBH. Similarly, if the spiral arm surface density $\tilde{\Sigma}_{\mathrm{s},0,i}$ is high and/or the number of spiral arm events $N_{\mathrm{s}}$ is large, then spiral perturbations can drive the unbinding of the binary from the GC. Therefore, the disruption rates do not always increase monotonically with these parameters.

3. The calculated disruption rates of massive binaries (binary masses of $\sim 8\,\mathrm{M}_{\odot}$ and higher) are typically a few times $10^{-7}\,\mathrm{yr^{-1}}$, which is a factor of $\sim 10$ higher compared to relaxation by stars only (cf. Table 3). The largest enhancement is a factor of $\sim 25$, assuming an inclination of $i=9.5^{\circ}$, $N_{\mathrm{s}}=1$ spiral arm events, and a surface density of $\tilde{\Sigma}_{\mathrm{s},0,i}=1000\,\mathrm{M}_{\odot}\,\mathrm{pc^{-2}}$. Our rates are similar to the rates found for relaxation by massive perturbers (Perets et al., 2007), indicating that nuclear spiral arms can act in conjunction with massive perturbers, in particular GMCs, to increase the injection rate of binaries to MBHs by orders of magnitude. The numbers of predicted HVSs and S-stars in our models are consistent with the observed numbers within factors of a few (assuming a binary fraction of 0.6 for HVS progenitors, and 0.8 for S-star progenitors).

4. The disruption rates are moderately dependent on the inclination of the progenitor binary, peaking around $10^{\circ}$ (cf. Fig. 16). The difference in the disruption rate between nearly coplanar and perpendicular orbits is a factor of $\sim 2$. Our model therefore provides a novel potential source for anisotropy in the distribution of HVSs, although further, more detailed study, is still warranted. There are no implications for S-star-like stars, since the reorientation time-scale for the latter is short compared to the typical S-star lifetime.

5. Our TDE rates are up to $20\%$ higher compared to relaxation by single stars only. Although the enhancement by nuclear spiral arms is not as significant as for binaries, it complements other models for enhanced TDE rates.

6. In addition to enhancing TDE, S-star, and HVS rates, nuclear spiral arms could also accelerate the coalescence of binary MBHs. Other implications include the disruption of clusters, building up nuclear star clusters, and supplying the inner regions of nuclear stars cluster with gas, triggering high-mass star formation. Also, the disruption rate of planets and planetesimals could be enhanced.

Acknowledgements

We thank Scott Tremaine, Ben Bar-Or and Nicholas Stone for stimulating discussions and comments on the manuscript, and the referee for a thorough review. ASH and HBP acknowledge support from the Israeli I-CORE center for astrophysics under ISF grant 1829/12; HBP acknowledges support from the Asher space research institute in the Technion. ASH gratefully acknowledges support from the Institute for Advanced Study, and NASA grant NNX14AM24G.

Appendix A Loss-cone diffusion coefficients for relaxation driven by stars

As discussed in Section 4.1, the orbit-averaged angular-momentum diffusion coefficient that describes the diffusion rate in angular-momentum space by stars, $\bar{\mu}_{\mathrm{stars}}(\mathcal{E})$, can be computed from the distribution function $f(\mathcal{E})$ and the potential $\psi(r)$. The formal definition of the diffusion coefficient is

[TABLE]

where $\langle(\Delta\mathcal{R})^{2}\rangle$ is the second-order diffusion coefficient in $\mathcal{R}\equiv(J/J_{\mathrm{c}})^{2}$ describing angular-momentum relaxation by stars. Using standard expressions for $\langle(\Delta\mathcal{R})^{2}\rangle$ (e.g., Cohn & Kulsrud 1978), $\bar{\mu}_{\mathrm{stars}}(\mathcal{E})$ can be expressed in terms of $f(\mathcal{E})$ as

[TABLE]

where $\ln(\Lambda)=\ln[M_{\bullet}/(2M_{\star})]$ is the Coulomb logarithm, and the integrals are given by

[TABLE]

Appendix B Disruption rates and expected number of HVSs and S-stars

Table 4 gives a complete list of the disruption rates (cf. the top part of Table 3). A complete list of the expected number of HVSs and S-stars (cf. the bottom part of Table 3) in given in Table 5.