Nuclear Star Clusters in Cosmological Simulations

We use the ROCKSTAR halo finder (Behroozi et al. 2013) to identify the dark matter halos in the simulation outputs. Then we find the stellar center of the galaxy by doing an iterative centering process. To account for discreteness due to the small number of stellar particles in the center, we use kernel density estimation (KDE) where we smooth each star particle with a 3D Gaussian kernel, then sum the contributions to create a full stellar density field. We start by smoothing the star particles with a large kernel (σ = 2 kpc). We pick the location with the highest density, then recalculate the stellar density around that point using a kernel that is three times smaller. We repeat until we reach the resolution limit of the simulation, given by the cell size at the finest refinement level. The smallest kernel uses σ = L_cell,min = 3 pc. This iterative process has the effect of first picking out large-scale galactic structure, then focusing in on the locations with the highest stellar density within those larger scale dense regions. The location of the absolute highest stellar density sometimes is in a very massive young cluster in the outskirts of the galaxy, but the algorithm avoids these clusters in favor of true nuclear clusters. Using the large-scale galactic structure is essential to this process.

Next we define the plane of the galaxy. Following the formalism in Zemp et al. (2011), we use the inertia matrix to calculate the axis ratios of the galaxy. The inertia tensor can be written as

$$\begin{eqnarray}&&{I}_{{ij}}=\displaystyle \frac{{\displaystyle \sum }_{k}{M}_{k}\,{{\boldsymbol{r}}}_{k,i}{{\boldsymbol{r}}}_{k,j}}{{\displaystyle \sum }_{k}{M}_{k}},\end{eqnarray} \tag{ 1 }$$

where M_k is the mass of the k-th particle, and ${{\boldsymbol{r}}}_{k,i}$ is the i component of the position of the k-th particle. We use only the star particles to calculate the axis ratios, as we are interested in the stellar component. The eigenvalues of this tensor are a²/3, b²/3, and c²/3, where a, b, and c are the semi-principal axes, and a ≥ b ≥ c. We define the normal vector to the plane of the galaxy to be the eigenvector corresponding to the smallest eigenvalue (c²/3).

We project all the star particles onto the plane of the galaxy, then perform a two-dimensional KDE to obtain the surface density profile. For star particles within the inner 12 pc, we take the Gaussian kernel with a one-sided width of 3 pc (L_cell,min). Outside of 12 pc, we use a larger 6 pc kernel (L_cell,max) to reduce counting noise. The surface density is calculated by integrating the KDE density over an annulus to get the mass, then dividing by the area of the annulus. This KDE estimate is used only for the central 100 pc, while outside this radius we simply bin the star particles.

To find the nuclear cluster, we perform a two-component fit of the stellar density profile. Galaxies at these redshifts are expected to be mostly disky with an additional central component (nucleus or bulge), so we use the sum of an exponential disk and a Plummer sphere:

$$\begin{eqnarray*}&&{{\rm{\Sigma }}}_{\mathrm{disk}}(R)\propto \exp (-R/{a}_{d}),\\ &&{{\rm{\Sigma }}}_{\mathrm{cluster}}(R)\propto {[1+{(R/{a}_{c})}^{2}]}^{-2},\end{eqnarray*}$$

where a_d and a_c are the scale radii for the two components (Dejonghe 1987). We define the size of the nuclear region, ${R}_{\mathrm{NSC}}$, to be the radius where ${{\rm{\Sigma }}}_{\mathrm{cluster}}({R}_{\mathrm{NSC}})={{\rm{\Sigma }}}_{\mathrm{disk}}({R}_{\mathrm{NSC}})$. Figure 1 shows an example of this fitting process.

Zoom In Zoom Out Reset image size

To estimate the uncertainty of our determination of ${R}_{\mathrm{NSC}}$ after marginalizing over the other structural parameters, we construct a multivariate Gaussian distribution with dimensions given by the covariance of the parameters of the decomposition (normalization and scale radius of both components). We then sample from this multivariate Gaussian and recalculate the radius of the nuclear region for each realization. We take the range that encloses 68.3% of these radii to be the error of ${R}_{\mathrm{NSC}}$.

The mass of the nucleus then is the sum of the star particles within this spherical region, ${M}_{\mathrm{NSC}}\equiv {\sum }_{k}{M}_{k}({r}_{k}\lt {R}_{\mathrm{NSC}})$. We do not subtract off the disk component within this region. We also calculate the half-mass–radius of the nuclear cluster, defined as

