Simulating the Spatial Distribution and Kinematics of Globular Clusters within Galaxy Clusters in Illustris
Felipe Ramos-Almendares Laura V. Sales Mario G. Abadi Jessica E. Doppel Hernan Muriel Eric W. Peng


firstpage–lastpage
Differentially private sub-Gaussian location estimators
    Marco Avella-Medina Columbia University, Department of Statistics, marco.avella@columbia.edu Victor-Emmanuel Brunel ENSAE ParisTech, Department of Statistics, victor.emmanuel.brunel@ensae.fr
    December 30, 2023
==================================================================================================================================================================================================





  We study the assembly of globular clusters (GCs) in
  9 galaxy clusters using the cosmological simulation
  Illustris.   GCs are tagged to individual galaxies at their infall time. The tidal removal of GCs from their galaxies and the distribution of the GCs within the cluster is later followed self-consistently by the simulation. 
  
  
  
  The method relies on the simple assumption of a single power-law 
  relation between halo mass (M_ vir) and mass in GCs (M_ GC)
  as found in observations. We find that the GCs
  specific frequency S_N as a function of V-band magnitude naturally
  reproduces the observed “U”-shape due to the 
  combination of the power law M_GC-M_ vir relation and the 
  
  non-linear stellar mass (M_*) - halo mass relation from the simulation. 
  Additional scatter in the S_N values is traced back to 
  galaxies with early infall times due to the evolution of the 
  M_*-M_ vir relation with redshift.
  GCs that have been tidally removed from their galaxies form the present-day
  intra-cluster component, from which about ∼ 60% were brought
  in by galaxies that currently orbit within the cluster potential. The
  remaining “orphan” GCs are contributed by satellite galaxies
  with a wide range of stellar masses that are fully tidally disrupted
  at z=0. This intra-cluster component is 
  a good dynamical tracer of the dark matter potential. 
  As a consequence of the accreted nature of most
  intra-cluster GCs, their orbits are fairly radial with a predicted
  orbital anisotropy β≥ 0.5. However, local tangential motions 
  may appear as a consequence of localized substructure, providing a 
  possible interpretation to the β<0 values suggested in observations 
  of M87.



galaxies: clusters: general – galaxies: haloes – galaxies: star clusters: general








§ INTRODUCTION

Galaxy clusters are populated by tens of thousands of globular clusters (GCs)
that distribute around galaxies as well as in intra-cluster space. 
However, little is known about their assembly history and their connection to the
cluster build up within the cosmological framework. GCs are among the 
densest stellar systems in the Universe, with typical stellar masses in 
the range [10^4-10^6  M_⊙] and sizes of only a few parsecs <cit.>.
 Because of their inferred old ages <cit.> and 
clustered distribution around galaxies, they are believed to be surviving 
probes of the star formation in the early universe.



GCs are, however, more metal poor than the bulk of the stars in
their host galaxy. In fact, GCs display a wide
range of colors and metallicities that suggests different origins.
One possible interpretation is that the youngest and more metal rich GCs are formed in-situ
while the metal poor component is likely acquired hierarchically via
the accretion of smaller galaxies
<cit.>.





This is consistent with the exquisite GCs data around the Milky Way 
(MW), other galaxies in the Local Group <cit.> and in massive ellipticals
<cit.>. In terms of numbers, however, the rich 
environment of nearby galaxy clusters, due to their high matter density and 
wide range of galaxy masses and morphologies, offer
the best opportunity to understand the connection between
GCs, galaxies and ultimately their dark matter halos. 


Significant effort and resources in the community
have focused on the generation of comprehensive GC maps in nearby galaxy clusters 
such as Virgo <cit.>, Fornax <cit.> 
and Coma <cit.>.
Along with other studies, these data have raised a number of interesting discoveries that 
shed light on the formation of GCs. First, the relation between the total mass in 
GCs (M_ GC) and the estimated dark matter halo mass (M_ vir) for galaxies
follows a single and relatively tight power-law correlation <cit.> that contrasts the complex and non-linear 
correlation of galaxy stellar mass (M_*) and halo mass suggested by abundance matching models 
<cit.>. 
Second, the specific frequency of GCs (S_N),
defined as the number of GCs normalized to the galaxy luminosity, is also highly non-linear, 
with dwarf galaxies in clusters that show high S_N ≥ 5 values which are at odds with
dwarf galaxies of similar mass located in the field and that are characterized by S_N ∼ 2 <cit.>. 

Some of the observed trends in GCs can be explained in models where they trace the most turbulent, high
density and gas rich star formation episodes in galaxies <cit.>,
and analytical calculations support this view <cit.>. 
Other scenarios where GCs are placed at the centers of their own dark matter halo 
<cit.>
and form completely independent of their host galaxy are, although compelling, currently disfavoured 
due to the large degree of tidal stripping expected 
and the observationally constrained abundances and radial distribution of GCs around galaxies  <cit.>. 

While other mechanisms
such as baryonic flows in the early universe may suggest the possible formation of 
GC-like objects <cit.>, 
the  most dominant formation channel for GCs seems therefore to be related 
to galaxies undergoing high-density and violent star formation activity. 
Additional observational support in this direction comes from the observed 
correlation between star formation density and an increased fraction
of stars born in bound stellar clusters in nearby star forming galaxies
<cit.>. 

Due to the low mass and small sizes of GCs, resolving these
systems in hydrodynamical simulations of galaxy formation is a daunting task, 
especially when embedded within the cosmological framework. 
  Encouragingly, full cosmological zoom-in simulations are starting to resolve
 the formation of individual GC-like objects <cit.>. 
 However, due to the expensive calculations the runs are necessarily stopped at
