Hydrodynamical simulations of early-type galaxies - AMS Dottorato

Hydrodynamical simulations of early-type galaxies - AMS Dottorato

Alma Mater Studiorum – Università di Bologna DIPARTIMENTO DI FISICA E ASTRONOMIA Corso di Dottorato di Ricerca in Astronomia CICLO XXVII Hydrodynamic...

16MB Sizes 0 Downloads 5 Views

Recommend Documents

The Evolution of Central Group Galaxies in Hydrodynamical Simulations
Nov 13, 2009 - simulated groups this leads to a dramatic rejuvenation of the central group galaxy at z < 1, affecting it

tesi dottorato defitivita - AMS Dottorato
GABRIELE BALDUCCI. Prof. FABIO PEZZI ..... del prodotto (dato dal fatto della grande quantità di prodotto ammostato) e

Schema frontespizio teso Dottorato - AMS Dottorato
List of Figures. 6 .... African Union‟s Inter-african Bureau of Animal Resources ..... forced by destitution to depend

chapter 1 - AMS Dottorato
satellite images for flood mapping and monitoring, damage assessment and risk assessment. The contribution of satellite

Tesi Dottorato Stanzani-Maserati + FRONTESPIZIO - AMS Dottorato
MICHELANGELO STANZANI MASERATI. Coordinatore Dottorato. Relatore. PROF. PIETRO CORTELLI. PROF. PIETRO CORTELLI. Esame fi

Schema frontespizio teso Dottorato - AMS Dottorato
Emma, Middlemarch, David Copperfield. 3.0. Introduzione p. 137. 3.1. Le traduzioni italiane di Emma (1816) p. 143. 3.1.1

Facsimile del frontespizio della tesi di dottorato - AMS Dottorato
Mar 11, 2015 - of oil is equivalent to ten minutes of operation at low power or three ...... sidering the average per ca

Low-temperature thermochronological evolution of - AMS Dottorato
Low-temperature thermochronological evolution of the Menderes and Alanya massifs. (Turkey). Presentata da: Dott. Frances

engineering agent-oriented technologies and - AMS Dottorato
II Engineering Agent-Oriented Technologies for Programming Multi-. Agent Systems .... 7.7.1 Comparison with State-of-the

Facsimile del frontespizio della tesi di dottorato - AMS Dottorato
List of ..... of Molecular Biology, 1974-1980” Isis 92.3 (2001): 541-75; Arnold Thackeray (ed.) ... the grocery store,

Alma Mater Studiorum – Università di Bologna DIPARTIMENTO DI FISICA E ASTRONOMIA Corso di Dottorato di Ricerca in Astronomia CICLO XXVII

Hydrodynamical simulations of early-type galaxies: effects of galaxy shape and stellar dynamics on hot coronae

Dottorando:

Andrea Negri

Relatore:

Chiar.mo Prof. Luca Ciotti Correlatore:

Chiar.ma Prof.ssa Silvia Pellegrini Coordinatore:

Chiar.mo Prof. Lauro Moscardini

Settore Concorsuale: 02/C1 – Astronomia, Astrofisica, Fisica della Terra e dei Pianeti Settore Scientifico-Disciplinare: FIS/05 – Astronomia e Astrofisica Esame finale anno 2014

The scientist does not study Nature because it is useful; he studies it because he delights in it, and he delights in it because it is beautiful. J. H. Poincaré

Contents List of acronyms

xi

1 Introduction

1

2 The simulations 2.1 The galaxy models . . . . . . . . . . . . . . . . . . . . . . . 2.2 The input physics . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Passive stellar evolution: stellar winds and SNe Ia . 2.2.2 Star formation and SNe II heating . . . . . . . . . . 2.2.3 Radiative cooling . . . . . . . . . . . . . . . . . . . . 2.3 The contribution of stellar kinematics to the ISM energetics 2.4 The code . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 The hydrodynamical equations . . . . . . . . . . . . 2.4.2 Numerical methods . . . . . . . . . . . . . . . . . . . 2.4.3 Code setup and outputs . . . . . . . . . . . . . . . .

. . . . . . . . . .

3 The effects of stellar dynamics on the X-ray emission of S0 3.1 Setting the galaxy models . . . . . . . . . . . . . . . . . . . 3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . 3.2.2 The thermalization parameter . . . . . . . . . . . . . 3.2.3 The X-ray observables: LX , TX , and ΣX . . . . . . . 3.2.4 Comparison with observed X-ray properties . . . . . 3.3 Discussion and conclusions . . . . . . . . . . . . . . . . . . .

galaxies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

4 The effects of galaxy shape and rotation on the X-ray haloes ETGs - Numerical simulations 4.1 The setting up of representative galaxy models . . . . . . . . . . . 4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 The X-ray ISM luminosity LX . . . . . . . . . . . . . . . . . 4.2.3 The X-ray emission weighted temperature TX . . . . . . . . 4.3 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

7 8 10 10 12 13 14 15 16 16 19 23 24 26 26 41 43 45 46

of . . . . . .

53 54 55 55 65 68 70

5 Ongoing and future developments 5.1 Star formation in ETGs . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Star formation input physics and galaxy models selection . .

79 80 80

v

vi

CONTENTS

5.2

5.1.2 5.1.3 AGN 5.2.1 5.2.2 5.2.3

Results . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . feedback in ETGs . . . . . . . . . Code implementation and testing Preliminary results . . . . . . . . Gravitational viscosity . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

80 92 100 102 104 104

6 Conclusions

109

A Basic concepts of fluid mechanics A.1 From the first principles of conservation to the basic equations of fluid mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . A.1.2 Linear momentum conservation . . . . . . . . . . . . . . . . . A.1.3 Energy conservation . . . . . . . . . . . . . . . . . . . . . . . A.1.4 Angular momentum conservation . . . . . . . . . . . . . . . . A.1.5 Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Source and sink terms inside an elliptical galaxy . . . . . . . . . . . . A.2.1 From microscopic to macroscopic source terms . . . . . . . .

113 114 114 115 117 118 120 121 122

B Vector operators in cylindrical and spherical coordinates 125 B.1 The Cartesian, cylindrical and spherical coordinate systems . . . . . 126 B.1.1 Coordinates and vector transformations . . . . . . . . . . . . 127 B.2 Differential operators in cylindrical and spherical coordinates . . . . 127 C Viscosity module implementation 129 C.1 Discretized viscosity equations . . . . . . . . . . . . . . . . . . . . . . 130

List of Figures 1.1

LX versus LK of the ATLAS3D sample of Sarzi et al. (2013). . . . . .

2.1 2.2

The bolometric cooling function Λ(T ) . . . . . . . . . . . . . . . Spatial position of an arbitrary point on the computational domain the ZEUS-MP 2 code . . . . . . . . . . . . . . . . . . . . . . . . . Numerical grid resolution in the radial direction. . . . . . . . . . The 0.3–8 keV specific emissivity ΛX (T ) . . . . . . . . . . . . . .

2.3 2.4 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 4.1 4.2

4.3 4.4 4.5

. . in . . . . . .

Meridional sections of vϕ and of σϕ (bottom) for the ISi and VDi models Meridional sections of vϕ and of σϕ (bottom) for the CRi and RDi models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hydrodynamical evolution of the VDi model . . . . . . . . . . . . . . ISM temperature evolution for model VDi . . . . . . . . . . . . . . . Hydrodynamical evolution of the ISi . . . . . . . . . . . . . . . . . . ISM temperature evolution of the ISi model . . . . . . . . . . . . . . Hydrodynamical evolution of the counter-rotating CRi model . . . . Temperature evolution of the counter-rotating CRi model . . . . . . Hydrodynamical evolution of the velocity dispersion supported model with rotating inner disc . . . . . . . . . . . . . . . . . . . . . . . . . Temperature evolution of the RDi model . . . . . . . . . . . . . . . . Time evolution of the thermalization parameter γth . . . . . . . . . . Time evolution of LX and X-ray emission weighted temperature TX Edge-on 0.3–8 keV surface brightness of the ISM at 13 Gyr . . . . . Meridional section of temperature and heating over cooling time ratio theat /tcool for the E0250 model of the NFW set . . . . . . . . . . . . . Time evolution of the X-ray luminosity LX and X-ray emission weighted temperature TX for the family derived from the E0250 model with the NFW halo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Meridional sections of the temperature and heating over cooling time ratio for the EO7250 VD model of the NFW set . . . . . . . . . . . . . . Meridional sections of the temperature and heating over cooling time ratio for the EO7250 IS model of the NFW set . . . . . . . . . . . . . . Meridional sections of the temperature and heating over cooling time ratio for the low mass EO7200 VD model of the NFW set . . . . . . . . . vii

2 14 17 18 21 27 28 31 32 34 36 38 39 40 42 43 47 48 57

58 59 61 63

viii

LIST OF FIGURES

4.6 4.7 4.8 4.9

4.10 4.11 5.1

ISM X-ray luminosity LX in the 0.3–8 keV band at 13 Gyr for all models in the NFW and in the Einasto sets, as a function of σe8 , of the hot ISM mass, and of the galaxy blue optical luminosity . . . . . Angle-averaged profile of the hot ISM density at t = 13 Gyr for the same models as in Fig. 4.2 . . . . . . . . . . . . . . . . . . . . . . . . Fraction of escaped ISM mass (Mesc ) with respect to the total injected mass at t = 13 Gyr, as a function of σe8 , for the whole NFW and Einasto sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ISM emission weighted temperature TX in the 0.3–8 keV band at 13 Gyr for all the models in the NFW and in the Einasto sets as a function of σe8 , of the hot ISM mass, and of the galaxy blue optical luminosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . X-ray luminosity LX with respect to X-ray luminosity weighted temperature TX at t = 13 Gyr for the whole NFW and Einasto sets . . . Edge-on 0.3–8 keV surface brightness of the ISM (ΣX ) at 13 Gyr, for 250 E0250 , EO7250 VD and EO7IS models . . . . . . . . . . . . . . . . . . .

Meridional sections of the heating over cooling time ratio for the EO7200 IS models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Time evolution of the X-ray luminosity LX and X-ray emission weighted 200 300 temperature TX for the EO7200 IS , FO7IS , FO7IS models . . . . . . . 5.3 Meridional sections of the heating over cooling time ratio for the FO7200 IS models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Meridional sections of the heating over cooling time ratio for the EO7300 IS models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Meridional sections of the time integrated density distribution of the 200 300 new stars formed up to 13 Gyr, for the EO7200 IS , FO7IS , and EO7IS models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Time evolution of the SFR and hot mass contained in the numerical 200 300 grid for the EO7200 IS , FO7IS , FO7IS models . . . . . . . . . . . . . . 5.7 ISM X-ray luminosity LX in the 0.3–8 keV band at 13 Gyr for the selected sub-group of the NFW set, as a function of σe8 . . . . . . . 5.8 Toomre Q parameter and surface density Σ for the cold gaseous disc of EO7300 IS at 13 Gyr. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 SFR, BH mass and accretion rate as a function of time from the begin of the simulation, for few representative tests . . . . . . . . . . . . . 5.10 Black hole bolometric luminosity and Eddington ratio evolution with time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Time evolution of the gaseous and stellar mass for the E0250 , EO4250 VD and EO4250 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IS 5.12 AGN duty cycle as a function of Eddigton ratio for the models of Figure 5.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64 66 67

67 68 70 82 83 85 88 89 90 91 101 105 106 107 108

List of Tables 2.1 2.2

ZEUS numerical grid formalism . . . . . . . . . . . . . . . . . . . . . X-ray emissivity coefficients . . . . . . . . . . . . . . . . . . . . . . .

16 20

3.1 3.2

Stellar kinematics of the models . . . . . . . . . . . . . . . . . . . . . Main outputs at the end of the simulations . . . . . . . . . . . . . .

25 29

4.1 4.2 4.3

Fundamental galaxy parameters for the NFW and Einasto sets of models 56 Simulations results for the NFW set at t = 13 Gyr . . . . . . . . . . 74 Simulations results for the Einasto set at t = 13 Gyr . . . . . . . . . 76

5.1

Simulations results for the σe8 = 300 km s−1 IS family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 94 Simulations results for the σe8 = 300 km s−1 IS family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 95 Simulations results for the σe8 = 200 km s−1 IS family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 96 Simulations results for the σe8 = 200 km s−1 IS family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 97 Simulations results for the σe8 = 200 km s−1 VD family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 98 Simulations results for the σe8 = 200 km s−1 VD family of the NFW set with star formation at t = 13 Gyr. . . . . . . . . . . . . . . . . . 99 Test simulations of the mass, momentum and energy injection terms. 103

5.2 5.3 5.4 5.5 5.6 5.7

ix

List of acronyms AGN Active Galactic Nuclei CE Cooling episode CR Counter-rotating DM Dark matter EO Edge-on ETGs Early-type galaxies FO Face-on FTCS Forward time-centred space IS Isotropic rotator ISM Interstellar medium RD Rotating disc SF Star formation SFR Star formation rate SMBH Super massive black hole SNIa Type Ia supernovae SNII Type II supernovae VD Velocity-dispersion

xi

Chapter 1

Introduction Early-Type galaxies (ETGs) are embedded in hot (106 − 107 K), X-ray emitting gaseous haloes (Fabbiano 1989; O’Sullivan et al. 2001), produced mainly by stellar winds and heated by Type Ia supernovae (SNIa) explosions, by the thermalization of stellar motions and occasionally by the central super-massive black hole (SMBH). In particular, the thermalization of the stellar motions is due to the interaction between the stellar and the SNIa ejecta and the hot interstellar medium (ISM) already residing in the ETG (e.g. Parriott & Bregman 2008a). A number of different astrophysical phenomena determine the X-ray properties of the hot ISM, such as stellar population formation and evolution, galaxy structure and internal kinematics, Active Galactic Nuclei (AGN) presence, and environmental effects. With the aid of high-resolution hydrodynamical simulations, in this Thesis we focus on the effects of galaxy shape, stellar kinematics and star formation on the evolution of the X-ray coronae of ETGs. To put this Thesis in perspective, one of the empirical discoveries that followed the analysis of first X-ray data of ETGs was the sensitivity of the hot gas content to major galaxy properties as the shape of the mass distribution, and the mean rotation velocity of the stellar component (see Mathews & Brighenti 2003; Kim & Pellegrini 2012 for a full discussion of the most relevant observational and theoretical aspects concerning the X-ray haloes). From Einstein observations it was found that, on average, flatter ETGs and S0 galaxies have lower total X-ray luminosity LX and LX /LB (a measure of the galactic hot gas content), than rounder systems of similar optical luminosity LB (Eskridge et al. 1995). Moreover, galaxies with axial ratio close to unity span the full range of LX , while flat systems have LX . 1041 erg s−1 . This result holds across the whole range of LB where the two shapes, round and flat, coexist (Pellegrini 1999a), as shown in Figure 1.1. In addition, it was found that LX /LB can be high only in slowly rotating galaxies, and is limited to low values for fast rotating ones (Pellegrini et al. 1997; Sarzi et al. 2010). The relationship between LX and shape/rotation was reconsidered, confirming the above trends, for the ROSAT PSPC sample (Pellegrini 2012), and for the Chandra sample (Li et al. 2011a; Boroson et al. 2011; Sarzi et al. 2013). In particular, Sarzi et al. (2013), after confirming that slow rotators generally have the largest LX and LX /LB values, also found that their gas temperature values TX are consistent just with the thermalization of the stellar kinetic energy, estimated from σe (the luminosity averaged stellar velocity dispersion within the effective radius Re ). Fast rotators, instead, have generally lower LX and LX /LB values, and the more so the larger their degree of rotational

The X-ray Halos of F 2

Introduction

Figure 1.1: LX versus LK of the ATLAS3D sample of Sarzi et al. (2013). The Figure 5. symbols Same as Figs. and 3, but withand elliptical symbolsaccording coded by colour anddegree flattening elliptical are 2coded by now colour flattening to the ofaccording to the de ellipticity of the galaxies they are representing, in order to visually convey the dependence of the scatter in their tot rotational support and the apparent ellipticity of the galaxies they are representing. parameters, as quantified in Fig. 4. The grey labels indicate the galaxies that were excluded from Low and high values of λRe correspond to slow and fast rotators (Emsellem et al. the analysis of the galaxies labelled left and right are, contribution respectively, most likely face-on Lbut values, intrinsically flat and rotati 2007). The solid on linetherepresents the panels expected to the observed X slow rotators in our low X-ray resolution ATLAS3D subsample. Finally, in both panels we also label in grey the tw with uncertainties, from the unresolved emission of low mass X-ray binaries, as given NGC 4649 and NGC 5353. The kinematic classification of NGC 4649 is somehow uncertain, however, since this obj in Boroson et al. (2011).

between fast and slow rotators defined in Paper III, whereas the very flat NGC 5353 happens to be the most massive HCG 68, so that its total LX may also include emission from hot gas confined by the group potential.

support; the TX values of fast rotators keep below 0.4 keV and do not scale with σe . Summarizing, there seems be a dependence of the involving hot gas content 3.2 High X-ray Resolution from the Spearman rank to analysis the correlations λRe and temperature on the galactic shape and internal dynamics. The investigation of the origin of this in Fig. 4 become more significant, and the same holds when exThe previous conclusions o dependence is the one of the main goal of the present Thesis. cluding the face-on fast rotators from the correlations involving �e . rotating, or even more genera From we an observational of view, the determination of the origin of these Finally, note that evenpoint a modest contamination from a centain their halos of hot gas hav trends is notortrivial. fact, since or rotation is medium dependent only galaxies tral AGN the theInintra-group -cluster to on theflattening, total are that basedthe on X-ray data of r sufficiently flattened are expected rotate significantly; moreover, given X-ray in the least massive objectstocould significantly impact on the contribution from a cent observed in real is biased projection effects, intrinsically flattened for objects these object. In by turn, this would generthe totalflattening LX values medium or LMXBs to the to galaxies, if observed face-on, can be found in thebetween region occupied rounder galaxies ally weaken the significance of any correlation λRe or �by e andobserved which in the case of LM inand thethe ellipticity-X-ray properties diagrams, increasing the scatter in the LX /LX,discr or LX /LX,discr+diff ratios, since low-mass 2 trends. Thus, it is difficult disentangle rotationalconsiderable and purely scatter. It is th galaxies generally tend to to beobservationally fast rotators (Paper III). Thepurely possible results hold when such factor flattening effects. impact of AGNs or the intra-group or -cluster medium is most nouse of Chandra or XMM. Theoretically, differentbetween explanations have proposed, and explored both the been LX /L ticeable in the correlation λRe and X,discr+diff Recently, Boroson, Kim analytically (Ciotti & Pellegrini hereafter CP96; Posacki et al. 2013b, hereafter deficit (lower left panel of Fig. 1996, 4), where the majority of the most the largest samples of early-ty P13) and numerically (Brighenti & Mathews 1996; D’Ercole Ciotti 1998, hereafter values&near rotationally-supported galaxies with LX /LX,discr+diff measurements, whic LX,gas 2 2 −2 DC98). Thelog proposed be classified categories: the enervalues below can of 15.4 LK,� kmin stwo.broad At such unity have LK σe explanations sample and 2 hydrodynamical ones. Of course, the difference is not sharp, as theexcludes the dom getic ones the a low logand LK σ e regime the typical AGN X-ray emission LX,AGN of galaxies. For the objects th flow energeticsgalaxies affects the hydrodynamical and different of early-type could be comparable evolution, to the predicted values hydrodynamical ple of Boroson et al., Fig. 6 sh for LX,discr and LX,diff (Fig. 3). For instance, the median LX,AGN ues compare to each other, th value for the ATLAS3D sample galaxies that appear in the study of the hot-gas luminosity corres Liu (2011, 42 objects) is 2.6 × 1039 erg s−1 . ergy from stellar-mass loss. T To summarize the results of this subsection, in Fig. 5 we redraw the same LK − LX and LK σe2 − LX diagrams presented in Figs. 2 and 3, but now with elliptical symbols coded by colour and 2 If we consider that globular flattening according to the values of λRe and the �e of the galaxies

3

configurations lead to different redistributions of the energy available. Explanations based on energetic effects suppose that the ISM in flat and rotating ETGs is less bound than in rounder and non rotating galaxies of similar luminosity, so that in the former objects the ISM is more prone to develop a global/partial galactic wind, with the consequent decrease of LX . A subdivision of energetic explanations is given by flattening effects versus rotational effects: the former supposes that flat galaxies have shallower potential wells than rounder galaxies of similar luminosity, so the gas is less bound (independently of the galaxy kinematical support, ordered rotation or velocity dispersion); the latter assumes that the gas injected in the galaxy retains the stellar streaming motion, and so it is less bound in rotating galaxies than in velocity dispersion supported galaxies of similar shape. Note that, at variance with the stellar random kinetic energy Lσ which is always supplied to the ISM (equation 2.29), the thermalization of the stellar ordered kinetic energy Lv depends on the relative motion between the stellar population and the ISM (equation 2.30), therefore it cannot be predicted a priori from the galaxy structure and kinematics. The energetic scenario has been explored analytically in CP96 by using two-component Miyamoto-Nagai models, and in P13 by building a wide set of more realistic ETG models, with different flattened structures and kinematics. These works showed that the binding energy of the gas depends on the procedure adopted to “flatten” the galaxy models, and that rotation affects also the hot gas temperature. Explanations based on hydrodynamical effects are less direct. In this scenario the rotation of the gas injected in the galaxy leads to hydrodynamical configurations of the ISM that, for different but cooperating reasons associated with angular momentum conservation, are less X-ray luminous than in non rotating systems. For example, the gas density in the central galactic region is lower in rotating models, where a rotationally supported cold disc forms, than in non-rotating ones, where the gas flows straight to the centre, forming a hot, dense core. Thus, the X-ray faintness is not due to the onset of galactic winds, but to redistribution of the gas inside the galaxy. The rotation is the main driver of X-ray under luminosity, and correlation with galaxy flattening is a by-product. The dependence of the energetics on the kinematical status of the ISM is the reason why hydrodynamical simulations are required; the implementation of such simulations and their results are the main subject of this Thesis. Numerical simulations have been used in this Thesis also to explore the effects of counter-rotating structures in ETGs. Counter-rotating discs have been observed in the central parts of ETGs, e.g., in the cases of NGC7097 (De Bruyne et al. 2001), NGC4478 and NGC4458 (Morelli et al. 2004), NGC3593 and NGC4550 (Coccato et al. 2013), IC719 (Katkov et al. 2013), NGC 4473 (Foster et al. 2013) with a counterrotating disc summing up to the 30% of the total stellar mass, and other cases in Kuijken et al. (1996) and in Erwin & Sparke (2002). Moreover, Krajnović et al. (2011) identified 19 galaxies hosting a counter rotating stellar component over the 260 ETGs of the ATLAS3D sample (Cappellari et al. 2011). Two competing effects could be at work as the ISM flows towards the centre of a galaxy with a counter-rotating stellar structure. On one hand, as a consequence of cooling and conservation of angular momentum, the infalling gas increases its rotational velocity until it reaches the region where the counter-rotating disc lies, and then interacts with a counter-rotating mass injection, with the consequent additional heating contribution in Lv (equation 2.30). On the other hand, this interaction also causes a reduction of the specific angular

4

Introduction

momentum of the local ISM, which will reduce the local centrifugal support, and will favour central accretion. What of the two competing effects will dominate can be quantified only with high resolution numerical simulations. An additional, potentially relevant effect related to ordered rotation is the possibility of large-scale instabilities in the rotating ISM, as those revealed by the simulations of hot gas flows in S0 galaxies presented in DC98, where the formation of large scale cold filaments are found in rotating flattened galaxies. Renewed interest in these hot halo instabilities came after the observation in fast rotators of cold, molecular gas discs, almost undetected in slow rotators (Davis et al. 2011; Young et al. 2011; Sarzi et al. 2013; Davis et al. 2013). The observed molecular discs show a kinematical misalignment that vanishes in the case of massive fast rotators, thus suggesting an internal origin of the cold gas. Davis et al. (2011) proposed that the hot halo, formed by the secular mass losses of the ageing stellar population, may be unstable and then hosting cooling episodes that produce the detected cold gas. However, the physical explanation of the lack of molecular gas in slow rotators is still debated. In this context, the simulations have explored the formation of cold gaseous discs and the related star formation processes in both rotating and non-rotating galaxies. Summarizing, the main goal of the present Thesis is to investigate the dependence of the hot gas content in ETGs on the galaxy shape and stellar dynamics by using the ZEUS-MP 2 code to perform high-resolution, 2D hydrodynamical simulations of gas flows in realistic, state-of-the-art models of ETGs. In particular, we aim to (1) disentangle the effects of flattening and rotation in determining the X-ray properties of flat and rotating ETGs, (2) investigate both the energetic and the hydrodynamical explanations, so as to quantify their relative importance and their role in determining the X-ray under-luminosity of flat and rotating ETGs, (3) understand the effects of a counter-rotating stellar disc on the hot gas dynamical status, (4) probe the hot halo instabilities in both rotating and non-rotating galaxies, estimate the consequent star formation history and investigate the star formation effects on the halo X-ray properties. Given the complexity of the subject, we performed three different, but correlated, studies. In a first, exploratory work, we followed the evolution of the hot ISM in flat ETGs, modelled as two-component, self-consistent, axisymmetric realistic S0 galaxies, with different degrees of rotational support, dark matter and counterrotation amounts. The total gravitational field, the ordered and random stellar motions are computed by solving the Poisson and the Jeans equations with the code described in P13. The hydrodynamical simulations account for a passively evolving stellar population with the associated injection of mass, momentum and energy from Type Ia supernovae and stellar winds of giant stars, along with the thermalization of stellar streaming and random motions. In order to reduce the number of free parameters, we focussed only on the effects of rotation and counter-rotation on the hot halo by keeping constant the shape of the model galaxy, leaving the study of galaxy flattening dependence for the subsequent work. Motivated by the simulations results of the exploratory work, we expanded our research by performing hydrodynamical simulations of a large set of realistic, state-ofthe-art galaxy models constructed with the P13 code and tailored to reproduce the observed properties and scaling laws of ETGs. The galaxy models are characterised by different stellar mass, intrinsic flattening, distribution of dark matter, and internal

5

kinematics. In the two previous works we found an ubiquitous presence of cooling episodes and formation of cold discs in rotating systems, due to angular momentum conservation; these are absent in velocity dispersion supported galaxies. For this reason, we investigated the effects of the star formation in ETGs on hot flows in a third work, where the numerical code accounts for star formation and Type II SNe explosions, as well as the injection of energy due to the thermalization of the random and ordered motions of the new stellar population. We re-simulated a sub-group of the previous entire large set of galaxies with different star formation schemes and efficiencies, in order to provide an adequate covering of the parameter space. Finally, we are now moving to consider also the effects of AGN feedback in a joint work in collaboration with G.S. Novak and J.P. Ostriker. The Thesis is organized as follows. In Chapter 2 we describe the galaxy models, the physical processes considered in the simulations and the adopted numerical code. In Chapter 3 we present our exploratory work on the X-ray emission of S0 galaxies, while Chapter 4 is dedicated to the exploration of hot flows in a large set of ETGs with different shape and internal kinematics. The studies regarding the star formation, its effects on the X-ray coronae, and the ongoing work related with AGN feedback simulations are presented in Chapter 5. In Chapter 6 we draw our conclusions.

Chapter 2

The simulations In the last two decades, numerical simulations became an essential tool in astrophysics, allowing the study of an astonishing variety of physical phenomena on a different scales, ranging from the large scale structure formation, to cluster and galaxy evolution, down to accretion on BHs, even in presence of a complex physics (e.g. cooling, star formation, AGN feedback, accretion disc, nuclear reactions, etc). The main focus of this Thesis is simulate the hot, X-ray emitting halos of ETGs. Given the great variety of physical phenomena present in the ISM evolution (cooling, formation of cold discs and star formation, to mention a few), this study would have been impossible without numerical hydrodynamical simulations. We adopted the ZEUS-MP 2 code (Hayes et al. 2006), widely used in the astrophysical community, to perform 2D hydrodynamical simulations of ETGs X-ray haloes. Since ZEUS-MP 2 is a general purpose fluid dynamical code, the computational core of the original (public) version includes only the standard equations of fluid dynamics (i.e. without source or sink terms). Moreover, it lacks of all the astrophysical processes typical of a galactic environment, such as the baryon cooling or the evolution of the stellar population. Finally, the initial conditions (stellar model, dark matter halo and central black hole) must be set by the user. We modified the code to simulate realistic, axisymmetric, two-component stateof-the-art galactic models having self-consistent rotation, random motions and selfgravity, calculated with a dedicated Jeans code (P13; Posacki 2014). In addition, we incorporated new physical processes, such as the ISM radiative cooling, the self-consistent injection of mass, momentum and energy by stellar winds and type Ia SNe, star formation, type II SNe feedback, and thermalization processes of the mass ejected by stars when it hits the ISM in situ. The intent of this Chapter is to give a complete description of the new physical modules implemented in the code. The Chapter is organised as follows. In Section 2.1 we describe the structural and dynamical properties of the galaxy models. In Section 2.2 we present the detailed input physics of the numerical code, while Section 2.4 the code numerical methods are reported.

8

2.1

The simulations

The galaxy models

For the simulations we adopted axisymmetric two-components state-of-the-art galaxy models, produced by a dedicated code built on purpose (P13; Posacki 2014). Here we summarize the building process and the main features of the galaxy models. Each galaxy model consist of a stellar distribution, embedded in a gravitational field due to the different contributions of the stellar population self-gravity, a dark matter halo (DM) and a central BH. The code can produce various type of stellar profiles having different intrinsic shape and flattening, among which we employed the de-projection (Mellier & Mathez 1987) of the de Vaucouleurs (1948) law, generalized for ellipsoidal axisymmetric distributions, ρ∗ (R, z) = ρ0 ξ −0.855 e−ξ

ρ0 =

b12

M∗ , 16πqRe3 0 Γ(8.58)

ξ=

1/4

b4 Re 0

(2.1)

,

s R2 +

z2 , q2

(2.2)

where (R, ϕ, z) are the cylindrical coordinates and b ' 7.67. The flattening is controlled by the parameter q 6 1, so that the minor axis is aligned with the z axis. Re 0 is the projected half mass radius (effective radius) when the galaxy is seen √ face-on; for an edge-on view, the circularized effective radius is Re = Re 0 q. While we applied the de Vaucouleurs (1948) profile to describe the stellar surface brightness of elliptical galaxies in Negri et al. (2014b), we employed the Miyamoto & Nagai (1975) density profile ρ∗ (R, z) =

M∗ b2 aR2 + (a + 3ζ)(a + ζ)2 , 4π ζ 3 [R2 + (a + ζ)2 ]5/2

(2.3)

to √ simulate S0 galaxies (Negri et al. 2014a), where a and b are scalelenghts, ζ ≡ z 2 + b2 . For a = 0, equation (2.3) reduces to the Plummer (1911) model, while for b = 0 to the Kuzmin (1956) disc. As a DM halo, we adopted the spherical untruncated Navarro-Frenk-White profile (Navarro et al. 1997, hereafter NFW): ρh (r) =

ρcrit δc rh , r(1 + r/rh )2

Φ(r) = −4πGρcrit δc rh3

(2.4)

ln(1 + r/rh ) , r

(2.5)

√ where r = R2 + z 2 and ρcrit = 3H 2 /8πG are the spherical radius and the critical density for closure, respectively. The total mass of the NFW profile diverges, so we refer to the DM mass enclosed within r200 (the radius of a sphere of mean interior density 200ρcrit ), as to the halo mass Mh . The concentration parameter c ≡ r200 /rh and the δc coefficient are related as δc =

200 c3 , 3 ln(1 + c) − c/(1 + c)

c≡

r200 . rh

(2.6)

2.1 The galaxy models

9

Recently, the Einasto (1965) profile has been reconsidered as a suitable choice to model DM haloes (Navarro et al. 2004; Merritt et al. 2006; Gao et al. 2008; Navarro et al. 2010). The spherical Einasto density distribution is the three dimensional equivalent of the Sérsic law:  1/n r ρh (r) = ρc edn −x , x ≡ dn , (2.7) rh   GMh Γ(3n, x) xn Γ(2n, x) Φh (r) = − 1− + , r Γ(3n) Γ(3n)

(2.8)

where rh is the half mass radius, n is a free parameter, and ρc = ρh (rh ) = −dn /(4πnΓ(3n)r 3 ) is the density at the spatial half-mass radius. For d Mh d3n n n e h we use the asymptotic relation dn ' 3n −

1 8 + 3 1215n

(2.9)

(Retana-Montenegro et al. 2012). The stellar and DM halo profiles are discretized on a high-resolution numerical grid in cylindrical coordinates, then the stellar gravitational potential, the radial and vertical forces are calculated by solving the Poisson equation. Afterwards, the same code solves the Jeans equations ∂ρ∗ σ 2 ∂Φtot = −ρ∗ , ∂z ∂z

(2.10)

σ 2 − vϕ2 ∂ρ∗ σ 2 ∂Φtot + ρ∗ = −ρ∗ , ∂R R ∂R

(2.11)

where Φtot is the total gravitational potential (due to stars, dark matter and BH), and σ is the stellar velocity dispersion (see below). The stellar random and ordered motions are therefore calculated. As well known (Ciotti & Pellegrini 1996; Ciotti 2000), for a two-integral distribution function the radial and vertical velocity dispersion are equal, σR = σz ≡ σ, and the only non-zero streaming motion is in the azimuthal direction. However, the Jeans equations provides only vϕ2 , so that we adopted the Satoh (1980) k-decomposition to calculate the stellar streaming velocity   vϕ2 = k 2 vϕ2 − σ 2 , (2.12) and the azimuthal velocity dispersion σϕ = σ + (1 + k 2 )(vϕ2 − σ 2 ),

(2.13)

where 0 ≤ k ≤ 1. In the case k = 0, no streaming motions are present and all the flattening is due to the azimuthal velocity dispersion, whereas the case k = 1 represents the isotropic rotator. In principle, the k parameter is a function of the position on the (R, z) plane, and so more complicated rotational velocity fields can be realized (e.g. see Chapter 3 and Negri et al. 2014a). Finally, all the galaxy model physical fields are then interpolated from the original high-resolution cylindrical grid on the hydrodynamical mesh by a separate code,

10

The simulations

via bilinear interpolation (Press et al. 1992). In addition, in the case of a hydro simulation in spherical coordinates, the code transforms all the vector quantities1 (stellar and DM gravitational fields) from cylindrical to spherical components by means of equation B.8.

2.2

The input physics

In this section we present the input physics implemented in the ZEUS code, while the detailed numerical schemes are presented in Section 2.4. As usual in numerical works of hot flows in ETGs (D’Ercole & Ciotti 1998; Ciotti & Pellegrini 1996; Negri et al. 2014a,b), in order to simulate the evolution of ETGs over cosmic time-scales, we account for: (1) passive stellar evolution, (2) star formation and Type II supernovae explosions, (3) radiative cooling. Moreover, we included the treatment of the ISM viscosity.

2.2.1

Passive stellar evolution: stellar winds and SNe Ia