$$\begin{eqnarray}&&M({R}_{{\rm{h}}})\equiv \displaystyle \frac{1}{2}M({R}_{\mathrm{NSC}}),\\ &&M(R)={\displaystyle \int }_{0}^{R}{{\rm{\Sigma }}}_{\mathrm{KDE}}(R)\,2\pi {RdR},\end{eqnarray} \tag{ 2 }$$

where Σ_KDE(R) is the KDE surface mass density as described above. This half-mass–radius is calculated in projection, unlike the mass, for a more direct comparison with observations. The errors on the mass and half-mass–radius are estimated by perturbing ${R}_{\mathrm{NSC}}$ according to its error.

Figure 2 shows the masses and sizes of the galactic nuclei in the simulations, along with observed NSCs in nearby late-type galaxies of intermediate and low mass (${M}_{\star }$ < 10¹¹) from Georgiev et al. (2016). We excluded one nuclear region with mass above 10⁹${M}_{\odot }$, in the SFE50-3SNR run that had insufficiently strong feedback to suppress the formation of a very dense bulge. Simulated nuclear regions with half-mass radii near 3 pc are due to a single star particle dominating the central mass distribution. For these regions, where a single particle contributes more than half of ${M}_{\mathrm{NSC}}$, we take the ${R}_{{\rm{h}}}$ given by Equation (2) to be an upper limit because they just reflect the size of the smoothing kernel. At low masses, there are many of these galaxies where a single star particle dominates. These unresolved clusters are consistent with many of the compact clusters in the Georgiev et al. (2016) sample.

Zoom In Zoom Out Reset image size

However, the nuclei that are resolved in the simulations tend to have ${R}_{{\rm{h}}}$ larger than observed clusters of the same mass. This overestimate is particularly noticeable at low masses and can be an order of magnitude or more. At high masses the discrepancy is a factor of several.

As discussed above, we select nuclear regions by means of a two-component fitting procedure. While this procedure generally provided sizes larger than the observed radii of the nearby NSCs, any more complicated model fit would be unwarranted due to the irregular structure of our galaxies. Our fit also appears to be robust and nuclear regions clearly defined. The central NSC density exceeds the extrapolated central density of the disk component typically by a factor of 10 or greater.

However, the physical interpretation of these components is not always clear. Strong feedback gives our galaxies an irregular structure at the redshift z ≈ 1.5 of our last output, making the identification of a center unclear in some cases. Some galaxies had a roughly constant density in the center, with one massive star particle that was chosen by our algorithm to be the center. Others had a more extended central component without a strongly peaked center, where the middle was chosen to be the center. Even in these situations, however, a central component is present and can represent a nuclear region, even if it is too large to be a typical NSC. Visually inspecting the profiles indicated that the returned sizes were an acceptable description of this central component. These results indicate that higher resolution or improved stellar physics are required to properly model NSCs.

Interestingly, we find that although NSC sizes are generally overestimated in the simulations, the overall galaxy sizes match relatively well with the observed sample of van der Wel et al. (2014). For late-type galaxies with 9.0 < log ${M}_{\star }$ < 9.5 at redshift 1.5 < z < 2, the median effective radius is 2.1 kpc, compared to 2.8 kpc for our simulated galaxies in the same range of mass and redshift. The distribution of sizes in the observed sample (68% lie between 1.2 and 3.7 kpc) is broader than in the simulations, which range from 2.5 to 3.3 kpc. The sizes of our simulated galaxies are thus statistically consistent with the available observations.

With the nuclear region defined, we can assess the spatial and kinematic structure of the collection of stars in this region. Using the inertia tensor (Equation (1)), we calculate the axis ratios of the nuclear regions.

Figure 3 shows the axis ratios as a function of the NSC mass. Nuclear regions that have less than 10 star particles are not shown, as they have too few particles to reliably determine the axis ratios. Our nuclear regions are mildly triaxial, but all of them are less flattened than the Milky Way NSC (Feldmeier-Krause et al. 2017b).

Zoom In Zoom Out Reset image size

To investigate whether the flattening is caused by rotation, we calculate the rotation velocity and the residual velocity dispersion of the NSCs. We define the plane of the galaxy in the same way as we did above, based on the inertia tensor. We found that this definition works better than using the angular momentum vector of the nuclear regions, because the latter is relatively small and does not determine the structure of the NSCs. We find that the direction of the nuclear angular momentum is misaligned with the normal to the plane of the galaxy by between 20° and 90°, and has no dependence on cluster mass.