high redshift (z=5) or they focus on very low mass dwarfs, making them 
not fully suitable for predictions of GCs in nearby galaxies and groups in the local Universe. 
Several theoretical efforts are currently underway
that couple the formation and evolution of stellar clusters to the star particles formed in cosmological
hydrodynamical simulations <cit.>. However, the price to pay for such an approach
is high, as it requires resolving 
the necessary conditions for molecular clouds formation and their posterior
tidal evolution at the sub-pc level. Currently this is only possible at the scale of
individual ∼ L_* and dwarf galaxies, but not for larger systems
such as groups and clusters. To take advantage of the richest set of observational
data of GCs available to date, which is collected from all nearby galaxy clusters, 
one must therefore look for alternative and numerically more efficient techniques. 

One possibility is to study the abundance and chemical properties of GCs by following 
their formation, evolution and destruction via semi-analytical modeling built on top of
cosmological N-body simulations <cit.>. 
Such models are powerful since they make quick predictions over a wide range
of halo masses and assembly histories by simultaneously sampling the large parameter
space that is inherent to following the formation of GCs. But it is 
also desirable to make predictions on the spatial distribution or kinematics
of these simulated GCs. A possible avenue to achieve the latter is by 
implementing tagging
techniques, where GCs are “painted” on top of other particle types in the simulation,
i.e. dark matter or stars, and use their dynamical evolution to self-consistently follow
the phase space information expected for the GCs within the hierarchical assembly of
structures in ΛCDM.
Such an approach, reminiscent of the insightful <cit.> technique
used in stellar halo analysis, has been implemented to study the stripping of GCs 
within clusters <cit.>
and the large specific frequency of dwarf ellipticals in clusters <cit.>. 

In this work, we apply a similar technique to that presented in <cit.>,
but instead of using  a single N-body only zoom-in simulation of a galaxy cluster, we couple it to
9 galaxy clusters selected from the hydrodynamical cosmological simulation Illustris <cit.>. 
 
 For the first time, this allows us to make predictions on the connection between GCs
 and the properties of galaxies, for which the stellar mass build up and gas evolution
 is naturally followed by the hydrodynamical treatment
 in the simulations.
 
 Our work is different from that presented in <cit.>
 (which also used Illustris) in that the abundances of GCs are calculated from the
 dark matter halo mass by calibrating
 the model to the observed M_ GC-M_ vir relation presented in <cit.>
 instead of modeling the GC formation by using the individual star formation histories of 
 galaxies as given by the simulation. Furthermore, our model extends that of 
 <cit.> in that all galaxies ever entering the galaxy clusters 
 are tagged with GCs rather than looking only at surviving dwarf galaxies
 and
 
 that the red/blue GC colors are assigned following observational fractions of <cit.>. 
 
 Our technique delivers GC maps in galaxy clusters that can be directly compared
 to observations in Virgo and the Fornax cluster. In this paper, the first of a series, 
 we will introduce the technique (Sec. <ref>) and study the resulting specific frequency
 S_N of galaxies (Sec.<ref>) as well as the predictions for the radial distribution
 and kinematics of the intra-cluster GC component (Sec. <ref>). We conclude by summarizing
 our main results in  Sec. <ref>.  









§ METHODS



The Illustris suite of simulations consists of a series of 
cosmological boxes 106 Mpc on one side that follow
the assembly of tens of thousands of halos and galaxies from z=127 until present day
<cit.>. The simulation adopts a cosmology consistent
 with the  Wilkinson Microwave Anisotropy Probe (WMAP)-9 measurements 
<cit.> and is run using 
the arepo hydrodynamical code <cit.>. Cooling, star formation, 
stellar feedback and stellar evolution are modeled following the implementations described 
in <cit.>, while black hole growth and feedback is modeled following
<cit.>. Relevant to the science presented in this paper, Illustris
has been shown to reproduce several of the observables at the galaxy population
level, including the scaling relations and angular momentum content <cit.>, 
the estimated merger rates <cit.>, general morphology
diversity <cit.> and color evolution of satellite
galaxies <cit.>, among others. 

The Illustris suite is composed of three different resolution levels and sibling
runs with/without the inclusion of baryons. We use for our study the highest 
resolution run (Illustris-1 referred to as Illustris hereafter for simplicity) 
achieving a mass per particle
m_ p, g = 1.3 × 10^6   M_⊙ and 
m_ p, dm = 6.3 × 10^6   M_⊙ for gas and dark matter, 
respectively, and a gravitational softening that is kept always smaller than 700 pc
for the collisionless components (dark and stars). Note that the high resolution 
gas typically resolves much smaller scales due to the adaptive Voronoi mesh construction,
reaching a typical cell size ∼ 50 pc at the centers of galaxies. 

We use galaxy and halo catalogs from subfind <cit.>
which run over pre-identified Friends-of-Friends (FoF) groups <cit.>. We
will use the commonly accepted terminology of “centrals" vs. “satellites" to refer either
to the main galaxy siting at the center of the gravitational potential (a.k.a central)
or to the rest of the galaxies associated to the same FoF group (satellites).
Further details on the simulations data products can be found in <cit.>. 



 §.§ Galaxy sample and infall time


We analyze 9 galaxy clusters from Illustris selected to have present day virial mass
M_ vir≥ 8 × 10^13  M_⊙. In what follows, we use 200
times the critical density of the universe to define all virial quantities. 
Although 10 objects satisfy our criteria at z=0, group 5 is not included
in our sample due to problems in tracing the central galaxy of the group backwards
in time. Our sample then comprises 9 galaxy clusters with virial masses consistent with
those of Virgo and Fornax. 

