Properties of Globular Clusters in Galaxy Clusters: Sensitivity from the Formation and Evolution of Globular Clusters

We build the cluster mass assembly history using three sets of high-resolution DM-only cosmological simulations from z = 200 to 0 (Taylor et al. 2019) with the N-body code, gadget2 (Springel 2005). The original simulations were carried out with cosmological parameters Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.04, and h = 0.68 with a box of (140 Mpc h⁻¹)³. The mass of each particle is 1.7 × 10⁹ M_⊙ h⁻¹ with a softening length () of 5.469 kpc h⁻¹. The power spectrum is calculated by the camb package ³ (Lewis et al. 2000).

From these original simulations, Virgo cluster analog DM halos ⁴ are identified at z = 0. We choose three cluster DM halos (Targets 1–3), which have small Lagrangian volumes, among various Virgo cluster analog DM halos in the original simulations to reduce the calculation time. Next, we re-simulate Targets 1–3 with a zoom-in technique, which involves rerunning the interesting regions of the low-resolution simulations with higher resolution (Porter 1985; Katz & White 1993; Navarro & White 1994). Multiscale initial conditions with positions and velocities were generated by the music package ⁵ (Hahn & Abel 2011). The high-resolution particle mass is 3.32 × 10⁶ M_⊙ h⁻¹ with = 0.683 kpc h⁻¹. Targets 1–3 have an M_halo of 9.25 × 10¹³, 1.14 × 10¹⁴, and 1.26 × 10¹⁴ M_⊙ h⁻¹ , respectively. The virial ratio, 2T/∣U∣, of Targets 1–3 is 1.56, 1.19, and 1.31, respectively, where T is a kinetic energy and U is a potential energy so Target 1 is the most unrelaxed cluster among them.

We identify halo and subhalo structures with the modified six-dimensional phase-space DM halo finder, rockstar ⁶ (Behroozi et al. 2013b) and make merger trees using Consistent Trees (Behroozi et al. 2013c). To define the DM halo in rockstar, we set the minimum number to 20 particles so that the minimum M_halo is 6.64 × 10⁷ M_⊙ h⁻¹. Throughout this paper, however, we use DM halos more massive than an M_halo of 10⁹ M_⊙ h⁻¹, which consists of more than ∼300 DM particles. The equivalent M_stellar of M_halo ≃ 10⁹ M_⊙ h⁻¹ would be 5.97 × 10⁵ M_⊙ h⁻¹. Here, We convert M_halo to M_stellar, using the M_stellar and M_halo relation from weak lensing data (Hudson et al. 2015). This match is what was done in the observation so we can compare our results with the observations identically. The slope of stellar MF of the observation is −0.56 (Lan et al. 2016), and Targets 1–3 have slopes of stellar MF with −0.54, −0.54, and −0.57, respectively.

The hierarchical merging scenario is adopted to make the GCs in our model. We assume GCs form by DM halo mergers that are larger than a minimum mass ratio of γ_MR, which is a free parameter in our model. When GCs form, the initial M_GC in each DM halo has a linear relation with M_halo (Peng et al. 2008; Harris et al. 2013; Hudson et al. 2014; Harris et al. 2015):

$$\begin{eqnarray}&&{M}_{\mathrm{GC}}=\eta {M}_{\mathrm{halo}},\end{eqnarray} \tag{ 1 }$$

where η is a constant formation efficiency and a free parameter in our model.

We assume that one DM particle represents one GC, although the mass of the GC is not equal to the mass of the DM particle. The initial M_GCs is distributed into individual GCs with a power-law initial globular cluster mass function (GCMF) (Elmegreen 2018):

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dM}}={M}_{0}{M}^{-2},\end{eqnarray} \tag{ 2 }$$

where M₀ is a normalization constant. We determine the minimum mass of the GC (${M}_{\min }$) is 10⁵ M_⊙ because we assume that most GCs below 10⁵ M_⊙ are destroyed due to two-body relaxation. If we normalize the probability of the MF, Equation (2) is

$$\begin{eqnarray}&&1={\int }_{{M}_{\min }}^{\infty }{M}_{0}{M}^{-2}{dM},\end{eqnarray} \tag{ 3 }$$

which gives ${M}_{0}={M}_{\min }$. Then, using the transformation method, the initial mass of each GC, M_i , is selected by a random number, 0 < r < 1:

$$\begin{eqnarray}&&{M}_{i}=\displaystyle \frac{{M}_{\min }}{r}.\end{eqnarray} \tag{ 4 }$$

We repeat Equation (4) until the sum of the initial mass of each GC in individual DM halos reaches the required M_GCs value at the GC formation epoch, and in this way, the number of GCs is also determined. The same number of DM particles from each DM halo is then selected to trace out the location of the GCs. Here, we pick particles according to in order of their bound energies, which assumes GCs form in the deepest location of the potential potential of DM halos. Finally, a different GC initial mass is assigned to each of the tagged DM particles, allowing us to trace their position and velocity down to z = 0.

