Infalling Young Clusters in the Galactic Centre: implications for IMBHs and young stellar populations
J. A. Petts A. Gualandris




firstpage–lastpage 2015

Unifying microscopic and continuum treatments of van der Waals and Casimir interactions
    Alejandro W. Rodriguez
    Submitted TEMP
=======================================================================================







The central parsec of the Milky Way hosts two puzzlingly young stellar populations, a tight isotropic distribution of B stars around SgrA* (the S-stars) and a disk of OB stars extending to ∼ 0.5pc. Using a modified version of Sverre Aarseth's direct summation code NBODY6 we explore the scenario in which a young star cluster migrates to the Galactic Centre within the lifetime of the OB disk population via dynamical friction. We find that star clusters massive and dense enough to reach the central parsec form a very massive star via physical collisions on a mass segregation timescale. We follow the evolution of the merger product using the most up to date, yet conservative, mass loss recipes for very massive stars. Over a large range of initial conditions, we find that the very massive star expels most of its mass via a strong stellar wind, eventually collapsing to form a black hole of mass ∼ 20 - 400 M⊙, incapable of bringing massive stars to the Galactic Centre. No massive intermediate mass black hole can form in this scenario. The presence of a star cluster in the central ∼ 10pc within the last 15 Myr would also leave a ∼ 2 pc ring of massive stars, which is not currently observed. Thus, we conclude that the star cluster migration model is highly unlikely to be the origin of either young population, and in-situ formation models or binary disruptions are favoured.



stars: winds – stars: evolution – stars: black holes – stars: massive




§ INTRODUCTION




The central parsec of the Milky Way hosts almost two dozen He-1 emission-line stars <cit.> and a population of many other OB stars in a thin clockwise disk extending from ∼0.04-1.0pc <cit.>. <cit.> recently extended the range of observations up to ∼ 4pc^2 centred on SgrA*; the radio source associated with the supermassive black hole (SMBH) at the centre of the Milky Way. The authors show that the OB population is very centrally concentrated, with 90% projected within the central 0.5pc. The clockwise disk exhibits a top heavy mass function (α∼ 1.7, <cit.>). <cit.> estimate the He-1 stars to be only ∼3-7 Myr old, which is puzzling as the tremendous tidal forces in this region make it difficult for a giant molecular cloud (GMC) to remain bound long enough for gas to cool and fragment <cit.>.

There appears to be very few He-1 stars farther than the central parsec, other than inside/near the young Arches <cit.> and Quintuplet <cit.> clusters at ∼ 30 pc. This led <cit.> to postulate that efficient dynamical friction on star clusters forming a few parsecs from SgrA*, where GMCs can more easily cool and fragment, could bring a dense core of massive stars to the central parsec within the age of the He-1 population.

Another model suggests that in-situ formation of the clockwise disk is possible if a tidally disrupted GMC spirals in to form a small gaseous disk, which can be dense enough to become Jeans unstable and fragment into stars <cit.>. The infalling cloud needs to be ∼ 10^5 M⊙ in order to reproduce observations <cit.>. Two large gas clouds of mass ∼5×10^5M⊙, M-0.02-0.07 and M-0.13-0.08, are seen projected at ∼ 7 and ∼ 13 pc from the Galactic Centre, respectively <cit.>. The top heavy mass function can
be reproduced by the in-situ model so long as the gas has a temperature greater than 100K, consistent with observations of the Galactic clouds. The rotation axis of the clockwise disk shows a strong transition from the inner to outer edge <cit.>, suggesting that the disk is either strongly warped, or is comprised of a series of stellar streamers with significant variation in their orbital planes <cit.>. In-situ formation is currently favoured for the clockwise disk, as an infalling cluster would likely form a disk with a constant rotation axis <cit.>. A caveat of in-situ formation is that it requires near radial orbits incident upon SgrA*, perhaps requiring cloud-cloud collisions <cit.>. 

Interior to the disk lies a more enigmatic population of B-stars in a spatially isotropic distribution around SgrA*, with a distribution of eccentricities skewed slightly higher than a thermal distribution (<cit.>, <cit.>). These “S-stars” have semi-major axes less than 0.04 pc, with S0-102 having the shortest period of just 11.5±0.3 yrs, and a pericentre approach of just ∼260 AU <cit.>. The S-star population could potentially be older than the disk population, as the brightest star in this population, S0-2, is a main sequence B0-B2.5 V star with an age less than 15 Myr <cit.>. The other stars in this population have spectra consistent with main sequence stars <cit.>, and observational limits require them to be less than 20 Myr old in order to be visible.

The tidal forces in this region prohibit standard star formation, so the S-stars must have formed farther out and migrated inwards. A possible formation mechanism of the S-stars is from the tidal disruption of binaries scattered to low angular momentum orbits, producing an S-star and a hyper-velocity star via the Hills mechanism <cit.>. The captured stars would have initial eccentricities greater than 0.97 <cit.>, but the presence of a cusp of stellar mass black holes around SgrA* could efficiently reduce the eccentricities of these orbits via resonant relaxation within the lifetime of the stars <cit.>. Additionally, <cit.> show that if a binary is not tidally disrupted at first pericentre passage, the Kozai-Lidov (KL) resonance <cit.> can cause the binary to coalesce after a few orbital periods, producing an S-star and no hyper velocity star.

