The role of rotation on the formation of second generation stars in globular clusters

ABSTRACT

1 INTRODUCTION

Spectroscopic and photometric studies have now clearly revealed that globular clusters host multiple stellar populations with different chemical compositions. A significant fraction of GC stars show, in fact, ‘anomalous’ chemical composition rarely found in the field, such as an enhancement in He, Al, and Na and a depletion in O, Mg (Piotto et al. 2005; Carretta et al. 2009; Milone et al. 2017; Gratton et al. 2019; Marino et al. 2019; Masseron et al. 2019). Stars are generally divided, depending on their chemical composition, into two or more groups: the first population (or first generation; hereafter FG) share the same chemical pattern with field stars, while the second population (or second generation; hereafter SG) display anomalous abundances. Multiple populations within the same GC have been detected not only in the Milky Way GCs, but also in external galaxies such as the Magellanic Clouds (Mucciarelli et al. 2009), Fornax galaxy (Larsen et al. 2014), and M31 (Nardiello et al. 2019).

So far, a complete picture explaining the origin of anomalous stars, and therefore the formation of multiple stellar populations within GCs, is still lacking. Several scenarios have been proposed to address this topic, nevertheless none of them is able to reproduce all observational constraints (Renzini et al. 2015; Bastian & Lardo 2018; Gratton et al. 2019). The main difference between the various scenarios resides in the identification of the stars providing the processed gas out of which SG stars formed. Possible sources of processed gas, suggested in the literature, include asymptotic giant branch (AGB) stars (D’Ercole et al. 2008; D’Ercole, D’Antona & Vesperini 2016), fast-rotating massive stars (Decressin, Charbonnel & Meynet 2007), massive stars (Elmegreen 2017), supermassive stars (Denissenkov & Hartwick 2014; Gieles et al. 2018), massive interacting binaries (de Mink et al. 2009; Bastian et al. 2013), black holes accretion discs (Breen 2018), and stellar mergers (Wang et al. 2020).

Among the aforementioned scenarios, the most thoroughly studied to date is the AGB one, which is also the one we are adopting in this work. In this scenario, after the formation of FG stars, FG massive stars clear the system from both the gas leftover by the FG formation and the one ejected by massive stars themselves (Calura et al. 2015). SG are then formed out of the ejecta of FG AGB stars plus pristine gas with same chemical composition of the gas from which FG were formed (D’Ercole et al. 2010, 2012; D’Antona et al. 2016).

D’Ercole et al. (2008) performed, for the first time, hydrodynamic simulations of a star-forming GC in the AGB framework finding that the gas ejected by AGB stars collects in a cooling flow towards the cluster centre. SG stars form in the central regions of the FG system out of this gas. The predicted central concentration of SG stars relative to the FG population has been observed in several clusters where some memory of the initial differences between generations has been preserved. Various studies have also shown that AGB ejecta alone do not reproduce the observed anticorrelations, such as the Na–O one (Carretta et al. 2009). In order to produce the observed chemical patterns, AGB ejecta should be diluted with pristine gas (namely gas with the same chemical composition as the one out of which FG were formed) while SG stars are formed (D’Antona et al. 2016, D’Ercole et al. 2012, 2016). Dilution of processed gas with pristine gas is required by many of the aforementioned scenarios (Gratton et al. 2019). Calura et al. (2019) performed the first series of 3D hydrodynamic simulations of a massive proto-GC in the AGB framework, following up the work of D’Ercole et al. (2008) which is instead based on 1D simulations. They take into account the radiative cooling and modelled, together with the AGB feedback, the accretion of pristine gas. They found, in addition to D’Ercole et al. (2008) results, that the most helium enriched SG stars, which are also the first to be formed, are characterized by a more compact spatial distribution than the less helium enriched ones. This trend is consistent with that found in observational studies of the SG populations spatial distribution (Johnson & Pilachowski 2012; Simioni et al. 2016).

In most studies on GC formation, clusters are modelled as non-rotating systems; however, in the last few years an increasing number of GCs have been found to show signatures of internal rotation (see e.g. Pancino et al. 2007; Bellazzini et al. 2012; Fabricius et al. 2014; Lardo et al. 2015; Boberg et al. 2017; Bianchini et al. 2018; Ferraro et al. 2018; Kamann et al. 2018; Lanzoni et al. 2018a,b; Sollima, Baumgardt & Hilker 2019; Vasiliev & Baumgardt 2021). The observed internal rotation is generally found to be quite moderate with typical ratios between the rotational amplitudes to the central velocity dispersion (V_rot/σ₀) ranging from 0.05 to about 0.7 (Bellazzini et al. 2012; Fabricius et al. 2014; Kamann et al. 2018). Present-day rotation, however, is most likely the remnant of a stronger early rotation (Hénault-Brunet et al. 2012; Mapelli 2017) which has been then lessened under the effects of two-body relaxation (Bianchini et al. 2018; Kamann et al. 2018; Sollima et al. 2019) and tidal forces exerted by the host galaxy. During the long term evolution of the cluster, the combined effect of angular momentum transport from the cluster centre outwards and of escaping stars carrying it away, leads to a loss of angular momentum and therefore to a decrease of the rotational amplitude (see e.g. Einsel & Spurzem 1999; Tiongco, Vesperini & Varri 2017; Livernois et al. 2022).

