The present-day globular cluster kinematics of lenticular galaxies from the E-MOSAICS simulations and their relation to the galaxy assembly histories Free

ABSTRACT

We study the present-day rotational velocity (V_rot) and velocity dispersion (σ) profiles of the globular cluster (GC) systems in a sample of 50 lenticular (S0) galaxies from the E-MOSAICS galaxy formation simulations. We find that |$82{{\ \rm per\ cent}}$| of the galaxies have GCs that are rotating along the photometric major axis of the galaxy (aligned), while the remaining |$18{{\ \rm per\ cent}}$| of the galaxies do not (misaligned). This is generally consistent with the observations from the SLUGGS survey. For the aligned galaxies, classified as peaked and outwardly decreasing (⁠|$49{{\ \rm per\ cent}}$|⁠), flat (⁠|$24{{\ \rm per\ cent}}$|⁠), and increasing (⁠|$27{{\ \rm per\ cent}}$|⁠) based on the V_rot/σ profiles out to large radii, we do not find any clear correlation between these present-day V_rot/σ profiles of the GCs and the past merger histories of the S0 galaxies, unlike in previous simulations of galaxy stars. For just over half of the misaligned galaxies, we find that the GC misalignment is the result of a major merger within the last |$10\, \mathrm{Gyr}$| so that the ex-situ GCs are misaligned by an angle between 0° (co-rotation) and 180° (counter-rotation), with respect to the in situ GCs, depending on the orbital configuration of the merging galaxies. For the remaining misaligned galaxies, we suggest that the in situ metal-poor GCs, formed at early times, have undergone more frequent kinematic perturbations than the in situ metal-rich GCs. We also find that the GCs accreted early and the in situ GCs are predominantly located within 0.2 virial radii (R₂₀₀) from the centre of galaxies in 3D phase-space diagrams.

1 INTRODUCTION

Lenticular (S0) galaxies are characterized by a bulge surrounded by a smooth disc with very little or no ongoing star-formation. The fraction of lenticular galaxies is found to increase within the high-density environments of rich galaxy clusters at the expense of spiral galaxies (Dressler 1980; Dressler et al. 1997; Fasano et al. 2000). For this reason, understanding how S0 galaxies form is important for constraining which physical processes play a dominant role in the transformation and quenching of galaxies in different environments.

Focusing on the kinematic properties of galaxies, previous studies have found that S0 galaxies in the field and small galaxy groups are more pressure-supported than S0 galaxies in clusters (Coccato et al. 2020; Deeley et al. 2020). Additionally, S0 galaxies in the field and small galaxy groups also showed more frequently misaligned stellar and gas kinematics than cluster S0s (Deeley et al. 2020). These results suggest that cluster S0 galaxies are consistent with spiral progenitors, i.e. faded spirals, which had their star formation quenched and their spiral arm structure suppressed as a result of the interaction with the environment, e.g. ram-pressure stripping, tidal interactions, and starvation (Gunn & Gott 1972; Larson, Tinsley & Caldwell 1980; Bekki 2009; Bekki & Couch 2011; Merluzzi et al. 2016). On the other hand, S0 galaxies in the field and small galaxy groups are likely the result of more complex formation histories that involved multiple mergers and accretion events (e.g. Tapia et al. 2017; Eliche-Moral et al. 2018; Dolfi et al. 2020, 2021).

Previous works have used simulations to study the formation of S0 galaxies through mergers and investigate how mergers influence their kinematic properties. Results showed that both major (mass-ratio >1:4) and minor (1:10< mass-ratio <1:4) mergers can produce S0 galaxies (Bournaud, Jog & Combes 2005; Naab et al. 2014; Wu et al. 2014; Tapia et al. 2017; Eliche-Moral et al. 2018; Schulze et al. 2020) and, specifically, major mergers are expected to spin-up the rotation of the stars out to large radii, while minor mergers are expected to decrease it (Naab et al. 2014; Wu et al. 2014; Schulze et al. 2020).

The commonly proposed two-phase formation scenario for massive early-type galaxies (ETGs) suggests that the galaxies had an early (⁠|$z$| ≥ 2) in-situ formation during which they formed most of their stellar mass, followed by the late (⁠|$z$| ≃ 0) accretion of dwarf galaxies in mini mergers (mass-ratio < 1:10) on to their outskirts (Oser et al. 2010; Damjanov et al. 2014; Zolotov et al. 2015; Rodriguez-Gomez et al. 2016). These late mini mergers are expected to decrease the rotation of the stars at large radii, where the galaxy transitions from a fast-rotating disc to a slowly rotating spheroidal component (Arnold et al. 2011, 2014; Guérou et al. 2016; Bellstedt et al. 2017; Dolfi et al. 2020).

These results suggest that it is important to investigate the full kinematic behaviour of the galaxies from the inner regions out to the outskirts, in order to identify the presence of transitions that can help us to constrain the galaxy-specific merger and accretion histories.

In a more recent work, Schulze et al. (2020) investigated the kinematic profiles of stars out to ∼5 effective radii (R_e) of a sample of ETGs with log (M_*/M_⊙) ≥ 10.3 from the Magneticum simulations. From the ratio of the rotation velocity to velocity dispersion profiles (V_rot/σ), they identified three distinct kinematic profile shapes: a peaked and outwardly decreasing (peaked, hereafter), an increasing, and a flat V_rot/σ profile out to |$\sim 5\, R_{\mathrm{e}}$| that were likely the result of the different merger histories of the galaxies. Specifically, they suggested that galaxies with peaked V_rot/σ profiles had an assembly history characterized by the late accretion of dwarf galaxies in minor and mini mergers that were disrupted beyond |$\sim 2\, R_{\mathrm{e}}$|⁠, enhancing the random motion of the galaxy outskirts without disrupting its central disc-like kinematics. The galaxies with flat and increasing V_rot/σ profiles were more likely to have experienced a late (i.e. |$z$| ≲ 1) major merger event. This event was likely more gas-rich for galaxies with increasing V_rot/σ profiles, leading to the re-formation of a disc component and the rising of the rotational velocity of the stars out to large radii. On the other hand, the late major merger event did not contribute to a significant gas fraction in the galaxies with flat V_rot/σ profiles and likely disrupted the central disc-like kinematics of the galaxy. Therefore, we expect that the present-day stellar V_rot/σ profiles of ETGs can be used to constrain their past assembly history.

In Dolfi et al. (2021), we studied the kinematic properties of a selected sample of nine S0 galaxies from the SAGES Legacy Unifying Globulars and GalaxieS (SLUGGS; Brodie et al. 2014) survey extending out to ∼4–|$6\, R_{\mathrm{e}}$|⁠. To reach out to |$\sim 5\, R_{\mathrm{e}}$|⁠, we combined the stellar kinematics (Arnold et al. 2014; Foster et al. 2016) with those of the globular clusters (GCs; Forbes et al. 2017a) and planetary nebulae (PNe; Pulsoni et al. 2018), which can be more easily probe out to large radii due to their brightness. Both the GCs and PNe are expected to trace the underlying stellar population of the ETGs. In fact, the PNe are stars with initial masses |$1\, \mathrm{M_{\odot }} \lt \mathrm{M}_{*}\lt 8\, \mathrm{M_{\odot }}$| that have evolved past the main sequence towards the red giant phase and should, thus, trace the kinematics of the stars out to large radii (e.g. Buzzoni, Arnaboldi & Corradi 2006; Romanowsky 2006). On the other hand, the GCs typically display a colour bimodality that reflects the two distinct sub-populations: the red, metal-rich (MR) and blue, metal-poor (MP) GCs. These two GC sub-populations are not necessarily equal to the in-situ and ex-situ (accreted) GCs, respectively (Forbes et al. 2018; Kruijssen et al. 2019), however they are expected to trace the galaxy kinematics differently (Trujillo-Gomez et al. 2021). The red, MR GCs are expected to trace the kinematics of the bulge and spheroid in ETGs, since they largely formed in-situ at roughly the same time as the bulk of the host galaxy stars. The blue, MP GCs are expected to trace the kinematics of the stellar halo of the galaxy, since they likely formed at earlier times than the red, MR GCs or were later accreted on to the galaxy (Forbes 1997; Strader et al. 2005; Brodie & Strader 2006). In a recent work, using the E-MOSAICS simulations, Reina-Campos et al. (2021) found that the MP GCs typically have more extended and shallower surface density profiles than the MR GCs, suggesting an accreted origin for a significant fraction of the former. The surface density profile of the MP GCs is also found to correlate with the slope of the dark matter halo density profile, therefore suggesting that the GC sub-populations are excellent probes for studying the outskirts of the galaxies. Additionally, simulations have also shown that major mergers are expected to spin-up the rotation of the GC population at large radii, as opposed to minor mergers that are expected to decrease it (Bekki et al. 2005).

From the study in Dolfi et al. (2021), we found that six of the nine galaxies had consistent stellar, GCs, and PNe kinematics with their rotation occurring along the photometric major-axis of the galaxy (aligned galaxies). In the remaining three galaxies, the PNe and GCs were rotating along a kinematic axis misaligned with respect to that of the stars, which was consistent with the photometric major-axis of the galaxy (misaligned galaxies). These results seem to suggest distinct assembly histories for the aligned and misaligned galaxies, with the latter having likely undergone more recent and multiple accretion events. Among the six aligned galaxies, we also found that four galaxies show a peaked V_rot/σ profile, while the remaining two galaxies show a flat V_rot/σ profile. Therefore, the comparison with the simulations of Schulze et al. (2020) would suggest an assembly from late minor and mini mergers for the peaked galaxies, while the flat galaxies were likely involved in a merger with higher mass-ratio.

However, these results assume that both GCs and PNe are tracing the kinematics of the underlying stellar population out to large radii. We find that this seems to be the case for the six aligned galaxies, but not for the remaining three misaligned galaxies that show twists in the kinematic position angle as a function of the galactocentric radius (see fig. 3 in Dolfi et al. 2021).

In light of this, the aim of this paper is to use the E-MOSAICS¹ simulations (Pfeffer et al. 2018; Kruijssen et al. 2019) to study the present-day kinematic properties of the GC systems of the S0 galaxies to understand how they are related to the past formation histories of the S0 galaxies. In a recent work, Deeley et al. (2021) have found that more than half (i.e. |$57{{\ \rm per\ cent}}$|⁠) of the S0 galaxies from the IllustrisTNG simulations have formed through significant merger events, while the remaining |$37{{\ \rm per\ cent}}$| have formed through gas stripping events as a result of the infall on to a galaxy group or cluster. Therefore, as one of the most dominant S0 formation pathways, it is important to further understand how the different merger events affect the present-day kinematic properties of S0 galaxies out to large galactocentric radii and, specifically, the kinematic properties of their GC systems. We also aim to further investigate the evolution of the aligned and misaligned galaxies to understand what physical processes are responsible for the present-day GC misalignment that we see in the observations from the SLUGGS survey (Dolfi et al. 2021).