We use the SubLink merger trees <cit.> to follow all
galaxies that ever interacted with each of our clusters. For GC tagging, we choose
all galaxies with stellar mass above M_* =10^8   M_⊙ at infall time, 
ensuring that our galaxies are numerically well resolved (≥ 60 stellar particles). 
We define the infall time as the last snapshot where a galaxy is a central of its own
FoF group. We note that such a definition does not necessarily correspond to the time
when they join the cluster but can instead be earlier if a galaxy joins a different
group before infalling onto the cluster. However, choosing a different infall time
will not impact our results significantly, as the dark matter halo of any galaxy
(which in our model is what determines the amount and distribution of GCs) stops
growing after they become satellites of any system.
 The central galaxy of each cluster requires a different treatment as their 
infall times are not well defined. 
We choose to tag them at the time when they reach 5% of their final virial mass 
in an attempt to capture their individual GC population separated from that resulting from the later hierarchical
assembly of the clusters. With this definition, we find z∼ 2.5-3.0 as the 
typical tagging redshifts for our centrals. We have explicitly checked that our 
results are not sensitive to this particular choice.
On average, our selection
criteria results in ∼ 300 galaxies with GCs tagged per galaxy cluster. 




 §.§ GC Tagging

   
  In order to select dark matter particles as tracers of globular cluster systems we use the  technique outlined by <cit.> and <cit.>. 
  This technique selects a subsample of particles that follow a chosen spatial density profile, and it has been used by the above authors to simulate stellar populations in dark matter only simulations.
  This method has been already successfully implemented to model the GC systems of galaxies in dark matter only simulations by our team  <cit.>. Briefly, the procedure is as follows.

 First we find the best Navarro-Frenk-White <cit.> fit that 
 describes the density profile of the dark matter in a given infalling halo:


    ρ_NFW(r)=ρ^0_NFW/(r/r_NFW)(1+r/r_NFW)^2

 where ρ^0_NFW and r_NFW are the dark matter halo characteristic density and scale length, respectively. We assume that r_NFW can be approximated by r_max=α r_NFW <cit.>, where α=2.1623 and r_max is the radius of the maximum circular velocity, a parameter provided by the subfind catalog in the Illustris database.

 As in <cit.> and <cit.>, we choose the density profiles of globular cluster systems to follow a <cit.> profile:

    ρ_HQ(r)=ρ^0_HQ/(r/r_HQ)(1+r/r_HQ)^3

where  ρ^0_HQ and r_HQ are the <cit.> characteristic density and scale length, respectively. 

 The radial distribution of GCs around galaxies have been studied in observations <cit.> and shown to be spatially more extended than the stars but less extended that the inferred underlying dark matter halo <cit.>. Furthermore, for the Milky Way galaxy there seems to be a dichotomy of GC properties: 
a more centrally concentrated distribution of red, younger and more metal rich GCs, and a more extended old and metal poor population of “blue" GCs with presumably accreted origins. Interestingly, this dichotomy or bi-modality may not hold true in other external galaxies, which instead may be better explained by a continuous gradient on their GC population <cit.>. 

Regardless of this dichotomy vs. gradient discussion – which is in itself worth of in-depth theoretical investigations that is beyond the scope of this paper–  observational evidence suggests that the redder GC component tends to be more centrally clustered compared to the bluer GCs. In this work, we attempt to model this spatial segregation by tagging two different GC populations, “blue” and “red", with different radial distributions. Notice, however, that this does not require the GC systems to be bi-modal, but instead can be interpreted as a rough description of the “redder" and “bluer" parts of a continuous distribution.

In Eq. <ref> take r_HQ=β r_NFW, with β=3.0 and β=0.5 for blue and red globular clusters, respectively.  The specific values for β were obtained as least square fits to the projected density profiles of red and blue GCs in observations at z=0 <cit.>.

 Once we have the right scale lengths for these density profiles, we compute the distribution function for dark matter (f_NFW) and globular cluster systems (f_HQ) using the equation


    f(ϵ)=1/8π[∫_0^ϵd^2ρ/dψ^2dψ/√(ϵ-ψ)+1/√(ϵ)(dρ/dψ)_ψ=0]


where ψ and ϵ are the relative gravitational potential and the relative energy, respectively <cit.>.

This leaves us in the position to select, in bins of relative energy (ϵ,Δϵ), a fraction f_HQ/f_NFW of particles as candidate GC particles.

To guarantee that the GC system is more concentrated than the dark matter halo as inferred from observations, 
we truncate the selection to a radius r_CutOff=r_50,h/3, where r_50,h is the half mass radius of the dark matter halo.  This choice has been used in previous work and it is inspired on the typical radius for the Milky Way GCs <cit.>. 







 §.§ GC Mass


In order to assign masses to particles tagged as GCs, we use <cit.> as observational constraints.
These authors found that, at z=0, the total mass of globular clusters inside a galaxy can be described by


    M_ GC,0 = a   M_ vir,0^b

 
where M_ GC,0 and M_ vir,0 are the mass of GCs and the halo mass of the galaxy at z=0 in solar mass units (M_⊙), respectively, and taking (a,b)= (4.9 × 10^-5, 0.96) for blue and (a,b)=(2.6 × 10^-8,1.2) for red GCs. <cit.> estimate M_ vir,0 from the stellar mass of the galaxies using the empirical star-halo mass relation from lensing results by <cit.> (see their Appendix C, eq. C1).

We use this result to compute the GC mass at infall in our simulated galaxies using:

    M_ GC, inf = a_ inf   M_ vir, inf^b_ inf