It has been shown that internal rotation can strongly affect GC long-term evolution, in particular both reducing the relaxation time-scale (or shortening the time of core-collapse), and therefore accelerating its dynamical evolution, and the mass-loss rate (Einsel & Spurzem 1999; Kim et al. 2002; Kim, Lee & Spurzem 2004; Ernst et al. 2007; Kim et al. 2008; Hong et al. 2013; Mastrobuono-Battisti & Perets 2021). Moreover, rotation affects also the present-day morphology (Hong et al. 2013), as well as the dynamics of multiple stellar populations (Mastrobuono-Battisti & Perets 2013, 2016; Hénault-Brunet et al. 2015; Tiongco, Vesperini & Varri 2019; Mastrobuono-Battisti & Perets 2021).

The observational study of the kinematics of multiple stellar populations is still in its early stages, but a few investigations have already revealed differences in the rotation amplitudes between stellar populations (Lee 2015, 2017, 2018; Cordoni et al. 2020a; Dalessandro et al. 2021; Szigeti et al. 2021), with the SG component rotating faster than the FG, with the exception of ω Cen where Bellini et al. (2018) found opposite results. Recently, Cordero et al. (2017) and Kamann et al. (2020) analysed the rotational patterns of two GCs subpopulations and found that extreme SG stars (Na-enhanced and extremely O-depleted) are characterized by a larger rotational amplitude than intermediate SG ones (moderately Na-enhanced and O-depleted). In some cases, instead, no difference as been detected in the rotational patters of distinct populations (Milone et al. 2020; Cordoni et al. 2020a,b; Szigeti et al. 2021). This rotational homogeneity between different populations could hint that these systems are dynamically evolved, where stars are already well kinematically mixed.

Rotation is, therefore, a fundamental physical process which needs to be included in theoretical models, however, very few studies have investigated its effects on the formation of multiple populations in GCs so far. Bekki (2010, 2011) and more recently McKenzie & Bekki (2021), studied the effects of rotation on the SG formation in the AGB framework, through 3D hydrodynamic simulations. They found that SG and FG stars display significant differences both in the rotational velocity (V_rot) and in the velocity dispersion (σ) profiles, with SG stars rotating faster but with a smaller dispersion. Moreover, SG stars show a flattened and compact spatial distribution which should be smoothed out during the cluster evolution in order to match the observations. These initial kinematical differences between various stellar populations are found to significantly decrease during the subsequent evolution of the cluster on two-body relaxation time-scales, which points towards a spatial and kinematical mixing of the populations (Mastrobuono-Battisti & Perets 2013; Vesperini et al. 2013; Hénault-Brunet et al. 2015; Mastrobuono-Battisti & Perets 2016; Tiongco et al. 2019; Vesperini et al. 2021).

In this work, we extend the study of Calura et al. (2019, hereafter C19) to investigate the effects of FG rotation on the formation of SG stars in a massive proto-GC, through a series of 3D hydrodynamic simulations. In our simulations, we will explore how the interplay between the kinematics of the AGB ejecta released by a rotating FG system and the accreted pristine gas affect the final kinematic properties of SG stars and the dependence of these properties on the chemical composition for different SG subpopulations.

The paper is organized as follows: in Section 2 we present the setup of our simulations focusing on the implementation of the rotating FG component. Section 3 contains the results of out models and describes the dynamics of both the gas and the stellar component. In Section 4, we discuss the outcomes of the simulations and compare them with previous results in the literature. Finally, we summarize our conclusion in Section 5.

2 SIMULATION SETUP

The initial setup of this work is similar to the one of C19 with the difference that we are here assuming a rotational FG. Therefore, all our simulations start |${t_{\rm AGB}=\rm 39\,Myr}$| after the FG formation, once all massive stars have already exploded as core-collapse supernovae (CC-SNe) and the system is completely cleared out of both SN-enriched ejecta and pristine gas (Calura et al. 2015). At this time, the intermediate mass stars are undergoing their AGB phase, returning mass and energy in the gas-free system.

We perform 3D hydrodynamic simulations using a customized version of the adaptive mesh refinement (AMR) code ramses (Teyssier 2002), which uses an unsplit second-order Godunov method to solve the Euler equations of hydrodynamics. The gas obeys the adiabatic equation of state for an ideal monoatomic gas P ∝ ρ^γ, where P and ρ are the pressure and density of the gas and the adiabatic index γ is set to 5/3. We take into account various astrophysical processes such as the mass and energy return from AGB stars, radiative cooling and star formation. For simplicity, the FG system is modelled as a static Plummer (1911) density profile with mass |$M_{\rm FG}=10^7 {\rm M_{\rm \odot }}$| and Plummer radius a = 23 pc. On the other hand, SG stars are instead modelled as collisionless particles; their dynamic evolution is derived by means of a Particle-Mesh solver. This paper is focused on the formation and early dynamics of multiple populations and an investigation of the origin of the various chemical patterns is beyond the scope of this work. Here, we only trace the evolution of the helium abundance of both the gas and the stars and use this abundance to identify the various SG subgroups. In all our simulations, we adopt a size of the computational box of L³ = (292 pc)³ uniformly divided into (512)³ cells, which corresponds to a spatial resolution of 0.6 pc.

2.1 Initial setup

