------------------------------------------------------------------------
sstars.tex ApJ, submitted
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-MailScanner-Information: Please contact the postmaster@aoc.nrao.edu for more information
X-MailScanner: Found to be clean
X-MailScanner-SpamCheck: not spam, SpamAssassin (not cached, score=0,
required 5, autolearn=disabled)
X-MailScanner-From: hagai.perets@weizmann.ac.il
X-Spam-Status: No
|%| arXiv:0903.2912||
|% |
%% LyX 1.6.1 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,twocolumn]{emulateapj}
%\documentclass[12pt,english,preprint]{emulateapj}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\setcounter{tocdepth}{3}
\usepackage{graphicx}
\usepackage{amssymb}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
%% Do not edit unless you really know what you are doing.
%\documentclass[english,twocolumn,dvips,usenames]{emulateapj}
\usepackage{array}
%\usepackage{dvipost}
%\usepackage{color}
\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
%\usepackage{dvipost}
%\usepackage{amsmath}
\usepackage{times}
\newcommand{\ns}{n_{\star}}
\newcommand{\Es}{E_{\star}}
\newcommand{\peryr}{\mathrm{yr}^{-1}}
\newcommand{\yr}{\mathrm{yr}}
\newcommand{\pc}{\mathrm{pc}}
\newcommand{\perpc}{\mathrm{pc}^{-1}}
\newcommand{\nbh}{n_{\bullet}}
\newcommand{\Md}{M_{\mathrm{disk}}}
\newcommand{\half}{{1\over 2}}
\newcommand{\scr}{s_{\mathrm{crit}}}
\shorttitle{Dynamical evolution of the S-stars}
\shortauthors{Perets et al.}
\makeatother
%\dvipostlayout
%\dvipost{osstart color push Red}
%\dvipost{osend color pop}
%\dvipost{cbstart color push Blue}
%\dvipost{cbend color pop}
\usepackage{babel}
\usepackage{babel}
\begin{document}
\global\long\global\long\def\Ms{M_{*}}
\global\long\global\long\def\LS{L_{*}}
\global\long\global\long\def\Rs{R_{*}}
\global\long\global\long\def\Mo{M_{0}}
\global\long\global\long\def\Lo{L_{0}}
\global\long\global\long\def\Ro{R_{0}}
\global\long\global\long\def\Mbh{M_{\bullet}}
\global\long\global\long\def\SgrA{Sgr\, A^{\star}}
\global\long\def\np{n_{p}}
\global\long\def\Np{N_{P}}
\global\long\def\Mp{M_{p}}
\global\long\def\Ns{N_{\star}}
\global\long\def\rMP{r_{\mathrm{MP}}}
%\centerline{\rule{\columnwidth}{1pt}}
\title{Dynamical evolution of the young stars in the Galactic center:
N-body simulations of the S-stars }
\author{Hagai B. Perets\altaffilmark{1}, Alessia Gualandris\altaffilmark{2},
Gabor Kupi\altaffilmark{1}, David Merritt\altaffilmark{2} and
Tal Alexander\altaffilmark{1} }
76100, Israel} \altaffiltext{2}{Department of Physics and Center
for Computational Relativity and Gravitation, Rochester Institute
of Technology, Rochester, NY 14623}
\begin{abstract}
We use $N$-body simulations to study the evolution of the orbital
eccentricities of stars deposited near ($\lesssim0.05$ pc) the Milky
Way massive black hole (MBH), starting from initial conditions motivated
by two competing models for their origin: formation in a disk followed
by inward migration; and exchange interactions involving a binary
star. The first model predicts modest eccentricities, lower than those
observed in the S-star cluster, while the second model predicts higher
eccentricities than observed. The $N$-body simulations include a
dense cluster of $10M_{\odot}$ stellar black holes (SBHs), expected
to accumulate near the MBH by mass segregation. Perturbations from
the SBHs tend to randomize the stellar orbits, partially erasing the
dynamical signatures of their origin. The eccentricities of the initially
highly eccentric stars evolve, in $20$ Myr (the S-star lifespan),
to a distribution that is consistent at the $\sim\!95\%$ level with
the observed eccentricity distribution. In contrast, the eccentricities
of the initially more circular orbits fail to evolve to the observed
values in 20 Myr, arguing against the disk migration scenario. We
find that $20\%$--$30\%$ of the S-stars are tidally disrupted by
the MBH over their lifetimes, and that the S-stars are not likely
to be ejected as hypervelocity stars outside the central 0.05 pc by
close encounters with stellar black holes.
\end{abstract}
\keywords{black hole physics --- galaxies: nuclei --- stars: kinematics }
\section{Introduction}
In recent years, high resolution observations have revealed the existence
of many young OB stars in the Galactic center (GC). Accurate measurement
of the orbital parameters of these stars gives strong evidence for
the existence of a massive black hole (MBH) which dominates the dynamics
in the GC \citep{ghe+98,eis+05,gil+08}. Most of the young stars are
observed in the central 0.5 pc around the MBH. The young star population
in the inner 0.05 pc (the {}``S-stars'') consists exclusively of
B-stars, in an apparently isotropic distribution around the MBH, with
relatively high eccentricities ($0.3\lesssim e\lesssim0.95$;
\citealt{gil+08}).
The young stars outside this region comprise O-stars in one or two
disks, and present markedly different orbital properties
\citep{lev+03,bar+08,lu+09}.
Since regular star formation in the region near the MBH is inhibited
by tidal forces, many suggestions have been made regarding the origin
of the S-stars. Many of these are probably ruled out by observations
and/or by theoretical arguments \citep[see][for a
review]{ale05,pau+06}. The various scenarios for the origin of the
S-stars predict very
different distributions for their orbits, which in principle could
be constrained by observations. However, it is not clear to what extent
relaxation processes can produce changes in the distribution of orbital
parameters after the stars have been deposited near their current
locations. Here we try to resolve this question. We use $N$-body
simulations to follow the evolution of stellar orbits around the GC
MBH for $20$ Myr, starting from various initial conditions that were
motivated by different models for the origin of the S-stars. However,
we do not discuss here the possibility that an intermediate mass black
hole was involved in the production and/or evolution of the S-stars,
which is discussed in details elsewhere \citep[see ][and references
therein]{mer+09}.
We find our N-body results to be consistent with our analytical predictions,
and compare them with current observations. We then discuss the implications
for the validity of the models for the production of the S-stars.
In addition we study the possible ejection of the S-stars outside
the inner $0.05$ pc and the contribution of such ejected stars to
the population of hypervelocity stars (\citet{bro+08})
and to the isotropic population of B-stars observed at distances of
up to $0.5$ pc from the MBH (i.e. at distances similar to those of
the young O/WR stellar disks, but outside of these disks at high
inclinations).
In \S2 we summarize the different models for the origin of the S-stars
and the predictions that they make for the stellar orbits at the time
when the stars are first deposited near the MBH. \S3 describes the
$N$-body simulations we carried out to follow the S-star orbital
evolution starting from these initial conditions. \S4 and \S5 present
the results of these simulations and discusses the implications for
the origin of the stellar populations of B-stars in the GC and for
the population of hypervelocity stars observed in the Galactic halo.
\S6 sums up.
\section{Models for the S-stars origin}
Many solutions have been suggested for the origin of the S-stars,
but many of these have been effectively excluded (see \citealp{ale05,pau+06}
for a review). Here we focus on two basic models which differ substantially
in their predictions for the initial orbital distribution of the S-stars
and/or the time passed since their arrival/formation at their current
location. These are (1) formation of the S-stars in a stellar disk
close to the MBH, followed by transport through a planetary-migration-like
scenario to their current positions \citep{lev07}; (2) formation
of the S-stars in binaries far from the MBH followed by scattering
onto the MBH by massive perturbers (e.g. giant molecular clouds) and
tidal disruption of the binaries \citep{per+07,per+08}, leaving a
captured star in a tight orbit around the MBH. Binary disruption scenarios
similar to (2) have been proposed, in which the S-stars formed in
a stellar disk (either the currently observed 6 Myr old disk, or an
older, currently not observed disk) and later changed their orbits
due to coherent torques through an instability of eccentric disks
\citep{mad+08}; or through the Kozai mechanism resulting from the
presence of two disks \citep{loc+08b}. The former of these two alternative
models may not be favored since the observed disk shows only moderate
eccentricities ($\sim0.34$; \citealp{bar+08}) and has too short
a lifetime for these processes to be effective \citep{mad+08}. The
latter alternative is also unlikely since the Kozai mechanism is quenched
in the presence of a massive enough cusp of stars such as exists in
the GC \citep{cha+08,mad+08}). Furthermore, the disk-scattering scenario
could be at work only in the last $\sim6$ Myr if originating in the
currently observed stellar disk(s), while the massive perturber scenario
could be at work at the last few tens of Myr (up to the life span
of the S-stars). In any case, these alternative binary-disruption
scenarios imply very similar initial distributions for the captured
S-stars.
In the following, we briefly discuss the initial distribution of the
eccentricities and inclinations of the S-stars expected from the different
scenarios for their production. The models are summarized in Table
\ref{t:models}.
%
\begin{table*}
\caption{\label{t:models}Models for the S-stars origin }
\centering{}\begin{tabular}{llccccl}
\hline
\# & Origin & Initial & Time & Model & Survival &
Refs.$^{c}$\tabularnewline
& & Eccentricity & (Myr) & Probability$^{a}$ & Fraction$^{b}$ &
\tabularnewline
\hline
1 & Capture following Binary disruption due & High & $6$ & $0.26$ &
$0.8$ & 1\tabularnewline
& massive perturbers or to disk instability & ($0.94\le e\le0.99$) &
& & & \tabularnewline
& in an old non-observed disk & & & & & \tabularnewline
2 & Capture following Binary disruption due & High & $20$ & $0.93$
& $0.7$ & 2\tabularnewline
& to disk instability in currently observed disk & ($0.94\le
e\le0.99$) & & & & \tabularnewline
3 & Disk formation + planetary like & Low & $6$ & $8\times10^{-3}$
& $0.9$ & 3\tabularnewline
& migration (currently observed disk) & ($e\le0.5$) & & & &
\tabularnewline
4 & Disk formation + planetary like & Low & $20$ & $0.06$ & $0.8$
& 3\tabularnewline
& migration (possible old disk) & ($e\le0.5$) & & & & \tabularnewline
\hline
\multicolumn{7}{l}{$^{a}${\footnotesize Probability for the samples of
the observed
and simulated S-stars}}\tabularnewline
\multicolumn{7}{l}{{\footnotesize to be randomly chosen from the same
distribution (see
text).}}\tabularnewline
\multicolumn{7}{l}{$^{b}${\footnotesize Fraction of S-stars not
disrupted by the MBH
during the simulation (see text).}}\tabularnewline
\multicolumn{7}{l}{$^{c}$ (1) \citet{mad+08} (2) \citet{per+07}
(3)\citet{lev07}}\tabularnewline
\hline
\end{tabular}
\end{table*}
\subsection{Binary disruption by a massive black hole}
\label{sec:Binary-disruption}
A close pass of a binary star near a massive black hole results in
an exchange interaction, in which one star is ejected at high velocity,
while its companion is captured by the MBH and left on a bound orbit.
Such interaction occurs because of the tidal forces exerted by the
MBH on the binary components. Typically, a binary (with mass, $M_{bin}$,
and semi-major axis, $a_{bin}$), is disrupted when it crosses the
tidal radius of the MBH (of mass $\Mbh$), given by
$r_{t}=a_{bin}(\Mbh/M_{bin})^{1/3}$.
One of the binary components is captured by the MBH \citep{gou+03}
on a wide and eccentric orbit while the companion is ejected with
high velocity \citep{hil88}.
The capture probability and the semi-major axis distribution of the
captured stars were estimated by means of numerical simulations, showing
that most binaries approaching the MBH within the tidal radius
$r_{t}(a_{bin})$
are disrupted \citep{hil91,hil92,gua+05,bro+06c}. The harmonic mean
semi-major axis for 3-body exchanges with equal mass binaries was
found to be \citep{hil91} \begin{equation}
\left\langle a_{cap}\right\rangle
\simeq0.56\left(\frac{\Mbh}{M_{\mathrm{bin}}}\right)^{2/3}a_{bin}\simeq\,0.56\left(\frac{\Mbh}{M_{\mathrm{bin}}}\right)^{1/3}r_{t},\label{e:afinal}\end{equation}
where $a_{bin}$ is the semi-major axis of the infalling binary and
$a_{cap}$ that of the captured star (the MBH-star {}``binary'').
Most values of $a_{cap}$ fall within a factor $2$ of the mean. This
relation maps the semi-major axis distribution of the infalling binaries
to that of the captured stars: the harder the binaries, the more tightly
bound the captured stars. The periapse of the captured star is at
$r_{t}$, and therefore its eccentricity is very high \citep{hil91,mil+05},
\begin{equation}
e\!=\!1-r_{t}/a_{cap}\!\simeq1-1.8(M_{\mathrm{bin}}/\Mbh)^{1/3}\!\simeq\!0.94-0.99\label{eq:capture-eccentricity}\end{equation}
for values typical of B-type main sequence binaries and the MBH in
the GC ($M_{bin}=6-30\,\Mo$; $\Mbh=3.6\times10^{6}\,\Mo$;
$a_{cap}=0.5-2\times\left\langle a_{cap}\right\rangle $).
Therefore, in order to study the evolution of S-stars from the binary
disruption scenarios we assume that the initial eccentricities of
S-stars are in the range $0.94-0.99$ (where most are close to the
mean value of $0.98$).
In principle the binary disruption scenario has specific predictions
for the semi-major axis distribution of the captured stars, which
could also be used for constraining the model. However such distribution
is highly sensitive to differences in the (unknown) binary distribution
in the GC region. The prediction of high eccentricities for the captured
S-stars, instead, is robust and has only a weak dependence on the
mass of the binary.
The inclinations of the captured S-stars in the massive perturbers
scenario \citep{per+07} are likely to be distributed isotropically
since the stars originate in an isotropic cusp. In the disk instability
scenario \citep{mad+08}, however, the progenitors of the captured
stars originate from a disk and therefore their inclinations should
resemble that of the stellar disk, i.e. have a very limited range,
of the order of the stellar disk angular thickness.
\subsection{Planetary like migration from the young stellar disk(s)}
\citet{lev07} suggested that the S-stars could have formed in the
currently observed stellar disk in the GC \citep{bar+08,lu+09}, and
then migrated inward in a way similar to planetary migration. The
migration time scale expected from such a scenario could be as short
as $10^{5}$ yr (for type I migration), which could be comparable
to (although possibly larger than; \citet{nay+07}) the lifetime of
the gaseous disk. Recent analytic work (\citet{ogi+03,gol+03b}; and
references therein) has shown that eccentricity is likely to be damped
during migration, unless eccentricity excitation occurs, which requires
the opening of a clean gap in the disk. In the latter case the migration
timescale might be larger \citep{lev07}, possibly inconsistent with
the lifetime of the gaseous disk \citep{nay+07}. It is therefore
more likely that the eccentricities of the stars are damped during
the migration. Even eccentricity excitation, if such took place, is
unlikely to excite very high eccentricities. The mean eccentricity
of the observed stars in the stellar disk is $0.34\pm0.06$. Therefore,
in order to study the evolution of the S-stars following their formation
in a stellar disk, we assume them to have low eccentricities, or,
being conservative, moderately high eccentricities ($e_{max}=0.5$;
where we use a thermal distribution of eccentricities, cut off at
$e_{max}$). These simulations include the less likely (\citealt{lev07})
possibility that the stellar disk extended inward to the current region
of the S-stars, in which case the S-stars were formed in-situ and
did not migrate.
\section{The N-body simulations}
To test these competing models, we carried out $N$-body simulations
of the inner Milky Way bulge using models containing a realistic number
of stars. All integrations were carried out on the 32-node GRAPE cluster
\texttt{gravitySimulator}
%\footnote{http://wiki.cs.rit.edu/bin/view/GRAPEcluster}
at the Rochester Institute of Technology which adopts a parallel
setup of GRAPE accelerator boards to efficiently compute gravitational
forces. The direct-summation code $\phi$GRAPE was used \citep{har+07}.
The simulations used a softening radius of $\sim4\, R_{\odot}$ comparable
to the radius of the S-stars $r_{\star}$, so as to be able to follow
even the closest encounters between stars.
Our initial conditions were based on a collisionally evolved model
of a cusp of stars and stellar remnants around the GC MBH
\citep[hereafter HA06; see also \citealt{fre+06}]{hop+06b}.
HA06 evolved the multi-mass isotropic Fokker-Planck equation representing
the stellar distribution in the region extending to $\sim1$ pc from
the MBH; in their models the contribution to the gravitational potential
from the distributed mass is ignored. HA06 fixed the relative numbers
of objects in each of four mass bins by assuming a mass function consistent
with continuous star formation. The $10M_{\odot}$ stellar-mass black
holes (SBHs) were found to follow a steep, $n(r)\sim r^{-2}$ density
profile near the MBH while the lower-mass populations (main-sequence
stars, white dwarfs, neutron stars) had $n\sim r^{-\alpha}$,
$1.4\lesssim\alpha\lesssim1.5$.
The SBHs were found to dominate the mass density inside $\sim0.01$
pc.
Based on this model, we constructed an $N$-body realization containing
a total of 1200 objects within 0.3 pc of the MBH: 200 {}``stars,''
with masses of $3\, M_{\odot}$, and 1000 {}``black holes'' with
masses $10\, M_{\odot}$, around a MBH of $3\times10^{6}\, M_{\odot}$.
We set $\alpha=2$ for the SBHs and $\alpha=1.5$ for the lower-mass
stars, and each density component was tapered smoothly to zero beyond
$0.1$ pc when computing the corresponding $f(E)$. Since the S-stars
may have masses as high as $\sim10\, M_{\odot}$, the higher mass
stars in our simulations could also be treated as S-stars. We did
not see any major differences in the evolution of the more massive
and the less massive stars, and we discuss the evolution of both together.
The number of SBHs contained within a radius $r$ in our $N$-body
models was \begin{equation}
N(0.94-0.99$; cf. section \ref{sec:Binary-disruption}).
Under these assumptions, the stars evolved for $6$ Myr (if formed
in the stellar disk) or longer (if the S-stars formed outside the
central pc, i.e. not in the young stellar disk). We evolved the models
for up to $20$ Myr, which is comparable to the lifetime of the most
massive observed S-star (S2). In the second scenario we assumed that
the S-stars formed in a gaseous disk and then migrated inwards (or
formed in situ in a disk extending close to the MBH). For this case
we assumed the S-stars to have low eccentricities ($<0.5$), as typical
of disk formation models, and to evolve for $6\,$Myr (the lifetime
of the observed stellar disk). In order to check both scenarios we
selected the stars with initially high eccentricity orbits ($0.94\leq
e\leq0.99$)
and low eccentricity orbits ($e\leq0.5$) and followed their evolution
for a time appropriate to their presumed origin.
In addition to these evolutionary scenarios we also studied the possibility
of ejection of S-stars as hypervelocity stars, following a close encounter
with a SBH in the vicinity of the MBH, as suggested by \citet{ole+07}.
Our high resolution simulations can accurately follow close encounters
between stars, and therefore track any resulting high velocity ejections
of stars. Motivated by recent observations of B-type main sequence
stars outside the inner $0.05$ pc of the Galaxy, which show evidence
of random orientations similar to those of the S-stars, we also looked
for stars scattered to larger distances by gravitational encounters.
In other words, we considered whether captured S-stars could dynamically
evolve to become the extended B-type stars population observed outside
the central $0.05$ pc.
Although some binaries could exist in the close regions near the MBH
\citep{per08b,per09}, these are likely to be rapidly disrupted and
not play an important role in the dynamical evolution of the S-stars
\citep{per09}. We therefore do not include any binaries in our initial
conditions for the S-stars and SBHs evolution.
\section{Results}
\subsection{Simulations vs. theory: resonant relaxation}
\label{sub:RR}
We applied the correlation curve method \citet{eil+08} to our simulations
to identify the relaxation process responsible for the dynamical evolution
of the stars (Fig. \ref{f:RR}). The method is able to detect and
measure relaxation in nearly Keplerian $N$-body systems. In the isotropic
system considered here, the angular momentum of the stars, $J$, evolves
both due to the slow stochastic two-body relaxation (e.g. \citealt{bin+87})
and to the rapid resonant relaxation (RR)
\citep{rau+96,rau+98,hop+06a,gur+07,eil+08}.
Two-body relaxation changes $J$ in a random walk fashion, $\left|\Delta
J\right|/J_{c}\!=\sqrt{\tau/\tau_{NR}}$
over the long two-body relaxation time-scale $\tau_{NR}\!\sim\!
Q^{2}/(N\log Q)$,
where $J_{c}$ is the maximal (circular orbit) angular momentum for
a given energy, $Q\!=\!\Mbh/\Ms$, $N$ is the number of enclosed
stars on the distance scale of interest, and time $\tau\!=\! t/P$
is measured in terms of the orbital period on that scale. Resonant
relaxation occurs when the symmetries of the potential act to constrain
the stellar orbits (e.g. closed ellipses in a Kepler potential, or
planar rosettes in a spherical one). As long as the symmetry is
approximately
maintained on the coherence timescale $t_{w}$, the stars experience
coherent torques, and $\Delta J/J_{c}\!\sim\!(\sqrt{N}/Q)\tau$. In
a nearly Keplerian potential, as is the case in the inner parsec of
the GC, RR can change both the magnitude of $J$ ({}``scalar RR'')
and its direction ({}``vector RR''). In the Newtonian context, the
coherence of scalar RR is limited by the precession of the apoapse
due to the enclosed stellar mass, on a time-scale $\tau_{M}\!\sim\! Q/N$.
On timescales $\tau\!>\!\tau_{M}$, the coherent change $\Delta J(\tau_{M})$
becomes the mean free path in $J$-space for a rapid random walk,
$\left|\Delta J\right|/J_{c}\!=\!\sqrt{\tau/\tau_{sRR}}$, where
$\tau_{sRR}\!\sim\! Q\!\ll\!\tau_{NR}$.
The coherence of vector RR is self-limited by the change in the orbital
orientation due to RR, and is even faster,
$\left|\Delta\mathbf{J}\right|/J_{c}\!=\!\sqrt{\tau/\tau_{vRR}}$,
where $\tau_{vRR}\!\sim\! Q/\sqrt{N}$.
Figure (\ref{f:RR}) shows the rms change in the scalar and vector
angular momentum of the stars in the simulation, as a function of
the time lag $\tau$ (correlation curves), up to the maximal time
lag for which the simulation can still be analyzed with high statistical
confidence, $\tau_{\max}\sim10^{4}$. In our simulation%
\footnote{For a spectrum of masses, $\tau_{M}\!\sim\! Q/N\left\langle
\Ms\right\rangle $
and $\tau_{sRR}\!\sim\!\Mbh/M_{\mathrm{eff}}$, where $\left\langle
\Ms\right\rangle \!=\!8.8\, M_{\odot}$
and $M_{\mathrm{eff}}\!=\!\left\langle \Ms^{2}\right\rangle
/\left\langle \Ms\right\rangle $.
$M_{\mathrm{eff}}\!\simeq\!9.7\, M_{\odot}$ for our simulation.%
}, the scalar RR coherence time expected from theory is
$\tau_{M}\!\sim\!283$
and the scalar RR timescale is $\tau_{sRR}\!\sim\!3.1\times10^{5}$.
The behavior of the scalar correlation curve clearly indicates that
RR dominates relaxation in our GC model. The turn from the coherent
phase to the random-walk phase of RR is seen at $\tau\!\sim\!300$,
close to the predicted value of $\tau_{M}$. The random-walk growth
continues up to full randomization and saturation at $\sim\!\tau_{\phi}$,
as expected \citep{eil+08}. The scalar RR timescale can be estimated
directly from the simulated data by $\tau_{RR}=\tau/(\left|\Delta
J\right|/J_{c})^{2}$;
substituting $\left|\Delta J\right|/J_{c}\!\sim\!0.15$ at
$\tau\!\sim\!\tau_{\phi}/2\!\sim\!5000$
yields $\tau_{sRR}\!\sim\!2.2\times10^{5}$, in good agreement with
the predicted value.
%
\begin{figure}[t]
\noindent \begin{centering}
\includegraphics[width=1\columnwidth]{rr_autofit_10Mo_SSBH1c}
\par\end{centering}
\caption{\label{f:RR}The correlation curves $\left|\Delta J\right|/J_{c}$
and $\left|\Delta\mathbf{J}\right|/J_{c}$ as a function of the time
lag $\tau$ for all the stars in the simulation. The simulation data
(points) are compared to the best fit theoretical curves (lines).
Four regimes are seen \citep[see
text and also]{eil+08}. At very short timescales, non-resonant
two-body relaxation ($\propto\!\sqrt{\tau}$)
dominates over RR ($\propto\!\tau$); At $\tau\!<\!\tau_{M}$ the
curves display a coherent, linear rise; At
$\tau_{M}\!<\!\tau\!<\!\tau_{\phi}$,
the curves display a random-walk growth and at $\tau\!>\!\tau_{\phi}$
they begin to saturate. }
\end{figure}
Furthermore, since $J/J_{c}\!=\!\sqrt{1-e^{2}}$, the RR-driven eccentricity
evolution relates the final eccentricity, $e_{f}$, after time lag
$\tau$ to the initial eccentricity, $e_{i}$. We therefore expect
the eccentricities at some given time lag to have a typical spread,
$e^{-}(\tau)\lesssim e_{f}(\tau)\lesssim e^{+}(\tau)$, where
\begin{equation}
e^{\pm}(\tau)=\left[1-\left(\sqrt{1-e_{i}^{2}}\mp\sqrt{\frac{\tau}{\tau_{RR}}}\right)^{2}\right]^{1/2}.\end{equation}
The magnitude of the predicted change in eccentricity agrees well
with that observed in the simulations. For example, an S-star initially
captured by a tidal exchange event on a $P=500$ yr ($a\!\sim\!0.04$
pc), $e_{i}=0.97$ orbit can evolve by RR to an $e_{f}\!=\!0.80$
orbit in 20 Myr. The short vector RR timescale $\tau_{vRR}\!\sim\!10^{4}$
($t_{vRR}\!=\!5\times10^{6}$ yr at 0.04 pc) implies full randomization
of the orbital planes after 6 Myr, throughout the S-cluster volume,
as is observed in the simulation. We therefore conclude that RR is
the dominant mechanism responsible for the dynamical evolution of
the S-stars and other stars close to the MBH in the GC.
\subsection{S-star eccentricities and inclinations}
In Figure~\ref{f:eccentricities-cum} we show the final cumulative
eccentricity distribution of the S-stars for the different origin
models (Table \ref{t:models}). These are compared to the the orbits
of the observed S-stars (taken from \citealp{gil+08}). The probabilities
for the samples of the observed and simulated S-stars to be randomly
chosen from the same distribution (calculated using a two-sample $\chi^{2}$
test) are given in Table \ref{t:models}. We find that the binary
disruption model taking place at least $20$ Myr ago is much favored
over all other models tested here. We find $93\%$ chance for the
observed S-stars and the simulated stars in such model to originate
from the same distribution (with the shorter timescale of $6$ Myr
consistent only at the $26\%$ level) . In contrast, the disk migration
scenarios seem to be excluded (for the given assumptions), since they
have major difficulties in explaining the large fraction of eccentric
orbits observed for the S-stars in the GC.
We find that the inclination distribution of the captured stars, initialized
with the same inclination, is rapidly isotropized, to resemble a random
distribution of inclinations (consistent with random at the $\sim65\%$
and $\sim20\%$ level, after evolution of $20$ and $6$ Myr, respectively).
This is expected from the RR process, as discussed in the previous
section (see also fig. \ref{f:RR}). We conclude that the observed
isotropic distribution of the S-stars angular momentum direction is
consistent with all the S-stars production models studied here, and
can not be used to discriminate between them, although it constrains
the lifetime of the S-stars system to be at least $\sim4$ Myr at
the $95\%$ level, assuming all S-stars were initially put on the
same plane.
%
\begin{figure}[!tbp]
\includegraphics[scale=0.3]{eccs_cum_dist}
\caption{{\footnotesize Cumulative distribution of observed and
simulated S-stars
eccentricities, for various models (see legend).}}
\label{f:eccentricities-cum}
\end{figure}
\subsection{Survival of the S-stars: tidal disruption, ejection and
hypervelocity
stars}
As already discussed above, the S-stars can change their orbits due
to their dynamical evolution. A star could therefore be scattered
very close to the MBH and be disrupted by it, if its pericenter distance
from the MBH becomes smaller than the tidal radius of the star
$r_{t}=r_{\star}(M_{\bullet}/m_{\star})^{1/3}$.
Many of the S-stars could therefore not survive for long close to
the MBH. We followed the orbits of stars in our simulations and calculated
the fraction of stars that have been disrupted. For the tidal disruption
calculations all stars were assumed to have the typical main sequence
radius according to their mass. We consider a star as being disrupted
if its pericenter became smaller than twice the tidal radius during
the simulation (i.e., when it is strongly affected by the MBH tidal
forces or even totally disrupted in one pericenter passage). We find
that most of the S-stars survived to current times in all the models
(see survival fractions in \ref{t:models}). The S-stars population
in the GC therefore gives a good representation of all the S-stars
formed/captured in this region. The production rate of the S-stars
required to explain current observations is therefore only slightly
higher ($1.2-1.3$ times higher) than that deduced from current number
of S-stars observed.
In principle the S-stars could be ejected by strong encounters to
orbits with larger semi-major axes, putting them outside the 0.05
pc region near the MBH \citep{mir+00}, or even ejecting them as unbound
hypervelocity stars \citep{ole+07}. The softening radius used in
our simulation was $r_{soft}=4R_{\odot}$, comparable to the radius
of observed S-stars, allowing us to follow even very close encounters.
Nevertheless, we find that only $\sim10\%$ of the stars were ejected
outside of the central $0.05$ pc, and even those had maximal semi-major
axis in the end of the simulation not extending beyond 0.1 pc. We
therefore conclude that such ejected S-stars can not explain recent
observations of many B-type stars outside the central $0.05$ pc (H
Bartko, private communication). Moreover, none of the $3\, M_{\odot}$
stars in our simulations have been ejected as a hypervelocity star,
suggesting that ejection of hypervelocity stars through encounters
with SBHs is not an efficient mechanism (see also \citealp{per09}
for the related constraints on this mechanism). We note that since
our simulations can be rescaled, we can probe much higher stellar
densities and smaller softening radius (see next section). When such
rescaling is used, in which all the 1200 stars are distributed between
$3\times10^{-4}$ pc to $0.1$ pc (rescaling the radii by half), we
find $5$ stars ejected beyond $0.1$ pc (to have final semi-major
axis of $0.1,\,0.16,\,0.16,\,0.21$ and $0.37$ pc), but none ejected
as hypervelocity stars or even as slow unbound stars.
\section{Dependence of the results on the assumed density of SBHs}
As discussed above, the density of the SBHs that are responsible for
the evolution in our $N$-body models is not well determined. Scaling
the $N$-body results to different assumed values of $N$ is complicated
by the fact that scalar resonant relaxation has two regimes, coherent
and random walk. In our models, the transition occurs at \begin{equation}
t_{M}\simeq
QP/N\sim3\times10^{4}\left(\frac{N_{\mathrm{SBH}}}{10^{3}}\right)^{-1}\left(\frac{P}{100\,\mathrm{{\rm
\textrm{yr}}}}\right)\,{\rm \textrm{yr}\,},\end{equation}
where the $N_{\mathrm{SBH}}^{-1}$ scaling is for scattering by SBHs
of a given mass.
In the coherent regime, $\Delta t\lesssim t_{M}$, orbital angular
momenta grow as
\begin{equation}
\frac{\Delta J}{J_{c}}\sim10^{-4}\left(\frac{N_{{\rm
SBH}}}{10^{3}}\right)^{1/2}\frac{\Delta t}{P}\,.\end{equation}
In the diffusive regime, $\Delta t\gtrsim t_{M}$, \begin{equation}
\left|\frac{\Delta
J}{J_{c}}\right|\sim\sqrt{\frac{\tau}{\tau_{sRR}}}\approx2\times10^{-3}\sqrt{\frac{\Delta
t}{P}}\,,\end{equation}
independent of $N_{\mathrm{SBH}}$.
We are interested in the orbital evolution of the S-stars over timescales
of $\Delta t\sim O(10^{7})$ yr. In our simulations and those with
larger $N$, $\Delta t\gg t_{M}$. In this large-$N$ regime, changes
in eccentricity are dominated by the diffusive relation and are therefore
expected to be nearly independent of $N$ for time scales of interest,
at least up to $N$-values of $\sim10^{5}$ where the distributed
mass begins to approach the mass of the MBH and resonant relaxation
is no longer effective. Only for $N_{\mathrm{SBH}}\lesssim10$ does
$t_{M}$ approach $10^{7}$ yr and our results start depending significantly
on $N_{\mathrm{SBH}}$. However, such a small number of SBHs in the
volume of interest is highly unlikely, and therefore our results are
robust to the details of the SBH cusp model.
The $N$-body simulations can also trivially be rescaled by \begin{equation}
r\rightarrow Ar,\ \ \ \ t\rightarrow A^{3/2}t\end{equation}
at fixed mass. This corresponds to placing the same number of SBHs
into a smaller (larger) region and integrating for a shorter (longer)
time. For instance, if we rescale our simulations to three times smaller
distances (in which case the stars are distributed between $3\times10^{-4}$
and $0.05$ pc), the integration time becomes $\sim5$ Myr.
\section{Summary}
The approximately thermal eccentricity distribution of the S-stars
near the massive black hole (MBH) in the galactic center (GC),
$N(