The distribution of stellar orbits in eagle galaxies – the effect of mergers, gas accretion, and secular evolution

ABSTRACT

The merger history of a galaxy is thought to be one of the major factors determining its internal dynamics, with galaxies having undergone different types or mergers (e.g. dry, minor, or major mergers) predicted to show different dynamical properties. We study the instantaneous orbital distribution of galaxies in the eagle simulation, colouring the orbits of the stellar particles by their stellar age, in order to understand whether stars form in particular orbits (e.g. in a thin or thick disc). We first show that eagle reproduces well the observed stellar mass fractions in different stellar orbital families as a function of stellar mass and spin parameter at z = 0. We find that the youngest stars reside in a thin disc component that can extend to the very inner regions of galaxies, and that older stars have warmer orbits, with the oldest ones showing orbits consistent with both hot and counter-rotating classifications, which is consistent with the trend found in the Milky Way and other disc galaxies. We also show that counter-rotating orbits trace galaxy mergers – in particular dry mergers, and that in the absence of mergers, counter-rotating orbits can also be born from highly misaligned gas accretion that leads to star formation.

1 INTRODUCTION

A galaxy’s assembly history is thought to be one of the major factors determining its internal kinematic structure (e.g. White 1979; Fall & Efstathiou 1980; Park et al. 2019). It is thus expected that the internal orbital structure of a galaxy should contain important clues about a galaxy’s past. In particular, stellar dynamics provide a fossil record of the formation history of galaxies. Stars that were born and remain in quiescent environments tend to be on regular rotation-dominated orbits, so their galaxies are likely dominated by fast-rotating, flat stellar discs with recent star formation (e.g. Fall & Efstathiou 1980; Lagos et al. 2017). On the other hand, stars born from turbulent gas or that have been dynamically heated after birth will be on warmer orbits with more random motions (e.g. Leaman et al. 2017). Dynamical heating mechanisms that can affect stars include major mergers (e.g. Benson et al. 2004; House et al. 2011; Few et al. 2012; Helmi et al. 2012; Ruiz-Lara et al. 2016; Lagos et al. 2018a) and long-term secular heating of the disc via internal instabilities (e.g. Jenkins & Binney 1990; Aumer, Binney & Schönrich 2016; Grand et al. 2016), for example. Major gas-poor mergers, in particular, can play a key role in destroying stellar discs and in creating dynamically hot spheroidal components in the remnant galaxies (Barnes & Hernquist 1996; Cox et al. 2006; Bois et al. 2010, 2011; Hoffman et al. 2010; Naab et al. 2014; Pillepich, Madau & Mayer 2015; Lagos et al. 2018a).

Stellar kinematics are also expected to be systematically correlated with stellar ages (e.g. van de Sande et al. 2018; Trayford et al. 2019). Observations have found that old stars dominate the light of pressure-supported bulges, while younger stars live on thinner discs. In general, stars born together tend to live on similar orbits (e.g. Bird et al. 2013; Stinson et al. 2013). If we study the phase-space distribution of stars or stellar components in a galaxy, we can potentially infer information on whether these components have formed in the galaxy itself or are the results of merging events. It is expected that both merging events and star formation episodes play an important role in setting a galaxy’s structure and its stellar populations. Thus, the more information we can gather about the different orbital structures in a galaxy, the closer we can get to unveiling its evolutionary history.

Dynamical models can be used to infer the orbital distributions within galaxies. One way to parametrize this distribution is through the ‘circularity plane’; the distribution of orbits in the λ_z−R plane, where λ_z represents the z-axis component of the angular momentum of the stars normalized by the maximum angular momentum a star with the same binding energy could have (i.e. that of a circular orbit), and the time-averaged radius, r, as a proxy for the binding energy. The distribution of stellar orbits in a galaxy can thus be represented by the probability density of orbits in the phase–space of λ_z versus r. The different dynamical structures of galaxies can also be identified within this phase–space plane with, for example, the subdivision of the stellar component in a galaxy into cold (λ_z > 0.8), warm (0.25 < λ_z < 0.8), hot (|λ_z| < 0.25), and counter-rotating (λ_z < −0.25) orbits (e.g. Zhu et al. 2018a, b; Jin et al. 2019, 2020; Poci et al. 2021, 2022; Santucci et al. 2022).

In recent years, Integral Field Spectroscopy (IFS) surveys such as SAURON (Spectroscopic Areal Unit for Research on Optical Nebulae; de Zeeuw et al. 2002), ATLAS^3D (Cappellari et al. 2011), CALIFA (Calar Alto Legacy Integral Field Array survey; Sánchez et al. 2012), SAMI (Sydney-Australian-Astronomical-Observatory Multi-object Integral-Field Spectrograph) Galaxy Survey (Croom et al. 2012; Bryant et al. 2015; Croom et al. 2021), MASSIVE (Ma et al. 2014), MaNGA (Mapping Nearby Galaxies at Apache Point Observatory; Bundy et al. 2015), and the Fornax 3D survey (Sarzi et al. 2018) have provided us with rich resolved views of galaxies, allowing their structure and evolution to be investigated in detail through the mapping of stellar kinematics across individual galaxies. These IFS surveys have made possible the use of techniques such as the Schwarzschild orbit-superposition method (Schwarzschild 1979) to estimate the internal mass distribution, intrinsic stellar shapes, and orbit distributions of galaxies across the Hubble sequence. The stellar orbit circularity distribution within 1 effective radius (R_e, radius of the isophote containing half of the total luminosity; de Vaucouleurs 1948) has been obtained for 300 CALIFA galaxies across the Hubble sequence (Zhu et al. 2018a, b) and for a few hundred early-type galaxies (ETGs) in MaNGA (Jin et al. 2020), SAMI (Santucci et al. 2022), and ATLAS 3D (Thater et al. 2023). These works show that the fraction of dynamically hot orbits increases from low-mass to high-mass galaxies in the local Universe, consistently with an increase in the mass fraction of a classical bulge (i.e. those with triaxial shape, high Sérsic index, and dominated by velocity dispersion; Weinzirl et al. 2009). However, even though the studies of galaxies’ orbital components in large IFS surveys have already started, the formation of these orbital components has not yet been quantitatively connected to the merger history of galaxies.

The details of merger histories have been (partially) uncovered in the nearest galaxies, where either single stars in the stellar haloes can be resolved, for example in the case of the Milky Way (Belokurov et al. 2018; Helmi et al. 2018; Helmi 2020), the Andromeda galaxy (D’Souza & Bell 2018b, a), and NGC 5128 (Rejkuba et al. 2011), or where individual globular clusters in the haloes are detected (Forbes et al. 2016; Beasley et al. 2018; Kruijssen et al. 2019), or where streams have been detected (Merritt et al. 2016; Harmsen et al. 2017). Recently, studies have attempted to constrain the global ex situ fractions – the amount of stellar mass that formed in galaxies other than the main galaxy studied (Boecker et al. 2020; Spavone et al. 2020; Davison et al. 2021; Angeloudi et al. 2023), or the typical mass of the accreted satellites (Pinna et al. 2019; Martig et al. 2021) combining stellar populations and orbital distributions in edge-on lenticular galaxies obtained from IFS data (Poci et al. 2019, 2021), as well as spiral galaxies in the CALIFA survey (Jin et al. 2023). However, so far the assembly and merger histories of most galaxies remain hidden, including in the nearby Universe.