The rotation velocity is calculated as the mass-weighted tangential velocity component of the star particles within the nuclear region. Figure 4 shows the ratio of rotational velocity to 3D velocity dispersion as a function of ellipticity, defined as ≡ 1–c/a. The black line shows the expected values for a system with an isotropic velocity dispersion that is flattened only by rotation. This relation can be approximated as

$$\begin{eqnarray*}&&\displaystyle \frac{{v}_{\mathrm{rot}}}{{\sigma }_{3{\rm{D}}}}\approx \sqrt{\displaystyle \frac{\epsilon }{1-\epsilon }}\end{eqnarray*}$$

(e.g., Mo et al. 2010). The simulated nuclei are predominately below this line, indicating that their nonsphericity is an intrinsic shape, not due to rotation. The accretion of particles composing the NSCs must have proceeded anisotropically. Indeed, we find that for nearly all clusters the radial velocity dispersion is larger than the tangential dispersion by up to a factor of two.

Zoom In Zoom Out Reset image size

As described above, our suite of simulations tracks enrichment from SNe Ia and SNe II, but not AGB winds. We calculate the mass-weighted total metallicity as the sum of all metals divided by the mass of all star particles in the nuclear regions.

To calculate the spread in metallicity, we consider two sources of variance: the internal spread within a star cluster particle due to self-enrichment (recorded in the simulation run time as the cluster is forming), and the dispersion of final metallicity among star particles.

To calculate the internal spread σ_Z,k² within the k-th particle, we transform the spread in overall metallicity from SN II. The details of how this is derived from the simulation outputs is described in H. Li & O. Gnedin (2018, in preparation), but we summarize the key points here. The star particles in the simulation contain a variable ${M}^{{ZZ}}=\sum {m}_{i}{Z}_{i,\mathrm{II}}^{2}$ that can be used to calculate the metallicity variance of a single star particle without needing its full accretion history. Here we are summing over all accreted mass elements m_i of a specific cluster particle as it is forming at time steps i. The final particle mass is ${M}_{k}=\sum {m}_{i}$. We only consider the spread in metallicity from SNe II, since the delay time for SN Ia is significantly longer than the formation timescale of star clusters and therefore enrichment from SNe Ia will not be relevant for this spread. The variance of the metallicity of a single star particle k is defined as:

$$\begin{eqnarray}{\sigma }_{Z,k}^{2} & = & \displaystyle \frac{{\displaystyle \sum }_{i}{m}_{i}{Z}_{i,\mathrm{II}}^{2}}{{\displaystyle \sum }_{i}{m}_{i}}-{\left(\displaystyle \frac{{\displaystyle \sum }_{i}{m}_{i}{Z}_{i,\mathrm{II}}}{{\displaystyle \sum }_{i}{m}_{i}}\right)}^{2}\\ & = & \displaystyle \frac{{M}_{k}^{{ZZ}}}{{M}_{k}}-{Z}_{k,\mathrm{II}}^{2}.\end{eqnarray} \tag{ 3 }$$

We then combine this with the dispersion in metallicity among cluster particles to obtain the total metallicity spread:

$$\begin{eqnarray}&&{\rm{\Delta }}{Z}^{2}=\displaystyle \frac{\sum {M}_{k}\,{\sigma }_{Z,k}^{2}}{\sum {M}_{k}}+\displaystyle \frac{\sum {M}_{k}{({Z}_{k}-\overline{Z})}^{2}}{\sum {M}_{k}},\end{eqnarray} \tag{ 4 }$$

where $\overline{Z}$ is the mass-weighted average metallicity of the nuclear region:

$$\begin{eqnarray}&&\overline{Z}\equiv \displaystyle \frac{\sum {M}_{k}{Z}_{k}}{\sum {M}_{k}}.\end{eqnarray} \tag{ 5 }$$

Figure 5 shows the metallicities of the nuclear regions plotted against their mass. We represent the metallicity spread with errorbars that cover the interval $\mathrm{log}(\overline{Z}-{\rm{\Delta }}Z)$ to $\mathrm{log}(\overline{Z}+{\rm{\Delta }}Z)$. Our nuclear regions have a wide spread of metallicities in the range of −2 < [Fe/H] < −0.25. Other then the most massive nuclear regions typically having higher metallicity, there is no correlation between metallicity and cluster mass.

