Internal R-process Abundance Spread of M15 and a Single Stellar Population Model

2.1.1. Cosmologial Simulations

2.1.2. Star Formation and Feedback

A gas cell bears star particles when the hydrogen number density of the gas is higher than 10⁴ per 1 cm³. Star formation is allowed only in cool (T < 10000 K) gas cells. This temperature condition is imposed to avoid star formation in dense gas shocked by supernovae (SNe).

We adopt a momentum injection method for the SN feedback as in Marinacci et al. (2019). Briefly, if the mass resolution of a simulation is not as good as ∼1000 M_⊙, the momentum enhancement during the Sedov–Taylor phase is taken into account (Kimm & Cen 2014; Hopkins et al. 2018; Marinacci et al. 2019). Taking ${P}_{\mathrm{SN}}=\sqrt{2{E}_{\mathrm{SN}}{M}_{\mathrm{SN}}}$ as the default momentum, we boost it by a factor of ${f}_{\mathrm{enhance}}=\sqrt{1+{M}_{\mathrm{cell}}/{\rm{\Delta }}M}$. Here the E_SN and M_SN are energy of supernova and mass of supernova ejecta, which we assume to be 10⁵¹ erg and 15.2 M_⊙, and ΔM is the mass imparted to a gas cell. There is a maximum in the boosted momentum p_terminal from each SN (Cioffi et al. 1988):

$$\begin{eqnarray}&&{p}_{\mathrm{terminal}}=4.8\times {10}^{5}{E}_{\mathrm{SN}}{\left(\displaystyle \frac{{n}_{{\rm{H}}}}{1\ {\mathrm{cm}}^{-3}}\right)}^{-1/7}f{\left(Z\right)}^{3/2}\,{\text{}}{M}_{\odot }\mathrm{km}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 1 }$$

where $f(Z)=\min [{\left(Z/{Z}_{\odot }\right)}^{-0.14},2]$.

The importance of radiation feedback has been pointed out in previous studies (Stinson et al. 2013). We implement a simple model for radiation feedback. We assume that each star particle keeps injecting energy for the first 4 Myr after formation. The energy emission rate is 2.3 × 10³⁶ erg s⁻¹ per solar mass. This rate corresponds to the number of UV photons of 8.0 × 10⁴⁶ s⁻¹ with an average energy of 〈h ν〉 = 18 eV (e.g., Fukushima et al. 2020). We turn off radiative cooling for 4 Myr in cells affected by the radiation feedback.

The kernel function for both the radiation feedback and the SNe feedback is proportional to the solid angle subtended by the cell viewed from the SN site (star particle). Namely,

$$\begin{eqnarray}&&{w}_{i}=\displaystyle \frac{{\rm{\Delta }}{{\rm{\Omega }}}_{i}}{4\pi }.\end{eqnarray} \tag{ 2 }$$

The kernel is calculated for the closest 64 cells around the star particle.

2.1.3. Identification of GCs

GCs are typically very compact and self-bound objects with an approximately constant [Fe/H]. We identify candidate star clusters as follows. First, we locate star particles in a high stellar density environment with at least 3 × 10⁵ M_⊙ of stars within 10 pc. Next, we group the stars using the k-means clustering algorithm. We obtain 19 star particle groups that are sufficiently dense and well separated from each other. Two of them are very massive (>10⁸ M_⊙) and contain star particles with various ages and metallicities. These are likely nuclear clusters in the centers of galaxies, and we eliminate them from our sample. The other 17 clusters have had short (≲10 Myr) star formation duration and have homogeneous iron abundances. The range of stellar masses of the 17 GCs is 4.0 × 10⁵ M_⊙–2.8 × 10⁶ M_⊙. Finally, we regard stars within 89 pc (300 cpc h⁻¹) from the cluster center as members of the GC. We check if the clusters survive in the final snapshot at z = 4. These star clusters likely remain bound until the present epoch, to be observed as "old" GCs. We note that, if we change the distance criterion to 5 and 20 pc, the number of the identified GC varies to 10 and 21.

Figure 1 shows the projected gas density at redshift 4. Among the simulated GCs, we select one as our target GC, which is marked by a large blue circle. Any of these GCs would give similar results and we pick up this sample as a typical example. The positions of other GCs are shown as black filled circles. The right panel of Figure 1 shows the distributions of the formation times and metallicities of the member stars of our target GC. The first star particle is formed at a cosmic age of 861.4 Myr, and star formation continues for 5.8 Myr. The mean formation time is 865.2 Myr, which corresponds to z = 6.5. The mean of iron abundance [Fe/H] is −2.59, and its standard deviation is 0.04. The total stellar mass of the target GC at our final snapshot is 6.4 × 10⁵ M_⊙. The small Fe abundance spread is consistent with that of a GC in observations (∼0.06 dex; see Carretta & Gratton 1997).

