<figure><img src="image_1.png" /></figure>

HAL
open science

# The impact of AGN on stellar kinematics and orbits in simulated massive galaxies

Matteo Frigo, Thorsten Naab, Michaela Hirschmann, Ena Choi, Rachel

Somerville, Davor Krajnovic, Romeel Davé, Michele Cappellari

► To cite this version:

Matteo Frigo, Thorsten Naab, Michaela Hirschmann, Ena Choi, Rachel Somerville, et al.. The impact of AGN on stellar kinematics and orbits in simulated massive galaxies. Monthly Notices of the Royal Astronomical Society, 2019, 489 (2), pp.2702-2722. 10.1093/mnras/stz2318. hal-03989608

HAL Id: hal-03989608

https://hal.science/hal-03989608

Submitted on 16 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_2.png" /></figure>

# The impact of AGN on stellar kinematics and orbits in simulated massive galaxies

Matteo Frigo<sup>1,★</sup> Thorsten Naab,<sup>1</sup> Michaela Hirschmann,<sup>2,3</sup> Ena Choi,<sup>4</sup>
Rachel S. Somerville,<sup>5,6</sup> Davor Krajnovic,<sup>7</sup> Romeel Davé<sup>8,9</sup>
and Michele Cappellari<sup>10</sup>

Accepted 2019 August 19. Received 2019 August 12; in original form 2018 November 27

## ABSTRACT

We present a series of 10 × 2 cosmological zoom simulations of the formation of massive galaxies with and without a model for active galactic nucleus (AGN) feedback. Differences in stellar population and kinematic properties are evaluated by constructing mock integral field unit maps. The impact of the AGN is weak at high redshift when all systems are mostly fast rotating and disc-like. After z ~ 1 the AGN simulations result in lower mass, older, less metal rich, and slower rotating systems with less discy isophotes – in general agreement with observations. 2D kinematic maps of *in situ* and accreted stars show that these differences result from reduced *in-situ* star formation due to AGN feedback. A full analysis of stellar orbits indicates that galaxies simulated with AGN are typically more triaxial and have higher fractions of x-tubes and box orbits and lower fractions of z-tubes. This trend can also be explained by reduced late *in-situ* star formation. We introduce a global parameter, $\xi_3$, to characterize the anticorrelation between the third-order kinematic moment $h_3$ and the line-of-sight velocity ($V_{avg}/\sigma$), and compare to ATLAS<sup>3D</sup> observations. The kinematic correlation parameter $\xi_3$ might be a useful diagnostic for large integral field surveys as it is a kinematic indicator for intrinsic shape and orbital content.

**Key words:** methods: numerical – galaxies: evolution – galaxies: formation – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

The kinematics of stars within galaxies are a powerful tool to understand a galaxy’s formation history. Since their dynamics are to good approximation collisionless (see Binney & Tremaine 1987), stars preserve information about their origin for a long time. Stars can either originate from the galaxy itself, forming *in situ* through the gravitational collapse of gas, or can get accreted from different galaxies during merger events. These two types of origin lead to

different dynamical properties: newly *in-situ*-formed stars follow the dissipational dynamics of the gas which formed them, and this can lead to the formation of orderly rotating stellar discs (e.g. Robertson et al. 2004). The process of accreting stars from other galaxies is instead purely collisionless, and tends to result in round dispersion supported systems (e.g. Naab, Johansson & Ostriker 2009). Because of this connection, many studies attributed the difference in dynamical properties of present-day galaxies to different stellar origins. The more massive and bright early-type galaxies (with absolute magnitude $M_V \lesssim -21.5$), whose stellar component has been for a significant part accreted, tend to have smaller angular momentum (Emsellem et al. 2011) and more

<footnote n="1"><sup>1</sup>Max-Planck-Institut fuer Astrophysik, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany</footnote>

<footnote n="2"><sup>2</sup>DARK, Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, 2100 Copenhagen, Denmark</footnote>

<footnote n="3"><sup>3</sup>Institut d' Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Department of Astronomy, Columbia University, New York, NY 10027, USA</footnote>

<footnote n="5"><sup>5</sup>Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA</footnote>

<footnote n="6"><sup>6</sup>Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Road, Piscataway, NJ 08854, USA</footnote>

<footnote n="7"><sup>7</sup>Leibniz-Institut fur Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote n="8"><sup>8</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="9"><sup>9</sup>Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa</footnote>

<footnote n="10"><sup>10</sup>Sub-Department of Astrophysics, Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>* E-mail: mfrigo@mpa-garching.mpg.de</footnote>
complex kinematics (Krajnović et al. 2011), while intermediate and low-mass galaxies ($M_v \gtrsim -21.5$), which have formed most of their stars in-situ, tend to be simple, flattened, fast-rotating systems (see Cappellari 2016 for a review). Using cosmological simulations, Naab et al. (2014) and Röttgers, Naab & Oser (2014) linked the present-day kinematics of galaxies to the type of galaxy mergers they experienced during their formation: minor or major, and with or without gas. They found that gas-rich formation histories are more likely to produce fast-rotating systems, while gas-poor galaxy mergers tend to lead to slow-rotating systems. They however also found major mergers with particular configurations which result in an increase of the galactic angular momentum, displaying a complex picture. This picture would however be incomplete without including the energy feedback from active galactic nuclei (AGN). The connection between AGN and their host galaxies has been subject of research for more than two decades. Soon after the discovery of supermassive black holes (SMBH) in the centres of early-type galaxies, correlations have been found between their mass and galactic properties such as galactic bulge mass and velocity dispersion (Dressler 1989; Kormendy 1993; Gebhardt et al. 2000). This connection has been in the focus of theoretical work, with the conclusion that the energy feedback from accreting black holes could be necessary to reproduce these scaling relations, as well as to obtain the correct masses and abundances of early-type galaxies in cosmological simulations (see e.g. Di Matteo, Springel & Hernquist 2005; Croton et al. 2006; Vogelsberger et al. 2014; Schaye et al. 2015 and reviews by Kormendy & Ho 2013, Somerville & Davé 2015 and Naab & Ostriker 2017). The impact of AGN might also go beyond affecting global properties. The cosmological simulations of Choi et al. (2015, 2017) showed that in low-redshift galaxies the fraction of stars that form in situ is much lower when including AGN feedback. This has strong repercussions on the morphological and kinematic properties of the host galaxy. Dubois et al. (2016) and Penoyre et al. (2017a) showed that including AGN feedback in cosmological simulations increases the abundance of slow-rotating systems. Bryan et al. (2012) showed that AGN feedback affects the orbital composition of simulated galaxies up to large radii.

In this paper we analyse a small sample of high-resolution ‘zoom’ simulations (see Fig. 1) for a more in-depth look at the impact of AGN feedback on the kinematic and stellar-population properties of galaxies, but also extending the analysis to higher-order kinematics and orbital structure. We compare our simulated galaxies to observations by mocking the images produced by integral field unit (IFU) spectrographs. These instruments collect a spectrum for each of their spatial pixels, so that one can observe the spatial distribution of spectrum-derived quantities, such as line-of-sight velocity, metallicity, and age. Recently a number of large galaxy surveys have been performed with IFU spectrographs, such as MaNGA (Bundy et al. 2015), SAMI (Croom et al. 2012), and CALIFA (Sánchez et al. 2012), resulting in the mapping of thousands of galaxies. The MUSE spectrograph (Bacon et al. 2010) also delivered detailed 2D maps of galactic properties (e.g. Emsellem, Krajnović & Sarzi 2014; Krajnović et al. 2018), including at high redshift (Guérou et al. 2016). For the study of AGN feedback, this means that there is a huge library of data that can be used to look for signatures of the effect of AGNs, and in this paper we want to understand the nature of generic signatures through cosmological simulations. We do this by running each of our simulations twice, once with and once without our AGN feedback implementation, in order to find and analyse the differences between the two cases.

In Section 2 we present the set of cosmological simulations analysed in this paper. In Section 3 we describe how our mock observational maps are created, and how the other values we present are calculated. In Section 4 we look at the effect of AGN feedback on one exemplary simulated galaxy, through our mock integral field maps. In Section 5 we analyse the full simulation sample to get an idea of the general impact of AGN feedback. In Section 6 we discuss and summarize our conclusions.

# 2 SIMULATION DETAILS

## 2.1 Cosmological ‘zoom’ simulations

In this work we analyse a set of 20 prototypical cosmological zoom simulations of massive galaxies for the impact of feedback from accreting SMBH. Each of our 10 initial conditions is simulated once with and once without the AGN feedback model. Throughout the paper the two cases will be labelled as AGN and NoAGN. The initial conditions for the simulations were constructed from a (100 Mpc)$^3$ dark matter only simulation with a WMAP3 cosmology (Spergel et al. 2007): $h = 0.72$, $\Omega_b = 0.044$, $\Omega_{dm} = 0.216$, $\Omega_\Lambda = 0.74$, $\sigma_8 = 0.77$, $n_s = 0.95$. All details on the construction of the zoom initial conditions are presented in Oser et al. (2010, 2012). The same initial conditions were used in e.g. Naab et al. 2014 and Hirschmann et al. 2012, 2013; Naab et al. 2014 and Hirschmann et al. 2017, but here we simulate at higher resolution. Dark matter particles have a mass of $m_p = 3.62 \times 10^6$ M$_\odot$ h$^{-1}$ and gas particles initially have mass of $m_p = 7.37 \times 10^5$ M$_\odot$ h$^{-1}$. The simulations are run from $z = 43$ to $z = 0$ with gravitational softening lengths of 0.2 kpc for gas, star, and black hole particles and 0.45 kpc for dark matter particles at the highest resolution level.

The simulation code is the same as the one used in Hirschmann et al. (2017). We used an improved version of GADGET-3 (Springel 2005), SPHGAL, which overcomes the numerical limitations of the classic smoothed particle hydrodynamics (SPH) implementation. All details of SPHGAL are given in Hu et al. (2014).

### 2.1.1 Star formation and feedback

The simulation code includes a model for the formation of stellar populations from gas particles, representing star formation and for feedback. The stellar populations provide thermal and kinetic feedback as well as metals to the interstellar medium. The chemical enrichment was originally described by Scannapieco et al. (2005, 2006) and later improved by Aumer et al. (2013) and Núñez et al. (2017) with an updated feedback model. Gas particles are stochastically converted into star particles depending on the density of the gas, in a way that reproduces the Kennicutt–Schmidt relation (Kennicutt 1998). To be eligible for conversion into stars, SPH particles need to have density higher than a threshold density value $\rho_{th}$, given by:

$$ \rho_{th} = \rho_0 \left( \frac{T_{gas}}{T_0} \right)^3 \left( \frac{M_{gas}}{M_0} \right)^{-2}, \qquad (1) $$

where $\rho_0 = 1 \text{ cm}^{-3}$, $M_0 = 1.3 \cdot 10^6 \text{ M}_\odot$, $T_0 = 6000 \text{ K}$, and $M_{gas}$ and $T_{gas}$ are the mass and temperature of the gas particle, respectively (Núñez et al. 2017). In each time-step the eligible particles are assigned a probability of converting into star particles, so that the time-averaged star formation rate is:

$$ \frac{d\rho_*}{dt} = \epsilon_{SFR} \frac{\rho_{gas}}{t_{dyn}}, \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Mock luminosity images of our sample of simulated galaxies, run without (top two rows) and with (bottom two rows) AGN feedback. All galaxies are viewed at an inclination of 60 degrees (where 90 would mean edge-on). Stars are colour-coded by V-band weighted age based on Bruzual & Charlot (2003). Several of the galaxies simulated without AGN feedback show the presence of young stellar discs, which is unusual for galaxies of this mass.</figcaption></figure>

where $t_{dyn} = \sqrt{4\pi G \rho_{gas}}$ is the dynamical time of the gas particle and $\epsilon_{SFR} = 0.025$ is the star formation efficiency (chosen as in Núñez et al. 2017 and Hirschmann et al. 2017). The newly created star particles are then treated as collisionless. Each particle represents a single stellar population assuming a Kroupa (2001) initial mass function, with a given age and the metallicity of the original gas particle. This stellar population then exerts feedback to the surrounding gas. This takes the form of type Ia and II supernovae and of winds from asymptotic giant branch (AGB) stars. Type II supernovae happen at a given time $\tau_{SNII} = 3$ Myr after the creation of the star particles. This is on the short end of the typical delay time distributions for type II supernovae. Type Ia supernovae and AGB winds are added every 50 Myr after the star particle creation. Each event provides momentum and thermal energy to the surrounding gas. The total feedback energy is given by:

$$E = \frac{1}{2} m_{\text{ejected}} v_{\text{out}}^2, \qquad (3)$$

where $m_{ejected}$ is the mass ejected by the stellar population and $v_{out}$ is the assumed ejecta velocity. These are determined depending on the mass, age, and metallicity of the particle and on the type of event. We assume $v_{out} = 4500$ km s$^{-1}$ for type Ia and type II supernovae and $v_{out} = 10$ km s$^{-1}$ for AGB stars. The ejecta mass is taken from Woosley & Weaver (1995) for type II and from Iwamoto et al. (1999) for type Ia supernovae. This energy and mass is then added to the surrounding gas both as thermal (heating) and as momentum feedback (pushing). The relative fraction depends on the density and distance between the supernova-undergoing stellar particle and the 10 neighbouring gas particles, mimicking the evolution of blast waves (a simplified version of the three-phase model adopted in Núñez et al. 2017; see Hirschmann et al. 2017). The feedback events also distribute metals to the surrounding gas. Eleven elements are tracked for every gas and star particle (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), and their abundances are used to compute the cooling rate of the gas with the yields from Karakas (2010), Iwamoto et al. (1999), and Woosley & Weaver (1995) for AGB winds and supernovae type Ia and II, respectively. All details can be found in Aumer et al. (2013), Aumer, White & Naab (2014), and Núñez et al. (2017).
## 2.1.2 AGN feedback

AGN feedback is represented through the model developed by Choi et al. (2012) and used in Choi et al. (2014, 2015, 2017). This model includes both a radiative and a kinetic (wind) component (see Naab & Ostriker 2017 for a discussion of alternative numerical implementations). While it does not explicitly include relativistic jets, this model has been proven to produce galaxies with realistic sizes, star formation rates, black hole mass relations, and in particular hot gas (X-ray) properties without separating into different feedback modes (Choi et al. 2012, 2014, 2015, 2017; Eisenreich et al. 2017; Choi et al. 2018).

Here we summarize the most important elements used for this study. Black holes are first seeded at the centre of haloes exceeding a mass of $10^{11}$ M$_{\odot}$, with an initial mass of $M_{BH} = 10^5$ M$_{\odot}$, and they are repositioned on the potential minimum of their host galaxy at every time-step. They can then grow either by merging with other black hole particles or by accreting neighbouring gas particles according to a modified Bondi–Hoyle–Lyttleton (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952) accretion rate:

$$ \dot{M}_{\text{BHL}} = \left\langle \frac{4 \pi G^2 M_{\text{BH}}^2 \rho}{(c_s^2 + v^2)^{3/2}} \right\rangle, \qquad (4) $$

where $M_{BH}$ is the mass of the supermassive black hole, and $\rho$, $c_s$, and $v$ are the gas density, speed of sound and velocity with respect to the black hole. The calculation is done on an individual particle basis and then averaged over the SPH kernel, which is denoted by the angle brackets. This reduces the dependency on the number of SPH particles (see Choi et al. 2012). We do not explicitly enforce an Eddington limit cap on the accretion rate, as it is self-regulated by the inclusion of the Eddington force. We also do not add a ‘boost’ factor $\alpha$ (see Choi et al. 2017). Of the gas particles which could be accreted, 90 per cent are re-emitted as a wind parallel to the angular momentum of the gas next to the black hole (see Ostriker et al. 2010). This simulates the broad-line winds commonly emitted by AGN (de Kool et al. 2001; Somerville & Davé 2015; Naab & Ostriker 2017). The remaining 10 per cent are accreted, increasing the mass of the black hole particle. The model also includes radiative feedback in two forms. There is an Eddington radiation pressure force, which depends on the accretion rate and represents low energy photons providing momentum to the gas isotropically. We are also representing the higher energy X-ray photons, using the formulae from Sazonov et al. (2005) for Compton scattering. This component provides both momentum and thermal energy to the gas. As in Hirschmann et al. (2017), our simulation code differs from the one used in Choi et al. (2017) in not including metallicity-dependent heating, which was shown to have negligible impact (Choi et al. 2017).

# 3 ANALYSIS METHODS

## 3.1 Voronoi-binned kinematic maps

We study our simulations by analysing a series of mock IFU maps (e.g. Cappellari et al. 2011a) of the kinematics, metallicities, and ages of the stellar populations of the simulated galaxies. These maps are constructed with a PYTHON code developed for this work, following the analysis presented in Jesseit et al. (2007, 2009), Röttgers et al. (2014), and Naab et al. (2014). The code is included in the publicly available PYGAD analysis package.<fnref n="1" /> Positions and

velocities of the simulated galaxies are centred on the densest nuclear regions using a shrinking sphere technique on the stellar component. In the AGN simulations we centre the galaxies on their central supermassive black hole particles, which we define as the most massive black hole particle within 1 kpc of the stellar density centre. We then calculate the eigenvectors of the reduced inertia tensor (Bailin & Steinmetz 2005) of all stellar particles within 10 per cent of the virial radius $R_{vir}$, and use them to align the galaxies’ principal axes with the coordinate systems, such that the x-axis is the long axis and the z-axis is the short axis. To mimic seeing effects, each star particle in the simulation is split into 60 ‘pseudoparticles’, which keep the same velocity as the original particle and the positions are distributed according to a Gaussian with $\sigma = 0.2$ kpc centred on the original position of the particle (see Naab et al. 2014). In projection, the pseudoparticles are mapped on to a regular 2D grid, with pixel size 0.1 kpc (at $z=0$). Adjacent bins of this grid are then joined so that each resulting spaxel has a similar signal-to-noise ratio. We do this using the Voronoi tessellation method presented in Cappellari & Copin (2003). For the maps presented in this paper, the target signal-to-noise level of the spaxels has been set such that each of them contains 25 000 pseudoparticles (~400 regular particles). This ensures good statistics and results in maps similar to the ones from modern integral field surveys. The Voronoi grid is then used to construct the plots of stellar kinematics, metallicity, and age shown in this paper. For the age and metallicity maps, the value of every spaxel is calculated through a mass-weighted sample average. For the kinematic maps, we construct a histogram of the line-of-sight (LOS) velocity distribution of each spaxel, with the bin size determined by the Freedman & Diaconis (1981) rule. We then follow the classic approach of van der Marel & Franx (1993), Gerhard (1993), and fit the LOS velocity histogram with a Gauss–Hermite function:

$$ f(V) = I_0 e^{-y^2/2} (1 + h_3 H_3(y) + h_4 H_4(y)), \qquad (5) $$

where $y = (V - V_{\text{avg}})/\sigma$ and $H_3$ and $H_4$ are the Hermite polynomials of third and fourth order:

$$ H_3(y) = (2\sqrt{2}y^3 - 3\sqrt{2}y)/\sqrt{6} \qquad (6) $$

$$ H_4(y) = (4y^4 - 12y^2 + 3)/\sqrt{24}. \qquad (7) $$

The four fitting parameters $V_{\text{avg}}$ (average velocity), $\sigma$ (velocity dispersion), $h_3$ (skewness of the distribution), and $h_4$ (kurtosis of the distribution) are the ones plotted in the four panels of the kinematic maps.

To characterize the specific angular momentum of our galaxies we also employ the dimensionless $\lambda_R$ parameter (Emsellem et al. 2007), defined as:

$$ \lambda_R = \frac{\sum_i F_i R_i |V_i|}{\sum_i F_i R_i \sqrt{V_i^2 + \sigma_i^2}}, \qquad (8) $$

where the sum has been carried out over the spaxels of the kinematic maps, and $F_i$, $R_i$, $V_i$, and $\sigma_i$ are the flux, projected radius, average LOS velocity, and velocity dispersion of each spaxel, respectively. By limiting the sum to bins within a certain radius $R$, it is also possible to evaluate the cumulative radial $\lambda_R$ profile for every galaxy. The values given in the tables and kinematic maps are all calculated within the projected half-mass radius $R_e$.

<footnote n="1">¹https://bitbucket.org/broett/pygad</footnote>
## 3.2 Higher order kinematics

The higher order moments of the LOS velocity distribution, $h_3$ and $h_4$, can provide additional information on the orbital structure of our galaxies. In rotating systems the $h_3$ parameter has been observed to be anticorrelated to the average LOS velocity $V_{avg}$, or more specifically to the $V_{avg}/\sigma$ ratio (Gerhard 1993; Krajnović et al. 2011; van de Sande et al. 2017; Veale et al. 2017). This anticorrelation indicates that the LOS velocity distributions typically have a steep leading wing and a broad trailing wing. Simple axisymmetric rotating stellar systems show this property due to projection effects – stars are typically on circular orbits and those with lower LOS velocities projected into each spaxel produce a broad trailing wing. The slope of this anticorrelation is then about ~0.1 (Bender, Saglia & Gerhard 1994). However, if the galaxy is more complex, i.e not axisymmetric, it can also contain stars orbiting around different axes or radial orbits. This can make the trailing wing broader, as these stars have lower LOS velocities. The slope of the anticorrelation would then be steeper, and in some slow-rotating galaxies it can become extremely steep (see e.g. van de Sande et al. 2017). If the group of rotating stars becomes subdominant the correlation between $h_3$ and $V_{avg}/\sigma$ can change sign and become positive. Here the few fast rotating stars create a broad leading wing in the LOS velocity distribution (Naab, Jesseit & Burkert 2006; Hoffman et al. 2009; Röttgers et al. 2014). This unusual property is typically seen in simulated gas-poor mergers (Naab & Burkert 2001; Naab et al. 2014).

We characterize this variety of behaviours with a global parameter indicating the slope of the relation between $h_3$ and $V_{avg}/\sigma$ for all spaxels of one galaxy. This definition is inspired by the finding in Naab et al. (2014) that different slopes indicate varying formation histories and by the improved empirical classifications of the SAMI and MASSIVE galaxy surveys (van de Sande et al. 2017; Veale et al. 2017). We define the kinematic correlation parameter $\xi_3$ as:

$$ \xi_3 = \frac{\langle h_3 V_{avg}/\sigma \rangle}{\langle h_3^2 \rangle} = \frac{\sum_i F_i h_{3,i} (V_i/\sigma_i)}{\sum_i F_i h_{3,i}^2}, \quad (9) $$

where the sum is calculated over each spaxel out to $R_e$ from the centre, and the angle brackets indicate flux-weighted averaging over spaxels. When $h_3$ and $V_{avg}/\sigma$ are correlated, this parameter estimates the inverse of the slope of the correlation to reasonable accuracy with a simple fraction of weighted sums; negative values indicate a negative correlation, while positive values indicate a positive one. This can be seen by assuming $\langle h_3 \rangle = 0$ and $\langle V_{avg}/\sigma \rangle = 0$ (which is justified if the kinematic maps are symmetric) and rewriting the definition of the parameter as:

$$ \xi_3 = \rho_{V/\sigma, h_3} \frac{\sigma_{V/\sigma}}{\sigma_{h_3}}, \quad (10) $$

where $\rho_{V/\sigma, h_3}$ is the Pearson correlation coefficient of $V_{avg}/\sigma$ and $h_3$, and $\sigma_{V/\sigma}$ and $\sigma_{h_3}$ are the dispersion values of the two parameters (all flux-weighted). If $h_3$ and $V_{avg}/\sigma$ are linearly correlated then $\rho = \pm 1$, and $\xi_3$ becomes exactly the slope of the correlation.