The structure of the paper is as follows: In Section 2 we describe the E-MOSAICS simulations. In Section 3 we describe the selection criteria of the simulated S0 galaxy sample, as well as of the GC systems of each galaxy. In Section 4 we describe the kinematic analysis and results of the GC systems of our selected S0 galaxies, and we also investigate the kinematics of the MR and MP GCs separately for a sub-sample of the galaxies. In Section 5 we study whether there exists a correlation between the location of the GC systems of the galaxies on the 3D phase-space diagrams and their accretion time on to the galaxy. In Section 6 we discuss the assembly histories of our simulated galaxies by investigating the relation between the V_rot/σ profiles of the GCs with the merger histories of the galaxies derived from the E-MOSAICS simulations. In Section 7 we summarize our results and conclusions.

2 THE E-MOSAICS SIMULATIONS

E-MOSAICS is a suite of models designed to study the formation and evolution of star clusters and of their host galaxies (Pfeffer et al. 2018). E-MOSAICS is obtained from a combination of the MOSAICS (Kruijssen & Lamers 2008; Kruijssen & Mieske 2009; Kruijssen et al. 2011; Pfeffer et al. 2018) semi-analytic model for the star cluster formation and evolution, and the suite of cosmological, hydrodynamical simulations from EAGLE (Crain et al. 2015; Schaye et al. 2015) for the formation and evolution of galaxies in a standard Λ cold dark matter Universe.

EAGLE is run using a modified version of the N-Body Tree-PM smoothed particle hydrodynamics (SPH) code, GADGET3 (Springel et al. 2005), which includes modifications to the SPH code, time stepping scheme and sub-grid physics models implemented in the simulations. Some of the main relevant physical processes that are implemented in the simulations within these sub-grid models include radiative cooling, reionization, star formation, stellar mass loss and metal enrichment, feedback from star formation and active galactic nuclei, and black hole growth from gas accretion and mergers. The efficiency of the sub-grid models for feedback are calibrated in order to reproduce certain observables. Specifically, the efficiency from the AGN feedback is calibrated to ensure that the simulations are reproducing the observed present-day relation between the central black hole mass and galaxy stellar mass. On the other hand, the efficiency from stellar feedback is calibrated to reproduce the observed present-day galaxy stellar mass-to-halo mass relation (Schaye et al. 2015).

MOSAICS is a semi-analytic model for the formation and evolution of star clusters that is incorporated into the EAGLE simulations of galaxy formation as a sub-grid model (Kruijssen et al. 2011). Within this model, star clusters are not simulated as individual particles, but they are formed as a sub-grid component of a star particle that is produced as a result of a gravitationally unstable gas particle in the simulations. Therefore, a fraction of the mass of the newborn star particle (which is set by the cluster formation efficiency; Bastian 2008) is used to generate a sub-grid population of star clusters, of which the initial properties are derived from the properties of the local gas and star particle at the moment of their formation (Kruijssen 2012; Reina-Campos & Kruijssen 2017). Subsequently to their formation, star clusters evolve and undergo mass-loss as a result of stellar evolution, modelled following the implementation for the star particles in EAGLE (Wiersma et al. 2009), two-body relaxation (Lamers et al. 2005), and tidal shocks from the interstellar medium (ISM). Dynamical friction is also another important process that may be responsible for the mass-loss and disruption of star clusters, and it is included in post-processing in E-MOSAICS (Pfeffer et al. 2018). The time-scales of dynamical friction are calculated for each star cluster according to the definition by Lacey & Cole (1993) and star clusters are assumed to be completely disrupted by dynamical friction if their ages exceed the time-scale for dynamical friction. Dynamical friction is most effective in the disruption of the most massive star clusters and at small galactocentric radii.

The present-day (i.e. |$z$| = 0) stellar masses of the star clusters formed within the simulations can range between 10² < M_*/M_⊙ < 10⁸ (Pfeffer et al. 2018). However, the lowest mass star clusters with |$\mathrm{M_{*}} \lt 5 \times 10^{3}\, \mathrm{M_{\odot }}$| are immediately removed after their formation. On the other hand, the star clusters with |$\mathrm{M_{*}} \sim 5\times 10^{4}\, \mathrm{M_{\odot }}$| are typically disrupted due to tidal shocks over the course of a Hubble time, while the high-mass star clusters (i.e. |$\mathrm{M_{*}} \gt 10^{5}\, \mathrm{M_{\odot }}$|⁠) are removed by dynamical friction. The GC mass function in E-MOSAICS is found to agree well with that of observed galaxies at the high-mass end (i.e. |$\mathrm{M_{*}} \gt 10^{5}\, \mathrm{M_{\odot }}$|⁠). However, the simulations overproduce the number of star clusters with |$\mathrm{M_{*}} \le 10^{5}\, \mathrm{M_{\odot }}$|⁠. Indeed, star clusters with |$\mathrm{M_{*}} \sim 10^{4}\, \mathrm{M_{\odot }}$| are overestimated by a factor of ∼10–100, while star clusters with |$\mathrm{M_{*}} \sim 10^{5}\, \mathrm{M_{\odot }}$| are overestimated by a factor of 2–10 in E-MOSAICS (Pfeffer et al. 2018). This issue may be traced to the fact that EAGLE does not incorporate an explicit model of the cold ISM and, as a result, low-mass star clusters are not efficiently disrupted through tidal shocks from the interactions with the ISM.

The inefficiency of the simulations at disrupting star clusters, due to the lack of an explicit cold ISM model in EAGLE, is also responsible for the excess of MR GCs in E-MOSAICS. Specifically, the simulations overproduce the number of MR GCs by a factor of 5 in the metallicity range −0.5 < [Fe/H] < 0.0 and by a factor of 2.5 in the metallicity range −1.0 < [Fe/H] < −0.5 (Kruijssen et al. 2019). However, Kruijssen et al. (2019) found that the overproduction of the MR GCs in the metallicity range −1.0 < [Fe/H] < −0.5 in the simulations does not have a statistically significant impact on the derived properties and correlations describing the formation and evolution of the galaxies.

In Section 3.2 we will discuss the adopted selection criteria to minimize the impact of the overproduction of the low-mass star clusters in our kinematic results.

Overall, the E-MOSAICS simulations are found to be able to reproduce most of the key observed properties of both young and old star cluster populations in their host galaxies. For this reason, they have been used in previous works to reconstruct the merger and accretion history of the Milky Way by comparing the age, metallicity, and orbital properties of the observed GCs with those from the simulations (e.g. Kruijssen 2019; Kruijssen et al. 2020; Pfeffer et al. 2020; Trujillo-Gomez et al. 2021).

In addition to the initial suite of 25 zoom-in simulations of Milky Way-mass galaxies, the cosmological volume of the E-MOSAICS simulations spans |$34.4\, \mathrm{cMpc}$| (Crain et al., in preparation; first reported in section 2.2 of Bastian et al. 2020) and the simulations have an initial gas particle mass of |$2.26\times 10^{5}\, \mathrm{M}_{\odot }$|⁠, thus they can resolve galaxies with |$\mathrm{M}_{*} \gt 10^{7}\, \mathrm{M_{\odot }}$| with more than 100 stellar particles. However, E-MOSAICS includes very few massive galaxies with |$\mathrm{M}_{*} \ge 10^{11}\, \mathrm{M_{\odot }}$| that are mainly ellipticals (Correa et al. 2017). At the same time, the largest galaxy cluster included in the simulations has halo mass |$\mathrm{M_{200}} \sim 6 \times 10^{13}\, \mathrm{M_{\odot }}$|⁠, i.e. Fornax cluster-like, meaning that E-MOSAICS does not probe the dense environments of rich galaxy clusters, e.g. Coma cluster-like. Therefore, E-MOSAICS is limited to the low-density environments of the field and galaxy groups.

3 THE DATA

3.1 The S0 galaxy sample

From the simulations, we initially select 170 galaxies with total stellar masses |$\mathrm{M_{*}}\gt 10^{10}\, \mathrm{M_{\odot }}$|⁠. This lower limit on the total galaxy stellar mass is chosen to match the stellar mass distribution of the observed ETGs from the SLUGGS survey (Brodie et al. 2014; Forbes et al. 2017b).

We produce mock SDSS g-band images of the edge-on and face-on view of each galaxy, colour-coded by the surface brightness calculated for each stellar particle using the fsps stellar population model (Conroy, Gunn & White 2009; Conroy & Gunn 2010). The images are modified to enhance the spiral arms so that we could more easily identify the S0 galaxies for our final sample. Therefore, following Trayford et al. (2017), young stars (i.e. |$\lt 100\, \mathrm{Myr}$| old) are re-sampled from the star-forming gas particles at higher resolution (i.e. |$1000\, \mathrm{M_{\odot }}$|⁠), assuming a constant star-formation rate (SFR), in order to better trace the spiral structure of galaxies. The synthetic images were generated using SPHviewer (Benitez-Llambay 2015).

From the visual inspection of the edge-on and face-on images, which was performed by one of the authors of this paper Warrick J. Couch (WJC), we identify 51 S0 galaxies. The S0 galaxies are selected based on the identification of a disc component and absence of the spiral arms from the edge-on and face-on view images. However, despite the best effort to identify a representative sample of S0 galaxies, we note that there may still be few mis-classified S0s in our sample that may be prolate ellipticals. Fig. 1 shows the edge-on and face-on view of one example S0 galaxy from the E-MOSAICS simulations, which is included in our final sample. We also show the corresponding edge-on and face-on spatial distributions of the GC system of the galaxy colour-coded by the line-of-sight velocity after applying the selection criteria described in Section 3.2.1–3.2.4. We note the presence of luminous clumps around the galaxy in Fig. 1. However, these very disperse clumps are not necessarily a sign of spiral arm structure, but they are a result of the smoothing adopted to enhance the youngest stars with respect to the old stars.