Zoom In Zoom Out Reset image size

We test our inhomogeneous mixing model by using direct numerical simulations. We take an approach similar to our previous paper (Tarumi et al. 2020), where we model an r-process event, NSM, as a point explosion. We deposit Eu to the interstellar medium (ISM) that is going to form the target GC. We follow the dispersal and dilution of Eu to the surrounding ISM. The main stellar population of the GC is formed after the r-process element enrichment.

We use 20 snapshots dumped at times before the GC formation, separated by 3–5 Myr from the adjacent output time. On each snapshot, we first determine the formation site of the target GC. To this end, we first identify the IDs of gas cells at a moment when the GC is forming. We then look for the cells back in time. In each snapshot we identify the flagged gas cells, and determine the GC formation planned site as the median of all the flagged gas cells. We then deposit Eu into the gas cells around the planned site. To investigate the relation between the distance and the resulting Eu abundance spread, we also test off-center NSM. First we take four concentric spheres centered at the planned site. On each surface, we pick up 26 points, which are (θ, ϕ) = (0, 0), (π/4, i), (π/2, i), (3π/4, i), and (π, 0) with i runs every π/4 from 0 to 7π/4 in spherical polar coordinates. In total, we take 26 × 4 = 104 points in addition to the planned site as the explosion sites. The radii of the spheres are chosen considering the time duration between the element deposition and the GC formation. Let t be the deposition time before the GC formation. For snapshots with (t < 30 Myr, 30 Myr < t < 60 Myr, 60 Myr < t), the radii are (100, 200, 300, 400)pc, (200, 400, 600, 800)pc, and (300, 600, 900, 1200)pc, respectively.

We inject a scalar component that mimics Eu into gas cells around the explosion. Eu is injected into a gas cell if it is within the snowplow radius r_sp of an injection center. The snowplow radius r_sp is calculated as follows (Mo et al. 2010). At the end of the self-similar Sedov–Taylor phase, the blast-wave radius is