Zoom In Zoom Out Reset image size

To turn the total metallicity into abundances for individual elements, we use computed supernova yields in the literature. We take the SN II yields of Nomoto et al. (2006) and integrate them over the Kroupa (2001) IMF from 10 to 40 ${M}_{\odot }$, using a 50% hypernova fraction for stellar masses where hypernova yields are available. Since the Nomoto et al. (2006) yields are provided at several fixed progenitor metallicities Z = 0, 0.001, 0.004, and 0.02, we interpolate between these models to get the yield at any metallicity. We use the W7 model of Iwamoto et al. (1999) for SN Ia yields at all metallicities.

To get the elemental abundances of a given star particle, we use the SN yields and scale them to produce the appropriate mass of metals, as set by the Z_Ia and Z_II variables in the simulation.

The two sets of published SN yields depend on the metallicity of the supernova progenitor star. Since this information is not recorded in the simulation output, we are using the metallicity of the star particle itself as the metallicity of the progenitor. This may overestimate the progenitor metallicity. However, as we discuss later, we use elements whose yields do not vary strongly with metallicity of the progenitor, reducing the effect on our results.

By using these yields, we calculate a mass-weighted value of [Fe/H] for both the star particles in the nuclear regions and the stars in the whole galaxy. The average contribution from all star cluster particles is

$$\begin{eqnarray}\overline{[\mathrm{Fe}/{\rm{H}}]} & = & \mathrm{log}\left[\displaystyle \frac{{\displaystyle \sum }_{k}{M}_{k}({Z}_{k,\mathrm{Ia}}{f}_{\mathrm{Fe},\mathrm{Ia}}+{Z}_{k,\mathrm{II}}{f}_{\mathrm{Fe},\mathrm{II}})}{{\displaystyle \sum }_{k}{M}_{k}\,{X}_{k}}\right]\\ & & \ -\ \mathrm{log}\left[\displaystyle \frac{{Z}_{\odot }{f}_{\mathrm{Fe}\odot }}{{X}_{\odot }}\right],\end{eqnarray} \tag{ 6 }$$

where f_Fe is the mass fraction of the metals that is Fe. The numerator in each fraction is the mass of iron, while the denominator is the mass of hydrogen. For the Sun, we use the solar abundances of Grevesse & Sauval (1998), which are Z_⊙ = 0.0169 and f_Fe = 0.0757.

We then transform the metallicity spread (Equation (4)) into spread in [Fe/H] by using:

$$\begin{eqnarray}&&{\rm{\Delta }}{[\mathrm{Fe}/{\rm{H}}]}^{2}={\rm{\Delta }}{Z}^{2}{\left(\displaystyle \frac{d[\mathrm{Fe}/{\rm{H}}]}{{{dZ}}_{\mathrm{II}}}\right)}^{2}.\end{eqnarray} \tag{ 7 }$$

We calculate this derivative numerically using interpolation of the published SN yields. While this equation is written using [Fe/H], it will also be used to calculate spreads in other elemental abundances.

We calculate the mean and total variance in [Fe/H] for both the nuclear regions and the galaxy as a whole. The results of these calculations are shown in Figure 6, along with an empirical mass-metallicity relation (MMR) for dwarf galaxies from Kirby et al. (2013), shifted down by 0.36 dex to account for evolution in the MMR from z = 0 to z = 1.5, based on the interpolation of the results of Mannucci et al. (2009) for Lyman-break galaxies at z ≈ 3.

Zoom In Zoom Out Reset image size

The model galaxies lie near the expected z = 1.5 relation, without any fitting, which provides evidence that our simulations are modeling the chemical enrichment of galaxies correctly. The nuclei are systematically more metal-rich. This indicates that our galaxies have a metallicity gradient, consistent with observations (e.g., Moustakas et al. 2010). The predicted iron abundances of NSCs span a wide range −2.5 < [Fe/H] < −0.5, uncorrelated with their mass.

To examine the effect of an extended star formation history, we first determine the time it took to form the bulk of the stellar mass in each nuclear region. To minimize the effect of few particles born much sooner or later than the bulk of the particles, we pick the smallest time interval in which 90% of the mass was formed. This interval is calculated using only the birth time of the particles, not including the time each star particle took to form. To account for this additional time, we add the maximum duration of star particle formation within the cluster ($2\,\mathrm{Myr}\lesssim {\tau }_{\mathrm{dur}}\lesssim 10\,\mathrm{Myr};$ Li et al. 2017) to the initial time interval.