Next, we consider the specific star formation rate (sSFR) of the 51 visually selected S0 galaxies, and we include in our final sample only those galaxies with sSFR more than |$0.3\, \mathrm{dex}$| below the star-forming main sequence (SFMS). This constraint should ensure that we are reliably selecting a population of quenched galaxies, if the scatter of the SFMS is |$\sim 0.3\, \mathrm{dex}$|⁠, as found in previous studies (e.g. Speagle et al. 2014). Fig. 2 shows the distribution of our S0 galaxies from the E-MOSAICS simulations on the sSFR-M_* plane (blue dots), with the best-fitting line of the location of the SFMS (solid black line) from Renzini & Peng (2015). We note that in the EAGLE simulations the |$\mathrm{sSFR}\simeq 10^{-10}\, \mathrm{yr}^{-1}$| for the SFMS at |$\mathrm{M_{*}}\simeq 10^{10}\, \mathrm{M_{\odot }}$| (see fig. 5 in Furlong et al. 2015), which is consistent with the value obtained from the SFMS relation of Renzini & Peng (2015) for the same stellar mass, as shown in Fig. 2. We find that only one galaxy fails to satisfy our selection criteria based on the sSFR, as it lies above the SFMS (red square in Fig. 2). Therefore, we exclude this galaxy from our final sample. We note that there are 11 galaxies in the simulations with very low SFRs that are practically equal to zero. In Fig. 2, we plot these galaxies with a fixed |$\mathrm{sSFR}=10^{-14}\, \mathrm{yr}^{-1}$|⁠. After the sSFR cut, we obtain a final sample of 50 S0 galaxies, with quenched star-formation, from the E-MOSAICS simulations.

In Fig. 3, we show the distribution of our final sample of 50 simulated S0 galaxies in the half-light radius (R_e)-stellar mass (M_*) plane (blue dots). Here, we assume that mass follows light in the simulations and, therefore, that the 2D half-mass radius of our E-MOSAICS S0 galaxies is equal to the 2D half-light radius, R_e. For comparison, we also show the nine S0 galaxies (red stars), observed from the SLUGGS survey (Brodie et al. 2014), that we have studied in our previous works (Dolfi et al. 2020, 2021). For the observed galaxies, the circularized radius (R_e) and the total stellar mass (M_*) are taken from Forbes et al. (2017b). The solid line represents the R_e-M_* single power-law relation of ETGs from Lange et al. (2015), which is obtained using the best-fitting parameters in the K-band for the ETGs with |$\mathrm{M_{*}} \gt 2 \times 10^{10}\, \mathrm{M_{\odot }}$| (see table B2 in Lange et al. 2015).

In Fig. 3, we note that there is quite a bit of scatter of both our simulated and observed S0 galaxies around the K-band single power relation of ETGs with |$\mathrm{M_{*}} \gt 2 \times 10^{10}\, \mathrm{M_{\odot }}$| from Lange et al. (2015). Additionally, there is also one simulated S0 with |$\mathrm{M_{*}} \gt 2 \times 10^{11}\, \mathrm{M_{\odot }}$| characterized by very large |$R_{\mathrm{e}}\sim 20\, \mathrm{kpc}$|⁠. In Lange et al. (2015), the authors do not report the 1σ scatter of the observed data around the R_e-M_* single power relation of ETGs. However, we note that the ETGs from Lange et al. (2015) show an overall similar scatter around the corresponding K-band single power-law relation, with ETGs having R_e ranging between ∼1 and |$3\, \mathrm{kpc}$| at |$\mathrm{M_{*}} = 10^{10}\, \mathrm{M_{\odot }}$| and between ∼2 and |$10\, \mathrm{kpc}$| at |$\mathrm{M_{*}} = 10^{11}\, \mathrm{M_{\odot }}$|⁠. At |$\mathrm{M_{*}} \gt 10^{11}\, \mathrm{M_{\odot }}$|⁠, some ETGs also show half-light radii as large as |$R_{\mathrm{e}}\sim 20\, \mathrm{kpc}$|⁠. This scatter in the R_e–M_* plane of ETGs slightly depends on the imaging band from which the R_e of the galaxies was calculated, being the largest in the u-band (Lange et al. 2015). In a recent work, de Graaff et al. (2021) produced mock r-band images of galaxies in the EAGLE simulations to compare the mass–size relation of simulated and observed galaxies. They found that using the r-band sizes, measured for the simulated galaxies, improves the agreement between the simulated and observed mass–size relation of quiescent and star-forming galaxies. Therefore, we conclude that our simulated and observed S0 galaxies are, overall, consistent with the R_e–M_* relation from Lange et al. (2015) in Fig. 3.

Fig. 3 also shows that our simulated S0 galaxies cover overall a similar stellar mass range as our observed S0 galaxies from the SLUGGS survey. However, we note a bias of the S0 galaxies from the SLUGGS survey to higher stellar masses (i.e. |$\gtrsim 10^{11}\, \mathrm{M_{\odot }}$|⁠) than the S0 galaxies from the E-MOSAICS simulations. This is due to the fact that the EAGLE simulations do not contain many high-mass galaxies with |$\mathrm{M_{*}}\gtrsim 10^{11}\, \mathrm{M_{\odot }}$| and that the majority of these massive galaxies in EAGLE are ellipticals (Correa et al. 2017). Finally, due to the limited volume of the E-MOSAICS simulations, the majority of the simulated galaxies are located in low-density environments (i.e. field and galaxy groups; see Section 2), similarly to the observed S0 galaxies from the SLUGGS survey. Therefore, the consistency between the environment of the S0 galaxies from the SLUGGS survey and E-MOSAICS simulations allow us to make comparisons between our simulated and observed S0 galaxy sample, when investigating how their kinematic properties relate to the merger histories of the galaxies. Additionally, in the next Section 4.1, we will consider the high-mass (i.e. 10.5 < log (M_*/M_⊙) < 11.6) S0 galaxy sub-sample in the simulations for the comparison with the SLUGGS survey.

3.2 The GC systems of the S0 galaxies in the E-MOSAICS simulations

3.2.1 GC selection: stellar mass cut

For each of the 50 S0 galaxies from the E-MOSAICS simulations in the previous Section 3.1, we select all of their star clusters that are gravitationally bound to the host galaxy and have stellar masses |$\mathrm{M_{*}}\gt 10^{5}\, \mathrm{M_{\odot }}$| at |$z$| = 0. This stellar mass cut ensures that our kinematic results will not be influenced by the overabundant population of the low-mass (i.e. |$\mathrm{M_{*}}\lesssim 10^{5}\, \mathrm{M_{\odot }}$|⁠) star clusters that were not destroyed from the interaction with the cold ISM (see Section 2). The same stellar mass cut (i.e. |$\mathrm{M_{*}}\gt 10^{5}\, \mathrm{M_{\odot }}$|⁠) was also used by Kruijssen et al. (2019), as they found that the GC mass function of their sample of simulated galaxies was consistent with that of the Milky Way after removing the lower mass GCs.

3.2.2 GC selection: magnitude cut

Fig. 4 shows the V-band GC luminosity function of one S0 galaxy from our selected sample, as an example. The grey histogram shows the V-band GC luminosity function prior to applying any stellar mass cut to the GC system, while the black histogram is obtained after removing the low-mass GCs with |$\mathrm{M_{*}}\lesssim 10^{5}\, \mathrm{M_{\odot }}$|⁠. We notice that the V-band GC luminosity function resembles more the shape of a Gaussian with a luminosity peak at M_V ≃ −7 when selecting only the GCs with |$\mathrm{M_{*}}\gt 10^{5}\, \mathrm{M_{\odot }}$| for typically all the GC systems of the S0 galaxies in our sample. This is more consistent with the V-band GC luminosity function of observed galaxies, which is characterized by a Gaussian-like shape with luminosity peak typically located at M_V ≃ −7.5 (e.g. Forbes et al. 2018). The non Gaussian-like shape with prominent peak below M_V ≃ −7 of the GC systems of our S0 galaxies seems likely to be the result of the overabundant population of the low-mass (⁠|$\mathrm{M_{*}}\lesssim 10^{5}\, \mathrm{M_{\odot }}$|⁠) star clusters that are not being efficiently disrupted in the simulations (see Section 2), as shown in Fig. 4 prior to applying any GC stellar mass cut. For this reason, we exclude these low-mass star clusters in order to avoid significant contamination from these objects. The red vertical-dashed line shows the V-band magnitude of the brightest GC that is included in the system (i.e. M_V = −10.5). This luminosity limit is chosen to include up to the V-band magnitude of ω-Centauri (M_V ≃ −10.3), which is the most luminous GC of the Milky Way, consistent with the observations of the GC systems from the SLUGGS survey (Pota et al. 2013).

3.2.3 GC selection: age cut

Since, in this work, we are mainly interested in studying the kinematic properties of the GCs (i.e. old star clusters), we then select only those star clusters with ages above |$8\, \mathrm{Gyr}$| in each galaxy. In the Milky Way, an age cut above |$8\, \mathrm{Gyr}$| would only exclude three of its GCs, as the bulk of its GCs are older than |$\sim 10\, \mathrm{Gyr}$| (see fig. 1 in Forbes 2020). While galaxies may have a larger spread in the ages of their GCs with some being |$\lt 5\, \mathrm{Gyr}$| old (see NGC 3377 in Usher et al. 2019), the bulk of the GCs has typically old ages (⁠|$\gt 8\, \mathrm{Gyr}$|⁠) in the remaining galaxies (Usher et al. 2019). Overall, we find that the fraction of the star clusters with ages below |$8\, \mathrm{Gyr}$| accounts for |$\lt 30{{\ \rm per\ cent}}$| of the total GC population in |$62{{\ \rm per\ cent}}$| of our simulated S0 galaxies. Therefore, the majority of our galaxies are dominated by an old GC population which is the main focus of this work.

3.2.4 GC selection: metallicity cut