The principal ingredient of the ETGs evolution is the initial stellar population, which evolves in a passive fashion. Two main contributions due to the primordial stellar population can be identified: gas injections from post main-sequence stars via stellar winds, which is the dominant source of mass in isolated ETGs, and Type Ia SNe explosions. Both the physical processes inject mass (due to stellar winds and SN ejecta mass), momentum (due to ordered streaming motions of the parent stars) and energy (by SNIa explosions and thermalization of random and streaming stellar motions) in the ISM. The rigorous derivation of the hydrodynamical equations with multiple isotropic source fields is given in Section A.2; by characterizing the mass and specific energy injection rates of the two physical processes we can explicitly write the source terms present in equations (A.48)-(A.50). In addition, the source terms are further simplified by the fact that the two sources belong to the same stellar population, thus they share the same kinematical properties. The mass injection rate due to stellar winds and SNIa explosions are parametrized respectively by ρ˙ SN = αSN (t)ρ∗ and ρ˙ ∗ = α∗ (t)ρ∗ (Pellegrini 2012), where αSN (t) =

1.4M RSN (t), M∗

−1.3 α∗ (t) = 3.3 × 10−12 t12 (yr−1 ).

(2.14) (2.15)

For the SNIa explosion rate RSN (t), we adopted the parametrization 2 −1 RSN (t) = 0.16 H0,70 × 10−12 LB t−s 12 (yr ),

(2.16)

where H0,70 is the Hubble constant in units of 70 km s−1 Mpc−1 , LB is the present epoch B-band galaxy luminosity in blue solar luminosities, t12 is the age of the stellar 1

In the case of axisymmetric galaxy models, it can be proved that the following relation holds: σR = σz = σr = σθ ≡ σ, therefore σ on the spherical mesh is simply calculated via bilinear interpolation, as all scalar fields, e.g. the stellar density.

2.2 The input physics

11

population in units of 12 Gyr, and s parametrizes the past evolution. Equation (2.15) holds for a Kroupa Initial Mass Function (Pellegrini 2012). The SNIa’s heating rate is obtained as LSN (t) = ESN RSN (t) erg yr−1 , where ESN = 1051 erg is the kinetic energy of one event. Following recent theoretical and observational estimates of the SNIa explosion rate (Mannucci et al. 2005; Greggio 2005, 2010; Sharon et al. 2010; Maoz et al. 2011), we adopted s = 1. The momentum injection term accounting for ˆϕ − u), both the sources follows directly from equation (A.49) and it is equal to ρ(v ˙ ϕe where we defined ρ˙ ≡ ρ˙ SN + ρ˙ ∗ . At variance with the mass and momentum explicit source terms, which have been calculated in an exact way, the energy injection requires an in-depth analysis. From equation (A.50), the exact source term for a generic isotropic source field is given by:   i u2s ρ˙ h 2 2 ˙ , kv − uk + Tr(σ ) + ρ˙ einj + E= 2 2

(2.17)

where ρ, ˙ v, u, einj , us , σ 2 are the mass injection rate per unit volume, the source streaming velocity field, the velocity of the ambient gas, the internal energy per unit mass of the injected gas, the modulus of the relative velocity of the injected material and the source (i.e., the velocity of the stellar winds and of the SNIa ejecta), and finally the velocity dispersion tensor of the source field. In the case of SNIa’s mass input, the thermalization p of random motions is usually neglected, due to the high velocity of the ejecta us = 2ESN /1.4M ' 8.5 × 103 km s−1 , far above the typical value of the velocity dispersion in ETGs (' 150 − 300 km s−1 ). In addition, the ejecta is thermalized via shocks in the ISM (due to its high velocity), so that internal energy einj of the ejected material is negligible with respect to u2s . The opposite applies to stellar winds: a typical red giant star injects mass in the ISM via winds with a speed of few 10 km s−1 (Parriott & Bregman 2008b), one order of magnitude lower than the velocity dispersion of a typical ETG, so that the contribution of the winds internal energy and kinetic energy is usually ignored. In our work we neglect einj term for both stellar wind and SNIa, the us term of stellar winds, but we consider that of SNIa ejecta and all the energy injection terms related to the stellar kinematic; this leads to the present form of the energy input from the ageing stellar population: i u2 ρ˙ h ˆϕ − uk2 + Tr(σ 2 ) , E˙ = ρ˙ SN ηSN s + kvϕ e 2 2

(2.18)

where we adopted a thermalization efficiency equal to 0.85 for the kinetic energy input from SNIa (e.g. Thornton et al. 1998; Tang & Wang 2005). An important point follows from the adopted trend of mass and energy injections. Given that LSN is typically larger than the integrated contribution from the thermalization of random and ordered motions, due the high velocity of SNIa ejecta, we can approximate the ˙ ρ˙ ' us αSN /α∗ ∝ t0.3 . This difference in the time dependence of specific heating as E/ the mass and energy inputs from the evolving stellar population produces a long-term time-increase of the specific heating of the input mass. In other words, the injected mass from the passively evolving stellar population is hotter and hotter as time increases.

12

2.2.2

The simulations

Star formation and SNe II heating

Star formation is often regarded as the most challenging problem in hydrodynamical simulations. A complete and self-consistent simulation of star formation in a galactic environment is nearly impossible, due to the extremely complex physics involved (magnetic fields, gas self-gravity, bipolar diffusion of the primordial cloud, molecular formation and dust shielding, to name a few; for a review see Larson 2011, 2012). For this reason, we adopted the Kennicutt (1998) recipe that reproduce the empirical Kennicutt-Schmidt law for star formation (Schmidt 1959; Kennicutt 1998), which links the star formation rate (SFR) to the local ISM density. Following Ciotti & Ostriker (2012) and references therein, the local star formation rate per unit volume is given by ρ ρ˙ SF = ηSF , tSF = max(tcool , tdyn ), (2.19) tSF where 0.01 6 ηSF 6 0.1 is the star formation efficiency, and the involved time-scales are the cooling (the ISM must be cold to form a star) and the Jeans one (the ISM must reach a critical density in order to collapse by self-gravity): r 3π tcool ≡ E/L , tdyn ≡ , (2.20) 32Gρ where G is gravitational constant. When star formation takes place in a computational cell, the cold ISM is removed from the numerical grid and a new stellar population appears. While the mass sink term related to the star formation episode is straightforward, being equal to the negative of ρ˙ SF , the momentum and energy sink are not trivial, since they depend on the cold cloud collapse process. We assumed that the cloud gas sink velocity is that of the local ISM, therefore the sink term in the momentum equation A.49 is identically zero, while the energy sink term is: ηSF E E˙ SF = , tSF

(2.21)

where E is the ISM internal energy density. For a given mass ∆M∗ = ρ˙ SF ∆t∆V of new stars formed in a computational cell of a given volume ∆V , by assuming a Salpeter IMF (x > 1, M ≥ Minf = 0.1 M ) dN ∆M∗ = (x − 1) 2 × dM Minf



M Minf

−1−x ,

(2.22)

the number of stars having a mass greater than MII = 8 M , that will explode as Type II SNe, are:    1 Minf x ∆M∗ ∆M∗ N = 1− ' 7 × 10−3 , (2.23) x MII Minf M where x = 1.35 (for details see Ciotti & Ostriker 2012). By assuming that each SNII leaves a remnant of 1.4 M , the total mass injected by SNIIe formed in a computational cell is ' 0.2∆M∗ , while they will injects a total amount of energy equal to N ESN ηSN εII = ' 3.9 × 10−6 ηSN . (2.24) ∆M∗ c2

2.2 The input physics

13

To recover the mass source term ρ˙ II for Type II SNe we have to consider that, on one hand, a given star formation episode generates SNIIe that will exponentially decline their mass and energy output on a time-scale τII = 2 × 107 yr. On the other hand, at a certain time t1 during the evolution of the very first SNIIe formed by the code, another star formation episode may take place, which forms younger SNIIe that in turn will eject into the ISM a mass of hot gas equal to 0.2ρ˙ SF (t1 ). This process can be formalized as: dρ˙ II ρ˙ II 0.2ρ˙ SF + . (2.25) =− dt τII τII By solving equation (2.25), the mass source term related to SNIIe is obtained, while the energy injection due to the ejecta thermalization is ρ˙ II εII c2 . As for the initial stellar population, the random and streaming motions of the SNIIe are thermalized via inelastic mixing of the ejecta in the ISM. However, in the case of the primordial stars, the stellar kinematics is determined by the initial conditions, whereas the velocity dispersion and azimuthal motions of the new stars are, in principle, a function of the pre-existing cold cloud that generated them, and they can even possess a non-vanishing velocity on the (R, z) plane (e.g. in the case of a cold clump that forms stars while falling towards the galaxy centre)2 . However, since the natural places for star formation are the rotationally supported discs formed by cooled gas (see Chapter 3 and 4), we assumed that the every new stellar population formed inherits all the kinematical features of the initial stellar component. Thus, the total heating due to Type II SNe is i ρ˙ II h ˆϕ − uk2 + Tr(σ 2 ) . E˙ II = ρ˙ II εII c2 + kvϕ e (2.26) 2

2.2.3

Radiative cooling

The radiative cooling is implemented by adopting a modified version of the cooling law reported in Sazonov et al. (2005). We neglect the Compton heating/cooling and the photo-ionization heating, allowing only for line and recombination continuum cooling. As usual in numerical simulations of hot flows in ETGs, we impose a lower limit for the ISM temperature of 104 K by modifying the cooling function at low temperatures, in order to obtain a smooth decline of the radiative emission as the temperature of a cooling element approaches the limit value. With these assumptions, our version of the cooling function, derived from equation (A32) of Sazonov et al. (2005), becomes L = n2H Λ(T ), where nH is the hydrogen number density and  2   104 K −23 Λ(T ) = S1 (T ) + 10 a(T ) 1 − (erg s−1 cm3 ), T

(2.27)

when T ≥ 104 K and Λ(T ) = 0 otherwise (Figure 2.1).√The coefficients in equation (2.27) are, in the cgs system, S1 (T ) = −3.8 × 10−27 T and a=− 2

18 2

e25(log T −4.35)



80 2

e5.5(log T −5.2)



17 e3.6(log T −6.5)

2

.

(2.28)

Even in the case of vanishing meridional velocity for a new stellar population, it is very unlikely that a star is formed directly on a circular orbit; thus the process of asymmetric drift may take place for every stars formed during the entire simulation.

14

The simulations

Λ(T) (erg s−1 cm3)

10−21

10−22

10−23

10−24 104

105

106 T (K)

107

108

Figure 2.1: The bolometric cooling function Λ(T ). For temperatures lower than 104 K the gas does not lose energy via radiative emission.

2.3

The contribution of stellar kinematics to the ISM energetics

One of the main goal of this Thesis is to study the effects of flattening and ordered rotation on the ETGs X-ray coronae. Analytical studies, based on global energetics arguments, showed that different and competitive effects should be taken into account (CP96, Pellegrini 2011, P13). However, the effects of rotation on the flow energetics are not trivial to predict from first principles. Indeed, the stellar population injects energy via SNIa explosions and thermalization of stellar motions (see equation 2.18). At variance of the stellar random kinetic energy Z 1 Lσ ≡ ρ˙ Tr(σ 2 ) dV, (2.29) 2 which is always supplied to the ISM, the thermalization of stellar ordered (streaming) motions, which provides an energy input per unit time to the ISM of Z 1 Lv ≡ ρ˙ kv − uk2 dV 2 Z Z (2.30) 1 1 2 2 2 = ρ(u ˙ R + uz ) dV + ρ(v ˙ ϕ − uϕ ) dV = Lm + Lϕ , 2 V 2 V depends on the relative motion between stars and ISM, thus it cannot be predicted a priori from the knowledge of the galaxy structure and kinematics. Note that in ˆϕ , and Lm and Lϕ are respectively the energy input rate our galaxy models v = vϕ e

2.4 The code

15

due to the ISM velocity in the meridional plane (R, z), and to the relative velocity of stars and the ISM in the azimuthal direction. To parametrize the amount of thermalization of the stellar streaming motions, and to compare the results of the hydrodynamics simulations with the global energetic estimates, we introduce the thermalization parameter γth ≡ Lv /Lrot ,

(2.31)

Z

(2.32)

where Lrot

1 ≡ 2

ρv ˙ ϕ2 dV

is the energy input per unit time that would be injected in a galaxy with an ISM at rest (i.e., u = 0), due to thermalization stellar streaming motions. Note that γth is undefined (formally, it diverges) for velocity dispersion supported models, and can be very large for slow rotators and/or for gas flows with large velocities in the meridional plane (as in the case of galactic winds). Using these definitions, the total energy supplied to the ISM due to thermalization of stellar motions can be written as Lkin ≡ Lσ + Lv = Lσ + γth Lrot .

(2.33)

Of course, Lσ decreases when increasing the galaxy rotational support, at fixed galaxy structure. Note that, if γth = 1, the virial theorem assures that at fixed galaxy structure, Lkin is independent of the level of rotational support. In the other extreme case, if γth = 0, the gas is injected everywhere with the same local velocity of the ISM, and then Lkin decreases for a larger rotational support. However, ordered rotation acts also in a competing way, i.e., it tends to unbind the gas; therefore, when rotation is unthermalized, the ISM is less heated but it is also less bound. Which of the two competing effects dominate can be inferred only with the aid of hydrodynamical simulations: a goal of this Thesis is to determine which of the two competing effects dominate, and to obtain estimates of γth that can be used in analytical works (P13). Finally, all the luminosities defined above can be converted into equivalent temperatures as Z Z µmp µmp Tσ = ρ∗ Tr(σ 2 ) dV, Trot = ρ∗ vϕ2 dV, (2.34) 3kB M∗ 3kB M∗ where kB is the Boltzmann constant, mp is the proton mass and M∗ is the total stellar mass, so that Tkin = Tσ + γth Trot ;

2.4

γth Trot = Tm + Tϕ .

(2.35)

The code

The ZEUS-MP 2 code (Hayes et al. 2006) is a widely used general purpose magnetohydrodynamical, Eulerian, operator splitting, fixed mesh, second order upwind code which operates in one, two and three dimensions in Cartesian, spherical and cylindrical coordinates. In the following we present the hydrodynamical equations solved by the code, and the adopted numerical schemes.

16

The simulations

Coordinate system

x1

x2

x3

dV1

dV2

dV3

g31

g32

Cylindrical Spherical

z r

R θ

ϕ ϕ

dz r2 dr

R dR sin θ dθ

dϕ dϕ

1 r

R sin θ

Table 2.1: ZEUS numerical grid formalism. In the coordinate independent formalism, the three spatial coordinates are denoted as (x1 , x2 , x3 ), and the total volume of a single cell is dV1 dV2 dV3 (for details see Stone & Norman 1992; Hayes et al. 2006).

2.4.1

The hydrodynamical equations

The code has been modified to take into account all the source and sink terms listed in the previous Section, thus the equations of the hydrodynamics solved by the code are ∂ρ + ∇·(ρu) = ρ˙ SN + ρ˙ ∗ + ρ˙ II − ρ˙ SF , ∂t ρ

∂u ˆϕ − u), + ρ (u · ∇) u = −∇p − ρ∇Φtot + (ρ˙ SN + ρ˙ ∗ + ρ˙ II )(vϕ e ∂t ∂E + ∇·(Eu) = −p ∇· u − L + E˙ + E˙ II − E˙ SF , ∂t

(2.36) (2.37) (2.38)

where ρ, u, E, p, Φtot , and L are respectively the ISM mass density, velocity, internal energy density, pressure, total gravitational potential, and bolometric radiative losses per unit time and volume. As usual, the gas is assumed to be an ideal monoatomic fully ionized plasma, so that p = (γ − 1)E, where γ = 5/3 is the adiabatic index. The chemical composition is fixed to solar (µ ' 0.62), and the gas self-gravity is neglected.

2.4.2

Numerical methods

The additional physics modules have been implemented with an operator-splitting, Forward Time Centered Space (FTCS) differencing scheme by adding a new sub-step in the ZEUS source step for each source/sink term, in conformance with the ZEUS numerical methods. The transport step is unmodified with respect to the public version. The discretization of equations (2.36)-(2.38) on the ZEUS grid (see Figure 2.2 and Table 2.1) with the FTCS is straightforward, e.g. in the case of the continuity equation: n ρn+1 i,j,k − ρi,j,k = ρ˙ SNi,j,k + ρ˙ ∗i,j,k + ρ˙ IIi,j,k − ρ˙ SFi,j,k , (2.39) ∆t where ρni,j,k and ρn+1 i,j,k are the old and new density values at (x1bi , x2bj , x3bk ). The same simple scheme applies also for the momentum and energy equations, except for the cooling term in L in equation (2.38), see below. In addition, the mass injection rate of SNII is calculated by discretizing equation (2.25) as ρ˙ n+1 IIi,j,k

 =

∆t 1− τII



ρ˙ nIIi,j,k +

0.2∆t n ρ˙ , τII ∗i,j,k

(2.40)

2.4 The code

17

x1a(i − 1)

x1a(i + 1)

x1a(i)

dx1a(i + 1)

dx1a(i)

dx1a(i − 1)

x1a(i + 2) x2a(j + 2)

x2b(j + 1)

dx2a(j + 1)

i, j + 1 i, j + 1 i, j + 1

x2a(j + 1)

dx2b(j + 1)

x2b(j)

i − 1, j i − 1, j i, j

dx2b(j)

i − 1, j

x2b(j − 1)

i, j

i + 1, j i + 1, j

i, j

i + 1, j

dx2a(j)

x2a(j)

dx2a(j − 1)

i, j − 1 i, j − 1 i, j − 1

x2a(j − 1)

dx1b(i + 1) dx1b(i) x1b(i + 1) x1b(i) x1b(i − 1) Figure 2.2: Spatial position of an arbitrary point on a (x1 , x2 ) slice. The discretized variables on the ZEUS computational domain are denoted by the staggered coordinate vectors x1ai or x1bi , and x2aj or x2bj . The a-mesh values are located at zone edges; the b-mesh values are located at zone centres. The dots denote the centring and indexing of zone-centred (scalars) and face-centred (vectors) variables on the mesh. As a consequence, the scalar quantities ρ, Φ, E, ρ∗ , uϕ , σ and σϕ are centred in (x1bi , x2bj ), whereas the first and second components of the velocity field are centred in (x1ai , x2bj ) and (x1bi , x2aj ), respectively. where ∆t is the simulation global time-step. Equation (2.40) allows to save computational storage, since it requires only the information at the previous time-step, thus avoiding the storing of the entire star formation history in every gridpoint at every time. Due to the ZEUS explicit scheme, the global hydrodynamical time-step ∆t takes into account the Courant–Friedrichs–Lewy stability condition imposing a minimum value ∆thyd (equation 60 in Hayes et al. 2006). Our input physics leads to the introduction of additional characteristic times, associated with the injection of mass and energy: ∆tρ =

ρ , ρ˙ SN + ρ˙ ∗ + ρ˙ II + ρ˙ SF

∆tc =

E , L

∆th =

E E˙ + E˙ II + E˙ SF

.

(2.41)

The momentum time-step requires a special treatment: if a non-rotating blob of gas

18

The simulations

400

dR (pc)

300

200

100

0.01

0.10

1.00 R (kpc)

10.00

100.00

Figure 2.3: Numerical grid resolution in the radial direction (the same is adopted for the z direction). Every dot is a numerical gridpoint of the computational mesh. hits a grid characterized by a rotating stellar population, a normal prescription as the ones above returns a null time-step. Thus, after several tests, we adopted a floor value ∆tmin = 1 Myr for the momentum time-step only, obtaining −2 −2 ∆t−2 m = ∆tmin + ∆tρ

(vϕ − uϕ )2 . u2ϕ

(2.42)

As a consequence, the global time-step is Ccfl , ∆t ≡ q −2 −2 −2 −2 + ∆t + ∆t + ∆t ∆t−2 + ∆t m c ρ h hyd

(2.43)

where Ccfl is the Courant coefficient (set to 0.5), and the minimum value of ∆t over the numerical grid is considered. While almost all the integration of equations (2.36)-(2.38) is performed by using an explicit temporal advancement (as prescribed by FTCS), for the integration of the cooling function we tested two different numerical algorithms: the fully explicit Bulirsch–Stoer method and the fully implicit Bader–Deuflhard method (Press et al. 1992). All the results presented in this work are based on the fully implicit algorithm, since it is far less computational time-consuming, and it allows to drop the cooling time-scale from equation (2.43), while giving the same global evolution of the hot gas flows (as proved with several tests).

2.4 The code

2.4.3

19

Code setup and outputs

The code is used in a pure hydro, 2.5D axisymmetric configuration (i.e. every physical quantity depends only on two spatial variables although the simulation is three dimensional) with a non uniform (logarithmic) computational mesh (R, z) of 480 × 960 gridpoints, having a resolution of ' 90 pc in the first 10 kpc from the centre, as shown in Figure 2.3. Reflecting boundary conditions were set along the z-axis, while on the outer edge of the simulated box the fluid is free to flow out of the computational grid. In this work we adopted a cylindrical grid in order to better resolve the regions near the equatorial plane, where a cold disc may form. Clearly, such a choice is quite expensive in terms of computational time, as more gridpoints than in the spherical case are needed, in order to maintain the shape of the grid reasonably regular with a logarithmic spacing. We verified by performing several tests that the code provides an excellent conservation of total mass and energy, where the former is given by II ∗ Mgas (t) + Mesc = Minj + Minj − Mnew ,

(2.44)

where Mgas and Mesc are the gas mass present into the simulation box at a given ∗ , M II time t and the cumulative mass escaped, Mnew inj and Minj are the cumulative mass of new stars, the cumulative mass injected by the passively evolving stellar population and by the SNIIe events, respectively; e.g. Z tZ Minj =

(ρ˙ SN + ρ˙ ∗ ) dV dt, t0

(2.45)

V

being t0 the simulation initial time. The conservation of energy can be written as Z

dEtot (E˙ − L + ρΦ) ˙ dV = + dt V

Z S

! p kuk2 ˆ dS, e+ + + Φ ρu · n ρ 2

(2.46)

where E˙ is the total injection/sink of energy, ρ˙ is the totalR injection/sink of mass, e = E/ρ is the ISM internal energy per unit mass, Etot = (e + kuk2 /2 + Φ)ρ dV , and the two integrals are extended over the whole numerical grid and its boundary, respectively. The hydrodynamical fields are saved every 100 Myr, while grid-integrated quan∗ , M II , the hot and cold gas mass (M tities, such as Minj , Mesc , Mnew hot , having inj 6 4 T ≥ 10 K and Mc ≤ 10 K), the integrated SNIIe luminosity LII , Lrot , Lσ , Lv , Lm , Lϕ , the SFR, the mean time of the SFR peak (hti∗ ), the X-ray emission in the 0.3–8 keV Chandra band LX , and the X-ray emission weighted temperature (TX ) are sampled with a time resolution of 1 Myr. The mean time of the SFR peak is calculated as Z t 1 hti∗ = ∗ (t0 − t0 )SFR(t0 ) dt0 dV, (2.47) Mnew (t) t0 while LX and TX are calculated using the thermal emissivity εX ≡ 10−14 ne ni (ΛX1 (T ) + ΛX2 (T )) erg s−1 cm−3

(2.48)

20

The simulations

Table 2.2: X-ray emissivity coefficients, see equations (2.49)-(2.50). Λ1X (T ) a b c d e f g h i l m n o p q r

100

36.6899 × 8.19020 × 100 3.69232 × 10−8 1.48406 × 10−9 1.13656 × 10−11 4.04378 × 100 3.43352 × 100 7.01936 × 10−9 5.4211 × 10−9 1.51675 × 10−9 1.07104 × 10−9 1.65655 × 10−9 1.41044 × 10−9 3.71927 × 10−2 1.20299 × 100 1.73516 × 100

Λ2X (T ) 1.42382 × 10−8 7.26800 × 100 2.71451 × 10−10 3.18379 × 100 6.20760 × 10−3 4.51332 × 100 2.12548 × 100 5.47970 × 10−4 5.00962 × 10−3 2.39339 × 10−1

over 0.3–8 keV emission of a hot, collisionally ionized plasma, using the spectral fitting package xspec3 (spectral model apec, Smith et al. 2001, see also Figure 2.4), where ne and ni are the electron and ion number density. The two functions Λ1X (T ) and Λ2X (T ) are expressed (in erg s−1 cm3 ) by the following equations:  2 aT − bT + 1   ,   cT 2 − dT + e     2 f T − gT + 1 Λ1X (T ) ≡ , 2 − iT + l  hT      mT 2 − nT + o   , pT 3 + qT 2 − rT + 1  b  aT ,     d Λ2X (T ) ≡ cT ,    4    eT + f T − g 10−10 , hT 4 + iT 2 + lT + 1

0.03 < T < 0.3 keV, 0.3 < T < 0.7 keV,

(2.49)

0.7 < T < 16 keV,

0.045 < T < 0.17 keV, 0.17 < T < 0.7 keV,

(2.50)

0.7 < T < 16 keV,

where all the temperatures are expressed in keV and the coefficients are listed in Table 2.2. Thus: Z Z 1 T εX dV, (2.51) LX = εX dV, TX = LX where the integration extends over the whole computational grid. Finally, in the data 3

http://heasarc.nasa.gov/xanadu/xspec/

21

ΛX(T) (erg s−1 cm3)

2.4 The code

10−23

10−24 105

106

107 T (K)

108

Figure 2.4: The 0.3–8 keV specific emissivity ΛX (T ) = 10−14 (Λ1X (T ) + Λ2X (T )). post-processing we generate the angle averaged profile of the density distribution4 Z π ρr (r) = ρ(r, θ) sin θ dθ, (2.52) 0

the X-ray surface brightness and temperature maps for an edge-on5 projection Z R=∞ p ΣX (x, z) = εX (R, z) d R2 − x2 , (2.53) R=x

1 Tp (x, z) = ΣX (x, z)

Z

R=∞

p T (R, z)εX (R, z) d R2 − x2

(2.54)

R=x

where we placed the y axis of the Cartesian reference system coincident with the line-of-sight, so that the x − z plane coincides with the plane of the sky. 2.4.3.1

Software engineering

The code is parallelized via MPI6 domain decomposition and has its roots in the venerable (and still used) ZEUS 2D (Stone & Norman 1992), written in FORTRAN 77 in the late ’80s. Indeed, ZEUS-MP 2 is the last born of a copious family of hydro codes, all sharing the same common ancestor, therefore a lot of old (although solid) 4

In the case of simulations performed in cylindrical coordinates, the field is pre-interpolated on a high-resolution spherical grid. 5 If necessary, the integrand is pre-interpolated on a high-resolution cylindrical grid. 6 http://www.mpi-forum.org/

22

The simulations

code is present and many modifications have been done over more that 20 years. During this time, the entire world of informatics has been changed many times, especially its High Performance Computing (HPC) branch, on both software and hardware sides. On one hand, workstations with tens of processors are now common, and state-of-art supercomputers are made up by ten thousands of cores, everyone with their dedicated hardware, optimizing compiler and numerical libraries. On the other hand, now numerical codes are required to be able to run in parallel on hundred of processors, in order to efficiently exploit the most recent supercomputing resources. The public version of the code is written in a mixture of FORTRAN 77, 90 and non-standard extensions, which are no longer supported by many compiler. A conspicuous part of the coding has been spent in removing such non-standard features, in order to (1) ensure an high level of portability among every architecture and (2) employ the new automatic features of the modern compilers, such as autovectorization. As the rest of the code, the new implementation written on purpose for this work strictly follows latest Fortran standards (the latest, and “modern” standards are Fortran 95, 2003 and 2008, although the code uses only two specific instruction of the latter). The code has been compiled and tested on an huge scale of hardware, from a simple netbook, to workstations, to a BluGene/Q architecture, spanning from one to thousands of processors, as well as with the Intel, GNU and IBM compilers. In addition, the code hydrodynamical output has been rewritten from scratch, by implementing a parallel output library based on the HDF5 file format7 (the de-facto standard on HPC systems), which allows multiple processors to write directly a single binary self-describing file with a single collective call, thus minimizing the I/O overload on HPC systems.

7

http://www.hdfgroup.org/

Chapter 3

The effects of stellar dynamics on the X-ray emission of S0 galaxies Negri A., Ciotti L., Pellegrini S., 2014, MNRAS, 439, 823 In this first exploratory work, we focus on the effects of galactic shape and internal kinematics on the X-ray emitting haloes of ETGs. In particular, our aim is to investigate the observed X-ray under-luminosity of flat rotating galaxies with respect to rounder objects of similar mass (see Chapter 1). The proposed theoretical explanations can be classified in two categories: the energetic ones, and the hydrodynamical ones. The former can be divided into flattening effects and rotational effects: the first explains the low LX of flat rotating galaxies with a low binding energy of the gas due to a shallower potential of such systems, while the latter assumes that the injected mass maintains the streaming motion of the stars, thus the ISM is less bound and it lacks the thermalization of the stellar streaming energy (Lv in equation 2.30). In the hydrodynamical explanations, LX is lower in rotating galaxies since the centrifugal barrier keeps the galaxy central regions at lower density with respect to non-rotating systems, where the ISM is free to flow directly towards the centre and a hot, denser and brighter atmosphere is formed. In order to investigate which of the above effects determine the under-luminosity of flat and rotating galaxies, we simulated the evolution of the hot ISM in flat ETGs, modelled as two-component, self-consistent, axisymmetric realistic S0 galaxies. Various degrees of rotational support, dark matter and counter-rotation amounts are considered. The stellar distribution consist of a Miyamoto & Nagai (1975) profile, whereas the dark matter halo is described by an Einasto (1965) profile. The total gravitational field, the ordered and random stellar motions are computed by solving the Poisson and the Jeans equations with the code described in P13. However, in order to reduce the number of free parameters, we focus only on the effects of rotation and counter-rotation on the hot halo by keeping constant the shape of the model galaxy and its stellar mass, leaving the study of galaxy flattening and stellar mass dependence in the following Chapter. The Chapter is organized as follows. In Section 3.1 the structural and dynamical tuning of galaxy models is described, while in Section 3.2 we expose our results, commenting the hydrodynamical features, the X-ray luminosity LX and the X-ray luminosity weighted temperature TX . Finally, Section 3.3 summarizes our conclusions.

24

The effects of stellar dynamics on the X-ray emission of S0 galaxies

3.1

Setting the galaxy models

For reference we adopt a galaxy model tailored to reproduce the main structural properties of the Sombrero galaxy (M 104, of morphological type Sa), taken as a representative case of a flat and rotating galaxy; at this stage, though, we are not concerned with reproducing in detail the properties of the X-ray halo of Sombrero (but see Section 3.2.4). The adopted stellar distribution consist of a Miyamoto & Nagai (1975) profile, while the dark matter halo is described by an Einasto (1965) profile, see equations (2.3)-(2.8), respectively, and Section 2.1 for the galaxies building procedure. In all models the stellar distribution is kept fixed, as we focus on the effects of galaxy rotation only. For a given stellar and DM halo mass model, we consider four different cases of kinematical support for the stellar component: the isotropic rotator (IS, k = 1), the fully velocity dispersion supported case (VD, k = 0), the counter-rotating disc (CR), and a velocity dispersion supported system with an inner rotating disc (RD). In order to build the CR and RD models, the following functional form for the Satoh parameter has been adopted k(R, z) = kext +

ρ∗ (R, z) (kint − kext ), ρ∗ (0, 0)

(3.1)

to calculate vϕ and σϕ in equations (2.12)-(2.13), where ρ∗ is given by equation (2.3), but with a = 18 kpc and b = 4 kpc. This choice of the Satoh parameter leads to a very flattened rotating structure in the central regions of the galaxy, with k(0, 0) = kint , while k = kext at large radii. In particular, the CR models are obtained for kint = −1 and kext = 1, while for the RD models for kint = 1 and kext = 0. Therefore, the RD and CR models at large radii are similar to the VD and IS models, respectively. The stellar mass of Sombrero is M∗ ' 2.3 × 1011 M (Tempel & Tenjes 2006). By adopting Sombrero’s apparent blue magnitude of 8.98 (de Vaucouleurs et al. 1991), and a distance of 9.8 Mpc (Jardel et al. 2011), the resulting B-band luminosity is LB ' 3.8 × 1010 LB, . In order to reproduce the major photometric and kinematical features of M 104 as given by Jardel et al. (2011), under the assumption that the galaxy is an isotropic rotator, we fixed a = b = 1.6 kpc in equation (2.3), and n = 4 and rh = 52.8 kpc in equations (2.7)-(2.8). The resulting DM mass of Mh = 2.5 × 1012 M and ρc = 4.96 × 10−26 g cm−3 . We call this model ISi , where the subscript “i” stands for “intermediate halo”, for reasons that will be clear in the following. The rotational velocity and the azimuthal velocity dispersion of ISi model are given in Figure 3.1 (left panels), superimposed on the isodensity contours of the stellar distribution. For reference, the maximum stellar streaming velocity, in the equatorial plane, is vϕ,max ' 396 km s−1 at R ' 3 kpc, and σ = 245 km s−1 at the centre. Starting from ISi , we built three more models characterized by a different internal kinematics, but with the same stellar and DM halo distributions. In model VDi all the galaxy flattening is supported by azimuthal velocity dispersion [k = 0 in equations (2.12)-(2.13); Figure 3.1, right panels]; in the equatorial plane σϕ,max ' 401 km s−1 at R ' 12 kpc. In the counter-rotating model CRi (kint = −1 and kext = 1 in equation (3.1)), the equatorial negative and positive rotational velocity peaks are −155 km s−1 and 377 km s−1 , reached at R ' 2 kpc and R ' 24 kpc, respectively,

3.1 Setting the galaxy models

Name

25

Table 3.1: Stellar kinematics of the models. k kint kext Mrot Mcrot vϕ,max σϕ,max

Jz

ISi VDi CRi RDi

1 0 – –

– – -1 1

– – 1 0

1.00 0.00 0.72 1.00

– – 0.28 0.00

396 0.00 378 200

245 401 370 400

31.2 0.00 24.4 3.44

ISl VDl CRl RDl

1 0 – –

– – -1 1

– – 1 0

1.00 0.00 0.72 1.00

– – 0.28 0.00

332 0.00 303 178

223 340 327 338

25.5 0.00 19.5 2.99

ISh VDh CRh RDh

1 0 – –

– – -1 1

– – 1 0

1.00 0.00 0.72 1.00

– – 0.28 0.00

513 0.00 499 239

284 516 445 515

40.3 0.00 31.9 4.18

Notes. For each model the columns give: the k Satoh parameter of IS and VD models, the parameters kint and kext in equation (3.1) for CR and RD models, the rotating (Mrot ) and counter-rotating (Mcrot ) stellar mass normalized to M∗ , the maximum values of the stellar streaming velocity and of the azimuthal velocity dispersion in km s−1 , and the total angular momentum of the stars in 1073 g cm2 s−1 .

while the circle of zero rotational velocity is at R ' 4.9 kpc. In practice, CRi is similar to ISi in the external regions, but has a thin counter-rotating stellar disc in the inner region (Figure 3.2, left panels). For this model σϕ,max ' 370 km s−1 at R ' 5.5 kpc. Finally, in the RDi model an inner stellar rotating structure is present (kint = 1 and kext = 0 in equation (3.1)) with vϕ,max ' 200 km s−1 at R ' 3 kpc, while at large radii the galaxy flattening is supported by the velocity dispersion, similarly to what happens for VDi (Figure 3.2, right panels). Note that in all these models, by construction, the velocity dispersion fields σ = σR = σz are the same as in model ISi (see Section 2.1), and coincide with the field σϕ of ISi . In addition to these four models, hereafter referred to as having an intermediate DM halo mass, we built two more groups of models, where the DM mass is doubled with respect to the intermediate halo ones (Mh = 5 × 1012 M ; hereafter “heavy halo” models ISh , VDh , CRh and RDh ), and where the dark mass is halved (Mh = 1.25 × 1012 M ; hereafter “light halo” models ISl , VDl , CRl and RDl ). In the heavy and light halo models again the stellar distribution is kept fixed, as also n and rh ; in addition the four choices for the k(R, z) field corresponding to the IS, VD, CR and RD internal kinematic pattern are maintained (see Table 3.1). Summarizing, we followed the ISM evolution for a set of 12 models; a few more models with “ad hoc” modifications in the input physics have been also run, in order to test specific issues, as discussed in the following Sections. As usual in similar studies, the galaxy structure and dynamics are kept fixed during the simulations, and the initial conditions assume the galaxy is devoid of gas, as expected after the period of intense star formation, giving birth to the galaxy, is ended by the strong feedback from type II supernovae. In this way, simulations do

26

The effects of stellar dynamics on the X-ray emission of S0 galaxies

not start with an equilibrium ISM configuration, but instead the hot ISM distribution builds-up from stellar mass losses with time increasing. All the simulations start at an initial galaxy age of 2 Gyr, and the evolution of the gas flow is followed for 11 Gyr. The equations of hydrodynamics and the adopted numerical methods reported in Section 2.4. As we account only for a passively evolving stellar population, the terms related to star formation in equations (2.36)-(2.38)(ρ˙ II , ρ˙ SF , E˙ II and E˙ SF ) are not present. Merging and gas accretion from outside are not considered; black hole feedback is also ignored.

3.2

Results

We present here the main results of our investigation, focussing on a representative selection of models. The detailed features of each simulated flow of course depend on the specific galaxy model and input physics. While the parameter space is too large for a complete exploration, fortunately, the global behaviour of the gas is quite robust against minor changes of the input parameters; thus, a reasonable amount of computational time is sufficient to capture the different behaviour of the flows resulting from major variations in the structural parameters of the parent galaxy.

3.2.1

Hydrodynamics

All models, independently of their internal dynamics, evolve through two well defined hydrodynamical phases. Initially, all the ISM properties are characterized by an almost perfect symmetry with respect to the galaxy equatorial plane (z = 0). As time increases, the specific heating of the stellar mass losses increases (Section 2.2.1), and the velocity field becomes increasingly structured, in a way that is related to the particular internal kinematical support of the stellar component, as described below. A time arrives when the reflection symmetry is lost, and from this moment on it is never restored. In the following we present the main characterizing features of the gas flows in VD, IS, CR and RD models. A summary of the relevant integrated quantities at the end of the simulations is given in Table 3.2. 3.2.1.1

VD models

We present here the time evolution of the ISM of the non-rotating, fully velocity dispersion supported models VDi , VDl and VDh . Snapshots of various flow properties in the meridional plane (R, z) for the VDi model, for a selection of 9 representative times, are shown in Figs. 3.3 and 3.4. In particular, in Figure 3.3 the colours map the ratio ∆th /∆tc (equation 2.41), with green and purple corresponding to cooling and heating regions, respectively, while in Figure 3.4 we show the ISM temperature field for the same times. In both figures, the arrows represent the ISM meridional velocity field (uR , uz ). The major feature characterizing the flow of VDi is present from the beginning of the evolution: this is the degassing along the galaxy equatorial plane, due to the concentrated heating there, and accretion on the galaxy centre along the z-axis. Above the plane, on a scale of ' 10 kpc, the flow is characterized by large-scale

3.2 Results

27

Figure 3.1: Meridional sections of the galaxy rotational field vϕ (top) and of the stellar azimuthal velocity dispersion σϕ (bottom) for the ISi (left) and VDi (right) models. Note that, by construction, the velocity dispersion components σR = σz = σ of the two models coincides with σϕ of ISi . The stellar isodensity contours (solid lines) correspond to a density of 1 M pc−3 at the innermost contour, and to a density value decreasing by a factor of ten on each subsequent contour going outwards.

28

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.2: Analogue of Figure 3.1 for the counter-rotating CRi model (left), and for the velocity dispersion supported model with an inner rotating stellar disc, RDi (right). The rotating inner stellar disc is apparent in the top panels.

2.21 2.40 2.53 2.50 4.28 4.50 4.90 4.71 1.28 1.37 1.42 1.43

1.17 1.35 1.38 1.36 1.27 1.17 1.01 1.05 1.16 1.68 1.52 1.54

ISi VDi CRi RDi

ISl VDl CRl RDl

ISh VDh CRh RDh

2.46 4.74 3.55 4.48

1.54 2.60 2.10 2.47

1.85 3.31 2.58 3.13

Lσ (1040 erg s−1 ) (4)

2.28 0.00 1.20 0.26

1.06 0.00 0.50 0.14

1.46 0.00 0.73 0.18

Lrot (1040 erg s−1 ) (5)

0.24 6.40 1.60 2.63

0.41 3.08 0.14 1.21

0.58 3.38 1.03 1.99

LX (1040 erg s−1 ) (6)

0.55 0.70 0.55 0.77

0.32 0.39 0.37 0.42

0.40 0.50 0.40 0.55

TX (keV) (7)

Notes. (1) Name of the model (2) Hot ISM mass within the computational grid. (3) Escaped mass from the grid boundary. (4)–(5) Lσ (equation 2.29) and Lrot (equation 2.32). (6)–(7) X-ray ISM 0.3–8 keV luminosity LX and luminosity weighted temperature TX (equation 2.51). For reference, at 13 Gyr the total mass injected in the galaxy from the beginning by the evolving stellar population (stellar winds plus SNIa ejecta) is Minj = 2.23 × 1010 M , and the SNIa’s heating rate is LSN = 1.5 × 1041 erg s−1 .

(1)

Mesc (109 M ) (3)

Mhot (109 M ) (2)

Name

Table 3.2: Main outputs at the end of the simulations (13 Gyr)

3.2 Results 29

30

The effects of stellar dynamics on the X-ray emission of S0 galaxies

regular vortices (a meridional circulation). Due to the lack of centrifugal support, cold gas accumulates at the centre from the beginning. Loss of reflection symmetry of the flow occurs at t ' 4.5 Gyrs. After this time little evolution takes place, and overall the gas velocity field slowly decreases everywhere. The flow remains decoupled kinematically: the axial inflow-equatorial outflow mode persists in an essentially time-independent way for the entire run, with heated and outflowing ISM in the disc, and almost stationary gas above and below the galactic disc. The ISM temperature, after an initial phase in which the gas is hotter in the outflowing disc, quickly establishes on a spherically symmetric structure (Figure 3.4). From the beginning, at the centre (within ' 500 pc), the gas cools and forms a dense cold core. Outside this region, the temperature is steeply increasing, forming (within ' 1 Gyr) a spherical, hot region (T ' 9 × 106 K at the peak), of radius r ' 5 kpc; at larger radii, the temperature is slowly decreasing outward, keeping a spherical distribution. Outside the central cool core, the temperature is everywhere slowly increasing with time, due to the secular increase of the specific ISM heating due to the adopted SNIa time evolution, see equations (2.14)-(2.16). The major features described above for the VDi model are qualitatively independent of the DM halo mass, although important trends with Mh are clearly detected. For example, the time of loss of reflection symmetry in the flow properties increases from ' 4 Gyr (VDl ) to ' 4.5 Gyr (VDi ) to ' 5.2 Gyr (VDh ). In addition, a more massive DM halo tends to “stabilize” the ISM velocity field, in the sense that in the light-halo model VDl the meridional vortices are more pronounced, and the temperature maps (while still showing a spherically symmetric structure on average) are more structured and less regular. In particular, the average values of the ISM velocity are higher (at any given time) for lighter DM haloes, and the equatorial purple region in figures analogous to Figure 3.3 (not shown) is less symmetric. Finally, at any time the average ISM temperature is larger for increasing Mh . 3.2.1.2

IS models

The evolution of the flow in the family of isotropic rotators is more complicated than in VD models, as already found in DC98, albeit for different galaxy models and different input physics. Figures 3.5 and 3.6 show the flow properties of the ISi model, at the same epochs of Figs. 3.3 and 3.4 for the VDi model. The only similarities with VDi are the loss of reflection symmetry, that happens at now at ' 3.3 Gyr, and a systematic decline in the flow velocity for increasing time. Noticeable differences are instead apparent. First of all, in ISi there is the formation, since the beginning, of a cold and thin gaseous rotating disc in the inner equatorial galaxy region (with size ' 5 kpc), due to angular momentum conservation. This cold disc is quite stable, even though cooling instabilities from time to time lead to the formation of cold blobs detached from it. In general, the cooling (green) regions are significantly more rich in substructures than in VDi . In particular, a second major difference with respect to VDi is given by the presence of a cooling V -shaped region containing the equatorial plane whose vertex matches the outer edge of the rotating cold disc. In this region the gas is colder than in the rest of the galaxy (except for the cold rotating disc; see for example the snapshots at 2.7, 8.6 and 13 Gyr in Figure 3.6). This V-region becomes cyclically more or less prominent during the evolution; when