where M_ GC, inf and M_ vir, inf are the mass of GCs and the halo mass of the 
galaxy at infall time in solar mass units, respectively. 
Since our definition of infall time ensures that the galaxy 
is a central galaxy of a FoF group at this time, we take M_ vir, inf = M_ vir.

The values of a_ inf and b_ inf in Eq. <ref> need to be calibrated so that, 
at redshift z=0 the resulting correlation between M_ vir,0 and M_ GC,0 reproduces
the observed normalization from <cit.>. This exercise gives us
a_ inf=3.5×10^-4  (2.0×10^-7) and b_ inf=0.9   (1.15) for blue (red) GC, respectively.  This calibration is needed to account for the average fraction of GCs lost by 
each galaxy to tidal stripping in the cluster potential.
We note that because at z=0 the
simulated galaxies are satellites within the cluster (and therefore their halo/virial
mass is ill defined) to make a fair comparison to <cit.>, we compute
the z=0 halo mass using the same procedure as the authors and implement the 
star-mass relation from <cit.>
using the stellar mass of the galaxies at z=0 as provided by the simulation. 






  For simplicity, 
we have chosen an infall relation that is independent of redshift
(a_ inf and b_ inf in Eq. <ref> 
are the same for all infalling galaxies). Theoretical models
give different predictions regarding a possible evolution
of the M_ GC-M_ vir relation with redshift, 
from nearly no evolution <cit.> to
changes by factor ∼ 10 <cit.> between z ∼ 3
and present day. The lesson learned from our model is that such
evolution is not needed based on available observations, 
besides a calibration that accounts for the average fraction of GCs 
lost to tidal removal by the cluster. 
However, this does not mean that normalization of the M_ GC-M_ vir
relation cannot change with time 
<cit.>.
It only means that, given the fair amount of free parameters in our model, 
we choose to adopt an unevolving infall relation and defer
an exploration of changing this hypothesis to future work and plausible 
observational evidence that may justify better such assumption.

Blue (red) points of figure <ref> shows the resulting blue (red) GCs mass associated to 
galaxies in sample as a function their halo mass at z=0.
 We define a GC as
“associated” to a galaxy if it lays within three times the stellar
half mass radius of a galaxy, i.e., r< 3r_h,*. This choice is inspired by the typical values used in observations <cit.>. However, replicating the exact procedure done on real galaxies is non-trivial and involves the fitting of parametrized profiles to account for an uniform background of not-associated GCs. Instead, here we adopt a simpler definition that scales the size with that of the central galaxy. We have explicitly checked that our main conclusions are not changed significantly by this choice by varying the criteria of association between 2-4 r_h,*. 


Blue and red lines in Fig. <ref> indicate the observational
relation by <cit.>. For comparison, we also show in green the relation reported by
<cit.> for a combination of all GCs, and shows that our results are in the ballpark of
those observations as well. Once M_ GC is known for each of our galaxies, we take the number of
candidate dark matter particles associated to each object, N_ tag, and assign a mass 
M_ GC, inf/N_ tag to each tagged particle. In this way, all GC particles with the same 
original host galaxy have the same mass.


With this method, each particle tagged as a GC has usually a mass of ∼ 10^3 M_⊙, 
far below the typical mass of a observed GCs. This happens because the number of tagged 
candidate particles that are consistent with the right distribution function of the GCs distribution
is greater than the expected number of real GCs in galaxies. Although this procedure allows us 
to have a better sampling of the phase-space of GCs around each galaxy, there are certain
applications where the number of GCs is important (for instance, specific frequency estimations). 
In those cases, it is not adequate to use all the candidate GCs to compare to observations, 
but instead we re-sample. This is done by adding back all the GC mass (previously distributed among all
tagged particles) and dividing by the typical GC mass chosen to be 
m_ GC=10^5 M_⊙ <cit.>.
This gives us the number of particles that need to be considered per galaxy as “realistic" GC numbers.
We then randomly select from all candidate GC particles, the needed number of realistic GCs. 
Fig. <ref> shows as an example the most massive galaxy cluster in our simulations
indicating in gray scale the stellar density and highlighted in magenta the tagged GC candidates.








 One of the caveats of our method is that the early destruction of GCs due to
tides, evaporation and infant mortality is not explicitly followed. Instead, our
model is calibrated to reproduce the z=0 content of GCs around galaxies, meaning that
those GCs that were disrupted after birth are never part of our sample. This is a necessary compromise, 
since the modeling of such destruction processes can be not only cumbersome, requiring many free parameters
<cit.>, but also requires a numerical resolution that is currently prohibitive for the study of GCs on the scale of galaxy clusters. Consequently, our method is not able to inform us about the early destruction mechanisms of GCs and it has been designed, instead, as a tool to study the predictions of the surviving population at z=0 within the cosmological assembly of structures. 









§ GCS AROUND GALAXIES IN CLUSTERS



 After tagging the corresponding GC candidate particles, our technique 
 self-consistently follows the orbits and tidal stripping experienced
 by each galaxy and their population of GCs after they join their
 cluster. As a result of that, the number of GCs that remain
 associated to a given galaxy today may differ from the one that was
 originally assigned. We start by quantifying the distribution and number
 of GCs  around galaxies in our 9 simulated galaxy clusters in Illustris-1.

Fig. <ref> shows a zoom-in of the most massive
galaxy cluster within half the virial radius. The background gray
scale map corresponds to the stellar distribution. Superimposed in blue
and red are the tagged candidate GCs. A close inspection of this
figure indicates that most GCs are found today heavily clustered
around the central and satellite galaxies. As introduce in Sec. <ref>, GCs are considered associated to a galaxy if they lay within 3r_h,*. We indicate this criteria with a green circle. 

