Abstract
We study the relation between stellar ages and vertical velocity dispersion (the age–velocity relation, or AVR) in a sample of seven simulated disc galaxies. In our simulations, the shape of the AVR for stars younger than 9 Gyr depends strongly on the merger history at low redshift, with even 1:10–1:15 mergers being able to create jumps in the AVR (although these jumps might not be detectable if the errors on stellar ages are of the order of 30 per cent). For galaxies with a quiescent history at low redshift, we find that the vertical velocity dispersion rises smoothly for ages up to 8–9 Gyr, following a power law with a slope of ∼0.5, similar to what is observed in the solar neighbourhood by the Geneva-Copenhagen Survey. For these galaxies, we show that the slope of the AVR is not imprinted at birth, but is the result of subsequent heating. By contrast, in all our simulations, the oldest stars form a significantly different population, with a high velocity dispersion. These stars are usually born kinematically hot in a turbulent phase of intense mergers at high redshift, and also include some stars accreted from satellites. This maximum in σz is strongly decreased when age errors are included, suggesting that observations can easily miss such a jump with the current accuracy of age measurements.
1 INTRODUCTION
In the solar neighbourhood, the vertical velocity dispersion (σz) of stars is correlated with their age. Stars with an age of 9–10 Gyr have a σz of about 40 km s−1 (or more, depending on the selection criteria, and in particular on the inclusion or exclusion of halo stars, see for instance Casagrande et al. 2011). By contrast, 1 Gyr-old stars have a σz of 12–16 km s−1, and the youngest stars are rather cold: Aumer & Binney (2009) find that the bluest stars have a σz of only 6 km s−1.
The exact shape of the age–velocity relation (or AVR) is however debated, partly because measuring stellar ages is difficult (Soderblom 2010), and partly because the selection criteria differ between different surveys (Haywood et al. 2013): we show in Fig. 1 two examples of observational results giving AVRs with very different shapes. For instance, Quillen & Garnett (2001) propose a rapid increase of σz for ages smaller than 3 Gyr, followed by a plateau corresponding to a saturation of σz at ∼20 km s−1, and a final jump for ages greater than 9 Gyr. Soubiran et al. (2008) find a similar saturation, but occurring for stars older than 5 Gyr. By contrast, others argue in favour of a continuously rising σz with age, following a power law σz ∝ tα, where α is close to 0.5 (e.g. Wielen 1977; Nordström et al. 2004; Holmberg, Nordström & Andersen 2007; Aumer & Binney 2009). However, Seabroke & Gilmore (2007) as well as Aumer & Binney (2009) have shown that the current data do not favour a power-law evolution over a saturation of σz for ages greater than 4–5 Gyr. The saturation after 3 Gyr that was proposed by Quillen & Garnett (2001) from a sample of 189 stars is excluded by the current data.
The AVR observed in the solar neighbourhood: blue triangles from Holmberg et al. (2007), red dots from Soubiran et al. (2008). The green square corresponds to the bluest stars from Aumer & Binney (2009), arbitrarily placed at a very small age. This illustrates the two competing ways of seeing the shape of the AVR: with a saturation after a few Gyr (as measured by Soubiran et al.), or with a smooth increase of σz with age (as measured by Holmberg et al.).
A possibility is that the shape of the AVR is determined by the properties of stars at birth, i.e. stars are born with a decreasing σz with decreasing redshift. This is supported by observations of turbulent gas discs at high redshift (e.g. Förster Schreiber et al. 2009), so that stars are naturally born hotter at higher redshift (note however that typical clumpy galaxies might be too massive to be the progenitors of Milky Way-type galaxies, e.g. van Dokkum et al. 2013; Inoue & Saitoh 2014), and by the idea that turbulence (and thus σz) decreases with redshift because of decreasing gas accretion rates (Bournaud, Elmegreen & Martig 2009; Forbes, Krumholz & Burkert 2012; Bird et al. 2013).
Most theoretical works have assumed that stars were born with a roughly constant σz over the last 8–9 Gyr, and that their present-day σz is the result of gradual heating. Spitzer & Schwarzschild (1951, 1953) first proposed that the source of that heating could be the scattering of stars by giant molecular clouds (GMCs) with masses of 106 M⊙. It was later shown that the number of GMCs in the Milky Way is too small to reproduce the observed heating, and that ratio of vertical to radial velocity dispersions is too high compared to observations (Lacey 1984).
By contrast, spiral arms are very inefficient at vertical heating but can provide some radial heating provided that they are transient and stochastic (Carlberg & Sellwood 1985), or that multiple density waves coexist (Minchev & Quillen 2006). It was thus suggested that heating by a combination of spiral arms and GMCs could match the observed AVR (Carlberg 1987; Jenkins & Binney 1990). However, Hänninen & Flynn (2002) showed that heating by GMCs is not as efficient as previously thought, and instead proposed a combination of heating by GMCs and halo black holes (see also Lacey & Ostriker 1985 for an early study of heating by halo black holes).
Another important mechanism for disc evolution is the radial migration of stars under the action of non-axisymmetric perturbations (e.g. Sellwood & Binney 2002; Minchev & Quillen 2006; Minchev & Famaey 2010). In recent years, the idea has emerged that radial migration could be a source of disc heating: the analytical model of Schönrich & Binney (2009a,b) combines recipes for radial migration and chemical evolution and shows that migration alone could be responsible for the chemo-dynamical trends observed in the Milky Way. As a result, disc heating has been attributed to radial migration in the simulations analysed by Loebman et al. (2011). However, the model by Schönrich & Binney (2009a,b) assumes that the vertical energy is conserved during migration, while simulations have later shown that it is the vertical action that is conserved (Minchev et al. 2012; Solway, Sellwood & Schönrich 2012). The important consequence is that the contribution of radial migration to disc heating is negligible, except in the outer regions where it induces some flaring (Minchev et al. 2012). This result is the consequence of the balance between outward migrating stars, which heat the disc, and inward migrators, which cool it (note that Roškar, Debattista & Loebman 2013 find that migration induces heating but only consider outward migrators).
It thus seems that while internal secular evolution clearly induces disc thickening, it is often difficult to figure out which mechanisms are acting in realistic simulations that include gas, spiral arms, bars and live dark matter haloes (see for instance Saha, Tseng & Taam 2010, for simulations where heating in the inner disc is due to a bar, but where the causes of heating in the outer disc are less clear). In any case, it would seem that radial migration is not a dominant contributor to secular heating.
For real galaxies however a final potentially important source of heating is the interaction with satellite galaxies, mostly in the form of minor mergers (e.g. Quinn, Hernquist & Fullagar 1993; Walker, Mihos & Hernquist 1996; Villalobos & Helmi 2008; Kazantzidis et al. 2009). The amount of heating provided depends not only on the merger ratio, but also on many other parameters like the orbit of the satellite Velazquez & White (1999) or the gas content of the main galaxy (Moster et al. 2010).
House et al. (2011) have shown from a sample of cosmological simulations performed with different techniques that mergers create jumps in the AVR, and thus point towards a quiescent merger history for the Milky Way, for which no such jumps are seen for ages lower than 9 Gyr.
This paper is the second of a series that uses a sample of seven simulated galaxies with different merger histories to understand the effect of mergers on the structure of galactic discs. In Paper I, we present in detail our sample and the analyses we perform. We dissect the discs in so-called mono-age populations (in 500 Myr-wide age bins), and compare their structure to Milky Way observations by Bovy, Rix & Hogg (2012a) and Bovy et al. (2012b,c). In this paper, we focus solely on the AVR. We present briefly our simulations in Section 2, and then discuss the velocity dispersion of young stars in Section 3. In Section 4, we compare the AVR in galaxies with quiescent or active merger histories. In Section 5, we study the origin of the slope of the AVR in quiescent galaxies, and discuss potential numerical limitations in Section 6 before summarizing our results in Section 7.
2 SIMULATIONS AND ANALYSIS
The simulations we use in this paper are the same as in Paper I, to which we refer the reader for a detailed description of the simulation technique. The sample consists of seven galaxies simulated using a zoom-in cosmological re-simulation technique (see Martig et al. 2009, 2012). This technique consists in first running a ‘dark matter only’ cosmological simulation, in extracting merger and accretion histories for a set of haloes, and in re-simulating these histories at higher resolution, starting at z = 5 with a seed galaxy. We follow the evolution of that seed galaxy down to z = 0 by including mergers as well as accretion of gas and dark matter as prescribed by the cosmological simulation.
The re-simulations use a Particle-Mesh code with a grid size of 150 pc, where gas dynamics is modelled with a sticky-particle algorithm. The mass resolution is 1.5 × 104 M⊙ for gas particles, 7.5 × 104 M⊙ for star particles (1.5 × 104 M⊙ for star particles formed from gas during the simulation), and 3 × 105 M⊙ for dark matter particles. For the galaxies studied here, there are typically 2–6 million star particles in the discs at z = 0.
We model star formation following a Schmidt law with an exponent of 1.5 above a volume density threshold of 0.03 M⊙ pc−3 (i.e. 1 H cm−3). We take into account stellar mass loss (see Martig & Bournaud 2010), as well as kinetic feedback from supernovae explosions, where 20 per cent of the energy of the supernovae is distributed to neighbouring gas particles.
The seven chosen galaxies have a total stellar mass at z = 0 of between 3 × 1010 and 2 × 1011 M⊙, and bulge-to-total ratios between 0.04 and 0.5 (see table 1 in Paper I). They were selected to have a range of formation histories as follows:
g37 and g92 have quiescent histories: over their last 9 Gyr of evolution they only have interactions with mass ratios smaller than 1:50.
g47 is also mostly quiescent, but undergoes a 1:15 interaction with a satellite at t ∼8.5–9 Gyr, the satellite does not merge with g47 but significantly disturbs its morphology.
g22 undergoes a 1:10 merger at t = 10 Gyr.
g48 undergoes a 1:4 merger at t = 8 Gyr.
g102 undergoes a 1:15 merger at t = 8.5 Gyr.
g106 undergoes a 1:5 merger at t = 7 Gyr.
Heating index for the four most quiescent galaxies. The heating index is measured as the slope of a linear fit to the logarithm of the AVR for ages between 2 and 7.5 Gyr (6.5 Gyr for g106). We show both the actual value of the index, and the value measured on the AVR when a 20 or 30 per cent age error is first taken into account (the black lines in Fig. 2).
| Name | Actual value | With 20 per cent age error | With 30 per cent age error |
|---|---|---|---|
| g37 | 0.53 ± 0.02 | 0.52 ± 0.02 | 0.46 ± 0.02 |
| g47 | 0.47 ± 0.03 | 0.46 ± 0.02 | 0.39 ± 0.01 |
| g92 | 0.52 ± 0.06 | 0.54 ± 0.04 | 0.49 ± 0.04 |
| g106 | 0.51 ± 0.04 | 0.62 ± 0.06 | 0.61 ± 0.07 |
| Name | Actual value | With 20 per cent age error | With 30 per cent age error |
|---|---|---|---|
| g37 | 0.53 ± 0.02 | 0.52 ± 0.02 | 0.46 ± 0.02 |
| g47 | 0.47 ± 0.03 | 0.46 ± 0.02 | 0.39 ± 0.01 |
| g92 | 0.52 ± 0.06 | 0.54 ± 0.04 | 0.49 ± 0.04 |
| g106 | 0.51 ± 0.04 | 0.62 ± 0.06 | 0.61 ± 0.07 |
Heating index for the four most quiescent galaxies. The heating index is measured as the slope of a linear fit to the logarithm of the AVR for ages between 2 and 7.5 Gyr (6.5 Gyr for g106). We show both the actual value of the index, and the value measured on the AVR when a 20 or 30 per cent age error is first taken into account (the black lines in Fig. 2).
| Name | Actual value | With 20 per cent age error | With 30 per cent age error |
|---|---|---|---|
| g37 | 0.53 ± 0.02 | 0.52 ± 0.02 | 0.46 ± 0.02 |
| g47 | 0.47 ± 0.03 | 0.46 ± 0.02 | 0.39 ± 0.01 |
| g92 | 0.52 ± 0.06 | 0.54 ± 0.04 | 0.49 ± 0.04 |
| g106 | 0.51 ± 0.04 | 0.62 ± 0.06 | 0.61 ± 0.07 |
| Name | Actual value | With 20 per cent age error | With 30 per cent age error |
|---|---|---|---|
| g37 | 0.53 ± 0.02 | 0.52 ± 0.02 | 0.46 ± 0.02 |
| g47 | 0.47 ± 0.03 | 0.46 ± 0.02 | 0.39 ± 0.01 |
| g92 | 0.52 ± 0.06 | 0.54 ± 0.04 | 0.49 ± 0.04 |
| g106 | 0.51 ± 0.04 | 0.62 ± 0.06 | 0.61 ± 0.07 |
The properties of each galaxy are analysed at z = 0. In this paper, we present the AVR for particles in age bins of 500 Myr. We first remove counter-rotating stars to minimize the contamination of our sample by halo stars. Then, we measure the vertical velocity dispersion for each population within a 2-kpc-wide annulus centred at 2 Rd (Rd is the scalelength obtained from an exponential fit to the total stellar mass distribution), and up to a height including 95 per cent of stars of that population (limited to 5 kpc above the disc plane). We estimate σz from the second moment of the velocity distribution; the resulting AVR is shown in Fig. 2 for each of the seven galaxies.
Vertical velocity dispersion as a function of age at 2Rd. The dots represent the actual value of σz, while the lines correspond to σz measured once a Gaussian error of 20 per cent (solid line) or 30 per cent (dashed line) has been first added to the age of each particle. Galaxies in the top panels have quiescent merger histories for z < 1.5, while galaxies in the bottom panels have more active histories. The small vertical lines mark the time of coalescence for the last merger undergone by each galaxy for z < 1.5 (the dashed line for g47 marks the end of a fly-by). We find a smoothly rising AVR for the quiescent galaxies, while mergers imprint jumps in the relation. The age errors do not change the shape of the AVR for quiescent galaxies, but smooth jumps for galaxies with mergers.
Since errors on ages are a significant observational issue, we perform an additional test: we add a Gaussian error on the age of each stellar particle before re-computing the AVR using the blurred ages. We test assumed age errors of 20 and 30 per cent, which are typical of current observational measurements (Holmberg et al. 2007; Soderblom 2010). The resulting AVRs are shown as black lines in Fig. 2.
In the next section, we compare the overall shape of the AVR for quiescent and active galaxies, and discuss the values of σz for young stars.
3 THE SHAPE OF THE AVR
The shape of the AVR is very different for galaxies with different merger histories (Fig. 2). Quiescent galaxies show a smooth rise of the AVR, while mergers leave well-defined jumps. An exception is g106: it undergoes a 1:5 merger, but this merger happens early enough so that it leaves little trace in the AVR. In the following, we thus study g106 together with the quiescent galaxies.
3.1 Quiescent galaxies
Apart from the oldest populations, which have a high σz in all cases, the quiescent galaxies show a very regularly rising σz with age. Interestingly, the fly-by undergone by g47 does not leave any trace in the AVR at a radius of 2Rd, only leaving signatures in the outer disc of g47.
In the Milky Way, the AVR is often fitted by a power law corresponding to σz ∝ tα in intermediate ages range (excluding young unrelaxed stars, and hot old stars). The heating index α is difficult to measure but usually found close to 0.5 (Holmberg et al. 2007; Seabroke & Gilmore 2007), and its value has been used to discriminate between disc heating mechanisms (e.g. Hänninen & Flynn 2002). We similarly fit the simulated AVRs with a power law for ages between 2 and 7.5 Gyr (6.5 Gyr for g106), and find heating indices close to 0.5 (see Table 1). This is a good match to the AVR in the Milky Way, suggesting that the heating in our simulations is reasonable.
When we add an error of 20 or 30 per cent to the age of each stellar particle, for ‘blurred ages’ smaller than 8 Gyr, we do not find significant changes to the shape of the AVR (black lines on Fig. 2). This is because the number density in this age range generally does not vary much and σz increases smoothly with age. The heating indices are mostly unchanged with 20 per cent age errors, and are slightly reduced with 30 per cent age errors, which means that the AVR is slightly flattened (see Table 1). An exception is g106, for which the stars with a true age of ∼8 Gyr (hotter because of the 1:5 merger), are now mixed also with younger stars so that the AVR is overall steeper, even when only 20 per cent age errors are included: the heating index is increased from 0.51 ± 0.04 to 0.62 ± 0.06. With such large error bars, though, g106 might still be consistent with the Milky Way.
In all galaxies, for the oldest stars, age errors introduce significant changes to the AVR: σz is lowered because the truly hot and old stars are now mixed with younger and colder stars, and because these young stars are more numerous. With age errors, the jumps in σz seen for the oldest stars are then considerably smoothened, and might be missed by surveys.
To illustrate that point, we show in Fig. 3 a comparison between the AVR in the solar neighbourhood from the Geneva-Copenhagen Survey (GCS; Nordström et al. 2004; Holmberg et al. 2007) to one of our simulated galaxies. The observed AVR shows a smooth increase of σz with age, with no jump for old stars, and a σz of about 30 km s−1 for the oldest stars. This observation might seem in contradiction with other studies finding σz ∼ 45 km s−1 for the oldest stars (e.g. Quillen & Garnett 2001), or similarly high values for populations with high [α/Fe] ratios (Lee et al. 2011; Bovy et al. 2012c; Haywood et al. 2013). However, the GCS sample is limited to a small volume around the Sun, and contains very few metal-poor stars (typically [Fe/H] <−0.7), which are the stars for which Bovy et al. (2012c) find σz > 40 km s−1. In addition, when selecting only stars close to the disc mid-plane (|z| < 200 pc) in our g92 simulation, we show in Fig. 3 that even if a jump in σz is present for old stars, this jump is totally erased when 30 per cent age errors are included, and the simulated galaxy matches the observed AVR very well. In this simulation, we find that the absence of a jump is both the consequence of the age errors and of the selection of a sample close to the disc plane, dominated by kinematically cold, thin disc stars (in Fig. 2, which includes all stars of a given age, adding age errors is not enough to totally erase the jump).
AVR for g92, compared with Milky Way data from Holmberg et al. (2007 – grey triangles). We also show a power-law fit to the simulated AVR, the resulting exponent is 0.47±0.03 for g92, showing a good match to the observed slope. We selected only stars close to the disc plane to replicate the selection criterion of the GCS. We find that this results only in a small overall number of old stars in the sample, and that, when age errors are added, the jump in σz found for old stars is totally washed out, making the simulation consistent with GCS data.
This shows that in this sense, our simulations are consistent with observations of the solar neighbourhood, and that even if a jump in σz were present in the Milky Way for old stars, it might not be detected by surveys limited to a small volume and with large age errors.
3.2 The effect of mergers
By contrast, the galaxies undergoing a merger show jumps in their AVRs, which appear at ages corresponding to the end of the merger. This is clearly seen for g22 and g102. For g48, where the merger mass ratio is about twice higher than in g22, the jump is not so sharp, and a lot of time is needed for the gas to cool down again: it takes about 2 Gyr for the velocity dispersion to decrease after the galaxies merge.
It is interesting that even the small 1:10–1:15 mergers leave a well-defined signature on the AVR (Fig. 2, bottom row). We also note here that mergers affect disc kinematics at all radii (as also found by Velazquez & White 1999), while the fly-by in g47 only affects outer regions. Whether this is a general feature of mergers versus fly-bys would need to be verified with larger samples of simulations.
These results are consistent with those of House et al. (2011), where mergers were similarly found to create jumps in the AVR, and are also consistent with idealized studies of the effect of mergers on disc kinematics (e.g. Quinn et al. 1993; Velazquez & White 1999; Kazantzidis et al. 2009).
However, contrary to quiescent galaxies, blurring the ages has a greater effect on the shape of the AVR for galaxies with mergers (black lines in Fig. 2). The jumps are considerably smoothened, so that they do not trace exactly the end of mergers anymore (especially for 30 per cent errors). Note also that for g22 the jumps actually become far less noticeable for age errors of 30 per cent, which suggests that errors on ages of 20 per cent or less are critical to detecting the traces of minor mergers in the observed AVR.
3.3 The velocity dispersion of young stars
For our seven galaxies, we find that the youngest stars (age < 500 Myr) have a vertical velocity dispersion that varies between 10 and 33 km s−1 at 2 Rd. In all galaxies, this value declines at increasing galactocentric distances, as shown in Paper I (fig. 6). For instance at 3 Rd, σz for young stars varies between 7 and 23 km s−1 depending on the galaxy.
It is likely that stars in our simulations are born too kinematically hot, because we do not fully resolve the structure of the interstellar medium. In the solar neighbourhood, Aumer & Binney (2009) find that the bluest (i.e. youngest) stars have a σz of only 6 ± 1 km s−1. This value is similar to the velocity dispersion of molecular clouds: Stark & Brand (1989) measured a cloud–cloud velocity dispersion of 7.8 ± 0.6 km s−1 for clouds within 3 kpc of the Sun. This dispersion depends on cloud mass, and the most massive GMCs have a cloud–cloud velocity dispersion of only 4 km s−1 (Stark & Lee 2005, 2006). If stars are born with such low velocity dispersions, they are then heated to σz ∼ 10–15 km s−1 within a few hundred million years (see also Fig. 1). This heating can be very fast: in an idealized simulation of a Milky Way type galaxy evolving in isolation, Renaud et al. (2013) find that stars initially born with a total velocity dispersion of 10 km s−1 are heated to 15 km s−1 in only 10 Myr. Assuming an isotropic velocity dispersion, this translates into a σz at birth of 5.8 km s−1. To obtain such a cold stellar component, they used a sub-parsec resolution (down to 0.05 pc) that we cannot achieve in cosmological simulations (our resolution is 150 pc).
It seems however that even if stars are not formed very cold, and even if we do not resolve that phase of intense early heating, the values of σz we obtain at 0.5–1 Gyr are reasonable compared to the Milky Way: in four of our simulations, young stars have a σz between 10 and 17 km s−1.
Our three other galaxies (g47, g48, and g22) have a higher σz for their young stars, up to 33 km s−1 for g22. These reflect a higher gas velocity dispersion for these galaxies. As shown in Fig. 4, what they have in common is both a higher stellar surface density, as well as higher star formation rate surface density, compared to the kinematically colder galaxies.
Relation between the vertical velocity dispersion σz for young stars at 2Rd and the stellar surface density at that same radius (top panel) as well as the local star formation rate surface density (bottom panel) for all seven galaxies. Both σz and the SFR are defined as corresponding to the last 500 Myr of evolution before z = 0.
A higher surface density has been shown to increase gas velocity dispersion through gravitational instabilities (Wada, Meurer & Norman 2002; Agertz et al. 2009; Bournaud et al. 2010). Another source of turbulence is the feedback from supernovae explosions, which increases with the star formation rate. Dib, Bell & Burkert (2006) showed that the atomic gas velocity dispersion does not depend on the SFR surface density up to a threshold corresponding to a starburst regime: for SFR densities greater than 0.5–1 × 10−2 M⊙ yr−1 kpc−2, the atomic gas velocity dispersion increases with the SFR density (see also Agertz et al. 2009; Tamburro et al. 2009). Interestingly, this threshold value is similar to the SFR density we find in our kinematically hottest galaxies, suggesting that feedback could be the source of the high velocity dispersion in the gas disc of these galaxies. Finally, it is also possible part of the heating is provided by satellite galaxies, either merging with the discs, or simply orbiting around the main galaxies.
It has also been observed that for molecular gas, the cloud–cloud velocity dispersion decreases with radius, and differs among local galaxies (Wilson et al. 2011); it is roughly twice smaller than the atomic gas velocity dispersion (Tamburro et al. 2009). Plausibly, different galaxies could have a different σz for their young stars.
In simulations, though, the gas velocity dispersion is also the result of numerical effects (e.g. Pilkington et al. 2011; Agertz et al. 2013). House et al. (2011) compared the vertical velocity dispersion in various simulations performed with different codes, resolutions and recipes. They find a ‘floor’ in σz of 15–20 km s−1 even in high-resolution simulations (50 pc), except when the star formation threshold is increased from 0.1 to 100 cm−3.
Our simulations are difficult to directly compare to the above results, because σz varies as a function of radius and from one galaxy to another. What we find is that an increase of the resolution by a factor of 2 does not change our results (see Section 6). This could mean that a much higher (sub-parsec) resolution would be necessary to successfully reproduce a very low σz for stars at birth (as in Renaud et al. 2013). We also find that increasing the star formation threshold from 1 to 5 H cm−3 decreases the overall normalization of the AVR (because we did not increase the star formation efficiency, so that the total mass is lower), but does not change its shape, in particular there is no fast heating for young stars.
4 THE ORIGIN OF THE AVR
To understand the origin of the continuous increase of σz with age in quiescent galaxies, we study how σz varies with time for each mono-age population. In the top panel of Fig. 5, we show σz at birth for a sample of mono-age populations for the three quiescent galaxies, as well as for g106, which shows similar characteristics to these galaxies (we here select stars being at 2Rd at z = 0, and trace them back to their time and place of birth). We find in all cases that 8–9 Gyr old stars form a different population, and are born significantly hotter than younger stars. They trace the end of a violent, high-redshift phase of star formation. For g47 and g106, this initial phase extends a bit longer, with stars still being born hot at ages of 7–8 Gyr. For instance, for g106, around 10 per cent of stars of that age and at that radius are accreted from satellites, and while the interactions contribute to disc heating, the total σz for these hot populations is dominated by the high σz at birth for stars born within the galaxy itself. After that phase of mergers, the younger stars are then generally born with a relatively constant σz, even though a slight decrease of the birth σz with time can usually be noticed (Fig. 5).
Investigating the origin of the slope in the AVR. The top panel shows σz at birth for a sample of mono-age populations in four of our simulated galaxies. The bottom panel shows the heating rate D for the same populations, where D is defined by |$\sigma _{z}^2(t)=\sigma _{z,0}^2+D (t - t_0)$|. While the oldest populations are born hot, for stars younger than ∼7 Gyr the σz at birth is relatively constant with time, and the slope of the AVR is due to subsequent heating of the stars.
We study how each mono-age population is heated by following the time evolution to z = 0 of the vertical velocity dispersion of its stars (again, for stars being at 2Rd at z = 0).
Disc heating can be modelled as a diffusion process in velocity space, in which case we can write |$\sigma _{z}^2(t)=\sigma _{z,0}^2+D (t - t_0)$| (with t0 and σz, 0 the values at birth). D is the diffusion coefficient, or the heating rate, and is plotted on the bottom panel of Fig. 5 for a sample of mono-age populations (we simply measure it from the difference between |$\sigma _{z}^2$| at birth and at z = 0).
These heating rates are generally below 200 km2 s−2 Gyr−1. They differ between galaxies: the lowest values are achieved for g92 and g106 (with D ∼ 40 km2 s−2 Gyr−1 for most mono-age populations), while g47 and g37 undergo more heating (D ∼ 100–150 km2 s−2 Gyr−1). While part of the heating could be due to spurious numerical effects (see discussion in Section 5.1), such numerical effects should be present in all simulations, still allowing comparisons between simulations. For instance, the much larger heating rate found in g37 compared to g92 suggests that this larger heating has a physical and not numerical origin.
For quiescent galaxies, the slope of the AVR in our simulations is thus governed by two separate effects: the oldest stars are born hot, while the rest of the disc (with ages under 8 Gyr) is mostly born with a similar σz and subsequently heated. These effects combine to produce a continuously rising σz with increasing age. This differs from the picture presented in Bird et al. (2013), where the vertical structure of the disc is mostly imprinted at birth (with also some subsequent heating). This could be the result of a different simulation technique (potentially for instance a stronger supernova feedback in their simulations), with young stars maybe born too hot so that there is not much room for subsequent heating.
5 INVESTIGATION OF HEATING MECHANISMS
5.1 Numerical effects
In the simulations presented in this paper, the relatively high σz for young stars, and the increase of σz with age for stars in galaxies with quiescent merger histories could both be numerical artefacts. In the Milky Way, the exact shape of the AVR is debated, but some studies find a low σz for young stars, and a rapid increase of σz in the first Gyr. If young stars are too hot in the simulations, we might just simply miss that first phase of intense heating and just have slower, more regular heating.
The velocity dispersion of stars at birth can be affected by the resolution of the simulation, by the recipe for star formation (and in particular the value of the gas density threshold above which stars are formed), by the recipe for supernova feedback (affecting the gas velocity dispersion), and also by our gas dynamics algorithm itself. In turn, subsequent heating is probably only affected by the resolution.
Low resolution can result in significant two-body heating, leading to equipartition of energy between rotational and random motion, thus thickening discs (e.g. Quinn et al. 1993). Massive dark matter particles are efficient at scattering stellar disc particles (Governato et al. 2004; Mayer 2004; Kaufmann et al. 2007), and also at increasing the gas thermal energy (Steinmetz & White 1997). The effects on gas should however be negligible for dark matter particle masses below 109 M⊙ (Steinmetz & White 1997), which is the case for all modern zoom-in cosmological simulations. The heating of the stellar disc should also be limited for haloes containing a few million dark matter particles. For instance, Widrow & Dubinski (2005) show discs undergoing only some minor heating in simulations using 2 million particles with each a mass of 2.9 × 105 M⊙ (similar to the dark matter particle mass in our simulations). Note also that Hänninen & Flynn (2002) show that halo black holes with a mass of 106 M⊙ are not efficient at heating discs. Compared to Hänninen & Flynn (2002), our use of a Particle-Mesh technique should reduce heating even further (e.g. Sellwood 2013).
For modern simulations, Sellwood (2013) argues that the most important source of numerical heating is actually the interactions between disc particles themselves, not between disc and halo. This effect is strongly dependent on the number of particles in the disc, and a few million disc particles are required to limit relaxation in simulated discs. While this result is not directly relevant for our simulations because we use a different simulation technique, the fact that all our simulated discs contain a few million stellar particles is encouraging.
Finally, in addition to relaxation effects, resolution also influences the growth of non-axisymmetric features in discs (Kaufmann et al. 2007; Khoperskov et al. 2007), and the properties of infalling satellites, both potential sources of heating. In the end, however, House et al. (2011) found that resolution did not have a strong impact on the AVR in their simulated galaxies, while increasing the star formation threshold from 0.1 to 100 H cm−3 strongly decreased the velocity dispersion of young stars.
We have tested the effect of resolution and star formation threshold on one of the simulated galaxies from the parent sample of Martig et al. (2012). Part of the results from these tests were already presented in appendix A of Martig et al. (2012), where it was found that the increased resolution did not modify the stellar mass and radial density profile of that simulated galaxy. We here extend these tests to the AVR.
The resolution test consists in a twice higher spatial resolution (corresponding to 75 pc), and a mass resolution increased by a factor of 6 (the mass of a gas particle is then 2500 M⊙). The test of the star formation threshold consists in increasing this threshold by a factor of 5 in the simulation at standard resolution (this moves the threshold from 1 to 5 H cm−3).
A twice higher resolution only has a very small impact on the AVR (Fig. 6). We find that the σz of young stars is similar to values at standard resolution, and the evolution of σz with age also has the same shape. It just seems that σz is consistently 1–2 km s−1 lower at high resolution, but this does not reconcile our simulations with observations of a fast rising AVR for small ages in the Milky Way. Note also that the constant offset of σz values at higher resolution indicates that whatever the heating mechanisms are, they are not strongly affected by resolution (at least, for the resolutions tested here).
Tests of the effect of an increased resolution (dashed blue line) or an increased star formation threshold (dotted red lines) on the shape of the AVR.
Increasing the threshold for star formation has a stronger effect on the simulated galaxy, partly because we did not increase the star formation efficiency so that the overall number of stars formed is lower, and the total stellar mass is only 80 per cent of the mass at standard resolution. With the higher threshold, the overall values of σz are thus shifted to lower values. In particular, σz for young stars becomes 14 km s−1 instead of 20 km s−1. However, the shape of the AVR is not affected, in particular there is no fast heating from 0 to 1 Gyr. If we arbitrarily re-scale all the values of σz to match the values of σz for young stars in the standard simulation, we then find a good match between both simulations.
5.2 Effect of radial migration
Recently, radial migration has been proposed as a source of disc heating (Schönrich & Binney 2009a,b; Loebman et al. 2011; Roškar et al. 2013), but we find no evidence for this in our simulations (as also discussed in Minchev et al. 2012; Minchev, Chiappini & Martig 2013, and is consistent with Bird et al. 2013). In Fig. 7, we show the fractional contribution of different types of stars (born in situ, formed in the inner disc, formed in the outer disc) to the total σz of each mono-age population. We find that for stars younger than 7–8 Gyr, the stars born in situ mostly determine the total σz (the total σz only differs from the σz of in situ stars by less than 10 per cent in all cases). Similarly to the findings by Minchev et al. (2012, their fig. 5), stars coming from the inner disc tend to have a slightly larger σz, and the opposite is true for stars from the outer disc. The net effect of migration is however to slightly cool the disc, because the heating from inner-disc stars is more than cancelled by the cooling from outer disc stars. Only in the very outer disc does the heating from outward migrators dominate, creating some flaring (Minchev et al. 2012). We emphasize here the importance of considering both inward and outward migrators when studying the effect of migration on disc heating or thickening.
Contribution of different populations of stars to the total vertical velocity dispersion as a function of age. For each mono-age population, we measure the fraction of its total σz (measured at 2Rd as previously) that is due to stars born in situ (i.e. in a 2 kpc ring centred on 2Rd), stars born in the inner disc, and stars born in the outer disc. For stars younger than 7–8 Gyr, the total velocity dispersion is very similar to the dispersion for stars born in situ, meaning that radial migration only has a minor effect on σz at that radius (if anything, it slightly cools the disc because of the combined effect of inward and outward migrating stars).
For the oldest stars, though, stars coming from the outer regions are the hottest. This is because the early mergers affect more significantly the outer disc, creating a flared morphology. The outward migration of cold but old stars from the inner disc is shown by Minchev et al. (2014) to be a key ingredient in reproducing some chemo-dynamical features observed in the RAdial Velocity Experiment (RAVE) and Sloan Extension for Galactic Understanding and Exploration (SEGUE) surveys (namely, the low vertical velocity dispersion for the most alpha-rich stars).
Ruling out migration as a source of heating at a radius of 2Rd still leaves a lot of possible sources of heating in the simulated discs. We investigate these sources of heating in the following sections.
5.3 Heating due to disc growth
To elucidate the mechanisms for heating, for each of the four galaxies of the previous section, we study the time evolution of one mono-age population, composed of the stars born at times between 8 and 8.5 Gyr, and that are found in a ring around 2Rd at t = 13.7 Gyr.
We show in Fig. 8 the time evolution of the baryonic surface density (i.e. including both gas and stars, which both contribute to the gravity in the disc plane) within an annulus at 2Rd (up to z of 3 kpc). We find that over the 5 Gyr that we consider, the surface density only increases by ∼10 per cent for g106 and g92, but it increases by 65–75 per cent for the other two galaxies. This could be an important source of heating. The exact contribution to the increase of σz actually depends on how the scaleheight of the population evolves with time (middle panel, Fig. 8). In all cases, z0 increases with time, with large differences between galaxies.
Left: time evolution of |$\sigma _z^2$| for one mono-age population in each of the four simulated galaxies (the slope of each line is the heating rate D as defined in Section 4). Middle: evolution of the scaleheight. Right: evolution of the baryonic surface density (gas+stars) at 2Rd. In thin dashed lines in the left-hand panel is shown the heating that can be attributed to the increase in surface density, computed using |$\sigma _z^2(t)=2 \pi G \Sigma (t) z_0$|(t = 8.5 Gyr). We find that disc growth only has a minor contribution to disc heating for g92 and g106, but is much more significant in g37 and g47.
This makes it difficult to estimate the contribution of gravity-driven heating to the increase of σz (and even more so because it is not the only heating mechanism present as explained in the next section). We provide a rough estimate of how much of the heating can be attributed to the increased surface density, by computing a gravity-driven σz, using the values of the surface density as a function of time, and assuming that the scaleheight remains fixed at its initial value. These are the dashed lines in Fig. 8, which show that the increase in surface density contributes to the heating rate by only 5 per cent for g92 and 14 per cent for g106 (their surface density evolves very little). On the other hand, it contributes by 37 and 57 per cent to the heating in g37 and g47, respectively.
Since these contributions are sub-dominant, in the next section, we investigate additional sources of heating.
5.4 Spiral arms, bar, and bending modes
A possible diagnostic of heating mechanisms is the ratio of vertical to radial velocity dispersions. In Fig. 9, we show the time evolution of σz, of σr and of σz/σr for the populations studied in the previous section (stars born between t = 8 and 8.5 Gyr). In all galaxies, we find a steady increase of σz with time, with the possible exception of g106 where σz seems to saturate to ∼22 km s−1. We similarly find that σr increases with time (with a saturation at ∼65 km s−1 for g37 over the last 3 Gyr of evolution). The resulting σz/σr ratio has a complex behaviour, suggesting different sources of heating from one galaxy to another (and confirming that most of the heating is not simply numerical).
Time evolution of the vertical velocity dispersion (σz, left-hand panel), radial velocity dispersion (σr, middle panel), and of the ratio between the two (σz/σr, right-hand panel), for the four galaxies with quiescent histories. We show for each galaxy the velocity dispersion for a population of stars born between 8 and 8.5 Gyr, and contained within an annulus at 2Rd at t = 13.7 Gyr (up to z = 3 kpc). The solid lines represent all stars, irrespective of their birth location, while the dashed line corresponds to stars born within the final 2Rd annulus (and thus corresponds to stars with minimal radial migration).
For g47, g92, and g106, σz/σr decreases over the 5 Gyr that we consider. For g37, σz/σr first decreases for 1.5 Gyr, and then returns to its initial value. These patterns are not consistent with numerical heating: Sellwood (2013) shows that σz/σr increases strongly for disc undergoing numerical heating (and in simulations with large number of particles, σz/σr evolves very little). A decrease of σz/σr is also inconsistent with the models of heating by GMCs and halo black holes by Hänninen & Flynn (2002), where σz/σr is either constant or slightly increasing with time.
Spiral arms and bars are more efficient at heating discs in the radial direction (Sellwood & Carlberg 1984; Carlberg & Sellwood 1985; Minchev & Quillen 2006). Vertical heating can then be obtained if the radial heating is redistributed in the vertical direction by GMCs or other massive objects (Carlberg 1987; Jenkins & Binney 1990; Jenkins 1992). Saha et al. (2010) also observe vertical heating linked with the presence of spiral arms, and speculate that even in the absence of GMCs the radial heating could be redistributed in the vertical direction by weak bending waves (as proposed by Masset & Tagger 1997).
A very likely source of heating in our simulations could then be the combination of spiral arms with either overdensities in the disc and/or weak bending waves. We show in Fig. 10 that such bending waves are indeed present in the simulations. We show the mean height above the mid-plane of stellar particles as a function of time for an annulus at 2Rd in each simulated disc, and we find low levels of vertical oscillations for all galaxies.
Time evolution on the average height above the plane for all stellar particles at 2Rvvd (up to z of 1 kpc). All galaxies show some oscillations (corresponding to bending modes), which are strongest for g37.
Interestingly, these bending waves are much stronger for g37, for which the evolution of σz/σr is also different. In this case, we suspect that the late increase of σz/σr could be due to the bending waves themselves (Khoperskov et al. 2010; Griv 2011). Khoperskov et al. (2010) also show that bending waves increase σz/σr mostly through an increase of σz, which is what we find for g37 where at late times σz increases while σr remains approximately constant (note that this was also observed by Saha et al. 2010). The bending waves themselves could be due to small satellites orbiting the simulated galaxies (Widrow et al. 2014), but could also naturally arise from internal instabilities (that are predicted to occur when σz/σr is below ∼ 0.4), or even be induced by the spiral arms (Faure, Siebert & Famaey 2014).
We conclude from the analysis in this section that disc heating in our simulations is due neither to numerical effects, nor to radial migration. The contributions to the heating vary from one galaxy to another, and include heating due to disc growth, and heating due to a combination of spiral arms and bars coupled with overdensities in the disc and vertical bending waves.
6 SUMMARY
We study the AVR in a sample of seven simulated disc galaxies. All galaxies first undergo an active phase of mergers at high redshift, creating in all cases a thick stellar component with a high velocity dispersion. Amongst the old stars, we find that most of them are actually born kinematically hot, with an additional contribution of accreted stars. At lower redshift (corresponding to stellar ages smaller than 8–9 Gyr), three of the simulated galaxies have a quiescent history while the others have various types of mergers. We find that in the quiescent galaxies σz increases smoothly with age. When fitted with a power law, we find heating indices close to 0.5, making our simulated galaxies consistent with results from the GCS for the solar neighbourhood (Nordström et al. 2004; Holmberg et al. 2007), rather than with observations suggesting that disc heating saturates after a few Gyr (e.g. Soubiran et al. 2008). Note that such a gradual heating was also found in other zoom cosmological simulations, as for instance shown by Brook et al. (2012) or Bird et al. (2013), even if they do not measure the heating index in their simulations.
By contrast with Bird et al. (2013), though, we mostly find that the slope of the AVR is not imprinted at birth, but is due to heating. While we cannot pin down the exact mechanisms responsible for that heating, we show that radial migration is not one of them. Because of the combined contribution of hotter stars coming from the inner disc and cooler stars coming from the outer disc, we even find that the total effect of migration is to very slightly cool the disc, at least for the radius we examined, and for stars younger than ∼8 Gyr. For older stars, the effect is reversed, because stars coming from the inner disc are cooler: this is due to early mergers creating a flared structure for the old stars.
Finally, we find that σz is very sensitive to mergers (see also e.g. Quinn et al. 1993; Velazquez & White 1999; Kazantzidis et al. 2009; House et al. 2011): even 1–10 to 1–15 mergers create jumps in the AVR. This can in theory be used to probe the merger history of the Milky Way. However, we find that errors of 30 per cent on stellar ages, even if they do not affect the shape of the AVR for quiescent galaxies (except for the oldest stars, for which age errors significantly decrease the measured σz), can significantly blur the signatures of mergers (errors of 20 per cent or less are necessary to properly detect these signatures). Stronger constraints on the history of the Milky Way can be obtained by combining the AVR with the structure of mono-abundance populations (Bovy et al. 2012b,c). As discussed in Paper I, our simulations, combined with current data, favour a very quiescent history for the Milky Way in the last 9 Gyr, with a potential contribution of mergers to the build-up of the thickest and oldest disc components.
We thank the referee for a constructive report. We also thank Hans-Walter Rix and Jo Bovy for useful comments. MM is supported by a Humboldt Research Fellowship. CF acknowledges financial support by the Beckwith Trust.
