Stellar collisions in globular clusters: Constraints on the initial mass function of the first generation of stars

Globular clusters display an anticorrelation between the fraction of the first generation of stars ($N({\rm G1})/N({\rm tot})$) and the slope of the present-day mass function of the clusters ($\alpha_{pd}$), which is particularly significant for massive clusters. In the framework of the binary-mediated collision scenario for the formation of the second-generation stars in globular clusters, we test the effect of a varying stellar initial mass function (IMF) of the G1 stars on the $(N({\rm G1})/N({\rm tot}))-\alpha_{pd}$ anticorrelation. We use a simple collision model that has only two input parameters, the shape of the IMF of G1 stars and the fraction of G1 stars that coalesce to form second-generation stars. We show that a variable efficiency of the collision process is necessary in order to explain the $(N({\rm G1})/N({\rm tot}))-\alpha_{pd}$ anticorrelation; however, the scatter in the anticorrelation can only be explained by variations in the IMF, and in particular by variations in the slope in the mass interval $\approx$ (0.1-0.5) M$_{\odot}$. Our results indicate that in order to explain the scatter in the $(N({\rm G1})/N({\rm tot}))-\alpha_{pd}$ relation, it is necessary to invoke variations in the slope in this mass range between $\approx -0.9$ and $\approx -1.9$. Interpreted in terms of a Kroupa-like broken power law, this translates into variations in the mean mass of between $\approx 0.2$ and $0.55$ M$_{\odot}$. This level of variation is consistent with what is observed for young stellar clusters in the Milky Way and may reflect variations in the physical conditions of the globular cluster progenitor clouds at the time the G1 population formed or may indicate the occurrence of collisions between protostellar embryos before stars settle on the main sequence.

The origin of the multiple stellar populations in GCs is highly debated. A number of scenarios have been proposed to explain the existence of two or more stellar populations. The oldest and most popular of these is the asymptotic giant branch (AGB) self-enrichment scenario, which is based on the formation of second-generation stars (hereafter G2 stars) from the ejecta of the first-generation (G1) AGB stars after gas has been cleared from the cluster by feedback from massive stars (Cottrell & Da Costa 1981;D'Ercole et al. 2016;Bekki 2017). Since the pure AGB scenario fails to reproduce the observed Na-O anticorrelation, a modified version of the model was introduced whereby pristine gas is accreted from outside the cluster and mixes with the AGB ejecta to form the G2 stars (Ventura & D'Antona 2009;Calura et al. 2019;Yaghoobi et al. 2022). Other variants of the model relied on the ejecta of fast rotating massive stars (FRMSs) instead of AGB stars (Decressin et 2007;Krause et al. 2013) or on the rapid efficient cooling of stellar wind material that A&A proofs: manuscript no. 43556corr_final2 falls back to the center of the cluster and mixes with leftover gas with a G1-like chemical composition (Wünsch et al. 2017). Other models invoked a dichotomy in GC mass whereby the efficiency of the AGB ejecta retention in the cluster would depend on the cluster mass (Valcarce & Catelan 2011). The AGB-and FRMS-based scenarios suffer from the same types of shortcomings. It is unclear whether gas ejecta of AGB stars (or FRMSs) can cool and form the second generation of stars due to uninterrupted heating by G1 stars, particularly from X-ray binaries, that heat the gas and prevent it from fragmenting (Conroy & Spergel 2011). Another issue is the so-called mass-budget problem. This problem reflects the inability of G1 stars, under the assumption of a Milky Way-like IMF, to produce enough chemically enriched material for the formation of significant populations of G2 stars. Furthermore, and owing to the GCs orbits in the galaxy, it is not clear whether accretion of external gas would be efficient (Conroy & Spergel 2011). Khalaj & Baumgardt (2015) investigated under which conditions the loss of G1 stars due to gas expulsion from the protocluster cloud can help reproduce the fraction of G2 stars and some of the clusters properties, such as their radial profiles and mass-half mass-radius relation. They find that reproducing the observations requires rapid gas expulsion timescales ( 10 5 yr) and a large number of supernova explosions, three to six times more than what would be found in a Galactic field-like IMF.
Another physical process that can lead to the formation of G2 stars over a short timescale, prior to gas expulsion from the clusters by the first supernovae, is stellar collisions (Sills et al. 2002;Sills & Glebbeek 2010;Jiang et al. 2014;Wang et al. 2020). Collisions can still occur at later times, albeit at a lower rate (Kravtsov & Calderón 2021;Kravtsov et al. 2022). Collisions not only set the relative fractions of the G1 and G2 stars, but they also alter the shape of the IMF of the G1 stars (Dib et al. 2007;Kravtsov et al. 2022). While more work is still needed in order to understand whether stellar collisions, particularly of low-mass stars, can lead to the chemical anticorrelations that are observed for GCs, Kravtsov et al. (2022) present observational evidence that supports the scenario in which stellar collisions of G1 stars in the mass range ≈ (0.1-0.5) M ⊙ can lead to the formation of G2 stars in the mass range (0.5-0.9) M ⊙ . In particular, they presented evidence for the existence of an anticorrelation between the fraction of G1 stars (N(G1)/N(tot)) and the slope of the present-day mass function of GCs in the stellar mass range (0.2-0.8) M ⊙ (α pd ). They also find that the fraction of G1 stars is anticorrelated with the present-day encounter rate measured in these clusters. The present-day encounter rates of the clusters are expected to scale with their primordial values when the clusters were denser (e.g., Maccarone & Peacock 2011). Using a simple collision model with a parametrized collision efficiency and a Kroupa-like IMF for the G1 stars (Kroupa 2001), Kravtsov et al. (2022) showed that it is possible to reproduce the (N(G1)/N(tot)) − α pd anticorrelation both in terms of slope and absolute values. Nonetheless, the (N(G1)/N(tot)) − α pd anticorrelation shows a level of scatter that cannot be simply explained by a fixed IMF nor by observational uncertainties. In this paper, and in the framework of this collision-based scenario, we focus on the role of the IMF in explaining both the trend and scatter in the (N(G1)/N(tot)) − α pd relation. In particular, we explore how variations in the shape of the IMF can impact the shape of the (N(G1)/N(tot)) − α pd anticorrelation. In Sects. 2 and 3 we briefly recall the observational data and the main elements of the model that is used to reproduce them, respectively. The results are presented in Sect. 4, and in Sect. 5 we conclude.