We assume that the metallicity of GCs is assigned by the M_stellar of their host DM halos and redshift at the GC formation epoch (Li & Gnedin 2014; Choksi et al. 2018; Choksi & Gnedin 2019b). We adopt the stellar mass–metallicity model (Ma et al. 2016; Choksi et al. 2018):

$$\begin{eqnarray}&&[\mathrm{Fe}/{\rm{H}}]=\mathrm{log}\left\{{\left(\displaystyle \frac{{{M}}_{\mathrm{stellar}}}{{10}^{10.5}\,{{M}}_{\odot }}\right)}^{{\alpha }_{m}}{(1+z)}^{-{\alpha }_{z}}\right\},\end{eqnarray} \tag{ 5 }$$

where α_m = 0.35 and α_z = 0.9 (e.g., Choksi et al. 2018; Chen & Gnedin 2022). Here, we define metal-poor (blue) GCs as [Fe/H] < −1.0 and metal-rich (red) GCs as [Fe/H] > −1.0.

After GCs form, they undergo mass loss and disruption by several internal and external processes: stellar evolution, two-body relaxation, tides from host galaxies, the orbital decay by dynamical friction, and tidal shocks (e.g., Spitzer 1987; Gieles et al. 2011; Shin et al. 2013; Webb et al. 2014; Madrid et al. 2017, and references therein). We include recipes for these processes in our PTM with the semianalytical approach. Note that we do not consider the GC disruption by tidal shocks because our simulations are DM-only simulations and we do not model the gas and stellar distributions in each galaxy.

2.3.1. Mass Loss of GCs by Stellar Evolution and Two-body Relaxation

The GC mass evolution is approximated to the first-order differential equation (Fall & Zhang 2001):

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=-({\nu }_{\mathrm{se}}+{\nu }_{\mathrm{rlx}}){M}_{i},\end{eqnarray} \tag{ 6 }$$

where ν_se and ν_rlx are time-dependent fractional mass-loss rates by the stellar evolution and two-body relaxation, respectively.

For ν_se, we use the stellar evolution model (Hurley et al. 2000), assuming the initial mass function (IMF) of each GC is distributed by a Kroupa IMF (Kroupa 2002). We tabulate the changes in the mass of stars with different masses as a function of time, ν_se(t), so we can calculate the amount of mass loss of each GC after they form. We assume a constant metallicity (0.1 Z_⊙), which is a typical metallicity of observed GCs. ⁷ Figure 1 shows the evolution of GCMF in the brightest cluster galaxy (BCG). The thin black line denotes the MF without a mass loss or disruption of GCs (original MF), and the green line shows the evolved MF, where we only apply the mass loss by the stellar evolution. The stellar evolution does not change the shape of the MF but it moves the original MF to the low-mass part.

Zoom In Zoom Out Reset image size

For ν_rlx, we adopt the formula in Spitzer (1987):

$$\begin{eqnarray}&&{\nu }_{\mathrm{rlx}}=\displaystyle \frac{{\xi }_{e}}{{t}_{\mathrm{rh}}},\end{eqnarray} \tag{ 7 }$$

where ξ_e of 0.033 is the normalization factor and t_rh is the half-mass relaxation timescale:

$$\begin{eqnarray}&&{t}_{\mathrm{rh}}=\displaystyle \frac{{M}_{i}^{1/2}{r}_{{\rm{h}}0}^{3/2}}{7.25\bar{m}{G}^{1/2}\mathrm{ln}\,{\rm{\Lambda }}},\end{eqnarray} \tag{ 8 }$$

where the $\bar{m}$ of 0.87 M_⊙ is the mean stellar mass of a Kroupa IMF, $\mathrm{ln}\,{\rm{\Lambda }}$ of 12 is the Coulomb logarithm, a typical value for GCs, and r_h0 is the initial half-mass radius of GCs. We treat r_h0 as a free parameter in our model.

The blue line in Figure 1 shows the evolved MF after applying the mass loss by the stellar evolution and two-body relaxation. The shape of the original MF is changed to a log-normal shape by two-body relaxation.

2.3.2. Disruption of GCs by Tides and Dynamical Friction

To consider tides from host DM halos where GCs belong, we divided GCs into those in the isolated regime and those in the tidal regime, using ${\mathfrak{R}}={r}_{{\rm{h}}}(t)/{r}_{{\rm{J}}}(t)$, where r_h(t) is a half-mass radius and r_J(t) is the Jacobi radius ⁸ with time (Gieles & Baumgardt 2008). The Jacobi radius is the distance from the cluster center to the Lagrange points L₁ and L₂ so it is used to define the boundary where stars dynamically belong to the GC. The variation of r_h is mainly driven by the internal dynamics of GCs, while r_J is mainly changed by the external dynamics of host galaxies.