Alternatively, <cit.> show that stars from the clockwise disk can be brought very close to SgrA* via global KL like resonances, if the clockwise disk of gas originally extended down to ∼ 10^-6pc (the lowest stable circular orbit around SgrA*). The authors also show that O/WR stars would be tidally disrupted within the region of the observed S-star cluster due to their large stellar radii, whereas B-stars could survive, in agreement with observations. Recently, <cit.> showed that a clockwise disk with 100% primordial binarity can produce ∼ 20 S-stars in less than 4Myr. KL oscillations can efficiently drag binaries close to SgrA*, producing an S-star and a hyper-velocity star. This mechanism produces S-stars with eccentricities lower than from the disruption of binaries originating from outside the disk. However, in order to thermalize the S-stars, ∼ 500 M⊙ in dark remnants are still required around SgrA* in order to match observations, consistent with Fokker-Planck models <cit.>. Three confirmed eclipsing binaries are observed within the clockwise disk, all being very massive O/WR binaries <cit.>. <cit.> estimate the present day binary fraction of the disk to be 0.3^+0.34_-0.21 at 95% confidence, with a fraction greater than 0.85 ruled out at 99% confidence. More recently, <cit.> predict that the binary fraction must be greater than 32% at 90% confidence.

An additional popular scenario is the transport of stars from young dense star clusters that migrate to the Galactic Centre via dynamical friction, with the aid of an intermediate mass black hole (IMBH). <cit.> showed that to survive to the central parsec from a distance ≥ 10pc, clusters either need to be very massive (∼10^6M⊙) or very dense (central density, ρc∼ 10^8 M⊙ pc^-3). <cit.> showed that including an IMBH in the cluster means the the core density can be lowered, but only if the IMBH contains ∼ 10 % of the mass of the entire cluster, far greater than is expected from runaway collisions <cit.>.

<cit.> (hereafter F09) revisited this problem using the tree-direct hybrid code, BRIDGE  <cit.>, allowing the internal dynamics of the star clusters to be resolved. The small tidal limits imposed by SgrA* meant the clusters had core densities greater than 10^7 M⊙pc^-3, leading to runaway collisions on a mass segregation timescale <cit.>. During collisions, the resulting very massive star (VMS) was rejuvenated using the formalism of <cit.>, and collapsed to an IMBH at the end of its main sequence lifetime, extrapolated from the results of <cit.>. The authors found that by allowing the formation of a 3-16 × 10^3 M⊙ IMBH <cit.>, some stars could be carried very close to SgrA* via a 1:1 mean resonance with the infalling IMBH. The orbits of these “Trojan stars” were randomised by 3-body interactions with the SMBH and IMBH, constructing a spatially isotropic S-star cluster. F09's simulation “LD64k” transported 23 stars to the central 0.1pc, however, the resolution of the simulation is ∼ 0.2pc, set by the force softening of SgrA*. The simulation also brought 354 stars within 0.5pc of SgrA*, 16 being more massive than 20M⊙, analogous to clockwise disk stars. The IMBH formed in LD64k is more massive than the observational upper limit of ∼ 10^4 M⊙, derived from VLBA measurements of SgrA* <cit.>. However, <cit.> state that an IMBH of 1500 M⊙ is sufficient for the randomisation of stars <cit.>.

Despite the successes of the F09 model, IMBH formation in young dense star clusters may be prohibited. VMSs of the order 10^3 M⊙ are expected to have luminosities greater than 10^7 L⊙ <cit.>, driving strong stellar winds. F09 assumed the mass loss rate of stars more massive than 300M⊙ to be linear with mass, however, recent work on VMS winds show steeper relations for stars that approach the Eddington limit <cit.>. F09's model also neglected the effect of the evolving chemical composition on the luminosity, and hence the mass loss, of the VMS <cit.>. We note that the initial mass function (IMF) used in F09, although employed due to numerical constraints, meant there were ten times more massive stars than expected from a full Kroupa IMF, leading to an increased collision rate and buildup of the VMS mass.

No conclusive evidence for the existence of IMBHs in star clusters has yet been found <cit.>. Sufficiently high mass loss could cause VMSs to end their lives as stellar mass black holes or pair-instability supernovae at low metallicity <cit.>. Pair-instability supernovae candidates have recently been found at metallicities as high as ∼ 0.1 Z⊙ <cit.>, with expected progenitors of several hundred solar masses <cit.>.

The most massive star observed, R136a1, is a 265^+80_-35M⊙ star in the 30 Doradus region of the Large Magellanic Cloud (LMC) <cit.>, with metallicity Z = 0.43 Z⊙. <cit.> suggest that it could be a very rare main sequence star, with a zero age main sequence mass of 320^+100_-40M⊙. However, it could be the collision product of a few massive stars. R136a1 has a large inferred mass loss rate of (5.1^+0.9_-0.8)× 10^-5M⊙yr^-1, ∼ 0.1 dex larger than the theoretical predictions of <cit.>. <cit.> predict that the evolution of all stars more massive than 300 M⊙ is dominated by stellar winds, with similar lifetimes of ∼2-3Myr. As such, it is not surprising that R136a1 is the most massive star currently observed, as more massive VMSs should be rare and short lived.

Whilst it may be unlikely for an IMBH to form at solar metallicity, a VMS could transport stars to SgrA* within its lifetime. In this paper we test the feasibility of the star cluster migration scenario as the origin of either young population in the Galactic Centre.