Data
We used the sample of multiple stellar populations in Galactic GCs obtained by Milone et al. (2017). Using uniform multiband Hubble Space Telescope (HST) photometry (Sarajedini et al. 2007;Piotto et al. 2015) in the central parts of a sample of 57 GCs, Milone et al. (2017) isolated red giant branch (RGB) stars and measured the fractions of the G1 population to the total number of RGB stars (N(G1)/N(tot)). These fractions were measured for a total of 54 GCs. For the slope of the present-day mass function, we used the compilation of Ebrahimi et al. (2020), which includes 32 GCs and for which the slope was measured by fitting a power law in the stellar mass range (0.2-0.8) M ⊙ . The cluster NGC 6584 was excluded from our sample since the data for this GC are incomplete. The data for four additional GCs were taken from Sollima & Baumgardt (2017), namely NGC 4833, NGC 6205, NGC 6397, and NGC 6656, and the measurement of α pd for these clusters was also performed by fitting a single power law over the same stellar mass range of (0.2-0.8) M ⊙ . In total, we have a sample of 35 GCs with measurements of (N(G1)/N(tot)) and α pd .
It is worth pointing out that the measurement of the slope of the present-day mass function by Ebrahimi et al. (2020) relies on the entire stellar census of the clusters. On the other hand, the fractions of G1 stars measured by Milone et al. (2017) were obtained within the HST/WFC3 (Wide Field Camera 3) field of view, and as such, these measurements are representative of the fractions in the central regions of the clusters. The fractions of G1 stars measured by Milone et al. (2017) underestimate their true values when measured for the entire clusters since G2 stars are often observed to be more concentrated toward the clusters' centers (e.g., Carretta et al. 2009;Kravtsov & Calderón 2021;see, however, Dalessandro et al. 2014 for a counterexample). Thus, the true values of the G1 star fractions are expected to be, on average, systematically higher than the values estimated by Milone et al. (2017).