Gieles & Baumgardt (2008) find that if ${\mathfrak{R}}\gt {{\mathfrak{R}}}_{c}$ star clusters can be treated in the tidal regime, while for ${\mathfrak{R}}\lt {{\mathfrak{R}}}_{c}$ they can be treated in the isolated regime, where ${{\mathfrak{R}}}_{c}$ is a criterion used to divide into the tidal regime and the isolated regime. Gieles & Baumgardt (2008) found that ${\mathfrak{R}}=0.05$, but we treat ${{\mathfrak{R}}}_{c}$ as a free parameter. In our model, we assume GCs staying in the tidal regime longer than the fixed number of t_rh (T_dur) are totally destroyed, where T_dur is the duration time when ${\mathfrak{R}}$ is larger than ${{\mathfrak{R}}}_{c}$. That is, T_dur determines how long GCs are affected by the tidal force from their host galaxies. We find that T_dur cannot change our results significantly so we use a fixed value of T_dur = 0.5 Gyr.

The evolution of r_h is described as follows. During the first t_rh, the stellar evolution is the main driver making GCs expand so the r_h0 of GCs expands ∼1.3 times during the first t_rh by the stellar evolution. During the pre-core-collapse stage of GCs, we assume the 1.3 r_h0 is not changed so 1.3 r_h0 is maintained until the first t_cc, where t_cc is the core collapse timescale of ≈ 10 t_rh (e.g., Spitzer 1987; Gürkan et al. 2004). After the first t_cc, because the binary-driven post-core-collapse expansion is dominant, r_h expands again (Goodman 1984; Spitzer 1987; Baumgardt et al. 2002; Shin et al. 2013):

$$\begin{eqnarray}&&{r}_{{\rm{h}}}{(t)={r}_{{\rm{h}}0}(t/{t}_{\mathrm{cc}})}^{(2+\nu )/3},\end{eqnarray} \tag{ 9 }$$

where ν ≈ 0.1 (e.g., Goodman 1984).

The evolution of r_J is described as follows. If we assume GCs are mainly affected by the tidal force of their host DM halos, r_J(t) is

$$\begin{eqnarray}&&{r}_{{\rm{J}}}{(t)=(M/2{M}_{g})}^{1/3}{R}_{g},\end{eqnarray} \tag{ 10 }$$

where M_g is the enclosed mass of host DM halos where GCs are located, and R_g is the galactocentric radius. Because we can trace the position of GCs in their host DM halos with time, we can calculate R_g as a function of time. To calculate M_g , we use the halo concentration and the virial radius in the results by the rockstar DM halo finder, which assumes the Navarro–Frenk–White profile of DM halos (Navarro et al. 1996; Jimenez et al. 2003; Mo et al. 2010). When we calculate r_J, we assume that only DM particles contribute to the M_g so we do not consider the mass of stellar and gas components. The contribution of stellar and gas components to M_g will be included in a follow-up study (e.g., Li & Gnedin 2014; Choksi et al. 2018).

Finally, we can calculate r_h and r_J with time using Equations (9) and (10). We assume GCs that satisfy ${\mathfrak{R}}\gt {{\mathfrak{R}}}_{c}$ longer than T_dur = 0.5 Gyr are totally destroyed. The thin red line in Figure 1 shows the MF after applying tides together with the stellar evolution and two-body relaxation. Compared with the blue line, the number of low-mass GCs is decreased due to the tidal force by their host DM halos.

Some GCs are gradually moving into the center of DM halos because of dynamical friction, described as orbital decay. We assume GCs that are destroyed by orbital decay can contribute to the formation of nuclear star clusters (NSCs) (e.g., Antonini et al. 2012; Gnedin et al. 2014; Sánchez-Janssen et al. 2019, and references therein). If GCs are in the circular orbit and only the DM contributes to the velocity dispersion of GCs, the dynamical friction timescale (Binney & Tremaine 2008; Gnedin et al. 2014) is

$$\begin{eqnarray}&&{t}_{\mathrm{df}}=0.65\,\mathrm{Gyr}{\left(\displaystyle \frac{{R}_{g}}{\mathrm{kpc}}\right)}^{2}\left(\displaystyle \frac{{V}_{c}({R}_{g})}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{M(t)}{{10}^{5}\,{M}_{\odot }}\right)}^{-1}{f}_{\epsilon },\end{eqnarray} \tag{ 11 }$$

where V_c (R) is a circular velocity of the DM halo at R_g , and f = 0.5 (Gnedin et al. 2014). Here, we assume V_c ∼ σ, where σ is the velocity dispersion of the host DM halo. We calculate the new t_df whenever the host DM halos of GCs are changed due to DM halo mergers. If t_df is shorter than the timescale when GCs are staying in their host DM halos, we assume GCs are destroyed and treated as NSCs.

