Forming young and hypervelocity stars in the Galactic Centre via tidal disruption of a molecular cloud

ABSTRACT

The Milky Way Galaxy hosts a four million solar mass black hole, Sgr A*, that underwent a major accretion episode approximately 3–6 Myr ago. During the episode, hundreds of young massive stars formed in a disc orbiting Sgr A* in the central half parsec. The recent discovery of a hypervelocity star (HVS) S5-HVS1, ejected by Sgr A* five Myr ago with a velocity vector consistent with the disc, suggests that this event also produced binary star disruptions. The initial stellar disc has to be rather eccentric for this to occur. Such eccentric discs can form from the tidal disruptions of molecular clouds. Here, we perform simulations of such disruptions, focusing on gas clouds on rather radial initial orbits. As a result, stars formed in our simulations are on very eccentric orbits (⁠|$\bar{e}\sim 0.6$|⁠) with a lopsided configuration. For some clouds, counterrotating stars are formed. As in previous work, we find that such discs undergo a secular gravitational instability that leads to a moderate number of particles obtaining eccentricities of 0.99 or greater, sufficient for stellar binary disruption. We also reproduce the mean eccentricity of the young disc in the Galactic Centre, though not the observed surface density profile. We discuss missing physics and observational biases that may explain this discrepancy. We conclude that observed S-stars, HVSs, and disc stars tightly constrain the initial cloud parameters, indicating a cloud mass between a few × 10⁴ and |$10^5\, {\rm M}_{\odot }$|⁠, and a velocity between ∼40 and 80 km s⁻¹ at 10 pc.

black hole physics, hydrodynamics, Galaxy: centre

1 INTRODUCTION

The central parsec of our Galaxy is home to a population of young stars, consistent with a star formation episode close to the supermassive black hole (SMBH) ∼3−6 Myr ago (Lu et al. 2013). The young stars may be divided into two apparently distinct kinematic structures: the S-stars and the clockwise disc. The S-star cluster (≲ 0.04 pc) has an isotropic and thermal eccentricity distribution (Gillessen et al. 2017) and the majority appear to be B stars. Considering the strong tidal forces in S-star cluster, it is thought that its stars did not form in place, but migrated from larger scales via binary disruption (Hills 1988; Gould & Quillen 2003) or interaction with a gas disc (Levin 2007). The clockwise disc (0.04 ∼ 0.5 pc) is so-called due to the common motion of the stars on the sky (Levin & Beloborodov 2003). The disc contains many young O and Wolf–Rayet stars in addition to B stars. Dynamical modelling yields an average stellar eccentricity of ∼0.3 (Yelda et al. 2014), though this measurement may be biased by binary stars as described by Naoz et al. (2018). Recent observations indicate that the age of the S-stars is consistent with the age of the clockwise disc (Habibi et al. 2017), suggestive of a common origin. Furthermore, Koposov et al. (2020) discovered a hypervelocity star (S5-HVS1) that was ejected from the Galactic Centre (GC) 4.8 Myr ago. S5-HVS1’s velocity vector is also consistent with a disc origin, though there is an |$\sim 10{{\ \rm per\ cent}}$| probability of a random alignment, considering the ∼7° half-opening angle of the disc (Lu et al. 2009).

A parsimonious explanation for these observations is that the formation of the disc produced an episode of binary disruptions, which would leave remnant stars in the S-star cluster. The disrupted binaries may come directly from the disc, or from an older background population. The spectroscopic age of S5-HVS1 (25-93 Myr; Koposov et al. 2020), suggests that it was a background star excited to disruption by the disc.

Madigan, Levin & Hopman (2009) and Generozov & Madigan (2020) showed that if the stars formed in a lopsided, eccentric disc, its binaries could have been excited to nearly radial orbits and tidally disrupted via the Hills (1988) mechanism. The excitation to radial orbits occurs as stellar orbits precess retrograde to their orbital angular momenta due to the influence of the surrounding star cluster. Orbits with higher eccentricities precess more slowly and end up behind the bulk of the disc, which then torques them to even higher eccentricities. The remnant bound stars’ orbits are initially highly eccentric, but approach a thermal distribution via resonant relaxation (Perets et al. 2009; Antonini & Merritt 2013; Generozov & Madigan 2020; Tep et al. 2021).

Hydrodynamical simulations show that formation of a lopsided, eccentric disc can arise from the infall of a turbulent |$10^4{-}10^5\, {\rm M}_{\odot }$| molecular cloud (Bonnell & Rice 2008; Gualandris, Mapelli & Perets 2012) or collisions of such clouds in the vicinity of the GC (Hobbs & Nayakshin 2009, HN09 hereafter). Infall of smaller molecular clumps with |$\sim 100\, {\rm M}_{\odot }$| can explain smaller compact groups of young stars in the GC like IRS 13N (Jalali et al. 2014).

Gualandris et al. (2012) initialized their N-body simulations of stellar dynamics with stellar orbits from the hydrodynamical star formation simulations of Mapelli et al. (2012). They found that eccentric disc instability indeed occurred and led to the increase of stellar eccentricities of some of the stars. However, the effect was insufficient to produce a significant number of stellar binary disruptions. Gualandris et al. (2012) noted that this result could have been due to a relatively low (e ∼ 0.2−0.4) eccentricities of stars formed in the hydrodynamical simulations in Mapelli et al. (2012), and that perhaps a molecular cloud on a more radial orbit could lead to a better quantitative agreement with observations. In this paper, we attempt to find cloud initial conditions that can produce more eccentric discs. We discuss the differences between our simulation setup and that in Mapelli et al. (2012) and Gualandris et al. (2012) in Section 5.5. We perform hydrodynamical simulations in which a single molecular cloud infalls on to the GC, shocks, and settles into an eccentric gaseous disc. The disc then fragments into stars. We use the orbits of the stars to initialize high-precision N-body simulations to track close encounters with the central SMBH and to make observational comparisons.

The remainder of this paper is organized as follows. In Section 2, we describe our overall approach to simulating gas cloud disruptions, using independent hydrodynamic and N-body simulations. In Section 3, we describe the setup and results for our hydrodynamic simulations. In Section 4, we summarize the setup and results of our N-body simulations, and present detailed observational comparisons. In Section 5, we discuss limitations of our simulations, and additional predictions. Finally, we summarize our findings in Section 6.

2 JOINT GAS DYNAMICS AND STELLAR N-BODY SIMULATIONS

The goal of this paper is to try and reproduce what we know about young stars in the GC from first-principle simulations. Ideally, we would like to model gas dynamics, star formation, and long-term evolution of stellar orbits with one code. However, this is currently impractical due to numerical limitations. Fortunately, the problem can be broken into two largely independent sub-problems: star formation and the subsequent orbital evolution. For both of these well-tested numerical approaches are available.

In the first part of our numerical investigation, we perform smoothed particle hydrodynamics (SPH) simulations of gas cloud infall on to the central parsec and the resulting star formation. These simulations are tailored to the relatively short gas-rich epoch of ∼(0.1−0.3) Myr in duration. Due to numerical expense in resolving close particle interactions, including via shocks, the SPH simulations cannot reach sufficiently close to Sgr A* to simulate highly eccentric orbits that would produce S-stars and hypervelocity stars (HVS) via the Hills (1988) binary disruption mechanism. However, these simulations are sufficiently robust to resolve star formation on scales R ≳ 0.01 pc. We recall that star formation inward of this region is not in fact expected (Nayakshin & Cuadra 2005) or observed (Paumard et al. 2006) because gas radiative cooling times are too long inside radius of ∼1 arcsec (0.04 pc) to permit disc fragmentation. Our gas dynamical simulations are therefore designed to capture the formation of the well-known ‘disc stars’ (e.g. Levin & Beloborodov 2003) from first principles, that is together with the gas disc formation.