Cosmological hydrodynamical simulations make it possible to trace back the formation of galactic structures in great detail (Lagos et al. 2018b; Obreja et al. 2019; Trayford et al. 2019; Du et al. 2021; Pulsoni et al. 2021; Lagos et al. 2022; Proctor et al. 2024), including adopting dynamical decomposition methods similar to those applied to observational data (Xu et al. 2019), finding quantitatively consistent results between simulations and observations. The number of studies connecting the dynamical structures of galaxies with their accretion histories is increasing. It is then essential to exploit the additional information we can derive from simulated galaxies in order to connect galaxy structures with their formation and evolution. Having a consistent method to measure dynamical structures will allow us not only to directly compare across simulation outputs and observations, but to also shed light on the connection between the orbital properties and the assembly histories of galaxies. For example, Zhu et al. (2022a) focused on the orbital distributions of galaxies in the IllustrisTNG (Pillepich et al. 2018; Springel et al. 2018; Nelson et al. 2019) and eagle (Crain et al. 2015; Schaye et al. 2015) simulations, identifying hot inner stellar halo components in a way that is also consistent and achievable with deep IFS observations. They show that IFS observations and dynamical models of the inner regions of galaxies can provide a way to quantitatively determine the mass and time of ancient massive mergers (Zhu et al. 2022b).

In this paper, we build on these studies and explore the effect of galaxy mergers on the orbital structure of galaxies that can be obtained for many galaxies with large IFS surveys, in order to understand whether stars form on particular orbits (e.g. in a thin or thick disc) or whether they migrate to their z ∼ 0 configuration. We employ the eagle simulations, as these have been shown to successfully produce galactic structures, such as bulge sizes, that are similar to observations in galaxies with M_⋆ ≥ 10¹⁰M_⊙ (Lange et al. 2016), as well as to reproduce the fundamental plane of elliptical galaxies (de Graaff et al. 2023) and the diversity of stellar rotation-to-pressure support observed in local Universe galaxies (Lagos et al. 2018b; van de Sande et al. 2019). We also explore how well Eagle reproduces current observational constraints on the orbital structures of galaxies from the local Universe.

This paper is organised as follows. In Section 2, we describe the sample of galaxies and the galaxy properties we use in this analysis. Section 3 presents the orbital fractions we derive for the eagle sample, comparing them with results from observations. The connection between the orbital fractions and a galaxy’s evolution history is discussed in Section 4. Our conclusions are given in Section 5.

2 DATA

2.1 The eagle simulation

The Eagle simulation suite (Schaye et al. 2015, hereafter S15; Crain et al. 2015) consists of a large number of cosmological hydrodynamic simulations with different resolutions, volumes, and subgrid physics models, adopting the Planck Collaboration XVI (2014) cosmology. S15 introduced a reference model, within which the parameters of the sub-grid models governing energy feedback from stars and accreting black holes (BHs) were calibrated to ensure a good match to the z = 0.1 galaxy stellar mass function, the size-mass relation of present-day star-forming galaxies and the black hole–stellar mass relation (see Crain et al. 2015 for details). In this work, we take the largest volume with the highest resolution presented in S15. Table 1 summarizes the parameters of this simulation.

A major aspect of the eagle project is the use of state-of-the-art sub-grid models that capture unresolved physics (i.e. processes happening below the resolution limit). The sub-grid physics modules adopted by Eagle are: (i) radiative cooling and photoheating (Wiersma, Schaye & Smith 2009a), (ii) star formation (Schaye & Dalla Vecchia 2008), (iii) stellar evolution and chemical enrichment (Wiersma et al. 2009b), (iv) stellar feedback (Dalla Vecchia & Schaye 2012), and (v) BH growth and active galactic nucleus (AGN) feedback (Rosas-Guevara et al. 2015).

The eagle simulations were performed using an extensively modified version of the parallel N-body smoothed particle hydrodynamics (SPH) code gadget-3 (Springel 2005). Among those modifications are updates to the SPH technique, which are collectively referred to as ‘Anarchy’ (see Schaller et al. 2015 for an analysis of the impact of these changes on the properties of simulated galaxies compared to standard SPH). We use SUBFIND (Springel et al. 2001; Dolag et al. 2009) to identify self-bound overdensities of particles within haloes (i.e. substructures). These substructures are the galaxies in eagle.

2.2 Merger trees

We identify mergers using the merger trees available in the eagle database (McAlpine et al. 2016). These merger trees were created using the d-trees algorithm of Jiang et al. (2014). Qu et al. (2017) described how this algorithm was adapted to work with eagle outputs. Galaxies that went through mergers have more than one progenitor, and for our purpose, we track the most massive progenitors of merged galaxies and define that as the main branch. The trees stored in the public database of eagle connect 29 epochs. The time span between snapshots can range from ≈0.3 to ≈1 Gyr. We use these snapshots to analyse the evolution of the orbital components in galaxies and the effect of galaxy mergers.

2.3 Stellar mass, luminosities, and the half-mass radius

The stellar mass of eagle galaxies is measured using a spherical aperture. This gives similar results to the 2D Petrosian aperture used in observational works, and provides an orientation-independent mass measurement for each galaxy. Following previous studies based on the eagle simulations, we take as stellar mass the mass of the stellar particles within a 30 pkpc (pkpc refers to proper kpc) aperture centred at the centre of potential (e.g. Schaye et al. 2015).

We compute the projected stellar half-mass radius for each galaxy as the projected physical radius enclosing half of the stellar mass in the subhalo, averaged over three orthogonal projections (referred to as r₅₀ hereafter).

To compute the r-band light, we model each stellar particle as a single stellar population (SSP) with the age and metallicity of the particle, using the EMILES population synthesis models (Vazdekis et al. 2016) and adopting a Chabrier (2003) initial mass function. This provides a spectrum for each stellar particle. The r-band light is then the result of the convolution of the spectra with the Sloan Digital Sky Survey r-band filter (Fukugita et al. 1996). We use these r-band luminosities to calculate r-band luminosity weighted properties of galaxies, which are more directly comparable to observational results from IFS surveys (e.g. Santucci et al. 2022, 2023).

2.4 eagle Sample Selection

We select galaxies from the eagle Ref-L100N1504 simulation, spanning a stellar mass range of |$10^{9} - 10^{12}\, \rm M_{\odot }$| at all snapshots between z = 0 and z = 2. We note that due to the resolution of the simulation, dynamical properties of galaxies are not reliable for stellar masses |$\lesssim 10^{9.5}\, \rm M_{\odot }$| (see Appendix A in Lagos et al. 2017, and Ludlow et al. 2023). Moreover, Ludlow et al. (2023) showed that a stellar mass |$\gtrsim 10^{10}\, \rm M_{\odot }$| is needed for stellar kinematics to be resolved at the 3D half-mass radii of eagle galaxies; at lower masses, gravitational scattering between stellar and DM particles affects their kinematics and spatial distribution. It is important to note that DM mass resolution and gravitational scattering can also affect the circularity distributions of stellar disc particles by transforming galactic discs from flattened structures into rounder spheroidal systems, causing them to lose rotational support in the process (Wilkinson et al. 2023). We present a comparison between the eagle ‘Reference’ simulation used here and the high-resolution DM simulation from Ludlow et al. (2023) in Appendix  A.