Fig. 2 shows an example of the $h_3 - V_{avg}/\sigma$ spaxel values within $R_e$ for four simulated galaxies with different LOS velocity distribution properties. The lines indicate the simple slope given by $h_3 = (1/\xi_3) \cdot V_{avg}/\sigma$. Purely rotating systems are expected to have $\xi_3 \sim -10$, while rotating systems with non-negligible fractions of different orbit types are expected to lie in the $-3 < \xi_3 < -6$ range. When there is no correlation, or when the slope is almost vertical (both of which are observed in slow-rotating galaxies), the value of $\xi_3$ comes close to zero. Because of this, $\xi_3$ can potentially be a good way of distinguishing different types of galaxies. van de Sande et al.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Four simulation examples for the relation between $h_3$ and $V_{avg}/\sigma$ for the Voronoi bin within $R_e$ with the corresponding $\xi_3$ values. The lines indicate $h_3 = (1/\xi_3) \cdot V_{avg}/\sigma$. The upper left-hand panel shows a typical axisymmetric fast-rotator (0204 AGN), with a large fraction of prograde z-tube orbits ($f_{z-tube}^{\text{pro}} = 60$ per cent). It has the most negative $\xi_3$ value. The more complex fast-rotator 0501 AGN (upper right-hand panel, $f_{z-tube}^{\text{pro}} = 26$ per cent) shows a steeper slope with a higher $\xi_3$. Slow rotators like 0215 AGN (bottom left-hand panel, $f_{z-tube}^{\text{pro}} = 13$ per cent) have $\xi_3$ close to zero. Unusual non-axisymmetric rotating systems like 0227 AGN ($f_{z-tube}^{\text{pro}} = 20$ per cent) have a weak positive correlation between $h_3$ and $V_{avg}/\sigma$ resulting in a slightly positive value of $\xi_3$ (bottom right-hand panel).</figcaption></figure>

(2017) have used best-fitting elliptical Gaussians with a maximum likelihood approach to characterize the slope of the relation, which is slightly more complicated than our procedure. Veale et al. (2017) perform linear least-squares fit to calculate the slopes directly. Using the inverse of the slopes highlights the difference between slow rotators and the slow rotators get values around zero. The $h_4$ parameter is known to relate to orbit anisotropy (Gerhard 1993; van der Marel & Franx 1993; Thomas et al. 2007), with negative values indicating the dominance of tangential orbits and positive values corresponding to radial orbits. We do not further analyse $h_4$ other than showing the projected maps.

## 3.3 Isophotal shape and triaxiality

Our analysis involves the calculation of photometric quantities, such as the ellipticity $\epsilon$ and the isophotal shape parameter $a_4/a$, using stellar mass as a proxy for stellar light. The ellipticity values are calculated by fitting the galaxy isophotes with ellipses. The isophotes are constructed as lines with constant stellar surface mass density. For each galaxy we use 10 isophotes between 0.25$R_e$ and $R_e$ and average their ellipticity values to obtain $\epsilon$. The $a_4/a$ parameter represents the deviation of the shape of the isophotes of the galaxy from a perfect ellipse. It is used to discriminate between galaxies with ‘boxy’ or ‘discy’ isophotes (Lauer 1985; Bender & Moellenhoff 1987). We calculate it by applying a Fourier expansion to the deviation of the actual isophotes from their best-fitting ellipse:

$$ \delta R(\theta) = R(\theta) - R_{\text{ell}}(\theta) = \sum_n (a_n \cos(n\theta) + b_n \sin(n\theta)), \quad (11) $$

where $\theta$ is the azimuthal angle (Jedrzejewski 1987). The first-, second-, and third-order coefficients are negligible if the ellipse is centred correctly and has the correct ellipticity and orientation angle. The fourth-order coefficient $a_4$, normalized to the zeroth coefficient $a$, represents the deviation of the isophote from a pure ellipse. A positive value of $a_4/a$ means that there is an excess of light along
the major axis of the ellipse, causing the real isophote to be more 'discy'. A negative value instead means that the shape of the real isophote is more 'boxy' (see e.g. Naab, Burkert & Hernquist 1999; Springel 2000).

We also look at the 3D shape of our galaxies by computing the triaxiality parameter:

$$T = \frac{1 - (b/a)^2}{1 - (c/a)^2}, \qquad (12)$$

where b/a and c/a are the ratios between the main axes. When $T = 0$ the galaxy is perfectly oblate, while when $T = 1$ the galaxy is perfectly prolate. We calculated the axis ratios through the reduced inertia tensor (Bailin & Steinmetz 2005) of all particles within the effective radius $R_e$:

$$\tilde{I}_{i,j} = \sum_{\text{Particles } k} m_k \frac{r_{k,i} r_{k,j}}{r_k^2}, \qquad (13)$$

where $m_k$ and $\tilde{r}_k$ are the masses and positions of the particles. If the eigenvalues of $I_{i,j}$ are $\tilde{a}^2 > \tilde{b}^2 > \tilde{c}^2$, the ratios of their square roots $\tilde{b}/\tilde{a}$, $\tilde{c}/\tilde{a}$ are related to the real axis ratios by:

$$b/a = (\tilde{b}/\tilde{a})^{\sqrt{3}} \quad \text{and} \quad c/a = (\tilde{c}/\tilde{a})^{\sqrt{3}}. \qquad (14)$$

The $\sqrt{3}$ exponent is an empirical fit accounting for the bias of the axis ratios towards unity due to measuring the reduced inertia tensor within a sphere (Bailin & Steinmetz 2005).

## 3.4 Orbit analysis

We analyse the orbital composition of each of our simulated galaxies following the approach of Jesseit, Naab & Burkert (2005) and Röttgers et al. (2014). This procedure starts by freezing the potential of the simulated galaxy at $z = 0$ and representing it analytically using the self-consistent field (SCF) method (Hernquist & Ostriker 1992): the density and potential are expressed as a sum of bi-orthogonal basis functions, which satisfy the Poisson equation. There are such multiple density-potential pairs (e.g. Vasiliev 2013 and Lilley et al. 2018). We used the one from Hernquist & Ostriker (1992), in which the zeroeth-order element is the Hernquist (1990) profile:

$$\rho_{000} = \frac{M}{2\pi a^3} \frac{1}{\frac{r}{a}\left(1+\frac{r}{a}\right)^3} \qquad (15)$$

$$\Phi_{000} = -\frac{GM}{r+a}, \qquad (16)$$

where a is the scale parameter of the Hernquist profile. Higher order terms then account for both radial and angular deviations. The maximum radial and angular orders are $n_{max} = 18$ and $l_{max} = 7$, respectively. In the AGN simulations, the central SMBH are represented with a point mass (with the same softening length as the original simulations), which is a simpler and more precise representation. This point mass potential is summed to the SCF potential calculated with all the other particles.

We then integrate the orbits of each stellar particle within this fixed analytical potential for about 50 orbital periods. This is enough for identifying the orbit type, but not so much that quasi-regular orbits diverge from regular phase-space regions forcing us to classify them as irregular. The integration is done with an eighth-order Runge-Kutta integrator. The orbit classification itself is then done using the code by Carpintero & Aguilar (1998), which distinguishes different orbit families by looking for resonances between their frequencies along different axes (up to $m < 29$, $n < 10$, if $m:n$ is the resonance).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Age distribution of star particles in the case study galaxy (0227) for the run with AGN feedback (orange) and the one without (blue). The top x-axis shows the corresponding redshift at which the stars have formed. Star formation proceeded at a similar rate up $z = 2$. Then it is rapidly terminated in the presence of AGN feedback. Without AGN feedback star formation continues all the way to $z = 0$.</figcaption></figure>

In this paper we consider three main families of orbits: z-tubes (orbits that rotate around the z axis), x-tubes (orbits that rotate around the x-axis), box orbits ($\pi$ boxes and boxlets), and irregular orbits. In addition to these, we computed the fraction of prograde z-tube orbits $f_{z-tube}^{pro}$, by only selecting z-tubes with angular momentum along the z-axis of the same sign as the overall galaxy. Our orbit classification scheme makes sure that the orbit type of each particle is solid by analysing three different sections of the integrated orbit independently. If the three segments are classified differently, then we consider the classification to have failed, and we exclude these particles from the analysis (but they are still counted for computing the fraction of the other orbit types).

# 4 A TYPICAL GALAXY SIMULATED WITH AND WITHOUT AGN FEEDBACK

Our study involves a small sample of 20 massive galaxies. As a test case, in this section, we first discuss the formation history, global galaxy properties, stellar kinematics, stellar age and metallicity, morphology, and redshift evolution for one prototypical galaxy.

## 4.1 Formation history and global properties

Galaxy 0227 is an early-type galaxy, with an effective radius of 4.0 kpc and a stellar mass of $2 \times 10^{11} M_{\odot}$ in the AGN case and $5 \times 10^{11} M_{\odot}$ in the NoAGN case. Its formation history is characterized by a major merger at redshift $z \sim 0.25$, with mass ratio of 1: 1.7 and 1: 1.2 in the NoAGN and AGN cases. The presence of AGN has a strong influence on the evolution after the merger. Fig. 1 shows a mock V-band image of this galaxy with and without AGN feedback. In the absence of AGN feedback (left-hand panel) the galaxy is still forming new stars in an extended disc. Instead, in the case with AGN feedback (right-hand panel) the system is spheroidal with a very old stellar population.

Fig. 3 shows the age distribution of stars in galaxy 0227 simulated with and without AGN feedback. The oldest stars (age > 10 Gyr) have very similar age distributions, with the bulk forming around $z \sim 2$. Towards lower redshifts, star formation gets quenched in the AGN case; a behaviour found in all our simulations. While in the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Edge-on 2D line-of-sight stellar kinematics ($V_{\text{avg}}-\sigma-h_3-h_4$ from left to right) of galaxy 0227 simulated without AGN feedback at $z=2$, $z=1$, and $z=0$ (from top to bottom). The maps show typical features of systems with a disc-like component: high LOS velocity in the mid-plane, dumb-bell shaped velocity dispersion with a suppression in the mid-plane disc region, anticorrelation of line-of-sight velocity, and $h_3$ negative $h_4$ along the disc. These features become strongest at $z=0$, when the disc is most prominent and can clearly be seen in the surface density contours (black lines).</figcaption></figure>

AGN case not many stars form after $z \sim 1$, in the NoAGN case star formation continues throughout the simulation, including a starburst at $z \sim 0.25$ during the major merger.

## 4.2 LOS kinematics

In order to identify features in the stellar kinematics originating from the impact of AGN feedback, we construct 2D maps visualizing kinematic properties, as detailed in Section 3.1. Specifically we show the stellar line-of-sight velocity, dispersion, and the higher order moments $h_3$ and $h_4$ in Figs 4 and 5 for galaxy 0227 without and with AGN feedback at $z=2$, $z=1$, and $z=0$. Initially (at $z=2$ and $z=1$) there are only moderate differences between the AGN and NoAGN simulations. The AGN and NoAGN galaxies (in brackets) have similar stellar masses of $M_* = 0.59 \times 10^{11} M_\odot$ ($M_* = 0.54 \times 10^{11} M_\odot$) at $z=2$, while at $z=1$ they are $M_* \sim 1.16 \times 10^{11} M_\odot$ ($M_* \sim 1.97 \times 10^{11} M_\odot$). The effective radii are $\sim 0.18 \text{ kpc}$ ($\sim 0.35 \text{ kpc}$) at $z=2$ and $0.95 \text{ kpc}$ ($1.53 \text{ kpc}$) at $z=1$. Down to $z=1$, the galaxies are supported by rotation. The average stellar line-of-sight velocities reach values of $\sim 200 \text{ km s}^{-1}$, and the velocity dispersion values around $300 \text{ km s}^{-1}$. The velocity increases only slightly from $z=2$ to $z=1$, but the rotating component becomes more extended for both cases. The $h_3$ parameter is anticorrelated with the LOS velocity – a typical signature for axisymmetric rotating systems (Krajnović et al. 2011; Naab et al. 2014).