3.2 Results

31

Figure 3.3: Hydrodynamical evolution of the VDi model, for a selection of representative times. Arrows show the meridional velocity field (uR , uz ) of the ISM; their lenght is proportional to the modulus of the gas velocity, and is normalized to the same constant value in each one of the nine panels, and in all panels of the subsequent Figs. 3.4-3.10; thus, the evolution of the velocity field as a function of time can be followed for a single model, and compared to that of the other models in Figs. 3.4-3.10. For reference, the longest arrow in the top left panel corresponds to 171 km s−1 . Colours map the ratio of heating over cooling times, ∆th /∆tc , as defined in equation (2.41); green and purple colours indicate cooling and heating regions, respectively.

32

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.4: ISM temperature evolution for model VDi , at the same times as in Figure 3.3. The arrows indicate the velocity field in the meridional plane, and are normalized as in Figure 3.3. The colour-bar indicates the temperature values in K.

3.2 Results

33

it is more prominent, the gas in it is almost at rest in the meridional plane (i.e., it is fully supported by its rotational velocity uϕ ). Inside the V -shaped region, the gas is outflowing along the equatorial disc, while outside the ISM velocity field is organized in large meridional vortices. Note how this V -shaped region nicely maps the region of similar shape in Figure 3.1 (bottom left panel), where it is clear how the heating contribution from the thermalization of the stellar azimuthal velocity dispersion is missing with respect to the VDi model. A third major difference between ISi and VDi is represented by the long-term time evolution of the heated (purple) regions in Figure 3.5. In fact, in ISi it is apparent the fading of the heated equatorial disc region, accompanied by the appearance of a central heated region. In the VDi model, instead, cooling always prevails over heating in the centre. This difference is due to the lower gas density in the central regions of ISi with respect to VDi , which is produced by the angular momentum barrier of the ISi model, that prevents the gas from falling directly into the central galactic region. Moreover, in ISi , the infalling gas accumulates on the cold disc, and this further decreases the hot gas density in the central galactic region with respect to VDi . Thus, in the central galactic region, the cooling time keeps shorter in VDi than in ISi , during their secular evolution. This difference in the central gas density between rotating and non-rotating models, with the consequent secular heating of the central gas in ISi 1 , and the constant cooling of that in VDi , is at the base of a fourth major difference in the respective gas evolutions: the evolution is quite smooth in VD models, while it shows a cyclic behaviour in IS ones, as apparent from the panels relative to 7.3 − 8.6 Gyr in Figure 3.5, which describe a full cycle (a new cycle starts at t = 8.9 Gyr; in Section 3.2.3 we describe the evolution of other gas properties during a cycle). At the beginning of a cycle (t = 7.3 Gyr in Figure 3.5), the ISM in a central and almost spherically symmetric region becomes hotter and hotter, which produces a pressure increase in this region. This pressure increase causes an outflow from the centre, along the disc; as a consequence, the gas residing at R ' 10 kpc is compressed, increasing its density and lowering its temperature (Figure 3.6). The regular shape of the V -region is disrupted: some of the centrally outflowing gas breaks into its vertex, and succeeds in reaching out along the disc (t = 7.6 and 7.8 Gyr); some other gas circulates above and below the disc, in complex meridional vortices, compressing the gas there. In ISi the heating is not strong enough to establish a full and permanent degassing; the compression increases the cooling, until a maximum in the extension of the cooling (and of the low temperature) regions is reached, after which the flow reverts to the “original” state (e.g., shown by t = 8.9 Gyr in Figs. 3.5 and 3.6). These periodic changes in the ISM structure are mirrored in the evolution of LX , as discussed in Section 3.2.3. The ISM temperature distribution of ISi is shown in Figure 3.6. The major characteristic features of the hydrodynamical evolution are apparent. The cold thin rotating disc is visible, together with the embedding spherical region of hotter gas: it is interesting to note how the disc size corresponds to the radial extent of the hot region. The cold disc is dense (number density n & 10 cm−3 ) and azimuthally supported by ordered rotation, with peak values of uϕ = 420 km s−1 . The presence of the hot, spherical region, with radius matching that of the disc, is due to the 1

Recall that in all these models the specific heating is increasing with time.

34

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.5: Hydrodynamical evolution of the ISi model at the same representative times of model VDi .

3.2 Results

35

efficient way in which gas cools and joins the disc; this depletes the central galactic region of gas, thus the heating of the remaining gas is more efficient. In fact, within ' 3 kpc from the centre, the average gas density is ≈ 10−2 cm−3 (or less) in the ISi model, and ≈ 10−1 cm−3 (or more) in the VDi model, at least an order of magnitude larger. Also, within the same radius, excluding the cold core (for the VDi ) and the cold disc (for the ISi ), the average temperature is ' 107 K for the ISi , and lower (5 × 106 K) for the VDi (Figs. 3.4 and 3.6). Outside the central hot sphere, though, the temperature of ISi is everywhere lower than that of the VDi model, and the ISM is centrifugally supported by ordered azimuthal velocity. Maps of the Mach number show that the ISM velocity field is in general subsonic over the whole galaxy body. As a consequence, the X-ray emission is not associated with shocks; instead, inhomogeneities in the ISM usually cools more effectively (as those in the V-shaped region), and contribute to the total X-ray emission2 (see also Section 3.2.3 where LX and TX of all models are discussed). As for the VD models, a decrease in DM halo mass causes the ISM temperature overall to decrease, and the density and velocity fields become more and more rich in substructures. Due to the lower importance of angular momentum (a consequence of the reduction of the ordered stellar streaming velocity), the size of the cold disc decreases from ' 10 kpc (ISh ), to ' 5 kpc (ISi ), to ' 4 kpc (ISl ). Remarkably, the size of the hot spherical region is always the same as the size of the cold disc. The ISM rotational velocities in the V -shaped region also decrease for decreasing DM halo mass; instead, at late times, ISl presents a polar outflow, with velocities of the order of uz ' 250 km s−1 at ' 10 kpc above the equatorial plane. These polar outflows are significantly reinforced in test models in all similar to ISi , but with doubled SNIa rate. A change in DM halo mass has also some complex consequences, coming from the interplay between heating and binding energies at the galactic centre. The main characteristics of the global evolution of ISi and ISl are very similar, while the cyclic behaviour is far less prominent in ISh , and becomes almost absent after a few Gyr of evolution (see also Section 3.2.3). In fact, being the central potential well deeper in ISh than in the other models, the outflow velocity of the disc gas is lower, the effect of compression of the surrounding gas is also lower, and major cooling episodes, with the associated substructure in the flow density and velocity patterns, are absent. 3.2.1.3

CR models

We now focus on the first of the two special families of rotating galaxy models, namely the counter-rotating ones (CR). As described in Section 3.1, counter rotation is introduced in the IS family adopting a convenient functional form for the coordinatedependent Satoh parameter; in particular, we constructed the counter rotation so that vϕ = 0 at R ' 5 kpc, i.e., at the edge of the region where ISi develops the cold rotating disc. As anticipated in Chapter 1, this choice maximizes the possible effects of counter rotation, both from the energetic and the angular momentum points of view. As can be seen from Figure 3.7, the global behaviour of the CRi model is 2

Empirical evidence for the absence of shock-related X-ray emission is given in Section 3.2.3, where VD models are shown to be the most X-ray luminous, while showing the less structured velocity pattern. Cooling inhomogeneities instead affect the secular trend of the X-ray emission in IS models, never being able, though, to make rotating models more X-ray luminous than non-rotating ones (Section 3.2.3).

36

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.6: ISM temperature evolution of the ISi model, at the same times as in Figure 3.5.

3.2 Results

37