Following Kruijssen et al. (2019), we remove the most MP GCs with [Fe/H] < −2.5. Specifically, Kruijssen et al. (2019) restricted the analysis to the GC metallicity range −2.5 < [Fe/H] < −0.5 to mimic the metallicity distributions spanned by the observed GCs of the Milky Way and M31, while reducing the contamination from the overproduced MR GCs (see Section 2). In this work, the adopted low metallicity cut is also observationally motivated. In fact, previous spectroscopy analyses have measured the metallicity of the GCs of the observed ETGs from the SLUGGS survey and found that all the GCs had metallicity [Fe/H] > −2.5 (Pastorello et al. 2015; Usher et al. 2019). Additionally, we find that the MP GCs with [Fe/H] < −2.5 are very low in numbers accounting for less than |$10{{\ \rm per\ cent}}$| of the total GC population in each one of our simulated S0 galaxies, in accordance with the results discussed in the previous works of Pfeffer et al. (2018) and Kruijssen et al. (2019). In fact, previous works also predict the existence of a metallicity floor at [Fe/H] = −2.5 as a result of the galaxy mass–metallicity relation, below which high-redshift low-mass galaxies could not form GCs that were massive enough (i.e. |$\mathrm{M_{*}}\gt 10^{5}\, \mathrm{M_{\odot }}$|⁠) to survive until the present-day (Beasley et al. 2019; Kruijssen 2019). Additionally, these low-mass galaxies are also not resolved in the E-MOSAICS simulations, since they are expected to have stellar mass |$\mathrm{M_{*}}\sim 10^{5}\, \mathrm{M_{\odot }}$| at [Fe/H] = −3 (Kruijssen 2019), which is the particle resolution of the simulations. For this reason, MP stars with [Fe/H] < −2.5 are not reliable as they are not modelled in the simulations. However, these results do not rule out the existence of GCs with metallicity [Fe/H] < −2.5, as suggested by the discovery of a very MP GC with [Fe/H] = −2.9 in M31 (Larsen et al. 2020). Therefore, while GCs with [Fe/H] < −2.5 may likely exist in the local Universe, they are expected to be rare and low in numbers as seen in the simulations. For this reason, their removal should not significantly influence the results in this work.

On the other hand, we do not apply any cut to the overproduced GC population in the metallicity range −1.0 < [Fe/H] < 0.0, since our tests show that the inclusion of this GC population does not significantly impact the kinematic properties of the GC systems of our simulated S0 galaxies. Additionally, recent results have found that the high metallicity GCs are not overproduced in massive galaxies (i.e. M_* ≳ 10¹¹M_⊙), but mainly in galaxies with stellar masses between |$10^{9.5}\, \mathrm{M_\odot } \lt \mathrm{M_{*}} \lt 10^{10.5}\mathrm{M_{\odot }}$| (Pfeffer et al., in preparation). Overall, the metallicity range −2.5 < [Fe/H] < 0.5 of the GC systems of our sample of selected S0 galaxies from the E-MOSAICS simulations is consistent with the metallicity measured for the GC systems of the observed ETGs from the SLUGGS survey (e.g. Pastorello et al. 2015; Usher et al. 2019), which we have studied in our previous work (Dolfi et al. 2021).

3.2.5 GC selection: surface density cut

We find that a very small fraction of our S0 galaxies are characterized by very extended GC systems that reach out to galactocentric radii as large as ∼0.5–|$1\, \mathrm{Mpc}$|⁠. Most of these galaxies also show evidence of shells or tidal tails from the visual inspection of their images with Mpc-scale field-of-views (FoV), suggesting that they may have recently experienced interactions. Therefore, we apply an observationally motivated cut to the surface density profile of the GC systems of our S0 galaxies, in order to avoid contamination from potential outliers GCs that may belong to the outer regions of the satellite galaxies orbiting around the central galaxy. Pota et al. (2013) found that the surface density profiles of the red and blue GCs of the ETGs from the SLUGGS survey reach levels as low as log (Σ_GCs/arcmin²) ≃ −2, corresponding to log (Σ_GCs/kpc²) ≃ −2. After we apply the same density cut to the GC surface density profiles of all our simulated S0 galaxies, the GC systems extend out to galactocentric radii ranging between ∼10 and |$200\, \mathrm{kpc}$|⁠. This is overall consistent with the radial extension of the observed S0 galaxies from the SLUGGS survey, with the most extended galaxy reaching out to |$\sim 100\, \mathrm{kpc}$|

4 KINEMATIC ANALYSIS AND RESULTS

In this section, we calculate the 1D kinematic profiles of the GC systems of our selected sample of 50 S0 galaxies from the E-MOSAICS simulations after applying the GC selection criteria described in Section 3.2. Here, all the galaxies are rotated such that the angular momentum vector of the stars aligns with the |$z$|-axis (i.e. the disc lies on the xy-plane). We carry out the kinematic analysis focusing on the edge-on projection of the galaxy. We note that the viewing angle may have an influence for the comparison with the observations from the SLUGGS survey that are characterized by random viewing angles, as well as for the comparison with other simulations.

To calculate the 1D kinematic profiles of the GC systems of our simulated S0 galaxies, we use the kinemetry² method developed by Krajnović et al. (2011). Within this method, we bin our data in elliptical annuli centred on the galaxy and we recover the best-fitting moments of the line-of-sight-velocity-distribution (LOSVD) along each annulus. The best-fitting LOSVD moments are calculated by performing a least-square minimization between the disc model, described by equations (3 and 4) shown in Dolfi et al. (2020), and our kinematic data. The details of the application of this method to sparse and non-homogeneous kinematic data sets are given in Proctor et al. (2009), Bellstedt et al. (2017), and Dolfi et al. (2020), while the details describing the χ²-minimization are given in Foster et al. (2011) and Foster et al. (2016).

As previously done in Dolfi et al. (2020, 2021), we calculate the first two moments of the LOSVD of the GC systems of our S0 galaxies, i.e. rotation velocity (V_rot) and velocity dispersion (σ) profiles, from which we derive the corresponding V_rot/σ profiles. The velocity dispersion of the GCs is calculated by performing a nearest neighbour binning, as described in our previous work (Dolfi et al. 2020). We also fit for both the kinematic position angle (PA_kin) and kinematic axial-ratio (q_kin), where they are well constrained and do not excessively vary as a function of the galactocentric radius. If this is not the case, then we fix the PA_kin to the photometric position angle (PA_phot)³ of the galaxy and we fix the q_kin to the mean value obtained from the unconstrained fit. The 1σ errors on the PA_kin and q_kin profiles are calculated by running kinemetry on 100 bootstrapped samples obtained by sampling with replacement of the original data set. On the other hand, the 1σ errors on the V_rot and 1σ profiles are the standard errors of the mean, calculated from the standard deviation of the velocity and velocity dispersion measurements of the data points in each elliptical annulus divided by the square root of the number of points in each bin.

4.1 1D kinematic profiles of the GC systems of the S0 galaxies in the E-MOSAICS simulations

Fig. 5 shows the 1D kinematic profiles of the GC systems of our selected sample of 50 S0 galaxies from the E-MOSAICS simulations. Since the colour-magnitude diagrams of the GC systems do not show evidence of any colour bimodality, as it was the case for the S0 galaxies from the SLUGGS survey (Pota et al. 2013), we do not split here between the red and blue GC sub-populations. However, in the next Section 4.2, we will investigate whether the MR and MP GC sub-populations show different kinematic behaviours by adopting the same fixed metallicity cut at [Fe/H] = −1 for all the GC systems of our S0 galaxies.

We classify the galaxies as aligned (Section 4.1.1) if the GCs are rotating along the photometric major axis (PA_phot) of the galaxy within the 1σ errors. Otherwise, we classify the galaxies as misaligned (Section 4.1.2).

Finally, we divide both aligned and misaligned galaxies into two stellar mass bins: low-mass 10 < log (M_*/M_⊙) < 10.5 and high-mass 10.5 < log (M_*/M_⊙) < 11.6. The high stellar mass bin is consistent with the stellar mass range covered by all the observed S0 galaxies from the SLUGGS survey that we have studied in Dolfi et al. (2021), with the exception of one very low-mass S0 (i.e. NGC 7457).

4.1.1 Aligned galaxies

Following the simulations from Schulze et al. (2020), we classify the V_rot/σ profiles of the GCs of the aligned galaxies (see Fig. 5, top) as peaked, flat, or increasing based on the V_rot/σ gradient calculated between the outer and inner regions of the GC V_rot/σ profiles. Specifically, we calculate the V_rot/σ gradient (i.e. ΔV_rot/σ) as the difference between the mean value of the V_rot/σ in the outer (i.e. 2.0 < R/R_e < 3.5) and inner (i.e 0.5 < R/R_e < 2.0) radial bins. We classify the galaxy as peaked, flat, or increasing if the V_rot/σ gradient is ΔV_rot/σ < −0.04, −0.04 < ΔV_rot/σ < 0.04 or ΔV_rot/σ > 0.04, respectively, as defined by Schulze et al. (2020).

However, we note here that the stellar V_rot/σ profiles of Schulze et al. (2020) do not extend beyond |$\sim 5\, R_{\mathrm{e}}$|⁠. On the other hand, the V_rot/σ profiles of the GC systems of our simulated S0s also typically extend out to |$\sim 10\, R_{\mathrm{e}}$|⁠, with 10 galaxies also reaching beyond |$10\, R_{\mathrm{e}}$|⁠. Additionally, the peaked V_rot/σ profiles of the GC systems of our S0 galaxies can have rotational velocities that reach their peak value at larger radii (i.e. beyond |$2\, R_{\mathrm{e}}$|⁠) than found by Schulze et al. (2020) for their peaked galaxies. Indeed, some of our S0 galaxies show V_rot/σ profiles peaking at |$\sim 5\, \mathrm{R_{\mathrm{e}}}$| (or beyond) and decreasing outwardly, so they would be classified as increasing according to the definition by Schulze et al. (2020) out to |$\sim 5\, R_{\mathrm{e}}$|⁠. Therefore, for these galaxies whose V_rot/σ profiles extend beyond |$\sim 5\, R_{\mathrm{e}}$|⁠, we also look at their kinematic behaviour at larger radii in order to properly classify them into one of the three categories from Schulze et al. (2020).

Fig. 5 (top) shows the 41 aligned galaxies (⁠|$82{{\ \rm per\ cent}}$| of total) split into the peaked, flat, and increasing V_rot/σ profile shapes from the top to the bottom panels, respectively. Each sub-figure of Fig. 5 (top) shows the 1D V_rot, σ, and V_rot/σ profiles of GC systems of the S0 galaxies in the low (10 < log (M_*/M_⊙) < 10.5; left column) and high (10.5 < log (M_*/M_⊙) < 11.6; right column) stellar mass bins.

We find that that the peaked, flat, and increasing V_rot/σ profile shapes occur in similar numbers, i.e. 10, 7, and 8, respectively, in the low stellar mass bin, while the peaked V_rot/σ profile shape is the dominant one in the high stellar mass bin. In fact, 10 galaxies show a peaked V_rot/σ profile shape as opposed to the flat and increasing V_rot/σ profile shapes that contain only three galaxies each at the high stellar masses. We also investigated the kinematics of the GC systems of the central and satellite galaxies. However, we did not see any clear differences, unlike previous works that found that S0 galaxies in clusters are typically more rotationally supported than S0 galaxies in the field ( Coccato et al. 2020; Deeley et al. 2020). We suggest that the lack of clear differences between the central and satellite galaxies in this work may be due to the fact that the E-MOSAICS simulations do not include the highest density environments of clusters, but they are limited to the field and group environments (i.e. Fornax cluster-like), as mentioned in Section 2.