We evolve direct N-body models of star clusters in the Galactic Centre, using the GPU-accelerated code NBODY6df, a modified version of Sverre Aarseth's NBODY6 <cit.> which includes the effects of dynamical friction semi-analytically <cit.>. In section <ref> we describe the theory behind our dynamical friction and stellar evolution models. In section <ref> we describe the numerical implementation. Section <ref> discusses prior constraints on the initial conditions and describes the parameters of the simulations performed. In sections <ref> and <ref>, we present our results and discuss their implications for the origin of the young populations. Finally, we present our conclusions in section <ref>.



§ THEORY





 §.§ Dynamical friction


The dynamical friction model used in this paper is a semi-analytic implementation of Chandrasekhar's dynamical friction <cit.>, described in <cit.>, which provides an accurate description of the drag force on star clusters orbiting in analytic spherical background distributions of asymptotic inner slope γ = 0.5,..,3 <cit.>. The novelty of our model is the use of physically motivated, radially varying maximum and minimum impact parameters (bmax and bmin respectively), which vary based on the local properties of the background. The dynamical friction force is given by:


    dvcl/dt = -4π G^2 Mclρlog(Λ) f(v* < vcl) vcl/v^3cl

where vcl is the cluster velocity, Mcl is the cluster mass, ρ is the local background density and f(v* < vcl) is the fraction of stars moving slower than the cluster; assuming a Maxwellian distribution of velocities, valid in the cuspy models explored here. The Coulomb logarithm is given by:


    log(Λ) = log(bmax/bmin) = log(min(ρ(Rg)/|∇ρ(Rg)|,Rg)/max(rhm, GMcl/vcl^2)),

where Rg is the galactocentric distance of the cluster and rhm is the half mass radius of the cluster. When coupled with the N-body dynamics, rhm is the live half mass radius, and  Mcl is well represented by the cluster mass enclosed within its tidal radius, including stars with energies above the escape energy.



 §.§ Evolution of very massive stars



<cit.> present similarity theory models of VMSs, for which the stellar properties can be calculated by solving a set of differential equations <cit.>. VMSs are predicted to have large convective cores containing more than 85% of the mass, surrounded by a thin extensive radiative envelope. In such stars the opacity becomes larger than the electron scattering value, and can be considered to come from Thomson scattering alone. Utilising such approximations, the authors  provide simple formulae to calculate the core mass and luminosity, as functions of stellar mass and chemical composition.

The luminosity of stars with μ^2M ≥ 100 can be found by substituting equation 36 of <cit.> into their equation 34:


    L ≈64826 M (1 - 4.5/√(μ^2 M))/1 + X,

where L is the luminosity, M is the mass of the VMS, X is the core hydrogen abundance, and μ is the mean atomic mass of the core. Assuming a fully ionised plasma, μ takes the form:

    μ = 4/6X + Y + 2,

where Y is the core helium abundance. Equation <ref> shows that at very large masses L ∝ M. However, unlike the F09 model, this formulation of the luminosity explicitly includes an L ∝ (1+X)^-1 dependence. As the mass loss rate depends on L, this leads to an increased mass loss rate in the late stages of evolution (see section <ref>).

<cit.> (hereafter B07) modelled the evolution of VMSs with zero age main sequence (ZAMS) masses of up to 1000M⊙, assumed to have formed via runaway collisions in a young dense star cluster. The authors numerically evolve the chemical composition of the star through the Core Hydrogen Burning (CHB) and Core Helium Burning (CHeB) phases via conservation of energy and mass loss from the stellar wind. In this section we briefly outline the model of <cit.> and describe how we include stellar collisions and their effect on VMS evolution.

As VMSs have large convective cores, one can reasonably approximate them as homogeneous (verified to be a good approximation down to 120 M⊙, B07). Applying conservation of energy, the hydrogen fraction in the core during CHB evolves as equation 1 of B07:


    Mcc(μ,M)dX/dt = - L(μ,M)/ϵH,

where Mcc is the mass of the convective core and ϵH is the hydrogen burning efficiency (i.e. the energy released by fusing one mass unit of hydrogen to helium).