Finally, the thick black line in Figure 1 shows the final GCMF at z = 0. The orbital decay makes the MF slightly lower than the one without it (the thin red line), maintaining a log-normal shape with a peak mass at ∼5.5 × 10⁵ M_⊙. To compare our final GCMF with the observation, we overplot the present GCMF of M49, whose M_halo ≃ 10¹⁴ M_⊙ is similar to that of Targets 1–3 (the average M_halo ≃ 10¹⁴ M_⊙). The GCMF of M49 is converted from the GC luminosity function with the constant mass-to-light ratio of 2.69 in the g band (Jordán et al. 2007b; Willmer 2018). We find that the peak mass and the log-normal shape match well between our model and the observation. Note that, if tidal shocks that we do not consider for the GC disruption are applied, the GCMF can move to the low-mass part while leaving the log-normal shape nearly invariant (e.g., Fall & Zhang 2001; Prieto & Gnedin 2008; Shin et al. 2008; Pfeffer et al. 2018; Li & Gnedin 2019; Reina-Campos et al. 2022a).

Our fiducial parameter set is determined to reproduce the overall feature of the observations, the M_GCs–M_halo relation, GC occupancy, and f_blue, but is not tuned to match the observations exactly. Our fiducial values for free parameters are as below. GCs are formed until z_c = 1, which corresponds to the lookback time of 8 Gyr, with an estimated minimum age of GCs in the Virgo cluster (e.g., Ko et al. 2022). We assume GCs form from DM halo mergers of γ_MR = 0.1 with the initial $\mathrm{log}\eta =-3.3$. We assume that r_h0 is 3 pc, which is the median half-light radii of typical GCs in the Virgo cluster (e.g., Jordán et al. 2005) and GCs with ${{\mathfrak{R}}}_{c}\gt 0.05$ for T_dur = 0.5 Gyr are totally destroyed. The values of the fiducial parameter set are shown in the first row in Table 1.

Table 1. Parameter Setting

	zc	γMR	Initial $\mathrm{log}\eta$	rh0	${{\mathfrak{R}}}_{c}$	χ2	PKS	χ2
				pc		MGCs–Mhalo	GC Occupancy	fblue
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
Fiducial	1.0	0.1	−3.3	3.0	0.05	0.057	0.96	0.041
	0.0	0.1	−3.3	3.0	0.05	0.053	1.00	0.043
	4.0	0.1	−3.3	3.0	0.05	0.212	0.00	0.051
	1.0	0.01	−3.3	3.0	0.05	0.132	0.81	0.033
	1.0	0.3	−3.3	3.0	0.05	0.103	0.63	0.047
	1.0	0.1	−3.0	3.0	0.05	0.045	0.94	0.033
	1.0	0.1	−4.0	3.0	0.05	0.204	0.04	0.065
	1.0	0.1	−3.3	2.0	0.05	0.115	0.10	0.051
	1.0	0.1	−3.3	10.0	0.05	0.102	0.02	0.056
	1.0	0.1	−3.3	3.0	0.5	0.053	0.98	0.035
	1.0	0.1	−3.3	3.0	0.005	0.355	0.00	0.254

Note. Column (1): the lowest redshift for the GC formation. Column (2): the minimum DM halo merger mass ratio. Column (3): the initial mass fraction of the GC system and their host galaxies. Column (4): the initial half-mass radius. Column (5): the criteria of the ratio between the half-mass radius and the Jacobi radius. Column (6): the χ² value for the M_GCs–M_halo relation. Column (7): the K-S probability for GC occupancy. Column (8): the χ² value for the number fraction of blue GCs. Bold numbers show the changed values.

Download table as:  ASCIITypeset image

Figure 2 shows the M_GCs–M_halo relation of our fiducial results and the observation of the ACS Virgo cluster survey (Peng et al. 2008). We convert M_stellar of observed galaxies to M_halo, using the M_halo–M_stellar relation (Hudson et al. 2015). The observation and our fiducial results have a χ² of 0.025 and 0.057, respectively, but the overall trend is very similar to each other. The mean M_GCs/M_halo of our fiducial results is 3.5 × 10⁻⁵, which is similar to the observation of M_GCs/M_halo ≈ 10⁻⁵ (Harris et al. 2013). The gray lines are fitting lines of the observation (dashed) and our fiducial results (solid). The slopes of fitting lines of our fiducial results and the observation are 1.04 and 1.14, respectively.

Zoom In Zoom Out Reset image size

Figure 3 shows the GC occupancy (the number ratio of galaxies that contain GCs among entire galaxies) of our fiducial results compared to the observation of NGVS (Sánchez-Janssen et al. 2019, see their Table 4). In both our fiducial results and the observations, most galaxies more massive than M_stellar of 10⁹ M_⊙ contain GCs so the overall match is good (P_K-S ∼ 0.96). The number of low-mass galaxies that contain GCs is smaller than massive galaxies because GCs that form in low-mass galaxies are easily destroyed due to their low mass and GCs in low-mass galaxies can move to massive galaxies by galaxy mergers.