The comparison with Schulze et al. (2020) simulations suggests that low-mass S0 galaxies (Fig. 5, left-hand side) can form through a range of different merger and accretion histories. Specifically, galaxies with peaked V_rot/σ profiles are expected to have typically formed through late (i.e. |$z$| < 1) mini mergers that mainly influence the kinematic properties of the outer regions of the galaxies without destroying their central disc-like kinematics. On the other hand, galaxies with flat V_rot/σ profiles are more likely to have experienced a late major merger that destroyed the central disc-like kinematics of the galaxies. This major merger may have also occurred in the form of multiple minor mergers, as a single major merger event is typically expected to spin-up the rotation of the GCs at large radii (Bekki et al. 2005). Finally, galaxies with increasing V_rot/σ profiles are suspected to have likely formed through a late gas-rich major merger that may have re-formed the central kinematically cold disc component of the galaxies, thus enhancing the rotation of the GCs out large radii (Bekki et al. 2005). Using the E-MOSAICS simulations, Trujillo-Gomez et al. (2021) also found that a late major merger is likely to produce more eccentric spatial distributions of the MR GCs compared to those of the MP GCs (i.e. increasing V_rot/σ profiles). On the other hand, galaxies characterized by a larger spread in GC eccentricities (low V_rot/σ profiles) are likely to have assembled |$50{{\ \rm per\ cent}}$| of their mass at earlier times (Trujillo-Gomez et al. 2021), consistent with the assembly history predicted for the peaked galaxies (i.e. two-phase formation scenario)

On the other hand, the high-mass S0 galaxies (Fig. 5, right-hand side) seem to be mostly dominated by a formation history from late mini mergers, as shown by the larger number of the peaked galaxies than the flat and increasing galaxies in the high stellar mass bin.

In Dolfi et al. (2021), we have studied the kinematic profiles of an observed sample of nine selected S0 galaxies from the SLUGGS survey extending out to |$\sim 5\, R_{\mathrm{e}}$| by combining the kinematic of the stars, GCs, and PNe. As mentioned above, eight out of the nine S0s from the SLUGGS survey cover similar stellar masses to the high-mass S0s from the simulations, i.e. 10.5 < log (M_*/M_⊙) < 11.6. From the observations, we find that six out of nine S0 galaxies (⁠|$67{{\ \rm per\ cent}}$|⁠) have GCs and PNe with consistent kinematics with respect to the underlying stars that rotate along the PA_phot of the galaxies (and, thus, are referred to as aligned galaxies)

Overall, the kinematic results of the GCs of these aligned galaxies from the SLUGGS survey seem to be consistent with those of the simulated S0 galaxy sample from E-MOSAICS. Specifically, in the range of stellar masses in common between the observations and simulations (i.e. 10.5 < log (M_*/M_⊙) < 11.6), we find that galaxies with peaked V_rot/σ profiles dominate with a |$\sim 62{{\ \rm per\ cent}}$| fraction in the simulations and |$\sim 60{{\ \rm per\ cent}}$| fraction in the observations. These are followed by the galaxies with flat V_rot/σ profiles that occur with a |$\sim 19{{\ \rm per\ cent}}$| fraction in the simulations and |$\sim 40{{\ \rm per\ cent}}$| fraction in the observations. We do not find evidence of any increasing V_rot/σ profile shapes in the observations, but this could simply be a result of sample selection effects or low number statistics. However, according to the simulations, galaxies with increasing V_rot/σ profile shapes should occur with similar frequencies as the galaxies with flat V_rot/σ profile shapes, i.e. |$\sim 19{{\ \rm per\ cent}}$|⁠.

Finally, across the entire stellar mass range covered by the simulated S0 galaxies (i.e. 10.0 < log (M_*/M_⊙) < 11.6), we find that the peaked V_rot/σ profile shape is still the dominant one among the aligned galaxies with |$\sim 49{{\ \rm per\ cent}}$| fraction as compared to the flat and increasing galaxies that are found to occur with a |$\sim 24{{\ \rm per\ cent}}$| and |$\sim 27{{\ \rm per\ cent}}$| fraction, respectively.

4.1.2 Misaligned galaxies

Fig. 5 (bottom) shows the 1D kinematic profiles of the GCs of the nine misaligned galaxies split into the low (10 < log (M_*/M_⊙) < 10.5; left column) and high (10.5 < log (M_*/M_⊙) < 11.6; right column) stellar mass bins, as for the aligned galaxies. Each sub-figure of Fig. 5 (bottom) also shows the 1D V_rot, σ and V_rot/σ profiles of GC systems as for the aligned galaxies, but this time we are also plotting the absolute difference between the PA_kin of the GCs and the PA_phot of the S0 galaxy to show the degree of misalignment of the GCs from the PA_phot of the galaxy.

As we see in Fig. 5 (bottom), the GC systems of the misaligned galaxies are rotating along an axis different from the PA_phot of the galaxy. Additionally, the V_rot/σ profile also show variations as a function of the galactocentric radius, with several galaxies being characterized by a second V_rot/σ peak at large radii. This suggests that the misaligned galaxies are characterized by a more complex assembly history that could likely be the result of multiple merger events occurring up to the most recent times in the formation history of the galaxies. Such an assembly history may have contributed to a significant fraction of the ex-situ GCs that have still not acquired the overall underlying rotation of the galaxy.

In Dolfi et al. (2021), we found that three out of the nine S0 galaxies (⁠|$33{{\ \rm per\ cent}}$|⁠) from the SLUGGS survey showed kinematic misalignment and twists in the rotation of the GCs and PNe with respect to the underlying stars as a function of radius. Specifically, we found that one galaxy (i.e. NGC 4649) was characterized by a similar ‘double-peaked’ V_rot/σ profile as we see for most of the misaligned simulated S0s. Therefore, these ‘double-peaked’ galaxies may be the results of a late major merger event, which spun up the rotation of the GCs at large radii (Bekki et al. 2005). On the other hand, the remaining two misaligned galaxies from the SLUGGS survey could be the result of multiple, late minor, and mini mergers that decreased the rotation of the GCs out to large radii.

In Section 6, we aim at investigating the presence of differences between the past merger histories of the aligned and misaligned galaxies to understand what physical processes are driving the GC misalignment in our S0 galaxies.

4.2 Kinematics of the MR and MP GC sub-populations

The colour-magnitude diagrams of the GC systems of our simulated S0 galaxies do not show evidence of any clear colour bimodality. For this reason, in Section 4.1, we did not adopt any colour cut to separate between the red and blue GC sub-populations. However, it is not clear whether this may be an issue of the simulations that are unable to reproduce the commonly observed colour-bimodality of the GC systems of observed galaxies.

In this section, we investigate whether the MR and MP GC sub-populations show any differences in their kinematic behaviours by adopting a fixed metallicity split at [Fe/H] = −1, which is the same for all the GC systems of our simulated S0 galaxies. This metallicity value is consistent with that adopted for the ETGs from the SLUGGS survey to separate between the two GC sub-populations (e.g. Usher et al. 2019). We note here that the issue of the overproduction of the MR GCs with −1.0 < [Fe/H] < 0 in the simulations could influence the comparison with the observations, specifically for the low-mass galaxies with 9.5 < log (M_*/M_⊙) < 10.5, which are found to be the most affected (see Section 3.2.4).

In order to ensure that each simulated S0 galaxy has a large enough sample of MR and MP GCs for reliably deriving the corresponding 1D kinematic profiles, we only consider here the 26 galaxies with a total number of MR and MP GCs greater than 40. This threshold is selected based on the comparison with the GC systems of the observed S0 galaxies from the SLUGGS survey, whose GC sub-populations contain more than 40 objects (see table 1 in Dolfi et al. 2021).

Following the classification in Section 4.1, we differentiate between aligned and misaligned galaxies. For the aligned galaxies, we show the kinematics of the MR and MP GCs by separating between the galaxies previously classified as peaked, flat, or increasing. Figs 6–9 show the 1D kinematic profiles of the MR and MP GCs for the 26 S0 galaxies, as selected above. Each sub-figure shows the PA_kin (where fitted), V_rot, σ, and V_rot/σ profiles from the top to the bottom panel for the MR (red line) and MP (blue line) GCs of each galaxy.

From the 1D kinematic profiles, we see that the MR GCs are rotating along the PA_phot of the galaxy within the 1σ errors in 9 out of the 11 peaked galaxies (see Fig. 6), as well as in all of the three increasing galaxies (see Fig. 7). Additionally, in the three increasing galaxies, the rotation of the MP GCs is, overall, consistent with the MR GCs, along the PA_phot of the galaxy. On the other hand, the MP GCs show little rotation compared to the MR GCs in seven of the peaked galaxies (see Figs 6a, c, e, g, h, i, and k) and, in two of these galaxies, their rotation is not consistent with the PA_phot of the galaxy (see Figs 6e and h). In one of the peaked galaxies, for which the PA_kin is not well constrained, the 2D kinematic maps do not show evidence of clear rotation for either the MR and MP GCs (see Fig. 6j).

These results may suggest different origins for the MR and MP GCs. In fact, the seven peaked galaxies with MP GCs characterized by low rotation suggest that this may be the result of an accretion history involving multiple low-mass (i.e. minor and mini) mergers from random directions, which would be consistent with the expected formation scenario of the peaked galaxies proposed by Schulze et al. (2020). On the other hand, the three peaked (see Fig. 6b, d, and f) and three increasing galaxies, with MR and MP GCs showing consistent rotation may have experienced fewer accretion events that preserved the overall rotation of the galaxy and of its GC sub-populations.

In Fig. 8, the flat galaxies have MR and MP GCs showing typically low rotation at all radii, i.e. V_rot/σ ≲ 0.6, as compared to the peaked galaxies. One galaxy (see Fig. 8c) also shows misaligned MR GC kinematics with respect to the PA_phot of the galaxy. Overall, this would seem to suggest a more violent and complex accretion history, possibly from a late major merger, for these galaxies that enhanced the random motion of the GC sub-populations out to large radii.