When the core is depleted of hydrogen, the VMS burns helium via equation 4 of B07 (see also <cit.>):



    Mcc(μ,M) dY/dt = -L(μ,M)/ϵratio,
    ϵratio = [(BY/AY - BO/AO)+(BC/AC - BO/AO) C'(Y)],


where ϵratio accounts for the fact that C and O are produced in a non-constant ratio, affecting the energy production per unit mass of helium burnt. Here, A and B are the atomic weights and binding energy of nuclei; with subscripts Y, C and O representing helium, carbon and oxygen respectively. C'(Y) is the derivative of the C(Y) fit from <cit.> with respect to Y (see B07 for the derivation of Equation <ref>). During CHeB, μ is defined as <cit.>:


    μ = 48/36Y + 28C + 27O,
    which by assuming Y+C+O=1 and using the fit to C(Y) by <cit.>, can be rewritten solely as a function of Y as:μ = 48/19Y + C(Y) + 27.


Subsequent stages of evolution are rapid and explosive. We assume that after core helium burning the remnant collapses to a black hole with no significant mass loss.



  §.§.§ Mass loss


The chemical evolution of the VMS is coupled to the mass evolution, as the luminosity of the star sets the wind strength. <cit.> (hereafter V11) show that the wind strength is heavily dependent on the proximity to the Eddington limit, when gravity is completely counterbalanced by the radiative forces, i.e. grad/ggrav = 1, where grad and ggrav are the radiative and gravitational forces, respectively. For a fully ionised plasma, the Eddington parameter, Γe, is dominated by free electrons and is approximately constant throughout the star (V11):


    Γe = grad/ggrav = 10^-4.813 (1 + Xs) (L/L⊙) (M/M⊙)^-1,

where Xs is the surface hydrogen abundance of the star. V11's fig. 2 shows that the logarithmic difference between the empirical <cit.> (here after V01) rates and the VMS rates follow a tight relation with Γe, almost independent of mass. The authors find that the mass loss rate is proportional to:


    Ṁ∝Γe^2.2,    if  0.4< Γe < 0.7
    Γe^4.77,    if  0.7<Γe < 0.95.


During the CHB phase, we model the stellar wind of the VMS using the formulae from V01, whilst correcting for the proximity to the Eddington limit by fitting on the data from table 1 of V11. In this way we obtain a coefficient that allows us to convert the V01 rate to the Γe enhanced rates of stars approaching the Eddington limit <cit.>. V11 modelled stars up to 300 M⊙, however, as the logarithmic difference between the V11 and V01 rates shows little dependence on mass, we extrapolate this approach to higher masses. V11 state that their predicted wind velocities are a factor 2–4 less than derived empirically. The effect of rotation is also neglected. It should be noted that due to these two effects, and our extrapolation of the V11 models, we likely underestimate the mass loss of our VMSs. Therefore the masses of our VMSs and their resulting remnants should be taken as a conservative upper limit at solar metallicity.

During CHeB VMSs are depleted of hydrogen and are expected to show Wolf-Rayet like features. We follow the approach of <cit.> and extrapolate the mass loss formula of <cit.>:

    log(Ṁ) = -11  + 1.29 log(L) +1.7 log(Y) + 0.5 log(Z).


B07 explored models with Wolf-Rayet like mass loss rates (arbitrarily) up to 4 times weaker, which only left a remnant twice as massive. The uncertainty arising from extrapolation of this formula should be of little significance to the transport of young stars to the central parsec, as post main sequence VMSs are not massive enough to experience substantial dynamical friction after the cluster is disrupted (B07). However, if sufficiently chemically rejuvenated, a CHB VMS may be capable of bringing stars to the central parsec before losing most of its mass. Thus the evolution during the CHB stage is of most interest.

We make sure that in both burning phases the predicted mass loss never exceeds the photon tiring limit, the maximum mass loss rate that can theoretically be achieved using 100% of the stars luminosity to drive the wind <cit.>:

    Ṁtir = 0.032 (L/10^6L⊙) (R/R⊙) (M/M⊙)^-1.

Here, the radii, R, of stars are taken from the mass-radius relation of <cit.>, which is in excellent agreement with <cit.>'s similarity theory models of VMSs, but requires less computational resources to calculate. The OB disk population is less than 7Myr old. Hence, we assume approximately solar abundances such that X0=0.7, Y0=0.28 and Z0=0.2 <cit.>.



  §.§.§ Rejuvenation following collisions


<cit.> show that VMSs have nearly all of their mass in their large convective cores. Repeated collisions can efficiently mix the core and the halo, keeping the star relatively homogeneous. The wind of the VMS also ensures homogeneity, as the loose radiative envelope is shedded by the stellar wind, leaving the surface with composition similar to the core.

We chemically rejuvenate a VMS following a collision with another star. We assume that stars colliding with the VMS efficiently mix with the convective core such that:


    Xnew = XstarMstar + XVMSMVMS/MVMS+star

Similarly for Y and Z. We approximate Xstar(t) and Ystar(t) for main sequence stars by interpolating the detailed stellar models of <cit.>. If a CHeB VMS collides with a hydrogen rich main sequence star, we assume that CHB is reignited. When two VMSs collide their composition is also assumed to be well mixed.



§ NUMERICAL METHOD



To model the effects of dynamical friction on self-consistent star cluster models we use the GPU-enabled direct N-body code NBODY6df <cit.>, which is a modified version of Aarseth's direct N-body code NBODY6 <cit.>. In this paper we model the background as an analytic stellar distribution with a central black hole (see section <ref>). In <cit.> we only tested our dynamical friction model for cases without a central black hole, however we discuss how to trivially add a black hole to the model in Appendix A. A validation of this approach via comparison with full N-body models computed with GADGET <cit.> is also given in the appendix.

We introduced an additional modification to the code to properly model the evolution of a VMS, as described in section <ref>. When a physical collision creates a star greater than 100M⊙ we flag it as a VMS and treat its evolution separately from the standard SSE package in NBODY6 <cit.> via the method described in section <ref>. As the mass loss can be very large for VMSs, fine time resolution is needed to prevent overestimation of the mass loss. We introduce a new routine which integrates the mass and composition of the star between the dynamical time steps using a time step of 0.1 years, sufficiently accurate to resolve the evolution. V11 predict terminal wind speeds of a few thousand kms^-1 for VMSs, as such we assume that the stellar wind escapes the cluster and simply remove this mass from the VMS. An arbitrary number of VMSs can potentially form and evolve simultaneously in the simulation.



§ INITIAL CONDITIONS







In NBODY6df the background potential is assumed static and analytic; an assumption valid over the short timescales considered here (less than 7Myr). We adopt a Dehnen model <cit.>, representing the central region of the Galaxy. We use a slope γ = 1.5, scale radius a = 8.625pc and total mass Mg = 5.9× 10^7 M_⊙, which closely reproduces the observed broken power-law profile obtained by <cit.> for the central region of the Galaxy, yet has simple analytic properties. We place a central fixed point mass of 4.3×10^6M⊙ to represent SgrA* <cit.>.



 §.§ Physical and numerical constraints on the initial conditions


There are two constraints on the initial conditions of the clusters. Firstly, they must reach the Galactic Centre within the age of the young populations. We therefore wish to model clusters that can potentially reach the Galactic Centre in less than 7Myr, so that we may test the migration model for both the clockwise disk and the S-stars. We obtain tight constraints on the initial orbital parameters by integrating the orbits of point masses in the Galactic Centre potential including dynamical friction. Fig. <ref> shows contours of equal inspiral time for different initial masses, apocentres and initial velocities. Initial conditions to the right of each line are such that the clusters can reach within 0.5pc in less than 7Myr. Arches like clusters <cit.>) could reach the Galactic centre in less than 7Myr if they formed at ∼5pc, or from 7-10pc if large initial eccentricities were assumed. More massive clusters can easily migrate ∼ 10pc in 7Myr. We note that these inspiral times are lower limits, as real clusters would lose mass from stellar winds and tides. We choose to model only those clusters for which a point mass object of the same mass can reach the Galactic Centre within ∼7Myr.

