Edinburgh Research Explorer
Central Dynamics of Multi-mass Rotating Star Clusters
Citation for published version:
Tiongco, M, Collier, A & Varri, AL 2021, 'Central Dynamics of Multi-mass Rotating Star Clusters', Monthly
Notices of the Royal Astronomical Society , vol. 506, no. 3, pp. 4488-4498.
https://doi.org/10.1093/mnras/stab1968
Digital Object Identifier (DOI):
10.1093/mnras/stab1968
Link:
Link to publication record in Edinburgh Research Explorer
Document Version:
Peer reviewed version
Published In:
Monthly Notices of the Royal Astronomical Society
General rights
Copyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)
and / or other copyright owners and it is a condition of accessing these publications that users recognise and
abide by the legal requirements associated with these rights.
Take down policy
The University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorer
content complies with UK legislation. If you believe that the public display of this file breaches copyright please
contact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately and
investigate your claim.
Download date: 10. Aug. 2022
Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 9 July 2021 (MN L
A
T
E
X style file v2.2)
Central Dynamics of Multi-mass Rotating Star Clusters
Maria Tiongco
1?
, Angela Collier
1
, and Anna Lisa Varri
2,3
†
1
JILA and Department of Astrophysical and Planetary Sciences, CU Boulder, Boulder, CO 80309, USA
2
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
3
School of Mathematics, University of Edinburgh, Kings Buildings, Edinburgh EH9 3JZ, UK
Accepted ?; Received ??; in original form ???
ABSTRACT
We investigate the evolutionary nexus between the morphology and internal kine-
matics of the central regions of collisional, rotating, multi-mass stellar systems, with
special attention to the spatial characterisation of the process of mass segregation.
We report results from idealized, purely N-body simulations that show multi-mass,
rotating, and spherical systems rapidly form an oblate, spheroidal massive core, unlike
single-mass rotating or multi-mass non-rotating configurations with otherwise identi-
cal initial properties, indicating that this evolution is a result of the interplay between
the presence of a mass spectrum and angular momentum. This feature appears to
be long-lasting, preserving itself for several relaxation times. The degree of flattening
experienced by the systems is directly proportional to the initial degree of internal
rotation. In addition, this morphological effect has a clear characterisation in terms of
orbital architecture, as it lowers the inclination of the orbits of massive stars. We offer
an idealised dynamical interpretation that could explain the mechanism underpinning
this effect and we highlight possible useful implications, from kinematic hysteresis to
spatial distribution of dark remnants in dense stellar systems.
Key words: methods: numerical – galaxies: star clusters: general – stars: kinematics
and dynamics
1 INTRODUCTION
The traditional picture of globular clusters as fully relaxed,
isotropic, non-rotating, spherical systems characterised by a
single old stellar population cannot last against the complex-
ities emerging from the new-generation data which are now
available and the theoretical ambitions that they stimulate.
In particular, after a few pioneering efforts (e.g., see Ander-
son & King 2003), there is now a convincing body of obser-
vational investigations mapping the internal kinematics of
several Galactic globular clusters (e.g., see Bianchini et al.
2013, Fabricius et al. 2014, Watkins et al. 2015, Ferraro et al.
2018, Kamann et al. 2018). Recent astrometric studies based
on Gaia Data Release 2 have further propelled the explo-
ration of the degree of anisotropy in the three-dimensional
velocity space (e.g., see Jindal et al. 2019), and confirmed
the growing evidence that the presence of internal rotation
in globular clusters is much more common than previously
assumed (e.g., see Bianchini et al. 2018, Sollima et al. 2019,
Vasiliev 2019).
This kinematic richness is now progressively being lever-
aged to attack a number of outstanding questions concerning
the internal dynamics of this class of stellar systems, from
?
E-mail: maria.tiongco@colorado.edu
the radial distribution of their angular momentum content
(e.g., see Bellini et al. 2017, Lanzoni et al. 2018a) to the
phase space properties of their present-day stellar popula-
tions (e.g., see Richer et al. 2013, Cordero et al. 2017, Cor-
doni et al. 2020).
Kinematic studies such as the those mentioned above
find an essential counterpart in detailed investigations of
the structural and morphological properties of star clusters,
which, after some early analyses (e.g., see Geyer et al. 1983,
White & Shawl 1987, Kontizas et al. 1989), unfortunately,
remain relatively scarce (e.g., see Chen & Chen 2010, Stet-
son et al. 2019). The exploration of the natural connection
between kinematics and morphology is indeed a crucial step
to fully understand the intrinsic phase space structure of
these stellar systems (e.g., see the informative studies con-
ducted by Davoust & Prugniel 1990, Han & Ryden 1994,
Ryden 1996, and, more recently van den Bergh 2008). Such a
joint approach can offer great insight into the physical origin
of their angular momentum content (e.g., see Frenk & Fall
1982, Fall & Frenk 1985), the importance of any tidal per-
turbation and, more generally, the different phases of their
dynamical evolution, as driven by the synergy of internal
and external processes.
Theoretical and numerical studies have indeed shown
that the total angular momentum content in collisional stel-
© 0000 RAS
arXiv:2107.03396v1 [astro-ph.GA] 7 Jul 2021
2 M. Tiongco, A. Collier, and A. L. Varri
lar systems is directly impacted by two-body relaxation
processes which determine its redistribution, transport and
eventual loss, especially in the case of tidally perturbed
systems (e.g., see Einsel & Spurzem 1999, Ernst et al.
2007, Hong et al. 2013, Tiongco et al. 2017). Therefore, the
strength and distribution of the angular momentum we mea-
sure in Galactic globular clusters today are remnant signa-
tures of the initial rotation content imprinted in such sys-
tems by their formation processes (e.g., see Lanzoni et al.
2018b for a comparison between the present-day rotation
curve of NGC 5904 and a long-term N-body simulation from
the survey by Tiongco et al. 2016).
More generally, even in the case of an isolated system,
the presence of non-vanishing total angular momentum may
lead to a more complex long-term dynamical evolution com-
pared to the one of a non-rotating system (e.g., the sugges-
tion by Hachisu 1979 of the existence of a “gravo-gyro catas-
trophe”, subsequently explored by several other authors).
Additional investigations have also highlighted interesting
effects of the interplay between bulk internal rotation and
a mass spectrum of stars in a stellar system. In particular,
Kim et al. 2004 showed that systems of that kind can pro-
duce an oblate core of fast rotating heavy masses, similar
to the one identified in the present study. Most recently,
Sz¨olgy´en et al. (2019) showed that in rotating star clusters,
the orbital inclinations of the heaviest stars decrease over
time, creating a mass segregation effect in the distribution
of orbital inclinations in addition to the well-known radial
(isotropic) mass segregation effect. This effect has also been
observed in simulations of stellar systems orbiting a mas-
sive black hole, such as nuclear star cluster simulations of
Sz¨olgy´en & Kocsis (2018), the eccentric nuclear disk simu-
lations of Foote et al. (2019).
So far, the explanation of why this effect occurs has
been attributed to resonant relaxation and resonant friction,
which was first introduced by Rauch & Tremaine (1996). In
a stellar system, coherent torques from stars on stable or-
bits enhance the rate of angular momentum relaxation. In a
stellar system dominated by a central mass, both the magni-
tude and direction of the angular momentum vectors change
in a random walk fashion, while in a stellar system without
a dominant central mass, only the direction of the angular
momentum vectors change stochastically. The latter process
is referred to as vector resonant relaxation and is considered
a more limited form of resonant relaxation. Meiron & Koc-
sis (2019) studied the effects of vector resonant relaxation
in globular cluster like systems. Rauch & Tremaine (1996)
also coined the term “resonant friction” that describes how
the orbital inclinations of massive objects in a stellar sys-
tem can be lowered by near-resonances. An application of
such concepts to the study of the statistical mechanics of
rotating systems with a central black hole has been recently
presented by Gruzinov et al. (2020).
In this work, we wish to concentrate on the investi-
gation of the evolutionary nexus between morphology and
kinematics of the central regions of collisional, rotating sys-
tems, with special attention to the process of mass segre-
gation. We perform and interpret a new series of N-body
simulations of rotating globular clusters with a spectrum of
stellar masses. In view of our specific interest in the central
dynamics and structural properties, we restrict our investi-
gation to a set of initial conditions characterised by spherical
Table 1. Summary of the properties of the N -body models pre-
sented in this study (in H´enon units). The 2nd column shows the
ratio of the initial rotational kinetic energy and total kinetic en-
ergy in the system, the 3rd column reports the initial value of the
spin parameter (Eq. 1) for each model, and the final two columns
denote the fraction of prograde particles for each mass bin (see
Section 2 for definitions). Rows 1-4 correspond to the primary set
of N-body models ordered by increasing spin; rows 5-9 refer to
the additional N-body experiments.
Model R
E
/R
KE
λ
Prograde Fraction (%)
H L
R000 0.001 0.001 50 50
R050 0.026 0.058 75 75
R075 0.056 0.084 87.5 87.5
R100 0.101 0.110 100 100
R
h
0.045 0.085 100 50
R
l
0.011 0.028 50 100
R100S50 0.131 0.111 100 100
R100S75 0.165 0.112 100 100
R
r
0.084 0.103 88.8 88.8
symmetry that distinguishes our study from previous simi-
lar investigations that used initially oblate, i.e., already flat-
tened, rotating models (commonly used models include the
ones from Lupton & Gunn 1987 and Varri & Bertin 2012).
Our parameter space explores different degrees of rotation
and velocity anisotropy, and we provide evidence that the
latter also plays a non-trivial role in the development of the
so-called ‘anisotropic mass segregation’. We explore in depth
such mass segregation process along different spatial direc-
tions within the cluster, and note that the central regions of
the cluster are flattened as an oblate spheroid that develops
and sustains itself for several relaxation times. We also of-
fer an idealised dynamical interpretation that could explain
the mechanism underpinning these effects and we present a
discussion of possible implication our results.
This article is structured as follows: numerical aspects
and initial conditions are described in Section 2 and our
results from numerical modeling are presented in Section 3.
Next, we discuss some analytical aspects of our results with
some additional numerical experiments to understand this
analysis further in Section 4. Conclusions drawn from this
work are reported in Section 5.
2 METHOD AND INITIAL CONDITIONS
In our survey of N-body simulations, the initial conditions
are first defined by a King (1966) spherical, isotropic, non-
rotating distribution function, with the concentration pa-
rameter chosen to be W
0
= 6. For the mass spectrum, we
adopted a power-law distribution with a slope of -2 and the
mass ratio of the heaviest star to the lightest star set to 100
(the resulting ratio of the heaviest star mass to the average
stellar mass is ≈ 21). Each N-body model has a total of
N = 65 536 particles. The initial conditions were generated
using the McLuster code (K¨upper et al. 2011).
Rotation is introduced into the system by randomly se-
lecting a fraction of particles and setting their tangential
velocities in the same direction, i.e., their tangential veloci-
ties are set as their absolute value. This implementation pre-
© 0000 RAS, MNRAS 000, 000–000
Central Dynamics of Rotating Star Clusters 3
serves the solution of the Boltzmann equation and also the
shape of the chosen equilibria; we refer to this change in the
velocities as the action of “Lynden-Bell’s demon” (Lynden-
Bell 1960, see also Rozier et al. 2019 for a recent applica-
tion). We note that the phase space invariance under veloc-
ity reversals is a direct result of the Jeans theorem (Jeans
1919). The models are named R for “rotating”, followed by
a number that denotes the fraction of particles that have
been selected for the application of the “demon”; we use the
terms “prograde” and “retrograde” to mean that v
φ
> 0 or
< 0, respectively, and the coordinate system is such that the
axis of rotation corresponds to the z direction). For example,
the R100 model has 100% of particles rotating in the same
direction. The models R000, R050, R075, and R100 com-
prise the primary models of our study, in order of increasing
rotation, and have identical initial positions of stars.
We measure the angular momentum content of each N-
body model via two distinct metrics. First, we calculate the
kinetic energy due to rotation and compare to the total ki-
netic energy in the model. The second parameter measured
to understand spin is the well-known “spin parameter” (e.g.,
Peebles 1969),
λ =
JE
1/2
GM
5/2
(1)
where J, E, and M are the total angular momentum, en-
ergy and mass respectively, and G is the gravitational con-
stant. Throughout the paper, we adopted the H´enon sys-
tem of units (H´enon 1971, Heggie & Mathieu 1986), where
G = M = 1 and E = −0.25. These two measures of an-
gular momentum content are listed in Table 1. The fiducial,
isotropic, R000 model has an initial λ of ∼ 0, and we increase
λ by increasing the fraction of prograde orbits.
We have also run some additional experiments to com-
plement our results, which we will describe in detail in Sec-
tion 4. To understand the roles of the different mass com-
ponents, we have created some additional N -body models
by dividing the given initial equilibrium into mass bins and
by rotating the bins individually (Experiment 1). The most
massive 1/3 of particles are denoted as “high mass” (H) and
the rest of the particles (the lower 2/3) are marked as “low
mass” particles (L).
In addition, we have also considered a subset of ini-
tial conditions (Experiment 2) characterised by some degree
of isotropic/radial mass segregation, as there is dynamical
and observational evidence suggesting the existence of pri-
mordial mass segregation in young star clusters (e.g., see
Bonnell & Davies 1998, de Grijs et al. 2002, McMillan et al.
2007). The prescription adopted to introduce such mass seg-
regation is based on the one featured in Baumgardt et al.
(2008), where, in short, by setting the segregation parameter
S closer to 1, the more likely a heavy particle is initialized
closer to the centre, with a segregation parameter of 1 being
fully segregated (with the heaviest particle at the shortest
radius from the centre, followed by the next heaviest particle
at the 2nd shortest radius, etc.).
Finally, we consider a more realistic rotation curve (Ex-
periment 3) defined by the rotation curve increasing from
zero from the centre of the cluster, peaking at approximately
the half-mass radius, then decreasing further out (see, e.g.,
Lanzoni et al. 2018a, Tiongco et al. 2017). This rotational
profile was generated by reversing the tangential velocities
of a different percentage of particles in each radial bin until
the desired profile is achieved as shown in Figure 1. Our sur-
vey of N-body simulations was performed using NBODY6
(Aarseth 2003) with GPU acceleration (Nitadori & Aarseth
2012). All N -body models are evolved in isolation and the
effects of stellar evolution are not included.
A reference time scale that we have adopted in all our
analyses is the initial half-mass relaxation time defined as
t
rh,i
=
0.138N
1/2
r
3/2
h
hmi
1/2
G
1/2
log(0.11N)
(2)
where hmi is the mean stellar mass, and r
h
is the 3D half-
mass radius, the radius enclosing half the mass of the cluster
(see e.g. Heggie & Hut 2003). In H´enon units, the t
rh,i
of all
of our models is ≈ 730 time units. For all of the models
featured in this study, the duration of the simulations is of
4000 time units.
We show relevant kinematic properties as a function
of radius for all of our models in Fig. 1, including ro-
tational velocity, velocity dispersion, and velocity disper-
sion anisotropy; the importance of the latter two proper-
ties are discussed in the next section. We acknowledge here
that while introducing rotation into the system via Lyden-
Bell’s demon is straightforward, the resulting initial rota-
tion curves do not resemble what is observed in globular
clusters. However, as these systems evolve over a short pe-
riod of time (i.e., over several dynamical times and much
shorter than a relaxation time), their angular momentum
distributions evolve into rotation curves which are compa-
rable to the ones observed in present-day star clusters. We
show the same kinematic properties in Fig. 1 after some
short evolution in Fig. A1. Overall, our main result shows
how increasing the amount of angular momentum affects the
spatial distribution of multi-mass stellar systems using ini-
tial conditions that do not change in physical structure when
increasing the amount of rotation, in contrast to initial con-
ditions realized from commonly used distribution function
based models with rotation, such as those from Lupton &
Gunn (1987) and Varri & Bertin (2012).
3 ANALYSIS OF PRIMARY MODELS
3.1 Central morphology
By the action of the “Lynden-Bell’s demon”, we have cre-
ated N -body models with initial conditions that differ in
their initial degree of rotation but maintain the original den-
sity and velocity distributions (i.e., their zeroth and second-
order velocity moments). Therefore, the difference between
the various N-body models appears in their first-order ve-
locity moment and corresponding velocity dispersion.
For clarity, we assume the following definition of the
velocity dispersion tensor (Binney & Tremaine 2008)
σ
2
i,j
= h(v
i
− hv
i
i)(v
j
− hv
j
i)i (3)
where i, j, k = r, θ, φ refer to conventional spherical coor-
dinates. While the initially non-rotating model (R000) is
isotropic, the process of introducing rotation lowers the com-
ponent σ
φ
, which, in turn, increases the associated degree of
velocity anisotropy. We depict the initial σ
φ
for each N -body
© 0000 RAS, MNRAS 000, 000–000
4 M. Tiongco, A. Collier, and A. L. Varri
0 0.1 0.2 0.3 0.4
v
φ
R000
R050
R075
R100
0 0.2 0.4 0.6 0.8
1
v
φ
σ
0
0
1 2
3
4
5 6
0 0.1 0.2 0.3 0.4 0.5
r r
h
σ
φ
R
r
R
h
R
l
R100S50
R100S75
0
1 2
3
4
5 6
0 0.2 0.4
r r
h
β
Figure 1. Top left: Initial rotational profile (mean velocity in the azimuthal/φ direction), as a function of cylindrical radius normalized
to the projected half-mass radius. Bottom left: Initial velocity dispersion σ
φ
profile. Top right: Initial rotational profile normalized to the
central velocity dispersion σ
0
=
q
1
3
(σ
2
r,0
+ σ
2
φ,0
+ σ
2
θ,0
). Bottom right: Initial velocity anisotropy parameter (β, see Eq. 4). In all panels,
solid lines denote the primary N-body models (Table 1, Row 1-4). Dashed lines represent the additional N-body models discussed in
Section 4 (Table 1, Row 5-9).
model in the top panel of Fig. 1. The other components of
the velocity dispersion tensor, σ
r
and σ
θ
, remain unchanged
under the action of the “Lynden-Bell’s demon” and are iden-
tical to the σ
φ
of the R000 model. The bottom panel of Fig.
1 illustrates the initial anisotropy parameter (β(r)) for each
N-body model; the definition adopted here is
β = 1 −
σ
2
jj
+ σ
2
kk
2σ
2
ii
. (4)
The initial differences in the amount of rotation and de-
gree of velocity anisotropy characterising our N -body mod-
els will subsequently determine, during their evolution, the
degree of flattening acquired by the central regions of the
systems. We stress that the flattening seen in our N-body
models should not be compared to the flattening observed
in present-day rapidly rotating systems, such as certain
classes of early-type galaxies, as the total amount of an-
gular momentum introduced via “Lynden-Bell ’s demon” is
much smaller than the one possessed by an oblate spheroid
whose morphology is completely shaped by rotation (e.g.,
see Chandrasekhar 1987, Binney 1978).
In all of our rotating cases, we observe the appearance
of a long-lasting oblate spheroidal structure in the core of
the N-body model, while, at larger radii, the system main-
tains a more spherical distribution. We illustrate this re-
sult in Fig. 2, showing the two-dimensional surface density
maps and a selection of demonstrative isodensity contours
of models R000 and R100 evaluated at 1.38t
rh,i
(i.e. a mo-
ment representative of a phase long after the anisotropic
mass segregation has peaked, discussed in the next section.)
At this time in the simulation, the R100 model shows a cen-
tral oblate shape in the (x, z) plane, while the R000 model
remains spherical throughout. The solid lines represent con-
tours that can be approximated by ellipses while the dashed
lines show a circular geometry; the contour containing 90%
of the total mass is denoted in red. We also note that rotat-
ing model has an oblate shape at this outer radii, while the
non-rotating model is more spherical.
To compare the morphology of all our N-body models,
we have measured the ellipticity of isodensity contours of the
© 0000 RAS, MNRAS 000, 000–000