Finally, in Fig. 9, the misaligned galaxies show mixed behaviours in the kinematics of their MR and MP GCs. In fact, in four of the misaligned galaxies, the rotation of the MR GCs is, overall, aligned with the PA_phot of the galaxy, while the MP GCs are not (see Figs 9a, d, e, and f) similar to the two peaked galaxies described above. On the other hand, two misaligned galaxies are more similar to the flat galaxies, as they show misaligned kinematics in both the MR and MP GC sub-populations that are characterized by less well-behaved kinematic profiles with several transitions as a function of radius (see Fig. 9b and g). Of the remaining two misaligned galaxies, one (see Fig. 9h) is characterized by low rotation in both the MR and MP GCs with unconstrained PA_kin, while the other (see Fig. 9c) shows low rotation in the MP GCs.

In Fig. 6–9, we have only selected those galaxies with a large enough sample of MR and MP GCs to calculate the 1D kinematic profiles. Here, we investigate the kinematics of the MR and MP GCs of all 50 S0 galaxies based on the 2D line-of-sight velocity maps (edge-on projection) to include in the analysis also those galaxies with less than 40 GCs in one or both of the GC sub-populations. In Fig. 10, we summarize the fraction of the galaxies showing the following four different kinematic behaviours in their GC sub-populations as follows:

• MR = MP: the MR and MP GCs are rotating along a similar axis (PA_kin), which is consistent with the photometric major axis of the galaxy (PA_phot);

• MR ≠ MP: the MR GCs are rotating along an axis consistent with the photometric major axis of the galaxy (PA_phot), while the metal poor GCs are misaligned;

• MP low rotation: the MR GCs are rotating along an axis consistent with the photometric major axis of the galaxy (PA_phot), while the MP GCs have little or no rotation;

• Miscellaneous: the MR and MP GCs are both misaligned from the photometric major axis of the galaxy (PA_phot) or have little rotation.

For the comparison between the simulations and observations, in Fig. 10, we also show the distribution of the eight out of the nine S0 galaxies from the SLUGGS survey characterized by high mass (i.e. log (M_*/M_⊙) > 10.5), which we studied in Dolfi et al. (2021).

We find that the high mass (log (M_*/M_⊙) > 10.5) simulated S0 galaxies are distributed with similar fractions in each of the four kinematic classes. On the other hand, the low mass (log (M_*/M_⊙) < 10.5) simulated S0 galaxies are more frequently characterized by MP GCs with low rotation velocity (MP low rotation|$\sim 50{{\ \rm per\ cent}}$|⁠), suggesting that these galaxies have mainly experienced isotropic accretion events from low-mass dwarf galaxies that mostly contributed to few MP GCs with low rotation.

For the high-mass S0 galaxies, we find that the simulations predict a fraction of galaxies with misaligned MP and MR GC kinematics (MR ≠ MP|$\sim 30{{\ \rm per\ cent}}$|⁠) that is consistent with the observations. On the other hand, the simulations underestimate the fraction of galaxies with aligned MR and MP GC kinematics (MR = MP|$\sim 30{{\ \rm per\ cent}}$|⁠) with respect to the observations (MR = MP|$\sim 50{{\ \rm per\ cent}}$|⁠). This discrepancy may be a result of sample selection effects as some of our simulated galaxies show signs of interaction from the Mpc-scale FoV images, while the observed galaxies from the SLUGGS survey are mainly undisturbed. Alternatively, it could be a result of viewing angle effects, as we are studying the edge-on view of the galaxies (see Section 5). Finally, we find an overall good agreement between the fraction of the simulated and observed S0 galaxies in the remaining two kinematic classes (MP low rotation and Miscellaneous) within the errors. Here, the discrepancy could also be the result of sample selection and viewing angle effects, as previously mentioned. Alternatively, it could also be due to differences in the splitting between the GC sub-populations that is colour-based in the observations and metallicity-based in the simulations.

5 PHASE-SPACE DIAGRAMS OF THE GC SYSTEMS

The aim of this section is to understand whether the location of the GC systems of the individual galaxies on the 3D phase-space diagrams correlate with their accretion times on to the galaxies. The results could subsequently be applied to the GC systems of observed galaxies for which we do not have explicit information about their in situ or ex situ origins as in the simulations.

In their study of simulated clusters, Rhee et al. (2017) separated the galaxies into three different accretion classes based on their accretion time on to the cluster: ancient infallers with t_inf(Gyr) > 6.45, intermediate infallers with 3.63 < t_inf(Gyr) < 6.45, recent infallers with t_inf(Gyr) < 3.63 and first infallers that have just fallen on to the cluster. They found that the ancient and recent infallers were mainly located within one virial radius from the centre of the cluster (i.e. |$\sim 1\, R_{\mathrm{200}}$|⁠) and that the recent infallers had higher orbital velocities than the ancient infallers. On the other hand, the intermediate and first infallers were mainly located between 1 and |$3\, R_{\mathrm{200}}$|⁠, with the intermediate infallers being characterized by lower orbital velocities than the first infallers.

Focusing on only the ex situ GC population of the sample of 50 S0 galaxies, which are defined by the E-MOSAICS simulations, we do not apply any surface density cut here, unlike in Section 3.2.5. Therefore, the spatial distributions of the GC systems will be more extended. Following Rhee et al. (2017), we separate the ex situ GCs based on the time they were accreted on to the host galaxy as defined by the simulations: ancient infallers (t_inf(Gyr) > 6.45), intermediate infallers (3.63 < t_inf(Gyr) < 6.45), and recent infallers (t_inf(Gyr) < 3.63), and we produce the 3D phase-space diagrams of each one of these accreted GC classes. The 3D phase-space diagrams are produced by plotting the 3D velocity (V_3D) as a function of the 3D galactocentric radius (R_3D) of each GC. We normalize the V_3D and R_3D of each GC by the velocity dispersion of the GC system, σ_3D, and virial radius of the galaxy, R₂₀₀, respectively, to which the GC belongs. For the central galaxies, R₂₀₀ is measured at |$z$| = 0. However, for the satellite galaxies, R₂₀₀ is not defined and can only be measured for the group as a whole. For this reason, for the satellite galaxies, we use the R₂₀₀ measured when the galaxy was last a central of its own group prior to becoming a satellite. The velocity dispersion, σ_3D, is calculated from the line-of-sight velocity dispersion along the x-, y-, and |$z$|-axes estimated as the standard deviation of the corresponding GC line-of-sight velocities.

Fig. 11 shows the 3D phase-space diagrams of the in situ GCs and ex situ GCs, split into the ancient, intermediate, and recent infaller GCs, of all our 50 S0 galaxies from the left to the right panel, respectively. In the left-most panel, we show together the in situ GCs and the ancient infaller GCs normalized by the total number of GCs and the total number of ex situ GCs of all our 50 S0 galaxies, respectively. The right-most panel shows together the ex situ GCs normalized by the total number of GCs of all our 50 S0 galaxies. The intermediate and recent infaller GCs are normalized the total number of ex situ GCs of all our 50 S0 galaxies, as for the ancient infaller GCs.

First of all, we see that the in situ and ex situ GCs, which are split into ancient, intermediate, and recent infallers GCs, show similar spatial distributions extending out to ∼1–|$2\, R_{\mathrm{200}}$|⁠. However, the majority of the in situ and ex situ GCs are contained within |$1\, R_{\mathrm{200}}$|⁠, with only a total of 10 GC systems among all our 50 S0 galaxies extending beyond |$1\, R_{\mathrm{200}}$|⁠.

Secondly, we see a trend that the GCs that were accreted more recently (t_inf(Gyr) < 6.45) are more likely found at larger radii (see also Kruijssen et al. 2020; Pfeffer et al. 2020). In fact, |$\sim 50{{\ \rm per\ cent}}$| of the ex situ GCs located within |$0.2\, R_{\mathrm{200}}$| are ancient infallers (t_inf(Gyr) > 6.45). The in situ GCs are also concentrated within |$0.2\, R_{\mathrm{200}}$|⁠, where they overlap with the ancient infaller GCs. On the other hand, only |$\sim 30{{\ \rm per\ cent}}$| of the ex situ GCs within |$0.2\, R_{\mathrm{200}}$| are intermediate infallers (3.63 < t_inf(Gyr) < 6.45) and |$\lt 30{{\ \rm per\ cent}}$| are recent infallers (t_inf(Gyr) < 3.63). The intermediate infallers are more numerous between 0.5 and |$1\, R_{\mathrm{200}}$| for V_3D/σ_3D < 1 with |$\ge 50{{\ \rm per\ cent}}$| fraction, while |$\ge 50{{\ \rm per\ cent}}$| of the recent infallers are found beyond |$1\, R_{\mathrm{200}}$| for V_3D/σ_3D < 1 or within |$1\, R_{\mathrm{200}}$| for V_3D/σ_3D > 1.

Overall, we see that the fraction of the in situ and ancient infaller GCs sharply drops beyond |$0.2\, R_{\mathrm{200}}$| as opposed to the ex situ GCs, whose fraction increases to |$\ge 70{{\ \rm per\ cent}}$| out to |$\sim 1\, R_{\mathrm{200}}$|⁠. Additionally, we note that the in situ GCs are also characterized by larger velocities on average than the ex situ GCs within |$0.2\, \mathrm{\mathit{ R}_{200}}$|⁠. This may be due to the nature of the in situ GCs, which are mainly MR in a fast-rotating disc configuration.

The results from Fig. 11 are overall consistent with those of the galaxies in individual galaxy clusters found by Rhee et al. (2017), suggesting that similar physical processes do apply to both GCs and galaxies accreting on to an individual galaxy or cluster of galaxies, respectively. However, we note that the boundaries between the different accreted classes of the GCs in the galaxies are less sharp than those of the accreted galaxies in galaxy clusters. In fact, as we previously mentioned, the ancient, intermediate, and recent GCs have overall similar distributions on the 3D phase-space diagrams. This may likely be a result of the degeneracy between the accretion redshift and the stellar mass of the accreted satellite galaxy. In fact, as shown in fig. 2 of Pfeffer et al. (2020), the apocentres of the accreted GCs decrease for increasing satellite stellar masses at fixed redshifts.

In summary, our results suggest that there exists a correlation between the location of the GCs of the individual galaxies on the 3D phase-space diagrams and the time they were accreted on to the galaxy. We also looked at the projected (2D) phase-space diagrams, along the three line-of-sights, i.e. x, y, and |$z$|⁠. We find that, overall, the boundaries delimiting the areas where the different accretion GC classes are mainly located are similar in the 3D and 2D phase-space diagrams. Therefore, the phase-space diagrams could be applied to the GC systems of the observed galaxies and we should expect that the GCs, found within |$0.2\, R_{\mathrm{200}}$| (⁠|$=20\, \mathrm{kpc}$| for |$R_{\mathrm{200}}=100\, \mathrm{kpc}$|⁠), have likely an in situ origin or they were accreted at early times (more than |$6.45\, \mathrm{Gyr}$| ago). On the other hand, the GCs found beyond |$0.2\, R_{\mathrm{200}}$| have a higher probability of being ex situ and accreted more recently (less than |$6.45\, \mathrm{Gyr}$| ago) on to the galaxy.