Secondly, the size of the clusters is limited by their small tidal limits when so close to SgrA*. Approximating the cluster as a point mass, the tidal radius is given by <cit.>:

    r^3t = GMcl/ω^2p +(d^2Φ/dR^2)p,

where ωp and (d^2Φ/dR^2)p are the angular velocity of a circular orbit and the second derivative of the potential at pericentre, respectively. The high mass requirement for fast inspiral, coupled with the small tidal limits, means that all models are inherently very dense and runaway mergers are expected. Although it is unknown whether such dense clusters are likely to form in the Galactic Centre, we explore these initial conditions in order to test the feasibility of the inspiral model.







 §.§ Initial Mass function





We sample stars from a Kroupa initial mass function (IMF) with an upper limit of 100 M⊙ <cit.>. A lower mass limit of 0.08 M⊙ would yield the most physically realistic results, but at a computational cost unfeasible for a parameter study of such massive clusters at the current time  (365k-730k particles for the most massive models explored). However, truncating the low end of the IMF means that one samples too many massive stars as compared with a full Kroupa IMF. To quantify the difference this has on VMS formation, we ran three test simulations at different mass resolutions, in the absence of a tidal field. Simulations 1lo, 1hi and 1kr have lower cutoffs of 1.0, 0.16 and 0.08 M⊙, respectively. We model star clusters as King models with dimensionless central potential, W0 = 6, and with no primordial mass segregation. The parameters of the isolated simulations are displayed in Table <ref>. Fig. <ref> shows the VMS mass as a function of time for simulations 1lo, 1hi and 1kr, showing that better sampling of the low end of the IMF inhibits the growth of the VMS. This occurs because primarily high mass stars build up the VMS, due to their short dynamical friction timescales and large cross sections for collision. In simulation 1kr, although half the cluster mass is comprised of stars less massive than 0.58 M⊙, only 37 stars less massive than 0.58 M⊙ are consumed throughout the entire lifetime of the VMS. The VMS initially grows very rapidly. However, the late main sequence evolution is dominated by the strong stellar wind of the helium rich VMS. Throughout its lifetime, the VMS in simulation 1kr removes 2244M⊙ of material from the cluster through its stellar wind, ∼ 2% of the cluster mass. During CHeB, simulations 1lo and 1hi reignite CHB via collision with a massive main sequence star, resulting in a lower remnant mass at collapse. The late evolution is very stochastic, however this is not important for the migration of young stars to the Galactic Centre, as the VMS only provides gravitational binding energy comparable to normal cluster stars during its CHeB phase.

Fig. <ref> shows that a lower limit of 0.16 M⊙ is sufficient to resolve the mass evolution of the VMS, and as we are only interested in the final distribution of OB stars, this IMF is sufficient for our simulations. We cannot evolve the most massive clusters at high mass resolution, as these models become too computationally expensive. As a compromise, we test a large range of initial conditions with a lower limit of 1 M⊙, and re-run a selection of initial conditions with a lower limit of 0.16 M⊙ to obtain more realistic results. We can simultaneously use the low resolution simulations to explore the possibility of an initially top heavy mass function for clusters forming close to SgrA*. A very top heavy function is observed for the clockwise disk <cit.>, however, it is unknown whether a top heavy IMF is expected from the collapse of GMCs at ∼5-10 pc from SgrA*.





 §.§ Binary fraction



Some simulations include a population of primordial binaries. Binaries are initialised as follows. Firstly all stars more massive than 5 M⊙ are ordered by mass. The most massive star is then paired with the second most massive star, and so on. This choice is motivated by observational data showing that massive OB stars are more likely to form in binary systems with mass-ratios of order unity <cit.>. Once all stars more massive than 5 M⊙ are in binaries, lower mass stars are paired at random until the specified binary fraction (the fraction of stars initially in a binary system) is reached <cit.>. For stars more massive than 5 M⊙, the periods and eccentricities are drawn from the empirical distributions derived in <cit.>, which show that short periods and low eccentricities are preferred in massive binaries. For lower mass stars the periods are drawn from the <cit.> period distribution and are assigned thermal eccentricities. The mass of a binary and its initial position in the cluster are assumed to be independent.



 §.§ Simulations