We note in Fig. <ref> that some small
stellar clumps are visible in the grayscale map and have no associated
green circle. These correspond to galaxies that remained below the
minimum mass for tagging (M_* < 10^8   M_⊙) throughout their
history and will consequently not be considered in the analysis that
follows. In a few exceptional cases a galaxy will be above our minimum
mass criteria but are identified by the merger tree as having been a
satellite of another system at all times. Because the tagging is done
at infall, defined as the last time that a galaxy was a central of its
own group, these objects cannot be tagged and therefore will also be
excluded from the analysis (an example of it can be seen near the
(-50,-250) kpc position in the figure, with a few blue GCs that were
acquired as some smaller companion merged to it). The fraction of such
objects in our sample is below ∼ 2%  of the galaxies and
therefore will not strongly impact our results.





 §.§ The GC specific frequency S_N


In observations of galaxy clusters, the specific frequency of GCs
S_N, which is defined as the number of GCs per unit M_V=-15 of
galaxy luminosity, shows a strong dependence on
galaxy brightness, with a characteristic “U”-shape indicating large
specific frequencies in bright as well as in dwarf galaxies
<cit.>. The origin of this behaviour is not well
understood. In particular, in light of the results of field and late type galaxies,
which suggest S_N ∼ 1-2 <cit.>. 
Although the large specific frequencies of the
central galaxies may not be surprising given their ability to capture GCs
from orbiting and merged galaxies in the cluster <cit.>; the
upturn of S_N in the dwarf scale regime is less well understood. 

We use our tagging technique to gain insight into the possible cause
for such trends. Fig. <ref> shows S_N as
a function of absolute V-band magnitude for our simulated galaxies.  S_N is
computed as,


    S_N = N_GC× 10^0.4(M_v + 15)


<cit.>, where N_GC is the number of
GCs within r< 3r_h,* of a galaxy (and assuming an average mass per
cluster 10^5   M_⊙, see
Sec. <ref>). Encouragingly, the GC specific frequency in
our simulated galaxies, shown in black small dots, follow a similar the familiar
trend found in observations. For comparison we include in
magenta starred symbols a set of observed S_N values for galaxies in
the Virgo cluster from <cit.>. Although our simulations tend
to underestimate S_N on the bright end, an effect driven by the
inefficiency of feedback to shut down star formation in high mass
halos for galaxies in Illustris <cit.>  artificially lowering our
S_N calculations, the good agreement with intermediate and low mass
dwarfs is a notable success of the GC model.

Our GC tagging method relies mostly on one specific assumption: a single
power-law relation between the mass in GCs and the halo virial
mass at infall. The S_N characteristic “U”-shape in our model
arises as a consequence of this assumption combined with the well known
non-linear relation between stellar mass (and, consequently,
luminosity) and halo virial mass.  Notice that for the latter we use
the information on stellar mass and luminosity directly from the
Illustris hydrodynamical run and it is not an input of our GC model. 
This conclusion is consistent with previous independent results 
using semi-analytical catalogs or analytic models of GCs formation
<cit.>

To guide the eye, we add in Fig. <ref>
yellow curves that are the results of a toy model, as follows. We vary
the halo mass M_ vir in the range
log(M_ vir/M_⊙)=[10, 15]. For halos in this mass range,
following the M_ GC - M_ vir relation from <cit.>, we
determine the mass in GCs. As a second step, assuming three different abundance
matching relations M_*-M_ vir <cit.> we estimate
their M_*. Adopting a uniform mass-to-light ratio γ=2 to go from
stellar mass to V-band magnitude and an average GC mass
M_ GC∼ 10^5  M_⊙, we compute the resulting S_N for
a galaxy with M_V magnitude. This simple model seems to capture the
essence of the simulated and observed trend; thus, it provides a possible
interpretation for the S_N shape. 


Our simulations also shed light on the origin on the scatter
of the S_N-M_V relation.  The S_N scatter increases 
substantially towards the low mass galaxies in our model, 
in good agreement with observational findings <cit.>. 
In particular, we note the existence of several dwarf galaxies that
at z=0 have no associated GCs. A closer inspection of those objects reveal
that their full GCs content has been tidal stripped by the gravitational 
potential of the host cluster. A more detailed analysis of the scatter
in the S_N-M_V relation suggests differences in the infall properties of
surviving galaxies at z=0. Fig. <ref> shows 
that for a given stellar mass today M_*, the corresponding virial
mass at infall can vary substantially (and therefore, the initial
mass of GCs assigned). As expected, objects with a larger halo mass at
a given present-day stellar mass should display larger specific frequencies
today, as confirmed by the color coding of the points. 

Interestingly, we find that the bulk of the scatter in this plot
is due to an evolving M_*-M_ vir relation with redshift in our
simulations. To illustrate this, we show the median M_*-M_ vir 
relation of those points selected in very narrow ranges of infall times:
t_ inf=[2.5,3.5] Gyr (green) and t_ inf=[13,13.5] Gyr (black). 
For halos infalling earlier into the cluster, the same halo mass is populated
by a smaller stellar mass content, resulting in an increased S_N value today.
Such evolution in the M_*-M_ vir relation seen in our simulations 
is also consistent with that determined via abundance matching studies
<cit.>. 

We have checked that, for those points well above the green curve in 
Fig. <ref>, tidal disruption plays a major role. 
Gray circles highlight objects where more than 90% of the stars 
have been tidally removed.
In most cases tides remove stars and GCs and are conducive to a
small S_N (blue points). However, there is a small subset of heavily
tidally stripped galaxies where the specific frequency increases
to S_N ≥ 8. These galaxies, located in the dwarf regime
in Fig. <ref>, have retained less than 10% of their stellar
budget but have held on to at least 10 GCs resulting in large specific
frequencies at z=0. 
 The distribution of stars follows naturally a shallower profile 