Zoom In Zoom Out Reset image size

The left panel in Figure 4 shows the number fraction of blue GCs in host galaxies (f_blue). We define blue GCs as [Fe/H] < −1.0 and red GCs as [Fe/H] > −1.0. For comparison, we overplot the observation of the ACS Virgo cluster survey (Peng et al. 2008). The gray-dashed line is the fitting line of the observed data points of the Virgo cluster (see Harris et al. 2015, their Equation (3)). Values of χ² of our fiducial results and the observation for f_blue are 0.041 and 0.022, respectively. Our fiducial results show higher χ² than the observation because there are many low-mass halos that contain only blue GCs (e.g., Peng et al. 2008; Harris et al. 2015). Red GCs can form in halos more massive than M_stellar of 2.6 × 10⁸ M_⊙ at z = 1 by Equation (5). In our simulations, the fraction of massive halos at z = 1 is only ∼0.01 so red GCs can form relatively less than blue GCs. However, the overall behavior is similar to what is observed (increasingly red GCs as the halo mass increases), even if the quantitative match is not very good.

Zoom In Zoom Out Reset image size

The right panel in Figure 4 shows the mean metallicity of blue and red GCs with each GC mass (M_GC) of our fiducial results and the observation of the ACS Virgo cluster survey (Peng et al. 2006; Choksi et al. 2018). In observations, the mean metallicity for blue GCs increases as M_GC increases (see the blue-dashed line), which is known as the blue tilt. Similar to the observation, our fiducial result (the solid line) also shows the blue tilt (e.g., Choksi et al. 2018; Usher et al. 2018), but the mean metallicity of blue and red GCs tends to have offsets to slightly lower and higher values of 0.1–0.3 dex.

In this section, we investigate how our model results are affected by the variation of free parameters. To try to understand the importance and impact of free parameters on our results, we systematically vary each free parameter from the fiducial values (see the first row in Table 1) while keeping the others fixed (bold numbers in Table 1 show the changed values). Free parameters that we can change in our model are the redshift cut (z_c), the merger ratio (γ_MR), the initial mass ratio of M_GCs and M_halo ($\mathrm{log}\eta$), the initial half-mass radius of GCs (r_h0), and the ratio of r_h and r_J (${\mathfrak{R}}$). Table 1 summarizes the free parameters and the altered values we consider for them. Note that the fiducial value is chosen to reproduce the overall observations, but is not finely tuned to match all observations exactly (see Section 3.1).

For comparison with our fiducial results, we use the M_GCs–M_halo relation, GC occupancy, and f_blue. Figures 5–7 show the results of each parameter set, compared to our fiducial parameter results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Redshift cut. First, we investigate how GC properties are changed with the redshift cut (z_c ), which limits the minimum GC formation epoch. The fiducial value of z_c is 1. However, young GCs exist in extragalactic systems (Glatt et al. 2008; Ko et al. 2018; Usher et al. 2019) and the mass and radii of young massive star clusters (YMCs) are similar to GCs in the local universe (Portegies Zwart et al. 2010; Longmore et al. 2014; Forbes et al. 2018). Thus, some simulation papers continue to make GCs until z_c = 0, assuming that YMCs could be one of the precursors of future GCs (e.g., Prieto & Gnedin 2008; Li & Gnedin 2014; Choksi et al. 2018; El-Badry et al. 2019; Chen & Gnedin 2022). Applying this assumption in our model, z_c = 1 is changed to z_c = 0. On the other hand, the average age of observed GCs in the Milky Way is old (10–13 Gyr) so we also test using z_c = 4.

The results are shown in the first and second panels in the first rows of Figures 5–7. In Figure 5, the height of M_GCs with z_c = 0 is 0.30 dex higher than that of the z_c = 1 case. The slope of the z_c = 0 case is 1.13 similar to that of the z_c = 1 case of 1.04. In Figure 6, the GC occupancy of the z_c = 0 case is similar to that of the z_c = 1 case (P_K-S = 1.0). Because both the z_c = 0 and z_c = 1 cases experience the epoch of the maximum number of galaxy merger frequencies at z = 3, this results in the similar trend of both the z_c = 0 and z_c = 1 cases. In Figure 7, the z_c = 0 case can make more red GCs than the z_c = 1 case because many red GCs form at low redshift. Thus, the z_c = 0 case has more galaxies that have lower f_blue than the z_c = 0 case.

The z_c = 4 case differs more strongly from the z_c = 1 case because galaxies stop making GCs before the epoch of the peak galaxy merger frequency (z = 3). It means insufficient GCs are made overall (see Figures 5 and 6) and the early formation of GCs makes galaxies have mostly blue GCs (see Figure 7 and Equation (5)). Therefore, the z_c = 4 case shows a worse match than the z_c = 0 case.