Model
In the framework of the collision-based scenario, our aim is to understand the effect of the shape of the IMF of the G1 stars on the observed (N(G1)/N(tot)) − α pd anticorrelation. Understanding the precise effects of collisions requires calculating the collision rates over the entire stellar mass range and a knowledge of the coalescence efficiency as a function of stellar mass (e.g., Dib et al. 2007). Here, we used a simpler model with two mass bins, similar to the one used in Kravtsov et al. (2022). We briefly recall its basic elements. We assume that stars that can coalesce are those exclusively found in the stellar mass range (0.1-0.5) M ⊙ . The collision products (G2 stars) are more massive stars, with masses that fall in the range (0.5-0.9) M ⊙ . If N 1 (G1) and N 2 (G1) are the total numbers of G1 stars that fall in the target and product mass bins, respectively, then the fractions of these two population of stars to the total number of stars (N ′ tot ) can be written as and If N 1,ext (G1) = f ext N 1 (G1) is the total number of stars that have experienced a collision, under the assumption that we can only have two stars colliding to form a new one, and with f ext being the fraction of N 1 (G1) stars that undergo a collision, then the new number of stars becomes After the collision process, a second generation of stars is formed in the mass range (0.5-0.9) M ⊙ (i.e., G2 stars). The new fractions of stars in the mass ranges mentioned above will be given by and The value of (N(G1)/N(tot)) for RGB stars is given by its analog for main sequence (MS) stars as For the GCs in our sample, we have access to (N(G1)/N(tot)) and α pd . In the model, and for any given functional form of the IMF, we can measure the values of f 1,in and f 2,in , and for a given value of f ext we can calculate the values of f 1, f l , f 2, f l , and g 1 . With f 1, f l and f 2, f l we can also derive the value of the post-collision value of the slope, α pc . We calculated the value of α pc as being α pc = (log( f 1, f l ) − log( f 2, f l ))/(log(0.7) − log(0.3)), where 0.3 and 0.7 (in units of M ⊙ ) are the mean masses in the intervals (0.1-0.5) M ⊙ and (0.5-0.9) M ⊙ , respectively. We note that the parameter g 1 for MS stars is formally equivalent to the ratio (N(G1)/N(tot)) deduced from the observations of RGB stars. However, some caution should be exercised. As already noted in Kravtsov et al. (2022) and in Sect. 2, the fractions (N(G1)/N(tot)) have been measured in the central parts of GCs. Since RGB G2 stars are typically observed to be more centrally located than their G1 counterparts, the real (N(G1)/N(tot)) fractions should be systematically higher than the values that are actually derived. In contrast, the parameter g 1 in the models is assumed to be measured for the clusters as a whole. Furthermore, g 1 is assumed to be a proxy for MS stars that are now found on the RGB. However, the observed fraction (N(G1)/N(tot)) could also contain, in addition to RGB stars whose masses fall in the range ≈ (0.8-0.85) M ⊙ , a contribution from more massive blue straggler stars that are either of collisional origin or have been formed by mass transfer in binaries. For the reasons discussed above, we prefer to formally distinguish between (N(G1)/N(tot)), which is derived from the observations, and g 1 , which is calculated in the models, even if in practice the two quantities are expected to be very close.