somewhat intermediate between those of the VDi and ISi models: in fact, although the V -shaped cooling region is still present, it is quite reduced with respect to that in ISi (even in regions where counter rotation does not have a direct effect), and a stronger galactic disc outflow along the equatorial plane takes place. Correspondingly, the cold gaseous central disc is smaller (with maximum size of ' 3 kpc), a consequence of the combined effects of a stronger local heating and the decrease of the local angular momentum of the ISM due to the mass injection of the counter-rotating stellar structure. Overall, however, the hydrodynamical evolution is similar to that of ISi , showing that the reservoir of angular momentum at large radii (where CRi and ISi are identical by construction) is the leading factor in determining the flow behaviour. In particular, the CRi velocity field is more rich in substructures than in VDi . The temperature evolution of CRi is presented in Figure 3.8. Again, as in ISi , the size of the cold disc strictly matches the size of the spherical region of hotter gas embedding the cold disc itself. In the azimuthal velocity field of the ISM, counter rotation is present at early times, when stellar mass losses are more important; this produces a region of hotter gas that is not present in the ISi model, and which is apparent at the radius of maximum stellar counter rotation in the first two temperature maps (t = 2.5 and 2.7 Gyr in Figure 3.8, cfr. with corresponding panels in Figure 3.6), as a lighter-coloured area along the equatorial plane, starting from ' 2 kpc. As time increases, however, the relative importance of the injected counterrotating gas decreases, and the rotational velocity of the ISM becomes dominated by the angular momentum of the ISM inflowing from the outer regions, and gas counter rotation is no longer present. A variation of the DM halo mass leads to the same systematic trends of the other families: the global ISM temperature decreases for decreasing Mh , the density and velocity fields become more structured, the linear size of the cold disc decreases (from 12 to 5 to 3 kpc in radius), and the gas rotational velocity in the V -shaped region decreases. Remarkably, a sustained polar outflow with uz ' 300 km s−1 at a height of ' 10 kpc above the equatorial plane develops in CRl by the present epoch, similarly to what happens for ISl . 3.2.1.4

RD models

We conclude with the family of the RD models, that are similar to the VD ones except for the presence of a rotating stellar disc in their inner equatorial region. Note that the rotating disc is the only source of angular momentum in this family. The hydrodynamical evolution of RDi is summarized in Figure 3.9, where the global similarities with VDi are apparent. In particular, at variance with ISi and CRi , the V -shaped region is now missing, while a small (' 2 kpc radius) cold disc is present, originated “in situ” by the stellar mass losses in the inner rotating stellar disc. The equatorial outflow is still present as in VDi , and the ISM velocities decrease as time increases. Also, in analogy with VDi and at variance with ISi and CRi , the central spherical heating region does not appear at late times. This shows again how the global evolution of a model is strictly linked to the amount of angular momentum stored at large radii, more than to the specific rotation in the central galactic regions. The temperature evolution of RDi is shown in Figure 3.10, where again the

38

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.7: Hydrodynamical evolution of the counter-rotating CRi model, at the same representative times of VDi .

3.2 Results

39

Figure 3.8: Temperature evolution of the counter-rotating CRi model, at the same representative times of Figure 3.7.

40

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.9: Hydrodynamical evolution of the velocity dispersion supported model with rotating inner disc, RDi , at the same representative times ofVDi .

3.2 Results

41

similarities with VDi are apparent. In particular, the temperature field is much less structured than in the rotating ISi and CRi models. As expected, almost no ISM azimuthal rotation is present in RDi , with the exception of some degree of rotation confined in the inner regions, where the small cold disc resides. The disc is embedded in a spherically symmetric region of hotter gas, similar to that present on a larger scale in the ISi model, and due to the same cause. The variation of the DM halo mass leads to the same overall changes as in the other families: at any given time, for decreasing Mh , the ISM temperature is lower and the hydrodynamical fields are systematically less regular, with higher average outflow velocities in the equatorial plane of the galaxy, while the extension and the rotational velocities of the small inner cold region decreases.

3.2.2

The thermalization parameter

As discussed in Section 2.3, one of the goals of this work is to measure the thermalization parameter γth , i.e., to estimate how much of the kinetic energy associated with ordered rotation of the stellar component is converted into internal energy of the ISM, equation (2.33). In fact, in addition to the obvious physical relevance of the question, reliable estimates of the value of γth as a function of the galaxy rotational status are useful in theoretical works (e.g., involving estimates of LX and TX based on energetic considerations, without simulations; CP96; Pellegrini 2011; Posacki et al. 2013a,b; Posacki 2014). The summary of the results for the three families IS, CR, and RD is given in Figure 3.11, where red, black, and green lines give the γth values for high, intermediate, and light DM haloes, respectively. The VD family is not considered, the R being 2 associated γth undefined (formally infinite, being Lrot = 0 and Lv = 0.5 ρ˙ kuk dV ). Note that in principle γth can be even larger than unity for galaxies with low rotation and thus low Lrot (as the RD models), or in cases of substantial counter rotation with high Lv (as for CR models). From Figure 3.11 a few common trends are apparent, that can be easily explained when considering the hydrodynamical evolution of the models. The first is that for all IS, CR and RD models the value of the thermalization parameter decreases for increasing DM halo mass. As γth = Lv /Lrot , this can be explained by the combination of two effects: the increase of Lrot with Mh , coupled with the decrease of the ISM velocity in the meridional plane and the increase of the azimuthal component of the ISM velocity in Lv = Lm + Lϕ (equation 2.30). The second common feature of all models is that fluctuations in the values of γth tend to increase for decreasing Mh , and this is due to the ISM velocity field becoming more rich in substructure for lighter DM haloes. Consistently with the hydrodynamical evolution, the fluctuations in the RD family are however smaller than in CR and IS models, as the complexity of the ISM velocity field is proportional to the amount of ordered rotation of the stellar component. In particular, the large fluctuations of γth in the intermediate and light halo models of the IS and CR families are due to the recurrent degassing events (with an increase of the velocity components uR and uz ) in the equatorial plane, as described in the previous Section. An important difference between the three classes of models is instead given by the average value of γth that, for the IS family, is much lower then for the CR and

42

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.10: Temperature evolution of the velocity dispersion supported model with an inner rotating stellar disc, RDi , at the same times of Figure 3.9.

3.2 Results

43

IS

CR

RD

γth

1.0

0.1

4

6

8 t (Gyr)

10

12

4

6

8 t (Gyr)

10

12

4

6

8 t (Gyr)

10

12

Figure 3.11: Time evolution of the thermalization parameter γth for the three families IS, CR, and RD: heavy, intermediate, and light DM haloes are shown with red, black, and green lines, respectively. RD families (for which, for the reasons explained above, γth can reach values larger than unity). Instead, the IS family is characterized by γth ' 0.1 − 0.3, a remarkably lower degree of thermalization. Being the large-scale kinematical support of the CR and IS families very similar, their significantly different γth must be explained by ˆϕ − uk2 term that originates in the counter-rotating disc of the CR the larger kvϕ e family. The range of values of γth for the IS family corresponds to the high-α models in P13, and it provides part of the explanation for the average lower X-ray luminosity of rotating models with respect to non-rotating ones of the same mass (a feature that is found in observed ETGs, see Chapter 1). In Section 3.2.3 and Section 3.3 we will return to this point for some additional considerations.

3.2.3

The X-ray observables: LX , TX , and ΣX

The time evolution of the observationally important ISM diagnostics LX and TX is shown for all models in Figure 3.12, and a list of LX and TX values at the end of the simulations (13 Gyr) is given in Table 3.2. Important similarities and differences, due to the different internal kinematical support and to the variable DM amount, are evident. Concerning similarities, in all models the ISM X-ray luminosity LX decreases with time, broadly reflecting the decrease of the hot gas content in the galaxies. Another similar behaviour is that, as already discussed in Section 3.2.1, within each family the X-ray luminosity weighted temperature TX increases when increasing the DM amount; in addition, TX increases with time (as more evident in the VDh and RDh models), due to the time evolution of the specific heating of the injected material (see Section 2.2.1). A major distinctive property of rotating models (IS and CR) with respect to non-rotating ones (VD and RD) is instead the presence of well defined oscillations in LX and TX . These oscillations are the result of the cyclic behaviour of the hydrodynamical evolution typical of IS and CR models (Sections 3.2.1.2 and 3.2.1.3). During each oscillation (see for example that corresponding to the large peak in LX at around 8.2 Gyr for ISi , mapped in Figure 3.5), LX and TX reach respectively a maximum and a minimum, all due to the onset of a cooling phase. At the beginning

44

The effects of stellar dynamics on the X-ray emission of S0 galaxies

of each cycle, the X-ray luminosity is low, the galaxy is filled with the gas coming from stellar evolution, and heated by the energy source terms. As the gas mass rises, the cooling becomes more and more efficient due to the compressional effect of the central outflow (Section 3.2.1.2), and LX increases too. This trend continues until a critical density is reached, such that the radiative losses dominate over the heating sources, and the gas catastrophically cools; at this point the peak in LX is produced, with the associated sharp decrement in TX . Finally, after the major radiative cooling phase has ended, the hot gas density and LX are low, and a new cycle starts. The global pattern of an oscillation is always the same, and governed by the mass injection and the radiative cooling rates. With time increasing, oscillations become more distant in time, since the refilling and the heating times become longer, due to the temporal decay of both the mass injection rate and the number of SNIa events (equations (2.14)-(2.16)). The second distinctive property of rotating models is that their LX and TX are always lower (for the same DM halo), than those of the models that are non rotating on the large scale (the VD and RD families). This is an important feature, also for its observational implications (Chapter 1), and thus deserves some consideration. This important effect of rotation, as that of causing a cyclic behaviour of the flow, originates in a different flow evolution that in turn is due not just to the different energetic input, but also to the different angular momentum of the gas at large radii. In fact, a first explanation of the X-ray under-luminosity and “coolness” of the IS family, that seems natural, lies in the lack of the σϕ -term in their Lσ (being σϕ replaced by the ordered rotational field of the stellar component, see Figure 3.1), and in the result that γth has low values (< 1, to be inserted in equation 2.18). However, this energy-based argument does not give the full explanation of the low LX values; a hint towards this conclusion is provided by the finding that LX is low also for the CR family, with similar γth values as for the RD one. We established definitively that the different energy input to the gas is not the sole explanation of the low LX and TX by performing some “ad hoc” experiments in the IS family. In practice, while retaining their internal dynamical structure, we modified the thermalization term in equation (2.18), replacing it with the full thermalization term of the VD family. Thus, from an hydrodynamical point of view, the ISM of these models still rotates as dictated by equation (2.37), but its energy injection is equal to that of the VD family3 . The results are interesting: on one side, LX and TX are higher than in the IS models; however, on the other side, LX and TX are still characterized by large oscillations (typical of rotating models, and absent in the VD family), and they are still lower than in the VD family. Having said this, one should also notice that, LX differs more, by comparing IS and VD models, than the whole of Lrot , which by itself shows that the energetic argument cannot account for the full gas behaviour (see Lσ and Lrot in Table 3.2). Therefore, these experiments prove that the X-ray under-luminosity and coolness of IS and CR models is not just due to a reduction of the injection energy in them, but -more importantly- to the global evolution of the ISM induced by ordered rotation. 2 In these tests the square brackets in equation (2.18) was substituted with kuk2 + vϕ + Tr(σ 2 ), 2 so that from the virial theorem the sum of the last two terms equals Tr(σ ) of VD models. Actually the heating in these modified IS models is even larger than in VD ones, because the ISM velocity of the former contains also a relevant rotational component uϕ . 3

3.2 Results

45

A few additional trends are shown by Figure 3.12. Time oscillations in LX and TX , for the rotating IS and CR models, become more and more important for decreasing Mh , reflecting the more structured density, velocity and temperature fields of the ISM for lighter DM haloes (see Sections 3.2.1.2 and 3.2.1.3). Another point is that the RD models always show the largest TX , for any DM halo; this is due to their central, small, hot region surrounding the small cold disc at their centres, a region that is not present in the VD class. Finally, in the rotating families (IS and CR), LX is systematically lower for increasing DM halo, while the opposite takes place for the two globally non-rotating families. This trend is due to the global angular momentum stored at large radii in IS and CR models, and its influence on the global behaviour of the flow: in rotating models, the increase of Mh corresponds to an increase of the total angular momentum, such that the galactic central region (where most of the LX comes from) is less dense of gas for larger Mh , due to the accretion of the gas on a larger cold disc. We finally discuss the edge-on appearance of the X-ray surface brightness maps ΣX at the end of the simulation (for the intermediate halo models; Figure 3.13). From the figure it is clear how large scale galaxy rotation leads to flatter and “boxy” X-ray isophotes, and to less concentrated X-ray emission, with respect to what is seen in non-rotating VD and RD models. Also the ISM density distribution is rounder in the VD, and more elongated along the equatorial plane in the IS models. In IS models, this elongation and the consequent “boxiness” (already found in Brighenti & Mathews 1996, DC98) are due to the rotational support on the equatorial plane, that prevents the gas from flowing inward. For what concerns a comparison with the optical surface brightness distribution Σ∗ of the parent galaxy (Figure 3.13), in non-rotating models ΣX is rounder than Σ∗ at all radii, reflecting the rounder shape of the total equipotential (mostly due to the DM). On the other side, rotation proves to be effective in determining a significantly flatter ΣX , that becomes more similar to Σ∗ , especially in the inner galactic region (within ' 10 kpc); however ΣX is never flatter than Σ∗ . ΣX becomes more spherical at large radii, and the boxiness decreases with radius, but it is still present at R ' 30 kpc.

3.2.4

Comparison with observed X-ray properties

As a general check of the reliability of the gas behaviour obtained from the simulations, we consider here a broad comparison with the observed X-ray properties of the Sombrero galaxy (whose structure was taken as reference), and of flat ETGs. In the Sombrero galaxy, diffuse hot gas has been detected in and around the bulge region with XMM-Newton and Chandra observations (Li et al. 2007, 2011b), extending to at least ' 23 kpc from the galactic centre, roughly as obtained here (Figure 3.13). The X-ray emission is stronger along the major axis than along the minor axis, and can be characterized by an optically thin thermal plasma with kT ' 0.6 keV, varying little with radius. The total 0.3–2 keV luminosity is LX = 2 × 1039 erg s−1 , and the hot gas mass is ' 5 × 108 M (Li et al. 2011b). The gas has a super-solar metal abundance, not expected for accreted intergalactic medium, thus it must be mostly of internal origin, as that studied in this work. In a simple spherical model for the hot gas, originating from internal mass sources heated by SNIa’s, a supersonic galactic wind develops for a galaxy potential as plausible for Sombrero, with an LX

46

The effects of stellar dynamics on the X-ray emission of S0 galaxies

far lower than observed (Li et al. 2011b). A flow different from a wind, and as found by our 2D hydrodynamical simulations, may provide the correct interpretation of the observed X-ray properties. Among the suite of models run in our work, an LX value comparable to that observed is reached at the present epoch by the IS models (Table 3.2); further, the TX of the ISh model is very close to the observed one, while ISi and ISl have a lower TX (note though that the LX and TX values in Table 3.2 refer to the whole computational grid, corresponding to a physical region larger than that used for the X-ray observations). Our work thus shows how it is crucial to account for the proper shape of the mass distribution (e.g., bulge, disc and dark matter halo), as well as for the angular momentum of the mass-losing stars, to reproduce the hot gas observed properties. For example, VD-like models predict LX larger by an order of magnitude, and inconsistent with those of Sombrero. Another feature clearly requiring angular momentum of the stars is provided by the observed X-ray isophotes, that show a boxy morphology in the inner regions (Li et al. 2011b), as obtained by our models only in case of rotation. Finally, note that the hot gas emission in Sombrero is lower than predicted by the best fit LX − LK correlation observed for ETGs (e.g. Boroson et al. 2011), as shown by Pellegrini (1999b, 2005, 2012). The X-ray luminosity could be reduced in Sombrero by the effects of rotation, as explained in Section 3.2.3. Moving to X-ray observations of S0 galaxies, a Chandra survey of their X-ray properties has been recently performed by Li et al. (2011b). They tend to have significantly lower LX than elliptical galaxies of the same stellar mass. While Li et al. (2011b) focussed on the possible cold-hot gas interaction to find an explanation (see also Pellegrini et al. 2012), we can suggest that rotation could have an important effect. A case S0 study is NGC5866 (Li et al. 2009), where the morphology of the hot gas emission appears rounder and more extended than that of the stars (the galaxy is seen edge-on), and again the X-ray isophotes have a boxy appearance in the inner region (see Fig. 1b in Li et al. 2009). However, the stellar mass is much lower than for Sombrero (M∗ ' 3 × 1010 M ), so we cannot directly compare LX and TX of our modelling with the observed hot gas properties of NGC5866 (the latter is just 7 × 1038 erg s−1 ).

3.3

Discussion and conclusions

In this chapter we studied the effects of the stellar kinematics on the ISM evolution of flat ETGs, focussing in particular on the effects of stellar rotation, leaving the study of rotation versus flattening effects in the following Chapter. We considered four different families of galaxies (based of a representative model tailored on the Sombrero galaxy), characterized by the same stellar distribution, and by a spherical DM halo of fixed scale-lenght, with three different total mass values. In the first family (IS), the galaxy flattening is entirely supported by ordered rotation of the stellar component, while the velocity-dispersion tensor is everywhere isotropic. In the second family (VD), the galaxy flattening is all due to stellar azimuthal velocitydispersion, i.e, the models are fully velocity-dispersion supported. The other two families are variants of the first two: in the CR family a thin counter-rotating stellar disc is placed in the central region of IS models; in the RD family, a rotating thin stellar disc is placed at the centre of the non-rotating VD models. The standard

3.3 Discussion and conclusions

IS LX (erg s−1)

1041

1040

0.8

0.7

0.7

0.6

0.6

0.5 0.4

0.5 0.4

0.3

0.3

0.2

0.2 3 4 5 6 7 8 9 10 11 12 13 t (Gyr)

3 4 5 6 7 8 9 10 11 12 13 t (Gyr)

CR LX (erg s−1)

1041

1040

0.8

0.7

0.7

0.6

0.6

0.5 0.4

0.5 0.4

0.3

0.3

0.2

0.2 3 4 5 6 7 8 9 10 11 12 13 t (Gyr)

RD

1040

0.8

TX (keV)

TX (keV)

LX (erg s−1)

1041

VD

1040

0.8

TX (keV)

TX (keV)

LX (erg s−1)

1041

47

3 4 5 6 7 8 9 10 11 12 13 t (Gyr)

Figure 3.12: Time evolution of LX and X-ray emission weighted temperature TX of the four families of models. Red, black, and green lines refer to the heavy, intermediate, and light DM haloes, respectively.

48

The effects of stellar dynamics on the X-ray emission of S0 galaxies

Figure 3.13: Edge-on 0.3–8 keV surface brightness of the ISM at 13 Gyr, for the intermediate halo models VDi , ISi , RDi and CRi ; the brightness values on the colourbar are given in erg s−1 cm−2 . Superimposed are the isophotes (Σ∗ ) obtained by projecting the galaxy stellar density distribution, starting from 104 M pc−2 on the innermost contour, and decreasing by a factor of ten on each subsequent contour going outwards. Note that the ΣX maps of VDi and RDi are dominated by a bright, round core, whereas the ISi and CRi maps show a boxy low-luminosity core.

3.3 Discussion and conclusions

49

sources of mass and energy due to a passively evolving stellar population are adopted, while we neglect star formation (and more in general diffuse mass sinks due to local thermal instabilities that are taken into account in Chapter 5), and feedback effects due to a central supermassive black hole. The simulations have been performed in cylindrical symmetry, and cover the ISM evolution for 11 Gyr. The main results can be summarized as follows. In all models, the ISM velocity and density at early times are symmetric with respect to the equatorial galactic plane, but soon this reflection symmetry is lost, and the hydrodynamical evolution of rotating models becomes more complicated than that of non-rotating ones. In VD and RD families the cooling gas tends to flow directly to the galaxy centre, while conservation of angular momentum leads to the formation of a cold rotating gaseous disc in IS ones. Moreover, the flow in rotating models is spatially decoupled (as already found in lower-resolution simulations, DC98), and shows large-scale meridional circulation outside a V -shaped region containing the galaxy equatorial plane, while inside this region the ISM is radially supported against the gravitational field by azimuthal rotation. These features are confirmed by a few runs with a number of gridpoints increased by a factor of four, reaching a spatial resolution of 15 pc in the central regions. In general, with increasing time the ISM velocity tends to decrease. The gas is not outflowing from the galactic outskirts in a significantly larger amount in rotating than in non-rotating families, as shown by the similar Mesc values in Table 3.2. Indeed, the major effect of ordered rotation in our models (that refer to a massive galaxy) is not that of making the gas less bound, but that of redistributing the ISM throughout the galaxy. A remarkable difference between globally rotating (IS and CR) and globally non-rotating (VD and RD) models is represented by large secular oscillations of LX and TX in the former (see also DC98, Negri et al. 2013). These oscillations are due to periodic cooling episodes in the V -shaped region of the rotating families, accompanied with degassing events along their equatorial plane (however limited to within few tens of kpc). Furthermore, in IS models a central hot region of radius comparable to that of the disc develops above and below the disc, due to the lower gas density there with respect to VD models. Outside of this hot sphere, the IS temperature is lower than that of VD models. The X-ray luminosity LX is largest for the velocity-dispersion-supported VD and RD families, and the highest X-ray emission-weighted temperatures are shown by the RD family, at any given DM halo mass. The strong rotators (IS and CR) are characterized by LX and TX significantly lower than for non-rotating models; for example, LX can be more than a factor of ten lower, at the same DM halo mass (see Table 3.2). In all cases, LX and TX are within the range of observed values for galaxies of similar optical luminosity and central velocity-dispersion. Sarzi et al. (2013), following CP96, suggested that the different LX of slow and fast rotators could be due to fast rotators, being on average flatter, and also more prone to loose their hot gas; here we find that the different LX is indeed due to a lower hot gas content of rotating systems, but this is not produced by a larger fraction of escaped gas, but instead by a larger amount of hot gas that has cooled below X-ray temperatures. In agreement with CP96 we also find that in low mass galaxies, generally tending to develop outflows, rotation favours gas escape. For increasing DM halo mass, the ISM velocity fields become more regular, with

50

The effects of stellar dynamics on the X-ray emission of S0 galaxies

less substructure; and TX increases, while LX behaves differently: in rotating models, LX decreases with increasing Mh , while it increases in non-rotating models. Rotating IS and CR models with light DM haloes at late times develop a polar wind. The (edge-on) X-ray isophotes are rounder than the stellar isophotes, as expected due to the round shape of the total gravitational potential. However, in rotating models the X-ray isophotes tend to be boxy in the inner regions, while in non-rotating models they are almost spherical. Note finally how the gas evolution and overall properties of the IS and CR families on one side, and of the VD and RD on the other, are remarkably similar, i.e., the presence of centrally rotating stellar discs does not alter significantly the global flow behaviour. In order to quantify the amount of galactic ordered rotation which is actually thermalized, we computed the thermalization parameter γth ≡ Lv /Lrot , that is the ratio of the heating due to difference between the streaming velocity of the stars and the ISM velocity, and the heating that would be provided by the stellar streaming if stars were moving in an ISM at rest. We found that γth is substantially less than unity in IS models, while it is of the order of unity or more in CR and RD ones; γth increases for lower DM contents. This shows that the different behaviour of VD and IS families is not entirely due to the different amount of thermalization of the stellar motions, but rather to the impact of angular momentum on the flow at large scales. In fact, despite of their different γth , all the main features of the IS family are still present in the CR one, including the well defined oscillations in LX and TX that are absent in VD and RD models. The parameter γth is used in works involving the global energy balance of the gas (e.g., CP96; Posacki et al. 2013b). Posacki et al. (2013b) showed that low values of γth go in the direction of accounting for the relatively low values of LX and TX in flat (and rotating) galaxies when compared to the values of their non-rotating counterparts. Also Sarzi et al. (2013) suggested that the kinetic energy associated with the stellar ordered motions may be thermalized less efficiently to explain why fast rotators seem confined to lower TX than slow rotators. The fact that γth is substantially less than unity in IS models could provide support to an energetic interpretation of the X-ray under-luminosity of flat and rotating galaxies, when compared to non-rotating ones of similar optical luminosity. The results for CR models, however, show that the lack of thermalization of ordered rotation cannot be the only explanation of the low values of LX : for these models γth is in fact of the order of unity, yet their LX is similar to that of the corresponding IS models, in which γth ' 0.1 − 0.2. Thus, even if the presence of a stellar counter-rotating thin disc can increase the thermalization of the ordered motions from 10-20 per cent in the isotropic rotators, up to 100 per cent or more, additional phenomena related to the global angular momentum (mainly stored at large radii) influence the behaviour of the flow, and produce a difference in LX and TX . Indeed, the conservation of angular momentum prevents the ISM to flow directly towards the centre, as in VD and RD models, so that the galaxy central regions are kept at low density with respect to the VD and RD models. This difference in the hot gas density distribution results in a systematic low LX and TX of rotating models, since the contribution to TX of the very hot gas in the central regions of rotating models is marginal (due to their low density) and more affected by colder gas located in the outer regions

3.3 Discussion and conclusions

51

(T ' 2 × 106 K). An additional glaring evidence for the important role of angular momentum is provided by the fact that, in rotating models (both IS and CR), LX decreases for increasing Mh , that is associated with an increase of galactic rotation. Summarizing the key points of this chapter, LX is significantly decreased not due to a larger degassing, or to a lower energetic input (due to a “missing” part in Lkin ), but by crucial angular momentum-related effects; pure energetic arguments cannot fully account for the changes in the overall gas properties (e.g., LX , TX ), and thus cannot solve the problem of the X-ray under-luminosity, and “coolness”, of rotating galaxies. In the following Chapter we expand our investigation by simulating a large set of ETGs characterized by variations of stellar mass, intrinsic flattening, dark matter distribution, in addition to various degree of rotational support. The present study can be relevant to the topic of black hole fuelling. One of the most debated aspects of SMBH accretion is how gas is carried to the centre of galaxies, especially in presence of rotation; another aspect is whether the source of fuel is a hot, roughly spherical atmosphere, from which accretion is almost steady, or it lies in cold material that sporadically and chaotically accretes (e.g., Novak et al. 2012; Russell et al. 2013; Werner et al. 2014). Related to these important issues, the present investigation shows that in velocity-dispersion-supported systems accretion is more hot and radial, with a large fraction of the total input from stellar mass losses flowing straight to the centre; in systems more supported by rotation, instead, the central density of hot ISM is lower, the mass accreted towards the centre is very small, and a cold rotating disc provides a large reservoir of cold gas, that can lead occasionally to clumpy multiphase accretion. Moreover, as anticipated in the Introduction, the presence of a counter-rotating structure affects the central feeding: the simulations show that, by reducing the amount of local angular momentum, accretion in the central grid is favoured with respect to what happens in pure isotropic rotators.

Chapter 4

The effects of galaxy shape and rotation on the X-ray haloes of ETGs - Numerical simulations Negri A., Posacki S., Pellegrini S., Ciotti L., 2014, MNRAS, 445, 1351 In this second work regarding the effects of galaxy shape and rotation on the X-ray emitting halo of ETGs, we expanded our previous results on the X-ray underluminosity of rotating galaxies (presented in Chapter 3). At variance with the previous investigation, where only variations of galaxies kinematical support have been considered, in this Chapter we take into account variations of galaxy shape, rotation and mass, by performing hydrodynamical simulations of a large set of self-consistent, state-of-the-art galaxy models characterised by different stellar mass, intrinsic flattening, distribution of dark matter, and rotational support. The dark matter haloes follow the NFW (Navarro et al. 1997) or the Einasto (1965) profile, while the stellar profile is a flattened de Vaucouleurs (1948). The galaxy flattening is supported by ordered rotation (isotropic rotators) or by tangential anisotropy. This work is a result of a joint research collaboration with Dr. Silvia Posacki, who built the entire set of axisymmetric galaxy models with the P13 code, and tailored their structural parameters to reproduce the observed properties and scaling laws of ETGs. The Chapter is organized as follows. Section 4.1 presents the procedure to produce and tune our flattened galaxy models, while in Section 4.2 we expose our results, commenting both the hydrodynamical features, the X-ray luminosity LX and the X-ray luminosity weighted temperature TX , in the same fashion as the previous Chapter. Finally, in Section 4.3 we discuss the results and the main conclusions are presented.

54

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

4.1

The setting up of representative galaxy models

In this large exploration, we adopted axisymmetric two-component galaxy models, where the stellar component can have different intrinsic flattening, while for simplicity the DM halo is kept spherical. In particular, the stellar profile is a flattened de Vaucouleurs (1948) (see Section 2.1, equations (2.1)-(2.2)), where we restrict the values of q in equation (2.2) to (1, 0.6, 0.3), corresponding to E0, E4 and E7 galaxies when seen edge-on. We produced two different sets of models, each of them characterised by the same DM profile: the Navarro-Frenk-White (Navarro et al. 1997) set and the Einasto (1965) one (again, see Section 2.1). In each of the two sets, we consider different families of models, built following the procedure described in P13 (Sections 3.3 and 3.4) and in Posacki (2014). Here we just recall the main steps to calculate their fundamental parameters, listed in Table 4.1. Each family is associated with a spherical galaxy, that we call the “progenitor”. The progenitor structural parameters are determined by assigning σe8 (the aperture luminosity-weighted velocity dispersion within Re /8), and then deriving the luminosity and effective radius Re of the galaxy from the Faber–Jackson and the size–luminosity relations (Desroches et al. 2007). Then, from a chosen stellar mass-tolight ratio, associated to a 12 Gyr old stellar population with a Kroupa initial mass function (Maraston 2005), the stellar mass M∗ is derived, and thus the photometric features of the progenitor are fixed. Finally, the parameters of the DM halo are determined in order to reproduce the assumed σe8 and fixing Mh /M∗ ' 20 (Behroozi et al. 2013). In the NFW set, these constraints produce rh ' 2Re , 22 . c . 37, and a DM fraction fDM within a sphere of radius Re of ' 0.6 for the spherical progenitors. For the Einasto set we fix n = 6, and we find that rh ' 20Re , and fDM ' 0.56 for the spherical progenitor. In each of the two sets we considered three values of σe8 for the spherical progenitors, i.e., 200, 250 and 300 km s−1 . Therefore, each of the two sets is made of 3 families of models, for a total of 6 spherical progenitors. Table 4.1 lists all the relevant parameters characterizing the progenitors galaxy models for both sets. The flattened descendants of each progenitor with intrinsic flattening of E4 (q = 0.6) and E7 (q = 0.3), are derived as follows. We produce two flattened models for each value of q. The first flattened model is called “face-on built” (hereafter FO-built), since, when observed face-on, its Re is the same as that of the spherical progenitor; this requires FO-built flattened models to be more and more concentrated as q decreases (ρ∗ ∝ q −1 ). The second flattened model instead, when seen edge-on, has the same circularized Re of the spherical progenitor, thus we call it “edge-on built” (hereafter EO-built); this property makes the EO-built models expand with √ decreasing q (ρ∗ ∝ q). Therefore, a spherical progenitor with a given value of σe8 produces four flat galaxies: two E4 models (FO and EO built), and two E7 models (FO and EO built). As a further step, in order to study the effects of galaxy rotation, we assume two kinematical supports for each flattened system: one corresponding to a velocity dispersion supported galaxy (VD models, k = 0), and the other one to an isotropic rotator (IS models, k = 1). In the flattening procedure the DM halo is maintained fixed to that of the progenitor. Note that our flattened models are representative of ETGs since they are consistent with their observed properties. We indeed checked for models lying outside the observed scatter of the scaling laws, but

4.2 Results

55

our adopted flattening procedure is quite robust in producing acceptable models, so that we retained all of them. Summarizing, from each spherical progenitor of given σe8 , eight flattened models are obtained (see Table 4.1), and we refer to this group of nine galaxy models as to a family. All models belonging to a family can be identified either by the σe8 value of the spherical progenitor, or by their stellar mass M∗ (or B luminosity), or DM halo mass; note however that while these last three quantities are kept constant within a family, the σe8 of the descendants varies. Indeed, the modification of stellar structure involves a change in the stellar kinematics, and so in the value of σe8 ; in particular, for our models σe8 decreases for increasing flattening (see P13 for a comprehensive discussion). Note that σe8 depends on the line-of-sight direction for non-spherical models; when quoting σe8 for the latter models, in the following, we refer to the edge-on projection. As in the previous chapter, we account only for a passively evolving stellar population, leaving the study of star formation effects in the following chapter.

4.2

Results

Here we present the main results of the simulations, focussing on the hydrodynamical evolution of a few representative models, and then on the global properties LX and TX for the two sets of models.

4.2.1

Hydrodynamics

For illustrative purpose, we present the hydrodynamical evolution of some selected EO-built models belonging to the NFW set. In particular, in the family derived from the E0 progenitor with σe8 = 250 km s−1 , we consider the velocity dispersion supported E7 model, and the corresponding E7 isotropic rotator. In Section 4.2.1.4 we summarize the main similarities and differences with the other models, as well as some considerations on the behaviour of the thermalization parameter γth . In general, as found in Chapter 3, the gas flows are found to evolve through two well defined hydrodynamical phases. At the beginning, all the ISM quantities (density, internal energy and velocity) are nearly symmetric with respect to the galactic equatorial plane. During the evolution, the velocity fields become more and more structured, until, after a certain time that depends on the specific model, the reflection symmetry is lost, and it is never restored. 4.2.1.1

The E0250 progenitor

The initial (t = 2.4 Gyr) and final (t = 13 Gyr) configurations of the ISM are shown in Fig. 4.1, where we show the meridional section of the ISM temperature (top panels), and the ratio of the heating and cooling time theat /tcool (bottom panels; green corresponds to a cooling dominated region while violet refers to a heating region). The arrows show the meridional velocity field. All the ISM physical quantities are stratified on a spherical shape, as a consequence of the galaxy spherical symmetry. A decoupled flow is soon established (t ' 2.4 Gyr), with an inflow in a round central region surrounded by an outflowing atmosphere.

56

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

Table 4.1: Fundamental galaxy parameters for the NFW and Einasto sets of models. Name

LB 11

Re

(10 LB, ) (2)

(kpc)

E0200

M∗

Mh

11

11

(3)

(10 M ) (4)

(10 M ) (5)

0.27

4.09

1.25

EO4200 IS

0.27

4.09

EO4200 VD EO7200 IS EO7200 VD FO4200 IS FO4200 VD FO7200 IS FO7200 VD

0.27

(1)

NFW σe8

(km s

−1

EIN σe8

)

(km s

−1

NFW fDM

EIN fDM

c

)

(6)

(7)

(8)

(9)

(10)

25

200

200

0.61

0.57

37

1.25

25

166

166

0.63

0.59

37

4.09

1.25

25

179

179

0.63

0.59

37

0.27

4.09

1.25

25

124

124

0.66

0.62

37

0.27

4.09

1.25

25

148

149

0.66

0.62

37

0.27

4.09

1.25

25

178

179

0.59

0.55

37

0.27

4.09

1.25

25

191

192

0.59

0.55

37

0.27

4.09

1.25

25

150

151

0.57

0.53

37

0.27

4.09

1.25

25

178

179

0.57

0.53

37

E0250

0.65

7.04

3.35

67

250

250

0.59

0.55

28

EO4250 IS

0.65

7.04

3.35

67

207

208

0.62

0.57

28

EO4250 VD EO7250 IS EO7250 VD FO4250 IS FO4250 VD FO7250 IS FO7250 VD

0.65

7.04

3.35

67

223

224

0.62

0.57

28

0.65

7.04

3.35

67

154

155

0.66

0.61

28

0.65

7.04

3.35

67

184

185

0.66

0.61

28

0.65

7.04

3.35

67

223

224

0.57

0.53

28

0.65

7.04

3.35

67

240

241

0.57

0.53

28

0.65

7.04

3.35

67

189

190

0.56

0.51

28

0.65

7.04

3.35

67

223

224

0.56

0.51

28

E0300

1.38

11.79

7.80

160

300

300

0.62

0.57

22

EO4300 IS

1.38

11.79

7.80

160

248

249

0.64

0.60

22

EO4300 VD EO7300 IS EO7300 VD FO4300 IS FO4300 VD FO7300 IS FO7300 VD

1.38

11.79

7.80

160

267

269

0.64

0.60

22

1.38

11.79

7.80

160

185

185

0.68

0.64

22

1.38

11.79

7.80

160

221

223

0.68

0.64

22

1.38

11.79

7.80

160

266

268

0.60

0.55

22

1.38

11.79

7.80

160

286

288

0.60

0.55

22

1.38

11.79

7.80

160

224

225

0.59

0.54

22

1.38

11.79

7.80

160

265

267

0.59

0.54

22

Notes. (1) Model name: E0 identifies the spherical progenitor, and the superscript is the value of σe8 . For the other models, the nomenclature is as follows: for example, FO4200 IS means a face-on flattened E4 galaxy, obtained from the E0200 progenitor, with isotropic rotation. (2) Luminosities in the B band. (3) Effective radius (for a FO view for FO-built models, and an EO view for EO-built models). For FO-built models, the √ edge-on effective radius is reduced by a factor q (Sect. 2.1). (4) Total stellar mass. (5) Total DM mass. (6) − (7) Stellar velocity dispersion, as the luminosity-weighted average within a circular aperture of radius Re /8, for the NFW and Einasto sets, respectively; for non-spherical models, σe8 is the edge-on viewed value. (8) − (9) DM fraction enclosed within a sphere of radius Re for the NFW and Einasto sets, respectively. (10) Concentration parameter for the NFW set.

4.2 Results

57

Figure 4.1: Meridional section of temperature (in K, top panels) and heating over cooling time ratio theat /tcool (bottom panels), for the E0250 model with NFW halo (Table 4.1), at the times specified in the boxes (in Gyr). We define tcool = E/L and theat as the ratio between E and the source terms in the r.h.s of equation (2.38). In the bottom plots, green regions refer to cooling gas, while purple indicate heating dominated regions, as indicated by the colour scale. Arrows show the meridional velocity field, with the longest arrows corresponding to 127 km s−1 .

58

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

1042

0.8 E0

E4

E7

1041

0.6 0.5

1040

TX (keV)

LX (erg s−1)

0.7

0.4 0.3 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr)

Figure 4.2: Time evolution of the X-ray luminosity LX and X-ray emission weighted temperature TX for the family derived from the E0250 model with the NFW halo. The red and black lines report the evolution of the VD models (solid), and of the IS models (dashed). The colours map the flattening: green, red and black correspond to the E0, E4 and E7 galaxies, respectively. At the same time cold gas accumulates into the centre, due to the lack of rotational support. Starting from the time of decoupling, the evolution appears to be nearly stationary. The evolution of the ISM temperature reflects the flow evolution: an hot atmosphere approximately isothermal (T ' 5 × 106 K) forms at the beginning, containing a cooling region of radius ' 5 kpc that leads to the formation of a cold core at the very centre (see the green region in the bottom panels). At the end of the simulation, a total of ' 2.6 × 1010 M of gas are cooled at the centre, while ' 5 × 109 M have been ejected as a galactic outflow. Overall, LX and TX of this model do not present significant fluctuations (Fig. 4.2, solid green line), with LX steadily decreasing and TX steadily increasing in pace with the time evolution of mass sources and specific heating (Section 2.2.1). 4.2.1.2

The EO7250 VD galaxy

The ISM evolution of the velocity dispersion supported EO7250 VD model presents important similarities with the spherical progenitor. This is not surprising, due to the absence of angular momentum, and to the fact that in general the gravitational potential is much rounder than the associated stellar density distribution (in addition, recall that the DM halo is kept spherical). Therefore, the only major differences between the E0250 progenitor and the EO7250 VD model are the different spatial regions where the gas is injected, and the different velocity dispersion field of the stars. A direct comparison of the evolution of the two models can be obtained by inspection of Fig. 4.3, analogous of Fig. 4.1. At early times the flow is kinematically decoupled, with an equatorial outflow due to the concentrated heating on the equatorial plane (purple region), associated with a polar accretion along the z-axis, evidenced by the

4.2 Results

59

Figure 4.3: Meridional sections of the temperature (top panels) and heating over cooling time ratio (bottom panels) for the EO7250 VD model of the NFW set, at the times specified in the boxes (in Gyr). Arrows are normalized to the same velocity as in Fig. 4.1.

60

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

green cooling region. As in the spherical progenitor, due to the lack of centrifugal support, the cooling material falls directly towards the galaxy centre, where a dense, cold core is formed during the first hundred of Myr. The early flow evolution is characterised by equatorial symmetry and large scale meridional vortexes. The symmetry is lost at t ' 3 Gyr, followed by a secular decrease of the velocity field. The flow velocities are larger that in the E0 progenitor, due to the weaker gravitational field (consequence of the edge-on flattening). The evolution of LX and TX are shown in Fig. 4.2 (black solid line). Compared to the E0 progenitor, the EO7250 VD model has a fainter LX and a lower TX , but a similar lack of significant fluctuations in LX and TX . At the end of the simulation, the cooled gas at the centre for this model is ' 1.8 × 1010 M , while the ejected ISM is ' 7 × 109 M . If this model is allowed to have the accretion physics activated for the central black hole, we expect to recover the complex AGN feedback phenomena described elsewhere (Ciotti & Ostriker 2012, see also Section 5.2), with significant reduction of the central accreted mass. 4.2.1.3

The EO7250 IS galaxy

Previous explorations (DC98; N14, see also Chapter 3), revealed that the evolution of gas flows in galaxies with significant ordered rotation of the stellar component is more complex than in velocity dispersion supported systems of similar structure. This is confirmed by the present study. The flow evolution of the EO7250 IS is shown in Fig. 4.4 (where more panels than the previous two models are shown, to better illustrate the more structured evolution of the ISM). The first major difference of the present model with respect to its VD counterpart is the formation, due to angular momentum conservation, of a rotationally supported, thin and dense cold disc, with a size of ' 5 kpc. The cold disc grows during galaxy evolution, reaching a final size of ' 10 kpc. A hot and rarefied zone that secularly increases in size surrounds the cold disc. At early times (t ' 2.1 Gyr), the ISM in the central regions cools and collapses, producing a low-density region that cannot be replenished by the inflowing gas, which is supported by angular momentum. As time increases, the combination of the centrifugal barrier, that keeps the centre at low density, and the secular increase of the specific heating produce the growth of the heating region (purple zone in Fig. 4.4, roughly extending as the cold thin disc). We stress that the time and spatial evolution of the theat /tcool ratio is more affected by cooling time variations than by the secular decrease of the heating time. Being the cooling time very sensitive to the ISM density, theat /tcool is strongly related to the density distribution evolution. 250 Another important difference between EO7250 IS and EO7VD concerns the ISM kinematics outside the equatorial plane. As apparent in Fig. 4.4, starting from t ' 8 Gyr the meridional velocity field develops a very complex pattern of vortexes above and below the equatorial plane. This behaviour is associated with the formation of a large cooling region (in green), and corresponds respectively to a peak and a drop in the evolution of LX and TX (Fig. 4.2, black dashed line). Note also that 250 LX and TX are the lowest of the three models E0250 , EO7250 VD and EO7IS . The cold mass accreted at the centre is now much smaller (' 2.9 × 103 M ), while the mass in the cold disc is ' 1.5 × 1010 M , and the mass ejected in the galactic wind is ' 1.1 × 1010 M . Note that a central black hole in this rotating model would produce

4.2 Results

61

Figure 4.4: Meridional sections of the temperature in K (six top panels) and heating over cooling time ratio (six bottom panels) for the EO7250 IS model of the NFW set, at the times specified in the boxes (in Gyr). Arrows are normalized to the same velocity as in Fig. 4.1.

62

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

a significantly weaker AGN activity than in the EO7250 VD model. Even if the models are structurally quite different from the S0/Sa models in N14, they show a similar ISM evolution, with a lower LX and TX in rotating models. However, the number of oscillations in the present set of models is much lower. 4.2.1.4

The thermalization parameter and an overview of all models

A useful global parameter that helps to quantify the heating of the ISM due to stellar ordered motions is the thermalization parameter γth (Section 2.3). For the EO7250 IS galaxy, and in general for medium-high mass models, we found that γth remains at low values (' 0.08 − 0.28, see Tables 4.2 and 4.3) over the evolution; this means that 1) the ISM almost co-rotates with the stellar population everywhere, and 2) there are not significant ISM velocities in the meridional plane. As a consequence, in the medium-high mass models Lkin ' Lσ (see equation 2.30). Note that γth can attain large values, even larger than unity (see the examples described later in this Section). This happens in general in low-rotation, low-mass systems, where γth is fully dominated by high velocity galactic winds (see equation (2.30)), so that Lm is large due to the thermalization of the strong meridional motions, even though Lϕ remains low. Remarkably enough, also in these high-γth cases, the azimuthal ϕ thermalization parameter γth ≡ Lϕ /Lrot (equation 2.30) remains low (see Tables 4.2 and 4.3), indicating that the ISM rotates almost as fast as the stellar population. One could be tempted to interpret the lack of thermalization of a significant fraction of ordered motions in all the IS models as the reason for the lower TX of the IS models with respect to their VD counterparts (Fig. 4.2, black and red dashed lines versus. the solid lines). However, even if this effect certainly contributes, it is not the main reason for the lower TX in rotating models. Indeed, we found that artificially adding the “missing” thermalization to the equations of hydrodynamics in dedicated test simulations of rotating models leads only to a negligible increase in TX (see also Chapter 3), showing that also other effects contribute to the low TX (see Section 4.2.3). We now discuss similarities and differences of the hydrodynamical evolution in galaxy models of different mass (i.e., derived from progenitors with different σe8 ). The main features of the family with the spherical progenitor of σe8 = 250 km s−1 are maintained in the σe8 = 300 km s−1 family. In particular, independently of the DM halo profile, increasing σe8 , both LX and TX increase. This is expected because more massive models can retain more and hotter gas independently of the flattening and kinematical support. In more massive models the LX and TX are less fluctuating with time, the outflow velocities of the galaxy outskirts are lower, and the complicated meridional circulation in the rotating models is reduced (as also found by N14). The final properties of all models are given in Tables 4.2 and 4.3. A full discussion of LX and TX is given in Sections 3.2 and 3.3. In general, the ISM temperature, luminosity, radius of the central cooling region, and inflow velocity are directly proportional to σe8 . In massive (σe8 ≥ 250 km s−1 ) models, at fixed galaxy mass, pure flattening does not affect significantly LX and TX , while a major reduction in LX and TX is obtained for the isotropic rotators. The situation is quite different for the families with low mass progenitors (σe8 = 200 km s−1 ). These are the only cases where a transition to a global wind can be

4.2 Results

63

Figure 4.5: Meridional sections of the temperature in K (six top panels) and heating over cooling time ratio (six bottom panels) for the low mass EO7200 VD model of the NFW set, at the times specified in the boxes (in Gyr). Arrows are normalized to the same velocity ad Fig. 4.1. Note the strong equatorial degassing established at late times.

64

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

LX (erg s−1)

1041

NFW

1040 1039 EO7200 VD

1038

EO7200 IS

10

37

EO4200 IS

1036

LX (erg s−1)

1041

EIN

1040 1039 1038 1037 1036 100

150 200 250 σe8 (km s−1)

300

108

109 Mhot (MO• )

1010

1011 LB (LB,O• )

Figure 4.6: ISM X-ray luminosity LX in the 0.3–8 keV band at 13 Gyr for all models in the NFW (top panels) and in the Einasto (bottom panels) sets, as a function of σe8 , of the hot (T > 106 K) ISM mass, and of the galaxy blue optical luminosity; spherical progenitors (green circles) with σe8 = (200, 250, 300) have been considered. The green, red and black colours refer to the E0, E4 and E7 models respectively. Filled and empty symbols indicate the fully velocity dispersion supported VD models, and the isotropic rotators IS models, respectively. Models in wind for the NFW set are labelled in the top left panel.

4.2 Results

65

induced by a change of shape or by rotation, in accordance with the energetic analysis of CP96. This is especially true for the less concentrated Einasto models. In these global wind cases, LX drops to very low values, due to the very low ISM density, and TX keeps larger than expected from the trend defined by non-wind models (see Section 4.2.3), due to the reduced cooling, and to thermalization of the meridional motions (see Sections 4.2.2 and 4.2.3 for a detailed discussion). The sensitivity of the flow phase for low-mass models near the transition to the outflow is shown for example by the EO7200 VD model with the NFW halo, that experiences two quite distinct evolutionary phases (Fig. 4.5). At the beginning, a significant equatorial degassing is apparent, coincident with the strong heating in that region. As time increases, the velocity field in the outflow region decreases and gas cooling becomes more and more important outside the V-shaped region around the equator. However, after ' 9 Gyr, the secular increase of the specific heating, coupled with the shallow potential well, induces again higher and higher velocities. The gas temperature increases again while LX decreases. The associated EO7200 IS model is in a permanent wind phase from the beginning, thus showing the additional effect of rotation in flattened, low-mass 200 galaxies. The differences between the EO7200 VD and EO7IS models are quantified by the associated values of the global quantities at the end of the simulation (see Table 4.2): Mhot = 0.66 × 109 M and 0.24 × 109 M in the VD and IS cases, respectively, where Mhot is the ISM mass having T > 106 K. Little accretion at the centre is present in the VD but not in the IS, and this shows how different AGN activity may be expected in rotating versus non rotating galaxies, also at low galaxy masses.

4.2.2

The X-ray ISM luminosity LX

We now move to describe the properties of LX for the whole set of galaxy models, as they would be observed at an age of 13 Gyr. The results are summarized in Fig. 4.6, where the top panels refer to the NFW set and the bottom panels to the Einasto set. LX is shown versus 3 different galaxy properties, i.e., σe8 (left panels), Mhot (central panels), and LB (right panels). Remarkably, the range of LX values spanned by the models matches the observed one (see for example the observed LX − LK and LX − σe8 trends in Figs. 2 and 5 in Boroson et al. 2011). The most interesting feature of Fig. 4.6 is the clear LX difference between flattened rotating models and models of similar σe8 but velocity dispersion supported. As described in the previous Section, the hydrodynamical simulations show that the under-luminosity of rotating galaxies with medium to large σe8 is due to a different flow evolution driven by the presence of angular momentum, which prevents the gas from accumulating in the central regions, leading to the creation of a very hot, low density atmosphere in the centre, and eventually resulting in a lower total LX . Instead, in VD models the ISM flows directly toward the central galactic regions, where a steep density profile is created. This difference in the hot gas density distribution is a major reason for the systematic difference of LX (see also Fig. 4.7). It nicely explains the lower LX observed for fast rotators than for slow rotators in the ATLAS sample (Sarzi et al. 2013). ETGs with the lowest σe8 , behave differently (see below). In the central panels of Fig. 4.6 LX is plotted against the hot gas content Mhot ; each rotating model is shifted to the left of the corresponding VD model, thus

66

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

100

ρhot (cm−3)

10−1

10−2

10−3

10−4

10−5 1

10

100

r (kpc) Figure 4.7: Angle-averaged profile of the hot ISM density at t = 13 Gyr for the same models as in Fig. 4.2. Solid lines refer to VD models, dashed lines refer to IS models. IS models have also a lower Mhot than VD models. This is due to the presence of recurrent cooling episodes driven by rotation, that further contribute to the lowering of LX . With the exception of the models with the lowest LB , the systematic differences in Mhot are not due to escaping ISM (Fig. 4.8). Finally, the right panels of Fig. 4.6 show how LX on average increases with the galaxy optical luminosity, however presenting at each LB a significant spread in LX , consistent with observations (Boroson et al. 2011). At fixed LB , round progenitors are found at high LX , while the dispersion is associated to a mix of flattening and rotation effects. At each LB , LX of the VD models is higher than that of IS ones by up to a factor of ' 40. The largest difference occurs for the more massive and flatter models, and it is much larger than the LX variation between a spherical progenitor and its most flattened VD version. Indeed, LX of VD models of identical LB with different flattening lies in a narrow range, with a weak trend for the X-ray luminosity to increase as the galaxy model gets rounder. The same behaviour occurs also among IS models with the same LB . This indicates that, at fixed LB and fixed internal kinematics, LX is only marginally sensitive to even large variations of the flattening degree of the stellar component. In Figure 4.6 it is apparent how the models with the σe8 = 200 km s−1 progenitor behave differently from the rest of the models; this is more evident for the EO flattening, when the galaxy potential well becomes shallower, and thus energetic effects of flattening and rotation are larger than for the FO flattening. For example, 200 the EO7200 VD model drops to low LX , at variance with the FO7VD model; this drop 200 happens also for the Einasto EO4VD model. This sharp LX difference is due to

4.2 Results

12

67

A. Negri, S. Posacki, S. Pellegrini & L. Ciotti NFW

1.0

EIN

Mesc /Minj

0.8 0.6 0.4 0.2 0.0 150

200 250 σe8 (km s−1)

300

150

200 250 σe8 (km s−1)

300

42

10 4.8: Fraction of escaped ISM mass (Mesc ) with respect to the total injected Figure mass (MinjNFW , see Tables 4.2 and 4.3) at t = 13 Gyr,EIN as a function of σe8 , for the whole 41 10 NFW and Einasto sets. The notation for the symbols is the same as in Fig. 4.6.

LX (erg s−1) TX (keV)

1040

1.0 1039 0.8

NFW

0.638 10 0.4 1037 0.2 36 10 0.0 1.0

0.2 EIN

0.4

0.6 TX (keV)

0.8

1.0

0.2

0.4

0.6 TX (keV)

0.8

1.0

TX (keV)

0.8 of escaped ISM mass (Mesc ) with respect to the total injected mass (Minj , see Tables A1 and A2) at t = 13 Gyr, as a function of Figure 8. Top panels: fraction σe8 , for the whole NFW and Einasto sets. The notation for the symbols is the same as in Fig. 6. Bottom panels: X-ray luminosity LX with respect to X-ray luminosity weighted temperature T X at t = 13 Gyr for the same models in the top panels. 0.6 0.4

of the flow phase to (even small) changes in the mass profile (e.g., with σe8 for the bulk of the models, a natural consequence of the flattening or mass concentration) and in the stellar kinematics (e.g., deeper potential well associated with larger σe8 . This leads to faster 0.2 rotation) is very high, and it is very difficult to predict systematic stellar (random and ordered) velocities, with the consequent larger 0.0that the VD and IS models in each panel are trends in LX . We stress energy input due to thermalization of the stellar velocity fields. In 9 100 by 150 200 250 300 108 so 10 1010 potential is more effective 1011 in retaining the hot characterized by construction the same gravitational potential, addition, a deeper −1 ) in turn can be further (km s rotation. ) (MO•that LB (L • ) by SNIa explosions. The that the difference in LX is only dueσtoe8galactic heated These cases Mhotgas, B,O point out the high sensitivity of the gas flow to the galaxy structure temperature range spanned by the models agrees well with that and kinematics at low galaxy galacticweighted winds are hosted. of real observed trend of higher T X for Figure 4.9:mass, ISMwhere emission temperature TXgalaxies; in the moreover, 0.3–8 keVtheband at 13 Gyr increasing σe8 is reproduced (e.g., see Fig. 6 in Boroson et al. 2011). for all the models in the NFW (top panels) and in the Einasto (bottom panels) sets An interesting trend can be noticed at low σe8 , marginally visible as a function of σe8 , of the hot (T > 106 K) ISM mass, andbutofclearly the galaxy bluein optical in the NFW set discernible the Einasto set, i.e. the 3.3 The X-ray emission weighted temperature luminosity. Symbols and colours are as in Fig. 4.6. of T X for decreasing σe8 , and this even more in rotating increase galaxies. This is the temperature analogous of the LX behaviour The second important diagnostic explored in this work is the lumidescribed in Fig. 6. In fact the transition to global winds in flattened nosity weighted ISM temperature in the band 0.3–8 keV. In analogy and rotating low-mass galaxies leads to a reduction in the LX but an with Fig. 6, the T X distribution of the whole set of models (both increase of T X , due to the thermalization of the resulting meridional NFW and Einasto) at the end of the simulations, as a function of flows (while the thermalization of galaxy rotation remains negligible, σe8 , Mhot and LB , is given in Fig. 10. see Sect. 3.1.4). A closer inspection of the T X −σe8 panels shows that From the left panels it is apparent how in general T X increases © 2014 RAS, MNRAS 000, 1–18

0.2 0.0 The effects shape rotation on 150 of galaxy 200 250 and 300 150the X-ray 200 haloes 250 of ETGs 300 −1 −1 68 Numerical simulations σe8 (km s ) σe8 (km s ) 1042 NFW

10

EIN

41

LX (erg s−1)

1040 1039 1038 1037 1036 0.2

0.4

0.6 TX (keV)

0.8

1.0

0.2

0.4

0.6 TX (keV)

0.8

1.0

Figure 8. Top panels: fraction of escaped ISM mass (Mesc ) with respect to the total injected mass (Minj , see Tables A1 and A2) at t = 13 Gyr, as a function of Figure 4.10: X-ray luminosity LXsymbols withisrespect toinX-ray luminosity weighted temperσe8 , for the whole NFW and Einasto sets. The notation for the the same as Fig. 6. Bottom panels: X-ray luminosity LX with respect to X-ray luminosity weighted temperature Gyr for the the same whole models inNFW the top panels. ature TX atT Xt at=t =1313Gyr for and Einasto sets. The notation for the

symbols is the same as in Fig. 4.6.

of the flow phase to (even small) changes in the mass profile (e.g., with σe8 for the bulk of the models, a natural consequence of the flattening or mass concentration) and in the stellar kinematics (e.g., deeper potential well associated with larger σe8 . This leads to faster fact flattening a flow transition to a global wind, in accordance rotation) is verythe high, andthat it is very difficult toproduces predict systematic stellar (random and ordered) velocities, with the consequent larger analysis, asindescribed in Section 4.2.1.4. Intothe NFW case, a stellar further trends in LX . Wewith stressthe that CP96 the VD and IS models each panel are energy input due thermalization of the velocity fields. In characterized byreduction constructionin byL the gravitational potential, so addition, a deeper is more effective in retaining the hot is attained when introducing rotation in thepotential EO7200 model, again in X same IS that the difference in LX is onlywith gas, that in turn be further heated does by SNIa due to CP96 galactic and rotation. These cases thermalization accordance P13, where of can ordered motions notexplosions. The point out the high sensitivity of the gas flow to the galaxy structure temperature range spanned by the models agrees well with that take place. Note how a transition to a very low LX value is also obtained for the NFW and kinematics at low galaxy mass, where galactic winds are hosted. of real galaxies; moreover, the observed trend of higher T X for EO4200 model, just by adding rotation. These findings point out the high sensitivity increasing σe8 is reproduced (e.g., see Fig. 6 in Boroson et al. 2011). IS An interesting trend (e.g., can be flattening noticed at low of the flow phase to (even small) changes in the mass profile orσmass e8 , marginally visible in(e.g., the NFW set but clearly discernible inmasses, the Einasto set, i.e. the concentration) and in the stellar kinematics rotation) at low galactic 3.3 The X-ray emission weighted temperature increase of T X for decreasing σe8 , and this even more in rotating for which then it is difficult to predict systematic trends in LX . We stress that the galaxies. This is the temperature analogous of the LX behaviour The second important diagnostic explored in this work is the lumiIS models in each family are characterized, by the same described in by Fig. construction, 6. In fact the transition to global winds in flattened nosity weightedVD ISM and temperature in the band 0.3–8 keV. In analogy gravitational potential, so that the difference in L is only due to galactic rotation. and rotating low-mass galaxies leads to a reduction in the LX but an with Fig. 6, the T X distribution of the whole set of models (both X increase of T X , due to the thermalization of the resulting meridional NFW and Einasto) at the end of the simulations, as a function of flows (while the thermalization of galaxy rotation remains negligible, σe8 , Mhot and LB , is given in Fig. 10. Sect. 3.1.4). A closerTinspection of the T X −σe8 panels shows that From the left panels itThe is apparent how inemission general T X increases 4.2.3 X-ray weightedseetemperature X © 2014 RAS, MNRAS 000, 1–18 The second important diagnostic explored is the 0.3–8 keV luminosity weighted ISM temperature TX . The distribution of the TX values for the whole set of models at the end of the simulations is given in Fig. 4.9, as a function of σe8 , Mhot and LB . In general TX increases with σe8 , a natural consequence of the deeper potential well associated with larger σe8 . This leads to faster stellar (random and ordered) velocities, with the consequent larger energy input from thermalization of the stellar motions. In addition, in a deeper potential the hot gas is retained at a larger TX . The temperature range spanned by the models agrees well with that of real ETGs, and the observed trend of TX with σe8 is reproduced (e.g., see Fig. 6 in Boroson et al. 2011, who measured TX of the pure gaseous component for a sample of 30 ETGs). At high σe8 , the observed TX values span a narrower range than in our models, likely because the models include very flat and highly rotating ETGs that

4.2 Results

69

are missing in the observed sample. Interestingly, instead, the low-σe8 end of the observed TX − σe8 relation shows an increase of dispersion in the TX values, and a hint for a flattening of the relation with respect to the trend shown at larger σe8 . These features are shown also by our models: at low σe8 the trend of TX flattens for NFW models, and the scatter around it increases considerably for the Einasto models. This is explained as the temperature counterpart of the LX behaviour at low σe8 in Fig. 4.6: the transition to global winds in flattened and rotating low-mass galaxies leads to a reduction in LX and an increase of TX with respect to the trend defined by more massive ETGs, or ETGs of similar mass but not in wind. The change in the relationship is due to the thermalization of the resulting meridional flows (while the thermalization of galaxy rotation remains negligible), and to the 200 lower cooling (see Section 4.2.1.4). For example, the EO4200 VD and EO4IS models in the Einasto set, have high TX as a consequence of the transition to the wind phase. The middle panels of Fig. 4.9 show the TX distribution as a function of Mhot . In the NFW set, there is a sequence of TX values clearly visible at Mhot > 2 × 109 M , with VD models hotter than the corresponding IS models. However, the three models with the smallest amount of hot ISM (Mhot < 109 M ) have higher temperatures than one would expect extrapolating the TX sequence to very low values of Mhot , as a consequence of the transition to the wind phase. A change in the trend is even more visible in the low mass Einasto models, where the stronger tendency to establish a global wind leads to an increase of TX at very low Mhot , reaching values even higher than in VD models with large X-ray haloes. In conclusion, at medium-high σe8 , TX of VD models tends to remain above that of rotating models; at low σe8 , in addition to the cooler branch of rotating models, another hotter branch of IS and VD models appears, made by models in wind. Finally, the right panels of Fig. 4.9 show again how TX of IS models is systematically lower with respect to that of VD ones of same LB , with the exception of those in the wind phase. As for LX , TX of VD models is dominated by the dense central luminous regions. In IS models, instead, the central region is hotter than in VD models, but it is also at a lower density, so that its contribution to TX is marginal, and TX is more affected by colder (T ' 2 × 106 K) gas located in the outer regions. Thus, the main reason of the lower TX in IS models of medium-high mass is not galaxy shape, but the importance of galaxy rotation, that drives the hydrodynamical evolution (Section 4.2.1.3). From the Jeans equations, the more a galaxy is flat, the more it can be rotating; thus the E7 IS models are cooler than their VD counterparts, and by a larger amount than for the analogous E4 pair, due to the stronger rotation in the E7 models. The trend of LX with TX for all models is shown in the lower panels of Fig. 4.10. Also in this figure the models behaviour is strikingly similar to that observed in the Boroson et al. (2011) sample, where a narrow correlation at high TX & 0.5 keV is broken into an almost vertical band of LX values spanning a large range (from 1038 to few 1041 erg s−1 ) for kT covering a small range (from 0.2 to 0.5 keV). This trend in the models is explained as the product of the effects described above, resulting in a high sensitivity of the flow phase to small variations in the galaxy structure at the lowest galaxy masses, that on average also have TX < 0.5 keV. The considerations on the thermalization parameter presented in the previous Section can be directly translated on the temperatures. Indeed, due the low values of

70

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

Figure 4.11: Edge-on 0.3–8 keV surface brightness of the ISM (ΣX ) at 13 Gyr, 250 for E0250 , EO7250 VD and EO7IS models, respectively; the brightness values on the −1 colour-bar are given in erg s cm−2 . Superimposed are the isophotes (Σ? ) obtained by projecting the galaxy stellar density distribution, with the innermost contour corresponding to 104 M pc−2 , and decreasing by a factor of ten on each subsequent contour going outwards. Note that the ΣX map of the EO7250 VD model shows a round shape and a luminous core very similar to the E0250 map, whereas EO7250 IS map presents a boxy shape and a low-luminosity core. γth , the simulations show that Tkin is almost coincident with Tσ in the medium-high σe8 models (i.e., models in a slow inflow, see Tables 4.2 and 4.3), where Tσ depend only on the galaxy structure, and do not contain contributions from gas cooling and SNIa heating (equation (2.35)). The low-σe8 wind models, instead, have Tkin > Tσ , and the temperature difference is due to thermalization of the strong meridional motions developed in the wind phase (Tkin ' Tσ + Tm , while Tϕ remains very small ϕ due to low values of γth ). Thus Tσ , except for wind cases, is a good proxy for Tkin . As a final comment, we note that, in general, at fixed σe8 , Einasto models tend to be slightly colder than the NFW models, both in TX and Tσ , due to the different dark matter profile. Relevant considerations on the global trend of Tσ with galaxy shape and kinematics are presented in Posacki (2014).

4.3

Discussion and conclusions

In this chapter, in a follow-up of a series of preliminary studies, we performed a large suite of high-resolution 2D hydrodynamical simulations, to study the effects of galaxy shape and stellar kinematics on the evolution of the X-ray emitting gaseous haloes of ETGs. Realistic galaxy models are built with a Jeans code, that allows for a full generality in the choice of axisymmetric galaxy shape and of the stellar and dark matter profiles, that can be tailored to reproduce observational constraints. The dynamical structure of the models obeys the implicit assumption of a 2-integrals phase-space distribution function. Stellar motions in the azimuthal direction are split among velocity dispersion and ordered rotation by using the Satoh (1980) decomposition. In particular, we explored two extreme kinematical configurations, the fully velocity dispersion supported system (VD) and the isotropic rotator (IS), in order to encompass all the possible behaviours occurring in nature. Of course,

4.3 Discussion and conclusions

71

the VD configuration applies only to a minor fraction of the flat galaxy population (e.g., Emsellem et al. 2011). Moreover, IS models approximate only to some extent the dynamical structure of flat and fast rotating galaxies, since the latter are more generally characterized by a varying degree of anisotropy in the meridional plane with intrinsic flattening (Cappellari et al. 2007). The source of gas is provided by secular evolution of the stellar population (stellar winds from ageing stars and SNIa ejecta). Heating terms account for SNIa events and thermalization of stellar motions. The main focus of this work is the explanation of long-standing and more recently observed trends of LX and TX with galaxy shape and rotation (as well as, of course, with fundamental galaxy properties as stellar velocity dispersion and optical luminosity). Evidences from previous exploratory theoretical (CP96; P13) and numerical works (DC98, Chapter 3) seem to point toward a cooperation of flattening and rotation in establishing the final X-ray luminosity and temperature of the ISM. However which of the two is the driving parameter, and what is the involved physical mechanism, had not been clarified yet. From the present investigation, we conclude that more than one physical effect is at play, and that the relative importance of flattening and rotation changes as a function of galaxy mass. We summarize the results discussing first the X-ray luminosity and then the emission-weighted ISM temperature. 1) In low mass galaxy models with a progenitor hosting a global wind, the effects of flattening and rotation are just to make the wind stronger, and all systems are found at the lowest values of LX . 2) In case of galaxies energetically near to the onset of a galactic wind, i.e., for ETGs with σe8 ≈ 200 km s−1 , flattening and rotation contribute significantly to induce a wind, in agreement with the energetic expectations discussed in CP96, with the consequent sharp decrease of LX . The transition to a global wind is favoured respectively by the facts that flattening can reduce the depth of the potential well, and that in rotating systems the ISM and the stellar component almost corotate; this reduces (in absolute value) the effective potential experienced by the ISM. 3) In models with σe8 > 200 km s−1 , galaxy shape variations, in absence of rotation, have only a minor impact on the values of LX , in the sense that fully velocity dispersion supported flattened models have LX similar to or just lower than that of their spherical progenitors. 4) In flat galaxies with σe8 > 200 km s−1 , rotation reduces significantly LX . Not only the thermalization parameter is low and part of the heating due to stellar motions is missing with respect to the corresponding VD models, but rotation acts also on the hydrodynamics of the gas flow: conservation of angular momentum of the ISM injected at large radii favours gas cooling through the formation of rotating discs of cold gas, reducing the amount of hot gas in the central regions and then LX . The effects of angular momentum are clearly visible in Figure 4.11, where we show the edge-on projected X-ray surface brightness maps. In conclusion, galaxy flattening has an important, though indirect effect for medium-to-high mass galaxies, in the sense that only flattened systems can host significant rotation of the stellar component. 5) The luminosity evolution and the luminosity values at the end of the simulations are similar for the NFW or Einasto dark matter haloes (at fixed stellar structure

72

The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

and similar values of the dark matter halo mass). The main results concerning the emission-weighted temperature TX can be summarized as follows: 6) As for LX , also for TX the response to a variation of shape and internal kinematics is different for low and high mass galaxies. TX does not change appreciably adding flattening and rotation to low mass progenitors that are in the global wind phase. Due to their low density and high meridional velocities, global winds are generally hotter than what expected by extrapolation of the TX of more massive systems. As described at point 2) above, adding flattening and rotation to ETGs energetically near to host a global wind leads to a transition to a wind phase, with the consequent increase of TX . 7) In the medium-high mass galaxies a change of shape produces small changes in TX . Adding rotation, instead, results in a much lower TX . This is because angular momentum conservation leads to the formation of a massive centrifugally supported cold disc and to a lower density of the hot ISM in the central regions above and below the equatorial plane, with respect to VD models. Then, the external colder regions weight more in the computation of TX . 8) Overall, for medium-high mass galaxies, TX increases with galaxy mass, independently of the specific dark matter halo profile. In general, in the Einasto haloes the hot gas is systematically cooler and with a larger scatter in TX , than in the NFW dark matter haloes of comparable mass. 9) In rotating models the ISM almost corotates with the stars, and so there is a corresponding reduction of the thermalization of the galaxy streaming velocity. At the same time the rotating ISM is less bound, due to the centrifugal support. With the exception of low mass galaxies in the wind phase, Tσ (the temperature associated with the thermalization of the stellar velocity dispersion) is a good proxy for Tkin , the true thermalization temperature of stellar motions, as computed from the simulations; for wind models instead Tkin > Tσ . A few important physical phenomena are still missing from the simulations. First, it is obvious that in rotationally supported models the massive and rotating cold discs are natural places for star formation. For observational studies it would be interesting to estimate age and mass of the new stars. From the point of view of the present investigation, the formation of stars, by reducing the amount of cold gas in the equatorial plane, could in principle also modify the evolution of the ISM. A second aspect missing here is the self-gravity of the gaseous cold disc. It is expected that self-gravity acts not only to promote star formation, but also to develop non axisymmetric instabilities, that lead to non-conservation of angular momentum of the gas. Phenomenologically, the effects of self-gravity can be viewed as a “gravitational viscosity” (e.g. Bertin & Lodato 2001), that favours accretion of cold gas toward the centre. Such a gas flow toward the centre is of great importance for feedback effects from a central massive black hole in rotating galaxies (Novak et al. 2011; Gan et al. 2014). In the next chapter we present our ongoing work of star formation in ETGs, by simulating a sub-group of the NFW set of galaxies shown in this chapter, with star formation activated. Moreover, we present also an improved numerical code, result of a joint collaboration with G.S. Novak, accounting for BH accretion and state-of-the-art AGN feedback in realistic ETGs, along with the the first simulations

4.3 Discussion and conclusions

in velocity supported galaxies.

73

74

Table 4.2: Simulations results for the NFW set at t = 13 Gyr. Minj

Mesc

Mgas

Mhot

LX

TX

LSN

Tkin

(109 M ) (109 M ) (109 M ) (109 M ) (1040 erg s−1 ) (keV) (1040 erg s−1 ) (keV)



Tv

Tm

(keV)

(keV)

(keV)

γth

ϕ γth

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

E0200

12.0

3.2

8.9

2.91

2.06

0.48

10.2

0.26

0.26

1.7E-3

1.7E-3





EO4200 IS EO4200 VD EO7200 IS EO7200 VD FO4200 IS FO4200 VD FO7200 IS FO7200 VD

11.9

9.8

2.4

0.15

8.29E-4

0.50

10.2

0.40

0.17

0.24

0.21

2.62

0.25

11.9

4.5

7.5

2.84

1.99

0.46

10.2

0.26

0.26

1.5E-3

1.5E-3





11.9

10.3

1.8

0.24

1.95E-3

0.49

10.2

0.29

0.08

0.22

0.18

1.24

0.22

11.9

9.4

2.8

0.66

1.89E-2

0.50

10.2

0.33

0.25

7.8E-2

7.8E-2





12.0

3.3

8.8

1.73

0.45

0.37

10.3

0.21

0.18

2.9E-2

9.0E-3

0.31

0.21

12.0

2.9

9.2

2.46

2.21

0.48

10.3

0.27

0.27

3.7E-3

3.7E-3





12.0

3.4

8.7

1.36

0.25

0.32

10.2

0.12

0.09

2.9E-2

1.1E-2

0.15

0.10

12.0

4.0

8.1

1.51

1.82

0.49

10.2

0.28

0.28

2.0E-3

2.0E-3





E0250

32.2

5.1

27.5

6.43

11.1

0.69

24.3

0.42

0.42

1.2E-3

1.2E-3





EO4250 IS EO4250 VD EO7250 IS EO7250 VD FO4250 IS FO4250 VD FO7250 IS FO7250 VD

31.7

8.3

23.8

4.02

0.76

0.55

23.9

0.29

0.27

1.8E-2

3.1E-3

0.13

0.10

31.7

6.9

25.2

6.17

9.50

0.67

23.9

0.42

0.41

1.5E-3

1.5E-3





30.6

11.4

19.7

3.42

0.33

0.55

23.1

0.18

0.13

5.5E-2

4.9E-3

0.20

0.18

30.6

12.5

18.7

3.83

4.87

0.62

23.1

0.41

0.40

1.8E-3

1.8E-3





32.2

6.5

26.2

3.80

0.87

0.56

24.3

0.30

0.29

1.9E-2

3.9E-3

0.13

0.10

32.2

5.2

27.4

5.62

11.1

0.72

24.3

0.43

0.43

1.6E-3

1.6E-3





32.2

6.4

26.0

2.91

0.43

0.50

24.2

0.19

0.15

3.7E-2

7.9E-3

0.13

0.10

32.2

7.4

25.3

3.82

10.3

0.72

24.2

0.45

0.45

1.7E-3

1.7E-3





E0300

71.3

7.6

64.7

14.70

43.3

0.94

49.1

0.65

0.65

1.0E-3

1.0E-3





The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

name

name

Minj

Mesc

Mgas

Mhot

LX

TX

LSN

Tkin

(109 M ) (109 M ) (109 M ) (109 M ) (1040 erg s−1 ) (keV) (1040 erg s−1 ) (keV)



Tv

Tm

γth

ϕ γth

(keV)

(keV)

(keV)

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

EO4300 IS

69.4

14.1

56.0

6.52

1.34

0.68

47.7

0.43

0.42

1.5E-2

1.9E-3

0.07

0.06

EO4300 VD EO7300 IS EO7300 VD FO4300 IS FO4300 VD FO7300 IS FO7300 VD

69.4

10.5

59.9

13.61

37.7

0.90

47.7

0.63

0.63

1.2E-3

1.2E-3





65.5

15.9

50.2

6.13

1.04

0.56

45.0

0.23

0.20

2.3E-2

1.5E-3

0.06

0.05

65.5

17.8

48.8

9.92

25.0

0.83

45.0

0.62

0.62

1.2E-3

1.2E-3





71.8

11.9

60.9

6.11

1.35

0.68

49.4

0.45

0.43

1.7E-2

4.0E-3

0.08

0.06

71.8

8.5

64.3

12.91

41.7

0.99

49.4

0.65

0.65

1.5E-3

1.5E-3





71.9

11.1

61.8

5.07

0.90

0.59

49.5

0.25

0.23

2.7E-2

4.0E-3

0.06

0.05

71.9

12.7

60.4

9.19

36.7

1.02

49.5

0.67

0.67

1.6E-3

1.6E-3





4.3 Discussion and conclusions

Table 4.2– continued

Notes. (1) Name of the model. (2) − (3) Total ISM mass injected into and escaped from the numerical grid, respectively. Differences in Minj for models of same LB are accounted for different sampling of ρ∗ over the numerical grid. (4) Total ISM mass retained within the galaxy at the end of the simulation. (5) − (7) ISM mass with T > 106 K, ISM X-ray luminosity in the 0.3–8 keV band, and ISM X-ray emission weighted temperature in the same band, at the end of the simulation. (8) SNIa heating rate at the end of the simulation. (9) − (12) Thermalization temperatures of stellar motions at the end of the ϕ simulation, defined accordingly to equations (2.35)-(2.34). By construction, Tkin = Tσ + Tv ; for rotating models Tv = γth Trot and Tϕ = Tv − Tm = γth Trot , while for velocity dispersion supported models Tv = Tm (see Chapter 1 and Section 2.4.3). (13) − (14) Thermalization parameter, and its azimuthal ϕ component γth = Lϕ /Lrot (see equation (2.30)), at the end of the simulation.

75

76

Table 4.3: Simulations results for the Einasto set at t = 13 Gyr. Minj

Mesc

Mgas

Mhot

LX

TX

LSN

Tkin

(109 M ) (109 M ) (109 M ) (109 M ) (1040 erg s−1 ) (keV) (1040 erg s−1 ) (keV)



Tv

Tm

(keV)

(keV)

(keV)

γth

ϕ γth

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

E0200

12.0

3.9

8.2

2.67

1.89

0.41

10.2

0.24

0.24

3.7E-3

3.7E-3





EO4200 IS EO4200 VD EO7200 IS EO7200 VD FO4200 IS FO4200 VD FO7200 IS FO7200 VD

11.9

12.1

0.1

0.09

5.25E-4

0.74

10.2

0.49

0.15

0.34

0.31

4.34

0.28

11.9

12.1

0.1

0.09

5.51E-4

0.73

10.2

0.54

0.23

0.31

0.31





11.9

12.0

0.2

0.19

1.25E-3

0.45

10.2

0.34

0.07

0.27

0.24

1.79

0.14

11.9

11.3

0.9

0.20

1.56E-3

0.45

10.2

0.42

0.22

0.20

0.20





12.0

6.0

5.7

1.52

0.97

0.21

10.3

0.21

0.16

5.2E-2

2.1E-2

0.63

0.37

12.0

3.6

8.5

2.31

1.94

0.45

10.3

0.25

0.25

2.1E-3

2.1E-3





12.0

3.9

8.2

1.69

0.16

0.25

10.2

0.14

0.09

4.8E-2

3.0E-2

0.28

0.11

12.0

4.9

7.3

1.60

2.10

0.44

10.2

0.26

0.26

2.8E-3

2.8E-3





E0250

32.2

6.6

26.0

6.47

10.1

0.63

24.3

0.37

0.37

1.4E-3

1.4E-3





EO4250 IS EO4250 VD EO7250 IS EO7250 VD FO4250 IS FO4250 VD FO7250 IS FO7250 VD

31.7

10.3

21.8

4.03

1.01

0.49

23.9

0.26

0.24

1.8E-2

3.4E-3

0.15

0.12

31.7

9.2

22.9

5.81

7.90

0.59

23.9

0.36

0.36

1.9E-3

1.9E-3





30.6

14.2

16.9

2.77

0.25

0.48

23.1

0.18

0.11

6.2E-2

9.4E-3

0.27

0.26

30.6

16.0

15.2

3.29

2.86

0.57

23.1

0.35

0.34

3.2E-3

3.2E-3





32.2

8.1

24.5

3.93

1.46

0.43

24.3

0.27

0.26

1.7E-2

3.8E-3

0.13

0.10

32.2

6.9

25.8

5.38

9.79

0.66

24.3

0.38

0.38

2.0E-3

2.0E-3





32.2

8.0

24.4

4.99

4.64

0.26

24.2

0.17

0.14

2.9E-2

1.1E-2

0.11

0.07

32.2

9.5

23.3

3.45

7.95

0.65

24.2

0.40

0.40

2.0E-3

2.0E-3





E0300

71.3

9.6

62.7

13.97

39.9

0.85

49.1

0.56

0.56

1.0E-3

1.0E-3





The effects of galaxy shape and rotation on the X-ray haloes of ETGs Numerical simulations

name

name

Minj

Mesc

Mgas

Mhot

LX

TX

LSN

Tkin

(109 M ) (109 M ) (109 M ) (109 M ) (1040 erg s−1 ) (keV) (1040 erg s−1 ) (keV)



Tv

Tm

γth

ϕ γth

(keV)

(keV)

(keV)

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

EO4300 IS

69.4

17.2

53.2

7.10

1.57

0.63

47.7

0.38

0.36

1.5E-2

1.6E-3

0.09

0.08

EO4300 VD EO7300 IS EO7300 VD FO4300 IS FO4300 VD FO7300 IS FO7300 VD

69.4

13.0

57.5

12.71

33.2

0.81

47.7

0.54

0.54

1.4E-3

1.4E-3





65.5

20.3

45.6

6.47

1.27

0.51

45.0

0.20

0.18

2.6E-2

2.2E-3

0.08

0.07

65.5

23.5

43.2

7.58

18.3

0.73

45.0

0.52

0.52

1.8E-3

1.8E-3





71.8

14.5

58.0

6.50

1.45

0.63

49.4

0.40

0.38

1.8E-2

3.6E-3

0.10

0.08

71.8

11.1

61.8

11.70

37.8

0.90

49.4

0.57

0.57

1.7E-3

1.7E-3





71.9

14.1

56.6

5.79

1.35

0.48

49.5

0.23

0.21

2.5E-2

3.5E-3

0.07

0.06

71.9

16.9

56.3

7.43

31.1

0.95

49.5

0.60

0.59

1.8E-3

1.8E-3





4.3 Discussion and conclusions

Table 4.3– continued

Notes. All quantities are as in Table 4.2.

77

Chapter 5

Ongoing and future developments In the previous studies on the effects of shape and rotation on the X-ray emitting halo of ETGs, we found that the ISM angular momentum plays a crucial role in determining the ETGs X-ray properties. A direct consequence of angular momentum conservation is the formation, in rotationally supported galaxies, of a massive, dense cold disc. These discs are the natural place for star formation that, in principle, can alter the ISM evolution, by removing cold gas in the equatorial plane. In addition, when a new stellar population is formed, SNII events take place on a short time-scale, injecting mass and energy via core-collapse stellar explosions. Thus, the X-ray features behaviour of the rotating models, studied in the previous Chapter, may be affected by the missing physics of star formation. Moreover, an interesting link between our finding of an ubiquitous formation of a cold disc in rotating systems and observed galaxy properties is given by the fact that, among ETGs, it is only in fast rotators that some degree of star formation is observed (Davis et al. 2011; Young et al. 2011; Sarzi et al. 2013; Davis et al. 2014). Motivated by the above considerations and observational results, we re-simulated a sub-group of the previous investigation NFW set of models, with the star formation physics activated, in order to study its effects on the ETGs hot haloes, and test the robustness of our previous results. The additional important physical phenomena currently missing in our simulations are the cold disc self-gravity and the AGN feedback, the latter widely studied in 1D and 2D numerical simulations with simple spherical galaxy models. Given the ubiquitous presence of BH at the centre of massive galaxies, we started a joint collaboration with Dr. Novak, with main goal to simulate the evolution of realistic, flat and rotating state-of-the-art galaxies in presence of AGN feedback, by including our galaxy models and the self-consistent feedback from the passive stellar population in the Novak et al. (2012) code. We present here the implementation details and test runs, along with the preliminary results of AGN feedback in both round and flattened E4 galaxies, supported either by rotation or tangential velocity-dispersion. As expected, the cold disc formation due to the centrifugal barrier inhibits the mass accretion onto the central BH, since the disc self-gravity is not considered. In order to mimic the angular momentum transport due to the three dimensional gravitational instabilities, we present here a numerical scheme based on the fluid viscosity formalism.

80

5.1 5.1.1

Ongoing and future developments

Star formation in ETGs Star formation input physics and galaxy models selection

In the present work we employed two different schemes for star formation (SF), in order to isolate and better study the effects of the cold disc in the previous simulations of rotating models. In the first scheme, only cold gas removal from the numerical mesh is allowed, following the Kennicutt (1998) recipe, see equations (2.19)-(2.21). In the second case, also the later injection of mass, momentum and energy from the new stellar population and SNIIe explosions is considered. In other words the former, labelled as “passive” SF, allows only for the terms ρ˙ SF and E˙ SF in equations (2.36)(2.38) in addition to the passive stellar population terms, while the latter, labelled “active” SF, solves the complete form of equations (2.36)-(2.38). This approach allows to investigate how much the massive cold discs of the previous study trigger (or reduce) the cooling episodes (CEs) observed in all rotating models. In order to explore the parameter space, each galaxy model is simulated both with passive and active SF schemes, where we adopted two different values of the star formation efficiency ηSF , respectively equal to 10−1 and 10−2 ; thus every model is re-simulated four times. Since a complete re-simulation of the two sets of galaxies analysed in the previous study would require an unaffordable cost in terms of computational time, we selected a sub-group of the NFW set of models, presented in Section 4.1. In order to explore the effects of SF on both the high and low galaxy mass models of the previous investigation, we focussed on the rotating models with the most massive and extended cold discs (the σe8 = 300 km s−1 family), and on the galaxies near to the transition from an inflow to a galactic wind (the σe8 = 200 km s−1 family). In addition, we simulated also the VD models of the σe8 = 200 km s−1 family so as to verify how much the removing of the cold gas collapsed in the central grid modifies the evolution of the X-ray halo.

5.1.2

Results

We present here the results of our simulations by focussing on LX , TX and the hydrodynamical evolution of a few representative models. All the models discussed here are gathered in groups of three, consisting of the “parent” model (i.e. without star formation), whose simulation results are presented in Table 4.2, and its passive and active SF counterparts, both with ηSN = 10−1 . For each group we describe the hydrodynamical evolution of each model, in order to describe the effects of SF on the hot gaseous halo. All the global quantities of the simulations at 13 Gyr are reported in Tables 5.2-5.5. As already found in the previous Chapters, at the beginning all the hydrodynamical quantities are nearly symmetric with respect the equatorial plane. During the evolution, the equatorial symmetry is lost, at a time that depends on the specific model, and it is never restored. In addition, the results of the previous study are confirmed, especially the low LX of the hot haloes in IS galaxies with respect the VD models at the same mass. Finally, as expected, the LX and TX time evolution of a given star forming model is more and more different with respect to its parent model at increasing ηSF .

5.1 Star formation in ETGs

5.1.2.1

81

The EO7200 IS models

The hydrodynamical evolution of the EO7200 IS models is presented in Figure 5.1. The early configuration of the parent model reflects the main characteristics of the σe8 = 250 km s−1 IS models presented in Section 4.2.1.3: a cold, rotationally supported disc is formed by the cooling of the hot gas injected in-situ, an extended cooling V -shaped region and an external equatorial degassing (green and purple regions in Figure 5.1) are also present. However, due to the combined effects of the secular increase of the specific heating and the shallower potential well of the EO7200 IS model, the very central purple region is more and more heated at time increasing, as a consequence it expands until it breaks the cooling V-shaped region at 3.5 Gyr (top second panel in Figure 5.1). After this moment, a galactic wind is established and it lasts until the end of the simulation (last panel of Figure 5.1), with the resulting decrease in LX and increase in TX , as shown in Figure 5.2 (see Sections 4.2.2 and 4.2.3 for a complete discussion of the effects of the galactic wind on LX and TX , and Section 3.2.1.2 for a complete hydro description of a cooling cycle). From Section 4.3, we remind that the σe8 = 200 km s−1 models lie on the onset of a galactic wind, i.e. a little variation of the flattening and/or rotational support may produce a significant difference in the flow evolution (e.g. the VD counterpart of the present model shown in Figure 4.5 switches to a global wind at 9 Gyr). At 13 Gyr the model presents 1.6 × 109 M of cold mass (i.e. gas having T < 2 × 104 K) confined in the disc (see Table 5.3). The introduction of passive star formation (i.e. of a gas sink in cold, dense regions) modifies the hydrodynamical evolution of the ISM, as it can be seen in the middle row of Figure 5.1. At 3 Gyr, the ISM configuration is extremely similar to that of the parent model, except for a reduced dimension of the cold disc from 2.5 to 1.3 kpc due to the SF sink. Nevertheless, instead of an early set up of a galactic wind, the V -shaped region strengthens, the radiative cooling becomes catastrophic and cold gas accretes onto the disc (dark green region in the middle second panel), prompting the LX peak and TX drop at ' 3.8 Gyr (Figure 5.2). The CE produces a central hot low density atmosphere where the radiative cooling is inefficient. This region is kept at low density by the centrifugal barrier, thus the central purple region is heated and it extends until a galactic wind is finally established (third panel, middle row of Figure 5.1), as in the parent model. The final configuration shows no cold disc, being the entire cold mass of the parent model converted into stars ∗ (Mnew ' 1.7 × 109 M ), while ' 1010 M of ISM has escaped due to the galactic winds. The bottom row of Figure 5.1 shows the hydrodynamical evolution of the active SF run. As in the passive SF case, at early times the ISM configuration is similar to that of the parent model, except for an additional heating due to the formation of a new stellar population in the cold disc (t = 3 Gyr in Figure 5.1). At variance of the passive SF run, where the SF sink prompts a major CE at 3.8 Gyr, the additional mass and energy returns due to SNIIe explosions produce a low velocity equatorial degassing. However, the energy injection rate is not sufficient to establish a wind, thus, the ISM density increases due to the mass injected from the initial stellar population (1 M yr−1 at 6 Gyr, while SNII injects ' 2 × 10−2 M yr−1 ), until the tcool becomes shorter than theat . At 8.2 Gyr a CE occurs, as evidenced by the dark