The initial conditions are described in Table <ref> and are referred to by the following naming convention: <M><mf><Ra>, where <M> is the cluster mass in units of ∼ 10^5 M⊙, <mf> is the mass resolution of the simulation, and <Ra> is the initial galactocentric distance in pc. For most simulations  we sample from a Kroupa IMF with an upper mass cut off, mup = 100M⊙. The “lo” resolution models have a lower mass cut off, mlow = 1M⊙ and mean mass, ⟨ m*⟩= 3.26M⊙. The “hi” resolution models have mlow = 0.16M⊙ and ⟨ m*⟩= 0.81M⊙. For simulations with <mf> = lu we use an IMF identical to the mass function of the clockwise disk <cit.>. The simulation name is followed by a suffix describing additional information about the simulation. The suffix W4 denotes that the dimensionless central potential, W0, is initially 4 instead of 6. The suffix vX indicates an eccentric orbit with initial velocity, 0.X vc (where vc is the circular velocity at the initial position). The Suffix “ms” indicates that the cluster is primordially mass segregated. Finally, the suffix “b” denotes the inclusion of primordial binaries (see section <ref>). Most models are Roche-filling at first pericentre passage, apart from runs marked with an asterisk, which are Roche-filling at their initial positions. The model with the suffix “d” is extremely Roche under-filling at its initial position.



§ RESULTS





In all models, the clusters are completely tidally disrupted in less than 7Myr. Massive clusters migrate farther in than lower mass clusters on the same initial orbits, due to shorter dynamical friction timescales and less efficient tidal stripping. However any cluster that reaches ∼ 3 pc is rapidly dissolved by its shrinking tidal limit as it approaches SgrA*. Clusters on eccentric orbits inspiral faster, as they pass through denser regions of the cusp periodically. However, clusters on very eccentric orbits (e.g. 2lo10_v2*, e ∼ 0.9) disrupt on the first few pericentre passages, depositing stars at large distances along the initial cluster trajectory.

Most simulations naturally form a VMS in less than 1Myr due to their high initial densities. However, the initial rapid mass accretion soon loses to the increasing mass loss rate and relaxation of the cluster, causing the VMS to collapse to a black hole of ∼ 20 - 250 M⊙ after 2-5 Myr (300-400 M⊙ for models with a <cit.> IMF), typically before their parent clusters completely disrupt. Table <ref> shows the maximum mass, remnant mass and lifetime of the VMS formed in each simulation. The clusters become completely unbound at ∼ 2-3 pc, and the IMBHs formed are not massive enough to experience significant dynamical friction and drag stars close to SgrA* (dynamical friction timescales for even a 400 M⊙ IMBH are longer than 100 Myr). Conversely, the evolution of the VMS does not appear to significantly inhibit the inspiral of the cluster, as only ≲ 2% of the initial cluster mass is typically expelled by the VMS throughout its lifetime.



For each model, Table <ref> shows the final distribution of stars after complete cluster dissolution and death of any VMSs. We show the distributions of semi-major axes for all stars and main sequence stars more massive than 8M⊙ at 7 Myr, as well as how many of these stars have final semi-major axes smaller than 1 pc. We use a 8M⊙ cut-off as these are the faintest main sequence stars spectroscopically observable in the Galactic Centre with current telescopes, K ≥ 15.5 <cit.>. Although photometric studies can see objects down to magnitudes of K<19-18 <cit.>, it is impossible to determine whether these stars are young or old. We also show the projected distributions of visible main sequence stars at 7 Myr and 15 Myr (as the S-star population may be older than the disk population, see section <ref>).



 §.§ Low Resolution Models








Fig. <ref> shows the final distributions of the semi-major axes and projected positions of stars for a representative selection of the models with a lower mass cutoff of 1 M⊙ (<mf> = lo). In all models the final distributions are broad, with a standard deviation of ∼ 2 pc. Other simulations show similar distributions, with less massive and less eccentric models dissolving farther out (see Table <ref>).

Models with very eccentric orbits (e.g. 2lo10_v2*) can bring stars close to SgrA*, however, very few stars have final semi-major axes smaller than 1pc. No stars more massive than 8 M⊙ are scattered to semi-major axes smaller than 1 pc in either 1lo10_v2* or 2lo10_v2*. This is likely due to the preferential loss of low mass stars, whereas high mass stars remain inside the cluster for longer, and end up tracing the final cluster orbit.

Simulation 2lo5 is the only non-radial model to bring stars to the central parsec, and the only model that brings a significant number of massive stars. However, one would expect to also see ∼ 3000 massive stars in the range 1-10 pc, about 10 times more than reach the central parsec. The right side of Fig. <ref> shows the distributions of projected distances of stars that are spectroscopically visible at 7 and 15 Myr. The amplitudes of the distributions are normalised to the expected number of stars had the simulation been run with a Kroupa IMF. The stars are projected to rotate clockwise in the sky. It can be seen that for all simulations, more than 1000 young stars are observed out to ∼ 10 pc. Considering current observational limitations, if a cluster were present in the central ∼ 10 pc within the last ∼ 15 Myr, a large number of stars would be observable up to ∼ 10 pc, suggesting it is unlikely that any clusters have inhabited this region in the last ∼ 15 Myr.