Results
We describe the IMF of the G1 stars using a multicomponent power-law function, which is given by (Kroupa 2001) and where A is a normalization constant that depends on the cluster mass. Its determination is not necessary since we are only interested in the fractions of stars in the target and product mass bins. When calculating the fractions, we assumed the minimum and maximum stellar masses to be M min = 0.01 M ⊙ and M max = 150 M ⊙ , respectively. In Eq. 7 the value of α 0 is the most uncertain and displays a significant amount of scatter amongst the Milky Way clusters. We fixed its value to the one A&A proofs: manuscript no. 43556corr_final2 derived for the Galactic field, namely α 0 = −0.3. In the first instance, we considered the effect of varying the value of the slope at the high mass end, α 2 . The Galactic field value is ≈ −2.3, and we considered a range of values going from −1.7 to −2.7. The value of α 1 was set to the Galactic field value of −1.3. The shape of the IMF for various values of α 2 , for the same arbitrary mass reservoir, is displayed in Fig. 1. For each set of free parameters (α 2 and f ext ), we calculated the values α pc and g 1 using the model described in Sect. 3. The results in the g 1 − α pc space are displayed in Fig. 2 and are compared to the observational data. Variations in α 2 , coupled to variations in f ext in the range (0.2-0.6), can help explain a fraction of the scatter that is observed in the (N(G1)/N(tot)) − α pd relation but fall short of explaining the values of (N(G1)/N(tot)) that are 0.4.
It is evident that variations in (N(G1)/N(tot)) (equivalently in g 1 for the models) would be more significant if the primordial fractions of G1 stars were different, as a consequence of variations in the physical conditions of the star formation process. Dib (2014) and Dib et al. (2017) show that the characteristic mass (i.e., the peak in the IMF when fitted with a tapered power-law function) varies in the range ≈ 0.1 to ≈ 0.8 M ⊙ when measured for young clusters in the Milky Way. It is not clear whether the same level of variation is present in the IMF of GCs at the time they formed. Since we are using a broken power-law function with no well-defined peak, we varied the intermediate slope of the IMF in the mass range (0.08-0.5) M ⊙ (i.e., α 1 ). Figure 3 displays a number of IMF realizations for various values of α 1 and where α 0 and α 2 are assigned the Galactic field values of −0.3 and −2.3, respectively. The derived g 1 − α pc relations for these different cases and for various values of f ext are displayed in Fig. 4. As already shown in Kravtsov et al. (2022), varying the value of f ext for a fixed set of the IMF parameters reproduces the anticorrelation between (N(G1)/N(tot)) − α pd . However, in order to reproduce the observed level of scatter in the observations, it is necessary to allow, for a fixed value of f ext , the existence of variations in the values of α 1 . This level of variation in α 1 translates into a mean mass that varies between 0.2 and 0.55 M ⊙ . This level of variation is consistent with what is observed for young stellar clusters in the Milky Way (Dib 2014;Dib et al. 2017). The simple models presented here reproduce the observations remarkably well. As discussed in Sect. 2, the ratios of (N(G1)/N(tot)) measured by Milone et al. (2017) are likely to systematically underestimate their true values. If that is indeed the case, it would bring our models into an even better agreement with the observations. It is also worth pointing out that variations in α 1 and α 2 can occur simultaneously. Furthermore, we have restricted the distinct power laws to mass ranges similar to those inferred for the Galactic field. There is, however, no indication that these limits hold for individual clusters; they may well vary from cluster to cluster, as is observed in young clusters in the Milky Way (Dib 2014).
Dynamical interactions in the clusters that lead to the preferential ejection of low-mass stars can also cause the slope of the mass function in GCs to flatten with time (e.g., Webb & Vesperini 2016;Webb et al. 2017;Baumgardt & Hilker 2018;Baumgardt et al. 2019). Sollima & Baumgardt (2017) as well as Ebrahimi et al. (2020) plotted the value of α pd as a function of the cluster age (t age ) expressed in units of the present-day half-mass relaxation time (t rh ). The positive correlation between α pd and the ratio (t age /t rh ) suggests that the slope becomes flatter with increasing values of (t age /t rh ), and this effect can be understood in terms of the dynamical evolution of the clusters. De Marchi et al. (2010) fitted the present-day mass function of many Galactic clusters, young and old, with a tapered power-law function and  deduced that the characteristic mass varies between 0.1 and 0.8 M ⊙ , which is roughly the same range of variation implied by our model. However, De Marchi et al. attributed the difference in the characteristic mass to the preferential evaporation of low-mass stars from the clusters, which over time leads to a shift in the characteristic mass to higher values.
In the context of the formation of a second generation of stars, N-body models of GCs' evolution show that dynamical evolution preferentially removes low-mass stars from the clusters' centers, thus decreasing their fractions with time (e.g., Vesperini et al. 2013Vesperini et al. ,2021Sollima 2021). We note, however, that all of these N-body studies adopted a fixed form of the IMF (i.e., usually the Kroupa 2001 functional form) and did not allow for stellar collisions. In Fig. 5 (upper subpanel), we reproduce the Fig. 5. Value of the present-day slope of the stellar mass function (α pd ) plotted as a function of the ratio (log(t age /t rh )), where t age are the ages of the clusters and the t h their present-day half-mass relaxation time (upper subpanel) and fraction of G1 stars in the clusters plotted as a function of log(t age /t rh ) (lower subpanel). In both subpanels, clusters are segregated by their mass, which is taken to be the mean value of the photometric and dynamical masses. α pd versus log(t age /t rh ) plot of Ebrahimi et al. (2020), which is expanded by the inclusion of the additional clusters (see Sect. 2) and where clusters are segregated by their mass, which we take to be the mean value between the photometric and dynamical mass. Figure 5 also displays the values of (N(G1)/N(tot)) plotted as a function of log(t age /t rh ) (lower subpanel). Here again, the clusters are segregated by their mass. While it is possible to observe a weak anticorrelation between (N(G1)/N(tot)) and log(t age /t rh ) for low-mass GCs in the sample (Spearman's ρ coefficient ≈ −0.69), this anticorrelation is not detected for highmass clusters (Spearman's ρ ≈ 0.22). It is possible that this "dichotomy" could indicate a stronger and weaker role for dynamical interactions in low-and high-mass clusters, respectively. This is confirmed by the results of Dalessandro et al. (2014), who find that G1 and G2 stars are entirely spatially mixed in the low-mass GC NGC 6362. Inversely, this would also indicate a more pronounced contribution from collisions in high-mass clusters and a weaker one in their low-mass counterparts. While it is probably too early to ascertain this fact from the current data, N-body models that follow the dynamical evolution of GCs while accounting for stellar collisions can help shed more light on this issue.