The origin of this effect is explained in detail in Section 3.2, as well as in Naab & Burkert (2001), Naab et al. (2006), Röttgers et al. (2014), and Naab et al. (2014) in the context of idealized models, merger simulations, and cosmological simulations. At redshift $z=0$, the situation is markedly different. In the NoAGN case the rotation signatures are significantly enhanced. The LOS velocities reach up to $320 \text{ km s}^{-1}$ in an extended disc. The velocity dispersion map shows a dumbbell feature with reduced velocity dispersion in the mid-plane, which is a signature of an edge-on rotation-supported disc embedded in a dispersion-supported spheroidal component. This can be seen by the isophotes (see Section 4.4). The LOS velocity distribution is asymmetric with anticorrelated $h_3$ values. The $h_4$ map shows characteristic features of disc rotation (bottom right-hand panel of Fig. 4). In the central kpc region, $h_4$ is positive, indicating a more peaked LOS velocity distribution with extended wings towards velocities lower and higher than the systemic velocity, as individual pixels cover significant fractions of the stars’ orbits. At larger radii (in the mid-plane), $h_4$ becomes negative indicating coherent rotation with very weak tails towards high and low velocities. As $h_4$ is known to roughly correlate with the velocity anisotropy (Gerhard 1993; Thomas et al. 2007), a negative $h_4$ indicates that tangentially biased orbits are dominating, which is to be expected in a rotating disc. Kinematic maps of this kind are regularly found in observational surveys like ATLAS<sup>3D</sup> (Cappellari et al. 2011b), CALIFA (Sánchez et al. 2012), or SAMI (Croom
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Fig. 4 for galaxy 0227 simulated including AGN feedback. The kinematics is qualitatively similar to the case without AGN at $z = 2$ and $z = 1$. By $z = 0$ however the strong rotational signatures are gone, and the galaxy looks more like a typical slow-rotator without kinematic disc signatures: low average LOS velocity, high velocity dispersion, no $h_3$ anticorrelation signal, positive $h_4$.</figcaption></figure>

et al. 2012). They are, however, more common for less massive galaxies. It is very unlikely to observe an elliptical galaxy of this high mass with such a prominent fast-rotating disc. The kinematic galaxy properties are very different in the AGN case (Fig. 5). By $z = 0$, there are no signatures of a prominent rotating stellar disc, as the AGN feedback prevents further gas accretion and in-situ disc formation (see e.g. Brennan et al. 2018). The galaxy is slowly rotating at ~80 km s⁻¹ and dispersion dominated, with only weak features in the higher order moments. Interestingly, $h_3$ is positively correlated with $V_{\text{avg}}$ in the central part of the galaxy. This is rare for observed galaxies, but relatively common in the simulated remnants of gas-poor mergers (see Naab & Burkert 2001; Naab et al. 2006; Röttgers et al. 2014). This positive correlation must originate from a particular orbital distribution, which will be analysed in Section 4.6. Also a core with negative $h_4$ is still visible. Values for $h_4$ are positive in most of the map indicating radially biased orbits. All of the above features are typical properties of massive early-type galaxies.

## 4.3 Age and metallicity distribution at $z = 0$

Fig. 6 shows a comparison of the projected stellar age (top panels) and metallicity (bottom panels) distributions for the NoAGN (left-hand column) and AGN (right-hand column) simulation at $z = 0$. At low redshifts the properties of the systems differ the most.

In the NoAGN case there is a distinct young $\lesssim 4$ Gyr stellar disc embedded in an older 7–9 Gyr stellar bulge. A moderate positive age gradient towards younger ages away from the centre is visible. The disc appears as a flattened metal-enriched region in the mid-plane, pretty much following the isophotes. These features indicate ongoing disc-like star formation and metal enrichment since $z = 1$. This is also consistent with the stellar age distribution in Fig. 3. In the AGN case (right-hand panels of Fig. 6) the stellar population is older (~10 Gyr, see also Fig. 3), less metal enriched – due to less ongoing star formation – with a shallower metallicity gradient. There is a mild positive age gradient with younger ages in the centre caused by residual nuclear star formation. The origin of age and metallicity gradients will not be discussed further in this paper (see e.g. Hirschmann et al. 2015; Rodriguez-Gomez et al. 2016).

## 4.4 Redshift evolution of kinematic and photometric properties

In this subsection we look at the evolution of three global parameters, $\lambda_R$, $\xi_3$, and $a_4/a$, through the whole formation history of our case-study simulation. We first use $\lambda_R$ (equation 8) to quantify the redshift evolution of specific angular momentum in the AGN and NoAGN cases. Fig. 7 shows the redshift evolution of $\lambda_R$ from $z = 2$ to $z = 0$. After a tumultuous phase at high redshift caused by mergers, at $z = 1$ $\lambda_R$ settles at around 0.3–0.4 in both cases.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Voronoi binned maps of the (mass-weighted) average stellar age (top) and metallicity (bottom) for our case-study galaxy (0227), in the *NoAGN* (left-hand panel) and *AGN* (right-hand panel) AGN cases, at $z = 0$. Without AGN feedback the higher star-formation rate at low redshift produces an overall much younger system, especially in the mid-plane, where a young stellar disc forms. Higher star-formation rate also results in high metallicities. With AGN feedback (right-hand panels) the galaxy is instead very old and its average metallicity is lower, closer to what we expect from observed early-type galaxies with this stellar mass (see e.g. McDermid et al. 2015 for a statistical study, or Guérou et al. 2016 for an example of spatially resolved age/metallicity maps).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Evolution of $\lambda_R$ for galaxy 0227, in the *NoAGN* and *AGN* cases. The values are indicated by inserted velocity maps out to the effective radius (isophote) of the galaxies. A major merger at $z \sim 0.25$ (vertical dotted line) strongly reduces the angular momentum of both systems. The *NoAGN* galaxy is less affected and can quickly regain angular momentum due to gas accretion and star formation. The *AGN* galaxy is instead unable to form new stars and remains a slow rotator.</figcaption></figure>

At $z = 0.25$ the angular momentum drops because of the major merger described in Section 4.1; the vertical dotted line marks the beginning of this merger. The subsequent evolution diverges for the two cases. In the *NoAGN* simulation the system is more gas rich, and thus loses less angular momentum and even regains it after the merger. This is a typical feature of gas-rich mergers and follow-up gas accretion (see review by Naab & Ostriker 2017). In

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution of $\xi_3$ for galaxy 0227, in the *NoAGN* and *AGN* cases. The values are indicated by small velocity maps out to the effective radius (isophote) of the galaxies. Up to $z = 0.25$ the value of $\xi_3$ is constant for both simulations and has a value as expected for a rotating system. However, after a major merger at $z \sim 0.25$ (vertical dotted line), the value for the *AGN* galaxy shifts towards zero and mildly positive values. This indicates that the galaxy lost its rotational support.</figcaption></figure>

the *AGN* case the system is already gas poor, without significant star formation before the merger (see Fig. 3). The merger then reduces the angular momentum significantly. Qualitatively this process for gas-poor mergers is discussed in detail in Naab et al. (2014). By $z = 0$ the two systems have very different rotation properties with a $\lambda_R$ value typical of fast rotators in the *NoAGN* case and a slow rotator value in the *AGN* case. This impact of AGN feedback on the rotation properties of massive galaxies has already been reported by Dubois et al. (2013), Martizzi et al. (2014), and Dubois et al. (2016) for cosmological RAMSES adaptive mesh refinement simulations with different AGN feedback models. We therefore assume this to be a generic feature of AGN feedback.

The major merger also affects the higher order kinematic features. We quantify them using the parameter $\xi_3$ defined in equation (9) and plot it as a function of redshift, as shown in Fig. 8. From $z = 1$ to $z = 0.25$ the two simulations show again the same behaviour, with the same degree of anticorrelation between $h_3$ and $V_{\text{avg}}/\sigma$: $\xi_3 \sim -7.5$ in both cases. As discussed in Section 3.2, this value is typical for a system dominated by tangential orbits, but higher than the one expected from a purely rotational system ($-10$). This indicates that a small amount of other orbit types contributes to skew the LOS velocity distribution. The major merger at $z = 0.25$ again makes the two cases diverge. In the *NoAGN* case the overall $\xi_3$ value stays the same. In the *AGN* case instead $\xi_3$ drops to 0 and the orbital structure of the system is more dispersion-supported – the correlation between $h_3$ and $V_{\text{avg}}/\sigma$ becomes weaker. The sign of $\xi_3$ oscillates a bit, but then settles to a weakly positive value, meaning that $h_3$ has the same sign as $V_{\text{avg}}$ as already pointed out.

We investigate the evolution of the isophotal shape parameter $a_4/a$, obtained by fitting the galactic isophotes at every snapshot (see subSection 3.3), with the galaxy seen edge-on. An example of these isophotes can be seen in the black lines of Figs 4 and 5. We show the evolution of $a_4/a$ since $z = 2$ in Fig. 9. Unlike in the previous cases, the *AGN* and *NoAGN* cases are already different at $z = 1$. The *NoAGN* case has systematically higher values of $a_4/a$ – more discy isophotes. This difference would however not be as
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Evolution of the isophotal shape of galaxy 0227, in the AGN and NoAGN case, quantified by $a_4/a$. The markers are surface brightness maps cut along the effective isophote. The red line represents perfectly elliptical isophotes. 'Boxy' galaxies have negative, 'discy' galaxies have positive $a_4/a$ values. Without AGN feedback the formation of a prominent disc results in discy isophotes at all times, despite the major merger at $z \sim 0.25$ (vertical dotted line). The AGN galaxy instead loses its discyness after the merger because further star formation is suppressed by the AGN.</figcaption></figure>

pronounced if the galaxy was not seen from an edge-on perspective. The value scatters due to minor mergers but drops to negative values after the major merger at $z = 0.25$. This is the common feature of major mergers destroying previously existing disc structures (see Naab et al. 1999; Naab & Burkert 2003). Subsequently a new stellar disc forms and the $a_4/a$ value becomes strongly positive again. In the AGN case the galaxy already lost its discyness at high redshift, because of the suppressed inflow of high-angular-momentum star-forming gas, and it keeps its elliptical or mildly boxy isophotes to $z = 0$. The effect of mergers and AGN feedback on the isophotal shape points in the same direction as the effect on $\lambda_R$ and $\xi_3$.

## 4.5 Kinematics of the accreted and in-situ-formed stellar components

Our kinematic maps can be generated for different stellar components of the galaxy to shed light on their respective kinematic structure. One might use the stellar age to distinguish different components; we show this example in the appendix. Perhaps even more interesting though, is to separate stellar particles according to their origin: either accreted from another galaxy or formed in situ in the main progenitor following the accretion of gas. Due to their intrinsically different origin, we can expect these two components to show very different kinematic (and stellar population) signatures (see e.g. Naab et al. 2014). To classify stars as in situ or accreted, we trace stars in the galaxies throughout the simulation from $z = 2$ to $z = 0$, and label them as in-situ stars when they form within 10 per cent of the virial radius (see Oser et al. 2010). All the remaining stellar particles are labelled as accreted. In the case of galaxy 0227 the in-situ fraction is $f_{in-situ} = 0.50$ and $f_{in-situ} = 0.17$ for the NoAGN and AGN cases, respectively. The values for the other galaxies are shown in Table 1.

Figs 10 and 11 show the stellar kinematic maps obtained for the separated in-situ and accreted components, in the NoAGN and AGN cases, respectively. The accreted components (upper panels of Figs 10 and 11) exhibit a very high velocity dispersion in both cases, but also have considerable net rotation, especially in the

NoAGN case. This larger net rotation is probably caused by the potential being more oblate-shaped in the NoAGN simulation ($T = 0.41$). In the AGN case the galaxy has a very triaxial, almost prolate shape ($T = 0.86$), which hinders the amount of z-tube orbits (more on this in Section 4.6) causing less rotation. The different shape is likely caused by the $z = 0.25$ merger being more gas-poor in the AGN case.

The in-situ components are very different in the two cases. In the AGN case (lower panel of Fig. 11), the in-situ stars follow the same kinematics as the accreted ones. Almost all of these stars formed before the major merger at $z = 0.25$, which means that their original orbits have been scrambled, resulting in a dispersion-supported system. In the NoAGN case the number of in-situ-formed stars is larger, both before and after the major merger, and the corresponding kinematic maps are more complex. There are two distinct features. The first is an orderly fast-rotating disc in the mid-plane, with low-velocity dispersion, a shallow $h_3-V_{avg}/\sigma$ trend, and strongly negative $h_4$. The second is a slow-rotating bulge with high velocity dispersion and a much steeper trend with $h_3$. The first component is mostly made of young stars which formed after the $z = 0.25$ major merger, hence the orderly motion. The surrounding bulge is instead older. These stars formed in situ at $z > 0.25$, and their orbits have been scrambled because of the major merger, resulting in less rotation. As the very high velocity dispersion suggests, there is also a counter-rotating component in this bulge, which explains why this component has a smaller net rotation than the accreted stars in the same potential.