As in C19, we follow the scenario of D’Ercole et al. (2016) and therefore assume that the cluster is located in the disc of a high-z star-forming galaxy (Kravtsov & Gnedin 2005; Bekki 2012; Kruijssen 2015). In addition, we assume that the cluster is orbiting in its host galaxy, which leads to an asymmetric accretion of gas from the side towards which the system is moving. In models, this is accomplished by placing the cluster at the centre of the computational box and, at time t_inf, gas starts to flow from one of the boundaries. Such accretion leads to a dilution of the AGB ejecta with pristine gas which, as discussed in the Section 1, is required in all the models to match the chemical trends derived by observations. The beginning of the infall of pristine gas, however, does not correspond to the beginning of the simulation. As in D’Ercole et al. (2016), we assume that FG CC-SNe explosions have carved a large cavity around the cluster, an hypothesis confirmed also by 3D hydrodynamic simulations of the expansion of SN driven bubbles located in galactic discs (Hopkins, Quataert & Murray 2012; Creasey, Theuns & Bower 2013; Walch et al. 2015; Kim, Ostriker & Raileanu 2017).

At the beginning of the simulation, the boundaries of all six faces of our computational box are set to be outflow. At t_inf, the pristine gas flows into the box though the yz plane at negative x.

Table 1 contains a summary of the main parameters assumed for our simulations, whereas Table 2 reports a description of the models we have performed with a particular emphasis on the rotational prescriptions we have adopted.

We point out that our model is aimed at exploring the formation of SG stars in a massive cluster with a present-day mass ≈10⁶ M_⊙. In a future investigation, we will explore the implications of initial rotation for a broader range of initial conditions including clusters with different masses and structural parameters (see e.g. Yaghoobi et al. 2022 for a study of the formation of multiple populations in non-rotating clusters with different initial masses).

2.2 Physical ingredients

We summarize here the main physical ingredients which are included in our simulations, with a particular focus on the implementation of rotation in the FG system (for more details on the basic setup, see C19 and Lacchin, Calura & Vesperini 2021).

The FG is instead modelled as a static component spatially distributed following the analytical Plummer mass density profile with M_FG = 10⁷M_⊙ and a = 23 pc. The system is evolved for 65 Myr to compare our results with the one of C19.

Here, we trace the helium composition both of the gas and of SG stars. We assume, for AGB stars the yields of Ventura & D’Antona (2011), therefore the helium mass fraction of the AGB ejecta spans between 0.36 for the gas released by the most massive AGBs to 0.32 for AGB ejecta at the end of the simulations. No iron is assumed to be produced inside AGBs, therefore the iron composition of their ejecta is the same as the one of the pristine gas.

The source terms |$\dot{\rho }_{\rm \star ,AGB}$| and S are added at each time-step to the density and energy of the fluid, respectively, as well as the cooling term which takes into account the loss of energy due to radiation (see Few et al. 2014). The pristine gas inflowing the system has a fixed temperature of T = 10⁴ K, a typical value for the warm photoionized ISM in a star-forming galaxy (Haffner et al. 2009).

2.3 Rotation

Signatures of internal rotation are found both in simulations of star cluster formation through the collapse of giant molecular clouds (Mapelli 2017; Ballone et al. 2020; Chen, Li & Vogelsberger 2021) and from observations of young clusters (see e.g. Fischer et al. 1992; Fischer, Welch & Mateo 1993; Hénault-Brunet et al. 2012).

In Appendix  A, we also report the results of models assuming a solid body rotation. Such profile has been adopted by some works in literature such as Bekki (2010, 2011) and McKenzie & Bekki (2021). The angular velocity ω has been derived assuming a velocity at the Plummer radius v_rot(r = a) = 2.5 km s⁻¹.

3 RESULTS

In this section, we present the results of the models assuming that the FG internal rotation follows the radial profile in Eq. 6. We have explored two different orientations for the rotational axis (parallel and perpendicular to the infall, see Tiongco, Vesperini & Varri 2018, 2022 for N-body studies of the long-term evolution of clusters with different orientations of the rotation axis) studying their effects for various pristine gas densities (⁠|$10^{-24}-10^{-23}\, {\rm g\, cm^{-3}}$|⁠). All model parameters are reported in Table 2.

The results obtained by adopting a solid-body rotation are very similar to those found assuming the profile in Eq. 6; for this reason, here we present only the results obtained for the models with the analytical profile in Eq. 6. We discuss the comparison with the models adopting a solid-body rotation in Section  A.

3.1 Rotation perpendicular to the infall

In this section, we present the results obtained assuming that the rotation axis of the FG system is perpendicular to the direction of the external gas infall and therefore lies along the z-axis.

3.1.1 Low-density model

Gas evolution

In Fig. 1, we show the 2D density, temperature and helium mass fraction slices for the gas component at 26 and 65 Myr, for the model with a low pristine gas density. To build the slices, we select all cells on the xy plane passing through the centre of the computational box.

Figure 1.

2D maps of the gas density (left-hand panels), of the temperature (central panels) and the helium mass fraction (right-hand panels) on the x–y plane for the LDanaz simulation. The corresponding evolutionary time of each set of panels is reported at the top. The white arrows in the helium mass fraction maps represent the gas velocity field.

Open in new tabDownload slide

We have displayed, on top of the helium mass fraction map, the gas velocity field as white arrows. For a comparison with C19, we show here the maps at similar evolutionary times.