Figure 7 shows that the iron spread in the nuclear regions does not correlate with the length of star formation history. In fact, several NSCs with the longest assembly time have relatively small spread, Δ[Fe/H] < 0.15, most of which is due to different metallicities of its constituent particles, rather than the intrinsic spread within each particle. The largest spread is seen in NSCs that took only 2–3 Myr to form, where all the spread is internal. These regions tend to be low metallicity ([Fe/H] < −1.8) and consist of only one or few star particles.

Zoom In Zoom Out Reset image size

The largest predicted iron spread is comparable with the spread in three massive globular clusters in the Milky Way, where it is clearly detected: ω Cen, M54, and Terzan 5. To make a fair comparison to our models, we define the observed spread as the rms scatter of the [Fe/H] content of the stars in each cluster. This value is calculated for M54 by Carretta et al. (2010a) and for Terzan 5 by Massari et al. (2014), while we calculate it for ω Cen using the data from Marino et al. (2011) In Figure 7, we do not show the observed age spreads, only the metallicity spreads. The age spreads are uncertain, with inconsistencies between different studies of ω Cen (Joo & Lee 2013; Villanova et al. 2014; Tailo et al. 2016). In our simulations, we are able to detect both age and metallicity spreads significantly below the observational limits, which remain relatively large. However, the simulated assembly timescales are within the inferred age spread of stellar populations in ω Cen of ∼500 Myr (Tailo et al. 2016), while being shorter than the age spreads in M54 and Terzan 5 (Siegel et al. 2007; Ferraro et al. 2016). As our simulations only reach z ≈ 1.5 (t ≈ 4.3 Gyr), we are not able to model the large age spreads present in these two clusters. These results lend support to the hypothesis that stripped NSCs could become progenitors of the anomalous globular clusters. There are still quantitative differences that need to be explored in future work: the nuclei with the largest iron spread all have lower metallicity than the three observed anomalous clusters and lie in the lower mass range below 10⁶${M}_{\odot }$.

When determining which elements can be modeled reliably by our simulations, we consider several aspects.

• (i)  
  Our simulations only model supernovae, not AGB stars or other sources of chemical enrichment, so we select elements without strong AGB contributions.
• (ii)  
  As discussed above, we are using the metallicity of the star particle itself as the metallicity of the supernova progenitor star. To reduce the impact on our results, we prefer elements whose yields do not vary strongly with metallicity of the progenitor.
• (iii)  
  Since the two sets of yields generally give different results, we prefer the elements with low spread between the different model sets.
• (iv)  
  We select the elements whose abundances match observations, whether in direct observation of SN ejecta, observations of low metallicity stars, or chemical evolution models.

Figure 8 shows how our two yield sets (Woosley & Weaver 1995; Nomoto et al. 2006) vary with metallicity, to address the second and third points. Table 1 describes in more detail several elements that are relevant for globular clusters and NSCs, the production source of those elements, and comments on the reliability of the yields. This table includes only elements that we deemed relevant for our analysis. We do not discuss elements such as C and N, which have significant AGB contributions, or other elements that are not strongly involved in the phenomenon of multiple populations.

Zoom In Zoom Out Reset image size

To evaluate the reliability of our use of supernovae yields to extract elemental abundances for these light elements, we surveyed the literature for studies examining this issue. Romano et al. (2010) and Andrews et al. (2017) put different sets of yields through the same chemical evolution model, then compared the result to observations, showing how different model choices affect the resulting abundances. These papers show both the scatter among different models and the differences between models and data, which are important for assessing reliability of the yields. Andrews et al. (2017) and Pignatari et al. (2016) describe the production sites of the various elements. Wiersma et al. (2009) implement the yields into cosmological smoothed-particle-hydrodynamics simulations and show in their Figure A4 the variations caused by different yields.

The synopsis of these studies is that oxygen is modeled reliably. It matches observations, has a weak metallicity dependence, and has a reasonably small spread between different yield sets. Magnesium is similar to oxygen. Aluminum and sodium have much stronger metallicity dependences and wider spreads, making their predicted abundances less certain.