This analysis implies that in-situ-formed stars and accreted stars tend to have intrinsically different kinematics from one another, at least until a major merger event scrambles their orbits. AGN feedback can thus significantly alter the present-day kinematics of galaxies by limiting in-situ star formation at late times, 'freezing' the kinematics at the most recent major merger, but also by affecting the merger history of the galaxy and thus the accreted stellar component.

## 4.6 Orbit distribution

It is also of interest to directly study the distribution of stellar orbits, and how it is affected by AGN feedback. We classify star particles into four global orbit types: z-tubes (rotating around the z-axis), x-tubes (rotating around the x-axis, including inner and outer major axis tubes), boxes (including $\pi$-boxes and boxlets), and irregular orbits. Fig. 12 shows the fraction of these orbit families as a function of radius. In both the AGN and NoAGN cases box orbits dominate within the central region, while tube orbits become dominant in the outer parts. In the NoAGN case, the fraction of z-tube orbits is larger at almost all radii. This is expected given the very prominent disc that has formed at low redshift. In the AGN case the z-tube fraction drops from 65 per cent to 49 per cent. The fraction of box orbits is lower in the centre and higher in the outskirts. The biggest change is however in the fraction of x-tube orbits, which went from an overall 5 per cent to 17 per cent. The likely reason for this is that the potential of the AGN galaxy has a more prolate shape ($T = 0.86$, instead of $T = 0.41$ for the NoAGN case), allowing for this kind of orbits. This change in the relative fraction of different orbit families also explains the positive correlation between $h_3$ and $V_{avg}/\sigma$; the bulk of the LOS velocity distribution is made of x-tube, box, and retrograde z-tube orbits, and the prograde z-tube orbits add a high-velocity tail to it. The fraction of irregular orbits is very low ($< 5$ per cent) in both cases. For about 20 per cent of the stellar particles our orbit classification scheme failed (see Section 3.4), and they are not included in the figure.
<table><caption>**Table 1.** General properties of our sample of simulated galaxies.</caption><thead><tr><th>GalID</th><th></th><th>M<sub>*</sub>(10<sup>10</sup> M<sub>☉</sub>)</th><th>R<sub>e</sub></th><th>Avg. age (Gyr)</th><th>f<sub>in-situ</sub></th><th>ε</th><th>a<sub>4</sub>/a</th><th>T</th><th>λ<sub>R</sub></th><th>ξ<sub>3</sub></th><th>f<sup>pro</sup><sub>z-tube</sub></th></tr></thead><tbody><tr><td>0175</td><td>NoAGN</td><td>26.73</td><td>1.86</td><td>7.85</td><td>0.23</td><td>0.24</td><td>0.017</td><td>0.72</td><td>0.08</td><td>-0.05</td><td>0.22</td></tr><tr><td>0175</td><td>AGN</td><td>18.93</td><td>2.57</td><td>10.72</td><td>0.05</td><td>0.35</td><td>0.001</td><td>0.51</td><td>0.12</td><td>-6.58</td><td>0.34</td></tr><tr><td>0204</td><td>NoAGN</td><td>19.59</td><td>3.09</td><td>8.17</td><td>0.63</td><td>0.78</td><td>0.196</td><td>0.29</td><td>0.46</td><td>-7.07</td><td>0.66</td></tr><tr><td>0204</td><td>AGN</td><td>16.41</td><td>2.06</td><td>9.50</td><td>0.18</td><td>0.37</td><td>0.026</td><td>0.15</td><td>0.36</td><td>-8.19</td><td>0.60</td></tr><tr><td>0215</td><td>NoAGN</td><td>27.79</td><td>1.76</td><td>9.61</td><td>0.27</td><td>0.38</td><td>0.042</td><td>0.48</td><td>0.37</td><td>-5.14</td><td>0.53</td></tr><tr><td>0215</td><td>AGN</td><td>7.38</td><td>1.70</td><td>11.26</td><td>0.09</td><td>0.31</td><td>0.018</td><td>0.50</td><td>0.02</td><td>-0.79</td><td>0.13</td></tr><tr><td>0227</td><td>NoAGN</td><td>48.46</td><td>3.27</td><td>7.72</td><td>0.50</td><td>0.49</td><td>0.124</td><td>0.41</td><td>0.47</td><td>-8.10</td><td>0.68</td></tr><tr><td>0227</td><td>AGN</td><td>22.24</td><td>2.60</td><td>9.95</td><td>0.15</td><td>0.15</td><td>0.000</td><td>0.86</td><td>0.10</td><td>0.50</td><td>0.20</td></tr><tr><td>0290</td><td>NoAGN</td><td>26.32</td><td>2.95</td><td>8.55</td><td>0.57</td><td>0.71</td><td>0.112</td><td>0.28</td><td>0.52</td><td>-11.49</td><td>0.74</td></tr><tr><td>0290</td><td>AGN</td><td>12.67</td><td>2.57</td><td>10.45</td><td>0.25</td><td>0.38</td><td>0.045</td><td>0.66</td><td>0.06</td><td>-2.62</td><td>0.23</td></tr><tr><td>0408</td><td>NoAGN</td><td>10.82</td><td>1.84</td><td>8.44</td><td>0.53</td><td>0.24</td><td>0.020</td><td>0.37</td><td>0.45</td><td>-5.96</td><td>0.53</td></tr><tr><td>0408</td><td>AGN</td><td>11.11</td><td>1.68</td><td>9.06</td><td>0.44</td><td>0.43</td><td>0.055</td><td>0.13</td><td>0.28</td><td>-6.48</td><td>0.50</td></tr><tr><td>0501</td><td>NoAGN</td><td>6.80</td><td>1.74</td><td>10.71</td><td>0.10</td><td>0.33</td><td>0.029</td><td>0.22</td><td>0.42</td><td>-6.74</td><td>0.54</td></tr><tr><td>0501</td><td>AGN</td><td>8.25</td><td>1.93</td><td>11.22</td><td>0.10</td><td>0.30</td><td>0.018</td><td>0.44</td><td>0.10</td><td>-3.55</td><td>0.26</td></tr><tr><td>0616</td><td>NoAGN</td><td>8.61</td><td>1.29</td><td>8.87</td><td>0.31</td><td>0.38</td><td>0.055</td><td>0.09</td><td>0.04</td><td>-1.02</td><td>0.13</td></tr><tr><td>0616</td><td>AGN</td><td>4.56</td><td>1.53</td><td>11.07</td><td>0.06</td><td>0.35</td><td>0.010</td><td>0.36</td><td>0.35</td><td>-7.27</td><td>0.50</td></tr><tr><td>0664</td><td>NoAGN</td><td>8.04</td><td>1.15</td><td>9.50</td><td>0.36</td><td>0.41</td><td>0.031</td><td>0.37</td><td>0.44</td><td>-5.74</td><td>0.24</td></tr><tr><td>0664</td><td>AGN</td><td>7.23</td><td>1.38</td><td>10.51</td><td>0.20</td><td>0.43</td><td>0.052</td><td>0.09</td><td>0.32</td><td>-7.61</td><td>0.67</td></tr><tr><td>0858</td><td>NoAGN</td><td>9.95</td><td>1.54</td><td>6.65</td><td>0.68</td><td>0.36</td><td>0.045</td><td>0.40</td><td>0.41</td><td>-8.85</td><td>0.58</td></tr><tr><td>0858</td><td>AGN</td><td>8.43</td><td>1.58</td><td>8.15</td><td>0.59</td><td>0.26</td><td>0.006</td><td>0.83</td><td>0.25</td><td>-6.88</td><td>0.19</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Stellar kinematics of galaxy 0227 (NoAGN) separated into its accreted component (above) and its *in-situ* formed one (below). The overall *in-situ* fraction is 50 per cent. The two components have strikingly different kinematics. The accreted component is mainly pressure-supported, but also rotates fast. The *in-situ* component shows two distinct features: a fast-rotating disc in the mid-plane with low velocity dispersion, and a slow-rotating bulge with very high velocity dispersion. The disc feature formed after a recent major merger, while the surrounding bulge is older, and its originally rotational orbits have been scrambled by the merger.</figcaption></figure>

The effect of AGN feedback on stellar orbits was also investigated in Bryan et al. (2012), who found a sharp increase in the fraction of box orbits. We do not see this in galaxy 0227, where the suppression of z-tube orbits is balanced by the increase in x-tube orbits. We however do see a trend towards more box orbits in the other galaxies of our sample, which will be discussed (and compared to Bryan et al. 2012) in Section 5.3.

## 5 RESULTS FROM THE SIMULATION SAMPLE

So far we focused on a single, example galaxy. In this section we show more general results for all 20 galaxies in our sample. This analysis cannot reveal the statistical kinematic properties of quiescent galaxy populations from recent cosmological simulations (Dubois et al. 2016; Penoyre et al. 2017b; Lagos et al. 2018; Schulze
<figure><img src="image_13.png" /><figcaption>Figure 11. Stellar kinematics of the accreted (above) and in-situ formed (below) stars of galaxy 0227 simulated with AGN feedback. The overall in-situ fraction is 17 per cent. In both cases the kinematics are pressure-supported, as no star formation happened since the last major merger at $z = 0.25$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12. Radial frequency of three different types of orbits in our case study galaxy: z-tubes, x-tubes, boxes, and irregular orbits. The dashed line shows the NoAGN case, and the full line shows the AGN one. In the latter, the fraction of z-tube orbits drops considerably due to the suppression of disc formation, and the fraction of x-tube orbits increases due to the more triaxial potential.</figcaption></figure>

et al. 2018). Instead, we would like to highlight the detailed impact of AGN feedback on massive galaxies for a few individual systems simulated at higher resolution. Table 1 shows for each galaxy in our sample the stellar mass $M_*$, effective radius $R_e$, the average stellar age, the in-situ formed fraction, the ellipticity $\epsilon$, the isophotal shape $a_4/a$, the triaxiality parameter $T$, $\lambda_R$, $\xi_3$, and the fraction of z-tube orbits $f_{z-tube}$. The galaxy IDs refer to a $z = 0$ halo in the parent simulation (the same as in Oser et al. 2012; Naab et al. 2014,...) and are ordered by halo mass. In general our galaxies have a lower stellar mass with AGN feedback due to the quenching of star formation, while the effective radius (for a given stellar mass) increases due to less dissipation (e.g. Crain et al. 2015; Choi et al. 2018). There are however notable exceptions: galaxies 0408 and 0501 have a larger stellar mass with AGN feedback. In both these cases the

NoAGN simulations go through strong stellar feedback events after high redshift starbursts, which result in a lot of gas leaving the galaxy. In the AGN case this does not happen, possibly because star formation happens more gradually. This causes the gas to stay in the galaxy and over time form more stars. Despite this, both systems experience less late star formation with AGN feedback, and at $z = 0$ have older ages and lower $\lambda_R$ values. In the following sections we will look at the distribution of kinematic ($\lambda_R$, $\xi_3$, orbit families) and morphological ($a_4/a$, triaxiality) properties at $z = 1$ and $z = 0$, and how AGN feedback affects them.

## 5.1 Angular momentum

In Fig. 13 we plot the $\lambda_R$ parameter of the sample galaxies versus their ellipticity $\epsilon$ for the simulations without (NoAGN, left-hand panels) and with AGN (AGN, right-hand panels) at redshift $z = 1$ (top panels) and $z = 0$ (bottom panels). The location of edge-on projections are indicated by the velocity maps. The blue/orange shaded regions indicate the typical distribution of these systems for 50 random orientations (projection effects for $\lambda_R$ based on simulations are discussed in e.g. Jesseit et al. 2009; Naab et al. 2014; Lagos et al. 2018). They were obtained by rotating the galaxies according to random rotation matrices, which were generated with the Stewart (1980) algorithm. The red line separates slow and fast rotators following to the definition by Cappellari (2016). A galaxy is considered a slow-rotator when

$$ \lambda_R < 0.08 + 0.25 \epsilon \quad \text{with} \quad \epsilon < 0.4. \qquad (17) $$