Simulations 2lo10 and 2lo10_W4 have the same initial orbit and mass, yet 2lo10_W4 is less concentrated. The lower density and longer relaxation timescale cause 2lo10_W4 to form a less massive VMS than 2lo10. However, the VMS in 2lo10_W4 lives longer as not all the most massive stars are consumed within ∼ 1 Myr. The models end up with similar final distributions of the resulting disk (see Table <ref>). The same trend is seen for the less massive analogues, 1lo10 and 1lo10_W4. The two most massive simulations 4lo10_W4 and 4lo10_W4v75, are massive enough to reach the central parsec from 15 pc in ∼ 7 Myr, but with central densities low enough to suppress the formation of VMSs. However, these simulations are more susceptible to tides, and are tidally disrupted at large radii.



 §.§ Higher resolution models


Fig. <ref> shows a comparison between simulations 1lo10_v5* and 1hi10_v5*, which have the same initial conditions, except 1hi10_v5* better samples the low mass end of the IMF. The panels on the left show the distributions of semi-major axes for all the stars and main sequence stars more massive than 8M⊙ at 7 Myr. The distributions are very similar, however 1hi10_v5* has a smaller ratio of spectroscopically visible stars to all stars due to differences in the IMF. The panels on the right show the projected distributions of main sequence stars visible at 7 and 15 Myr. For the projected distributions, the number of stars is re-normalised to the expected number of stars had the simulation been run with a Kroupa IMF from 0.08-100 M⊙. Although massive stars are consumed to construct the VMS, this is a small fraction of the population. The distributions look very similar in shape and magnitude, indicating that models run with a lower limit of 1 M⊙ produce similar final distributions to simulations that better sample the IMF. This verifies the validity of the normalisation approach used on the projected visible distributions in Fig. <ref>.

Fig. <ref> demonstrates how simulations 1hi5, 1hi10, 1hi10_v5* and 1hi10_v2* evolve with time. The top two panels show the evolution of the Galactocentric distance of the cluster and the mass enclosed within the tidal radius. The bottom two panels show the evolution of the VMS mass and the half mass radius of the cluster. Simulations 1hi5, 1hi10 and 1hi10_v5* quickly form a VMS and expand due to rapid two body relaxation in the dense core. The expansion lowers the core density and thus the collision rate. The reduced collision rate allows the VMSs to rapidly burn their fuel and collapse without significant hydrogen rejuvenation. Simulation 1hi5 forms a more massive VMS than the other simulations as it is initially ∼ 10 times as dense, however the resulting increased luminosity decreases its lifetime. In simulation 1hi10_v2*, the cluster becomes unbound before the massive stars can reach the centre of the cluster, however the initial density is high enough that a 235 M⊙ VMS forms by the first pericentre passage. The self-limiting nature of the VMS formation is discussed in section <ref>.









 §.§ Models with extreme initial conditions







The young clockwise disk population exhibits a top heavy mass function, with power law index α∼ 1.7 <cit.>. In the context of the cluster inspiral scenario this has been explained by mass segregation inside the cluster, with the most massive stars reaching the central parsec, and low mass stars being preferentially lost due to tides during inspiral. However, as we have shown in section <ref>, clusters lose massive stars as well as low mass stars throughout inspiral, via dynamical ejections and the shrinking tidal limits as the clusters approach SgrA*. In order to test the effect of mass segregation, we ran simulation 1hi5_ms, which we primordially mass segregated using the method described in <cit.>. For simulations 1hi5 and 1hi5_ms, Fig. <ref> shows the semi-major axes of all stars and main sequence stars more massive than 8M⊙ at T=7Myr, as well as the distributions of projected distances of spectroscopically visible stars at 7 and 15Myr. Their distributions look similar, as simulation 1hi5 has an initial dynamical friction timescale of tdf∼ 0.1 Myr for the most massive stars, causing the cluster to rapidly mass segregate. As such, primordial mass segregation does not significantly enhance the transport of massive stars to the central parsec, as clusters of high enough density become mass segregated before tides become important.





As star formation close to a SMBH is not well understood, we also test a model in which the cluster is born with the top heavy mass function derived in <cit.>. Fig. <ref> shows the evolution of simulations 1lo5, 2lo5, 1lu5 and 2lu5, where the two latter clusters have stars sampled from the <cit.> mass function. Models computed with the top-heavy mass function form VMSs of greater mass, as more massive stars are sampled, and their cross sections for collisions are larger. However, as the cluster mass is distributed amongst fewer stars, simulations 1lu5 and 2lu5 relax faster than 1lo5 and 2lo5 and dissolve more rapidly. A flatter mass function does not help bring stars closer to SgrA*, and leaves more visible stars spread across the central 10 pc.

As a final test of extreme initial conditions we re-run simulation 1hi5 with an initial size and density corresponding to being Roche filling at 1 pc. In simulation 1hi5_W4d, this cluster is placed initially on a circular orbit at 5 pc, so that it is initially very Roche under-filling. Fig. <ref> shows the evolution of Galactocentric distance, cluster mass, VMS mass and half mass radius of simulations 1hi5 and 1hi5_W4d as a function of time. Increasing the density shortens the relaxation time, and after ∼ 2 Myr the clusters in simulations 1hi5 and 1hi5_W4d are similar. The cluster is able to hold onto its mass for slightly longer in 1hi5_W4d, but after the cluster expands it ultimately gets disrupted by the tidal field in the same way as 1hi5. As such, making clusters arbitrarily dense does not help the cluster migration scenario.



 §.§ Models with primordial binaries