When we stop making GCs at z = 1, it matches the overall observed properties of GCs in the Virgo cluster. It implies that the GC formation in galaxy clusters after z = 1 might be rare. One possible explanation is that the cosmic star formation rate (SFR) has a peak at between z = 2 and 1 and the SFR is rapidly decreasing after z = 1 in the range of 10¹¹ M_⊙ < M_halo < 10¹³ M_⊙ (Behroozi et al. 2019). The other is that after galaxies fall into the galaxy cluster, they lose their gas by harassment and ram pressure stripping (Moore et al. 1996; Hester 2006; Boselli et al. 2019). In this case, although the galaxy mergers continue after falling into the galaxy cluster, GCs cannot form.

Merger ratio. We consider the minimum galaxy merger ratio to produce the GC system (γ_MR) is 0.1 as a fiducial value so both major mergers and minor mergers can make GCs. We use 0.01 rather than 0.1 to include the extreme minor mergers and use 0.3 for major mergers only.

The results are shown in the third and fourth panels in the first rows of Figures 5–7. In Figure 5, the γ_MR = 0.01 case makes the height of M_GCs increase 0.75 dex entirely because most galaxies undergo more frequent mergers than the γ_MR = 0.1 case. The slopes of the γ_MR = 0.01 and γ_MR = 0.1 cases are 1.04 and 0.97, respectively, so the slope stays quite similar. More frequent galaxy mergers of the γ_MR = 0.01 case make the GC occupancy higher than the γ_MR = 0.1 case (see Figure 6). When γ_MR is 0.3, the height of M_GCs and the GC occupancy slightly decrease in each galaxy (see Figures 5 and 6). This is because fewer GCs form in individual galaxies due to a smaller major merger frequency than the minor merger frequency. On the other hand, in Figure 7, we can see that f_blue is not significantly affected by γ_MR because the total number of GCs per galaxy is just determined by γ_MR while the number ratio of blue and red GCs is mainly determined by z_c .

The initial mass ratio between the GC system and the host galaxy. The initial mass ratio of the GC system and its host galaxies ($\mathrm{log}\eta$) is set to be −3.3 at the GC formation epoch as our fiducial value. We test two more initial $\mathrm{log}\eta$ of −4.0 and −3.0 (e.g., Choksi & Gnedin 2019a). These values assume a linear relation between $\mathrm{log}{M}_{\mathrm{GCs}}$ and $\mathrm{log}{M}_{\mathrm{halo}}$, but only the height is changed on the plot.

The panels in the second rows of Figures 5–7 show the results. In Figure 5, if the initial $\mathrm{log}\eta$ is high, more GCs can form in each galaxy so the height of M_GCs of the $\mathrm{log}\eta =-3$ case is 0.03 dex higher than the $\mathrm{log}\eta =-3.3$ case with slopes of 1.07 and 1.04, respectively. If the initial $\mathrm{log}\eta$ is low, insufficient GCs form so the height of M_GCs of the $\mathrm{log}\eta =-4$ case is 0.51 dex lower than the one of the $\mathrm{log}\eta =-3.3$ case. The GC occupancy is also affected by the initial $\mathrm{log}\eta$ because the $\mathrm{log}\eta =-3$ and $\mathrm{log}\eta =-4$ cases can make more and less GCs than the $\mathrm{log}\eta =-3.3$ case (see Figure 6). In Figure 7, the initial $\mathrm{log}\eta$ does not affect f_blue because the number ratio of blue and red GCs is mainly driven by z_c .

The initial half-mass radius. Next, we change the initial half-mass radius (r_h0). Our fiducial value of r_h0 is 3 pc. We use 2 and 10 pc rather than 3 pc to see the effect of r_h0.

The first and second panels in the third rows in Figures 5–7 show the results. In Figure 5, the slopes of the r_h0 = 2 pc and the r_h0 = 10 pc cases are 1.11 and 0.98, similar to the fiducial results of 1.04, but the heights of both cases are 0.22 and 0.41 dex lower than the r_h0 = 2 pc case. Both the r_h0 = 2 pc and r_h0 = 10 pc cases produce lower GC occupancy than the r_h0 = 3 pc case entirely (see Figure 6). This means too many GCs are being disrupted, but the main disruption process is different between the r_h0 = 2 pc and the r_h0 = 10 pc case. If r_h0 is 2 pc, two-body relaxation is the dominant process because t_rlx is shorter than the r_h0 = 3 pc case (see Equation (7)). But in the r_h0 = 10 pc case, the tidal force from host galaxies is the dominant process because most GCs are in the tidal regime by the criteria of ${\mathfrak{R}}\gt {{\mathfrak{R}}}_{c}$. Even though more GCs are being destroyed, they are not changing the overall M_GCs–M_halo trend. Instead, data points are simply being removed from the trend so the occupancy falls.