6 THE ASSEMBLY HISTORY OF THE SIMULATED S0 GALAXIES

One of the advantages of the simulations is that they contain information about the detailed merger and accretion histories of the galaxies. For this reason, we now investigate the presence of correlations between the distinct kinematic V_rot/σ profiles of the GCs, described in Section 4.1, and the corresponding assembly histories of the S0 galaxies. This adds to the previous work using E-MOSAICS, which has revealed correlations between the galaxy assembly history and GC age-metallicity distributions (Kruijssen et al. 2019) and GC kinematics (Trujillo-Gomez et al. 2021).

Furthermore, in Section 4.2, we have seen that the MR and MP GC sub-populations can have misaligned kinematics. Therefore, we aim to understand whether the observed kinematic misalignment in the GC sub-populations of our simulated S0 galaxies is the result of specific merger events.

6.1 The relation between the V_rot/σ profiles of the GCs and the accretion history of the S0 galaxies

From the E-MOSAICS simulations, we derive the merger trees of our simulated 50 S0 galaxies. These merger trees include all the accretion events with different mass ratio (i.e. |$\mathrm{M_{*}^{\mathrm{sat}}}/\mathrm{M_{*}}$|⁠, where |$\mathrm{M_{*}^{\mathrm{sat}}}$| is the satellite galaxy mass and M_* is the target galaxy stellar mass) experienced by our target S0 galaxy as a function of redshift.

Following Schulze et al. (2020), we classify the merger events into four categories depending on the mass ratio of the merger: major mergers (mass-ratio>1:4), minor mergers (1:10<mass-ratio<1:4), mini mergers (1:100<mass-ratio<1:10), and smooth accretion (mass-ratio<1:100). We focus here on the mergers that occurred within the last |$10\, \mathrm{Gyr}$| for the comparison with Schulze et al. (2020).

On the left-hand side of Fig. 12, we show the fraction of the peaked, flat, and increasing galaxies that experienced their last major, minor, and mini merger event within a time interval of |$\pm 1\, \mathrm{Gyr}$| for a given look back time in Fig. 12 up to a maximum of |$10\, \mathrm{Gyr}$|⁠, similarly to fig. 10 of Schulze et al. (2020). Each panel is normalized by the total number of the galaxies of a given V_rot/σ class that experienced their last major, minor, and mini merger within the last |$10\, \mathrm{Gyr}$|⁠, such that the sum data points of a given type is equal to one in each panel. However, overall, we do not see any clear trend between the merger histories of the different galaxy types as a function of redshift, as it was previously found by Schulze et al. (2020). In fact, Schulze et al. (2020) found that |$\sim 60{{\ \rm per\ cent}}$| of the peaked galaxies had no major mergers and that, if a major merger occurred, it was most probably more than |$5\, \mathrm{Gyr}$| ago. On the other hand, they found that only |$40{{\ \rm per\ cent}}$| and |$30{{\ \rm per\ cent}}$| of the flat and increasing galaxies experienced no major mergers, respectively, while the remaining ones more likely experienced at least a major merger between 3 and |$7\, \mathrm{Gyr}$| ago. Finally, Schulze et al. (2020) found that all galaxies experienced typically between 1 and 4 mini mergers with a probability increasing towards low redshifts.

On the right-hand side of Fig. 12, we show the fraction of the peaked, flat, and increasing galaxies that experienced a total number of major, minor, and mini mergers in the last |$10\, \mathrm{Gyr}$| from the top to the bottom panel, respectively. We see that the majority of the galaxies are likely to have experienced either zero or one major merger. A small fraction of the galaxies (∼10-|$20{{\ \rm per\ cent}}$|⁠) may have also experienced two major mergers, however the peaked and increasing galaxies are unlikely to have experienced more than two. On the other hand, a small fraction of the flat galaxies may have experienced up to four major mergers. Almost all of the increasing galaxies (⁠|$80{{\ \rm per\ cent}}$|⁠) did not experience any minor merger, while roughly half of the peaked and flat galaxies may have experienced up to two or three minor mergers, respectively. Approximately |$40{{\ \rm per\ cent}}$| of the peaked galaxies have had at least one mini merger and a similar fraction of the flat galaxies have had at least three. On the other hand, half of the increasing galaxies have had no mini mergers, while the remaining half is likely to have experienced between one and two mini mergers. These results show that the peaked, flat, and increasing galaxies are overall similar in terms of the number of the different types of mergers that they experienced in the last |$10\, \mathrm{Gyr}$|⁠. An exception is represented by the frequency of minor mergers that are more typical for the peaked and flat than the increasing ones.

The discrepancy with the results found by Schulze et al. (2020) may be due to several possible reasons. It could be a result of low sample statistics as we only have 50 galaxies compared to the ∼500 of Schulze et al. (2020). Additionally, the fact that we are also only selecting S0 galaxies in this work could add a certain bias. Secondly, it could be due to the fact that we are studying the V_rot/σ profiles of the GCs, while Schulze et al. (2020) looked at the stellar V_rot/σ profiles. Additionally, our V_rot/σ profiles extend beyond |$5\, R_{\mathrm{e}}$|⁠, while those from Schulze et al. (2020) were limited to |$5\, R_{\mathrm{e}}$|⁠. Finally, it could be a result of differences between the EAGLE and Magneticum simulations.

To test whether our results may partly depend on the radial extension of the V_rot/σ profiles, we perform an alternative classification of the V_rot/σ profiles restricted to |$5\, R_{\mathrm{e}}$|⁠. However, the results in Fig. 12 do not change significantly and we still do not see the trends found by Schulze et al. (2020) for the major and mini mergers specifically.

In a previous work, Pulsoni et al. (2021) studied the relation between the stellar V_rot/σ profiles out to |$15 R_{e}$|and the past accretion histories of the ETGs in the IllustrisTNG simulations. Interestingly, they found evidence that the stellar V_rot/σ profiles did not strongly depend on the mass-ratio of the mergers, but mainly on the presence or absence of gas. Therefore, the results from Pulsoni et al. (2021) specifically suggest that the peakedstellar V_rot/σ profiles are shaped by the in-situstellar component and recent cold gas accretion that helped preserving the central peak in rotation of the stars (and not by the late mini mergers as suggested by Schulze et al. 2020). Additionally, Pulsoni et al. (2021) have also found that the location of the kinematic transition between the fast-rotating and the slowly rotating spheroidal component in the peakedgalaxies does not necessarily correspond to the transition betweem the in-situand ex-situdominated regions of the galaxies, as previously suggested by Schulze et al. (2020).

The results in Fig. 12 could possibly suggest that the distinct present-day kinematic profiles of the GCs are not the result of the specific different merger histories of the S0 galaxies, in agreement with the results from Pulsoni et al. (2021). Other factors, such as the orbital configuration of the merging galaxies or the presence/absence of gas, may be playing a more dominant role. A possible follow-up project could involve investigating the kinematic properties of the GC systems back in time to understand whether and how the different physical processes (e.g. mergers) influence the GC kinematics in our galaxies during their evolution and increase the sample size beyond 50 galaxies.

6.2 The origin of the GC misalignment in S0 galaxies

In Section 4.1, we classified 9 out of the 50 simulated S0 galaxies as misaligned, since their GCs were rotating along a different PA_kin compared to the PA_phot of the galaxy.

In this section, we aim to investigate what is the cause of the misalignment in these galaxies and, specifically, whether it is a result of the specific merger histories of the misaligned galaxies with respect to the other aligned galaxies. For each one of the nine misaligned galaxies, we identify those accretion events that each contributed to a significant fraction of the present-day GC population (i.e. |$\ge 10{{\ \rm per\ cent}}$|⁠), and we study the kinematic properties of these accreted GCs, with the main focus being their 2D velocity map and 1D PA_kin and V_rot profiles.

We find that five of the nine misaligned galaxies (⁠|$56{{\ \rm per\ cent}}$|⁠) have experienced a major merger event within the last |$10\, \mathrm{Gyr}$| that contributed between 30 and |$80{{\ \rm per\ cent}}$| of the total present-day GC population. Three of these galaxies have also experienced a minor merger event between 4 and |$7\, \mathrm{Gyr}$| ago that contributed between |$10{{\ \rm per\ cent}}$| and |$20{{\ \rm per\ cent}}$| of the GC population. The ex situ GCs are largely MR (i.e. |$\gtrsim 40{{\ \rm per\ cent}}$|⁠) in these five galaxies and they show different rotation properties with respect to the in situ GCs. In fact, while the in situ GCs show typically strong rotation, which is overall well aligned with the PA_phot of the galaxy, the ex situ GCs show rotation which can be misaligned by angles varying between 0° < ΔPA < 180° (where 0° is co-rotating and 180° is counter-rotating) with respect to that of the in situ GCs. In Fig. 13, we show the edge-on projection of one of these five galaxies (i.e. GalaxyID =4362302) as an example, which experienced a major merger |$\sim 9\, \mathrm{Gyr}$| ago that contributed to |$\sim 70{{\ \rm per\ cent}}$| of the total present-day GC population with |$60{{\ \rm per\ cent}}$| of these ex situ GCs being MR. We see that the ex situ GCs are counter-rotating with respect to the in situ GCs in this galaxy.

We also find that few of the aligned galaxies in Section 4.1 (e.g. GalaxyID =13394165) have accreted a significant fraction of the ex situ MR GCs from a major merger (i.e. |$\sim 40{{\ \rm per\ cent}}$|⁠), but these ex situ MR GCs are co-rotating with the in situ MR GCs. These aligned galaxies have experienced the major merger event between 6 and |$9\, \mathrm{Gyr}$| ago, similarly to the five misaligned galaxies. Therefore, we suggest that the misalignment between the in situ and ex situ MR GCs in these galaxies is not necessarily connected to the time of occurrence of the major merger, but it, possibly, depends on other properties of the mergers, such as the orbital configuration of the merging galaxies. In this scenario, the misalignment between the in situ and ex situ GCs could then be the result of the change of the axis of accretion from the cosmic web during the galaxy formation process, such that the spin axis of the accreted GCs was aligned differently with respect to the spin axis of the host galaxy.