The distribution of galaxies at $z = 1$ is similar between the AGN and NoAGN cases, with most galaxies being flattened fast rotators with $\lambda_R$ in the range $0.2 < \lambda_R < 0.4$. The ellipticity values are a bit higher in the NoAGN case ($0.3 < \epsilon < 0.8$) than in the AGN one ($0.3 < \epsilon < 0.6$), but qualitatively the two populations are very similar. By $z = 0$ many (7 out of 10) of the NoAGN galaxies are still fast rotators with a similar ellipticity distribution. This trend is in agreement with results for massive galaxy populations from cosmological box simulations without AGN feedback (Dubois et al.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** $\lambda_R$ as a function of ellipticity $\epsilon$ for the galaxies at redshift $z = 1$ (top panels) and $z = 0$ (bottom panels), simulated without (NoAGN, left-hand panel) and with (AGN, right-hand panel) AGN feedback. The edge-on values are indicated by velocity maps. The coloured contours indicate the distribution of our galaxies when they are seen from 50 random orientations each. The dark red line marks the limit between slow- and fast-rotators according to Cappellari (2016). With AGN feedback the systems become rounder and rotate more slowly at $z = 0$.</figcaption></figure>

2016). Instead, in the AGN case by $z = 0$ the galaxies have become rounder ($\epsilon < 0.4$) and more slowly rotating, with $\lambda_R$ no larger than ~0.35. More than half of the galaxies would be considered bona fide slow rotators even in their edge-on projections. As discussed earlier, the trend towards slower rotation with AGN feedback is caused by the suppression of late *in-situ* star formation (see Brennan et al. 2018 for a discussion of ejective and preventative AGN feedback), which in most cases significantly reduces rotation observed at $z = 0$. The effect is strongest for the largest and most massive galaxies in our samples (lower numbers, like 0227), which – without AGN feedback – develop massive fast-rotating disc structures (in the case of galaxy 0175, this young disc structure is on a different plane, and thus does not increase $\lambda_R$ significantly).

We find a correlation between $\lambda_R$ and $\epsilon$, at least for the fast-rotators: faster rotating galaxies tend to be more flattened. Most of our slow-rotating galaxies exhibit a relatively high ellipticity, which is a trend found in other simulation studies as well (Bois et al. 2010; Naab et al. 2014), and possibly due to resolution limits. An interesting case is galaxy 0616, which contradicts our expectations by being a slow-rotator when simulated without AGN feedback but turns into a fast-rotator when simulated with AGN feedback. What happens here? In the NoAGN case gas infall triggers a starburst that forms a disc that counter-rotates with respect to the rest of the

galaxy. This lowers the projected $\lambda_R$ value, but leaves a relatively high ellipticity. In the AGN case the gas is kept from forming this new disc and the galaxy retains most of the (projected) angular momentum of the older stellar component.

In Fig. 14 we plot the $\lambda_R$ radial profiles for all galaxies, at $z = 1$ and $z = 0$. Typically, the values increase from the centre until they reach an asymptotic value, usually within $R_e$. This is consistent with previously published simulation data, even though we are missing systems with dropping $\lambda_R$ profiles (Naab et al. 2014; Wu et al. 2014; Lagos et al. 2018). At $z = 1$ there is not much difference between the AGN and NoAGN galaxies, while at $z = 0$ galaxies simulated with AGN feedback show once again systematically lower $\lambda_R$ values, even among the fast rotators. Many galaxies that would be rotationally supported without AGN, become pressure-supported when an AGN is present. Overall, AGN feedback results in more slow-rotating and dispersion-supported galaxies in agreement with previous simulations (Dubois et al. 2016).

## 5.2 Higher order kinematics

As discussed in Sections 3.2 and 4.4, rotating galaxies are expected to have anticorrelated $h_3$ and velocity fields, but the degree of this anticorrelation depends on the orbital structure of the galaxy, and
<figure><img src="image_16.png" /><figcaption>**Figure 14.** $\lambda_R$ radial profiles of our galaxy sample at $z = 1$ (top) and $z = 0$ (bottom). The AGN galaxies (orange) evolve towards lower $\lambda_R$ values than their NoAGN counterparts (blue).</figcaption></figure>

we can employ our $\xi_3$ parameter to evaluate this for our sample. In Fig. 15 we plot $\xi_3$ as a function of $\lambda_R$ at $z = 1$ and $z = 0$. The edge-on values are plotted with velocity maps, while the contours represent the location of the sample in the $\xi_3 - \lambda_R$ plane for random orientations. At $z = 1$ all galaxies have a negative of $\xi_3$ and $h_3$ is anticorrelated with the velocity, as expected for fast rotators. This is also true for the two galaxies which are slow rotators (according to $\lambda_R$) at $z = 1$. At $z = 0$ the sample splits into two groups: slow rotators with low values of $\lambda_R$ tend to have $\xi_3 \sim 0$ (very steep correlation or no correlation), while all fast rotators have $\xi_3 < -3$ (negative correlation). The specific value of $\xi_3$ for the fast rotators depend on their orbital structure; the galaxies where a disc feature is particularly prominent (0204, 0227, and 0290 in the NoAGN case) have the lowest values, reaching about $\xi_3 = -11.5$. In other words, more flattened and simple rotating systems have a less steep correlation between $h_3$ and $V_{\text{avg}}/\sigma$ than fast rotators with more complex kinematics. A similar behaviour was also observed in real galaxies by Veale et al. (2017). This results in a weak correlation between $\xi_3$ and $\lambda_R$ for the fast rotators that was not present at $z = 1$ when the kinematics of the galaxies were overall simpler. The bimodality of slow and fast rotators in the $\xi_3 - \lambda_R$ plane is seen in both the NoAGN and AGN cases, but with AGN feedback the group of galaxies with $\xi_3 \sim 0$ is larger. A few galaxies (0175 NoAGN and 0227 AGN) have a positive value of $\xi_3$ at $z = 0$ when viewed from certain orientations.

If we compare these results with observational IFU surveys, we find a small discrepancy. In Fig. 16 we plot the $\xi_3$ values of galaxies from the ATLAS$^{3D}$ survey (Cappellari et al. 2011b),<fnref n="2" /> compared with the contours of our AGN simulations seen at random inclinations. The ATLAS$^{3D}$ values also include a re-extraction of the kinematics from the subset of galaxies in the SAURON survey originally presented in Emsellem et al. (2004). To compute $\xi_3$ and $\lambda_R$ for the ATLAS$^{3D}$ sample, we only considered spaxels with $\sigma > 120 \text{ km s}^{-1}$, since the Gauss-Hermite moments can only be extracted from the data when the galaxy velocity dispersion is well resolved by the spectrograph (e.g. Cappellari & Emsellem 2004). The distribution of $\xi_3$ values is similar between observations and simulations, and can be divided in two groups: slow rotators with $\xi_3 \sim 0$ and fast rotators with $-3 < \xi_3 < -10$. However, at given $\lambda_R$ the ATLAS$^{3D}$ galaxies seem to have lower $\xi_3$ (in absolute value) than the simulations. We believe there are at least three reasons for this difference. Several of the ATLAS$^{3D}$ fast rotators have strong bar features, which are not present in our sample of simulations. In their presence the kinematic maps often show a positive correlation between $V_{\text{avg}}$ and $h_3$ (Chung & Bureau 2004), causing $\xi_3$ values closer to zero or sometimes even positive. In Fig. 16 galaxies with clear bars have been highlighted, and they tend to have slightly larger $\xi_3$, although this does not seem to account for all ATLAS$^{3D}$ fast rotators with $\xi_3$ close to zero. It is possible that hidden or weak bars are present in the other galaxies too, affecting the $h_3$ values. Secondly, as previously mentioned, constraining the $h_3$ value of each spaxel is harder in observations. The selection of spaxels with $\sigma > 120 \text{ km s}^{-1}$ limits this problem, but does not eliminate it. This results in more noisy $h_3$ maps, which makes the $h_3 - V_{\text{avg}}/\sigma$ trend less tight, and thus moves the $\xi_3$ value of observed galaxies closer to zero. At equal $\sigma$, this effect is stronger for slower rotating galaxies, as their LOS velocity distribution have lower $h_3$ values. Lastly, our (AGN) sample consists of only 10 massive galaxies, all of which have relatively low $\lambda_R$ values. This could be coincidental, but could also be caused by the limitations of AGN feedback models, which are known to be difficult to control and in this case might struggle to produce true fast rotators. If our AGN simulations did explore the $\lambda_R > 0.3$ regime, we would however expect most of them to have $-10 < \xi_3 < -5$, typical of rotation-dominated systems, matching the observations.

## 5.3 Orbit distribution and connection to kinematic correlation

The fraction of prograde z-tube orbits $f_{z\text{-tube}}^{\text{pro}}$ of the galaxies in our sample, measured in the same way as in Sections 3.4 and 4.6, are listed in Table 1. We generally find slightly lower prograde z-tube fractions in the AGN sample compared to the NoAGN one. A trend in this direction was also found by Bryan et al. (2012), who similarly derived the fractions of different orbit types in cosmological simulations with and without AGN feedback, although they focused on larger haloes and at large radii. If we compare our simulations at the same radius, we however find a weaker dependence on AGN feedback. In Bryan et al. (2012) the median fraction of tube orbits within 0.1 $R_{\text{vir}}$ goes from 45 per cent without AGN feedback (in their fiducial REF simulation) to 20 per cent with AGN. In our sample it goes from 63 per cent to 48 per cent. Similarly, the box orbit fraction of Bryan et al. (2012) goes from 50 per cent to almost 80 per cent, while in our simulations it stays much lower, from 20 per cent without AGN to 25 per cent with. This difference is likely due to

<footnote n="2">²Available from http://purl.org/atlas3d</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 15.** $\xi_3$ versus $\lambda_R$ at $z=1$ (top panels) and $z=0$ (bottom panels), simulated without (left-hand panel) and with (right-hand panel) AGN feedback. The kinematic map markers indicate the values when the galaxy is seen edge-on, while the density contours indicate the distribution when our galaxies are seen through 50 random orientations each. At $z=1$ all galaxies have values of $\xi_3$ in the anticorrelation regime, typical of fast rotators, while at $z=0$ many galaxies have $\xi_3 \sim 0$ or in a few cases even positive, and this effect is stronger with AGN feedback.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** $\xi_3$ as a function of $\lambda_R$ for the galaxies from the ATLAS$^{3D}$ sample (circle markers), compared with our AGN simulations (orange contours; same as Fig. 15). The ATLAS$^{3D}$ galaxies are distinguished in slow-rotators (red) and fast-rotators (light blue) according to the Cappellari (2016) definition. At equal $\lambda_R$, observed fast-rotators seem to have smaller $\xi_3$ (absolute) values than the simulation, possibly because of more complex kinematic features (bars) and of more noisy $h_3$ measurements. Slow-rotators have $\lambda_R \sim 0$ and $\xi_3 \sim 0$ in both observations and simulations.</figcaption></figure>

their more massive halo sample; their mean dark matter halo mass is $6\ 10^{13}\ h^{-1}\ M_{\odot}$, meaning that they contain galaxy groups rather than single galaxies. This explains the larger box orbit fractions, as in-situ star formation, which favours tube orbits, is less important at this scale. It also explains why in the Bryan et al. (2012) sample stars and dark matter particles roughly have the same orbital structure, which would not be true in our sample. The difference could however also be due to the different simulation code used in this work, particularly the stellar and AGN feedback models.