In Figure 7, if r_h0 is 2 pc, there are more red GCs than in the r_h0 = 10 pc case. This is because if r_h0 is 2 pc, blue GCs that form at high redshift are destroyed by a short t_rlx. In the case of red GCs, they have formed relatively recently compared to blue GCs, so red GCs can survive until z = 0 in spite of a short t_rlx. When r_h0 is 10 pc, most blue and red GCs have larger ${\mathfrak{R}}$ than ${{\mathfrak{R}}}_{c}$ so they are quickly destroyed by tides from host galaxies. Although GCs have r_h0 = 10 pc, some blue and red GCs that can have low ${\mathfrak{R}}$ might be isolated from the tides. Therefore, it makes isolated GCs have a long t_rlx so they can survive until z = 0. These results are also seen in N-body simulations of individual star clusters (Webb et al. 2014; Zonoozi et al. 2016; Park et al. 2018).

The ratio of the half-mass radius to the tidal radius. We use an ${{\mathfrak{R}}}_{c}$ of 0.05 (e.g., Gieles & Baumgardt 2008) to divide the GCs into the tidal regime (${\mathfrak{R}}\gt {{\mathfrak{R}}}_{c}$) and the isolated regime (${\mathfrak{R}}\lt {{\mathfrak{R}}}_{c}$), so GCs in the tidal regime are destroyed by the tidal force from host galaxies. To investigate how ${{\mathfrak{R}}}_{c}$ affects the disruption of GCs, we change ${{\mathfrak{R}}}_{c}$ from 0.05 to 0.5 and 0.005.

The results are shown in the third and fourth panels in the third rows in Figures 5–7. The trend of ${\mathfrak{R}}$ of each GC is increasing due to Equation (9), although there is a fluctuation of ${\mathfrak{R}}$ by r_J. That is, GCs that have ${{\mathfrak{R}}}_{c}=0.5$ have already experienced the ${{\mathfrak{R}}}_{c}=0.05$ regime. This results in similar trends of the M_GCs–M_halo relation and the GC occupancy between the ${{\mathfrak{R}}}_{c}=0.5$ and ${{\mathfrak{R}}}_{c}=0.05$ cases (see Figures 5 and 6).

In Figure 7, many galaxies have lower f_blue in the ${{\mathfrak{R}}}_{c}=0.5$ case than the ${{\mathfrak{R}}}_{c}=0.05$ case. Because the tides are weaker in the ${{\mathfrak{R}}}_{c}=0.5$ case compared to the ${{\mathfrak{R}}}_{c}=0.05$ case, more red GCs can survive in the ${{\mathfrak{R}}}_{c}=0.5$ case. However, most GCs with ${{\mathfrak{R}}}_{c}=0.005$ are destroyed because they are much more affected by a tidal force from their host galaxies. As a result, none of the results with ${{\mathfrak{R}}}_{c}=0.005$ are similar to the ${{\mathfrak{R}}}_{c}=0.05$ case. Therefore, increasing ${{\mathfrak{R}}}_{c}$ does not do much but decreasing this value is very significant for our model.

Because we use the PTM, we can trace the positions of tagged particles as GCs with time. In this section, we focus on additional observations for comparison, using GC position information. To compare our results with the observation, we add two more GC catalogs in this section: the point-source catalog in the Sloan Digital Sky Survey (SDSS) Sixth Data Release (Adelman-McCarthy et al. 2008) and the Canada–France–Hawaii Telescope Legacy Survey (CFHTLS; Hudelot et al. 2012). Lee et al. (2010) select the brightest GC candidates in a circular field with a radius of 9°, using the photometry of the point sources in the SDSS Sixth Data Release. They use the criteria for color and magnitude, 0.6 < (g − i)₀ < 1.3 and 19.5 < i₀ < 21.7 mag, with reddening correction, and the magnitude limit is i₀ < 21.7 mag. CHFTLS covers 155 deg² across four patches that comprise several fields with five filters: u*, ${g}^{{\prime} }$, ${r}^{{\prime} }$, ${i}^{{\prime} }$, and ${z}^{{\prime} }$. For a point source, the limiting magnitude in the i band is ∼24.7.

First, we investigate the size of the GC system in each galaxy. Figure 8 shows the median radii of the GC system ($\mathrm{log}{R}_{0.5}$) with the M_stellar of their host galaxies. The observed effective radii of the GC system are taken from nearby galaxy groups (Hudson & Robison 2018). To measure the median radii of the GC system, we use galaxies that contain more than 10 GCs inside 0.1R_vir (orange solid dots) or 1.0R_vir (brown open circles). Adopting galaxies that have more than 10 GCs can improve the statistics because using 10 GCs results in higher confidence in calculating $\mathrm{log}{R}_{0.5}$. When we increase the radial cut from 0.1R_vir to 1.0R_vir, $\mathrm{log}{R}_{0.5}$ slightly increases but the overall trend is almost the same. The trend of $\mathrm{log}{R}_{0.5}$ is found to increase with M_stellar and it is similar to the observations of galaxies with M_stellar > 10¹⁰ M_⊙, although there are large vertical scatters in the observation. From M_stellar = 10⁸ M_⊙ to M_stellar = 10¹⁰ M_⊙, the trend of $\mathrm{log}{R}_{0.5}$ is almost constant. Thus, the trend of $\mathrm{log}{R}_{0.5}$ shows a broken-power law with a broken point at ∼5 × 10¹⁰ M_⊙.