In the first snapshot, at 26 Myr, the infall has already started. The shock front has just crossed the system creating two ‘wings’ which are both visible in all the three maps. Contrary to C19, the wings are not perfectly symmetric because of the rotation of FG stars and therefore of their ejecta on the xy plane. The presence of the infalling gas is even more clearly visible in the maps of the helium mass fraction where the helium-rich gas ejected by FG stars, marked in red, is located downstream to the cluster and in the central regions while the blue helium-poor pristine gas of the infall is dominant upstream to the system. In the central part, the system is composed of several cold (T ∼ 10³ K) and dense (⁠|${\rm \rho =10^{-19} \, g \ cm^{-3}}$|⁠) clumps. These structures are the result of the fragmentation of a torus which has been recently described in detail by Inoue, Yoshida & Hernquist (2021) through a linear perturbation analysis. In codes which exploit the AMR technique, the computational box is discretized in multiple grids, which in turn leads to a discretization of the hydrodynamic equations. Such modellization can lead to the so called ‘artificial fragmentation’ which takes place when small perturbations grow to form fragments. Truelove et al. (1997) studied the conditions that give rise to this phenomenon, by means of a 3D hydrodynamic AMR code, finding that, to avoid artificial fragmentation, the Jeans length λ_J (Jeans 1902) should be at least four times greater than the resolution element. For the gas properties in our simulation, we end up with a λ_J is equal to 2.6 pc and therefore satisfies the Truelove criterion for avoiding artificial fragmentation. The clumps found in our simulation are thus not numerical artifacts but rather represent the result of the actual dynamics of the gas out of which SG stars form. A detailed investigation of the formation and evolution of these clumps is beyond the scope of this work and will be further investigated in a future study.

At the end of the simulation, ∼65 Myr, a tail of cold and dense gas is pointing towards the cluster centre, the so called ‘accretion column’. This helium-poor gas is generated out of the infall event through the Bondi–Hoyle–Litterton accretion mechanism (Bondi & Hoyle 1944; Moeckel & Throop 2009).

As shown in Lacchin et al. ( 2021), the structure and dynamics of the accretion column may be significantly affected by the feedback of SNe Ia.

Evolution of the stellar component

In Fig. 2, we plot the 2D surface density, helium mass fraction and line-of-sight velocity maps on the xy and yz planes for the stellar component at the same evolutionary times shown in Fig. 1.

Figure 2.

2D maps of the stellar component at two evolutionary times for the LDanaz simulation on the x–y plane and y–z plane at two evolutionary times (reported on the left-hand panels). The first column shows the surface density, the second represents the helium mass fraction Y, while the third the line-of-sight velocity of the stars.

Open in new tabDownload slide

In the 26 Myr snapshots, stellar particles are preferentially lying on a disc in the xy plane; this is the expected distribution of stars forming out of the gas released from an FG system rotating about the z-axis. The surface density map in the xy plane clearly reveals the presence of five clumps in the very central part of the system. Their average mass is 4 × 10⁴M_⊙ and more than |$90{{\ \rm per\ cent}}$| of all the particle mass is confined in these clumps. The five clumps are formed and orbit the centre of the system, and eventually merge after a few Myr.

All stars at this stage are strongly helium-enhanced, given that the infall has just started and has not had enough time to significantly dilute the AGB ejecta. These stars are concentrated in the cluster’s innermost regions (their half-mass radius is r_h = 3.7 pc) and distributed in a toroidal structure; they are characterized by a significantly more rapid rotation than that of the FG population. Specifically, the SG peak rotational velocity is equal to about V_peak = 16 km s⁻¹ at a distance from the cluster’s centre of R_pk = 6.7 pc.

At 65 Myr the clumps have already merged and collected in the centre. The disc is still clearly visible in all maps and contains |$80{{\ \rm per\ cent}}$| of the total mass of the SG. It extends for ∼8 pc and is composed of helium-enhanced stars that are still significantly rotating. It is interesting to note the presence of a small region in the innermost part of the disc characterized by a helium abundance smaller than in the rest of the disc. This feature is a consequence of the early dynamics of the SG formation: as discussed above, the most He-rich AGB ejecta are first collected in a clumpy toroidal structure, which later forms a stellar disc; the non-rotating pristine gas, on the other hand, flows directly in the centre and creates this small region with slightly lower helium abundance. This is also illustrated by the density profiles presented in Fig. 3 showing a slightly lower central density of the most He-rich population. In order to link these fine structural and chemo-dynamical features within the SG system with the present-day properties of old GCs, it is necessary to carry out further explorations of the long-term dynamical evolution with initial conditions similar to those emerging from our simulations.

Figure 3.

Total SG density profile at the end of the simulation and SG density profiles for three intervals of the helium mass fraction Y, for the low-density model LDanaz. The density profile of FG stars is also shown (see the legend for the details).

Open in new tabDownload slide

In Fig. 4, we show, for the LDanaz model, the angular variation of the Cartesian components of the velocity of the SG stars (left-hand panels) and the radial profiles of the mean velocity and velocity dispersions in cylindrical coordinates (right-hand panels) measured at the end of the simulation (t = 65 Myr). This figure further illustrates the kinematical properties of the SG population and the various SG subgroups. SG stars rotate more rapidly, with an amplitude of ∼12 km s⁻¹, than the FG population, whose amplitude is equal to 2.5 km s⁻¹. Such difference in the rotational intensity between the two populations has been observed in several present-day clusters (Lee 2017, 2018; Dalessandro et al. 2021; Szigeti et al. 2021), although the amplitudes are significantly smaller than the once we obtain, since our simulations do not include the long-term evolution of the system. In addition, more helium-enriched stars are found to rotate faster than weakly enriched ones, qualitatively in agreement with recent observational works on M13 (Cordero et al. 2017) and M80 (Kamann et al. 2020) (with the exception, already discussed above, of the rotation in the very inner regions, where the intermediate SG group rotates slightly faster). The small angular variation in the mean value of V_z is a consequence of the fact that the SG disc is slightly tilted (see also Fig. 2). Moreover, we found that FG and SG stars are rotating in phase, a feature often found in present-day GCs (Milone et al. 2018; Cordoni et al. 2020b).