than the cuspy Hernquist profile selected for GCs. Therefore, 
once a galaxy has been significantly depleted of stars, it may enter a 
regime where more stars are lost than GCs, explaining the rise in S_N.
However, this might be a direct consequence of
our specific choice of a cuspy radial profile for the tagged GCs
and we prompt the reader to consider S_N ≈ 10 values and above 
with caution.










Although our assumption of an almost constant slope in the 
M_ GC-M_ vir relation with redshift partially drives 
some of these results, it is reassuring to see that
estimates of GC formation as a result of mergers might support such an 
assumption <cit.>. Moreover, our model makes a
testable prediction: dwarf galaxies with the largest S_N values
within galaxy clusters will correspond to the earliest infalling halos
and therefore should be located preferentially in the inner regions of
clusters. We note that our results are in line with those presented in
<cit.> with an independent GCs analysis in the Illustris
clusters. <cit.> also conclude that
earlier infalling dwarfs have larger S_N values, an effect partially
driven by their larger dark matter halo mass –in agreement with our
results– and that more bursty star formation may increases the
chances of GCs formation. Encouragingly, observational evidence seems
to suggest an increase of S_N for dwarfs near the cluster center
<cit.>, in good agreement with our
predictions.

We hasten to add that our prediction for field dwarfs would place them
at S_N ∼ 5 at M_V=-16 (i.e. near the green curves in
Fig. <ref>) but they will not scatter upwards. We
defer the comparison to field objects to future work.









§ DISTRIBUTION & KINEMATICS OF INTRACLUSTER GCS



GCs that lay beyond 3r_h,* from the center of a galaxy (and are
therefore not considered in the S_N calculations of the previous
section) constitute a diffuse and extended intracluster GC (ICGC)
component. They are seemingly not connected to any galaxy today but
are instead linked to the gravitational potential of the host galaxy
cluster. In observations, hints to such intracluster GC population can
be found in early work <cit.> 
and confirmed in more recent observations of nearby
galaxy clusters such as Virgo <cit.>, Fornax
<cit.> and Coma <cit.>. We use our simulated
GCs catalog to study this diffuse intracluster component, understand
its origin and determine their potential as dynamical tracers of the
dark matter.








 §.§ Building the intracluster GCs component


Fig. <ref> shows the projected radial profile of red and
blue intracluster GCs. Thin lines indicate individual simulated galaxy
clusters in Illustris whereas thick solid curves represent the
median. ICGCs have a radial distribution that roughly follows that of
the intracluster light (cyan curves) and that are steeper than the
dark matter (shown in gray), with a surface density profile scaling
Σ_ GC∼ r^-2.5. A complete census of this diffuse
component in observations is extremely challenging, and as a result
only a few estimates are available to date. For reference, we compare in
Fig. <ref> with the SDSS measurement for ICGCs in Virgo by
<cit.> shown with starred symbols. 

The blue GCs intracluster component in our simulated clusters seems in
reasonable agreement with <cit.> although we note that our
model predicts a shallower profile for the red (or metal rich)
distribution than observed. This is likely due to an insufficient
initial segregation between red and blue components already set at
infall time of individual galaxies in our GCs tagging method. However,
given the non-linearity required for building the ICGC and the
simplicity of our tagging method, the good agreement between the
radial distribution of observed and simulated GCs is  encouraging. To be
conservative, in what follows, we will not distinguish between red and
blue ICGC components and consider it as a single population.





What is the origin of this ICGCs? We find two main contributors to
building the ICGC component, the GCs that were tidally stripped from
(surviving) galaxy satellites in the cluster as well as GCs that were
brought in by galaxies that had fully merged to the central
potential. We will refer to them in what follows as “orphan” GCs
since their progenitors not longer exist today. Fig. <ref>
allows a closer look to this contribution as a function of radius
from the cluster center. We show the 3D median density profile of the
ICGC component in a stacking sample of our simulated
clusters. ICGCs are shown in dashed green and it is split into
contributions from orphans (blue dashed-dotted) and those that were
originally assigned to the central galaxy (dotted red). The remaining
contribution to the ICGCs, accounting for about 62% of this component,
originates from GCs initially associated to galaxies that still survive today 
as satellites in the cluster[The exact contribution of surviving satellites 
to the ICGC component is only weakly dependent on the adopted cutoff of 3r_h,* 
to consider a GC as not associated to a galaxy. For example, this changes to 
54% if the cutoff is extended to 5r_h,* instead]. Vertical
arrows indicate the half number radius split by components.

GCs initially associated to the central galaxy dominate only in the
inner region and their contribution is not significant beyond
r/r_ vir≥ 0.1. From the rest of the intracluster component,
orphans are more centrally concentrated compared to those from
surviving galaxies. This is consistent with orphan GCs having having had
their progenitor galaxies totally disrupted, which explains their
segregation towards the center of the cluster. For completeness, we
show in Fig. <ref> the total 3D median profile of GCs
(black). The difference between ICGCs and the total curve indicates,
as expected, that in the outskirts of the galaxy clusters most GCs are
still associated to their host galaxies.






Because the chemical composition of GCs is linked to that of their
progenitor galaxy <cit.>, 
it is of interest to identify the spectrum
of galaxy masses contributing GCs within a galaxy
cluster. The top panel in Fig. <ref> shows the 
added distribution of progenitor galaxy masses for GCs in 
all our simulated clusters, divided according to the different populations. 
Taken as a whole, most GCs today in clusters are linked to massive galaxies with
M_* ∼ 10^10.5  M_⊙ and above (black solid curve), 
which seems to also apply for the intracluster component (green line). 
We have checked that this is the case when we look individually to each of the simulated clusters. 