$$\begin{eqnarray}&&{r}_{\mathrm{sh}}=23\times {\left(\displaystyle \frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{-19/45}\times {\left(\displaystyle \frac{E}{{10}^{51}\mathrm{erg}}\right)}^{13/45}\,\mathrm{pc},\end{eqnarray} \tag{ 3 }$$

and the velocity is

$$\begin{eqnarray}&&{v}_{\mathrm{sh}}=200\times {\left(\displaystyle \frac{n}{1{\mathrm{cm}}^{-3}}\right)}^{2/15}\times {\left(\displaystyle \frac{E}{{10}^{51}\mathrm{erg}}\right)}^{1/45}\mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 4 }$$

where n is the number density in units of cm⁻³. After that, the shell expands conserving the momentum (snowplow phase). When the shell expansion velocity becomes comparable to the turbulent velocity of the ISM (∼10 km s⁻¹), the snowplow phase ends and it blends into the ISM. The radius at the end of the snowplow phase is

$$\begin{eqnarray}&&{r}_{\mathrm{sp}}={r}_{\mathrm{sh}}\times {\left(\displaystyle \frac{{v}_{\mathrm{sh}}}{10\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{1/3}.\end{eqnarray} \tag{ 5 }$$

This approach typically distributes Eu to ∼10⁵ M_⊙ of gas cells.

The amount of injected Eu is 5 × 10⁻⁵ M_⊙, which is consistent with GW170817 (Cowperthwaite et al. 2017) and with estimates from Reticulum II observation (Ji et al. 2016). The amount of Eu assigned to each gas cell is set proportional to the volume of the gas cell. Observations suggest that the r-process elements in M15 is produced with more than one event, with different ratio of light-to-heavy elements (Otsuki et al. 2006). We assume two kinds of r-process events, lanthanide-poor and lanthanide-rich, exist as the observations suggest. The natal cloud is enriched with both of these events, and it has an abundance of [Eu/H] = −2.1. This value is motivated by the fact that the metallicity of M15 is [Fe/H] ∼ −2.4 and Milky Way's halo stars with similar metallicity have [Eu/Fe] ∼0.3. This constant "pre-enrichment," though we do not specify its origin, makes it possible to directly compare our simulation results with M15 observations. Further discussion on the origin of first-peak r-process elements and lanthanides is presented in Section 4.3.

In Figure 2, we show Eu abundance distributions for four different snapshots. The red lines in the four panels show the kernel density estimation of Eu abundance distribution for M15 stars (Worley et al. 2013). The blue histograms are the results of our numerical simulations. The top left panels show the result of one of our "best" models in which Eu is deposited 19.0 Myr before the GC formation. The explosion site is 100 pc away from the formation site. The observed Eu abundance distribution is reproduced well. Other realizations with the same model parameters with Δt = 19.0 Myr and R = 100 pc have comparable Eu abundance scatter, although with somewhat different distributions. The top right panel shows a slightly different model with Δt = 30.7 Myr and R = 100 pc. Although the shape of the distribution function is different, it still has a significant abundance spread within one system. Interestingly, this particular realization yields an apparent double-peak distribution. The low-abundance peak is at [Eu/H] = −2.1; some stars are barely enriched by the NSM. It is interesting that the apparent bimodality can be reproduced by our inhomogeneous ISM model without invoking multiple stellar formation epochs.

Zoom In Zoom Out Reset image size

The bottom two panels of Figure 2 are results for models that do not reproduce the observed distribution. In the model with Δt = 56.1 Myr and R = 0 pc (bottom left), even if an NSM occurs at the center, the ejecta dilutes with a large amount of ISM. We find a small Eu abundance spread with the standard deviation of 0.1 dex, which is comparable to observational uncertainties. Another model is with Δt = 73.7 Myr and R = 0 pc (bottom right). In this long duration model, Eu from the last NSM is almost completely mixed within the ISM, and the overall Eu abundance and its scatter are inconsistent with observation. Based on the simulation results, we argue that the observed Eu abundance can be reproduced if the r-process enrichment (NSM) occurs in the outer part of the natal cloud, and if a short period of time is available for ISM mixing.

Figure 3 shows the distribution of Eu abundance spread for Δ t = 10.3 and 30.7 Myr models. In Δ t = 10.3 Myr, R = 100 pc models, more than 80% show abundance spread comparable to or more than M15. The fraction decreases as we increase the distance, and none of the R ≳ 300 pc models show observable spread. In longer duration models (Δ t = 30.7 Myr), the inhomogeneous probability for R = 100 pc models decreases to 35%. This is because the NSM ejecta is diluted to a larger volume. In this long duration model, however, 15% of R = 300 pc models show significant Eu spread. Time evolution decreases the probability for small R models, but gives chances for relatively large R models to show Eu abundance spread.

Zoom In Zoom Out Reset image size

We can calculate the volume in which an NSM ejecta remains inhomogeneous. If a GC is formed within, it would have an Eu abundance spread σ of more than σ > 0.15 dex. The volume is ∼0.02 kpc³, and it is roughly constant for the first 50 Myr. This reflects the fact that the probability decrease for low R models and increase for high R models have compensated with each other. After ∼50 Myr has passed, the volume quickly decreases to zero due to the dilution.

In Figure 4, we show the observed distribution of M15 stars on the [Eu/H]–[Na/H] plane and the corresponding histograms. We see a marginal bimodality in [Eu/H] with two peaks at [Eu/H] ≃ −2.1 and [Eu/H] ≃ −1.7 (see also Section 3.1). We do not find correlation between Eu and Na. Based on these facts, we discuss a few scenarios for the formation and r-process enrichment of M15.

Zoom In Zoom Out Reset image size

A popular formation model of GCs assumes multiple star formation epochs. During the formation of the FG stars, some BNSs are formed. One (or more) of the BNSs merges before the second-generation (SG) star formation and enriches the cluster. Then, the SG stars have higher [Eu/H] than the FG stars. The observed Eu abundance spread can be reproduced in this scenario (Bekki & Tsujimoto 2017; Zevin et al. 2019). An important feature is that the r-process enrichment affects only the SG stars. In the bottom right panel of Figure 4, we show schematically the expected distribution function. The FG stars have abundances similar to those of halo stars, and we do not expect Eu spread among them. The FG stars are not affected by the r-process event nor by the neon-sodium chain (NeNa-chain). On the other hand, the SG stars are highly enriched with both the r-process event and the NeNa-chain, therefore appearing on the top right. For SG stars there could be a spread in Eu abundance, depending on the timescale of star formation, similarly to our model studied in Section 3. In this case, a clear bimodal distribution of [Eu/H] can be seen. Another prediction is the correlation of Eu to Na (and other light elements). It is known that GCs have Na–O anticorrelation (see, e.g., Bastian & Lardo 2018), which is a sign of hydrogen-burning between the FG and the SG star formation epochs. SG stars are likely enriched with both Na and Eu. Comparing the FG and SG stars, we expect a positive correlation between [Eu/H] and [Na/H]. The observed abundances of M15 do not show this feature.

In another scenario, an r-process abundance spread in the natal cloud is the origin of the abundance spread in M15 stars. In this case, the number of star formation epochs in M15 would significantly affect the [Na/H]–[Eu/H] relation. If all the stars are formed in a single burst, Eu and Na are independent of each other. Therefore, the distribution would be random, as in the bottom left panel (red colormap and histogram) of Figure 4. On the other hand, if M15 stars are formed in multiple star formation bursts, only the FG stars would have Eu spread. The SG star, being formed out of the mixture of stellar winds from the first-generation stars, would not show Eu abundance spread. The expected distribution is shown in the bottom middle panel (green colormap and histogram).

Carretta et al. (2009) investigate the Na–O anticorrelations of 15 GCs and discuss that stars with ${\rm{[Na/Fe]}}\gt {{\rm{[Na/Fe]}}}_{\min }+0.3$ can be safely regarded as SG stars. Taking [Na/H] = −2.09 as the criterion motivated by the lowest Na abundance in W13 stars ([Na/H] = −2.39: [Fe/H] = −2.32 and [Na/Fe] = −0.07), it is clear that the SG stars also show significant [Eu/H] spread comparable to that of the FG stars. Therefore, we conclude that the single starburst is strongly favored as the formation scenario of GC stars.

We have seen that a GC can have a significant Eu abundance spread if the last NSM occurred ≲50 Myr before the GC formation. Each NSM creates a region of ∼0.02 kpc³ that allows formation of a GC with Eu spread. Here we estimate how often this happens.

One NSM creates a region of 0.02 kpc³ for 50 Myr. Assuming the stellar mass of a galaxy to be 10⁷ M_⊙ at the cosmic age of 1 Gyr, the average star formation rate is 10⁻² M_⊙ yr⁻¹. Since a core-collapse supernova (CCSN) occurs in 100 M_⊙ of stars formed, the CCSN rate is ∼10² Myr⁻¹. Assuming the NSM rate is 10⁻³ of that of CCSN, NSM rate is ∼10^-1 Myr⁻¹. Therefore, ∼5 NSMs occur in 50 Myr, making 0.1 kpc³ of the region with significant Eu spread. In our simulations, most GCs are formed in the central 1 kpc. ⁸ Therefore, assuming the size of the star-forming region is ∼1 kpc³, the fraction of GCs with significant Eu spread is 10%. This is consistent with one system in ∼seven metal-poor GCs with [Eu/H] ≲ −2 (M68, NGC6287, NGC6293, M92, NGC6397, M15, and M30 to our knowledge; see also Cohen 2011; Roederer 2011). Note that the fraction is highly dependent on the metallicity (or equivalently, the baseline Eu abundance). A GC is more likely to have Eu abundance spread if the homogeneous Eu background is assumed to be low. The estimate on the volume of the GC forming cloud bears significant uncertainty. Considering the low metallicities of GCs, we argue that the assumption on the formation site may be valid only for early, immature galaxies with low stellar masses of ≲10⁷ M_⊙. Note that the low [Fe/H] of our target GC does not help realize Eu abundance spread, because we assume a [Eu/H] = −2.1 background that is consistent with the metallicity of M15.

We assume that NSMs are the origin of r-process elements. However, other processes such as collapsars (Siegel et al. 2019) and magnetorotational supernovae (Nishimura et al. 2015) are also proposed as the origin. Although there is evidence for the delay of the r-process in very metal-poor stars (Tarumi et al. 2021), the origin is still under debate.

For M15 stars, Otsuki et al. (2006) point out that the abundances of light r-process elements such as yttrium (Y) and zirconium (Zr) show smaller spread compared to that of Eu, as can be found in Figure 5. There, the lines show the enrichment by NSMs with two different yields. For the top two panels, we draw the model curve assuming that (i) the homogeneous Y, Zr abundances are [Y/H] = −2.7, [Zr/H] = −2.3, (ii) the homogeneous Eu abundance is [Eu/H] = −2.1, and (iii) the NSM ejecta has [Y/Eu] = −1.0, [Zr/Eu] = −0.9. The last conditions are motivated by the abundance ratio of r-II stars (Montes et al. 2007; Ji et al. 2019; Tsujimoto et al. 2020). Note that this ratio is lanthanide-rich by ∼0.5 dex compared to the solar r-process abundance. The NSM (or an r-process event) that has dominantly enriched M15 has lower [light-r/Eu] than the solar r-process ratio, which is consistent with that of r-II stars. Interestingly, recent multiwavelength observations of a kilonova event suggest a low opacity for the ejecta (Waxman et al. 2018). The estimated lanthanide fraction is also low, but with a significant uncertainty. Future observations and modeling of kilonova will help further clarify the detailed abundances of r-process elements produced in an NSM. Different origins for light and heavy r-process elements are supported as in previous works (e.g., Honda et al. 2006; Otsuki et al. 2006). Note that the number of lanthanide-rich r-process events that have shaped the Eu spread in M15 is likely one, considering the small physical scale of M15. The lanthanide fraction of the event and r-II stars in MW halo are similar, suggesting that the numbers of events that have contributed to the r-process element abundances of r-II stars are likely one. These stars may have formed in a narrow, isolated environment that is not much affected by lanthanide-poor r-process events.

Zoom In Zoom Out Reset image size