The remaining four of the nine misaligned galaxies (⁠|$44{{\ \rm per\ cent}}$|⁠) did not experience any significant accretion events and, thus, they have low ex situ GC fractions (i.e. |$\lesssim 30{{\ \rm per\ cent}}$| of the total present-day GC population). The ex situ GCs, which are mostly MP in these galaxies, do not show evidence of any clear rotation. For this reason, we investigate the 2D velocity maps and the 1D PA_kin and V_rot profiles of the in situ MR and in situ MP GCs in these galaxies. In two of these galaxies, we find that the rotation of the in situ MR GCs is, overall, aligned with the PA_phot of the galaxy, while the in situ MP GCs are rotating along a kinematic axis offset by ΔPA ≃ 90° with respect to the PA_phot of the galaxy. In other two galaxies, the in situ MP GCs do not show evidence of any net rotation, while the in situ MR GCs show rotation that oscillates around the PA_phot of the galaxy, but it is overall consistent within the large 1σ errors. In Fig. 13, we show one of these four galaxies (i.e. GalaxyID =6451189) as an example, whose in situ MP GCs are misaligned by ∼90° with respect to the in situ MR GCs.

From the 2D kinematic maps of all 50 S0 galaxies, we find that the majority of the galaxies have the in situ GCs co-rotating with the stars about the photometric major-axis of the galaxy (i.e. |$70{{\ \rm per\ cent}}$|⁠). On the other hand, a small fraction of the galaxies have in situ GCs that are not co-rotating with the stars (i.e. |$16{{\ \rm per\ cent}}$|⁠) or that do not show evidence of any net rotation (i.e. |$14{{\ \rm per\ cent}}$|⁠). Additionally, we find that |$20{{\ \rm per\ cent}}$| of the S0 galaxies have the in situ MP GCs that are rotating about an axis different from that of the in situ MR GCs, which is overall consistent with the photometric major axis of the galaxy (e.g. GalaxyID =13394165 and 17827697 in Fig. 6).

We do not find any clear differences in the merger histories of the galaxies with misaligned in situ MR and in situ MP GC kinematics. Therefore, these results suggest that the misalignment of the in situ GC sub-populations is not related to the specific merger histories of the galaxies.

The in situ MP GCs likely formed at early times, i.e. |$z$| > 2, in the turbulent and high-pressure discs of low-mass gas-rich galaxies. Prior to their formation, the in situ MP GCs suffered kinematic perturbations from merger events (more frequent at high redshifts) that redistributed these GCs in the host galaxy halo where they survived until the present-day (Kruijssen 2015; Keller et al. 2020). The in situ MR GCs likely formed more recently after the host galaxy built up more of its stellar mass. For this reason, the in situ MR GCs are less likely to have suffered kinematic perturbations from merger events (less frequent and with smaller masses) and are generally co-rotating with the host galaxy stars.

We find that the in situ MR GCs are typically characterized by a more centrally concentrated and flattened spatial distribution than the in situ MP GCs, and they are rotating consistent with the stars about an axis consistent with the photometric major axis of the galaxy in a disc-like configuration (i.e. |$70{{\ \rm per\ cent}}$|⁠), as commonly found in the observations (e.g. Forbes 1997; Brodie & Strader 2006; Arnold et al. 2011; Pota et al. 2013; Dolfi et al. 2020, 2021). On the other hand, the in situ MP GCs typically have a more extended spatial distribution.

This in-situ GC formation scenario could also explain the misalignment of the in situ MR GCs seen in |$16{{\ \rm per\ cent}}$| of the galaxies if the merger events had a large enough mass ratio to perturb the orbits of these GCs and the subsequent star formation produced a stellar disc that was aligned differently from these previously formed in situ MR GCs.

Additionally, previous observational works have also found evidence of ETGs showing ‘kinematically distinct cores’ (KDC). These types of galaxies are generally not very common (i.e. |$\sim 7{{\ \rm per\ cent}}$| identified in ATLAS^3D by Krajnović et al. 2011), but some previous works have shown that the stars and GC sub-populations in these galaxies may have distinct and misaligned rotation (e.g. Blom et al. 2012), suggesting different origins. Therefore, the distinct rotational properties of the in situ MR and in situ MP GCs could also have their origins in a similar decoupling event that leads to KDCs.

7 SUMMARY AND CONCLUSIONS

In this work, we have studied the kinematic profiles of the GCs in a selected sample of 50 S0 galaxies from the E-MOSAICS simulations (Pfeffer et al. 2018; Kruijssen et al. 2019), with the aim of finding the link between the present-day kinematic properties and the past merger histories of the S0 galaxies in low-density environments (i.e. field and small galaxy groups).

From the 2D kinematic maps and 1D kinematic profiles of the GCs (see Section 4.1), we find that |$82{{\ \rm per\ cent}}$| of the galaxies have GCs that are rotating about an axis that is consistent with the photometric major axis of the galaxy (aligned), while the remaining |$18{{\ \rm per\ cent}}$| of the galaxies do not (misaligned), in general agreement with the observations from the SLUGGS survey (Dolfi et al. 2021)

Among the aligned galaxies, we find that |$49{{\ \rm per\ cent}}$|⁠, |$24{{\ \rm per\ cent}}$|⁠, and |$27{{\ \rm per\ cent}}$| show a peaked, flat, and increasing V_rot/σ profile, respectively (see Section 4.1). From the comparison of these three distinct V_rot/σ profiles with the past merger histories of the S0 galaxies derived from the E-MOSAICS simulations (see Section 6.1), we do not find evidence of any clear correlation between the present-day V_rot/σ profile shape and the past merger histories of the galaxies, unlike in the Magneticum simulations of the stellar profiles by Schulze et al. (2020). A future work could investigate the kinematic properties of the GC systems back in time to understand whether and how the different physical processes can influence the GC kinematics in our galaxies during their evolution and increase the sample size beyond 50 galaxies.

Among the misaligned galaxies, we find that (see Section 6.2):

• |$56{{\ \rm per\ cent}}$| of the galaxies experienced at least one major merger event that contributed between |$30{{\ \rm per\ cent}}$| and |$80{{\ \rm per\ cent}}$| of the total present-day GC population. In these galaxies, a large fraction of the ex situ GCs are MR and they are rotating about an axis that is different from that of the in situ MR GCs (which is consistent with the photometric major axis of the galaxy). The misalignment angle between the ex situ and in situ MR GCs can vary between 0° (co-rotation) and 180° (counter-rotation), while the MP GCs do not show evidence of any clear rotation in these galaxies (see Fig. 13; top panel).

• |$44{{\ \rm per\ cent}}$| of the galaxies did not experience any major merger event and have low ex situ GC fraction (i.e. |$\lesssim 30{{\ \rm per\ cent}}$|⁠). In these galaxies, the ex situ GCs are predominantly MP with no evidence of any clear rotation. In two of these galaxies, the in situ MP GCs are rotating about an axis offset by 90° from that of the in situ MR GCs (which is consistent with the photometric major axis of the galaxy).

For the total sample, from the 2D kinematic maps, we find that the majority of the galaxies (i.e. |$70{{\ \rm per\ cent}}$|⁠) have the in situ GCs that are co-rotating with the stars about an axis that is consistent with the photometric major axis of the galaxy, as commonly found in the observations (e.g. Forbes 1997; Brodie & Strader 2006; Arnold et al. 2011; Pota et al. 2013; Dolfi et al. 2020, 2021). A minority of the galaxies (i.e. |$16{{\ \rm per\ cent}}$|⁠) have the in situ GCs that are rotating about an axis that is different from the photometric major axis of the galaxy, while the remaining |$14{{\ \rm per\ cent}}$| of the galaxies have the in situ GCs that do not show evidence of any net rotation. Finally, for the total sample, we find that |$20{{\ \rm per\ cent}}$| of the galaxies have the in situ MP GCs that are rotating about an axis different from that of the in situ MR GCs. Therefore, we suggest that there are likely two possible origins for the GC misalignment in our S0 galaxies:

• in situ GC formation during the evolution of the host galaxy. The in situ MP GCs formed early (⁠|$z$| > 2) and suffered kinematic perturbations from the frequent merger events, which redistributed these GCs in the host galaxy halo (Kruijssen 2015; Keller et al. 2020). The in situ MR GCs are less likely to have suffered kinematic perturbations since they formed more recently when mergers became less frequent and the galaxy built up more of its stellar mass.

• a major merger event that contributed to a large fraction of the ex situ MR GCs that can be misaligned by angles varying between 0° (co-rotation) and 180° (counter-rotation) with respect to the in situ MR GCs (which are rotating consistent with the photometric major axis of the galaxy). The degree of misalignment will depend on the orbital configuration of the merging galaxies.

Finally, we have studied the distribution of the GC systems of our S0 galaxies in 3D phase-space diagrams (see Section 5). We find that there is a correlation between the accretion time of the GCs on to the galaxy with its location on the 3D and 2D phase-space diagrams, with the in situ and ancient infaller GCs (accreted more than |$6.45\, \mathrm{Gyr}$| ago) being predominantly located within |$0.2\, \mathrm{\mathit{ R}}_{\mathrm{200}}$|⁠. On the other hand, the GCs that were accreted more recently (less than |$6.45\, \mathrm{Gyr}$| ago) are most likely found at larger virial radii, i.e. beyond |$0.2\, \mathrm{\mathit{ R}}_{200}$| (see Fig. 11). Therefore, in future work, we can use the 2D phase-space diagrams of observed GC systems to assign a probability of the GCs of being in situ or ancient/recent infallers.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their very constructive comments and suggestions that helped improving this paper. We thank the E-MOSAICS team for providing access to the data in the simulations used in this paper. DF, WC, KB, and AD acknowledge support from the Australian Research Council under Discovery Project 170102344. JP acknowledges support from the Australian Research Council under Discovery Project 200102574. AJR was supported by National Science Foundation grant AST-1616710 and as a Research Corporation for Science Advancement Cottrell Scholar. JMDK gratefully acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) in the form of an Emmy Noether Research Group (grant number KR4801/1-1), as well as from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme via the ERC starting Grant MUSTANG (grant agreement number 714907).

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is a part of the National E-Infrastructure. The work also made use of high-performance computing facilities at Liverpool John Moores University, partly funded by the Royal Society and LJMU’s Faculty of Engineering and Technology.

DATA AVAILABILITY

The data are available upon request from the E-MOSAICS project (https://www.astro.ljmu.ac.uk/~astjpfef/e-mosaics/).

Footnotes

REFERENCES