We include a primordial binary population in three of our simulations 1hi10_b, 1hi10_v2b* and 1hi5_b. The inclusion of primordial binaries is interesting as the clockwise disk has three confirmed eclipsing binaries <cit.>, with a total binary population estimated to be greater than 30% <cit.>. Secondly, a popular formation scenario for the S-stars is from the tidal disruption of binaries by SgrA* via the Hills mechanism <cit.>, where one star is captured  and the other is ejected as a hyper-velocity star.


Fig <ref> shows the final projected distances of binaries in simulations 1hi10_b, 1hi10_v2b* and 1hi5_b. In these models 5% of the stars are initially in binary systems, the properties of which are described in section <ref>. A large number of binary systems survive, despite some being consumed during the formation of the VMS. The final distributions of binaries with main sequence primaries more massive than 8 M⊙ are very similar to the distribution of single stars more massive than 8 M⊙ in models with no primordial binaries.



In all three simulations, no binaries end up with semi-major axes less than 1pc. In 1hi10_v2b* one massive binary of total mass 68.9 M⊙ came within 0.1 pc of SgrA*. Fig. <ref> shows the orbit of this star, which came 0.09 pc from SgrA* at its third pericentre passage. However, the binary remained bound, as its tidal disruption radius by SgrA* was equal to ∼ 10 AU, ∼ 2000 times smaller than its distance. The binary coalesced at apocentre. Due to the scarcity of binaries that approach SgrA*, and the fact that many binaries would be observed beyond the disk, we conclude that if the binary breakup scenario is the origin of the S-stars, it is unlikely that the progenitors originated from nearby star clusters.



§ DISCUSSION





The formation and evolution of a VMS appears to have little effect on cluster inspiral as compared with the collisionless models of <cit.>, yet the suppression of IMBH formation strongly inhibits the radial migration of a sub population of massive stars towards the central parsec (F09). However, even in the case of IMBH formation, one would still observe a broad distribution of massive stars out to ∼ 10 pc, making the scenario unlikely even if IMBH formation were efficient.

All simulations form a large disk of massive stars from ∼ 1 - 10pc, contradicting observations. This implies that no cluster has been present in this region in the past ∼ 15 Myr, as a large population would still be visible with current telescopes. Two ∼ 10^5 M⊙ gas clouds, M-0.02-0.07 and M-0.13-0.08, are seen projected at ∼ 7 and ∼ 13 pc from SgrA* <cit.>, suggesting the presence of GMCs in this region is commonplace. However, the absence of young stars in this region suggests that perhaps GMC collapse is suppressed, only triggering from a tidal shock at close passage to SgrA* <cit.>. Verifying this hypothesis is beyond the scope of this work, and would require further study of GMC collapse in the close vicinity of a massive black hole.



§ CONCLUSIONS




We ran N-body simulations of young dense star clusters that form at distances of 5-15 pc from SgrA* and inspiral towards the Galactic Centre due to dynamical friction. Most models are dense enough that runaway collisions are inevitable, forming a very massive star in less than 1 Myr. However, careful treatment of the evolution of this very massive star shows that it is likely to lose most of its mass through its stellar wind and end its life as a ∼ 50-250 M⊙ black hole. As no intermediate mass black hole can form in this model, clusters dissolve a few pc from SgrA*, leaving a population of bright early type stars that would be observable for longer than the age of both the clockwise disk and S-star population, contradicting observations. It is therefore unlikely that a cluster has inhabited the central 10 pc in the last ∼ 15 Myr, as such the S-stars are unlikely to have formed via disrupted binaries originating from star clusters. Instead, the clockwise disk likely formed in-situ, perhaps from a gas cloud on a radial orbit incident on SgrA* <cit.>, and the S-star cluster is likely to be populated either by dynamical processes in the clockwise disk <cit.>, or through the binary breakup of scattered binaries <cit.>.



§ ACKNOWLEDGEMENTS

JAP would like to thank the University of Surrey for the computational resources required to perform the simulations used in this paper. 


mn2e,natbib





§ APPENDIX A: DYNAMICAL FRICTION COMPARISON WITH GADGET





In <cit.>, we only tested our dynamical friction formulation against N-body models of single component Dehnen profiles. In this paper we add an additional central massive black hole. In the Maxwellian approximation, valid for cuspy distributions, this comes into Chandrasekhar's formula via the black hole's contribution to the velocity dispersion of the stars. The addition of a black hole to the model is described in the appendix of <cit.>.

Here we briefly show that our dynamical friction formulation in the vicinity of a black hole is accurate by means of two N-body models of point mass clusters, of mass 10^5M⊙, orbiting the potential described in section <ref>. The N-body models are computer using the mpi-parallel tree-code GADGET2 <cit.>. The stellar background is comprised of 2^24 particles of mass 3.5M⊙, with a central black hole of 4.3×10^6M⊙. The softening of the cluster potential is ϵ = 0.0769pc, corresponding to rhm∼ 0.1pc. The same softening length is used for the background particles. The black hole is given a softening length, ϵ = 0.2pc to reduce numerical inaccuracies resulting from the large mass ratio of the black hole to background particles. In GADGET2 the force is exactly Newtonian at 2.8ϵ, so the semi-analytic and N-body models should agree to ∼0.56pc. The cluster is initially 5pc from SgrA*.

Fig. <ref> shows the orbital evolution of the circular and eccentric cases computed semi-analytically and with GADGET2. Our formalism agrees very well with the N-body models in both cases. In the eccentric case the circularisation appears to be slightly under-predicted. This is likely because we neglect the drag force from stars moving faster than the cluster <cit.>.