The kinematic correlation parameter $\xi_3$ that we introduced is related to the orbital structure of galaxies. Here we would like to study this connection using our orbit analysis. In Fig. 17 we plot $\xi_3$ as a function of the fraction of prograde z-tube orbits within $R_e$, $f_{z-tube}^{\text{pro}}$, at $z=0$. The plotted $\xi_3$ values are the average for 50 random views of each galaxy, and the error bars mark the dispersion (negligible for galaxies 0616 NoAGN and 0215 AGN, because they have projected velocity maps close to zero; see the kinematic map markers in Fig. 15). Most galaxies with high values of $f_{z-tube}^{\text{pro}}$ have a $\xi_3 < -3$ as expected, and there is a rough correlation between the two quantities. The galaxy with the highest $f_{z-tube}^{\text{pro}}$ (0290 NoAGN) is also the one with the lowest value of $\xi_3$: $\sim -11.5$ when seen edge-on and $\sim -6.5$ when averaging between many different viewing angles. The reason for this is that when the system is dominated by orbits
<figure><img src="image_19.png" /><figcaption>**Figure 17.** $\xi_3$ as a function of the fraction of prograde z-tube orbits $f_{z-tube}^{\text{pro}}$ for our sample of simulated galaxies at $z=0$. The $\xi_3$ values of each galaxy are an average over 50 random views, and the error bars are their standard deviation. Galaxies with high $f_{z-tube}^{\text{pro}}$ tend to have $\xi_3 < -5$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** $\xi_3$ as a function of the triaxiality parameter for our sample of simulated galaxies at $z=0$. The $\xi_3$ values of each galaxy are an average over 50 random views, and the error bars are their standard deviation. There is a correlation between the two parameters: prolate galaxies have positive values of $\xi_3$, while oblate galaxies have negative values.</figcaption></figure>

that rotate (progradely) around the z axis, these stars form the bulk of the LOS velocity distribution, and all other orbit types make the $h_3$ signal stronger for that given $V_{\text{avg}}/\sigma$. When non-rotational orbits are dominating ($f_{z-tube}^{\text{pro}} \sim 0$), then $V_{\text{avg}}/\sigma \sim 0$ and consequently $\xi_3 \sim 0$.

A few galaxies (0175 NoAGN and 0227 AGN) have a positive correlation between $h_3$ and $V_{\text{avg}}/\sigma$ in large parts of their kinematic maps, resulting in a positive value of $\xi_3$. This is likely connected to the fact that these galaxies have a prolate potential. We investigate this by plotting $\xi_3$ as a function of the triaxiality parameter $T$ in Fig. 18. There seems to be a rough correlation between the two quantities in our sample. The most prolate galaxies ($T \sim 1$) have positive values of $\xi_3$, while almost all oblate galaxies ($T < < 1$) have negative values. An exception is galaxy 0616 NoAGN, which as already discussed is made of two counter-rotating components and looks like a 'fake' slow rotator.

This connection between morphology and kinematics likely arises because different potential shapes allow different kinds of orbits; specifically, x-tubes are more common in prolate potentials. We see this by plotting $\xi_3$ as a function of the fraction of x-tube orbits

<figure><img src="image_21.png" /><figcaption>**Figure 19.** $\xi_3$ as a function of the fraction of x-tube orbits $f_{x-tube}$ for our sample of simulated galaxies at $z=0$. Colours and markers are the same as in Figs 17 and 18. The $\xi_3$ values of each galaxy are an average over 50 random views, and the error bars are their standard deviation. Galaxies with high $f_{x-tube}$ tend to have $\xi_3 > 0$.</figcaption></figure>

$f_{x-tube}$ in Fig. 19. There is again a rough correlation, meaning that galaxies with higher $f_{x-tube}$ are more likely to display a positive correlation between $h_3$ and $V_{\text{avg}}/\sigma$ in their kinematic maps. This follows from the correlation between $f_{x-tube}$ and the triaxiality $T$, which has previously been observed in isolated (Jesseit et al. 2005) and cosmological simulations (Röttgers et al. 2014). It should however be noted that in a pure prolate system only x-tube orbits and box orbits are allowed, and if there is net rotation around the long axis $h_3$ and $V_{\text{avg}}/\sigma$ become anticorrelated again. We do not see this in our sample because none of our galaxies is dominated by x-tube orbits (at most $f_{x-tube} = 0.25$, for 0227 AGN).

## 5.4 Isophotal shape

In Fig. 20 we plot the $a_4/a$ parameter of all our galaxies versus their ellipticity $\epsilon$ at $z=1$ and $z=0$. Like for Figs 13 and 15, we also added contours to show the distribution of values for smaller inclinations. At $z=1$ the panels with and without AGN feedback look qualitatively very similar. All galaxies have discy isophotes when viewed edge-on. When viewing the galaxies from different points of view both the ellipticity and the $a_4/a$ values tend to become smaller. At $z=0$, the cases with and without AGN behave as expected. The NoAGN galaxies show systematically higher $a_4/a$ values, due to the formation of embedded stellar discs at low redshift. In the AGN case the $a_4/a$ values are lower, meaning that the isophotes are less discy and closer to elliptical. Even though we do not have a clearly boxy galaxy in our sample, two galaxies (0175 and 0227) have almost perfectly elliptical isophotes.

We also computed the 3D shape of our galaxies using the triaxiality parameter $T$, defined in Section 3.3. The values of $T$ for our galaxies are found in Table 1, or in Fig. 18. We found that with AGN feedback a bigger fraction of our galaxies (five out of ten, instead of two out of ten) has a triaxial or almost prolate shape ($T > 0.5$). A prolate shape is more common for massive ellipticals, as found in both observations (Tsatsi et al. 2017; Graham et al. 2018; Krajnović et al. 2018) and simulations (Li et al. 2018). Without AGN feedback more of our galaxies are oblate ($T \sim 0$) despite their larger mass, which makes them more similar to the significantly
<figure><img src="image_22.png" /><figcaption>**Figure 20.** *a4/a* versus galaxy ellipticity *ϵ* at *z* = 1 (top panels) and *z* = 0 (bottom panels), simulated without (left-hand panel) and with (right-hand panel) AGN feedback. The edge-on locations are indicated by the isophotal maps, while the density contours indicate the distribution of our galaxies when they are seen from 50 random orientations each. The black line indicates elliptical isophotes and separates boxy (*a*4 < 0) from discy (*a*4 > 0) galaxies. Galaxies with AGN feedback are rounder and have more elliptical – in one case even boxy – isophotes at *z* = 0.</figcaption></figure>

less massive fast rotators we observe (Krajnović et al. 2011; Cappellari 2016).

## 5.5 Kinematic properties of the accreted component

The differences we analysed between the *AGN* and *NoAGN* simulations are most easily understood as a product of the prevention of late *in-situ* star formation by AGN feedback. As pointed out in Section 4.5 however, this is not the full story.

In the top panel of Fig. 21 we plot the $\lambda_R$ parameter as a function of ellipticity $\epsilon$, for our *AGN* and *NoAGN* sample, considering only their accreted component. This is achieved by filtering out all the particles which formed *in situ* since $z=2$, prior to the kinematic analysis. The *AGN* galaxies have a similar distribution to the one of the full (accreted + *in situ*) systems (Fig. 13) and to what we expect from observational surveys. This is a consequence of these systems having formed in large part through accretion. The *NoAGN* galaxies are instead different from full counterparts having generally lower angular momentum and ellipticity. This is again expected, since the missing *in-situ* component is predominantly flatter and faster rotating. However, the accreted components of the *NoAGN* simulations are also different than the accreted *AGN* components; despite the lower $\lambda_R$ values, they are in the range of fast rotators.

This effect is also visible in the $\xi_3-\lambda_R$ plane: the accreted component of the *NoAGN* galaxies has different kinematics from the full galaxies and from their *AGN* counterpart, mostly lying in-between slow and fast rotators ($\lambda_R \sim 0.1, \xi_3 \sim -2$).

The fact that the accreted components of our two samples have different kinematics implies that the impact of AGN feedback is not limited to the suppression of *in-situ* star formation, but also affects the merger history of these galaxies, specifically making it more gas poor. Dryer mergers tend to further decrease the net angular momentum of galaxies and result in less oblate systems, which support a wider variety of orbits (most notably, x-tubes).

## 6 DISCUSSION AND CONCLUSIONS

From the analysis of these simulated galaxies emerges a clear picture, which confirms the previous studies on the subject and adds new insights. The energy output of AGNs heats up and pushes away the interstellar gas, effectively suppressing the *in-situ* formation of stars. This affects the kinematics and morphology of the systems with a stronger impact at later cosmic times, when the central black holes become more massive. In our simulations AGN feedback results in realistic early-type galaxy properties at $z=0$. From our detailed stellar assembly, stellar population, mock IFU, isophotal shape, and stellar orbit analysis we get the following general picture:
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Top: Angular momentum proxy $\lambda_R$ as a function of ellipticity $\epsilon$, for our NoAGN (blue) and AGN (red) samples at $z=0$, when considering only their accreted stellar component. Bottom: Same for $\xi_3$ versus $\lambda_R$. The fact that in both cases there is no overlap implies that the effect of AGN feedback on galaxy kinematics is not simply due to the prevention of late in-situ star formation.</figcaption></figure>

(i) The stellar kinematics of massive early-type galaxies is significantly affected by AGN feedback, as seen both in the mock observational kinematic maps and in the orbit analysis of our simulation. Without AGN feedback massive early-type galaxies would develop young fast-rotating stellar discs even at low redshift, giving them kinematic signatures typical of less massive fast rotators. With AGN feedback massive early-type galaxies are instead more likely to become slow rotators due to the suppression of late in-situ star formation, in agreement with previous studies (Dubois et al. 2013; Martizzi et al. 2014; Penoyre et al. 2017b; Lagos et al. 2018).

(ii) As shown in Fig. 7, AGN feedback's effect of lowering the rotational support of galaxies is more pronounced in, but not limited to, late major mergers. Apart for some cases where mergers can cause a spin-up of the galaxy thanks to a favourable orbital configuration (Naab et al. 2014), most of the time mergers tend to disrupt the orbits of stars, reducing the angular momentum of the galaxy. However, without AGN feedback the further accretion of gas can produce a new rotating stellar disc and make the galaxy recover its angular momentum. With AGN feedback the in-falling star-forming gas is heated up and blown away. The origins of this

mechanism lie in the different spatial and kinematic properties of in-situ-formed and accreted stars (Rodriguez-Gomez et al. 2016).

(iii) AGN feedback starts having a significant impact on the stellar angular momentum only after $z=1$, and is stronger for more massive galaxies. With some exceptions, like galaxy 0616 in our sample which without AGN feedback develops a counter-rotating core, having AGN feedback always decreases the angular momentum of the galaxies in our sample.

(iv) While the main way that AGN feedback affects galaxy kinematics is by suppressing late in-situ star formation, its impact goes beyond that. By lowering the gas mass fraction of mergers it also affects the dynamics of the accreted stellar component.

(v) We compute the ellipticity $\epsilon$ and the $a_3/a$ isophotal shape parameter and follow their evolution through cosmic time. By suppressing the formation of discs, AGN feedback makes galaxies less flattened and their isophotes significantly less discy (more elliptical or even boxy), especially when seen edge-on. Like for the angular momentum, this difference starts arising at $z \sim 1$, and its effect is again stronger for the most massive galaxies of our sample.

(vi) We introduce a new global parameter, the kinematic correlation $\xi_3$, to quantify the anticorrelation between the LOS-velocity and $h_3$ from 2D kinematic maps. Slow and fast rotators have different typical values of this parameter owing to their different orbital structures. AGN feedback pushes the $\xi_3$ value towards the slow-rotator regime ($\xi_3 \sim 0$, meaning a very steep anticorrelation between $V_{avg}/\sigma$ and $h_3$ or lack of such a correlation).

(vii) We perform a full orbit analysis for all simulated galaxies and find that systems with AGN feedback have a higher fraction of x-tube and box orbits and a lower fraction of z-tubes. This is consistent with them being more triaxial due to the lack of late in-situ star formation and the more stellar accretion dominated assembly history. We find that the $\xi_3$ parameter is well correlated to the fractions of prolate z-tubes and x-tubes, as well as with the triaxiality of the galaxy.

(viii) We compared the $\xi_3$ values of our simulations with observed galaxies from the ATLAS$^{3D}$ sample finding an interesting discrepancy. At equal $\lambda_R$, observed fast rotators seem to have values of $\xi_3$ closer to zero and sometimes even positive; this could be because many of these galaxies show bar features, which cause a positive correlation between $h_3$ and LOS velocity, and/or possibly because of noise in the observed $h_3$ values. Our AGN sample also lacks galaxies with high $\lambda_R$ values, which are instead very common in the ATLAS$^{3D}$ sample.