On the other hand, the orphan population seems quite diverse.
Although adding all clusters hints to a highly skewed 
distribution towards massive progenitors for the orphan population
(blue line, top panel), when analyzed separately, each cluster
may show a different behavior. We demonstrate this in the bottom panel of
Fig. <ref>, which shows for each cluster the progenitor
stellar mass distribution for their orphan population. Because the number
of orphans may vary from object to object, each histogram
has been normalized such that the bins add up to 1 to ease the comparison. 
We find that about half of the galaxy clusters are heavily dominated by
orphan GCs brought in by one or two very massive M_* ≥ 10^11  M_⊙
progenitors (major merger events). However, for the other half of the clusters, 
the contribution to the orphan GCs population is dominated instead by 
dwarf-like objects, with median masses M_* ∼ 10^9.5  M_⊙. 
Such variations should be imprinted in the chemical signatures of 
the ICGC component, providing clues to unravel the past accretion
histories of their cluster hosts. 

Therefore, our simulations indicate an interestingly broad distribution of galaxy masses 
that contributed to the build up of the ICGC component. 
As such, a heterogeneous range of metallicities and
ages are expected to be found in observations, albeit modulated by the radial
trends shown in Fig. <ref>. We note that histograms in
the top panel of Fig. <ref> are number-weighted, meaning that each
candidate GC counts individually and the same galaxy progenitor is
counted as many times as GCs it contributed. For GCs still associated
to galaxies (contributing to the all curve in black), we
consider the present-day mass of their host galaxy. In the case of
orphans GCs, we count the mass of the galaxy progenitor at their
infall time, since they have totally merged to the cluster by z=0 by
definition.



 §.§ Intracluster GCs as kinematic tracers of the dark matter






This population of free floating GCs provides also a unique opportunity
to trace kinematically the gravitational potential of the galaxy
cluster in regions where stars or other luminous tracers are scarce.
Obtaining spectroscopic data of several dozens to hundreds of GCs in
the intra-cluster medium of nearby clusters like Virgo or Fornax
is, although expensive, within current to near-future capabilities of
observational campaigns (e.g., Subaru/PFS, Maunakea Spectroscopic Explorer). 
In the MW itself, 3D information of GCs have been used to estimate the enclosed
dark matter mass out to ∼ 100 kpc, well beyond the stellar and
gaseous disk <cit.>. However, for systems outside our own Galaxy, the
kinematical information will be projected. And how well GCs are
expected to trace the gravitational potential will then depend on
relative biases between dark matter and GCs as well as the orbital
structure predicted for this ICGC component.


Fig. <ref> compares the global kinematics of the ICGCs
(green) to that of the underlying dark matter halo (black) for our 9
simulated clusters. We show the three spherical components of the
velocity in a spherical coordinates system centered at each cluster. For
this, we have stacked the information of all individual clusters by
scaling the distances and velocities by their respective virial
quantities. The results from individual clusters are not different
from the overall ensemble. We find that the velocity dispersion for
the ICGCs in all three directions is comparable to that of the dark
matter, albeit systematically smaller by ∼ 20-25%. The more
concentrated radial distribution found for GCs compared to the dark
matter in Fig. <ref> is a possible explanation for this
difference. Values for the r.m.s dispersion in all coordinates are quoted individually in each panel for both components to ease the comparison. In all cases, the distributions are well described by a Gaussian function with the measured r.m.s dispersion (indicated by the thin lines).

The radial velocity component for the ICGC as well as for the dark
matter exceeds that of the θ and ϕ components, suggesting
the prominence of radial orbits over tangential ones. This is shown in
more detail in Fig. <ref> by means of the orbital anisotropy
parameter β = 1 - 2σ_r^2/σ_t^2, where 
σ_t^2 = σ_ϕ^2 + σ_θ^2 <cit.>. 
Thick solid lines show the median β profiles for the ICGC 
separated in red and blue component (although no significant 
difference is seen between both populations). 

We find that the orbits of GCs are expected to be highly radially
biased (β≥ 0.3) at all radii within the cluster. 
This agrees with the accreted nature of the intracluster population, where galaxies in radial orbits are preferentially tidally disrupted, thus donating their GCs to the cluster potential. 
Dark matter halos in
ΛCDM galaxy clusters have long been known to be radially 
biased <cit.>.
Fig. <ref> also shows the median β profile of the host
dark matter halos (black curve) and confirms that the ICGC component
is comparably and even slightly more radially biased than the dark
matter in the clusters. Large β values therefore seem an
unavoidable prediction for any scenario where structure grows
hierarchically and tidal stripping plays a mayor role
<cit.>. 

A few observational estimates of β are currently available in
the literature and offer important constraints to theoretical models
of GC assembly. 







In particular, the kurtosis of the line
of sight velocity dispersion in combination with an estimate of the
circular velocity can be used as an indication of the
orbital anisotropy of a luminous tracer (stars, PNe, GCs) in extragalactic objects
<cit.>. 
Such measurements are particularly challenging due to projection effects and the availability of a discrete number of GCs.

Applying the kurtosis technique to GCs around M87 suggests
that their orbits are radially biased (β > 0), consistent with our predictions (magenta triangles in
Fig. <ref> compared to solid lines).  






Interestingly, an
independent analysis of the GCs around M87 presented in
<cit.> argues for a β<0 beyond r ≥ 100 kpc.
Such preferentially tangential orbits for GCs are difficult to
reconcile with the accretion origin for the ICGC component explored
here <cit.>.