Figure 4.

Stellar rotation profiles for the model LDanaz at 65Myr. On the left: Rotation amplitude of the Cartesian velocity components for SG stars as a function of |$\theta$|⁠, defined as the angle between the direction of each star projected on the |$x-y$| plane and the x axis, for three bins of the helium mass fraction Y (see the legend for more details). The FG rotation amplitudes are also shown by the black line. On the right: Radial profiles of the cylindrical velocity components (solid lines) and velocity dispersions |$\sigma$| (dashed lines) for SG stars for three bins of the helium mass fraction Y.

Open in new tabDownload slide

Regarding the velocity profiles in Fig. 4, it is not surprising that the peak of the SG rotational velocity V_t is located approximately at the SG half-mass radius equal to 1.5 pc, as suggested by observational studies of very young stellar clusters (Hénault-Brunet et al. 2012), but also of older globular clusters (Bianchini et al. 2018).

As for the radial profiles of the velocity dispersion shown in the right-hand panels of Fig. 4, it is interesting to note the presence of a few bumps indicating that the SG system is still in the process of settling into a dynamical quasi-equilibrium. Along the tangential direction, | V_rot/σ| is greater than unity in the disc region, with a peak of 3.5 at the edge of the stellar disc, confirming that the system is rotationally supported. This ratio is much larger than the observed one, but, as for the rotational amplitudes of the two populations, the long-term evolution has been shown to lead to a significant decrease of this quantity (Tiongco et al. 2017). Finally, we point out that here we have explored a single case for the initial FG rotational profile; an extension of this study to explore different initial FG rotational properties will be carried out in the future. The early study on the formation of SG in rotating clusters by Bekki ( 2010) suggests that the strength of the SG rotation emerging at the end of the formation phase is correlated with that of the FG system.

3.1.2 High-density model

Gas evolution

Fig. 5 represents the 2D density and temperature slices at 5 and 62 Myr for the HDanaz model. At 5 Myr, the system is encountering the infalling gas which creates a shock front that is, as in Fig. 1, not perfectly symmetric because of the rotation of the already present AGB ejecta. The presence of a bow shock is also visible in the temperature map, at a distance comparable to the Plummer radius of the cluster (see also Naiman, Ramirez-Ruiz & Lin 2011). In the central 10 pc, the gas is rotating counter-clockwise due to the imposed rotation of the FG stellar component. In particular, the gas is mainly confined in a torus due to the balancing effects of gravity and rotation.

Figure 5.

2D slices of the gas density (left-hand panels) and temperature (right-hand panels) on the x–y plane for the HDanaz simulation. The corresponding evolutionary time is reported on the top of each set of panels. The white arrows in the density maps represent the velocity field of the gas.

Open in new tabDownload slide

At 62 Myr the gas is not rotating anymore and, similarly to the results obtained by C19, the accretion column is carrying pristine gas towards the cluster centre.

We do not report here the maps of the gas helium mass fraction since, at 5 Myr, the map is similar to the one of the LDanaz model at 26 Myr and, at the end of the simulation, no feature is visible since the entire computational box is dominated by gas with moderate helium enhancement.

Evolution of the stellar component

Contrary to the LDanaz model, here the earlier incoming of the infall prevents the clump formation while the disc is still formed as shown in Fig. 6 in the snapshot at 5 Myr. Not surprisingly, SG stars are mainly located in a torus resembling the gas distribution.

Figure 6.

2D surface density maps of the stellar component at two evolutionary times for the HDanaz simulation on the x–y plane, on the left, and y–z plane, on the right.

Open in new tabDownload slide

At 62 Myr the stellar component is much more extended than at the previous snapshot and a significant number of stars form along the accretion column. Small clumps form along the accretion column, move towards the cluster centre and merge with the central component. Instead, a minority of them is kicked away.

In contrast to the LDanaz model, in the HDanaz one the disc does not survive once the infall starts, both because of the higher density of the accreted gas, and, in turn, its stronger ram pressure, but also of because of its earlier arrival when the disc is still forming. The system is dominated by moderately helium enriched stars which are formed out of the AGB ejecta and a large amount of the infalling pristine gas. Nevertheless, a small fraction of highly enhanced helium stars preserves the rotation acquired at birth.

Fig. 7 shows the velocity profiles as a function of the θ angle (left, defined as in Fig. 4) and radius (right) for the HDanaz model. SG stars rotate with an amplitude similar to that of the FG with peaks in the x component. This feature is due to the stars born along the accretion column and consequently possessing a significant velocity towards −x. Along z, the velocity is not zero due to the gas motions produced by the infall and the subsequent formation of the accretion column.

Figure 7.

Stellar rotation profiles of the model HDanaz at 65 Myr. The reported quantities are as in Fig. 4.

Open in new tabDownload slide

At variance with the low-density model, here the dispersion velocities of all the various SG subgroups are similar. The maximum value of |V_rot/σ| is 0.2, much smaller than in the low density models confirming that the system is dispersion supported.

3.2 Rotation parallel to the infall

We present here the results for the model where the FG is assumed to rotate about the x-axis and therefore with a rotational axis parallel to the direction of the infall. We describe here only the low density model (LDanax), given the similar results obtained for the HDanaz and HDanax models. In both cases a disc is formed both in the gas and stellar components, but, once the infall starts, it is disrupted and no significant rotation is found.

3.2.1 Low-density model

Gas evolution