Figure 5.1: Meridional sections of the heating over cooling time ratio for the EO7200 IS models. From top to bottom, each row refers to the −1 run without SF, passive and active SF simulations (ηSF = 10 ), respectively.

82 Ongoing and future developments

5.1 Star formation in ETGs

83

0.8 EO7

200 IS

41

10

no SF A SF P SF

0.7

0.5

TX (keV)

LX (erg s−1)

0.6

0.4

1040

0.3 0.2 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr) 0.8 FO7

200 IS

41

10

no SF A SF P SF

0.7

0.5

TX (keV)

LX (erg s−1)

0.6

0.4

1040

0.3 0.2 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr) 0.8 FO7 41

10

300 IS

no SF A SF P SF

0.7

0.5 1040

TX (keV)

LX (erg s−1)

0.6

0.4 0.3 0.2 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr)

Figure 5.2: Time evolution of the X-ray luminosity LX and X-ray emission weighted 200 300 temperature TX for the EO7200 IS (top panels), FO7IS (middle panels), FO7IS (bottom panels) models. The blue, orange and magenta lines refer to the parent model (without SF), passive and active star formation simulations, respectively, with ηSF = 10−1 . 300 The EO7300 IS model is not shown since it is almost identical to FO7IS model.

84

Ongoing and future developments