As discussed in <cit.>, a possible origin for the
negative anisotropy values would be the
preferential tidal disruption of GCs in radial orbits, which would
create an orbital distribution that favors preferentially tangential
orbits. Our method is not well suited to test such hypothesis since
our model is calibrated to reproduce the final GC content in
halos and the destruction of individual GCs is not explicitly
modeled. Instead, our simulations suggest an alternative explanation
for negative β values.  

Thin blue lines in Fig. <ref>
show the individual anisotropy profile for each of our simulated
clusters and confirm that local deviations toward lower β values
are not uncommon. 




We have associated such dips in the profiles with the presence of infalling galaxies that have not yet virialized into the cluster. The GCs of these infalling objects mostly trace the substructure of the cluster (and hence the tangential motion), and they are not directly tracing the cluster potential.





The substructure explanation suggested by our simulations seems well aligned with GCs on possible
tangential orbits in M87, as the Virgo cluster is known to be not
fully relaxed. Moreover, substructure in the kinematics of planetary
nebulae in the Virgo cluster have already been identified <cit.>.


Our results indicate that localized
tangentially biased orbits are not inconsistent with the accretion
build up of the ICGC component. However, such β<0 values should
be local, not extending for more than a few hundred kiloparsecs, offering
a path to proving or disproving such explanation. 








§ CONCLUSIONS



We present a novel method to follow the evolution of GCs within galaxy
clusters in cosmological simulations. The method relies on the
observational single power-law relation between mass in GCs and halo
mass to determine the number of GCs associated to a galaxy. By
following the cluster assembly, we trace back all galaxies at infall
time, and by means of a particle tagging technique, we are able to follow the
posterior dynamical evolution of its associated GCs within the
cluster potential. We apply this technique to 9 of the most massive cluster
halos with M_ vir∼ 10^14  M_⊙ in the cosmological
hydrodynamical simulation Illustris, providing useful constraints to
observations of GCs in systems such as the Virgo cluster.

The successes of this model include reproducing (i) the large
scatter in the specific frequency of GCs, S_N for dwarf galaxies, helping build
the observed “U”-shape in the S_N- stellar mass or magnitude relation;


and (ii) the build up of an intracluster population of
GCs that is spatially distributed following a steeper profile than the
dark matter, in agreement with observations. We find that S_N ∼ 5
values are naturally expected for dwarf galaxies with M_V > -18 as a
result of the non-linear relation between halo mass and stellar
mass. Furthermore, dwarfs that scatter upwards in the S_N-M_V
relation correspond to the early accreting subhalos onto the cluster,
and their large specific frequency stands from the lower galaxy
formation efficiency at fixed halo mass expected at higher
redshifts. This result is in line with observational evidence of large
S_N values distributed preferentially near the center of the Virgo
cluster <cit.>.

We use our simulations to
make predictions about the galaxy progenitors that are expected to be
the main contributors to the intracluster population of GCs. We find
that in numbers the ICGC population is dominated by galaxies
like the MW and above, with masses M_* ∼ 10^10.5  M_⊙
contributing about 80% of the GCs in this extended and diffuse
component. Furthermore, we find that the tidal origin for the ICGC
translate, on average, into a very radial (β∼ 0.5) orbital 
structure for the GCs. However, the anisotropy profile can present localized deviations
towards tangentially biased orbits as a result of surviving
substructure of galaxies within the otherwise virialized cluster. It
would be interesting to explore this possibility in light of the
tangentially biased orbital structure measured by the galaxy group
centered at NGC 1407 <cit.> and for M87 according to
<cit.>.

 An important caveat of our analysis is that the GC content and distribution
for each galaxy is not computed ab-initio, but instead parametrized
based on empirical relations following Eq. <ref> and <ref>. 
This means that only surviving GCs are tagged, providing no valuable insights on
the formation, dynamics and destruction mechanisms of GCs at birth.
Fully self-consistent models of GC formation and their evolution
within the local interstellar medium embedded within the cosmological
context are the ultimate goal to understand the GC connection to halos
and galaxies. However, the wide dynamical range needed to be resolved
(from sub-pc to several mega-pc scales) are outside the reach of current
numerical simulation capabilities. Important efforts to follow the
formation of stellar clusters are starting to shed light on such
connections, but they are limited to the scale of individual galaxies
<cit.>, 
higher redshifts <cit.> or idealized galaxy mergers set-ups
<cit.>.

We argue that tagging techniques like the one presented here, when
coupled to fully hydrodynamical simulations that follow the formation
of galaxies within ΛCDM such as Illustris, offer a unique
opportunity to make predictions on the properties of GCs on the scales
of galaxy groups and clusters  until more sophisticated models become
available. This approach is complementary to the
aforementioned theoretical strives on smaller scales and can
potentially be informed by them if more complex modeling is required.
Our current predictions cover the need for a theoretical framework to
interpret the wealth of existing observational data of GCs in nearby
groups and clusters such as Virgo <cit.>, 
Coma <cit.>, Fornax <cit.> and CenA 
<cit.> among others.




§ ACKNOWLEDGEMENTS

We would like to also thank
the anonymous referee for a detailed and constructive referee report 
that helped clarify and improve the first version of this manuscript.
Authors are grateful for insightful comments and discussions with
Oleg Gnedin, Thorsten Naab, Hui Li and Aaron Romanowski. 
LVS acknowledges support from NASA through the HST Program AR-14583 and 
from the Hellman Foundation. FRA, MGA and HM acknowledge financial support from 
Secyt-UNC and Conicet, Argentina. EWP acknowledges support from the 
National Natural Science Foundation of China under Grant No. 11573002.











mnras
 




