In Fig. 8 we show the 2D density, temperature, and helium mass fraction slices at almost the same times as in Fig. 1. At 28 Myr, the shock front due to the infalling event has already overcome the cluster and the accretion column has formed. However, the shape of the accretion tail is very different from the one in Fig. 1 since here the rotation leads to the formation of a disc on the yz plane.

Figure 8.

Two-dimensional maps of the gas density (left-hand panels), of the temperature (central panels) and helium mass fraction (right-hand panels) on the |$x-y$| plane for the LDanaz simulation. The time is reported on the top of each set of panels. The white arrows in the helium mass fraction maps show the gas velocity field.

Open in new tabDownload slide

The infalling gas is colliding with the disc face-on and is forced to circulate around the disc giving rise to two accretion columns divided by a very narrow tail, still not polluted by the pristine gas. This feature is even more clear in the helium mass fraction map, where a tight tail of extremely helium enhanced gas is present downstream of the system. The two helium-enriched blobs near the centre are instead due to the bounce of the helium-rich gas first pushed by the infalling gas downwards. to explain

At 64 Myr, the tail of extremely helium-enhanced gas has disappeared, while the accretion column has acquired its standard shape. In the central regions, the gas is still enriched in helium, whereas in the outskirts and along the accretion column the AGB ejecta have been strongly diluted by the infalling gas and have an helium abundance very similar to that of the pristine gas.

Evolution of the stellar component

Fig. 9 shows the 2D maps for the stellar component at the same evolutionary times of Fig. 8.

Figure 9.

2D maps of the stellar component for the LDanax model. The reported quantities are as in Fig. 2.

Open in new tabDownload slide

At 28 Myr, stars with low helium enhancement start to form both downstream of the system along the two cold and dense tails clearly seen in Fig. 8 and upstream of it at around 20pc from the cluster centre. At this distance, a shock is formed due to the interaction of the infalling gas with the potential well of the cluster, which induces the formation of new stars. Four clumps of extremely helium enhanced stars are formed in the very centre and located in the disc similarly to what has been obtained for the LDanaz model.

At 64 Myr, the stellar disc is still present and particularly evident in all the maps in Fig. 9. Similarly to the LDanaz model, we find that the very inner regions host stars slightly less enriched in helium. The disc is slightly thicker than in the LDanaz model; this is due to the dynamical effect of the gas infalling along the accretion column in the direction perpendicular to the plane of the disc. In addition, the disc appears slightly tilted by about 10 deg relative to the x–y plane. The left-hand panel of Fig. 10 shows that the SG component has a higher rotational amplitude than the FG one, with a peak value similar to the one found for the LDanaz model, and, in addition, SG stars with high helium composition are rotating faster than those with moderate helium enrichment. The different orientation of the rotation axis with respect to the direction of the infall implies that the external gas motions do not affect the rotation plane of the disc, leading to a smoother rotation profile as shown in the left-hand panel of Fig. 10.

Figure 10.

Stellar rotation profiles of the model LDanax at 65 Myr. The reported quantities are as in Fig. 4.

Open in new tabDownload slide

Similarly to the Ldanaz, | V_rot/σ| is increasing inside the disc region with a peak of 2, therefore the system is rotationally supported even at large radii.

4 DISCUSSION

Our simulations have explored the role of cluster rotation on the formation of SG stars in a massive proto-GC. We have shown how different inclinations and pristine gas densities affect the kinematical and structural evolution not only of the FG and SG, but for the first time, of different subpopulations. We here discuss our results, connecting them with the relevant literature both from the observational and theoretical side.

4.1 On the chemical composition of SG stars

While cluster rotation has a significant effect on the morphological and kinematical properties of SG stars, it does not have a strong impact on the overall chemical composition of the system and its final stellar mass. In Fig. 11, the helium distribution functions of both the low and high density models are compared. As shown in this figure, for both the low-density and high-density models, the helium distribution is not significantly affected by the orientation of the rotation axis relative to the direction of motion of the cluster in the external gas environment. The two panels are also in general agreement with the helium distribution found in the non-rotating models by C19.

Figure 11.

The mass distributions of Y in SG stars at the end of the simulations (⁠|$\sim 65$|Myr) for different models. On the left: low-density models: LDanaz (blue), LDanaz (orange), LDsbz (purple) and LD (green). On the right: high-density models: HDanaz (blue), HDanax (orange) and HD (green). Each distribution is normalised to its maximum value. The pristine gas helium mass fraction Y is shown by the black dashed line.

Open in new tabDownload slide

4.2 Comparison with literature

Recently, McKenzie & Bekki (2021) studied the formation of MPs in proto-GCs within the context of their parent galaxy, through 3D hydrodynamic simulations. In all their models, the cluster is assumed to rotate as a solid body. At variance with our simulations, the infalling gas is modelled as clumpy. For their fiducial model, which assumes a FG of mass M_FG = 10⁶M_⊙, an order of magnitude lower than in our case, they found a flattened distribution of the SG, similarly to Bekki (2010, 2011) and this work. Discs, however, are not a common feature in GCs, which are generally described as spheroidal systems, although some degree of ellipticity has been reported by Frenk & Fall (1982) who found an age–ellipticity relation in Galactic and LMC GCs. It is however verified by N-body simulations that SG stars born in a disc mix with FG ones through angular momentum exchange, ending up with a flattened system after 12 Gyr of evolution (see also Mastrobuono-Battisti & Perets 2016; Tiongco, Collier & Varri 2021). At variance with our results, instead, McKenzie & Bekki (2021) found a less concentrated SG compared with the FG, a configuration found only in few GCs. On the other hand, in both our studies, FG and SG have been derived to rotate in phase, with the SG rotating significantly faster than the FG. Even though we are exploring here the evolution of a much more massive system, our resulting SG rotational amplitude is in agreement with the mass–amplitude relation found by McKenzie & Bekki (2021).