To test whether our nuclear regions may have properties similar to the observed multiple populations in globular clusters, we examine their detailed elemental abundances. First, we compare the overall normalization of the abundances as a function of metallicity. Figures 9 and 10 show [O/Fe] and [Mg/Fe] for our nuclear regions and compare them to the compilation of 202 red giants from 17 globular clusters presented in Carretta et al. (2009a). Because of the limited sample size of the available measurements, we combine the available data from all clusters and plot the regions enclosing 50% and 90% of all stars. We compare this total range of observed abundances with our model predictions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each plot, a solid line shows the track the objects would follow if their metallicity came entirely from SNe II. Points below this line indicate some enrichment from Type Ia supernovae, which contribute Fe, but not O or Mg.

Additionally, we examine the elemental anticorrelations (Na–O and Al–Mg) that are key features of the multiple population phenomena in globular clusters. We compute the [Na/Fe] and [O/Fe] ratios for each star particle. Since the underlying metallicity dispersion is responsible for the elemental dispersion of a star particle, the elemental dispersions will be correlated. A given metallicity spread will cause a spread along a certain direction in the [Na/Fe]–[O/Fe] plane. We smooth each star particle by a Gaussian whose width along this direction will reproduce the elemental spreads (see Equation (7)). For visual clarity, a lower threshold on this width was chosen to be 0.005 dex, which was also used as the width perpendicular to the direction of the spread. We mass-weight the contributions of all star particles in a cluster, then plot contours that enclose 50% and 90% of the mass for each cluster. Figure 11 shows the results of this for the Na–O anticorrelation, while Figure 12 shows the Mg–Al anticorrelation. Black lines show the tracks that objects would follow if their metallicity came from a given percentage of SNe II, for metallicities −3 < [Fe/H] < 0. Nuclear regions consisting of a single particle will exhibit spreads along these lines, while regions with multiple particles can have spreads in any direction due to potential variations between particles. It is apparent that nearly all of our simulated nuclear regions do not have a significant spread in any elemental abundance, unlike stars in globular clusters. It is also apparent that the yield lines are not capable of producing the full range of elemental variations, especially for Al, even though they include metallicities up to solar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen in Equation (7), elements whose yields change more with metallicity will have larger spreads. Na and Al yields change more than O and Mg (see Figure 8), explaining their larger abundance spreads in these figures.

As shown in Section 3.5, we are unable to reproduce the full spread in light element abundances that are present in globular clusters. For [O/Fe] and [Mg/Fe], the one-sided spread is less than 0.02 dex for all clusters, while for [Na/Fe] and [Al/Fe] the spread is 0.1 dex at most. The observed one-sided spreads are significantly larger: about 0.2 dex for [O/Fe] and [Na/Fe], 0.1 dex for [Mg/Fe], and 0.3 dex for [Al/Fe]. The higher spread in Na and Al in our clusters is due to the stronger variation in the yields of these elements with metallicity (Figure 8). Even with this large variation, even the yields themselves are not able to produce as much Al as is present in enriched second-generation globular cluster stars.

The spread we calculate includes both the internal spread within a star particle due to a metallicity variation in the accreted gas during its formation, as well as the spread among different star particles within a given nuclear region. Even with both of these sources, the elemental spreads are too small. Our simulations would be able to reproduce the distinct stellar populations by having multiple star particles with differently enriched material, but we do not see this. Many nuclear regions consist of only one star particle, but those that have multiple do not show the necessary variations. In addition, a substantial age spread cannot be responsible for the elemental spreads, as none of our nuclear regions reach the required light element spread, regardless of their formation timescale.

As described in Section 2.1, the sizes of our nuclear regions are an acceptable description of the central component of the central density component, even though they are typically larger than observed NSCs. This indicates that the particles that are included within each NSC are a reasonable selection. The main conclusion of our calculation is that the extended star formation histories, without internal pollution, are insufficient to reproduce the observed abundance trends in globular clusters. This conclusion would only be strengthened if our selection of NSC particles is an overestimate. If we reduce the sizes and masses of model NSCs, their abundance spread would only decrease.