green region in Figure 5.1, expanding the gaseous disc in size and mass. After this moment, following the same dynamics of the passive SF model, a global wind takes place. Thus, due to this last CE, the active SF run has a larger cold disc (R ' 5 kpc and Mc ' 5.7 × 107 M ) with respect to the passive SF simulation, but again the cold disc present in the parent run has been consumed by SF processes, and a stellar disc of 3.8 × 109 M has been produced (visible in the top left panel of Figure 5.5), which is more massive than the parent model gaseous cold disc thanks to the low total mass escaped via galactic winds in the active SF simulation. The different evolution of LX and TX due to the introduction of SF can be appreciated from Figure 5.2. It is clear that the introduction of the cold gas removal induces an additional CE at early times (cf. blue and orange lines). Interestingly, the active SF run does not behave as one may expect on the ground of simple energetic considerations: the additional heating due to SNIIe explosions avoids an early CE, thus allows the refilling of the galaxy halo mainly with mass losses of the old stars, after which a major cooling occurs at later times. 5.1.2.2

The FO7200 IS models

We now consider the FO7200 IS models, whose hydrodynamical evolution is presented in Figure 5.3. The hot gas flow evolution is much more complicated with respect to the previous cases, due to recurring CEs in every models. This is not surprising, since the FO flattening produces a deeper potential well than the EO one (Section 4.1). Thus, the increased capacity of retaining the hot gas produces in the parent model an inflow initial configuration that evolves in three major cooling events at t = 2.8, 3.8 and 10.2 Gyr (the last two shown in the first and fourth panels of Figure 5.3, top row), unlike the EO7200 IS parent model, with the associated oscillations in LX and TX (Figure 5.2); the hydrodynamical features of every CE is the same as the ones occurring in the precedent models. Additional differences include the presence of a larger cold disc (Mc ' 7.2 × 109 M and radius 8.5 kpc at 13 Gyr, due higher stellar streaming velocity than the EO counterpart), and the final kinematical configuration of the hot gas, which is in a partial inflow state. In the passive SF run, the hydrodynamical evolution changes substantially. In the initial stages, the flow kinematical configuration is almost the same as for the parent model: they both experiment a major CE at t = 3.8 Gyr (dark green regions in both first and second row of Figure 5.3). However, the reduction of cold gas due to the passive SF reduces the pressure support in the central regions, and thus another CE develops at t = 5.8 Gyr (second panel). After this moment, the halo is emptied of hot ISM to a point such that a global wind is established at 7.5 Gyr and it lasts until 13 Gyr, with the consequent reduction in LX (orange line in Figure 5.2), while the parent model exhibits a long refilling phase of the hot gas (cf. the last three panels of Figure 5.3). At 13 Gyr the disc, which is not subjected to the galaxy wind but consumed by the passive SF, has a radius of ' 2.5 kpc and contains ' 2 × 107 M of cold gas, while 5.8 × 109 M of hot gas has escaped and 5.8 × 109 M has been converted into stars. The presence of the SNIIe explosions further modifies this picture, as displayed by the last row in Figure 5.3. Indeed, the active SF run does not suffer major CEs as in the previous models, due to the additional injection of energy. After the formation

Figure 5.3: Meridional sections of the heating over cooling time ratio for the FO7200 IS models. From top to bottom, each row refers to the −1 run without SF, passive and active SF simulations (ηSF = 10 ), respectively.

5.1 Star formation in ETGs 85

86

Ongoing and future developments

of the cold disc at early times by the cooling of the in-situ ISM (a feature common to all the rotating models), a new stellar disc forms and SNIIe began to explode, injecting additional energy with respect to the passive SF run. As a consequence, a global cooling of the hot ISM is prevented and a decoupled flow is then established, as shown in the last three panels of Figure 5.3, with the purple region extending due to the time increasing specific heating (Section 2.2.1). However, a galactic wind does not take place, instead until the end of the simulation the halo is filled by the stellar mass losses, presenting a reduced heating central zone and an increased cooling region. Presumably, if the run had continued, the ISM would have irradiated its internal energy and then fallen towards the centre to form a new cold disc. Thus, as in the EO7200 IS models, the active SF delays the CE but in this case it does not prompt a wind. This is reflected in the evolution of LX and TX in Figure 5.2, where they do not present strong oscillations, instead LX steadily increases during the last 4 Gyr. The different evolution of the cold gaseous disc in this group of simulations deserves some considerations. Indeed, a major difference between the parent model and the active SF one is the outstanding reduction of the cold disc, from 8.5 to ' 0.5 kpc in radius and from ' 7.2 × 109 M to ' 2.6 × 106 M in mass (see also Figure 5.5). This is due both to the SF processes, which completely consume the cold disc formed at early times, and to the fact that no major CEs develop thanks to the additional heating of SNIIe. On the contrary, in the parent simulation a CE take place at 10.2 Gyr, thus doubling the cold disc size (cf. the fourth and fifth panels of Figure 5.1, first row). From Figure 5.2 it is clear that the enhanced cooling of the passive SF produces a strong deviation from the parent model history, while the active SF, although without ample oscillations and with a higher final luminosity, produces a LX and TX evolution less divergent from that of the parent one. 5.1.2.3

The EO7300 IS models

Finally, the hydrodynamical evolution of the EO7300 IS models are presented in Figure 5.4; note the different scale of each panel, in order to show the large cold disc in their entirety. As it is clear from the figure, the ISM behaviour is fairly independent from the adopted scheme for the SF. After the initial formation of the cold disc, the green cooling region expands and moves outwards, as more and more ISM cools onto the disc, while the central heating region (in purple) grows in size and intensity. The entire process is almost stationary, without the CEs that characterise the lighter models, and the mass loss via galactic wind is nearly independent of SF. The parent model disc is massive (' 4.5 × 1010 M ) and it extends until a radius of 28 kpc. The passive and active SF models behave in the same way, although the intense star formation in the disc reduces its size and mass to 25 kpc and ' 1.3 × 109 M in the case of passive SF, and to 23 kpc and ' 1.5 × 109 M for the active SF run (note the deep purple regions at R = 12 kpc in the last bottom panel of Figure 5.4, a clear sign of a new stellar population). The bottom panels of Figure 5.2 show an evolution of LX and TX almost independent of SF, where the blue, orange and magenta lines nearly overlap each other. Indeed, due to both the strong gravitational potential well and the large galaxy mass,

5.1 Star formation in ETGs

87

LX is completely dominated by the secular decrease of the passive stellar population heating term, without presenting any sign of CEs. 5.1.2.4

Star formation rates

In analogy with Figure 5.2, in Figure 5.6 we present the evolution of both the SFR 200 300 and the hot mass content for the EO7200 IS , FO7IS and EO7IS models. It is evident the different evolution of the SFR at low and high galaxy masses. The EO7200 IS and 200 FO7IS models present recurrent increases of the SFR (with maximum values of 1.34 and 1.77 M yr−1 ), followed by a secular decline due to the cold gas depletion. The peaks are correlated with drops in the hot gas mass, thanks to the CEs that, in this case, fuel the SF. Instead, the EO7300 IS model shows a smooth declining evolution of both the SFR and Mhot . Therefore, in high-mass models SF is sustained by the large cold gas reservoir in the disc, formed at early times, while in low-mass galaxies it is almost consumed on a time-scale of few Gyr (see the discussion of the EO7200 IS model), and then the SF process is fuelled mainly by the cooling of the X-ray emitting halo. 5.1.2.5

A review of all the models

We discuss here similarities and differences between the effects of star formation on the hydrodynamics, LX and TX for the entire set of simulations. As a general trend, all the models evolution are more and more similar to the respective parent galaxy at decreasing ηSF , due to the reduced capacity of transforming the cold, dense ISM into stars. Moreover, the new stellar population mean formation time hti∗ (equation 2.47) for the entire set of models lies in the range ' 2.8–6.6 Gyr (Tables 5.1 and 5.5), thus the global stellar component (passive stellar population plus stellar disc) at t = 13 Gyr is old, being the stellar disc dominated by giant stars. The main features of the EO7300 IS model are maintained in the hydrodynamical evolution of the entire σe8 = 300 km s−1 IS family, whose models present a smooth evolution and form at early times a cold, massive gaseous disc that is consumed by the SF. In particular, for high star formation efficiencies the cold disc is almost entirely consumed at 13 Gyr (only ' 109 M of cold gas remains), whereas a gaseous disc of ' 1010 M is still present for low ηSF (see Tables 5.1 and 5.2). The today rate of SF is ' 1.7–3 M yr−1 , well above the mean SFR of 0.14 M yr−1 found in Davis et al. (2014) for the ATLAS3D sample of ETGs; however we point out that this family poorly represents real ETGs due to the high stellar and DM mass content. The time evolution of LX and TX is not modified by star formation processes, although a large stellar disc is formed with a mass ranging from 50 to 90 percent of the total mass losses of the initial stellar population. The situation is different in the low galaxy mass family having σe8 = 200 km s−1 . The common trend for all the models is the lack of the cold central disc (or core for non-rotating models) at 13 Gyr, being consumed by SF during the entire time evolution. In addition, the final SFR ranges from 0.01 to 0.69 M yr−1 (with a mean of 0.14 M yr−1 for the IS models, see Tables 5.3-5.6), well in accordance with similar galaxies in the Davis et al. (2014) sample. The VD models behave as their respective parent models: the almost radial accretion allowed by the lack of angular momentum forms a cold dense core at early

88

Ongoing and future developments

Figure 5.4: Meridional sections of the heating over cooling time ratio for the EO7300 IS models. From top to bottom, each row refers to the run without SF, passive and active SF simulations (ηSF = 10−1 ), respectively. Note the different scale adopted with respect to Figures 5.1 and 5.3, in order to show the large cold discs in their entirety.

5.1 Star formation in ETGs

89

Figure 5.5: Meridional sections of the time integrated density distribution of the new stars formed up to 13 Gyr (in units of M pc−3 ), for the EO7200 IS (left panels), 300 models (right panels). Top and bottom panels FO7200 (middle panels), and EO7 IS IS refers to passive and active SF simulations (ηSF = 10−1 ), respectively. Note the different scale adopted with respect to Figures 5.1 and 5.3, in order to show the large cold discs in their entirety. times, which is entirely converted into stars1 , with a final SFR mean and maximum value of 0.27 and 0.45 M yr−1 . Despite the early intense SF at the centre of VD galaxies, LX and TX evolutions do not show any variation in presence of star formation. On contrary, the evolution of the IS counterparts can be substantially altered. This is not surprising, since the rotating models in this range of galaxy masses are known to be extremely sensible to variation of shape and rotation (see Section 4.3). However, although the hydrodynamical evolution of the IS models can be modified, in general the final galactic wind configuration of the parent EO-built IS models2 is maintained, while for the more gravitationally bound FO-built models the final status of the hot corona is altered, switching from inflow to a wind configuration. The effects of SF on the final value of the X-ray luminosity are even clearer from Figure 5.7, where we report LX at 13 Gyr for the parent models and for all the simulations with SF. In the most massive galaxies, along with the lighter VD 1

Note that if the black hole had been activated in these simulations, only a part of the cold mass would have been converted into stars, and a part accreted onto the BH. 2 At exception of the EO4200 IS galaxies in presence of passive SF, in which a global wind is establishing at 13 Gyr.

90

Ongoing and future developments

EO7200 IS

1010

100 10-1 109

Mhot (MO•)

SFR (MO• yr-1)

101

10-2 10-3 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr) FO7200 IS

1010

100 10-1 109

Mhot (MO•)

SFR (MO• yr-1)

101

10-2 10-3 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr) EO7300 IS

1010

100 10-1 109

Mhot (MO•)

SFR (MO• yr-1)

101

10-2 10-3 3 4 5 6 7 8 9 10 11 12 13 3 4 5 6 7 8 9 10 11 12 13 t (Gyr) t (Gyr)

Figure 5.6: Time evolution of the SFR and hot mass contained in the numerical grid 200 300 for the EO7200 IS (top panels), FO7IS (middle panels), FO7IS (bottom panels) models. The blue, orange and magenta lines refer to the parent model (without SF), passive and active star formation schemes, respectively, with ηSF = 10−1 . As a reference, the black dotted line represent the mass injection (in M yr−1 ) of the passively evolving stellar population (note that each model in a given row shares the same initial stellar population).

5.1 Star formation in ETGs

1041

91

No SF

LX (erg s−1)

1040 1039 1038 1037 1036

Passive SF ηSF = 10−1

1041

1040

1040

1039

1039

1038

1038

1037

1037

1036

1036

1041

LX (erg s−1)

Passive SF ηSF = 10−2

300

Active SF ηSF = 10−2

Active SF ηSF = 10−1

LX (erg s−1)

LX (erg s−1)

1041

200 250 σe8 (km s−1)

1041

1040

1040

1039

1039

1038

1038

1037

1037

1036

LX (erg s−1)

150

1036 150

200 250 σe8 (km s−1)

300

150

200 250 σe8 (km s−1)

300

Figure 5.7: ISM X-ray luminosity LX in the 0.3–8 keV band at 13 Gyr for the selected sub-group of the NFW set (see Section 5.1.1), as a function of σe8 . The same colour and symbol scheme of Figure 4.6 is adopted. For reference, the top panel refers to the original models (i.e. without star formation), while middle and bottom rows concern the passive and active SF simulations, respectively. Left and right panels correspond to low (10−2 ) and high (10−1 ) star formation efficiencies.

92

Ongoing and future developments

ones, LX is only marginally sensible to the presence of SF processes, while for the IS models the variation of the X-ray luminosity can be remarkable, thus increasing the LX spread at low galaxy masses.

5.1.3

Conclusions

In this chapter, we employed numerical hydrodynamical simulations to study the star formation in ETGs, in particular to test the robustness of our previous results regarding their the X-ray properties. The previous exploration of a large set of galaxies revealed the crucial role of the angular momentum in determining the X-ray under-luminosity of rotating ETGs. In addition, a major consequence of angular momentum conservation is the creation of massive, cold and rotationally supported discs. These dense discs are the ideal place for star formation processes that, in principle, can alter the hot ISM evolution, by acting as a sink of cold gas in the equatorial plane. Moreover, an eventual new stellar population will inject mass and energy into the ISM via SNII explosions, thus adding a new heating term to those associated with a passively evolving stellar population. As a consequence, the entire evolution of our sample of ETGs may be altered, and our previous conclusions may be dependent from the over-simplistic assumptions of no star formation made in Chapters 3 and 4. In order to investigate the effects of star formation on the X-ray halo of ETGs, we re-simulated a sub-group of the previous study NFW set of models, with two different SF schemes, the first consisting in a cold gas removal from the numerical grid (passive SF), while the second accounts also for the mass and energy injection by SNIIe explosions (active SF). Two different SF efficiencies has been also adopted to better explore the parameter space. We selected two different families of galaxies for this study: the high mass σe8 = 300 km s−1 family of isotropic rotators, which hosts the most massive rotating cold discs, and the low mass σe8 = 200 km s−1 family, consisting of both IS and VD galaxies, due to its sensitivity to variation of shape and kinematics, and to the fact that rotating models in this mass range present major recurring cooling episodes that can fuel SF processes. We summarize the results by discussing first the SFR and then the SF effects on the X-ray emission. 1) In low mass (σe8 = 200 km s−1 ) VD galaxies, the cold central core formed at early times is the only reservoir of cold material, due to the lack of CEs, and it is completely converted into stars, with a rate at 13 Gyr between 0.3–0.5 M yr−1 . The hot halo hydrodynamical evolution is not affected by the SF. 2) In low mass IS models, the hot atmosphere evolution can be substantially altered by the introduction of the SF. Cold discs are formed at early times, and completely converted into stars on a time-scale of few Gyr. Afterwards, SF is mainly fuelled by CEs of the X-ray halo, which prompt peaks of the SFR followed by a a secular decline, with a final value of 0.01–0.69 M yr−1 , with a mean of 0.14 M yr−1 . Overall, the hot halo final configuration is very susceptible to the action of the SF that, in general, contributes to induce a galactic wind, thus confirming the sensitivity of the σe8 = 200 km s−1 galaxies to small variations of structural and physical parameters. 3) In the case of high mass (σe8 = 300 km s−1 ) isotropic rotators, the X-ray corona do not present major CEs and its global evolution is only marginally sensible

5.1 Star formation in ETGs

93

to SF. The massive, cold discs formed at early times represent the only source of material for the star formation processes: they are partially or totally consumed during an Hubble time, depending of the star formation efficiency, presenting a smooth decline of the SFR with time, with a final value up to 3 M yr−1 . 4) The SFRs are well in accordance with the observed range of ≈ 0.01−3 M yr−1 (Davis et al. 2014) for the ATLAS3D sample of ETGs, where galaxies similar to our low mass family present SFRs in the 0.02–0.14 M yr−1 range. In addition, all the simulations suggest that the mean formation time for new stellar population is relatively short, ranging from 2.8 to 6.6 Gyr. Therefore, this excludes the possibility that young stellar discs are detected at the present epoch. 5) The effects of SF on the hydrodynamical evolution of the hot ISM are reflected on the ISM X-ray luminosity LX and luminosity weighted temperature TX . For VD models of low galactic mass, the evolution of LX and TX is not affected by the presence of SF, due to the fact that the cold ISM is accumulated in the central region and the SNIIe released energy is irradiated outside the X-ray band. In the case of IS models of the same mass, the additional CEs induced by the SF produce large oscillations in LX and TX , thus increasing the LX spread at low σe8 . 6) Overall, for high mass rotating models, LX and TX are only marginally sensible to the effects of SF. Their evolution is smooth as the respective parent model, since every system in this mass range is energetically dominated by the heating from the massive initial stellar population.

0.1

0.1

0.1

0.01

0.01

0.01

0.1

0.1

0.1

0.01

0.01

0.01







(2)

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

×

×

×

×

(1)

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

EO4300 IS

EO7300 IS 300 FO4IS FO7300 IS

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

71.93

71.81

65.46

69.38

73.01

72.94

66.56

70.53

73.01

72.94

66.56

70.53

73.01

72.94

66.56

70.53

73.01

72.94

66.56

70.53

(10 M ) (4)

9

Minj

11.14

11.93

15.87

14.10

10.83

11.75

15.68

13.85

10.87

11.71

15.77

13.89

10.86

11.79

15.71

13.83

11.04

11.89

15.89

13.95

61.82

60.92

50.24

56.04

16.13

14.15

17.02

14.87

6.86

7.73

7.82

8.23

17.79

15.06

18.53

15.98

7.40

8.04

8.23

8.55

9

Mgas

5.07

6.11

6.13

6.52

5.07

6.23

6.50

6.86

5.30

6.40

6.45

6.90

5.36

6.27

6.22

6.95

5.72

6.61

6.58

7.05

9

Mhot

56.7

54.7

44.5

49.8

10.9

7.85

10.4

7.93

1.47

1.22

1.31

1.24

12.4

8.73

12.3

8.97

1.59

1.36

1.58

1.42

9

Mc









46.06

47.04

33.86

41.81

55.29

53.50

42.96

48.41

55.45

57.60

40.38

50.89

68.21

66.25

53.03

60.03

9

∗ Mnew

























11.08

11.51

8.07

10.17

13.63

13.24

10.60

12.00

9

II Minj

hti∗









4.98

4.90

5.34

5.02

4.19

4.25

4.42

4.36

5.07

4.94

5.47

5.12

4.23

4.29

4.50

4.43

(10 M ) (10 M ) (10 M ) (10 M ) (10 M ) (10 M ) (Gyr) (5) (6) (7) (8) (9) (10) (11)

9

Mesc









9.25

9.6

6.34

8.33

13.2

12.59

9.72

11.10

10.94

11.66

7.38

9.94

16.13

15.44

11.78

13.55

(M yr (12)

−1

SFR









2.54

2.47

2.10

2.24

2.03

2.12

1.66

1.99

3.09

3.21

2.67

3.08

2.49

2.65

2.07

2.51

) (M yr−1 ) (13)

∗ Mnew /hti∗

Notes. (1) Name of the model. (2) Star formation scheme adopted (Active, Passive, not present). (3) Star formation efficiency. (4) − (5) Total ISM mass injected into and escaped from the numerical grid, respectively. Differences in Minj for models of same LB are accounted for different sampling of ρ∗ over the numerical grid. (6) Total ISM mass retained within the galaxy at the end of the simulation. (7) − (8) ISM mass with T > 106 K, and T < 2 × 104 K, respectively. (9) Mass of the new stellar population. (10) Cumulative mass ejected by type II SNe. (11) − (13) mean formation time with respect to the simulation beginning, when the old initial stellar population is already 2 Gyr old, mean star formation rate and today star formation rate of the new stellar population.



0.01

0.1

0.01

0.1

(3)

SF

Name

ηSF

Table 5.1: Simulations results for the σe8 = 300 km s−1 IS family of the NFW set with star formation at t = 13 Gyr. 94 Ongoing and future developments

0.1

0.1

0.1

0.01

0.01

0.01

0.1

0.1

0.1

0.01

0.01







(2)

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

×

×

×

×

(1)

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

EO4300 IS

EO7300 IS 300 FO4IS FO7300 IS

EO4300 IS

EO7300 IS 300 FO4IS FO7300 IS

EO4300 IS

EO7300 IS

FO4300 IS FO7300 IS

EO4300 IS

EO7300 IS FO4300 IS FO7300 IS

(10 (4)

erg s

0.900

1.35

1.04

1.34

0.885

1.38

1.07

1.47

0.915

1.50

1.04

1.54

0.944

1.44

0.991

1.56

1.09

1.65

1.13

1.65

40

LX −1

)

(10 (5)

erg s

49.49

49.40

45.04

47.74

50.23

50.18

45.79

48.53

50.23

50.18

45.79

48.53

50.23

50.18

45.79

48.53

50.23

50.18

45.79

48.53

40

LSN −1

) (10 (6)

erg s

























11.82

11.91

10.15

11.49

9.42

10.11

7.76

9.83

40

LII −1

) (10 (7)

erg s

























1.91

4.55

1.38

4.11

1.57

3.81

1.08

3.46

40

LII kin −1

) (10

−1

9.74E-2

1.72E-1

7.85E-2

1.60E-1

9.71E-2

1.72E-1

7.91E-2

1.60E-1

9.68E-2

1.72E-1

7.89E-2

1.60E-1

9.71E-2

1.71E-1

7.89E-2

1.60E-1

9.65E-2

1.71E-1

7.90E-2

1.59E-1

(8)

erg s

Lkin 40

0.59

0.68

0.56

0.68

0.55

0.68

0.58

0.69

0.58

0.68

0.59

0.68

0.57

0.68

0.58

0.69

0.57

0.70

0.58

0.69

(9)

) (keV)

TX

0.25

0.45

0.23

0.43

0.25

0.45

0.23

0.43

0.25

0.45

0.23

0.43

0.25

0.45

0.23

0.43

0.25

0.45

0.23

0.43

(10)

(keV)

Tkin

0.23

0.43

0.20

0.42

0.23

0.43

0.20

0.42

0.23

0.43

0.20

0.42

0.23

0.43

0.20

0.42

0.23

0.43

0.20

0.42

(11)

(keV)



0.06

0.08

0.06

0.07

0.06

0.07

0.06

0.06

0.06

0.07

0.06

0.07

0.06

0.07

0.06

0.07

0.06

0.07

0.06

0.06

(12)

γth

Notes. (1) Name of the model. (2) Star formation scheme adopted (Active, Passive, not present). (3) Star formation efficiency. (4) ISM X-ray luminosity in the 0.3–8 keV band. (5) SNIa heating rate at the end of the simulation. (16) Type II SNe luminosity. (7) − (8) Kinematical heating of the new and primordial stellar population (under the assumpion that the new stellar population has the same kinematical feature of the primordial one). (9) ISM X-ray emission weighted temperature in the same band, at the end of the simulation. (10) − (11) Thermalization temperatures of stellar motions at the end of the simulation, defined accordingly to equations (2.35)-(2.34). By construction, Tkin = Tσ + Tv ; ϕ for rotating models Tv = γth Trot and Tϕ = Tv − Tm = γth Trot , while for velocity dispersion supported models Tv = Tm (see Chapter 1 and Section 2.4.3). (12) Thermalization parameter as defined in equation (2.30).



0.01

0.01

0.1

0.01

0.1

(3)

SF

Name

ηSF

Table 5.2: Simulations results for the σe8 = 300 km s−1 IS family of the NFW set with star formation at t = 13 Gyr.

5.1 Star formation in ETGs 95

0.1

0.1

0.1

0.01

0.01

0.01

0.1

0.1

0.1

0.01

0.01

0.01

(2)

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

×

×

×

×

(1)

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS 200 FO4IS FO7200 IS

EO4200 IS

EO7200 IS

FO4200 IS

FO7200 IS

11.96

12.00

11.93

11.92

12.08

12.14

12.09

12.07

12.08

12.14

12.09

12.07

12.08

12.14

12.09

12.07

12.08

12.14

12.09

12.07

(10 M ) (4)

9

Minj

3.41

3.32

10.34

9.76

2.98

6.65

10.43

4.32

5.79

3.75

10.07

4.65

3.09

6.77

9.70

9.40

3.82

6.00

8.27

7.77

8.69

8.80

1.84

2.39

3.12

0.68

0.34

3.28

0.49

3.64

0.32

2.00

3.77

0.98

0.47

0.27

3.82

1.48

0.80

0.29

9

Mgas

1.36

1.73

0.24

0.15

1.57

0.43

0.27

1.63

0.40

2.88

0.29

1.77

1.42

0.72

0.24

0.17

3.28

1.19

0.62

0.28

9

Mhot

7.16

6.97

1.57

2.24

1.40

2.19E-1

3.99E-2

1.57

1.98E-2

8.98E-4

0.00

7.31E-2

2.19

1.86E-1

2.09E-1

9.92E-2

2.55E-3

8.70E-4

5.66E-2

7.34E-4

9

Mc









5.99

4.80

1.32

4.47

5.81

4.74

1.71

5.42

6.53

5.47

2.40

3.00

5.57

5.81

3.78

5.02

9

∗ Mnew

























1.30

1.09

0.48

0.60

1.11

1.16

0.76

1.00

9

II Minj

hti∗









5.39

4.20

3.84

6.61

3.43

3.28

2.37

5.74

5.27

4.11

4.39

3.90

2.92

3.10

5.34

3.72

(10 M ) (10 M ) (10 M ) (10 M ) (10 M ) (10 M ) (Gyr) (5) (6) (7) (8) (9) (10) (11)

9

Mesc

Notes. All quantites are as in Table 5.1.









0.01

0.1

0.01

0.1

(3)

SF

Name

ηSF









1.11

1.14

0.34

0.68

1.69

1.45

0.72

0.94

1.24

1.33

0.55

0.77

1.91

1.88

0.71

1.35

(M yr (12)

−1

SFR









0.37

0.08

0.02

0.45

0.04

0.08

0.00

0.12

0.69

0.09

0.06

0.05

0.11

0.01

0.06

0.00

) (M yr−1 ) (13)

∗ Mnew /hti∗

Table 5.3: Simulations results for the σe8 = 200 km s−1 IS family of the NFW set with star formation at t = 13 Gyr.

96 Ongoing and future developments

0.1

0.1

0.1

0.01

0.01

0.01

0.1

0.1

0.1

0.01

0.01

0.01







(2)

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

×

×

×

×

(1)

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS FO4200 IS FO7200 IS

EO4200 IS

EO7200 IS 200 FO4IS FO7200 IS

EO4200 IS

EO7200 IS 200 FO4IS FO7200 IS

(10 (4)

erg s

LX −1

2.48E-1

4.49E-1

1.95E-3

7.29E-4

1.63E-1

5.77E-3

3.00E-3

1.10E-1

9.73E-3

5.87E-1

2.82E-3

8.64E-2

4.60E-1

1.09E-2

2.15E-3

1.01E-3

9.82E-1

4.08E-2

1.15E-2

1.61E-3

40

)

Notes. All quantites are as in Table 5.2.



0.01

0.1

0.01

0.1

(3)

SF

Name

ηSF (10 (5)

erg s

10.20

10.25

10.19

10.17

10.31

10.36

10.32

10.30

10.31

10.36

10.32

10.30

10.31

10.36

10.32

10.30

10.31

10.36

10.32

10.30

40

LSN −1

) (10 (6)

erg s

LII

























2.58

0.36

0.24

0.18

0.40

0.05

0.24

0.01

40

−1

) (10 (7)

erg s

























0.15

0.06

0.02

0.03

0.04

0.01

0.02

0.00

40

LII kin −1

) (10 (8)

erg s

−1

9.70E-3

1.63E-2

2.31E-2

3.17E-2

1.12E-2

2.30E-2

2.43E-2

1.57E-2

1.93E-2

1.74E-2

2.44E-2

1.70E-2

9.50E-3

2.14E-2

2.25E-2

3.06E-2

1.29E-2

2.04E-2

1.58E-2

3.04E-2

40

Lkin

0.32

0.37

0.49

0.50

0.36

0.40

0.50

0.39

0.50

0.43

0.51

0.40

0.29

0.43

0.49

0.54

0.41

0.44

0.44

0.53

(9)

) (keV)

TX

0.12

0.21

0.29

0.40

0.14

0.29

0.31

0.20

0.25

0.22

0.31

0.22

0.12

0.27

0.29

0.39

0.16

0.26

0.20

0.39

(10)

(keV)

Tkin

0.09

0.18

0.08

0.17

0.09

0.18

0.08

0.17

0.09

0.18

0.08

0.17

0.09

0.18

0.08

0.17

0.09

0.18

0.08

0.17

(11)

(keV)



0.15

0.31

1.24

2.62

0.25

1.22

1.33

0.36

0.80

0.44

1.34

0.53

0.14

1.00

1.19

2.47

0.37

0.86

0.70

2.44

(12)

γth

Table 5.4: Simulations results for the σe8 = 200 km s−1 IS family of the NFW set with star formation at t = 13 Gyr.

5.1 Star formation in ETGs 97

0.1

0.1

0.1

0.01

0.01

0.01

0.01

0.1

0.1

0.1

0.1

0.01

0.01

0.01

0.01









(2)

A

A

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

P

P

×

×

×

×

×

(1)

E0200

EO4200 VD

EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD 200 FO4VD FO7200 VD

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

11.96

12.00

11.93

11.92

12.00

12.08

12.14

12.09

12.07

12.14

12.08

12.14

12.09

12.07

12.14

12.08

12.14

12.09

12.07

12.14

12.08

12.14

12.09

12.07

12.14

(4)

(109 M )

Minj

Mgas

Mhot

Mc

∗ Mnew II Minj

hti∗

4.00

2.94

9.37

4.51

3.19

3.90

2.91

8.89

4.48

3.25

3.77

2.96

9.35

4.58

3.27

3.89

2.92

8.64

4.49

3.21

3.87

3.05

9.29

4.68

3.42

(5)

8.12

9.19

2.80

7.54

8.92

1.84

2.63

1.08

2.99

2.96

1.83

2.57

0.65

2.90

2.94

1.79

2.64

1.30

3.00

3.05

1.92

2.62

0.78

2.95

2.95

(6)

1.51

2.46

0.66

2.84

2.91

1.61

2.53

0.97

2.92

2.90

1.66

2.54

0.56

2.88

2.93

1.56

2.53

1.16

2.93

2.96

1.74

2.59

0.69

2.93

2.93

(7)

6.44

6.68

2.05

4.69

5.99

6.39E-2

7.13E-2

2.29E-3

6.40E-2

6.27E-2

1.35E-2

1.37E-2

4.04E-4

1.17E-2

9.66E-3

8.04E-2

8.36E-2

3.22E-3

7.38E-2

8.91E-2

1.25E-2

1.57E-2

4.74E-4

1.30E-2

1.45E-2

(8)











6.35

6.60

2.11

4.60

5.93

6.48

6.60

2.09

4.59

5.92

8.00

8.22

2.69

5.72

7.34

7.88

8.08

2.54

5.55

7.21

(9)































1.60

1.64

0.54

1.14

1.47

1.57

1.62

0.51

1.11

1.44

(10)











4.30

4.72

3.00

4.89

4.92

4.31

4.60

2.76

4.74

4.78

4.37

4.74

3.14

4.90

4.91

4.22

4.59

2.78

4.70

4.76

(11)

(109 M ) (109 M ) (109 M ) (109 M ) (109 M ) (109 M ) (Gyr)

Mesc

Notes. All quantites are as in Table 5.1.



0.01

0.1

0.01

0.1

0.1

(3)

SF

Name

ηSF

SFR











1.48

1.40

0.70

0.94

1.21

1.50

1.43

0.76

0.97

1.24

1.83

1.73

0.86

1.17

1.49

1.87

1.76

0.91

1.18

1.51











0.28

0.33

0.00

0.28

0.28

0.34

0.36

0.00

0.31

0.22

0.43

0.41

0.00

0.36

0.45

0.34

0.40

0.00

0.32

0.37

(M yr−1 ) (M yr−1 ) (12) (13)

∗ /hti∗ Mnew

Table 5.5: Simulations results for the σe8 = 200 km s−1 VD family of the NFW set with star formation at t = 13 Gyr.

98 Ongoing and future developments

0.1

0.1

0.1

0.1

0.01

0.01

0.01

0.01

0.1

0.1

0.1

0.1

0.01

0.01

0.01

0.01







(2)

A

A

A

A

A

A

A

A

A

A

P

P

P

P

P

P

P

P

P

P

×

×

×

×

×

(1)

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD FO4200 VD FO7200 VD

E0200

EO4200 VD EO7200 VD 200 FO4VD FO7200 VD

E0200

EO4200 VD

EO7200 VD FO4200 VD FO7200 VD

(10

erg s (4)

−1

1.82

2.21

1.89E-2

1.99

2.06

1.86

2.21

3.61E-2

1.96

1.80

2.18

2.24

1.41E-2

1.64

1.51

2.46

2.18

5.60E-2

1.93

3.01

1.83

2.26

1.93E-2

1.63

1.74

40

LX )