For completeness, we present our results down to stellar masses of |$10^9\, \rm M_{\odot }$| (13199 galaxies), but we highlight the regime in which our results are likely affected by resolution. When looking at the connection between orbital fractions and mergers, we focus on galaxies with stellar masses |$\gtrsim 10^{10}\, \rm M_{\odot }$| (3561 galaxies), since this is the regime where kinematic quantities can robustly be measured in eagle.

During our analysis we often distinguish between galaxies that have or have not had mergers. To reliably trace minor and major mergers for all galaxies with stellar masses |$\gt 10^9\, \rm M_{\odot }$|⁠, we study all potential progenitors with stellar masses |$\ge 10^8\, \rm M_{\odot }$|⁠. We make this distinction based on the merger history of galaxies from z = 2 to z = 0: galaxies with ‘no mergers’ will be those that did not have a minor or major merger during that period, with minor and major mergers being those with a stellar mass ratio between the second most massive and most massive progenitors between [0.1−0.3[ and [0.3,1], respectively. Note that we can trace mergers with mass ratios <0.1 but the sample would not be complete for all galaxies in our sample given the limitation of the progenitors needing to have a stellar mass |$\ge 10^8\, \rm M_{\odot }$|⁠. Hence, we opt to only track minor and major mergers.

2.5 Orbital decomposition of eagle galaxies

To get to the orbital distribution and decomposition of stellar particles in eagle galaxies, we first calculate the ratio between the specific angular momentum of a stellar particle relative to the specific angular momentum of a circular orbit of the same energy, λ_z (a.k.a. circularity parameter). For this, we follow the method introduced by Abadi et al. (2003) described below:

• We first define the z-axis of a galaxy as that parallel to the stellar angular momentum vector of the galaxy computed with all stellar particles of the subhalo with the spatial reference frame being centred on the position of the most bound particle, and the velocity reference being the mean mass-weighted velocity of all stellar particles in the subhalo.

• Compute the specific angular momentum of each stellar particle in the z direction, j_z.

• We then compute the potential and kinetic energy of all stellar particles in the subhalo and bin particles by their total energy, with bins containing 100 particles.

• The specific angular momentum of a circular orbit with the total energy of a given stellar particle, j_c, is then taken as the maximum value of the specific angular momentum of all particles j_max among the 100 particles in the bin of total energy of that stellar particle.

• We calculate a λ_z = j_z/j_c for each stellar particle.

With the calculation above, we construct a two-dimensional (2D) distribution of λ_z versus projected distance (in the |$\rm xy$| plane of the galaxy) for each galaxy. We have verified that the resulting orbital distributions are insensitive to reasonable variations in the number of binding energy bins used to calculate j_c. The top-left panel in Fig. 1 shows examples of these distributions for three eagle galaxies at z = 0, revealing the diversity of orbital structures that galaxies can have. We provide a detailed analysis of these example galaxies later in this section.

We note that for massive galaxies the central black hole may have an impact on the central orbital structure if its sphere of influence is comparable to r₅₀ (e.g. Gebhardt et al. 2003; Thomas et al. 2014, 2016; Thater et al. 2017; Frigo et al. 2021; Merrell et al. 2023; Wheeler et al. 2023). However, the central black holes in galaxies in our sample have spheres of influence that are well below their r₅₀ values with an average |$r_{\rm SOI}\, \sim$| 0.05 r₅₀, with a standard deviation of 0.015, so we do not expect any significant effect on our results.

We separate orbits into four different components as follows:

• Cold orbits selected as λ_z > 0.8, are those with near-circular orbits;

• Warm orbits selected as 0.25 < λ_z < 0.8, are somewhat disc like, but with a considerable amount of angular momentum out of the disc (z) plane;

• Hot orbits selected as −0.25 < λ_z < 0.25, are those with near radial orbits;

• Counter-rotating orbits are those with λ_z < −0.25.

These thresholds were adopted to follow the analysis of Zhu et al. (2018a), Jin et al. (2020), Santucci et al. (2022) and allow for easier comparison with observations.

The fractional stellar mass and r-band light in each orbital component are calculated within 1r₅₀, by summing the mass or r-band luminosity, respectively, of each stellar particle within 1r₅₀ that is in the corresponding orbital family, and dividing by the total stellar mass or r-band luminosity within 1r₅₀. We also compute these fractions for the whole galaxy, without imposing a radial distance restriction on particles (i.e. using all the stellar particles of the subhalo).

We show the orbital distributions for three typical example galaxies – one low-mass (A), one intermediate–mass (B), and one high–mass galaxy (C) – in Fig. 1 (with stellar masses at the top of each column). These galaxies were visually selected to be representative of galaxies in their respective mass bins. For each of them, we show the binned orbital distribution of the stellar particles in the top panels. The colour indicates the median stellar age of the particles in each bin. The middle panels show the radial distribution of λ_z of the orbits for stellar particles in three different stellar age bins, as labelled. The bottom panels show the probability density function (PDF) of λ_z for particles in three different age bins. Galaxies in the low-mass range have a large fraction of their stellar mass distributed in cold and warm orbits. The cold component is also the youngest. This is typical for galaxies with stellar masses |$10^{9.5} \lt M_{\star }/\rm M_{\odot } \lt 10^{10.3}$|⁠, which in eagle are generally supported by rotation and are actively forming stars in their disc. Intermediate-mass galaxies (⁠|$10^{10.3} \lt M_{\star }/\rm M_{\odot } \lt 10^{11}$|⁠) show a higher contribution from the hot component, but with the warm and cold components still present. These galaxies are still forming stars on average, as shown by the young population that we see concentrated in cold orbits. High-mass galaxies (⁠|$M_{\star }/\rm M_{\odot } \gt 10^{11}$|⁠) are dominated by a dynamically hot and old component, extending to their outskirts. They have little contribution from cold orbits. On average, for all masses, particles in cold and warm orbits are younger than the hot component. Fractions of orbits for these galaxies can also be seen in the left-hand panels of Fig. 2 as indicated by the labels.

3 THE STELLAR ORBITAL DISTRIBUTIONS IN GALAXIES IN THE LOCAL UNIVERSE

3.1 Stellar orbital families of Eagle galaxies at z = 0

We analyse the distribution of the fraction of orbits (mass-weighted) within 1 r₅₀ for all galaxies at z ∼ 0, as a function of stellar mass (at z ∼ 0) in the top panels of Fig. 2.

We find that the warm and hot orbital components dominate at all masses within 1r₅₀, while the median contribution of the cold and counter-rotating components is always below ∼30 per cent. At |$\log M_{\star }/\rm M_{\odot } \lesssim 10.4$|⁠, the fraction of hot orbits decreases with increasing stellar mass, while warm and cold components become more important with increasing stellar mass, meaning that more massive galaxies are more rotationally supported than galaxies at lower masses. However, there is a clear change in the contributions of the different types of orbits at around a stellar mass of |$\log M_{\star }/\rm M_{\odot } \approx 10.4$|⁠, where we see a peak in the fractions of cold and warm orbits and a minimum in the fractions of hot and counter-rotating orbits. The transition mass we find is consistent with the transition mass reported in Clauwens et al. (2018), where the fraction of ex-situ stars starts to become more dominant in the eagle simulations. This transition mass is similar to what has been reported in observations to be the mass above which galaxy mergers become the dominant channel of galaxy growth in observations (Kauffmann et al. 2003; Cappellari et al. 2013; Robotham et al. 2014) and eagle (Lagos et al. 2018a). This transition is more evident when we divide our sample into galaxies that have evolved without experiencing minor and/or major mergers in the last 10 Gyr (between z ∼ 0 and z ∼ 2; middle panels), and galaxies that have undergone mergers (right-hand panels). Galaxies with no recent mergers in their evolutionary history show fractions of cold and warm orbits that are steadily increasing, while the fractions of hot and counter-rotating orbits decrease with increasing stellar mass. We see evidence of a turning point at around |$\log M_{\star }/\rm M_{\odot } \sim 10.75$|⁠, but the number of galaxies with no mergers at those high stellar masses is very small (112 galaxies). Even though these galaxies have not experienced any recent minor or major mergers, they have undergone very minor mergers (mass ratio ≤ 0.1). Excluding all the galaxies that also had very minor mergers removes the change in trend that we see at around |$\log M_{\star }/\rm M_{\odot } \sim 10.75$|⁠. Note that for the sample of galaxies with mergers, we still see the increase of the contribution of cold orbits with stellar mass up to |$\log M_{\star }/\rm M_{\odot } \sim 10.5$|⁠, but the peak contribution is smaller than what we see for galaxies that have not experienced mergers. This is expected, if we assume that mergers have an average heating effect on cold orbits (this will be explored in detail in Section 4).

The trends we find when looking at the fractional contribution of different orbital families in the whole galaxies are qualitatively similar to what we find within 1r₅₀. However, the total fraction of cold orbits shows an overall higher contribution than those within 1r₅₀. This is not surprising since the outer regions of a galaxy are generally disc-dominated (as seen in Fig. 2).

In observations, stellar orbit distributions have only been derived explicitly for three large (N>100) samples of galaxies, in the CALIFA (Zhu et al. 2018a, b), MaNGA (Jin et al. 2020), and SAMI (Santucci et al. 2022, 2023) surveys, with the results being generally consistent with each other, as shown in Santucci et al. (2022). We show the orbital fractions derived for SAMI passive galaxies, as well as the results from this work in Fig. 3. We only select passive galaxies (with specific star formation rates |$\lt 10^{-11} \, \rm yr^{-1}$|⁠) from our eagle sample, stellar mass-matching it with the SAMI one. Since SAMI orbital fractions are r-band light-weighted, we derive r-band light-weighted orbital fractions for our sample of eagle galaxies, to ensure a robust comparison. We show the median values of the fraction of orbits within 1 R_e (derived from the r-band light) for four different mass bins reported in Santucci et al. (2022) (shown as triangles, with the 1σ scatter represented with errorbars), while median values of the light-weighted fractions in eagle are shown as squares, and the 1σ scatter in Eagle is represented by the filled regions. The trends of the fractions of orbits that we observe with stellar mass and λ_Re are consistent with SAMI within the uncertainties. The inclusion of star-forming galaxies in the eagle mass-matched sample only has the effect of increasing the scatter in the orbits, in particular in the cold orbits (due to more prominent disc components being added to the sample).

We show the median age of the stellar particles in each orbital component in Fig. 4, top panels. We find that median ages of different components show qualitatively similar trends with stellar mass; stellar particles are older with increasing stellar mass regardless of their orbital family. All ages increase with increasing stellar mass, so that the most massive galaxies have the oldest populations. When looking at the different components, stellar particles on cold orbits are generally the youngest, with the oldest particles being on counter-rotating or hot orbits (which show a very similar stellar mass dependence), at all masses. However, there are some differences depending on whether galaxies had mergers or not. For example, the difference in age between the counter-rotating and hot components is very small, but slightly larger in the sample of galaxies with no merger; stellar particles on cold orbits in galaxies with no mergers are slightly younger than those in galaxies with mergers (about 1 Gyr at |$\log M_{\star }/\rm M_{\odot } \sim 10.4$|⁠); the opposite is seen for counter-rotating orbits, where we see older stellar particles in galaxies with no mergers compared to those with mergers. This is not surprising since counter-rotating orbits in the sample with mergers likely result from the mergers themselves (see Section 4 for a demonstration of that), which can happen at later times in a galaxy’s life. Moreover, looking at the median ages of the total stellar particles in the galaxy (see Fig. B1, top panels, in Appendix  B), we find that the total cold component is generally younger than the other components, in particular at low stellar masses. This difference is not as prominent when we look at the ages within 1 r₅₀, indicating that continuing star formation has a larger impact in the outskirts of a galaxy than within 1 r₅₀. Galaxies with mergers have inner cold components with slightly younger ages than galaxies with no mergers, but similar trends with mass.

The concept of ‘downsizing’ in observations refers to more massive galaxies forming their stars earlier but also quicker than lower mass galaxies (Thomas et al. 2010). With this in mind, we also analyse the standard deviation of the ages of the stellar particles (σ Age) in each component, shown in Fig. 4, bottom panels. Even though the underlying distribution of the ages of the stellar particles is not always symmetric (discs, in particular, tend to have longer tails towards young ages), σ Age is a simple way of comparing the spread in the ages of the stellar particles in different components. We find that σ Age of the cold component shows no dependence on stellar mass, while the warm, hot, and counter-rotating components have lower values of σ Age with increasing stellar mass in agreement with the ‘downsizing’ scenario. The cold component almost always shows a broader range in ages compared to the warm component – indicating a more prolonged star formation history. Hot and counter-rotating components show very similar variations in ages, and display on average less spread in ages than what we see in the cold and warm orbital components at fixed stellar mass. We find that galaxies with no mergers (on the middle panels) show steadily decreasing σ Age for particles in all components with increasing stellar mass. Our results do not change if we consider log  Age and σ log  Age instead of Age and σ Age.

Galaxies that have experienced minor or major mergers (right-hand panels) have higher values of σ Age for the particles in the cold component, causing the increase in σ Age that we see in the whole sample (left-hand panel). Trends for σ Age within 1r₅₀ and for the whole galaxies are generally similar. The only noticeable difference is in the σ Age of the cold component at low stellar masses. Within 1r₅₀, the cold component has the highest values of σ Age at all masses. However, if we consider stellar particles in the whole galaxy (see Fig. B1, bottom panels, in Appendix  B), cold orbits have lower σ Age than hot and warm orbits at low stellar masses. Mergers seem to lead to slightly higher σ Age for all components.

3.2 Stellar orbital families of eagle galaxies across cosmic time

We now explore how the orbital structures of galaxies change with redshift. We trace the main progenitors of each z = 0 galaxy in our sample, and we calculate the fraction of their orbital components at z = 0.5, z = 1, and z = 2 the same way as we did for our z = 0 sample. We show the total fraction of orbits in each component as a function of the stellar mass that the galaxies have at z = 0 for different redshifts in Fig. 5, left-hand panels. Note though, that the fractions themselves are calculated relative to the stellar mass of the progenitor galaxy at the respective redshift. We only show the total fractions of orbits since the behaviour within 1r₅₀ is very similar to that in the whole galaxy. The idea of showing the fraction of orbits as a function of the stellar mass at z = 0 is to isolate the evolution in the orbital fractions that happens in the y-axes. We note that, if we consider all the galaxies in our sample, regardless of whether they survived to z = 0, we find similar trends, albeit with slightly weaker overall evolution. Similarly, if we split the sample between galaxies that by z = 0 had mergers or not, we see overall qualitatively the same differences discussed in the context of the z = 0 distributions (Figs 2 and 4), and hence here we limit to showing the full sample without splitting by merger history.

Galaxies with |$M_{\star }/\rm M_{\odot }\gt 10^{10.5}$| at z = 0, have a peak in the fraction of cold and warm orbits at z ∼ 1, with the most massive galaxies reaching their peak earlier (z ∼ 2), and then steadily decrease with time, presumably due to the merging activity in these galaxies which leads to an increase in the fraction of hot and counter-rotating orbits. In the case of hot and counter-rotating orbits, we obtain the opposite trends with time, where the orbital fractions are higher at z = 0 for high-mass galaxies. Most of the evolution in the fractions of orbits (cold, warm, hot, and counter-rotating) in massive galaxies (≳10^10.5 M_⊙, happens between z = 0.5 and z = 0, presumably due to their active merger history. The MUSE MAGPI survey (Foster et al. 2021) is the ideal survey to test these predictions, since it focuses on studying galaxy transformations between z = 0 and z ∼ 0.4 in massive galaxies.

For galaxies with |$M_{\star }/\rm M_{\odot }\lt 10^{10.5}$| at z = 0, we see that their fraction of cold orbits steadily increases from z = 2 to z = 0, albeit with little evolution between z = 0 and z = 0.5. This is reminiscent of the evolutionary pathways identified by Lagos et al. (2017) for galaxies that at z = 0 have high specific angular momentum in Eagle: most of their spin-up happens at z ≲ 1. Similar behaviours are also found for the fractions of warm orbits, albeit with an overall weaker evolution at |$M_{\star }/\rm M_{\odot }\lt 10^{10.5}$| than what we see for cold orbits. Hot and counter-rotating orbits, on the other hand, decrease with decreasing redshift. For these galaxies (with |$M_{\star }/\rm M_{\odot }\lt 10^{10.5}$| at z = 0) we see that most of the increase in their orbital fractions happens between z = 2 and z = 1, with little evolution at z < 1. In fact, these galaxies display almost no evolution in their orbital fractions from z = 0.5 to z = 0.

These results suggest that at z ∼ 2 only the progenitors of very massive z = 0 galaxies (⁠|$M_{\star }\gtrsim 10^{11}\, \rm M_{\odot }$|⁠) had their disc components in place. For galaxies with |$10^{10.3}\lesssim M_{\star }/\rm M_{\odot }\lesssim 10^{11}$|⁠, the disc primarily forms between z = 2 and z = 1, while for galaxies of stellar mass |$\lt 10^{10}\, \rm M_{\odot }$|⁠, discs are formed preferentially at z < 1.

We show the median age of the stellar particles in each component as a function of stellar mass at z = 0, z = 0.5, z = 1, and z = 2 in the middle panels of Fig. 5 and  σ Age for each component in the right-hand panels. At fixed redshift, the ages of the components increase with increasing stellar mass. The degree of increase with stellar mass changes with redshift, becoming steeper at lower redshifts. Similarly, at fixed stellar mass, the ages increase with decreasing redshift. The trends still hold if we normalize the ages of each component by the age of the universe at each redshift (as seen in Appendix  C, Fig. C1). Cold orbits are always the youngest component in galaxies at all redshifts, although very close to the ages of the other components at z = 2. We also note that the counter-rotating component is the oldest at all redshifts, closely followed by the hot component. The warm component is always younger than the hot component, but older than the cold one. This is consistent with observations that found evidence of older stars being on hotter orbits than young stars in present-day systems such as the Milky Way (e.g. Bird et al. 2013), M31 (e.g. Quirk et al. 2018), and in MaNGA (Shetty et al. 2020) and SAMI (Foster et al. 2023) galaxies. Interestingly, in our sample, this is true not only at all redshifts, but also at all stellar masses (as seen in the top panels of Fig. 4). This indicates that the formation of stars in the different orbital components happens in similar ways for all galaxies: the counter-rotating and hot components are the first to cease star formation while the cold component is the last to form. Whether this means that the components are already dynamically in place, since stellar particles can be added (or migrate) to different components due to mergers and galaxy interactions, is still an open question. We also see a decrease in the scatter of the ages going towards high stellar masses for all components but the cold orbits, where we see a fairly stellar-mass independent scatter in stellar ages. This is consistent with stars in more massive galaxies forming earlier and over a shorter timescale (i.e. ‘downsizing’; Thomas et al. 2010). For z ≳ 1, σ Age does not show any trend with stellar mass, while we see a decrease in the σ Age with increasing stellar mass for the warm, hot, and counter-rotating components at z = 0. Cold orbits do not show any decrease in σ Age with stellar mass, even at z = 0.

4 DISCUSSION

We have analysed the orbital distributions of galaxies at z ∼ 0, with stellar masses above |$\log M_{\star }/\rm M_{\odot } \sim 9.5$| in the eagle simulations. We divided each galaxy into four different orbital components, depending on their orbital circularity: a cold component (with close to circular orbits), a warm component, a hot component, and a counter-rotating component. We calculated the fraction of stellar mass and the median stellar age in each component. We find that changes in the orbital distributions are mostly correlated with the stellar mass of the galaxies, as previously seen in the literature (e.g. Zhu et al. 2018a; Jin et al. 2020; Santucci et al. 2022). We also find differences in the ages and orbital fractions depending on the galaxy merger history.

4.1 Number of wet and dry mergers

In order to constrain the correlation between galaxy mergers and a galaxy’s orbital properties, we show in the left-hand panel of Fig. 6 the median values of the orbital fractions within (solid lines) and outside 1r₅₀ (dotted lines) as a function of the number of minor and/or major mergers each galaxy has undergone since z = 2, for galaxies of intermediate stellar masses (⁠|$10.3 \lt \log M_{\star }/ \rm M_{\odot } \lt 10.8$|⁠). Lower and higher mass galaxies show similar global trends, but with different absolute values (see Figs D1 and D2 in Appendix  D). We see a decrease in the fractions of warm and cold orbits with increasing number of mergers, while the fractions of hot and counter-rotating orbits increase with increasing number of mergers. The decrease of cold orbits with the number of mergers is stronger when we consider the region outside |$1\, r_{\rm 50}$|⁠, instead of only what is enclosed within |$1\, r_{\rm 50}$|⁠. The conjecture is that mergers tend to disrupt cold and warm orbits, which results in an increase in the fraction of hot and counter-rotating orbits after a galaxy has experienced a merger. The effect is cumulative, as galaxies with the highest fractions of counter-rotating and hot orbits are those with a cumulative number of mergers ≥3. This result is consistent with a two-phase formation scenario (White 1980; Kobayashi 2004; Oser et al. 2010), where more massive galaxies (which are also the galaxies showing high fractions of counter-rotating and hot orbits), start forming stars at earlier epochs than less massive ETGs (De Lucia et al. 2006; Dekel et al. 2009). After star formation is quenched, a second phase of slower evolution is dominated by mergers, leading to the increase of the hot and counter-rotating fractions. This is also what was shown to happen in Eagle for slowly rotating galaxies by Lagos et al. (2022). They show that the vast majority of slow-rotating galaxies experience quenching before they suffer the major morphological transformation that leads to their slowly-rotating nature.

Ages of stellar particles as a function of the number of major and/or minor mergers are shown in the middle panel of Fig. 6. Particles in galaxies that have undergone a higher number of mergers have younger ages, independently of the component they belong to, with the cold component having the youngest ages overall. This is likely a consequence of any gas present being funnelled to the galaxy centre and leading to new episodes of star formation by the galaxy merger. In addition, accreted stars from the secondary galaxy are likely to be younger than those in the primary due to the underlying positive correlation between stellar age and mass. This is also visible by looking at the age dispersion (Fig. 6, right-hand panel), where we see an increase in σ Age for the stellar particles in each component, but not for cold orbits, with increasing number of mergers.

Different types of mergers during the second phase can have different effects on a galaxy’s structure and composition, leading to the formation of slow- or fast-rotating galaxies. Cappellari (2016) described how these galaxies can form through two main formation channels: fast-rotating galaxies start as star-forming discs and evolve through a set of processes dominated by gas accretion, bulge growth, and quenching. By comparison, slowly rotating galaxies assemble near the centre of massive haloes via intense star formation at high redshift (Naab, Johansson & Ostriker 2009; Clauwens et al. 2018), and evolve from a set of processes dominated by gas-poor mergers, resulting in more triaxial shapes. Gas-poor and gas-rich mergers can also have very different effects on the spin of galaxies (as shown by Naab et al. 2014 and Lagos et al. 2022, for example).

To explore how orbital fractions are influenced by different types of mergers, we classify the mergers into gas rich (wet) and gas poor (dry) using the ratio of gas to stellar mass of the merger (f_{gas, merger}, derived by Lagos et al. 2018a). The latter is the sum of the gas masses of the two merging galaxies divided by the sum of the stellar masses of the two galaxies. In this paper, we classify mergers with a gas fraction below 0.2 as dry or gas poor, and mergers with gas fractions above this value as wet or gas rich. We note that changing the threshold does not impact the results, since the distribution of the gas fraction is highly bimodal: dry mergers have a median gas fraction of ∼0.10, while wet mergers have a median gas fraction of ∼0.98.

We show the median values of the fractions of orbits as a function of the number of wet and dry mergers in Fig. 7. Note that galaxies that experienced both wet and dry mergers appear in both panels; selecting galaxies with dry mergers only or wet mergers only does not change the trends we see here. In general, the fractions of orbits follow the same qualitative trends we see in Fig. 6. However, the effect of wet mergers is weak and leads to a small change in the orbital fractions of up to 20 per cent. Their effect, however, appears to accumulate with the number of wet mergers. Galaxies having undergone ≥3 wet mergers have ∼15 per cent more stellar mass in hot and counter-rotating orbits than galaxies with =1 wet merger.

Dry mergers, on the other hand, seem to have a great impact on the amount fractions change. In particular, one dry merger is enough to see a large increase (50–100 per cent) in the fractions of hot and counter-rotating orbits (and a decrease in the fractions of warm and cold orbits). This is not only true within 1r₅₀, but also in the outer regions, where the effects are stronger. Interestingly, we see that in general galaxies have fractions of hot orbits above 0.4 and counter-rotating orbits above 0.2 only if they have experienced at least one dry merger in their evolution (since z = 2). Similarly, dry mergers seem to be the main channel to reduce the contributions from warm and cold orbits (with the fractions of warm orbits below 0.3 and cold orbits below 0.1). We do not show the effects of dry and wet mergers on the ages of the components here, since we do not see significant differences for the different components (we find slightly younger ages, for all components, with increasing number of wet mergers, and slightly older ages for the cold component with increasing number of dry mergers). We do not find significant differences when we look at the fractions of orbits as a function of minor and major dry mergers.

4.2 Ex situ fractions

The internal orbital structures of the galaxies in our sample are also consistent with a scenario where a massive galaxy’s evolution is dominated by galaxy mergers, leading to their average spin-down. Previous results from simulations suggest that stars on different orbits have different formation paths. The cold components are mostly young stars formed in situ, the warm component likely traces old stars formed in situ or stars being heated from cold discs via secular evolution, and a small fraction of the warm component stars could be accreted (Gómez et al. 2017; Park et al. 2021). The stars on hot orbits in the outer regions should mostly be accreted via minor or major mergers (Gómez et al. 2017; Tissera et al. 2017), while stars on hot orbits in the inner regions are predicted to have formed in situ at high redshift (Williams et al. 2014; Du et al. 2021; Yu et al. 2021). However, the origin of the hot orbits, particularly in galaxies that have not experienced mergers, are predicted to be different for low- and high-mass galaxies: in low-mass galaxies, stellar particles are mostly born hot; while in high-mass galaxies, they are born in cold/warm orbits, but heated by bar-like secular evolution. For example, Yu et al. (2023) found that particles on hot orbits in Milky Way-like galaxies in the FIRE-2 cosmological zoom-in simulations (Hopkins et al. 2018) are all born hot. This is an interesting aspect that requires further investigation. However, the simple assumptions made in Eagle regarding the structure of the interstellar medium, specifically the fact that gas does not cool down below 8000 K, and therefore by construction is bound to have a significant dispersion component in galaxies of stellar masses |$\lesssim 10^{9.5}\, \rm M_{\odot }$|⁠, together with resolution limitations, prevent us from studying how orbits could transition from cold to warm in low-mass galaxies.

We find that massive galaxies have a substantial amount of stellar mass on hot orbits (over 40 per cent and increasing with increasing stellar mass) in the inner regions (within |$1\, r_{\rm 50}$|⁠), with stellar ages >10 Gyr. We also have high fractions of hot orbits when considering the whole galaxy. We also find that the most massive galaxies have smaller ranges in the ages of their particles (i.e. stars formed at around the same time) in the inner regions, and that these galaxies, being less supported by rotation, also generally have older mean ages than less massive (more rotationally supported) galaxies.

To better constrain the correlation between orbital components and ex situ fractions, we show in Fig. 8 the fraction of orbits in different components in the whole galaxy as a function of stellar mass, divided into hexabins colour coded by the median of the total ex situ fractions (derived by Davison et al. 2020) for the galaxies included in the bin. A minimum of five galaxies is required for each bin. Low-mass galaxies (⁠|$\log M_{\star }/ \rm M_{\odot } \lt $| 10.5) have very small ex situ fractions (with typical values being less than 0.1), except for galaxies with hot orbital fractions ≳0.5 where the ex situ fraction increases to ≈0.36. The ex situ fractions have a sharp increase around |$\log M_{\star }/ \rm M_{\odot } \sim 10.5$|⁠, as previously seen in eagle (e.g. Clauwens et al. 2018). Fractions of orbits within 1r₅₀ show similar, but weaker, trends and are shown in Appendix  E. This is not surprising since Fig. 6 shows that mergers and galaxy interactions have a stronger effect on the outskirts of galaxies (also seen by Lagos et al. 2018a; Karademir et al. 2019).

We find that low-mass (⁠|$\log M_{\star }/ \rm M_{\odot } \lt $| 10.5) in situ dominated galaxies with negligible ex situ fractions are characterized by high fractions of cold and, in particular, warm orbits. In galaxies with larger masses (⁠|$\log M_{\star }/ \rm M_{\odot } \gt $| 10.5) the galaxy structure is altered by the accreted stars and gas from merger events in a way that is strongly dependent on stellar mass. It is clear, however, that there are correlations between the ex-situ fraction and the fraction of mass in different orbital families at fixed stellar mass for galaxies with stellar masses above |$\log M_{\star }/\rm \rm M_{\odot } \sim 10.6$|⁠. In these galaxies, higher ex situ fractions are associated with lower fractions of cold and warm orbits and higher fractions of hot and counter-rotating orbits. Since the fraction of orbits and the stellar mass of a galaxy are measurable quantities (using advanced modelling techniques), we can use them to help us to put constraints on the galaxy’s ex situ fraction.

Our results suggest that galaxies with |$\log M_{\star }/ \rm M_{\odot } \lt 10.5$| in eagle are rotationally supported, dominated by young (∼6 Gyr) stellar particles in cold and warm orbits, which are driven by in situ star formation (>90 per cent of the stellar particles in these galaxies are formed in situ). We also find a small fraction of galaxies in this stellar mass range where stellar particles formed ex situ contribute to enhance the fraction of hot-orbits (as shown in Fig. 6).

Above |$\log M_{\star }/\rm \rm M_{\odot } \sim 10.5$| we see a change in the contributions of the different components (Fig. 2), so that galaxies transform towards a more pressure-supported configuration (dominated by old stellar particles in hot orbits). This transformation seems to be generally driven by dry mergers (Fig. 7), which is associated with accretion from stars formed ex situ (Fig. 8).

4.3 Counter-rotating orbits

The establishment of counter-rotating orbits, which are generally composed of old stellar particles, is thought to happen at a much earlier stage of galaxy evolution, during which both prograde and retrograde mergers contribute nearly equally to galaxy bulge growth. However, the total fraction of counter-rotating orbits should decrease with redshift, when the galaxy halo effectively accretes ambient cold gas, which cools and eventually forms a coherent gaseous and stellar disc (see e.g. Driver et al. 2013; Pillepich et al. 2019). For high-mass galaxies (⁠|$\log M_{\star }/ \rm M_{\odot } \gt $| 10.5), we see an increase in the fractions of counter-rotating orbits (Fig. 5), generally due to dry mergers (as seen in Fig. 7).

However, we have cases (∼8 per cent) where galaxies with no minor or major mergers in the last 10 Gyr (since z ∼ 2) have fractions of counter-rotating orbits comparable to those of galaxies that have experienced mergers, pointing to the need for an alternative path for the formation of counter-rotating orbits. Previous numerical studies (e.g. Scannapieco et al. 2009; Clauwens et al. 2018; Garrison-Kimmel et al. 2018) and cosmological simulations (e.g. Park et al. 2019) found that galaxies sometimes develop counter-rotating discs due to the gas infalling in a direction misaligned with the existing co-rotating disc plane. Therefore we suspect a possible way to form these orbits is through inflows of cold gas at particular angles. To test our hypothesis, we select a sample of central galaxies at z ∼ 0 with stellar mass above |$\log M_{\star }/ \rm M_{\odot } \sim 10.1$|⁠. We then calculate the median value of the fraction of counter-rotating orbits for these galaxies and select galaxies that are in the top (highest fractions of counter-rotation) and bottom 25th (lowest fractions of counter-rotation) percentiles. We exclude satellite galaxies since their tidal interactions may contribute to the counter-rotating orbits. The mass threshold is selected in order to have kinematic quantities robustly measured. Selecting a different mass limit – for example |$\log M_{\star }/ \rm M_{\odot } \sim 10.6$| – does not change our conclusions. From these two populations (top and bottom 25th percentiles), we exclude all the galaxies that have not experienced minor or major mergers since z = 2. For the remaining galaxies in these two subsamples, we show the median evolution of the mass of their counter-rotating as a function of time (in Gyr), as well as their total stellar mass and the cosine of the inclination angle of the accreted cold gas (as derived by Jiménez et al. 2022), in Fig. 9. To provide an example of the variation of the alignment of the gas inflow and the mass in the counter-rotating component with time in individual galaxies, we show values for an example galaxy in each population (green squares for the fraction of counter-rotating orbits and plum stars for the cos θ_acc). A value of cos θ_acc close to 1 means that the inflow of gas is closely aligned with the rotation axis of the galaxy. We find that the misalignment between the inflowing gas and the disc’s rotation axis, especially in the early Universe, can have a strong impact on the formation of counter-rotating orbits. Both populations have, in general, similar net accretion rates since their median total stellar mass distributions show a similar increase with time (black dotted lines in Fig. 9). Inflows with inclination angles highly misaligned (cos θ_acc <0.8) from the rotation axis of the disc, in particular the early Universe, seem to support the formation of stellar particles in counter-rotating orbits, as we can see from the increase in the mass of the counter-rotating component (which is increasing at a similar rate to the total stellar component) with redshift (left-hand panel). On the other hand, galaxies that do not show an increase in the mass of their counter-rotating component with redshift, had inflows of cold gas closely aligned with the rotation axis of the disc (right-hand panel). It is possible that both properties (counter-rotating orbits and misaligned gas inflows) could be tracing some other mechanism; however, Jiménez et al. (2022) found that the vertical velocity dispersion of cold gas correlates most strongly with the specific gas accretion rate onto the disc as well as with the degree of misalignment between the inflowing gas and the disc’s rotation axis. This evidence supports misaligned accretion to be a main driver of gas kinematics.

It would be interesting to check if similar results can be confirmed by observations. Counter-rotating components formed at high redshift can be distinguished from those formed through mergers by looking at the ages of the components. However, ages cannot help us distinguish between counter-rotation due to mergers and due to misaligned inflows, since in both cases we would have counter-rotating components with relatively younger ages than the components formed at high redshift. One possible way to determine their formation could be to study the metallicity of these counter-rotating components, or to look at the globular cluster population in the galaxy: metal-rich (red) clusters are born from an initial central burst, while metal-poor (blue) clusters are brought in by the later accretion of less massive satellites (Searle & Zinn 1978; Beasley et al. 2002; Tonini 2013; Leaman, VandenBerg & Mendel 2013; Kruijssen et al. 2019; Beasley et al. 2018). If the galaxy had no mergers, most of its globular clusters would be red, metal-rich globular clusters. If the galaxy had a merger, there should be evidence of the presence of blue, metal-poor globular clusters.

5 CONCLUSIONS

We analysed the orbital distributions of galaxies at z ∼ 0 in the eagle simulations in order to shed light on the connection between their orbital components and their evolution history. We divided each galaxy into four different components, depending on their orbital circularity: a cold component (with close to circular orbits), a warm component, a hot component, and a counter-rotating component. We calculated the fraction of stellar particles and their median age in each component. We find that:

• Changes in the orbital distributions of galaxies are correlated with the stellar mass of the galaxies, so that galaxies have higher fractions of hot and counter-rotating orbits and lower fractions of cold and warm orbits with increasing stellar mass (Fig. 2). We also find small differences when considering galaxies that did not experience any minor and/or major merger since z ∼ 2 and galaxies that have undergone minor and/or major mergers.

• Stellar particles are older with increasing stellar mass regardless of their orbital family. Stellar particles on cold orbits are generally the youngest, with the oldest particles being on counter-rotating or hot orbits, at all masses. We also find small differences depending on whether galaxies had mergers or not. Stellar particles on cold orbits in galaxies with no mergers are slightly younger than those in galaxies with mergers, while the opposite is seen for counter-rotating orbits, where we see older stellar particles in galaxies with no mergers compared to those with mergers. σ Age of the cold component shows no dependence on stellar mass, while the warm, hot, and counter-rotating components have lower values of σ Age with increasing stellar mass (Fig. 4).

• Galaxies with |$M_{\star }/\rm M_{\odot }\gt 10^{10.5}$| at z = 0, have a peak in the fraction of cold and warm orbits at z ∼ 1, with the most massive galaxies reaching their peak earlier (z ∼ 2), and then steadily decrease with time. Opposite trends with time are found for hot and counter-rotating orbits, so that the orbital fractions are higher at z = 0 for high-mass galaxies. For galaxies with |$M_{\star }/\rm M_{\odot }\lt 10^{10.5}$| at z = 0, the fraction of cold and warm orbits steadily increases from z = 2 to z = 0, albeit with little evolution between z = 0 and z = 0.5. Hot and counter-rotating orbits, on the other hand, decrease with decreasing redshift. The ages of the components increase with increasing redshift. σ Age decreases going towards high masses for all but the cold orbits, where it stays the same (Fig. 5).

• The fractions of orbits are closely connected to the number of mergers a galaxy experiences. In particular, dry mergers seem to be the main channel to reduce the contributions from warm and cold orbits and to increase the contributions from hot and counter-rotating orbits (Figs 6 and 7).

• At fixed stellar mass, for stellar masses |$\log M_{\star }/ \rm M_{\odot } \gt 10.5$|⁠, galaxies with higher fractions of hot and counter-rotating orbits show higher fractions of ex situ particles (Fig. 8).

• Counter-rotating components have three possible channels for their formation:

  • at high redshift, during a phase where both prograde and retrograde mergers contribute nearly equally to galaxy bulge growth;

  • through mergers, in particular dry mergers (Fig. 7).

  • in the absence of mergers, inflows with inclination angles highly misaligned from the rotation axis of the disc, can lead to an increase in the fractions of counter-rotating orbits (Fig. 9).

ACKNOWLEDGEMENTS

GS and KH acknowledge funding from the Australian Research Council (ARC) Discovery Project DP210101945. GS acknowledges funding from the UWA Research Collaboration Awards 2022. CL has received funding from the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) through project number CE170100013. CF is the recipient of an Australian Research Council Future Fellowship (project number FT210100168) funded by the Australian Government. CL and CF are the recipients of ARC Discovery Project DP210101945. RMcD acknowledges funding support via an Australian Research Council Future Fellowship (project number FT150100333). AP is supported by the Science and Technology Facilities Council through the Durham Astronomy Consolidated Grant 2020-2023 (ST/T000244/1). KLP acknowledges support from the Australian Government Research Training Program Scholarship. GvdV acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 724857 (Consolidator Grant ArcheoDyn).

DATA AVAILABILITY

The data underlying this article are available through http://icc.dur.ac.uk/Eagle. Data products (orbital distributions) created for this article are available upon request to the (Eagle orbital distributions: data request) corresponding author.

Measurements from Schwarzschild models of SAMI galaxies are available SAMI Schwarzschild models: data request contacting the corresponding author.

References

APPENDIX A: eagle HIGH-RESOLUTION DARK MATTER COMPONENT

We show in Fig. A1 the distribution of λ_z considering the total galaxy for the reference run (Ref100, in red) and the high-resolution dark matter run (Ref50 HR DM, in blue), for different mass bins. In general, the distributions for the two runs are very similar, with only small differences for stellar masses 10 < log M_⋆/M_⊙ <10.5 and in the scatter of the cold orbits.

If we calculate the fraction of orbits for both runs (Fig. A2), we find good agreement between the median values of all the components. Note that the stellar mass used for this figure is the total subhalo stellar mass, not the stellar mass within 30 kpc used in Fig. 2.

APPENDIX B: AGES AND σ AGE OF THE TOTAL ORBITAL COMPONENTS

In this Section, we show the median values of Age and σ Age for the orbits within the whole galaxy as a function of stellar mass. We find that median Ages and σ Age of the different components in the whole galaxy show very similar trends to those within 1 r₅₀.

APPENDIX C: STELLAR ORBITAL FAMILIES OF eagle GALAXIES ACROSS COSMIC TIME, NORMALIZED BY THE AGE OF THE UNIVERSE

We show the fraction of orbits, age, and σ Age of the stellar particles in each component as a function of stellar mass at z = 0, z = 0.5, z = 1, and z = 2 in Fig. C1. Ages and σ Age of each component are normalized by the Hubble time (t_H) at each redshift. The trends we see here are similar to those found in Fig. 5, with different absolute values due to the normalization, so that at fixed redshift the ages of the components increase with increasing stellar mass. The degree of increase with stellar mass changes with redshift, becoming steeper at lower redshifts. Similarly, at fixed stellar mass, the ages increase with decreasing redshift. σ Age decreases going towards high masses for all but the cold orbits, where it stays the same.

APPENDIX D: ORBITAL FRACTIONS AND NUMBER OF MERGERS

In this section, we show the median values of the orbital fractions within and outside 1₅₀ as a function of the number of minor and major mergers each galaxy has undergone, for galaxies with low mass (9.5 < log M_⋆/M_⊙ < 10.3; Fig. D1), and high mass (log M_⋆/M_⊙ < 10.8; Fig. D2). Low-mass galaxies and high-mass galaxies show similar global trends to those of galaxies with intermediate mass (shown in Fig. 6), although with different absolute values.

APPENDIX E: EX SITU FRACTIONS WITHIN 1r₅₀

We show in Fig. E1 the fraction of orbits in different components within 1r₅₀ as a function of stellar mass, colour coded by their total ex situ fractions. Low-mass galaxies have very little contribution from ex situ fractions, which have a sharp increase above log M_⋆/M_⊙ ∼ 10.5−10.6. These trends are consistent with what we see when considering the whole galaxy (Fig. 8), although less strong.