Discussion and conclusions
In this work we explore, in GCs that harbor multiple populations, the origin of the anticorrelation that is observed between the fraction of the first generation of stars (N(G1)/N(tot)) and the slope of the present-day mass function of the clusters in the stellar mass range (0.2-0.8) M ⊙ (α pd ). We compare the observations to the results of a simple model that is based on the idea of stellar collisions between G1 stars in the mass range (0.1-0.5) M ⊙ that lead to the formation of a second generation of stars with masses in the range (0.5-0.9) M ⊙ . The model has two input parameters, the shape of the IMF of the G1 stars and the fraction of G1 stars in the (0.1-0.5) M ⊙ mass range that collide to form G2 stars (parameter f ext ). The parameter f ext encapsulates much of the physics that governs the efficiency of the collision process. Its value depends on the compactness of the protocluster cloud, and hence on the mean distance between stars and also on other factors such as the initial levels of mass segregation in the clusters, the binary fraction, and the period distribution of binaries. We find that the appropriate range for f ext is ≈ (0.2 − 0.6). Smaller and larger values of f ext lead to values of the slope that are outside the range of observed values.
Our results show that variations in f ext are necessary to explain the anticorrelation between (N(G1)/N(tot)) and α pd for a fixed shape of the IMF. However, the large scatter that is observed in the (N(G1)/N(tot))−α pd anticorrelation can only be explained, in the framework of this collision-based model, by variations in the IMF of the G1 stars. In particular, we show that variations in the slope of the IMF in the intermediate-mass regime of a Kroupa-like IMF (≈ 0.08 − 0.5) in the range −1.9 to −0.9 can reproduce the scatter in (N(G1)/N(tot)) at a given value of α pd . This level of variation in α 1 corresponds to a range of ≈(0.2-0.55) M ⊙ in the mean stellar mass in the clusters. In this work we have restricted the target G1 stars to the mass range (0.1-0.5) M ⊙ . However, mergers between stars whose masses fall in the range (0.1-0.9) M ⊙ and that result in G2 stars with masses ≤ 0.9 M ⊙ can further decrease the values of g 1 and bring the models into a better agreement with the observations. Variations in the IMF of GCs have been pointed out by other groups. Zonoozi et al. (2016) show that a dependence of the high mass end of the IMF on metallicity as proposed by Marks et al. (2012) can help explain the mass-to-light versus metallicity anticorrelation that is observed for the population of GCs in M31.
While it is probably safe to say that the jury is still out concerning the relative importance of collisions versus dynamical interactions in modifying the shape of the IMF in GCs, the collision-based model presented here provides an explanation for the observed (N(G1)/N(tot)) − α pd anticorrelation, and it is not yet clear if dynamical evolution can do the same. The level of variations in the IMF inferred in this work, particularly for the slope in the stellar mass range (0.08-0.5) M ⊙ that is needed in order to explain the scatter in the (N(G1)/N(tot)) − α pd relation, is consistent with what is observed for young stellar clusters in the Milky Way (Dib 2014) and may reflect variations in the physical conditions of the GC parental clouds at the time the G1 population started to form, the occurrence of other processes such as gas accretion onto protostars (Dib et al. 2010), or collisions between protostars before they settle on the MS (Dib et al 2007).