When at least half of the initial gas mass is turned into stars or accreted on to Sgr A* we stop these SPH simulations and use the data for the stars formed in the simulations to initialize the N-body (stars-only) simulations. The latter follow a much smaller number of particles (stars), N ≲ 10³ but are able to address orbits that happen to pass much closer to Sgr A*. These simulations are also run for much longer duration, ∼5 Myr. Our N-body simulations are therefore streamlined to simulate – potential – transformation of the disc stars into S-stars via the development of a secular instability.

3 HYDRODYNAMICAL SIMULATIONS

3.1 Background potential

Sgr A* is surrounded by a massive nuclear star cluster (Launhardt, Zylka & Mezger 2002). When a massive gas cloud strikes the central parsecs of the Galaxy and gets bound to it, both Sgr A* and the cluster must respond dynamically. Readjustment of the stellar cluster and the SMBH-cluster interactions during that is a complex problem that is outside of the scope of our investigation due to numerical limitations. Observations indicate that presently Sgr A* is exactly in the centre of the Galaxy with a very small 3D velocity (Reid & Brunthaler 2020). Therefore, for simplicity we keep Sgr A* – an SMBH with the initial mass of M_BH = 4 × 10⁶M_⊙ – fixed in the centre of the Galaxy for the duration of our SPH simulations. The black hole is surrounded by a spherically symmetric static stellar cusp with a broken power-law profile, motivated by Schödel et al. (2018)’s fits to the surface brightness profile in the GC. In particular, the stellar density is