Notes. All quantites are as in Table 5.2.





0.01

0.1

0.01

0.1

(3)

SF

Name

ηSF (10

erg s (5)

10.20

10.25

10.19

10.17

10.24

10.31

10.36

10.32

10.30

10.36

10.31

10.36

10.32

10.30

10.36

10.31

10.36

10.32

10.30

10.36

10.31

10.36

10.32

10.30

10.36

40

LSN −1

) (10

erg s (6)































1.50

1.63

0.02

1.37

1.74

1.26

1.56

0.02

1.25

1.46

40

LII −1

) (10

erg s (7)































0.25

0.32

0.00

0.18

0.29

0.21

0.28

0.00

0.17

0.23

40

LII kin −1

) (10

erg s (8)

−1

2.23E-2

2.17E-2

2.59E-2

2.04E-2

2.09E-2

2.24E-2

2.16E-2

2.44E-2

2.04E-2

2.09E-2

2.23E-2

2.15E-2

2.68E-2

2.03E-2

2.08E-2

2.23E-2

2.16E-2

2.34E-2

2.04E-2

2.10E-2

2.23E-2

2.15E-2

2.59E-2

2.03E-2

2.09E-2

40

Lkin

0.49

0.48

0.50

0.46

0.48

0.49

0.49

0.49

0.46

0.48

0.49

0.47

0.50

0.47

0.49

0.47

0.49

0.49

0.46

0.45

0.49

0.47

0.50

0.47

0.48

(9)

) (keV)

TX

0.28

0.27

0.33

0.26

0.26

0.28

0.27

0.31

0.26

0.26

0.28

0.27

0.34

0.26

0.26

0.28

0.27

0.30

0.26

0.27

0.28

0.27

0.33

0.26

0.26

(10)

(keV)

Tkin

0.28

0.27

0.25

0.26

0.26

0.28

0.27

0.25

0.26

0.26

0.28

0.27

0.25

0.26

0.26

0.28

0.27

0.25

0.26

0.26

0.28

0.27

0.25

0.26

0.26

(11)

(keV)





















































(12)

γth

Table 5.6: Simulations results for the σe8 = 200 km s−1 VD family of the NFW set with star formation at t = 13 Gyr.

5.1 Star formation in ETGs 99

100

5.2

Ongoing and future developments

AGN feedback in ETGs

In Chapters 3 and 4 we showed the effects of the stellar dynamics on the evolution of the X-ray coronae of ETGs by means of 2D hydrodynamical simulation with realistic galaxy models. In particular, we found that the angular momentum plays a fundamental role in determining the hydrodynamical evolution of the hot rotating haloes and their X-ray under-luminosity. However, the massive rotationally supported cold discs that forms due to angular momentum conservation are not observed in ETGs. By including the fundamental physics of star formation, we showed in the last Section how the massive cold discs present in our simulated rotating galaxies are almost completely consumed by star formation processes, producing realistic values of the SFR. At this point of our investigation, the missing natural ingredient is the AGN feedback, given the ubiquitous presence of black holes (BHs) at the centre of massive galaxies (Gebhardt et al. 2000; Ferrarese & Merritt 2000; Tremaine et al. 2002; Novak et al. 2006; Gültekin et al. 2009; Fabian 2012; Kormendy & Ho 2013). Simulations of feedback in ETGs has been carried in a sequence of papers by means of 1D hydrodynamical simulations (Ciotti & Ostriker 1997, 2001, 2007; Ciotti et al. 2009, 2010; Ostriker et al. 2010; Shin et al. 2010; Jiang et al. 2010). Despite their spherical symmetry, the implemented micro-physics is extremely accurate, accounting for the radiation pressure due to the accretion luminosity, Compton heating, radiative transport and a physically motivated model of AGN mechanical feedback. Recently, an extension to bi-dimensional, axisymmetric simulations has been performed by Novak et al. (2011), and further developed in Novak et al. (2012) and Gan et al. (2014). The transition from 1D to 2D simulations of an axisymmetric galaxy is a crucial step to include galaxy flattening and rotation, and capture a wide spectrum of physical phenomena, e.g. the breaking of a cold, falling gas shell by means of the Kelvin-Helmholtz instability, which cannot be simulated in a spherically symmetric code. Although the great accuracy of the implemented micro-physics, in these 2D codes the galaxy is a spherically symmetric idealised model and galaxy rotation is not self-consistent. Moreover, the stellar streaming velocity is kept at low values and artificially tuned in order to allow the gas produced by the ageing stellar population to flow inside the first computational grid, without the formation of a rotationally supported disc due to the centrifugal barrier, as it occurs with realistic velocity fields; even specific studies on rotating accretion flows such as Li et al. (2013) adopt a constant specific angular momentum distribution built ad-hoc to allow accretion inside the Bondi radius. We started a joint scientific collaboration with Dr. Novak, with the principal goal to simulate the evolution of realistic, flat and rotating state-of-the-art galaxies in presence of AGN feedback. In order to reach this aim, we implemented our realistic galaxy models and passive stellar population feedback into the spherical code of Novak et al. (2012). The main reason relies in the fact that the entire treatment of the micro-physics is by far more technically and physically complicated than the stellar population feedback, therefore we adopted as a base the very robust Novak’s code. Preliminary results of flattened velocity dispersion supported galaxies with AGN feedback are presented in Section 5.2.2. By implementing galaxy models with realistic rotation curves, the formation of a

5.2 AGN feedback in ETGs

101

102

104

Q

102 101

100

Σ (MO· pc−3)

103

101

100 10-1 10-1

0

1

10

2

10

10

-1

10

R (kpc)

0

1

10

10

10−1 102

R (kpc)

Figure 5.8: Toomre Q parameter and surface density Σ for the cold gaseous disc of EO7300 IS at 13 Gyr. massive, cold rotating gaseous disc hosting star formation is expected, due to the angular momentum conservation of the non-viscous ISM. As a consequence, the BH accretion is strongly inhibited, and the cold material is completely consumed by star formation. However, the above argument does not take into account the ISM self-gravity: when such a rotating disc reaches the self-gravitating regime, it becomes prone to gravitational instabilities, which transport angular momentum from the disc inner regions to the outer edge, thus allowing the accretion of the rotating material. The local stability of a self-gravitating disc can be expressed with the aid of the Toomre (1964) Q parameter, defined as Q≡

cs κ , πGΣ

κ2 =

∂2 Φ + 3Ω2 , ∂R2

(5.1)

where cs is the sound crossing time, κ the epicyclic frequency, Σ the disc surface density, and Ω the circular angular velocity. For illustrative purposes, in Figure 5.8 we calculated the Q parameter of the cold disc for the EO7300 IS simulation without star formation, at t = 13 Gyr. It is apparent that in the disc denser region, the Toomre stability criterion Q > 1 is violated. Of course this is an extremely simplified picture, since we are implicitly assuming that the rotationally supported gaseous disc is completely dominated by its own gravity field. The actual stability criterion is more complicated, since also the external gravitational fields of stars and DM halo should be taken into account (e.g. see Jog 1996; Romeo & Falstad 2013; Jog 2014). In addition, Figure 5.8 refers to one of the most massive disc at the end of the simulation, when almost all the ISM has cooled down to the central regions. However, it is reasonable to expect that the cold disc would develop gravitationally unstable regions during its evolution. Unfortunately, a 2D axisymmetric code is intrinsically unable to capture the three-dimensional break of a rotating disc. Various methods has been proposed over the years in order to transport angular momentum in a gaseous disc (for a review, see Blaes 2014), each with their pros and cons. In particular, we adopted a numerical

102

Ongoing and future developments

scheme that mimics a gravitational instability via the inclusion of a “gravitational viscosity”, employing the same mathematical formalism of the fluid viscosity (already applied in studies with non-public versions of the ZEUS code in spherical coordinates, see Stone et al. 1999; Li et al. 2013). We are implementing a coordinate independent version of the gravitational viscosity formalism in order to include it in both the ZEUS-MP 2 code (i.e. the code described in Section 2.4, which can handle both cylindrical and spherical coordinates), and in the spherical customized AGN code presented in this Section. The Section is organized as follows. In Section 5.2.1 a series of runs built on purpose to test the new implemented physics is presented. Section 5.2.3 shows the detailed viscosity formalism and its coordinate independent implementation are shown. Finally, in Section 5.2.2 the preliminary simulations of flattened velocity dispersion supported galaxies are presented.

5.2.1

Code implementation and testing

In order to implement our realistic galaxy models and passive stellar population feedback into the spherical code of Novak et al. (2012), we split the implementation phase in two different steps: first, we interfaced the realistic galaxy models produced by the Jeans code of P13 with the hydro code. We took advantage of the numerical interpolation code presented in Section 2.1, which can handle both spherical and cylindrical staggered hydro meshes. Thus, the interface between the galaxy models and the hydrodynamical code is the same of our version of ZEUS-MP2. The porting of the routines that actually load the interpolated galaxy model did non present any particular problem, since it were already designed in a coordinate independent fashion. Afterwards, the passive stellar population mass (taken from Ciotti & Ostriker 2007), energy and momentum original injections terms has been replaced with the corresponding up-to-date self-consistent terms present in Section 2.4.1. We present here a set of tests performed with the purpose of reproducing the simulations of Novak et al. (2012) with the same physics, in order to track down implementation errors in the up-to-date feedback from the passive stellar population. In particular, we tested each of the mass, momentum and energy injection terms separately, by activating the new implementation for one of them, and leaving the remaining ones unaltered. In Table 5.7, simulations 1-8 test the stellar winds injecting terms without the AGN feedback (but having the central BH sink activated), since the large energy injections due to AGN and SNIa can mask errors in the stellar winds thermalization. In the same way, we performed the four additional runs 25-28 to exclude the effects of the AGN feedback over the SNIa energy injection. As initial conditions, we adopted the galaxy model of Novak et al. (2012), consisting of a Jaffe profile with a total mass of 3 × 1011 M and a projected halfmass radius of 6.9 kpc, while the total gravitational potential is a singular isothermal sphere plus the point-mass potential of the central BH. The stellar population is assumed to rotate as a solid body inside the first 10 pc, whereas a constant specific angular momentum is assumed at large radii. A trend common to every group of simulations is the lack of difference between the old and new momentum injection terms, due to the very low galaxy rotational velocity. We focussed on the evolution of three different quantities: the SFR, and the

5.2 AGN feedback in ETGs

103

Table 5.7: Test simulations of the mass, momentum and energy injection terms. Name θSN AGN ρ˙ E˙ v˙ ϕ (1) (2) (3) (4) (5) (6) 1 2 3 4 5 6 7 8

10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2

× × × × × × × ×

old new old new old new old new

old old new new old old new new

old old old old new new new new

9 10 11 12 13 14 15 16

10−2 10−2 10−2 10−2 10−2 10−2 10−2 10−2

X X X X X X X X

old new old new old new old new

old old new new old old new new

old old old old new new new new

17 18 19 20 21 22 23 24

1 1 1 1 1 1 1 1

X X X X X X X X

old new old new old new old new

old old new new old old new new

old old old old new new new new

25 26 27 28

1 1 1 1

× × × ×

old new old new

old old new new

new new new new

Notes. (1) Test name. (2) present day SNIa explosion rate. (3) AGN feedback presence (X = present, × = not present). (4)-(6) Implementation of mass, energy and momentum injection terms, respectively (old/new).

104

Ongoing and future developments

BH accretion rate and mass. Figure 5.9 reports the most relevant tests: it is evident that the new injection mass and energy terms do not change the BH accretion and the SFR in simulations 5-8 (first row in Figure 5.9). In the case of a significant energy input from SNIa (last row in Figure 5.9), a burst of SF occurs with a SFR more peaked than the old SNIa rate, although the final gas mass converted in stars is the same. The action of the AGN feedback produce a noisy evolution, but in general all the tests show that the newly implemented passive population feedback is correct.

5.2.2

Preliminary results

250 We simulated three different models: the E0250 , EO4250 VD and EO4IS of the Einasto set (Table 4.1). The input physics is described in Novak et al. (2011), except for our upgrade of the mass, momentum and energy return of the passive stellar population. Figures 5.10-5.12 report the BH Eddington ratio, hot and cold gas mass, SFR, stellar mass and AGN duty cycle. We do not present here a complete description of their time evolutions, however the difference between the velocity dispersion supported models (both E0 and E4) and the rotating one is apparent: E0 and E4 VD galaxies show little difference, while rotation strongly affects the IS model. Indeed, the BH accretes at a rate up to the Eddington limit in the case of VD models, while in the IS model the BH luminosity ratio is very sub-Eddington (Figure 5.10). This is due to the cold disc formation, where all the cooled mass of the IS model is confined (Mc ' 8 × 109 M ) and slowly consumed by star formation during the whole evolution (Figure 5.11). On the contrary, the E0 and E4 VD models present no cold gas, since it is directly accreted on the BH, and recurrent star formation bursts, coincident with the peaks of AGN activity. The difference between the VD and IS models is even more evident from Figure 5.12, where the AGN duty cycle is reported: all the simulations spent most of their time at very sub-Eddington accretion rates (the median in Figure 5.12 left panel is ' 10−8 ), but an half of the AGN energy has been emitted with an Eddigton ratio greater that 0.1 in the case of VD models, while the IS one emits almost all its energy below 10−3 . To sum up, for VD models the vast majority of the energy is emitted at high Eddington ratio (as observed in Kollmeier et al. 2006), while the angular momentum conservation inhibits the BH accretion in rotating models. As a consequence, a mechanism of angular momentum transport is needed in order to induce the accretion of the cold disc on the central BH. To this end, in the next Section we present an angular momentum transport algorithm based on the inclusion of a “gravitational viscosity” in the equations of hydrodynamics.

5.2.3

Gravitational viscosity

We present here the adopted angular momentum algorithm, based on the viscosity formalism (presented in Section A.1.5). The aim of this method is to induce a non-conservation of the specific angular momentum Jz ≡ Ruϕ , so as to mimic the effects of gravitational instabilities. Following Stone et al. (1999), we assumed that the only non-zero components of the viscous stress tensor are the azimuthal ones,

5.2 AGN feedback in ETGs

105

5

5

10 15 t (Myr)

-1 Mrate BH (MO • yr )

2 1

6 5

5

3 2 1

6 5

5

4 3 2 1 0

5

10 15 t (Myr)

5

10 15 t (Myr)

0

5

10 15 t (Myr)

0

5

10 15 t (Myr)

0

5

10 15 t (Myr)

102 101 100 10-1 10-2 10-3 10-4

10 15 t (Myr)

25 26 27 28

0 102 101 100 10-1 10-2 10-3 10-4

10 15 t (Myr)

21 22 23 24

4

102 101 100 10-1 10-2 10-3 10-4

10 15 t (Myr)

3

0

MBH ( 108 MO•)

SFR (MO• yr-1)

10 0

4

10 15 t (Myr)

100

5

13 14 15 16

0

MBH ( 108 MO•)

SFR (MO• yr-1)

10 5

6 5

10 15 t (Myr)

100

0

2 1 0

MBH ( 108 MO•)

SFR (MO• yr-1)

10 5

3

10 15 t (Myr)

100

0

4

-1 Mrate BH (MO • yr )

0

5 6 7 8

-1 Mrate BH (MO • yr )

10

6 5

-1 Mrate BH (MO • yr )

MBH ( 108 MO•)

SFR (MO• yr-1)

100

102 101 100 10-1 10-2 10-3 10-4

Figure 5.9: SFR, BH mass and accretion rate as a function of time from the begin of the simulation, for few representative tests (see Table 5.7). From top to bottom, each line refers to the group of four simulations listed in the legend.

106

Ongoing and future developments

100 10-2

1044

10-4

1042

10-6

1040

10-8

1038

10-10 10

2

3

4

5

6 7 t (Gyr)

8

9

10 2

3

4

5

6 7 t (Gyr)

8

9

Lbol BH/LEdd

-1 Lbol BH (erg s )

1046

100 10-2

1044

10-4

1042

10-6

1040

10-8

1038

10-10 10

2

3

4

5

6 7 t (Gyr)

8

9

10 2

3

4

5

6 7 t (Gyr)

8

9

Lbol BH/LEdd

-1 Lbol BH (erg s )

1046

100 10-2

1044

10-4

1042

10-6

1040

10-8

1038

10-10 10

2

3

4

5

6 7 t (Gyr)

8

9

10 2

3

4

5

6 7 t (Gyr)

8

9

Lbol BH/LEdd

-1 Lbol BH (erg s )

1046

Figure 5.10: Black hole bolometric luminosity and Eddington ratio evolution with 250 time. From top to bottom, each row refers to E0250 , EO4250 VD and EO4IS of the Einasto set (Table 4.1). The low accretion rate of the EO4250 IS model is apparent.

107

1010

109

109

108

Mcold (MO•)

Mhot (MO•)

5.2 AGN feedback in ETGs

E0 E4 VD E4 IS

108 2

3

4

5

6 7 t (Gyr)

8

9

10 2

3

4

5

6 7 t (Gyr)

8

9

100 109

M* (MO•)

1010

101 SFR (MO• yr-1)

107 10

10-1

10-2 2

3

4

5

6 7 t (Gyr)

8

9

10 2

3

4

5

6 7 t (Gyr)

8

9

108 10

Figure 5.11: Time evolution of the gaseous and stellar mass for the E0250 (green), 250 EO4250 VD (red) and EO4IS (black). Top row refers to hot and cold mass, while bottom row reports the SFR and the stellar mass evolution.

Ongoing and future developments

100

100

10−1

10−1

10−2

10−2 E0 E4 VD E4 IS

10−3 10−10 10−8

10−3 −6

−4

10 10 Lbol /L BH Edd

10

−2

0

−6

10 10

10

−5

−4

10

−3

−2

10 10 Lbol BH/LEdd

−1

10

10

cumulative energy above/below

cumulative time above/below

108

0

Figure 5.12: AGN duty cycle as a function of Eddigton ratio for the models of Figure 5.11 (the adopted colours are the same). The solid/dashed line show the cumulative time (left panel) and energy (right panel) above/below the given Eddington ratio. thus from equation (A.31) the components of T in spherical coordinates are ∂  uϕ  sin θ ∂  uϕ  Trϕ = ηr , Tθϕ = η , (5.2) ∂r r r ∂θ sin θ whereas in cylindrical coordinates assume the form ∂uϕ ∂  uϕ  TRϕ = ηR , Tzϕ = η , (5.3) ∂R R ∂z where we used the differential operators listed in Appendix B, and we imposed Trr = Trθ = Tθθ = TRR = TRz = Tzz = 0. Equation (A.16) reports the force on the fluid due to the viscosity as ∇· T ; given the assumed form of the viscous stress tensor, its non-vanishing component is the azimuthal one:   u  1 ∂  sin θ ∂  u  1 ∂ ϕ ϕ 2 ∂ ηr + η (∇· T )(ϕ) = r ∂r ∂r r r ∂θ r ∂θ sin θ ∂  uϕ  cos θ ∂  uϕ  + 2η + 2η 2 , (5.4) ∂r r r ∂θ sin θ     ∂uϕ ∂ ∂ ∂  uϕ  ∂  uϕ  (∇· T )(ϕ) = η + ηR + 2η ; (5.5) ∂z ∂z ∂R ∂R R ∂R R where the two equations express the force in spherical and cylindrical coordinates, respectively, and η is assumed to be a function of the position on the meridional plane only, i.e. it does not depends on the azimuthal angle ϕ. The effect of the viscosity is to reduce the fluid velocity gradients by transforming kinetic energy into internal. From equation (A.21), the viscosity internal energy input is T : ∇u that, under the above assumptions, becomes T : T /η. A complete description of the coordinate independent implementation is given in Appendix C. As a future development, we will employ the presented formalism to simulate rotating ETGs in presence of AGN feedback.

Chapter 6

Conclusions In this Thesis we studied the effects of galaxy shape, stellar kinematics and star formation on the evolution of the hot, X-ray emitting gaseous haloes of ETGs by means of high-resolution, axisymmetric hydrodynamical simulations of state-of-theart galaxy models, tuned to reproduce real galaxies. In particular, we focused our efforts on two different problems. The first main goal is to investigate the observed trend of LX and TX with galaxy shape and stellar kinematics, in particular the X-ray under-luminosity and coolness of flat and rotating galaxies. Precedent analytical studies (CP96; P13) proposed two different explanations based on energetic considerations. In both cases, the ISM in flat and/or rotating galaxies is less bound than in rounder and non-rotating objects of the same mass due to (1) the assumed shallower potential of flat galaxies (independently of the kinematical support), and (2) to the assumption that the injected mass in rotating galaxies retains the streaming motion of the stellar population, with the consequent minor heating due to the lack of thermalization of the stellar global rotation. Thus, flat and rotating galaxies are more prone to develop a global wind, with the consequent reduction of LX and hot mass content. Instead, the numerical study of DC98 seems to indicate the redistribution of the hot ISM due to angular momentum conservation as the main actor in determining the ETGs hot halo evolution. Furthermore, these simulations found a generalized coolness of rotating models, where cold filaments are formed. In order to understand the main physical mechanism governing the X-ray haloes evolution, their LX and TX , we performed two different investigations based on hydrodynamical simulations. Firstly we focussed on the sole effects of rotation and counter-rotation in S0 galaxies; subsequently, we explored variations of galaxy mass, shape, kinematical support and dark matter halo in a large set of realistic, stateof-the-art elliptical galaxy models, built with a dedicated Jeans code (P13; Posacki 2014). The hot haloes are formed by the mass losses of ageing stellar population (stellar winds from giant stars and SNIa ejecta), while ISM heating is provided by SNIa explosions and by thermalization of stellar motions. From the results of our investigations, whose numerical simulations reproduce the LX and TX of observed ETGs (Boroson et al. 2011; Sarzi et al. 2013), we conclude that the above explanations and the relative importance of flattening and rotation are functions of the galaxy mass. Indeed, at varying with galaxy mass, different physical mechanisms dominates in determining the hot halo evolution: 109

110

Conclusions

(1) At low galaxy masses (σe8 < 200 km s−1 ), spherical galaxies are in a permanent wind state, presenting low LX and high values of TX , due to the thermalization of the strong meridional motions; adding flattening and rotation has a minor impact on the halo evolution. (2) The value of σe8 ≈ 200 km s−1 represents a sort of threshold where galaxies are extremely sensitive to variations of shape and kinematics. In this mass range, where galaxies are energetically near to the onset of a wind, adding flattening or rotation contributes to induce a global wind, in agreement with the energetic explanations of CP96, with the consequent reduction of LX and increase of TX . (3) In galaxies having σe8 > 200 km s−1 , the angular momentum conservation becomes the main actor in determining the evolution of the X-ray coronae. Pure shape variations in velocity dispersion supported galaxies produce only a marginal effect; on the contrary, flat rotating galaxies shows low LX and TX . This is not a consequence of the missing thermalization of stellar streaming motions, due to the fact that ISM and stars co-rotate, but instead to the hydrodynamical configurations induced by angular momentum conservation. In fact, a cold disc forms and the centrifugal barrier keeps low the hot gas density in the central regions, thus producing a lower LX with respect to non-rotating galaxies, where a denser and brighter atmosphere forms. For the same reason, the ISM redistribution in rotating galaxies induces also recurring cooling episodes and a generalized lower TX than in non-rotating models. In this regime of mass, flattening is an indirect effect in the sense that only flattened systems can host significant rotation of the stellar component. The second main goal of this Thesis is to probe star formation in ETGs and its effects on the X-ray halo evolution. The previous explorations revealed the crucial role of the angular momentum in determining the halo evolution, with the creation of massive (up to 5 × 1010 M ), cold and rotationally supported discs by recurring cooling episodes. Such dense discs are the ideal place for star formation processes that, in principle, can alter the hot ISM evolution, by acting as a sink of cold gas and re-injecting mass and energy into the ISM via SNII explosions. Moreover, on the observational side, molecular discs and star formation regions are observed in fast rotators only (Davis et al. 2013, 2014): although the origin of this molecular gas has not been clarified yet, we investigated its link with the formation of cold discs due to instability of rotating X-ray coronae. We employed a sub-group of galaxy models of the previous studies to perform hydrodynamical simulations with different star formation schemes and efficiencies, in order to investigate its effects on the X-ray halo of ETGs. Feedback from SNII explosions and thermalization of ordered and random motions of the new stellar population was also considered. Simulations results quite agree with the observed present-day SFR in ETGs, which ranges from ' 0.02 to 0.14 M yr−1 (Davis et al. 2014), and show that star formation effects are a function of galaxy mass and stellar kinematics. Depending on galaxy mass, the entire reservoir of cold gas is partially or completely consumed by star formation: in the case of rotating low-medium mass galaxies, cold discs are formed and completely converted into stars in few Gyr, then star formation is sustained by the recurring cooling episodes of the rotating X-ray halo. In this mass range, the hot halo evolution in rotating galaxies is very susceptible to the action of the star formation that, in general, contributes to induce a galactic wind, thus increasing the

111

LX spread at low σe8 and confirming the sensitivity of the σe8 ≈ 200 km s−1 galaxies to small variations of structural and physical parameters. On the contrary, the halo evolution in velocity-dispersion-supported systems of similar mass is only marginally affected by star formation. In more massive rotating galaxies, star formation is fuelled entirely by massive cold discs, which have been partially or completely consumed at the present epoch, depending of the star formation efficiency. We conclude that the hot halo instabilities induced by angular momentum conservation, and thus present only in rotating galaxies, can explain both the origin of the molecular discs observed in fast rotators and their present-day SFRs. Finally, we dealt with the inclusion of AGN feedback in the presented set of numerical simulations. We implemented the realistic galaxy models of this Thesis and the self-consistent feedback from the passive stellar population in the numerical code of Novak et al. (2011, 2012), used to simulate AGN feedback with an accurate micro-physics but with simple spherical galaxy models. Preliminary results of nonrotating spherical and flattened galaxies present an AGN luminosity and duty cycle in good agreement with observations, where the majority of the AGN energy is emitted at high Eddington ratios. However, BH accretion is strongly inhibited in rotating models due to the formation of massive, cold, star forming discs. Cold disc fragmentation and bar instabilities cannot be captured by an axisymmetric code, thus we implemented a numerical scheme, based on the fluid viscosity formalism, that allows for angular momentum transport from the inner parts of the disc to the outer edge. As a future development, we will simulate the evolution of realistic stateof-the-art galaxy models in presence of AGN feedback, by means of the presented novel implementation.

Appendix A

Basic concepts of fluid mechanics In this Chapter we present the mathematical foundation of our work, leaving all the physical characterization in the following chapter. Starting from the first principles of mass, momentum and energy conservation, we derive the fluid mechanics fundamental equations in presence of an arbitrary number of sources and sinks in full generality. We then examine the total angular momentum conservation, its link with the momentum conservation and its implications on the symmetry of the stress tensor. Afterwards, a treatment of the viscosity for a Newtonian fluid is presented. Finally, we derive the general source and sink terms in the case of an elliptical galaxy.

114

A.1

Basic concepts of fluid mechanics

From the first principles of conservation to the basic equations of fluid mechanics

In this Section we derive the continuity, the momentum and energy equations from the first principles of mass, momentum and energy conservation, applied in their integral form. We also examine the total angular momentum conservation, its links with the momentum equation and its implication on the symmetry of the stress tensor of a fluid. To derive these basic equations of fluid dynamics we will employ the fundamental theorem of continuum mechanics, the Reynolds Transport Theorem, which describes the time evolution of a fluid quantity integrated on a material volume. The latter is a volume formed by the same fluid particles, so it is comoving with the fluid and at the same time is dragged and deformed by the fluid flow. Let V (t) be a material volume; integrating a quantity F (x, t) over V (t), the Reynolds Transport Theorem states that the time variation of this volume integral is equal to the net flux over the closed surface bounding V (t), plus the Eulerian change of F integrated over the material volume V (t). The quantity F can be a scalar, vector or tensor field, so that this is a kinematical result of wide application. In this circumstance we only enunciate the theorem, the proof is presented in Aris (1989). Theorem A.1 (Reynolds Transport Theorem). Let F (x, t) be a general scalar, vector or tensor field; V (t) a material volume, S(t) the close surface bounding it and ˆ the outward normal to S(t). The theorem states that: n   Z Z d DF F (x, t) dV = + F (∇· u) dV dt V (t) Dt V (t) Z I (A.1) ∂F ˆ dS. = dV + Fu·n V (t) ∂t S(t) Note the use of the d/dt sign (instead of ∂/∂t or D/ Dt), since an integration in space gives a function of time only. From this fundamental theorem we will derive the basic equation of fluid dynamics.

A.1.1

Mass conservation

Starting from the conservation of mass for a material volume V (t), we derive now the continuity equation. Denoting with ρ(x, t) the fluid density field, the total mass M of a general volume V is: Z M= ρ(x, t) dV. (A.2) V

If V is a material volume (usually denoted as V (t)), the mass rate of change of this volume is equal to Z dM d = ρ(x, t) dV. (A.3) dt dt V (t) Applying the Reynolds Theorem with F = ρ, we have   Z dM Dρ = + ρ(∇· u) dV. dt V (t) Dt

(A.4)

A.1 From the first principles of conservation to the basic equations of fluid mechanics 115

We assume now the presence of n sources or sinks of mass, each of them denoted with the term ρ˙ l , with l ranging from 1 to n. Of course, if the l-th field is a source (sink) of mass, ρ˙ l is defined positive (negative). With this additional assumption, equation (A.4) becomes:   Z Z n X Dρ + ρ(∇· u) dV = ρ˙ l dV, (A.5) V (t) Dt V (t) l=1

and since this equality is true for any choice of the material volume V (t), the integrand itself must vanish1 . As a consequence, we can write the continuity equation: X Dρ ∂ρ + ρ ∇· u = + ∇·(ρu) = ρ˙ l , Dt ∂t

(A.6)

where we omitted the index l from the summation symbol for the sake of notation simplicity. Combining the continuity equation with the Reynolds Theorem, we obtain a useful corollary:  Z Z X  DF d ρF dV = ρ (A.7) ρ˙ l dV. +F dt V (t) Dt V (t)

A.1.2

Linear momentum conservation

Here we derive the momentum equation from another first principle: the linear momentum conservation. The procedure is quite similar to the one adopted in the previous Section to derive the continuity equation. We start considering the Newton II law for the usual material volume V (t): Z d ρu dV = F , (A.8) dt V (t) where the l.h.s. is the total pulse rate of change which, using the Reynolds Transport Theorem and the continuity equation can be expressed as:  Z Z X  d Du ρu dV = ρ +u ρ˙ l dV. (A.9) dt V (t) Dt V (t) The term F present in the r.h.s. of equation (A.8) represents the sum of the forces acting on the material volume plus the net momentum injections. There are three type of forces acting on a fluid: body forces, surface forces and line forces. We shortly describe them here: 1

The continuity equation expresses the mass conservation in its strong form. Indeed, the weak and the strong form of a conservation law are respectively the integral and the differential form. This nomenclature arises from the fact that integral form admits the existence of solutions represented by discontinuous functions, whose derivatives do not exist in the classical sense of the term, but only in the sense of distribution theory. Therefore a derivative of one of these solutions is a distribution rather than a function. For this reason, a solution represented by a function whose derivative is a distribution is called a weak solution; since only the integral form admits this kind of solutions, the integral form is also denoted as the weak form. On the other hand, the differential form does not admit solutions other than the differentiable ones. This is a stronger condition compared to the integrability requested by the weak form, therefore the differential form is also denoted as the strong form.

116

Basic concepts of fluid mechanics

• Body forces: they arise from the immersion of the fluid into a force field (e.g. the gravitational field). They are distributed throughout the entire fluid and they are proportional to its mass. Since the body forces are generated by the interaction between the fluid and the force field (which is not a material object), no physical contact is needed in order to apply a body force on the fluid. We will denote a generic body force with the vector field g, having the dimensions of a specific force. Body forces can be conservative; in this case g = −∇Φ, and we call Φ the force potential. The total force applied by the field g on a generic volume V can be formulated as: Z ρg dV, (A.10) V

where V might not be a material volume. • Surfaces forces: they act on a surface through direct contact. They are proportional to the extent of the area and the resulting force on a generic surface S can be expressed as: Z t(n) dS, (A.11) S

where t(n) is the surface stress vector (force per unit area). This field can be represent in terms of the stress tensor σij in this way2 : t(n)i = σij nj ,

(A.12)

where nj is the j-th component of the outward normal to the surface (denoted ˆ therefore this component lie on the ej direction3 . as n), The stress tensor is a second order tensor, it depends only on the physical features of the fluid and it can be symmetric or antisymmetric, it depends on the evolution of the moment of linear momentum. In Section A.1.4 a detailed description of the link between the symmetry of σij and the conservation of angular momentum is given. Now we look for a moment at the product σij nj (here considered without the Einstein convention), present in equation (A.12). This term represents the force per unit area in the ei direction acting on a surface normal to the ej direction. When i = j, the force acts perpendicularly to the surface, while if i 6= j the force acts parallel to the surface. A force acting parallel to any surface is called a shear force. Inserting the stress tensor in equation (A.11), we have the net force produced by surface forces, expressed as a tensorial flux across the surface S 4 : Z Z t(n)i dS = σij nj dS. (A.13) S 2

S

Herefter, when not explicitly declared, the Einstein convention is present. In the entire Thesis we will often switch between the vector and the indicial notation of a physical quantity. The former is of course valid for every coordinate system, while for the latter we explicitly use Cartesian coordinates. Since our space of work is the standard Euclidean three-dimensional space, we can work in Cartesian coordinates without loss of generality, as long as the final expression is writable in vector notation. 4 Hereafter, where not specified, the Einstein summation convention is present. 3

A.1 From the first principles of conservation to the basic equations of fluid mechanics 117

If the surface S is closed, through the Gauss’ theorem5 we can express the r.h.s. of the last equation as the divergence of the stress tensor, integrated over the volume V enclosed by S, obtaining: I I Z ∂σij t(n)i dS = σij nj dS = dV. (A.14) S S V ∂xj • Line forces: they are generated by surface tension forces, they act along a line and have an intensity proportional to the extent of the line. These forces are present between two immiscible liquids, or between a gas and a liquid. They do not appear directly in the equation of momentum, but only in its boundary conditions. Inserting the Equations A.9, A.10 and A.14 in the momentum conservation A.8, considering the net injection (sink) of momentum, and imposing the resulting relation holds for any arbitrary material volume V (t), we obtain (in vector notation) ρ

 X Du = ρg + ∇· σ + P˙ l − ρ˙ l u , Dt

(A.15)

where P˙ l is the momentum of the fluid when it is injected (subtracted) by the l-th source (sink); it is defined with the same sign of ρ˙ l . At this stage of the discussion, we must specify the nature of the fluid we want to treat and what kind of physics we want to consider (e.g. viscosity, presence of magnetic fields). In this work we deal with a Newtonian, perfect, viscous fluid, in the presence of a conservative body force field represented by a potential Φ. We neglect the effects of magnetic field and thermal conduction. For such a fluid, the stress tensor reduces to σ = −pI + T , where p is the thermal pressure of the classical thermodynamics, I is the identity tensor and T is the viscous stress tensor. In addition, the assumption of a Newtonian fluid implies for the viscous stress tensor a very specific analytic form, as a function of the velocity field. This specific relation will be presented in Section A.1.5; until then we will denote the viscous stress tensor simply as T . Then, putting the above assumptions into equation (A.15), the latter reduces to the well-known form of the momentum equation in presence of source and sink terms:  X Du ρ = −∇p − ρ∇Φ + ∇· T + P˙ l − ρ˙ l u . (A.16) Dt fluid.