This indicates that our simulations do not contain the enrichment sources that are the true cause of the spreads. As supernovae are the only enrichment source we include, these results provide support to the idea that the polluter responsible for globular cluster abundances must lie within the clusters during their formation epoch. A number of models have been proposed to explain these light element abundance variations. Most scenarios call for multiple generations of star formation within the globular cluster's formation epoch, with enriched material from a first generation of stars seeding a second metal enriched generation. Possible versions suggested in the literature include supermassive stars (Denissenkov & Hartwick 2014), rotating massive stars (Decressin et al. 2007; Krause et al. 2013), AGB stars (Cottrell & da Costa 1981; Renzini et al. 2015), and interacting binaries (de Mink et al. 2009). However, all of these current scenarios have some shortcomings in matching the full set of available observations, as discussed by Bastian & Lardo (2018). These models fail to produce correct elemental abundances of the globular clusters, both for the elements listed above and helium. Additionally, to reproduce the large fraction of enriched stars in a globular cluster, most models require huge amounts of mass loss (of only first generation stars), which does not appear to be realistic (Larsen et al. 2012; Cabrera-Ziri et al. 2015).

One important aspect of our simulations is that they only reach z ∼ 1.5. As later star formation would make the masses, metallicity spreads, and age spreads increase, the quantities we report for those properties are best interpreted as lower limits. This complicates the comparison of our NSCs to observed NSCs and GCs, which are all at z ≈ 0. For example, the later bursts of star formation in Terzan 5 and M54 would not be present in our simulation.

With this caveat in mind, we will discuss the predicted iron spread in model NSCs. The majority of clusters have a spread of less than 0.1 dex, which is roughly the limit of the observational sensitivity (e.g., Gratton et al. 2004; Willman & Strader 2012). However, there are some clusters with [Fe/H] spreads larger than this that could be the progenitors of the iron-complex GCs. These model NSCs are very metal poor, with the metallicity distribution dominated by internal dispersion. Our simulations also include two nuclear regions that have a longer formation timescale and a [Fe/H] spread dominated by variation among the star particles. These regions may be analogous to clusters like ω Cen or M54, which may share a similar origin (Carretta et al. 2010b). Our objects formed on a timescale consistent with that of ω Cen, which appears to have an age spread of at least 500 Myr (Tailo et al. 2016). The metallicity distribution of the object with the largest spread between particles has a central peak at a metallicity consistent with the peak of the metallicity distribution of ω Cen and M54, but the tail extends to lower metallicities rather than higher. Our objects have lower mass than the observed clusters as well. These quantitative differences need to be explored in future work, but our results lend support to the hypothesis that objects like ω Cen could be made through tidal stripping of NSCs.

If this model is correct, it has implications for the location of the polluter responsible for the light element spreads of the massive GCs. If the NSC formed by in situ formation (e.g., Hartmann et al. 2011), the polluter must be present for each new burst of star formation. However, if the NSC formed by the inspiral of globular clusters (e.g., Antonini et al. 2012; Gnedin et al. 2014; Tsatsi et al. 2017), the presence of already existing light element spreads in each cluster would provide a natural explanation for the fact that light element spreads are found within each peak in the metallicity distribution of the iron-complex GCs. No polluter would be needed in the central region of the galaxy. Future observations may be able to test these hypotheses by looking for light element abundance spreads in NSCs. If not detected, both the model of massive GC formation discussed in this paper and the GC inspiral model of NSC formation would need to be discarded. If light element abundance variations are present in NSCs, these models would remain plausible.

The amount of dark matter contained within the extent of NSCs is important for their ability to represent progenitors of massive globular clusters. Even the largest globular clusters can have only a dynamically subdominant amount of dark matter within their optical radii (e.g., Conroy et al. 2011; Ibata et al. 2013). Since NSCs form at the bottom of the dark matter potential well, they could retain some of the dark matter even if the rest of the galaxy is tidally stripped.

We find that only 1 in every 11 NSCs contains any dark matter particles within ${R}_{\mathrm{NSC}}$. This does not necessarily imply that no dark matter is present. The mass and spatial resolution for dark matter is relatively low in our simulations and may prevent a reliable estimate of the dark matter content within the compact clusters. The dark matter particle mass is about 10⁶${M}_{\odot }$, so that even one particle can be more massive than the whole NSC. Also, the spatial force resolution for dark matter is kept separate from the gas and star particle resolution, at 50–100 pc, so that the discreteness of massive dark matter particles does not affect the dynamics of gas and stars. This resolution makes it difficult to put strict constraints on the amount of dark matter present, but we can still conclude that dark matter is dynamically subdominant within the model nuclear regions. If these central regions turned into an object resembling a massive globular cluster after being tidally stripped by the host galaxy, the remaining amount of dark matter would likely not violate observational constraints.