As in McKenzie & Bekki (2021), we have tested different inclinations of the GC rotational axis with respect to the one of the host galaxy (which is here identified by the direction of the infalling gas). However, in the low-density models, where the disc survives, the inclination of the disc is not significantly affected by the infalling gas and, therefore, it does not align with the host galaxy as it instead happens in their simulations. The models presented here and those of McKenzie & Bekki (2021) differ in various aspects (FG initial mass, external infalling gas, and time interval explored). Further investigation is necessary to clarify the origin of the differences between the final SG system orientation found in the two works.

4.3 Dynamical and kinematical features of present-day globular clusters

Our simulations follow the formation and very early evolutionary phases of a cluster composed by multiple stellar populations. For a close comparison with the properties of present-day old and intermediate-age GCs, it would be necessary to take into account the effects of the subsequent long-term dynamical evolution. A number of studies have shown that the long-term evolution gradually modifies the properties imprinted by the formation processes. Although, in some dynamically old clusters, the differences between the dynamical properties of SG and FG stars may be completely erased during the clusters’ long-term evolution, in dynamically younger ones some memory of these differences may still be preserved (see e.g. Mastrobuono-Battisti & Perets 2013, 2016; Vesperini et al. 2013; Hénault-Brunet et al. 2015; Tiongco et al. 2019; Sollima 2021; Vesperini et al. 2021; Hypki et al. 2022). Indeed, many of the observational studies cited in the Section 1 have found that SG stars are more centrally concentrated than the FG population, in general agreement with what found in this paper and in previous theoretical studies (e.g. D’Ercole et al. 2008; Bekki 2010; Calura et al. 2019).

As for the kinematical properties on which this paper is focused on, a number of observational studies (see the discussion in the Section 1) revealed that in some clusters the SG system rotates more rapidly than the FG, in agreement with what found in our work. The observational study of the kinematics of multiple populations is still in its early stages and much work remains to be done to build a comprehensive observational picture of the kinematical differences between FG and SG stars and, for clusters with multiple SG groups, between the various SG groups. Concerning the possible difference between various SG subgroups, our study predicts that the more extreme SG population is initially rotating more rapidly than the intermediate (moderately enriched) SG groups; this result is consistent with the trend found by Cordero et al. ( 2017) and Kamann et al. ( 2020) in M13 and M80, respectively. As for the morphology of different stellar populations, a few early studies have found that in some clusters the SG subsystem is more flattened than the FG one (see Lee 2017; Cordoni et al. 2020a,b; see also Lee 2018 for a cluster showing instead an opposite trend) providing possible examples of clusters retaining some memory of the initial differences emerging from our simulations. Additional observational studies and numerical simulations exploring the long-term evolution of clusters with initial conditions informed from our study will be necessary to further explore these aspects of the structure and kinematics of multiple stellar populations.

Finally, we conclude this section with a brief remark on the mass of SG stars forming in our simulations. The total mass of SG stars formed by the end of our simulation is about 0.1 the mass of the FG system, in the low-density models, and about equal to the mass of FG stars, in the high-density models. As shown in a number of previous studies (see e.g. D’Ercole et al. 2008; Sollima 2021; Vesperini et al. 2021), the fraction of the total cluster mass in SG stars may significantly increase as a result of the preferential loss of FG stars during the cluster’s early and long-term evolution but, as already pointed out in Calura et al. (2019), the extent of mass-loss required to reach the values of the SG mass fraction observed in present-day globular clusters (see e.g. Milone et al. 2017) is much less extreme than sometimes reported in the literature. More extended surveys of simulations exploring the formation of multiple populations are necessary to test the role of various physical ingredients (e.g. FG structural properties and initial stellar mass function, stellar feedback, external environment) and shed further light on the initial mass of the SG population and its dependence on the initial FG properties (see e.g. McKenzie & Bekki 2021; Yaghoobi et al. 2022 for two recent studies addressing these issues).

5 CONCLUSIONS

An increasing number of observational studies of the kinematics of stars clusters are revealing that internal rotation is a common kinematic feature of these stellar systems. It is therefore important to include rotation in theoretical studies of star clusters and explore its role on their formation and dynamical evolution.

In this paper, by means of 3D hydrodynamic simulations, we have studied the role of rotation on the formation and early dynamics of multiple stellar populations in GCs.

Our models follow the formation of SG stars out of the ejecta released by AGB stars in a rotating FG system and the external pristine gas accreted by the cluster. We have explored the resulting structural and kinematical properties of the SG subsystem, and studied the differences between the dynamical properties of FG stars and those of the various SG subgroups.

In our simulations, we have modelled a massive proto-GC with an FG mass of 10⁷ M_⊙, moving through a uniform gas distribution for which we have assumed two different density values: a low-density one characterized by ρ_pg = 10⁻²⁴ g cm⁻³ and a high-density one with ρ_pg = 10⁻²³ g cm⁻³. In order to explore the interplay between the internal rotational dynamics of the AGB ejecta and that of the external infalling gas, we have considered two different configurations: one in which the FG rotational axis is perpendicular to the motion of the cluster through the external medium and one in which it is parallel to it.

The main results of the paper are the following:

• Our simulations have revealed the complex hydrodynamics/stellar dynamics of the SG formation phase in the presence of a rotational FG system. As derived in previous investigations, we find that the SG forms concentrated in the innermost regions of the FG cluster. Both the SG morphology and kinematics are significantly affected by rotation and the interaction between rotating AGB ejecta and non-rotation infalling external gas. The AGB ejecta initially collect in a disc of gas in the inner regions of the FG cluster and form a rotating disc of helium-enhanced SG stars. The disc survives in the model of a cluster moving in a low-density external medium where the infalling, non-rotating gas has a minor effect on the overall evolution of the system. An inner SG disc also forms in the high-density model but, as a consequence of the earlier arrival of the infalling gas and its higher density, the disc is disrupted before the end of the simulation. Although the long-term dynamical processes may gradually alter the SG disc and drive it towards an increasingly spherical spatial distribution, massive clusters forming a helium-enhanced population in a low-density external environment are those where some memory of an initial flattened SG subsystem might be found.

• The SG populations forming in a rotating FG cluster embedded in low-density pristine gas are characterized by a rotational amplitude larger than that of the FG population. The differences between the SG and the FG rotational kinematics we find in our simulations are generally consistent with the findings of previous theoretical (Bekki 2010, 2011; McKenzie & Bekki 2021) and observational investigations (Cordero et al. 2017; Cordoni et al. 2020a; Dalessandro et al. 2021; Szigeti et al. 2021).

• The more He-rich SG subgroups forming earlier, mainly out of AGB ejecta, rotate around the cluster’s centre more rapidly than SG stars formed later, out of a mix of rotating AGB ejecta and non-rotating infalling pristine gas. However, for a comparison with the present-day properties of old GCs, it is necessary to take into account the effects of the cluster long-term dynamical evolution (Mastrobuono-Battisti & Perets 2013, 2016; Hénault-Brunet et al. 2015; Tiongco et al. 2019; Mastrobuono-Battisti & Perets 2021; Vesperini et al. 2021; Sollima et al. 2022). The findings of our simulations are generally consistent with those of the first observational studies that have explored the kinematics of different SG subgroups and found that the more helium-enhanced SG stars rotate faster than the helium-poor ones (Cordero et al. 2017; Kamann et al. 2020). When SG stars are formed from a rotating FG component embedded in high-density pristine gas, no significant differences are found both between the SG subgroups and the FG and SG systems.

• Very minor differences have been found between models assuming a FG solid-body rotational profile and those assuming the analytical profile of Eq. 6 or changing the orientation of the rotational axis with respect to the direction of the infalling external gas.

• The more complex hydrodynamics of SG formation in a rotating FG cluster does not significantly affect the final distribution of the helium abundances of the SG populations.

In future studies, we will further expand the investigation presented here to fully explore the dependence of our results on the initial properties of the FG clusters (e.g. initial mass, structure, and strength of initial rotation) and those of the external environment (e.g. clumpy ISM).

ACKNOWLEDGEMENTS

We are grateful to the anonymous referee for the useful suggestions. AMB acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 895174. FC acknowledges support from grant PRIN MIUR 2017- 20173ML3WW 001, from the INAF main-stream (1.05.01.86.31) and from PRIN INAF 1.05.01.85.01. EV acknowledges support from NSF grant AST-2009193. We acknowledge PRACE for awarding us access to Discoverer at Sofia Tech Park, Bulgaria. We acknowledge the computing centre of Cineca and INAF, under the coordination of the ‘Accordo Quadro MoU per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione’, for the availability of computing resources and support. We acknowledge the use of computational resources from the parallel computing cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at the Physics and Astronomy Department in Bologna. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU is also supported in part by Lilly Endowment, Inc.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

APPENDIX A: MODEL WITH SOLID BODY ROTATION

Until now, all presented results have been obtained assuming the analytical velocity profile of Eq. 6 for the FG component. For the low-density model, which is the one where the rotation is much stronger, we have performed, for a comparison, an additional run where we have adopted a solid body rotation profile (LDsbz, see Table 2 for the model description). Such velocity profile was earlier adopted in the 3D hydrodynamic simulations of Bekki (2010, 2011) and McKenzie & Bekki (2021).

Even in this model, several clumps are formed at around 10 Myr, and clearly visible also at 27 Myr in Fig. A1, however their stellar surface density is 3 times lower than in the simulation assuming the analytic velocity profile (LDanaz model). In addition, stars in the disc possess a stronger line of sight velocity in the yz plane, as seen comparing it to Fig. 2, a consequence of the higher velocity imparted to FG ejecta in the outskirt in the solid body model. At the end of the simulation, the SG velocity is overall larger than the one obtained for the LDanaz model due to the slightly larger velocity imposed to the AGB ejecta with the solid body profile. Helium-enhanced stars are, in fact, rotating faster in the outskirts, a consequence of the more extended disc (∼12 pc), obtained for the model with a solid body profile. In particular, Fig. A2 shows the rotational amplitudes along the three Cartesian coordinates with an average value of ∼15 km s⁻¹ along both the x and y axes. The z component is instead lower than in the LDanaz model, even though the disc is tilted of about 15°.

Figure A1.

2D maps of the stellar component for the LDsbz model. The reported quantities are as in Fig. 2.

Open in new tabDownload slide

Figure A2.

Stellar rotation profiles of the model LDsbz at 65 Myr. The reported quantities are as in Fig. 4.

Open in new tabDownload slide

It is worth noticing that, even though the FG is set to rotate with a solid body profile, the SG has a pattern much more similar to Eq. 6, with a peak at around the half-mass radius, which is compatible with the one obtained for the LDanaz model.