A.1.3

Energy conservation

We begin considering the usual material volume V (t), and let e(x, t) and E(x, t) = ρe be the internal energy per unit mass and the internal energy per unit volume. The 5 Let V be a volume, S the close surface bounding it and F a scalar, vector or tensor field; the Gauss’ theorem states: Z I ∇· F dV = F dS. V

S

118

Basic concepts of fluid mechanics

Second Law of Thermodynamics states that the rate of change of the fluid kinetic and internal energy is equal to the net work rate done by volume and surface forces on V (t), plus every source or sink of energy: d dt

Z V (t)

ρ kuk2 E+ 2

!

Z

Z

t(n) · u dS

ρg · u dV +

dV =

S(t)

V (t)

+

XZ

(El + Kl ) dV −

Z

L dV, (A.17)

V (t)

V (t)

where El and Kl are the sources (sinks) of internal and kinetic energy, respectively, and L is the bolometric radiative losses per unit time and volume. By means of the Reynolds Transport Theorem, the l.h.s. can be expressed as # " Z DE ρ kuk2 D(ρ kuk2 ) dV, (A.18) + E ∇· u + ∇· u + Dt 2 2 Dt V (t) and by inserting equations A.6-A.15 in it, the integrand can be recasted as   X DE ρ˙ l ˙ + E ∇· u + u · ρg + ∇· σ + (A.19) Pl − u . Dt 2 The r.h.s. of equation (A.17) requires a little of math. By applying the Gauss’ Theorem on the second term and expressing it in indicial notation, as in equation (A.14), we have I I Z ∂(σij ui ) t(n)i ui dS = σij ui nj dS = dV ∂xj S S V (t)   (A.20) Z ∂σij ∂ui + σij , = ui ∂xj ∂xj V (t) which can be represented in vector notation6 as the integral of u · (∇· σ) + σ : ∇u. We can now express the energy conservation by collecting equation (A.19), the r.h.s. of equation (A.17) and equation (A.20), obtaining: X DE + (E + p) ∇· u = T : ∇u + u · Dt



ρ˙ l u − P˙ l 2

where we performed the substitution σ = −pI + T .

A.1.4

 +

X

(El + Kl ) − L , (A.21)

Angular momentum conservation

In the previous Section we introduced the stress tensor σij , and we pointed out that it can be symmetric or antisymmetric. This feature of the stress tensor basically depends on two factors: the conservation of total angular momentum and the type of source of angular momentum we want to treat. 6 Hereafter we use the double dot product A : B between two tensors A, B to denote the scalar quantity Aij Bij .

A.1 From the first principles of conservation to the basic equations of fluid mechanics 119

If a fluid is such that the torques arise only from macroscopic forces (i.e. they are moments of body and surface forces), it is denoted as nonpolar. For this class of fluids we will prove that the conservation of the moment of linear momentum per unit mass, x × u, implies the symmetry of the stress tensor. On the other hand, a fluid having intrinsic body torques and stress couples in addition to the moments of body and surface forces is called a polar fluid. The total angular momentum of this kind of fluid is composed by the moment of linear momentum per unit mass, x × u, plus an intrinsic angular momentum per unit mass, l. We will not treat here the case of a polar fluid. We examine now the conservation of total angular momentum for a nonpolar fluid. In this case the total angular momentum per unit mass reduces to x × u (the moment of linear momentum per unit mass), thus its conservation is expressed by the following equation: Z I Z   X d x × ρg + ρ˙ l wl dV + x × t(n) dS. (A.22) ρx × u dV = dt V (t) V (t) S(t) The l.h.s. can be written as: Z

 X  Du x× ρ +u ρ˙ l dV, Dt V (t)

(A.23)

where we have used corollary A.7 and the identity u = Dx/ Dt. Equation of angular momentum conservation can be written in indicial notation as:  Z Z   X  X Duk ρ˙ l dV = εijk xj ρgk + ρ˙ l wk,l dV εijk xj ρ + uk Dt V (t) V (t) I + εijk xj σkl nl dS, (A.24) S(t)

where, as usual, V (t) is a material volume, S(t) is the close surface bounding it, σkl is the stress tensor and εijk is the Levi-Civita symbol. Through Gauss’ theorem, the last term of equation (A.24) is: I Z ∂ εijk xj σkl nl dS = (εijk xj σkl ) dV S(t) V (t) ∂xl (A.25)   Z ∂σkl = εijk xj + σkj dV. ∂xl V (t) Substituting the previous term in equation (A.24) and rearranging it, the conservation of angular momentum assumes this form:   Duk ∂σkl X εijk xj ρ − ρgk − − ρ˙ l (wk,l − uk ) = εijk σkj . (A.26) Dt ∂xl The l.h.s. of this equation vanishes by means of equation (A.15): this implies that εijk σkj = 0. Noting that the components of εijk σkj are (σ23 − σ32 ), (σ31 − σ13 ), (σ12 − σ21 ), thus εijk σkj = 0 implies that σkj = σjk , therefore the stress tensor must be symmetric.

120

Basic concepts of fluid mechanics

At this point, we can make some considerations regarding the link between the symmetry of the stress tensor and the conservation of angular momentum for a nonpolar fluid. If the stress tensor is general (i.e. it is not specified by the physics of the problem), the conservation of angular momentum states that it must be symmetric. On the other hand, if σ is known and it is symmetric (e.g. σ = −pI), the equation of angular momentum conservation is reduced to the moment of the Cauchy’s equation. As a consequence, for a nonpolar fluid the symmetry of σ is a necessary and sufficient condition for the conservation of angular momentum.

A.1.5

Viscosity

In this Section we deal with the mathematical formulation of the viscosity in the case of a Newtonian fluid, especially for what concerns the links between the viscous stress tensor and the kinematical quantities of a fluid. Indeed, in the last Section A.1.2 we made specific assumption on the nature of the fluid we want to treat, without specify the form of the viscous stress tensor T . In this kind of analysis, it is useful to examine the velocity gradient of a fluid ∇u, which can decomposed in two different parts as     ∂uj ∂uj 1 ∂ui 1 ∂ui ∂ui = + − + , (A.27) ∂xj 2 ∂xj ∂xi 2 ∂xj ∂xi where the first term is second order symmetric tensor named strain tensor e, while ˆk with i = the second term is clearly antisymmetric, and it is equal to (∇×u)· e 6 k 6= j. A necessary and sufficient condition for the conservation of the angular momentum is the symmetry of stress tensor (see Section A.1.4). The hypothesis of a Newtonian fluid implies a much stronger condition: the viscous stress tensor must be directly proportional to the strain e: Tij = Aijkl ekl and the fluid must be isotropic, i.e. the fluid properties are the same in all directions (Aris 1989). As a consequence the Aijkl must be symmetric in i and j, and must be an isotropic fourth order tensor. The most general form that satisfies these requirements is Aijkl = ηδik δjl + αδil δjk + βδij δlk ,

(A.28)

as a consequence the viscous stress tensor becomes Tij = ηeij + αeij + βδij ekk = 2ηeij + βδij ekk  ekk  δij + ζekk δij , = 2η eij − 3

(A.29)

where in the second equality we take advantage of the symmetry of T by imposing α = η. The last expression shows the viscous stress tensor decomposed in a traceless tensor, called the shear tensor, and a spherical tensor having as a trace the divergence of the velocity field (where we defined ζ ≡ β + 2η/3). The shear tensor describes the rate at which the fluid is deformed by shearing, ignoring any changes of volume, while the last tensor describes the stress (i.e. dissipation) when a fluid is expanded or compressed. By explicitly showing the dependence of the strain tensor from the velocity we obtain  h i  2 T T = η ∇u + (∇u) + ζ − η (∇· u) I, (A.30) 3

A.2 Source and sink terms inside an elliptical galaxy

121

where η is called (shear or kinematical) viscosity, while ζ is usually denoted as bulk viscosity (or second viscosity, see Aris 1989). The importance and the value of the bulk viscosity is still debated today, since direct measurements are difficult to carry out, however it has been verified that for a common gas the bulk viscosity is negligible, although important in some physical processes as the sound absorption in fluids (Landau & Lifshitz 1959). On the theoretical ground, Stokes claimed that for a generic Newtonian fluid ζ = 0 (Stokes’ hypothesis), and this assumption has been proved to be true for a low density monoatomic gas (for a review on the importance of the bulk viscosity, see Gad-el-Hak 1995). In the case of this work, the viscosity is not directly related to the fluid micro-physics, but instead we aim to use its formalism to create a sort of “gravitational viscosity”, that will mimic a three-dimensional disc instability (see Section 5.2.3). For simplicity, we therefore accept the Stokes’ hypothesis, and we treat the viscous stress tensor as h i 2 T = η ∇u + (∇u)T − η (∇· u) I. (A.31) 3

A.2

Source and sink terms inside an elliptical galaxy

In the previous Section we derived the fundamental equations for a viscous Newtonian fluid with an arbitrary number of general source and sink terms, without specify neither the nature of the source/sink fields, nor their physical parametrization. However, the goal of this Thesis is to simulate the evolution of an elliptical galaxy over cosmic times, so as to we have to derive appropriate macroscopic source and sink terms to account for various physical phenomena, e.g. the shed of mass, momentum and energy from an ageing stellar population. This task is far from being simple, even if we simply consider the mass injection process only. Indeed, consider a single star, i.e. point-like source, moving through the ISM, which in turns has a non-vanishing velocity field. The way this source injects mass, by shedding material in every direction, is complex and can be anisotropic7 , for example in fast-rotating massive O and B stars. In addition, the single star ejects also momentum, internal (due to the non-vanishing temperature of the ejected mass) and kinetic energy (the stellar ejecta pushes away the in-situ ISM) that will be thermalized in part. To sum up, since an elliptical galaxies contains one hundred billion stars, it is not immediate to write a macroscopic, continuous field ρ˙ that accounts for all this phenomena. Luckily, we are in possession of a powerful tool that is the distribution function f (x, w, t), which describes at any time the dynamical state of every star in a stellar population. Indeed, it provides the number of stars dN with a position in volume d3 x centred on x and with a velocity included the range w + d3 w and w: dN = f (x, w, t) d3 x d3 w.

(A.32)

The basic idea of this procedure is to take the microscopic source terms (valid for every stars) and averaging them of the velocity space, thus obtaining their macroscopic counterparts, i.e. associated not to a single star but to the entire stellar population. 7 Even if this injection was isotropic with respect to the star itself, its effects would be anisotropic in the ISM, due to the relative velocity between the stars and the ISM.

122

Basic concepts of fluid mechanics

We treat here the the simple case of an ageing single stellar population that injects mass, momentum and energy by means of processes as Type Ia supernovae and stellar winds from red giant stars. In fact, in this case the distribution function is the same for every source terms, and the number of microscopic sources does not change with time. On the other hand, in the case of source fields related to a forming stellar population, e.g. Type II supernovae, the number of sources varies continuously (and so the distribution function) due to the continuous formation and death of high-mass stars. We treat such a special case in Section 2.2. The detailed calculation is presented in Posacki (2011), here we briefly recall the main points.

A.2.1

From microscopic to macroscopic source terms

From the definition of the distribution function it follows that the number density of the stellar sources is: Z n(x, t) = f (x, w, t) d3 w. (A.33) R3

We introduce now a set of functions, represented by m = m(x, w, t, n),

(A.34)

p = m [w + us (x, w, t, n)n] ,

(A.35)

einj = einj (x, w, t, n), 1 k = m kw + us (x, w, t, n)nk2 , 2

(A.36) (A.37)

that parametrize the return in the ISM of mass, momentum, specific internal and kinetic energy rates per unit solid angle, respectively, where n is a unit vector and us is the velocity of the ejected material. Note that every function depends not only on the position in the phase-space, but also on the direction towards the material is ejected, since in general the ejection process is not isotropic. We can write the macroscopic source terms as: Z Z ρ(x, ˙ t) = mf d3 w d2 n, (A.38) 3 Ω R Z Z P˙ (x, t) = pf d3 w d2 n, (A.39) 3 Ω R Z Z E (x, t) = meinj f d3 w d2 n, (A.40) 3 Ω R Z Z K (x, t) = kf d3 w d2 n, (A.41) Ω

R3

where we performed the integration over the velocity space and the whole solid angle. Equations (A.38)-(A.41) exhibit the macroscopic mass, momentum, internal and kinetic energy source terms for a general source field, so that they can be inserted directly in the equations of fluid dynamics. However, these terms are still quite general. Indeed, they represent a non-isotropic source field, whose injection functions (m, p, us , e, k) depend on the star velocity v. By consider the special case in which

A.2 Source and sink terms inside an elliptical galaxy

123

the injection functions are isotropic and do not depend on the source velocity, the above equations can be written8 as (A.42)

ρ(x, ˙ t) = 4πmn, P˙ (x, t) = ρv, ˙

(A.43)

E (x, t) = ρe ˙ inj , i 1 h K (x, t) = ρ˙ kvk2 + Tr(σ 2 ) + u2s , 2

(A.44) (A.45)

where v and σ 2 are the source streaming velocity and is the velocity dispersion tensor, defined as Z Z 1 1 3 2 v(x, t) = wf d w, σij = (wi −vi )(wj −vi )f d3 w; (A.46) n(x, t) R3 n(x, t) R3 and it follows that 1 Tr(σ ) + kvk = n(x, t) 2

2

Z

kwk2 f d3 w.

R3

(A.47)

Finally, we can insert equations (A.42)-(A.45) into eqs. A.6, A.16 and A.21, recovering: X Dρ + ρ ∇· u = ρ˙ l , Dt X Du ρ = −∇p − ρ∇Φ + ∇· T + ρ˙ l (v l − u) , Dt X DE + (E + p) ∇· u =T : ∇u + (ρ˙ l einj,l ) − L Dt i X ρ˙ l h + kv l − uk2 + Tr(σ 2l ) + u2sl . 2

(A.48) (A.49) (A.50)

Equations (A.48)-(A.50) are the governing equations of a Newtonian viscous fluid in presence of an arbitrary number of isotropic source/sink terms. It is apparent that every source (or sink) must be physically characterized, by specifying the rate of injection or subtraction of mass (ρ˙ l ) as a function of position and time.

8

Note that the integral of n over the whole solid angle vanishes.

Appendix B

Vector operators in cylindrical and spherical coordinates In this Appendix we present the vector basis of the cylindrical and spherical reference system and its relations with the standard Cartesian basis. Afterwards, we recall the general form of the most common differential vector operators expressed in cylindrical coordinates.

126

B.1

Vector operators in cylindrical and spherical coordinates

The Cartesian, cylindrical and spherical coordinate systems

In this Appendix we illustrate three types of coordinate systems, the Cartesian, the cylindrical and the spherical ones, distinguishing between the coordinates of a particular coordinate system and its vector basis (which is composed by the unit vectors). In the case of the Cartesian system the coordinates are indicated by (x, y, z), ˆx , while the unit vectors constituting the Cartesian basis are usually denoted as e ˆy , e ˆz . As a consequence the position vector is expressed as x = xˆ e ex + yˆ ey + zˆ ez . The main feature of the Cartesian reference system is the independence of its unit vectors from the position in space, therefore any spatial derivative of the Cartesian unit vectors vanishes. On the other hand, in the cylindrical coordinate system the vector basis is ˆR , e ˆϕ , e ˆz , while the coordinates are denoted as (R, ϕ, z), respectively formed by e the cylindrical radius, the azimuth and the vertical coordinate. The position vector in this coordinate system is x = Rˆ eR + zˆ ez . ˆr , Finally, in the vector basis of the spherical coordinate system is denoted as e ˆθ , e ˆϕ , whereas the coordinates are indicated as (r, θ, ϕ), called the spherical radius, e the polar angle and the azimuth (note that the azimuth in cylindrical and spherical systems is the same). The cylindrical and spherical systems are called curvilinear systems because their vector bases depend on the position, unlike the Cartesian one. In fact, if we express the unit vectors of the cylindrical and spherical reference system as a function of the Cartesian one, we have:  ˆR = e ˆx cos ϕ + e ˆy sin ϕ,  e (B.1) ˆϕ = −ˆ ˆy cos ϕ, e ex sin ϕ + e   ˆz = e ˆz , e  ˆr  e ˆθ e   ˆϕ e

ˆx sin θ cos ϕ + e ˆy sin θ sin ϕ + e ˆz cos θ, =e ˆx cos θ cos ϕ + e ˆy cos θ sin ϕ − e ˆz sin θ, =e ˆy cos ϕ. = −ˆ ex sin ϕ + e

(B.2)

It is clear that the unit vectors depend on the azimuthal angle ϕ and on the polar angle θ (for the spherical basis only), so their derivatives are:  ∂ˆ eR   ˆϕ ,  ∂ϕ = e (B.3)  ∂ˆ e   ϕ = −ˆ eR . ∂ϕ  ∂ˆ er   ˆϕ sin θ, =e   ∂ϕ    ∂ˆ eθ (B.4) = −ˆ eϕ cos θ,  ∂ϕ     ˆθ  e ˆr . =e ∂θ

B.2 Differential operators in cylindrical and spherical coordinates

127

Of course the cylindrical and the spherical coordinate system are right-handed, so that the external product between the unit vectors gives:  ˆR × e ˆϕ  e ˆϕ × e ˆz e   ˆz × e ˆR e

ˆz , =e ˆR , =e ˆϕ . =e

(B.5)

 ˆr × e ˆθ  e ˆϕ × e ˆr e   ˆθ × e ˆϕ e

ˆϕ , =e ˆθ , =e ˆr . =e

(B.6)

A general treatise of the curvilinear coordinate systems can be found in Aris (1989), whereas Arfken & Weber (2005) present a lot of examples of curvilinear coordinate systems, among which the cylindrical and the spherical ones.

B.1.1

Coordinates and vector transformations

It is useful deriving the coordinates and vector transformations between the cylindrical and the spherical basis. We treat here the case of vector expressed in terms of the spherical basis but with cylindrical components. Indeed, let A be a generic vector, ˆi (where i can be every coordinate of any coordinate and by defining Ai ≡ A · e system) it can be expressed in the cylindrical and spherical system as ˆ R + Aϕ e ˆ ϕ + Az e ˆz = Ar e ˆr + Aθ e ˆ θ + Aϕ e ˆϕ . A = AR e

(B.7)

In order to express Ar and Aθ in terms of the cylindrical components AR and Az , ˆr and e ˆθ , obtaining: we multiply equation B.7 by e ( Ar = AR sin θ + Az cos θ, Aθ = AR cos θ − Az sin θ; while the coordinates transformations formulae are: ( ( √ R = r sin θ, r = R2 + z 2 , z = r cos θ, θ = atan2(R, z),

(B.8)

(B.9)

where atan2 is the arctangent function with two arguments.

B.2

Differential operators in cylindrical and spherical coordinates

We present now a list of the most common differential operators, expressed in the cylindrical and spherical coordinate system. The operators listed below are, respectively, the gradient, the divergence, the laplacian of a scalar and the curl operators.

128

Vector operators in cylindrical and spherical coordinates

Here, f indicates a generic scalar field while A represents a generic vector field. We first present the operators in the cylindrical system: 1 ∂f ∂f ∂f ˆϕ ˆz , +e +e ∂R R ∂ϕ ∂z

(B.10)

1 ∂(RAR ) 1 ∂Aϕ ∂Az + + , R ∂R R ∂ϕ ∂z

(B.11)

ˆR ∇f = e ∇· A =

1 ∂ ∆f = R ∂R  ∇×A=

∂Aϕ ∂Az − ∂ϕ ∂z



  ∂f 1 ∂2 f ∂2 f R + 2 2 + 2, ∂R R ∂ϕ ∂z

 ˆR + e

∂AR ∂Az − ∂z ∂R

 ˆϕ e 

+

(B.12)

1 ∂(RAϕ ) 1 ∂AR − R ∂R R ∂ϕ



ˆz , (B.13) e

whereas in the spherical system the assume the following form: ∇f =

∇· A =

1 ∂f 1 ∂f ∂f ˆr + ˆθ + ˆϕ , e e e ∂r r ∂θ r sin θ ∂ϕ

1 ∂(sin θAθ ) 1 ∂Aϕ 1 ∂(r2 Ar ) + + , 2 r ∂r r sin θ ∂θ r sin θ ∂ϕ

1 ∂ ∆f = 2 r ∂r

  1 ∂2 f ∂2 f 2 ∂f r + 2 + 2, ∂r r sin θ ∂θ ∂z

(B.14) (B.15) (B.16)

    ∂(rAϕ ) 1 ∂(Aθ sin θ) ∂Aθ 1 1 ∂Ar ˆr + ˆθ ∇×A= − e − e r sin θ ∂θ ∂ϕ r sin θ ∂ϕ ∂r   1 ∂(rAθ ) ∂Ar ˆϕ . (B.17) + − e r ∂r ∂θ

Appendix C

Viscosity module implementation In this Appendix we describe the implementation details of the viscosity terms presented in Section 5.2.3.

130

Viscosity module implementation

C.1

Discretized viscosity equations

The new viscosity module has been implemented as a separate sub-step of the momentum and energy equations. While equations (5.4)-(5.5) show the explicit analytical formulae in the case of spherical and cylindrical coordinates, in the code the equations are implemented following the covariant formalism, so that a single equation can handle different coordinate systems (Stone & Norman 1992; Hayes et al. 2006). We present here the covariant version of equations (5.4)-(5.5), discretized on the ZEUS numerical grid (see Figure 2.2). From Stone & Norman (1992) and the physical assumptions above, we have T13 = η [(∇u)13 + (∇u)31 ] ,

T23 = η [(∇u)23 + (∇u)32 ] ,

T : ∇u = η [(∇u)13 + (∇u)31 ]2 + η [(∇u)23 + (∇u)32 ]2 , (∇· T )(ϕ) = g31

2T32 ∂g32 ∂ (g2 T31 ) g32 ∂T32 2T31 ∂g31 + + + , ∂V1 g2 ∂V2 g31 ∂x1 g2 g32 ∂x2

(C.1) (C.2) (C.3)

where we indicated with index numbers the components of a tensor with respect to the ZEUS coordinates, see Table 2.1. By writing the necessary ∇u components in the ZEUS covariant formalism, we obtain (∇u)13 = g2 g31 (∇u)23 =

∂u3 , ∂V1

g32 ∂u3 , g2 ∂V2

(∇u)31 = −

u3 ∂g31 , g31 ∂x1

(C.4)

(∇u)31 = −

u3 ∂g32 . g2 g32 ∂x2

(C.5)

By indicating the centring of a variable on the (x1, x2) grid as a superscript, e.g. (∇u)ab 32i,j,k is centred on (x1ai , x2bj ), we discretize the viscous stress tensor, given by equation (C.1), on the numerical grid1 : h i bb bb bb T31 = η (∇u) + (∇u) (C.6) i,j,k 31i,j,k 13i,j,k , i,j,k ab T31 = i,j,k

i ηi,j,k + ηi−1,j,k h ab (∇u)ab 31i,j,k + (∇u)13i,j,k , 2

h i bb bb bb T32 = η (∇u) + (∇u) i,j,k 32i,j,k 23i,j,k , i,j,k ba T32 = i,j,k

i ηi,j,k + ηi,j−1,k h ba (∇u)ba + (∇u) 32i,j,k 23i,j,k . 2

(C.7)

(C.8) (C.9)

Equation (C.2) can be immediately written as h i2 h i2 bb bb bb bb (T : ∇u)bb = η (∇u) + (∇u) + η (∇u) + (∇u) , i,j,k i,j,k i,j,k 13i,j,k 31i,j,k 23i,j,k 32i,j,k (C.10) 1

As every scalar field, η is centred on (x1bi , x2bj ).

C.1 Discretized viscosity equations

131

whereas equation (C.3) becomes (∇· T )bb (ϕ)i,j,k = g31bi

ab ab g2ai+1 T31 − g2ai T31 i+1,j,k i,j,k

dV 1ai +

bb 2T31 i,j,k

g31bi

+

ba ba g32bj T32i,j+1,k − T32i,j,k g2bi dV 2aj

dg31bd1i +

bb 2T32 i,j,k

g2bi g32bj

dg32bd2j . (C.11)

In the expressions above we used the following quantities: (∇u)bb 31i,j,k = −

u3i,j,k dg31bd1i , g31bi

(C.12)

(∇u)bb 32i,j,k = −

u3i,j,k dg32bd2j , g2bi g32bj

(C.13)

(∇u)bb 13i,j,k

g2bi g31bi = 2

(∇u)bb 23i,j,k

g32bj = 2g2bi





u3i,j,k − u3i,j−1,k u3i,j+1,k − u3i,j,k + dV 2bj dV 1bj+1

(∇u)ab 31i,j,k = −dg31ad1i (∇u)ab 13i,j,k = g2ai g31ai

u3i,j,k + u3i−1,j,k , 2g31ai

u3i,j,k − u3i−1,j,k , dV 1bi

(∇u)ba 32i,j,k = −dg32ad2j (∇u)ba 23i,j,k =

u3i+1,j,k − u3i,j,k u3i,j,k − u3i−1,j,k + dV 1bi dV 1bi+1

u3i,j,k + u3i,j−1,k , 2g2bi g32aj

g32aj u3i,j,k − u3i,j−1,k . g2bi dV 2bj

 ,

 ,

(C.14) (C.15) (C.16) (C.17) (C.18) (C.19)

Given the explicit nature of the ZEUS code, the inclusion of the fluid viscosity substep requires the additional time-scale ∆tv = ρη −1 min(dV1 , dV2 )2 in equation (2.43), where the volume elements are listed in Table 2.1.

Bibliography Arfken G. B., Weber H. J., 2005, Mathematical methods for physicists 6th ed., Arfken, G. B. & Weber, H. J., ed. Aris R., 1989, Vectors, Tensors and the Basic Equation of Fluid Mechanics. Dover Books on Engineering Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57 Bertin G., Lodato G., 2001, A&A, 370, 342 Blaes O., 2014, Space Sci. Rev., 183, 21 Boroson B., Kim D.-W., Fabbiano G., 2011, ApJ, 729, 12 Brighenti F., Mathews W. G., 1996, ApJ, 470, 747 Cappellari M. et al., 2007, MNRAS, 379, 418 Cappellari M. et al., 2011, MNRAS, 413, 813 Ciotti L., 2000, Lecture notes on stellar dynamics. Pisa, Scuola Normale Superiore Ciotti L., Ostriker J. P., 1997, ApJ, 487, L105 Ciotti L., Ostriker J. P., 2001, ApJ, 551, 131 Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038 Ciotti L., Ostriker J. P., 2012, in Astrophysics and Space Science Library, Vol. 378, Astrophysics and Space Science Library, Kim D.-W., Pellegrini S., eds., p. 83 Ciotti L., Ostriker J. P., Proga D., 2009, ApJ, 699, 89 Ciotti L., Ostriker J. P., Proga D., 2010, ApJ, 717, 708 Ciotti L., Pellegrini S., 1996, MNRAS, 279, 240 (CP96) Coccato L., Morelli L., Pizzella A., Corsini E. M., Buson L. M., Dalla Bontà E., 2013, A&A, 549, A3 Davis T. A. et al., 2013, MNRAS, 429, 534 Davis T. A. et al., 2011, MNRAS, 417, 882 133

134

Bibliography

Davis T. A. et al., 2014, MNRAS, 444, 3427 De Bruyne V., Dejonghe H., Pizzella A., Bernardi M., Zeilinger W. W., 2001, ApJ, 546, 903 de Vaucouleurs G., 1948, Annales d’Astrophysique, 11, 247 de Vaucouleurs G., de Vaucouleurs A., Corwin, Jr. H. G., Buta R. J., Paturel G., Fouqué P., 1991, Third Reference Catalogue of Bright Galaxies D’Ercole A., Ciotti L., 1998, ApJ, 494, 535 (DC98) Desroches L.-B., Quataert E., Ma C.-P., West A. A., 2007, MNRAS, 377, 402 Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87 Emsellem E. et al., 2011, MNRAS, 414, 888 Emsellem E. et al., 2007, MNRAS, 379, 401 Erwin P., Sparke L. S., 2002, AJ, 124, 65 Eskridge P. B., Fabbiano G., Kim D.-W., 1995, ApJS, 97, 141 Fabbiano G., 1989, ARA&A, 27, 87 Fabian A. C., 2012, ARA&A, 50, 455 Ferrarese L., Merritt D., 2000, ApJ, 539, L9 Foster C., Arnold J. A., Forbes D. A., Pastorello N., Romanowsky A. J., Spitler L. R., Strader J., Brodie J. P., 2013, arXiv:1308.3531 Gad-el-Hak M., 1995, Journal of Fluids Engineering, 117, 3 Gan Z., Yuan F., Ostriker J. P., Ciotti L., Novak G. S., 2014, ApJ, 789, 150 Gao L., Navarro J. F., Cole S., Frenk C. S., White S. D. M., Springel V., Jenkins A., Neto A. F., 2008, MNRAS, 387, 536 Gebhardt K. et al., 2000, ApJ, 539, L13 Greggio L., 2005, A&A, 441, 1055 Greggio L., 2010, MNRAS, 406, 22 Gültekin K. et al., 2009, ApJ, 698, 198 Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188 Jardel J. R. et al., 2011, ApJ, 739, 21 Jiang Y.-F., Ciotti L., Ostriker J. P., Spitkovsky A., 2010, ApJ, 711, 125 Jog C. J., 1996, MNRAS, 278, 209

Bibliography

135

Jog C. J., 2014, AJ, 147, 132 Katkov I. Y., Sil’chenko O. K., Afanasiev V. L., 2013, ApJ, 769, 105 Kennicutt, Jr. R. C., 1998, ApJ, 498, 541 Kim D.-W., Pellegrini S., 2012, Hot Interstellar Matter in Elliptical Galaxies, Astrophysics and Space Science Library, Vol. 378. Springer Kollmeier J. A. et al., 2006, ApJ, 648, 128 Kormendy J., Ho L. C., 2013, ARA&A, 51, 511 Krajnović D. et al., 2011, MNRAS, 414, 2923 Kuijken K., Fisher D., Merrifield M. R., 1996, MNRAS, 283, 543 Kuzmin G., 1956, AZh, 33, 27 Landau L. D., Lifshitz E. M., 1959, Fluid mechanics Larson R. B., 2011, in IAU Symposium, Vol. 270, Computational Star Formation, Alves J., Elmegreen B. G., Girart J. M., Trimble V., eds., pp. 1–5 Larson R. B., 2012, in American Institute of Physics Conference Series, Vol. 1480, American Institute of Physics Conference Series, Umemura M., Omukai K., eds., pp. 38–42 Li J., Ostriker J., Sunyaev R., 2013, ApJ, 767, 105 Li J.-T., Wang Q. D., Li Z., Chen Y., 2009, ApJ, 706, 693 Li J.-T., Wang Q. D., Li Z., Chen Y., 2011a, ApJ, 737, 41 Li Z. et al., 2011b, ApJ, 730, 84 Li Z., Wang Q. D., Hameed S., 2007, MNRAS, 376, 960 Mannucci F., Della Valle M., Panagia N., Cappellaro E., Cresci G., Maiolino R., Petrosian A., Turatto M., 2005, A&A, 433, 807 Maoz D., Mannucci F., Li W., Filippenko A. V., Della Valle M., Panagia N., 2011, MNRAS, 412, 1508 Maraston C., 2005, MNRAS, 362, 799 Mathews W. G., Brighenti F., 2003, ARA&A, 41, 191 Mellier Y., Mathez G., 1987, A&A, 175, 1 Merritt D., Graham A. W., Moore B., Diemand J., Terzić B., 2006, AJ, 132, 2685 Miyamoto M., Nagai R., 1975, PASJ, 27, 533 Morelli L. et al., 2004, MNRAS, 354, 753

136

Bibliography

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 Navarro J. F. et al., 2004, MNRAS, 349, 1039 Navarro J. F. et al., 2010, MNRAS, 402, 21 Negri A., Ciotti L., Pellegrini S., 2014a, MNRAS, 439, 823 (N14) Negri A., Pellegrini S., Ciotti L., 2013, Mem. Soc. Astron. Italiana, 84, 762 Negri A., Posacki S., Pellegrini S., Ciotti L., 2014b, MNRAS, 445, 1351 Novak G. S., Faber S. M., Dekel A., 2006, ApJ, 637, 96 Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26 Novak G. S., Ostriker J. P., Ciotti L., 2012, MNRAS, 427, 2734 Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642 O’Sullivan E., Forbes D. A., Ponman T. J., 2001, MNRAS, 324, 420 Parriott J. R., Bregman J. N., 2008a, ApJ, 681, 1215 Parriott J. R., Bregman J. N., 2008b, ApJ, 681, 1215 Pellegrini S., 1999a, A&A, 351, 487 Pellegrini S., 1999b, A&A, 351, 487 Pellegrini S., 2005, MNRAS, 364, 169 Pellegrini S., 2011, ApJ, 738, 57 Pellegrini S., 2012, in Hot Interstellar Matter in Elliptical Galaxies, Astrophysics and Space Science Library, Vol. 378, pp. 21–54, Kim, D.–W. and Pellegrini, S., eds Pellegrini S., Held E. V., Ciotti L., 1997, MNRAS, 288, 1 Pellegrini S., Wang J., Fabbiano G., Kim D.-W., Brassington N. J., Gallagher J. S., Trinchieri G., Zezas A., 2012, ApJ, 758, 94 Plummer H. C., 1911, MNRAS, 71, 460 Posacki S., 2011, Gas Flows in Galaxies: Treatment of Sources and Sinks from First Principles, Master Thesis, Bologna University, Unpublished Posacki S., 2014, The dynamics of early-type galaxies as a tool to understand their hot coronae and their IMF, PhD Thesis, Bologna University, Unpublished Posacki S., Pellegrini S., Ciotti L., 2013a, Mem. Soc. Astron. Italiana, 84, 766 Posacki S., Pellegrini S., Ciotti L., 2013b, MNRAS, 433, 2259 (P13) Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing. Cambridge: University Press, 2nd ed.

Bibliography

137

Retana-Montenegro E., van Hese E., Gentile G., Baes M., Frutos-Alfaro F., 2012, A&A, 540, A70 Romeo A. B., Falstad N., 2013, MNRAS, 433, 1389 Russell H. R., McNamara B. R., Edge A. C., Hogan M. T., Main R. A., Vantyghem A. N., 2013, MNRAS, 432, 530 Sarzi M. et al., 2013, MNRAS, 432, 1845 Sarzi M. et al., 2010, MNRAS, 402, 2187 Satoh C., 1980, PASJ, 32, 41 Sazonov S. Y., Ostriker J. P., Ciotti L., Sunyaev R. A., 2005, MNRAS, 358, 168 Schmidt M., 1959, ApJ, 129, 243 Sharon K. et al., 2010, ApJ, 718, 876 Shin M.-S., Ostriker J. P., Ciotti L., 2010, ApJ, 711, 268 Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91 Stone J. M., Norman M. L., 1992, ApJS, 80, 753 Stone J. M., Pringle J. E., Begelman M. C., 1999, MNRAS, 310, 1002 Tang S., Wang Q. D., 2005, ApJ, 628, 205 Tempel E., Tenjes P., 2006, MNRAS, 371, 1269 Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998, ApJ, 500, 95 Toomre A., 1964, ApJ, 139, 1217 Tremaine S. et al., 2002, ApJ, 574, 740 Werner N. et al., 2014, MNRAS, 439, 2291 Young L. M. et al., 2011, MNRAS, 414, 940