Zoom In Zoom Out Reset image size

We also overplot the initial $\mathrm{log}{R}_{0.5}$ with black dots in Figure 8 to show how the GC system size has changed from their initial size. In our model, it is difficult to define the initial $\mathrm{log}{R}_{0.5}$ in each galaxy because GCs are added whenever galaxies experience mergers. Thus, we just overplot the initial $\mathrm{log}{R}_{0.5}$ at the first merger in each galaxy for simplicity. We find that the initial $\mathrm{log}{R}_{0.5}$ is an extension of the current $\mathrm{log}{R}_{0.5}$.

Next, we investigate the distribution of ICGCs in galaxy clusters. Various galaxy cluster surveys have detected ICGCs and found that the number density profile of blue ICGCs is more extended than that of red ICGCs (Lee et al. 2010; Peng et al. 2011; Durrell et al. 2014; Madrid et al. 2018; Harris et al. 2020). Figure 9 shows the projected density profiles of the blue and red ICGCs of our fiducial results and the observation of the Virgo cluster (Lee et al. 2010). In the observation, ICGCs are defined by the masking method (e.g., Lee et al. 2010; Ko et al. 2017, 2018; Harris et al. 2020). ⁹ In our simulation, we define ICGCs that are actual members of the galaxy cluster, rather than members of satellite galaxies using the binding energy calculated by the rockstar DM halo finder.

Zoom In Zoom Out Reset image size

Our fiducial results show that the blue ICGCs generally have a more extended surface density profile as the observation, although there is an increasing trend of red ICGCs at the outskirt of the galaxy cluster due to red ICGCs that are coming from the recent mergers of Milky Way-size galaxies. Because blue GCs are older than red GCs, there is a high probability that blue GCs can escape their host galaxies by accretion or merger. It makes blue ICGCs have a higher mass density profile than red ICGCs.

In Figure 10, we investigate the mean metallicity (left panel) and the mean age (right panel) of GCs as a function of a clustercentric radius. For comparison, we overplot the NGVS observation as open circles (Ko et al. 2022). In the left panel of Figure 10, both our fiducial results and the observation show a decreasing mean metallicity of blue and red GCs with the clustercentric radius. Especially, the gradient and the height of the mean metallicity of blue GCs match the observation well. However, in the case of red GCs, there is an offset of the mean metallicity between our fiducial results and the observation. We revisit this issue in the discussion (Section 5). In the right panel of Figure 10, the mean age of blue GCs is decreasing with the clustercentric radius in both our fiducial results and the observation but there is no overlap. However, there is a possibility that the age of GCs is estimated as younger (1–2 Gyr) because integrated stellar populations from the observations can be affected by changing the morphology of the horizontal branch (Ko et al. 2022). In this case, there might be an overlap of blue GCs between our fiducial results and the observation but the age of red GCs in our fiducial model is still younger than the observation.

Zoom In Zoom Out Reset image size

To understand how properties of the GC system are changed by the clustercentric radius, we investigate the M_GCs–M_halo relation and the specific mass (S_M = 100M_GCs/M_stellar; Peng et al. 2008). We divide galaxies into inner (r < 1 Mpc) and outer galaxies (r > 1 Mpc).

The left panel in Figure 11 shows the M_GCs–M_halo relation of our fiducial results. The slopes of the inner and outer galaxies are 1.06 and 0.94, respectively, and the height of the inner galaxies is 1 dex higher than the outer galaxies. The right panel in Figure 11 shows the S_M of our fiducial results and the observation of the ACS Virgo cluster survey (Peng et al. 2008). The inner galaxies have higher S_M than the outer galaxies in both our fiducial results and the observation (e.g., Peng et al. 2008; Harris et al. 2013). There are insufficient high-mass and low-mass galaxies in both our simulation and the observation, respectively, so it results in a significant difference between our fiducial model and the observation at high-mass and low-mass regions of galaxies (e.g., Mistani et al. 2016; Carlsten et al. 2022). However, we can still see a U-shape of the S_M of the inner galaxies in both our fiducial results and the observation.

Zoom In Zoom Out Reset image size

Overall, Figure 11 shows that the inner galaxies have more GCs than the outer galaxies. This results from the inner galaxies generally experiencing more frequent mergers than the outer galaxies due to an early infall to galaxy clusters.