$$\begin{eqnarray} \rho (r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho _o \left(\frac{r}{3}\right)^{-1.16} & r \le 3 \,{\rm pc}\\ \rho _o \left(\frac{r}{3}\right)^{-2.9} & r \gt 3 \,{\rm pc}, \end{array}\right. \end{eqnarray}$$

(1)

where r is the distance from the GC. The density is normalized so that the mass enclosed within 3 pc is |$9.1\times 10^6\, {\rm M}_{\odot }$| (the mass within 1 pc is |$\sim 1.2\times 10^6\, {\rm M}_{\odot }$|⁠). We model the effects of the stellar cusp by adding an extra acceleration to each gas or star particle, viz.

$$\begin{eqnarray} a= -\frac{G M_{\rm *} (r)}{r^2}, \end{eqnarray}$$

(2)

where M_*(r) is the stellar mass within the star’s distance from the origin, r. We use this potential for both SPH and N-body simulations.

3.2 Simulation physics

Our star formation simulations are performed with the well-known cosmological hybrid N-body/SPH code Gadget-2 (Springel 2005) that includes gas hydrodynamics, cooling, gravity, N-body dynamics, and conversion of dense self-gravitating regions of gas into stars. As in HN09, we use sink particle prescription to model both Sgr A* and newly born stars. In particular, the SMBH interacts with gas and stars via gravity but also accretes SPH and star particles if they approach it closer than the capture radius r_c = 5 × 10⁻³ pc and if they are gravitationally bound to it. Hence, the SMBH mass increases during the simulations. Stars are treated in a similar vein except their accretion radius is much smaller, r_c = 10⁻⁴ pc, and they are of course not fixed in space, free to follow their orbits. When stars accrete SPH particles, they accrete not only mass but also momentum of the latter.

In HN09, gas radiative cooling was modelled with an idealised ‘β-cooling’ approach standard to simulations of self-gravitating discs (Rice, Lodato & Armitage 2005; Nayakshin, Cuadra & Springel 2007). This approach fixes the gas cooling time to be proportional to the local orbital period, which is useful in differentiating between fragmenting and non-fragmenting self-gravitating discs Gammie (2001) but is less suitable to spherical geometries. Here, we use more self-consistent approach following Lombardi, McInally & Faber (2015) in which radiative cooling rate is a function of the local gas density and temperature. We use the gas/dust opacity from Zhu, Hartmann & Gammie (2009).

When gas cools and contracts due to its own gravity sufficiently, dense self-gravitating gas condensations form. These have densities up to multiple orders of magnitude higher than the surrounding gas, and they continue to collapse, numerically speaking to a point. Star particles are introduced when gas density exceeds ρ_sf = 10⁻⁸ g cm⁻³. We note that ρ_sf is much higher than the local tidal density for relevant distances from SMBH, ensuring that star particles are introduced only when local gravitational collapse is well underway. Unlike HN09, we do not allow mergers between star particles to avoid potentially significant modification of their orbital elements if two-star particles on very different orbits merge.

Our unit mass is M_u = 1 M_⊙, unit length L_u = 1 pc, and unit time is t_u = 1.49 × 10⁷ yr.

3.3 Initial conditions

We do not have direct observational constraints on what the initial gas configuration looked like. When choosing initial conditions here, we are motivated by (1) the need to bind a significant quantity of gas to the SMBH in the centre of the Galaxy (this we would not have the young stars we observe in the Galaxy’s centre Paumard et al. 2006); (2) the hint provided by the eccentric disc instability mechanism for producing HVS and S-stars. For the instability to work, the initial stellar orbits need to be rather eccentric (Madigan et al. 2009; Gualandris et al. 2012; Generozov & Madigan 2020). This obviously requires a non-axisymmetric deposition of gas into the central parsec, and strong shocks in which a large fraction of gas kinetic energy can be lost to radiation. For a single cloud falling on to Sgr A*, this favours a cloud on a rather radial orbit, so that different sides of the cloud pass the SMBH on opposite sides, collide behind it, and then get bound to the central parsec.

Beyond these general principles,there are only very rough constraints. The cloud probably originated within the inner ∼200 pc since the OB disc stars appear coeval within ∼1−2 Myr (Paumard et al. 2006; Yelda et al. 2014); the free-fall time of gas from distances larger than that would likely be spread over a longer time interval and hence result in a star formation event distributed over a longer time period. Exploring such initially distant gas clouds is however prohibitively expensive for us here because we wish to resolve regions as close as ∼0.01 pc to the SMBH.

We started our investigation with a gas cloud positioned 30 pc away from the GC. This resulted in stellar discs much larger than the ∼0.5 pc observed disc of young stars (Paumard et al. 2006). More reasonable results were found for a gas cloud initially a distance R₀ = 10 pc away from the GC; this is the scenario we pursue here. For simplicity, we assume a uniform density spherical cloud with radius R_c = 8 pc and mass M_c = 10⁴, 3 × 10⁴, or 10⁵M_⊙ with a velocity vector that would take the centre of the cloud to pericentre distance R_p = 0.21 pc, if the gas self-interaction (shocks and gravity) were neglected. We found that none of the lowest mass gas clouds, |$M_{\rm c} = 10^4\, {\rm M}_{\odot }$|⁠, resulted in a sufficient mass of young stars to explain either the observed clockwise disc stars, or the S-stars, so we do not present these experiments here. Furthermore, we explored a large range of the initial magnitude of cloud velocity, v₀, while keeping the pericentre passage distance the same for all of the simulations. This corresponds to a positive fixed y-velocity component for all of our initial gas cloud configurations (see Fig. 1), and an anticlockwise motion of the cloud centre.

Figure 1.

Snapshots of gadget Cloud 1 simulation at selected times. Gas density is shown with the colour map, with the colour bar in code units (M_⊙ pc⁻³ ≈ 7 × 10⁻²³ g cm⁻³). Cyan pluses mark positions of the stars formed in the simulation. The position of Sgr A* is indicated with the ‘+’ symbol at (x, y) = (0, 0). Note that gas orbits evolve from clockwise oriented at early times to anticlockwise at late times. See the text for discussion of this effect.

We found that all of the simulations that started with a relatively large initial velocity, v₀ ≳ 100 km s^-1, produced young stars on orbits distinct from those seen in the GC. Such velocities are still somewhat lower than the local escape velocity (∼170 km s⁻¹), so the gas is bound to the GC. However, we find that the first passage tends to bind too small a fraction of the cloud to the SMBH, and so a second approach is needed to form a strongly bound disc. During this evolution, star formation begins in shocked filaments, rather than in the disc, and results in stars on eccentric orbits with semimajor axis typically greater than 1 pc. This is in contrast to the observed young stars that are contained to the central ∼0.5 pc (Paumard et al. 2006). We also found that clouds on orbits with smallest initial velocities, v₀ ≲ 20 km s⁻¹, resulted in too strong cancellation of the initial angular momentum via shocks. This led to discs too small in radial extent to explain the clockwise disc. Below, we therefore focus on the middle range in v₀ that led to star formation most reminiscent of the observed young stars in the GC.

3.4 Overview of SPH simulation results

The parameters of the SPH simulations that we chose to follow up with N-body simulations are shown in Table 1. Here, while describing the resulting gas dynamics and star formation, we centre our discussion around two numerical experiments specifically, which we refer to as Cloud 1 and Cloud 2. These two appear to reproduce observations best and are listed in the top two rows of Table 1.

Table 1.

Cloud simulations that form well-defined stellar discs, and are chosen for follow-up N-body simulations. The first and second rows correspond to Clouds 1 and 2, respectively, the most promising initial conditions for reproducing observations of the GC.

Cloud mass	Cloud velocity	Snapshot timea	Disc mass	Accreted massb	<e>	f(j/jlc < 10)c	Percent disruptedd	Percent Disruptede
(M⊙)	(km s−1)	(yr)	(M⊙)	(M⊙)				(high resolution)
(C1) 1 × 105	82	2.2 × 105	2.2 × 104	2.65 × 104	0.55	0.5	7–8 [7.8]	7.8
(C2) 3 × 104	41	3.0 × 105	8.1 × 103	4.71 × 103	0.55	0.55	4–8 [5.6]	5.5
1 × 105	54	3.0 × 105	5.5 × 103	4.11 × 104	0.63	0.045	0–0 [0.0]
1 × 105	54	3.7 × 105	1.1 × 104	4.13 × 104	0.56	0.05	0-1 [0.2]
1 × 105	54	4.5 × 105	1.3 × 104	4.15 × 104	0.53	0.094	0–0 [0.0]
1 × 105	54	5.4 × 105	1.3 × 104	4.16 × 104	0.48	0.08	0–0 [0.0]
3 × 104	54	4.5 × 105	2.9 × 103	1.15 × 104	0.62	0.036	0–0 [0.0]
3 × 104	54	6.0 × 105	3.4 × 103	1.17 × 104	0.58	0.096	0–0 [0.0]
3 × 104	68	4.2 × 105	5.1 × 103	1.03 × 104	0.58	0.27	1–2 [1.2]
3 × 104	82	3.0 × 105	6.6 × 103	8.84 × 103	0.57	0.47	2–7 [4.2]
3 × 104	82	3.7 × 105	6.7 × 103	8.94 × 103	0.54	0.48	1–3 [2.4]

^aStarting time for N-body simulations.

^bMass accreted by Sgr A* during SPH simulations.

^cFraction of particles that start with angular momentum within a factor of 10 of the loss cone.

^dPercent of disc particles that become disruptees over 5 Myr from five different 100-particle simulations. The number in bracket is the mean.

^ePercent of disc particles that become disruptees in high-resolution simulations (see the text for details).

Fig. 1 shows several snapshots exemplifying the complex dynamics of gas and its gradual conversion into stars in the Cloud 1 simulation. The last snapshot shown in the figure is used to initialize N-body simulations in the next section. We see that early on (in the t = 0.06 Myr snapshot in Fig. 1) gas with clockwise rotation dominates disc formation in the central ∼0.3 pc. The disc is very eccentric (with typical eccentricity e ∼ 0.6) and is not smooth because it forms from gas falling towards SMBH along a shocked filament located outside the figure, at y ≈ 0 along large negative x. As gas infall into the centre continues, some of the filaments making up the disc become dense enough to be self-gravitating. This results in the formation of self-bound, high-density clumps inside the filaments, some of which can be seen as nearly white coloured clumps in the t = 0.1 Myr snapshot. When the star formation density threshold is exceeded in these clumps, star particles are introduced. These are shown with the cyan crosses (the very central cross is Sgr A* itself).

Since stars are born from the gas, they inherit the orbital parameters of the parent SPH particles, and initially follow very similar orbits. However, new gas arriving in the central region constantly mixes with the gas already there. Thus, gaseous orbits evolve under the effects of both gravity and shocks and/or gas pressure forces whereas star evolve mainly due to gravity.¹ This leads to a gradual divergence in orbits of stars and gas.

One can also see in the top row of snapshots shown in Fig. 1 some apsidal orbital precession for both gas and stars due to the stellar potential in the GC. However, the bottom row of the snapshots shows strikingly different orbital structure for stars and gas. There is a depletion of gas in the central region, with most of the mass there being in the young stars (not counting the background potential, of course). Furthermore, as time progresses, the sign of the angular momentum of the dominant inflow of gas into the central region flips. By the time of the last snapshot shown the gas rotates predominantly counterclockwise in the snapshot. This very rapid evolution is due to star formation, strong shocks, and angular momentum cancellation, an effect that we shall discuss in greater detail below.

Fig. 2 shows gas and star dynamics in the central region for Cloud 2 simulation. The orientation of angular momentum vector for the gas dominating the infall at all times is counterclockwise, so we do not observe a flip in the angular momentum vector as in the Cloud 1 simulation. The forming disc is still eccentric but somewhat less so, and is also more uniform. This indicates that shocks are less important in the central ∼0.2 pc for the Cloud 2 simulation than they were in Cloud 1 simulation. As a result, gravity dominates the dynamics of both gas and stars in the central region, and thus gaseous and stellar orbits in the Cloud 2 simulation are much closer to one another than in Cloud 1 simulation. The last snapshot shown in Fig. 2 is used to start our N-body simulations.

Figure 2.

Same as Fig. 1 but for Cloud 2. Note that the orientation of the dominant component of gas in the very centre is opposite to that of Cloud 1 at all times. Therefore, the stars that form there orbit Sgr A* in the counterclockwise direction.

Fig. 3 compares gas dynamics of Cloud 1 (top row of panels) and Cloud 2 (bottom row of panels). In the former case, gas entering the inner regions (∼0.25 pc) has a clockwise angular momentum at early times (the leftmost panel). By the time t = 0.1 Myr (middle panel), there are two gas streams entering the central region, one with clockwise angular momentum and one with anticlockwise angular momentum. Finally, at later times (the right-hand panel), there is only one dominant orientation of angular momentum for the incoming gas: it is now anticlockwise. Note that the anticlockwise direction is the dominant one in the cloud, so it ‘wins’ eventually in Cloud 1 simulation. In the Cloud 2 simulation, however, we observe that the clockwise filament feeding the central region is never dominant. Instead, it is the anticlockwise one that always dominates gas supply to the central disc.

Figure 3.

Cloud 1 (top panel) and Cloud 2 (bottom panel) at selected times (0.052, 0.1, and 0.13 Myr from the left- to the right-hand panel, respectively) on larger physical scales than shown in 1 and 2. Note that in Cloud 1 simulation the gas inflowing into the very central region is dominated by the filament oriented along the negative x-axis and has a clockwise angular momentum.

Fig. 4 presents angular momentum distribution for SPH particles in the simulations for Cloud 1 and Cloud 2 at several times. The dashed vertical line marks purely radial orbits (j_z = 0). Note that the mean of j_z at t = 0 is exactly the same for two clouds, but because the 3D initial velocity of Cloud 1 is larger, the initial distribution of angular momenta for Cloud 1 is broader than that for Cloud 2. The legend lists the fraction of SPH particles that have ‘clockwise’ orbits (with negative j_z), F_cw.

Figure 4.

Specific angular momentum distribution of gas in Cloud 1 and Cloud 2 at a few different snapshots. For each cloud and snapshot, we indicate ‘F_cw’–the fraction of particles on clockwise orbits with negative angular momentum.

The higher F_cw for Cloud 1 is very important. At t = 0.06 Myr, we observe low angular momentum peaks for both Clouds 1 and 2, but they have opposing signs. These peaks are formed due to shocks. As streams with positive and negative angular momentum collide, the absolute value of angular momentum, |j_z| is reduced for both streams. The exact location of the peak in j_z is, however, not always trivially predicted since it depends on which of the streams gets into Sgr A* vicinity first, their relative densities and the amount of mass and momentum present in both. We see from Fig. 4 that in Cloud 1 the negative j_z material has the upper hand but only at early times. Eventually, as more and more gas with positive angular momentum falls in, we see that the negative shoulder of j_z distribution disappears in both Clouds 1 and 2.

However, this – temporary – victory of negative j_z gas in the inner disc has a lasting effect on the resulting stellar distribution. Fig. 5 shows the distribution of angular momenta for stars formed in these two simulations at several different times. We can see that the stars in Cloud 1 are mainly rotating clockwise, whereas those in Cloud 2 all rotate anticlockwise. Note that if we continued the SPH simulations for longer, more stars in Cloud 1 with anticlockwise orbits could have formed, but they are likely to be outside of the central ∼0.5 pc that we focus on here given their relatively large angular momentum.

Figure 5.

Specific angular momentum distribution of star particles formed from Cloud 1 and 2 (thin lines). For reference, we also plot the distribution of gas at the same times with thick lines.

3.5 Accretion on to Sgr A*

As shown in the fifth column of Table 1, Sgr A* accretes between ∼5 × 10³ and |$4\times 10^4\, {\rm M}_{\odot }$| of stars and gas (16–|$42{{\ \rm per\ cent}}$| of the mass of the original cloud) over the course of our simulations.

This accretion can explain the Fermi Bubbles–kpc-scale gamma-ray structures, extending above and below the GC (Su, Slatyer & Finkbeiner 2010). The bubbles indicate Sgr A* was active within the last few Myr, and output 10⁵⁵ − 10⁵⁷ erg of energy during this period (Guo & Mathews 2012). All of the clouds in Table 1 can account for this energy, as long as the radiative efficiency exceeds |$\sim 0.1{{\ \rm per\ cent}}$|⁠.

4 N-BODY SIMULATIONS

It is necessary to switch to a high-accuracy N-body integrator to dynamically evolve the stellar disc going forward. Binaries that undergo disruption by the SMBH move along high eccentricity orbits which require short, adaptive time-steps to minimize energy and angular momentum error particularly near pericenter. We must also do away with the capture radius of the SMBH in the hydrodynamical code to retain high eccentricity orbits at pericenter.

We select snapshots from the hydrodynamic simulations that formed well-defined discs for follow-up N-body simulations. In particular, we select snapshots where roughly half of the initial gas mass has either been accreted by Sgr A* or turned into stars. The full set of cloud masses, velocities and snapshot times we used for follow-up N-body simulations are listed in Table 1.

The positions and velocities of the stars and central SMBH are used to initialize REBOUNDN-body simulations (Rein & Liu 2012). We exclude unbound stars and those with semimajor axes greater than 0.5 pc. (We do not expect the outer regions of the disc to be dynamically relevant, considering the short range of secular torques; Gürkan & Hopman 2007). For computational efficiency, we remove any binaries from the initial conditions by deleting one star from each binary pair. We ignore the remaining gas from the hydrodynamical simulations in our N-body simulations. This is in contrast to Gualandris et al. (2012), who include in their simulations an axisymmetric potential modelled on the remnant gas mass profile. The dynamics of the stars is dominated by the more massive spherically symmetric stellar cusp and we note that Gualandris et al. (2012) found similar dynamical evolution in simulations that excluded the gas disc potential. Šubr & Haas (2016) also ran N-body simulations of an eccentric disc in the GC, motivated by Mapelli et al. (2012), but did not include the stellar cusp which is crucial for comparing to observational data.

To quickly identify the most promising Clouds for producing binary disruptions and the S-stars, we use low-resolution simulations. Instead of using all of the stars, we initialize five REBOUND simulations with different random sets of 100 stars from each snapshot.² As discussed in Section 4.1, we then select the most promising Clouds (the aforementioned Clouds 1 and 2) for follow-up, higher resolution simulations. We find excellent agreement for the number of disruptions produced by high- and low-resolution simulations, though we find the eccentricity distribution and surface density profile of the disc are resolution dependent, as discussed in Section 4.2.

The mass of the central SMBH is |$\sim 4\times 10^6\, {\rm M}_{\odot }$|⁠, and is taken directly from the hydrodynamic simulations. The disc mass in the N-body simulation is the total mass of all bound stars with semimajor axes less than 0.5 pc in the corresponding Gadget snapshot (⁠|$2.2\times 10^4$| and |$8.1\times 10^3\, {\rm M}_{\odot }$| for Clouds 1 and 2, respectively). Stars are initialized with the same mass in the N-body simulations. Each disc is evolved forward in time for 6 Myr using the IAS15 integrator (Rein & Spiegel 2015). Perturbations from the surrounding stellar cusp are included with the REBOUNDX package (Tamayo et al. 2020).

Our N-body simulations are initialized with the centre of mass of the black hole–disc system at the origin. However, in the presence of the external cusp force the centre of mass can drift. To be consistent with our hydrodynamic simulations (where the central black hole is fixed), we keep the centre of mass of the black hole and disc fixed at the origin.

We also performed a few simulations, where the black hole and disc are allowed to drift. In this case, the SMBH feels the acceleration of the stellar cusp if it moves more than 5 × 10⁻⁴ pc from the origin, where the enclosed stellar mass is |$\approx 1\, {\rm M}_{\odot }$|⁠. The results are consistent with simulations where the centre-of-mass drift is negated.

Unlike, the simulations of Generozov & Madigan (2020), the simulations here do not include post-Newtonian corrections. This is because, simulations with post-Newtonian effects can have significant numerical artifacts. For example, the SMBH can receive spurious velocity kicks and start drifting from the centre. Such kicks occur when particles penetrate within tens of gravitational radii of Sgr A*, where the first-order post-Newtonian prescription in REBOUNDX is inadequate. In most post-Newtonian simulations, such kicks do not occur, and we find good agreement with Newtonian simulations (see also Wernke & Madigan 2019).

4.1 Results of N-body simulations: disruptions

Table 1 shows statistics on ‘binary disruptions’ in our simulations. The initial conditions do not have any binaries, as they are deleted for computational efficiency before the start of the simulation. However, we record a disruption whenever particles cross within 3 × 10⁻⁴ pc of the central SMBH.³ We do not delete such particles (‘disruptees’) from our simulations, though only a single disruption is recorded for each one. We also do not include disruptees in plots of orbital elements at late times. In reality, binaries would typically be split into two stars that are removed from the disc (with one deposited at small semimajor axes and one ejected from the system). To test the effect of this removal, we ran additional simulations, where disruptees were merged with the central black hole. We found this did not affect the disruption statistics.

We have identified two promising snapshots for producing binary disruptions, where |$6{-}8{{\ \rm per\ cent}}$| of the particles became disruptees, corresponding to the first and second row of Table 1 (Cloud 1 and Cloud 2). These two clouds produce enough disruptions to explain the bulk of the S-star population, as discussed below.

Following Generozov & Madigan (2020), we expect

$$\begin{eqnarray} N_s=f_{\rm mf} N_{\rm d} f_{\rm d} (1-f_{\rm c}). \end{eqnarray}$$

(3)

S-stars within the observed mass range, where N_d is the total number of binaries in the disc, f_d is the fraction that are disrupted, f_c is the fraction of binaries that collide instead of disrupting, and f_mf is the fraction of S-stars within the observed mass range (⁠|$\sim 3{-}15\, {\rm M}_{\odot }$|⁠). Assuming f_d ≈ 0.06−0.08, as it is for Clouds 1 and 2, and assuming N_d = 1000, f_c = 0.2, and f_mf = 0.33 as in Generozov & Madigan (2020), we would predict 21 (16) S-stars between 3 and |$15\, {\rm M}_{\odot }$| for Cloud 1 (Cloud 2). The collision probability, f_c, may be a factor of ∼2−3 higher once repeated close encounters with the SMBH are taken into account (Generozov & Perets 2021). On the other hand, N_d could also be higher by a factor of a few. Observationally, there are ∼20 early-type S stars with semimajor axes less than 0.03 pc (Gillessen et al. 2017).

The other simulations in Table 1 produced compact eccentric discs, but few disruptions. In some cases, this can be explained by the small mass of the stellar disc, as some snapshots have stellar discs with only a few|$\times 10^3\, {\rm M}_{\odot }$|⁠. We also find a correlation between the fraction of stars at small angular momentum and the number of disruptions, as shown in Fig. 6. The top panel of this figure shows the number of disruptions versus the harmonic mean of the angular momentum (in units of the loss cone angular momentum) for different clouds. (We use the clouds in Table 1, but fix the stellar mass to |$10^4\, {\rm M}_{\odot }$| in each case to better isolate the effect of the angular momentum distribution). As the fraction of stars near the loss cone increases, the harmonic mean and the number of disruptions also increase. In the bottom panel of Fig. 6, we show the number of disruptions as a function of the fraction of stars within a factor of 10 of the loss cone. There is a positive correlation between the number of disruptions and this fraction too.

Figure 6.

Mean number of disruptions for different clouds, as a function of the harmonic mean of the stars’ angular momentum (top panel) and the fraction of stars within a factor of 10 of the loss cone (bottom panel). The clouds are the same as those in Table 1, except the masses are fixed to |$10^4\, {\rm M}_{\odot }$| to better isolate the effect of the angular momentum distribution. Clouds with more stars at low angular momenta have more disruptions. These data are from low-resolution (100-particle) simulations. Error bars correspond to the standard error of the mean.

So far, we have shown results from stacked 100-particle simulations. To check whether small particle numbers are affecting our results, we ran higher resolution simulations of Clouds 1 and 2 (For Cloud 1 the higher resolution simulation had 900 stars after the binary deletion step, which is approximately half the number of bound stars within 0.5 pc in the corresponding hydrodynamic snapshot. For Cloud 2, there were 328 stars, close to the total number within 0.5 pc). As shown in the final column of Table 1, the fraction of disruptees is independent of resolution. Thus, we expect our prediction for the fraction of binaries disrupted to be robust.

However, as discussed in the following section, the semimajor axis and eccentricity distribution is resolution dependent, with the disc becoming more spread out and eccentric in lower resolution simulations, due to artificially strong two-body relaxation.

4.2 Eccentricity and semimajor axis distribution

Fig. 7 shows the initial semimajor axis and eccentricity distribution for the high-resolution Cloud 1 and 2 simulations. (Note that we use t_N to denote time from the start of our N-body simulations). Initially, stars formed in each cloud have a similar spread of semimajor axes and eccentricities.⁴ The outer stellar orbits in Cloud 1 are initially more eccentric than the inner orbits. This positive eccentricity gradient is caused by more efficient gas-driven circularization in the inner disc.

Figure 7.

Initial profile of stellar orbital eccentricity versus semimajor axis from high-resolution Cloud 1 (left-hand panel) and Cloud 2 (right-hand panel) simulations. Colours show inclinations.

Fig. 8 shows the eccentricity and semimajor axis distribution (excluding disruptees) of the discs formed by Clouds 1 and 2 after 5 Myr of evolution in N-body simulations. After 5 Myr, there is no longer any correlation between eccentricity and semimajor axis for Cloud 1, though a negative correlation is present for the Cloud 2 disc (i.e. orbits at smaller semimajor axes are, on average, more eccentric). Also, the discs have lost their apsidal alignment due to differential precession (see Fig. 9).

Figure 8.

Eccentricity versus semimajor axis after 5 Myr for Cloud 1 (top row) and Cloud 2 simulations (bottom row). The left-hand and centre panels show the distribution from low- and high-resolution simulations, respectively. In the former case, we use stacked data from five different simulations. The right-hand panels compare the eccentricity distributions from these simulations. Black points are individual particles from the simulations. Red contours show the particle density in this space, and are constructed by binning the black points and smoothing with a Gaussian kernel.

Figure 9.

Initial (top panel) and final (bottom panel) orbits for the high-resolution Cloud 1 (left-hand panel) and Cloud 2 (right-hand panel) simulations.

As shown in Table 2 and Fig. 8, the semimajor axis and eccentricity distribution of the disc is dependent on resolution. For example, the final mean eccentricity of the Cloud 1 stars is 0.48 in the low-resolution simulation and 0.38 in the high-resolution simulations. This is expected, as the disc heats up via two-body relaxation, increasing its mean eccentricity (Alexander, Begelman & Armitage 2007).

Table 2.

Mean eccentricity after 5 Myr in low- and high-resolution Cloud 1 and 2 simulations.

Simulation	<e>a	<e> (Low inc)b	|$t_{\rm rx}^{c}$|	|$t_{\rm rx, real}^{d}$|
			(yr)	(yr)
Cloud 1 (low resolution)	0.48	0.46	3.0 × 104	1.3 × 105
Cloud 1 (high resolution)	0.38	0.33	1.8 × 105	1.3 × 105
Cloud 2 (low resolution)	0.43	0.36	9.5 × 104	1.8 × 105
Cloud 2 (high resolution)	0.40	0.31	2.5 × 105	1.8 × 105

^aMean eccentricity of all bound particles in simulation

^bMean eccentricity of all low inclination orbits – those within the median inclination. Inclination is defined with respect to the mean orientation of the disc.

^cInitial relaxation time of a disc annulus at 0.1 pc in the simulation, estimated from equation (3) in Alexander et al. (2007).

^dEstimate for initial disc relaxation time with a realistic mass function (m^−1.7, extending from 1 to |$100\, {\rm M}_{\odot }$|⁠).

The relaxation time is proportional to (N|$\langle$|m²|$\rangle$|⁠)⁻¹, where N is the number of stars and |$\langle$|m²|$\rangle$| is the second moment of the mass function. We estimate the initial relaxation time at 0.1 pc for the discs formed from Cloud 1 and Cloud 2, in Table 2, using equation (3) from Alexander et al. (2007).⁵ For comparison, we also show the relaxation time for a realistic m^−1.7 mass function, extending from 1 to |$100\, {\rm M}_{\odot }$| (for the disc masses in Table 1).⁶ The relaxation times for the high-resolution simulations are within |$40{{\ \rm per\ cent}}$| of the realistic estimates. Therefore, we use the high-resolution simulations for observational comparisons. Also, we focus on low inclination stars for these comparisons, since high inclination stars would not be identified as members of the disc. Here, inclination is measured with respect to the mean of all the unit angular momentum vectors. Low inclination orbits are those with inclinations less than the median (∼13° in high-resolution simulations of both clouds after 5 Myr). Such orbits and the observed clockwise disc have a similar spread in inclination (e.g. Lu et al. 2009 infer the half-opening angle of the disc to be 7° ± 2°.)

The disc eccentricity inferred from observations is 0.27 ± 0.07 (Yelda et al. 2014), which is consistent with the mean eccentricity of low inclination stars in our high-resolution simulations (see the third column of Table 2).

The eccentricity distribution of Clouds 1 and 2 discs after 5 Myr is unimodal like the eccentricity distribution inferred by Yelda et al. (2014) from observations (see Fig. 10). In contrast, in Madigan et al. (2009) and Generozov & Madigan (2020), the eccentric disc instability leads to a prominent bimodal structure. In fact, the shape of the final eccentricity distribution depends on the background density profile. The final eccentricity distribution for the disc simulation with the flatter r^−1.1 background in Generozov & Madigan (2020) also lacks a bimodal structure. In this paper, we also assumed a relatively flat density profile (with a power-law index of −1.16); hence, the bimodal eccentricity distribution is absent at late times.

Figure 10.

Initial and final eccentricity distribution in (high resolution) N-body disc simulations initialized from Cloud 1 (top panel) and Cloud 2 (bottom panel).

The bimodal structure is also time-dependent. At ∼1 Myr, there is a second peak at high eccentricities. However, the eccentricity distribution only appears bimodal when high inclination stars are included. When these are filtered out, the peak at high eccentricities disappears. We find this is the case for various simulations, including those from Generozov & Madigan (2020).⁷ This means we would not expect a bimodal eccentricity distribution to be observed within the disc, as high inclination stars would not be identified as disc members.

Fig. 11 shows a comparison of the surface density profiles from Clouds 1 and 2 with the observed profile from Lu et al. (2009). We scale the simulated profile by a constant to approximately match the overall normalization of the observed profile, considering the observations are incomplete at low stellar masses. The discs in our simulations are significantly more compact than the observed young disc(s) in the GC.

Figure 11.

Surface density of low inclination stars after 5 Myr for Cloud 1 (top panel) and Cloud 2 (bottom panel), compared to the observed surface density profile from Lu et al. (2009). The surface density from simulations is scaled by an arbitrary constant.

Overall, we are able to reproduce the mean eccentricity of the disc, but not its surface density profile. Also, in our simulation, the inclination distribution of the disc becomes broader at small radii, whereas in the inclination distribution of observed clockwise disc becomes narrower at small radii (Naoz et al. 2018). However, our simulations lack several important physical effects as described in the following section and caution is warranted in interpreting such comparisons.

5 DISCUSSION

5.1 Additional effects

In this section, we discuss addition physical effects that could change the physical properties of the simulated discs, as well as observational effects that might bias the observed properties.

• Gas may persist for longer periods of time than we assume here, and circularize the stars via dynamical friction or gravitational torques. Alternatively, massive gas clumps could accelerate the two-body relaxation within the disc, which would increase the average eccentricity, and cause it to spread radially.

• Star formation formation proceeds from the inside out in our hydrodynamic simulations (i.e. star formation first occurs close to the Sgr A*). Thus, some star formation may occur towards the outskirts of the stellar disc after the start of ourN-body simulations. This could make the stellar discs less compact, and more consistent with observations.

• Our N-body simulations lack a mass spectrum. We would expect higher mass stars (that dominate the observed population) to have a lower average eccentricity and inclination than low-mass stars due to two-body relaxation (Alexander et al. 2007). Additionally, secular effects could lead to different mean eccentricities for different species (Foote, Generozov & Madigan 2020; Gruzinov, Levin & Zhu 2020), and lead to heavier species sinking towards the disc mid-plane (Szölgyén & Kocsis 2018). Finally, the heavier species would likely develop a steeper surface density profile than the light stars, potentially exacerbating the tension with the observed surface density profile.

• Unresolved binaries in the disc may bias the observed eccentricity distribution, as discussed in Naoz et al. (2018).

  We have used the procedure described in Naoz et al. (2018) to measure the apparent eccentricity distribution of Cloud 1 and Cloud 2 after 5 Myr, with a hypothetical, isotropic population of binaries. The orbital elements of the disc orbits are taken from our simulations, with the exception of the inclination and longitude of the ascending node, which are fixed to the observed mean values for the disc in the GC (130° and 96°, respectively; Yelda et al. 2014). We assume the binaries consist of two |$10\, {\rm M}_{\odot }$| stars on circular orbits of radius 0.05 or 0.5 au.

  In the latter case, the disc eccentricity distribution is virtually unaffected. In the former case, the mean disc eccentricity of bound stars for Cloud 1 increases from ∼0.38 to ∼0.5. Also, in this case, |$\sim 10{{\ \rm per\ cent}}$| of the stars appear to be unbound. Such a large population of unbound stars is not observed.

• We do not include the effects of discrete stars and compact objects in the nuclear star cluster surrounding the disc. As shown by Perets et al. (2018), this can have a dramatic impact on its evolution, leading to flips in orientation, transient spiral arms, and warping.

• In Section 4.2, we compared the two-body relaxation time in our simulations to that expected with a realistic mass function. In that comparison we assumed a mass function extending to 100 M_⊙. In reality, stars above |$\sim 40\, {\rm M}_{\odot }$| would have died over 5 Myr, increasing the two-body relaxation time by a factor of ∼3.

Finally, Madigan (2011) previously found that the slope of the background stellar density can affect the number of disruptions coming from the disc. We explore this in the next section.

5.2 Slope of the stellar density profile

We have assumed a relatively flat density profile with a power-law index of –1.16, based on Schödel et al.’s (2018) fits to the diffuse stellar light profile in the GC. In fact, the density profile may be significantly steeper, especially on the scales of interest (∼0.1 pc).

Fits to the resolved stellar population find a steeper profile close to r^−1.4 (Gallego-Cano et al. 2020). In fact, the diffuse light may be probing a somewhat different, younger stellar population that has not yet relaxed dynamically (see the discussion in section 8.3 of Schödel et al. 2020).

Theoretically, a relaxed stellar-density profile is expected to fall between r^−1.5 and r^−1.75 (Bahcall & Wolf 1976; Alexander & Hopman 2009). In fact, this relaxed cusp would only develop a factor of ∼5 times inside of the influence radius (or 0.4 pc in the GC – see Vasiliev 2017). Thus, even if the power-law index of the density profile is −1.16 near 1 pc, it would be significantly steeper on small scales for a relaxed population. Finally, the precession rate of the disc depends on the total extended mass profile, including stellar mass black holes, which are expected to accumulate in the central region via mass segregation (Morris 1993; Miralda-Escudé & Gould 2000; Freitag, Amaro-Seoane & Kalogera 2006; Hopman & Alexander 2006; Merritt 2010). This population can form a very steep cusp. For example, previous work has found that the black holes can have r⁻² over a large range of radii (Freitag et al. 2006; Alexander & Hopman 2009; Aharon & Perets 2016; Vasiliev 2017). We note that a significant population of stellar mass black holes is likely required to reproduce the eccentricity distribution of the S-stars (Antonini & Merritt 2013; Generozov & Madigan 2020). Furthermore, observations of X-ray binaries in the GC also suggest the presence of a black hole cusp (Muno et al. 2005; Generozov et al. 2018; Hailey et al. 2018).

To test the effects of a steeper density profile, we have added an additional external potential into (low resolution) simulations of Cloud 1 and Cloud 2, representing a black hole cusp. Specifically, we assume the black holes have an r⁻² density profile, with mass |$M(r)=2.4\times 10^4\, {\rm M}_{\odot } \left(r/0.1 {\rm pc}\right)$| inside of radius r, as in Antonini & Merritt (2013).⁸ The number of mean disruptions from Cloud 1 (Cloud 2) is 7.8 (5.6) with this steeper density profile. Thus, this black hole cusp does not strongly affect the results. We also tried a steeper, single component cusp with an r^−1.5 density profile and the same stellar mass inside of 1 pc (⁠|$\sim 1.2\times 10^6\, {\rm M}_{\odot }$|⁠). The mean number of disruption for Clouds 1 and 2 are 8.6 and 6.

Previously, Madigan (2011) found that the slope of the background density profile strongly affects the disc instability that pushes binaries to disruption. In order for disruptions to occur, the outer disc has to precess ahead of the inner disc. As the stellar density becomes flatter, differential precession and disruptions are suppressed.

However, in our case, the discs start with a strong primordial twist as illustrated in Fig. 12. For both Cloud 1 and Cloud 2, the bulk of the outer disc starts significantly ahead in its precession. This twist allows the eccentric disc instability to develop even in the presence of a steep cusp.

Figure 12.

Initial orientation of stellar orbits in Cloud 1 (top panel) and Cloud 2 (bottom panel) as a function of semimajor axis. Cloud 1 (2) has its angular momentum along the negative (positive) z-axis, and precesses towards increasing (decreasing) ϖ. Thus, the disc will precess towards the right (left) in the top (bottom panel) panel.

We recover a bimodal eccentricity distribution for Cloud 2 with the r^−1.5 background density profile. However, the bimodal structure is still absent for Cloud 1. This suggests that a steep background density profile is a necessary, but insufficient condition, for producing a bimodal eccentricity at late times. Cloud 1 produces the most massive stellar disc, with the fastest two-body relaxation. This enhanced two-body relaxation may be suppressing the peak at high eccentricities, by causing stellar orbits to thermalize.

The numerical experiments described in this subsection are not self-consistent. Although we changed the background stellar density in our N-body simulations, we are still using snapshots from hydrodynamic simulations with a flatter background density profile as initial conditions. Thus, the twist in Fig. 12 may be an artefact of strong differential precession induced by this flatter density profile. However, the initial conditions here are somewhat arbitrary to begin with, as discussed in Section 3.3. It is possible that a strong primordial twist could also be seeded by a non-uniform cloud, even with a steep background density profile.

5.3 Single star disruptions

We can infer the number of single star disruptions from our disc simulations. For a |$10\hbox{-} {\rm M}_{\odot }$| star, the tidal disruption radius is (M_bh/M_*)^1/3R_* ≈ 6.3 × 10⁻⁶ pc. We find that |$\sim 2{-}3{{\ \rm per\ cent}}$| of particles reach a smaller pericentre for both clouds. If the particles are in fact binaries, they would be tidally disrupted before their component stars could reach their tidal disruption radii, since they would be in the empty loss cone regime (Generozov & Madigan 2020). Thus, the single star disruption rate is |$\sim 2{{\ \rm per\ cent}} (1-f_{\rm bin})$|⁠, where f_bin is the binary fraction in the disc.

5.4 Properties of ejected stars

Recently Koposov et al. (2020) discovered a HVS that was ejected from the GC 5 Myr ago. Such stars would be produced by our discs, since binary disruptions typically eject one star from each disrupted binary from the GC (Hills 1988).

As shown in Fig. 13, most disruptions occur within the first Myr for both clouds, though there is a tail extending to late times. We follow the procedure described in Generozov (2020) to predict the angular distribution of HVSs produced by each cloud. In particular, for each disruptee, we initialize a hyperbolic orbit with the same orbital angles and pericentre, and assume a fixed semimajor axis of –0.01 pc.⁹ We then measure the velocity at infinity for this orbit. Fig. 14 shows an aitoff projection of the polar angles for these velocities. Stars ejected by Cloud 1 show weak clustering in azimuth. This will lead to a concentration of young HVSs in a particular region of the sky. For stars ejected within first 1 (5) Myr, the p-value that the azimuthal angle is drawn from a uniform distribution is 3 × 10⁻⁴ (0.002) according to the Kolmogorov–Smirnov test. The circular standard deviation of the azimuthal angle is 47 (69) degrees. For comparison, in the simulations of Generozov & Madigan (2020) and Generozov (2020), the standard deviation of stars ejected in first Myr was 30º−60º.

Figure 13.

Distribution of disruption times for N-body disc simulations initialized from Cloud 1 (top panel) and Cloud 2 (bottom panel). These results are from our high-resolution simulations.

Figure 14.

Aitoff projection of polar angles of ‘ejected stars’ velocity vectors. Here, the equator corresponds to the disc mid-plane. The top (bottom panel) panel corresponds to N-body disc simulations initialized from Cloud 1 (Cloud 2). Each point is coloured by the time disruption (and ejection) would occur.

The production rate of HVSs may be dominated by other channels like scattering by molecular clouds from large scales (Perets, Hopman & Alexander 2007; Perets & Gualandris 2010). Alternatively, some background stars in the Milk Way’s nuclear star cluster may be excited to disruption by the disc. In fact, S5-HVS1 could be such a star, considering its spectroscopic age (25–93 Myr; Koposov et al. 2020) exceeds the disc’s

5.5 Differences with previous work

Gualandris et al. (2012) performed high accuracy N-body simulations initialized from hydrodynamical simulations of molecular cloud infall and star formation by Mapelli et al. (2012). In many regards, our work is similar to that of the quoted authors, but results do bear quantitative differences. In particular, Mapelli et al. (2012) found discs that are significantly less eccentric (e ∼ 0.2−0.4). In contrast, the mean eccentricity of the discs at the start of our N-body simulations is ∼0.5 (with a maximum eccentricity up to 0.99).

The disparity in the initial eccentricity distribution may be driven by differences in the initial conditions, physical problem setup, and numerical techniques. In particular, the biggest differences are as follows:

• The Mapelli et al. (2012) gas cloud is 1.9 times larger; it starts at a galactocentric distances of 25 pc, whereas we focus on clouds that start at 10 pc. The initial gas velocity is the same everywhere in our clouds while the former authors seed their clouds with supersonic turbulence.

• The background stellar density profile is different. Mapelli et al. (2012) assume an r⁻² (r^−1.4) density profile outside (inside) of 0.39 pc, whereas, we use the broken power law in equation (1).

• Our clouds are on less radial orbits. The initial pericentre of the centre of our clouds is ∼0.2 pc. In Mapelli et al. (2012), the impact parameter of the cloud centre is 0.01 pc (and the pericentre would be even smaller).

• The thermodynamics of gas in the simulations is different. We impose a minimum (irradiation) temperature of 15 K, whereas Mapelli et al. (2012) impose a much higher minimum temperature of 100 K. These temperatures are all well below the local virial temperature, ∼10⁶−10⁷ K, so these choices are probably immaterial for the large-scale cloud infall dynamics but may affect the dynamics of gas and the stellar mass function in the disc. Radiative cooling is included in both our hydrodynamic simulations and the one used by Gualandris et al. (2012).

• The artificial viscosity of the SPH simulations in the two studies is different. The latter detail may affect shocks that form during the cloud infall and subsequent disc formation.

• Another significant difference is the prescription for star formation. While we introduce sink particles that are allowed to grow by further gas accretion, Mapelli et al. (2012) do not introduce sink particles.

Currently, it is not possible to pinpoint which of these factors drives the differences in our results.

Overall, more simulation work is needed to ascertain the robustness of our (and Gualandris et al. 2012) conclusions, with attention to both numerics and the uncertainties in the stellar potential and the initial configuration of the molecular cloud prior to its infall in the GC. In particular, the effect of the initial cloud orbit warrants attention, as counterintuitively we obtained a more eccentric disc with a less radial orbit.

6 SUMMARY

In this paper, we have simulated disruptions of molecular clouds in the GC. Such disruptions can lead to the formation of the lopsided, eccentric discs that would explain the young disc observed within the central parsec of the Galaxy. Moreover, these discs are subject to a secular gravitational instability that can excite binary stars in the disc to extreme eccentricities and disruption. The remnants from these disruptions would explain the S-star cluster within ∼0.04 pc of the GC. Our main conclusions are summarized as follows:

• The velocity of the initial cloud can be constrained from observations. Clouds with velocities ≳100 km s⁻¹ at 10 pc form long star-forming filaments that are not observed in the GC. Therefore, such high velocities can be ruled out. Additionally, velocities ≲ 20 km s⁻¹ can be ruled out, as they produce discs that are too small in radial extent to explain the observed clockwise disc. The most promising cloud velocities are ∼40−80 km s⁻¹.

  Similarly, the mass of the cloud can be constrained. Clouds with masses between a few|$\times 10^4$| and |$10^5\, {\rm M}_{\odot }$| are the most promising. Lower mass clouds produce too few stars to explain the observations. Higher mass clouds would likely produce too many.

• We identified two promising cloud initial conditions for reproducing the young clockwise disc in the GC, as well as the S-stars. Overall, these clouds can push |$\sim 6{-}8{{\ \rm per\ cent}}$| of their constituent binaries to disruption.

• We are able to reproduce the mean eccentricity of the clockwise disc. The final eccentricity profiles in our simulations is unimodal, although a bimodal distribution is present at early times. Previous work (Madigan et al. 2009; Generozov & Madigan 2020) found a bimodal eccentricity distribution due to a steeper background density profile.

• The eccentricity distribution is inclination dependent. For example, even if a bimodal eccentricity distribution is present, it will only show up if high inclination stars are included in the distribution. As such stars are not identified as disc stars in observations, we would not expect a bimodal eccentricity in the observed clockwise disc.

• We are not able to reproduce the observed surface density profile of the clockwise disc. In particular, the surface density falls off too steeply with radius.

• There are a number of simplifications we have made that may affect the aforementioned observational comparisons. In particular, our N-body simulations lack a mass spectrum, a live nuclear star cluster, and gas.

We have demonstrated that the tidal disruption of a molecular cloud is a viable way to explain the young disc stars in the GC, as well as the S-stars and HVSs.

The eccentricity distribution of non-disc S-stars can be used as a test of our model. In particular, we predict young stars outside of the disc should have higher average eccentricities. Additionally, the presence or absence of counterrotating discs as claimed in Genzel et al. (2003), Paumard et al. (2006) (though see Lu et al. 2009; Yelda et al. 2014) can be used as an additional constraint on the initial conditions.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments.

We thank Andreas Eckart, Stefan Gillessen, Yuri Levin, Avi Loeb, Hagai Perets, and Rainer Schödel for helpful comments.

AM gratefully acknowledges support from the David and Lucile Packard Foundation. This work utilized resources from the University of Colorado Boulder Research Computing Group, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. SN acknowledge the funding from the UK Science and Technologies Facilities Council, grant No. ST/S000453/1. This work made use of the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

DATA AVAILABILITY

Data used in this study will be made available upon reasonable request.

Footnotes

1

Dynamical friction and gas accretion on to stars are relatively weak compared to gravitational forces acting on stars in the central parsec.

2

In practice, we first select a random set of 150 stars from these snapshots, delete all binaries, and then delete more stars until the total remaining is 100.

3

This is approximately the tidal radius of a 20-M_⊙ binary with a semimajor axis of 1 au.

4

Here, we use osculating orbital elements about the central SMBH.

5

To accurately compute the relaxation time of a disc it is necessary to account for its flattened geometry, as done in this reference.

6

Interestingly, the combination N|$\langle$|m²|$\rangle$| is weak function of the minimum mass for an m^−1.7 mass function, changing by ∼15 per cent as the minimum of the mass function is decreased from 1 to |$0.1\, {\rm M}_{\odot }$| (at fixed disc mass).

7

In this case, filtering out the high inclination stars results in a flat eccentricity distribution.

8

This model can reproduce the S-star eccentricity distribution within these stars’ lifetimes. However, assuming a single power-law index for the black hole profile is likely unrealistic: as black holes begin to dominate relaxation on small scales, their density profile should flatten (Vasiliev 2017).

9

The orbital geometry depends on the ratio of pericentre to semimajor axis, which is approximately (M_bh/m_bin)^−1/3. This would have a broader range than assumed here. However, a realistic spread in this ratio does not significantly affect the angular distribution in Fig. 14.