Even though slow-rotating galaxies could also form without AGN feedback through particularly gas-poor formation paths, our simulations suggest that AGN feedback might be essential to produce the observed amount of quiescent, slow-rotating and non-discy early-type galaxies. The impact of AGN on the rotation properties are in line with earlier studies using different AGN feedback models and simulation codes (Dubois et al. 2013; Martizzi et al. 2014; Penoyre et al. 2017b; Lagos et al. 2018). In this study we demonstrate that also higher order properties in the isophotal shape and line-of-sight kinematics, as well as the underlying orbital content, are significantly affected by accreting SMBH. The effects typically result in a better agreement with observations. The newly introduced kinematic correlation parameter $\xi_3$ might provide a useful diagnostic for large integral field surveys, as it is a kinematic indicator for intrinsic shape and orbital content. This study is not statistically complete nor can the assumed AGN feedback model
be considered as an accurate description of the process. We can just give a model perspective on the observable effect of processes eventually happening in nature.

## ACKNOWLEDGEMENTS

This research was supported by the German Federal Ministry of Education and Research (BMBF) within the German-South-African collaboration project 01DG15006 'Ein kosmologische Modell für die Entwicklung der Gasverteilung in Galaxien'. TN acknowledges support from the German Research Foundation (DFG) Excellence Cluster 'Universe'. MH acknowledges financial support from the European Research Council (ERC) via an Advanced Grant under grant agreement no. 321323-NEOGAL.

## REFERENCES

* Aumer M., White S. D. M., Naab T., Scannapieco C., 2013, MNRAS, 434, 3142

* Aumer M., White S. D. M., Naab T., 2014, MNRAS, 441, 3679

* Bacon R. et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. SPIE, Bellingham, p. 773508.

* Bailin J., Steinmetz M., 2005, ApJ, 627, 647

* Bender R., Moellenhoff C., 1987, A&A, 177, 71

* Bender R., Saglia R. P., Gerhard O. E., 1994, MNRAS, 269, 785

* Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton University Press, Princeton.

* Bois M. et al., 2010, MNRAS, 406, 2405

* Bondi H., 1952, MNRAS, 112, 195

* Bondi H., Hoyle F., 1944, MNRAS, 104, 273

* Brennan R., Choi E., Somerville R. S., Hirschmann M., Naab T., Ostriker J. P., 2018, ApJ, 860, 14

* Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

* Bryan S. E., Mao S., Kay S. T., Schaye J., Dalla Vecchia C., Booth C. M., 2012, MNRAS, 422, 1863

* Bundy K. et al., 2015, ApJ, 798, 7

* Cappellari M., 2016, ARA&A, 54, 597

* Cappellari M., Copin Y., 2003, MNRAS, 342, 345

* Cappellari M., Emsellem E., 2004, PASP, 116, 138

* Cappellari M., Emsellem E., Krajnović D., McDermid R. M., Scott N., 2011a, MNRAS, 413, 813

* Cappellari M. et al., 2011b, MNRAS, 413, 813

* Carpintero D. D., Aguilar L. A., 1998, MNRAS, 298, 1

* Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125

* Choi E., Naab T., Ostriker J. P., Johansson P. H., Moster B. P., 2014, MNRAS, 442, 440

* Choi E., Ostriker J. P., Naab T., Oser L., Moster B. P., 2015, MNRAS, 449, 4105

* Choi E., Ostriker J. P., Naab T., Somerville R. S., Hirschmann M., Núñez A., Hu C.-Y., Oser L., 2017, ApJ, 844, 31

* Choi E., Somerville R. S., Ostriker J. P., Naab T., Hirschmann M., 2018, ApJ, 866, 91

* Chung A., Bureau M., 2004, AJ, 127, 3192

* Crain R. A. et al., 2015, MNRAS, 450, 1937

* Croom S. M. et al., 2012, MNRAS, 421, 872

* Croton D. J. et al., 2006, MNRAS, 365, 11

* de Kool M., Arav N., Becker R. H., Gregg M. D., White R. L., Laurent-Muehleisen S. A., Price T., Korista K. T., 2001, ApJ, 548, 609

* Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604

* Dressler A., 1989, in Osterbrock D. E., Miller J. S., eds, IAU Symposium Vol. 134, Active Galactic Nuclei. Kluwer, Dordrecht, p. 217

* Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, 433, 3297

* Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, MNRAS, 463, 3948

* Eisenreich M., Naab T., Choi E., Ostriker J. P., Emsellem E., 2017, MNRAS, 468, 751

* Emsellem E. et al., 2004, MNRAS, 352, 721

* Emsellem E. et al., 2007, MNRAS, 379, 401

* Emsellem E. et al., 2011, MNRAS, 414, 888

* Emsellem E., Krajnović D., Sarzi M., 2014, MNRAS, 445, L79

* Freedman D., Diaconis P., 1981, Z. Wahrscheinlichkeitstheor. Verwandte Geb., 57, 453

* Gebhardt K. et al., 2000, ApJ, 539, L13

* Gerhard O. E., 1993, MNRAS, 265, 213

* Graham M. T. et al., 2018, MNRAS, 477, 4711

* Guérou A., Emsellem E., Krajnović D., McDermid R. M., Contini T., Weilbacher P. M., 2016, A&A, 591, A143

* Hernquist L., 1990, ApJ, 356, 359

* Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375

* Hirschmann M., Naab T., Somerville R. S., Burkert A., Oser L., 2012, MNRAS, 419, 3200

* Hirschmann M. et al., 2013, MNRAS, 436, 2929

* Hirschmann M., Naab T., Ostriker J. P., Forbes D. A., Duc P.-A., Davé R., Oser L., Karabal E., 2015, MNRAS, 449, 528

* Hirschmann M., Charlot S., Feltre A., Naab T., Choi E., Ostriker J. P., Somerville R. S., 2017, MNRAS, 472, 2468

* Hoffman L., Cox T. J., Dutta S., Hernquist L., 2009, ApJ, 705, 920

* Hoyle F., Lyttleton R. A., 1939, Proc. Camb. Phil. Soc., 35, 405

* Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, MNRAS, 443, 1173

* Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439

* Jedrzejewski R. I., 1987, MNRAS, 226, 747

* Jesseit R., Naab T., Burkert A., 2005, MNRAS, 360, 1185

* Jesseit R., Naab T., Peletier R. F., Burkert A., 2007, MNRAS, 376, 997

* Jesseit R., Cappellari M., Naab T., Emsellem E., Burkert A., 2009, MNRAS, 397, 1202

* Karakas A. I., 2010, MNRAS, 403, 1413

* Kennicutt R. C., Jr, 1998, ApJ, 498, 541

* Kormendy J., 1993, in Beckman J., Colina L., Netzer H., eds, The Near-est Active Galaxies. Consejo Superior de Investigaciones Científicas, Madrid, p. 197

* Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

* Krajnović D. et al., 2011, MNRAS, 414, 2923

* Krajnović D., Emsellem E., den Brok M., Marino R. A., Schmidt K. B., Steinmetz M., Weilbacher P. M., 2018, MNRAS, 477, 5327

* Kroupa P., 2001, MNRAS, 322, 231

* Lagos C. del P., Schaye J., Bahé Y., Van de Sande J., Kay S. T., Barnes D., Davis T. A., Dalla Vecchia C., 2018, MNRAS, 476, 4327

* Lauer T. R., 1985, MNRAS, 216, 429

* Li H., Mao S., Emsellem E., Xu D., Springel V., Krajnović D., 2018, MNRAS, 473, 1489

* Lilley E. J., Sanders J. L., Evans N. W., Erkal D., 2018, MNRAS, 476, 2092

* Martizzi D., Jimmy, Teyssier R., Moore B., 2014, MNRAS, 443, 1500

* McDermid R. M. et al., 2015, MNRAS, 448, 3484

* Naab T., Burkert A., 2001, ApJ, 555, L91

* Naab T., Burkert A., 2003, ApJ, 597, 893

* Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

* Naab T., Burkert A., Hernquist L., 1999, ApJ, 523, L133

* Naab T., Jesseit R., Burkert A., 2006, MNRAS, 372, 839

* Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178

* Naab T. et al., 2014, MNRAS, 444, 3357

* Núñez A., Ostriker J. P., Naab T., Oser L., Hu C.-Y., Choi E., 2017, ApJ, 836, 204

* Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312

* Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63

* Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642

* Penoyre Z., Moster B. P., Sijacki D., Genel S., 2017a, MNRAS, 468, 3883

* Penoyre Z., Moster B. P., Sijacki D., Genel S., 2017b, MNRAS, 468, 3883
Robertson B., Yoshida N., Springel V., Hernquist L., 2004, ApJ, 606, 32

Rodriguez-Gomez V. et al., 2016, MNRAS, 458, 2371

Röttgers B., Naab T., Oser L., 2014, MNRAS, 445, 1065

Sánchez S. F. et al., 2012, A&A, 538, A8

Sazonov S. Y., Ostriker J. P., Ciotti L., Sunyaev R. A., 2005, MNRAS, 358, 168

Scannapieco C., Tissera P. B., White S. D. M., Springel V., 2005, MNRAS, 364, 552

Scannapieco C., Tissera P. B., White S. D. M., Springel V., 2006, MNRAS, 371, 1125

Schaye J. et al., 2015, MNRAS, 446, 521

Schulze F., Remus R.-S., Dolag K., Burkert A., Emsellem E., van de Ven G., 2018, MNRAS, 480, 4636

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Spergel D. N. et al., 2007, ApJS, 170, 377

Springel V., 2000, MNRAS, 312, 859

Springel V., 2005, MNRAS, 364, 1105

Stewart G. W., 1980, SIAM J. Numer. Anal., 17, 403

Thomas J., Saglia R. P., Bender R., Thomas D., Gebhardt K., Magorrian J., Corsini E. M., Wegner G., 2007, MNRAS, 382, 657

Tsatsi A., Lyubenova M., van de Ven G., Chang J., Aguerri J. A. L., Falcón-Barroso J., Macciò A. V., 2017, A&A, 606, A62

van de Sande J. et al., 2017, ApJ, 835, 104

van der Marel R. P., Franx M., 1993, ApJ, 407, 525

Vasiliev E., 2013, MNRAS, 434, 3174

Veale M. et al., 2017, MNRAS, 464, 356

Vogelsberger M. et al., 2014, Nature, 509, 177

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Wu X., Gerhard O., Naab T., Oser L., Martinez-Valpuesta I., Hilz M., Churazov E., Lyskova N., 2014, MNRAS, 438, 2701

## APPENDIX: ADDITIONAL FIGURES

Figs A1 and A2 show the kinematic maps of galaxy 0227 at $z=0$ separating the stars between different age groups: younger than 3 Gyr (formed after the major merger at $z=0.25$), between 3 and 10 Gyr old, and older than 10 Gyr. In the case with AGN feedback the galaxy has too few stars belonging to the first group, so we skipped it. Without AGN feedback there are instead many stars in this group, and they are almost all on rotational orbits in a thin disc, with few stars above or below its plane. This disc rotates very fast, at around $400 \text{ km s}^{-1}$, has very little velocity dispersion and shows quite extreme signatures in the $h_3$ and $h_4$ maps. Interestingly the intermediate age group shows an extended velocity dispersion signature, perhaps due to the presence of several non-aligned remnants of discs that were also quenched in the case with AGN feedback. The older stars behave in a relatively similar way with or without AGN feedback. Their rotational velocity is smaller and their dispersion much higher, similar to the maps for the accreted component (Top panel of Fig. 11). From the youngest to the oldest age group, the value of the $\lambda_R$ parameter is 0.84, 0.38, and 0.35 for the case without AGN, and 0.10 and 0.13 for the case with AGN. The old stars rotate faster in the case without AGN, likely because of the deep central potential.

<figure><img src="image_24.png" /><figcaption>**Figure A1.** Kinematic maps ($V_{\text{avg}}$, $\sigma$, $h_3$, $h_4$) for age-selected stellar components of galaxy 0227 simulated without AGN. From top to bottom, stars between 0 and 3 Gyr old, between 3 and 10 Gyr old, and older than 10 Gyr.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure A2.** Kinematic maps ($V_{\text{avg}}$, $\sigma$, $h_3$, $h_4$) for age-selected stellar components of galaxy 0227 simulated with AGN. From top to bottom, stars between 0 and 10 Gyr old and stars older than 10 Gyr.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.