The innate origin of radial and vertical gradients in a simulated galaxy disc
We examine the origin of radial and vertical gradients in the age/metallicity of the stellar component of a galaxy disc formed in the APOSTLE cosmological hydrodynamical simulations. Some of these gradients resemble those in the Milky Way, where they have sometimes been interpreted as due to internal evolution, such as scattering off giant molecular clouds, radial migration driven by spiral patterns, or orbital resonances with a bar. Secular processes play a minor role in the simulated galaxy, which lacks strong spiral or bar patterns, and where such gradients arise as a result of the gradual enrichment of a gaseous disc that is born thick but thins as it turns into stars and settles into centrifugal equilibrium. The settling is controlled by the feedback of young stars; which links the star formation, enrichment, and equilibration timescales, inducing radial and vertical gradients in the gaseous disc and its descendent stars. The kinematics of coeval stars evolve little after birth and provide a faithful snapshot of the gaseous disc structure at the time of their formation. In this interpretation, the age-velocity dispersion relation would reflect the gradual thinning of the disc rather than the importance of secular orbit scattering; the outward flaring of stars would result from the gas disc flare rather than from radial migration; and vertical gradients would arise because the gas disc gradually thinned as it enriched. Such radial and vertical trends might just reflect the evolving properties of the parent gaseous disc, and are not necessarily the result of secular evolutionary processes.
keywords:galaxies: kinematics and dynamics – galaxies: structure – galaxies: formation – galaxies: evolution
The ages and metallicities of disc stars in the Milky Way (MW) show well-established radial, vertical, and kinematic gradients. In the solar neighbourhood, for example, the velocity dispersion of stars increases with age (e.g., Wielen, 1977; Quillen & Garnett, 2001; Holmberg et al., 2009), and age correlates in turn with total metallicity and, in particular, with the abundance of elements at given [Fe/H] (e.g., Freeman & Bland-Hawthorn, 2002; Haywood et al., 2013). In the disc midplane, the average age and metallicity of stars decline with increasing cylindrical radius, (e.g., Boeche et al., 2013; Hayden et al., 2014). The radial age trend is preserved but the metallicity gradient appears to reverse when considering stars at fixed vertical distances, , away from the disc symmetry plane (e.g., Anders et al., 2014; Bergemann et al., 2014; Martig et al., 2016; Ness et al., 2016).
At fixed , the vertical scaleheight of stars increases with decreasing metallicity or increasing age (e.g., Mikolaitis et al., 2014; Casagrande et al., 2016); these gradients depend on and become steeper at smaller radii (Hayden et al., 2014). Finally, at fixed metallicity, the vertical scaleheight of stars increases (‘flares’) from the centre outwards, especially for stars with low Fe ratios; the flare depends on metallicity, becoming more pronounced at lower [Fe/H] (Bovy et al., 2016).
Although the detailed form of each of these correlations vary somewhat from study to study, their qualitative nature seems robust. Their physical origin, however, is still a matter of debate. One possibility is that, to first order, the differences between old (metal-poor) stars and young (metal-rich) stars are driven by the integrated effect of secular evolutionary processes, such as orbital scattering off molecular clouds, radial migration induced by non-axisymmetric disc features, or dynamical heating by accretion events, external perturbations, and disc instabilities. These processes naturally affect older stars more, inducing vertical gradients and systematic trends between age/metallicity and kinematics that resemble the observed ones (for a review, see Rix & Bovy, 2013, and references therein).
A contrasting view is that stellar gradients are largely imprinted at birth, and thus reflect the properties of the parent (gaseous) disc at the time of formation of each coeval stellar population. In this scenario, the star formation and enrichment history of the disc are linked to the evolution of its vertical structure and to the timescale of its dynamical equilibration (Brook et al., 2012; Bird et al., 2013; Stinson et al., 2013; Grand et al., 2015; Miranda et al., 2016; Ma et al., 2017; Grand et al., 2017). At early times the gaseous disc was thick and metal poor, but it enriched as it thinned down, leaving behind a strong vertical gradient in the age and metallicity of the stars at each fixed radius, in a process reminiscent of the classic Eggen et al. (1962) dissipative collapse picture. Radial gradients arise in this model simply because the gaseous disc assembles inside out and is denser near the centre, leading to faster transformation of gas into stars and more rapid enrichment than in the outskirts.
Each of these views has its pros and cons. For example, secular evolution models have difficulty accounting for the radially increasing scaleheight (‘flare’) of stars, since kinematic perturbations are expected to be weaker at larger radii on account of lower densities, and to proceed more slowly there because of longer orbital times. Birth models, on the other hand, have difficulty explaining why the vertical settling and enrichment timescales needed to explain the data are so much longer than the local dynamical timescale, which would naively be expected to set the collapse and equilibration timescale of the disc.
In the Milky Way, a further complication is the presence of an apparent ‘gap’ in the vs Fe abundances of stars at fixed , which has often been interpreted as signalling the presence of two disc populations of possibly distinct origin: a relatively old, -rich ‘thick’ disc mixed with a younger, -poor ‘thin’ population of overlapping [Fe/H] (see; e.g., Venn et al., 2004; Bensby et al., 2005; Navarro et al., 2011; Recio-Blanco et al., 2014, and references therein).
These complexities explain why ‘chemodynamical’ models resort to features from both formation scenarios in order to reproduce the structure of the Milky Way disc(s). Radial migration; external accretion; self-enrichment; gas outflows and inflows; these are some of the mechanisms invoked to reproduce the peculiar kinematic and chemical structure of the Galactic thin and thick discs (see; e.g., Matteucci & Francois, 1989; Chiappini et al., 2001; Minchev et al., 2013, 2014; Schönrich & McMillan, 2017, and references therein). The specificity of such models, however, hinders their general applicability to understanding the relative importance of various physical mechanisms in a typical disc galaxy.
Direct cosmological hydrodynamical simulations of disc galaxy formation offer a complementary approach, where, modulo the algorithmic choices for simulating star formation and feedback, the roles of evolutionary effects and of conditions ‘at birth’ may be contrasted and assessed in statistically significant samples. The price to pay is that very few, if any, of the models will reproduce the detailed structure of the Milky Way, implying that the physical origin of observed Galactic trends must be inferred from judicious application of lessons learned from simulated galaxies that might not necessarily resemble our Galaxy in detail.
We adopt this approach here, where we analyze the formation of a disc galaxy in the APOSTLE111A Project Of Simulating The Local Environment suite of cosmological hydrodynamical simulations of the Local Group (Sawala et al., 2016; Fattahi et al., 2016). Our analysis focusses on the physical origin of vertical and radial gradients of disc stars, on how they reflect the properties of the gaseous disc from which they descend, and on the effect of secular evolutionary processes that operate after their formation. Our paper is organized as follows. We describe briefly the APOSTLE numerical simulations in Sec. 2. Our main results are presented in Sec. 3. A toy model meant to illustrate how these results may be used to interpret the physical origin of vertical gradients is presented in Sec. 4. We end with a brief summary of our main conclusions in Sec. 5.
2 Numerical Simulations
2.1 The Apostle project
The APOSTLE suite of zoomed-in CDM cosmological hydrodynamical simulations evolves 12 volumes from a large cosmological box selected to resemble the main characteristics of the Local Group (LG). Each volume is centred on a pair of dark matter halos of combined virial222The virial properties of a dark matter halo are computed within spheres of mean density the critical density for closure, . Virial quantities will be indicated with a ‘200’ subscript. mass of order . The pairs at lookback time 333We shall hereafter use lookback time, , to quote simulation results, in order to prevent confusion between ‘redshift’ and the vertical coordinate, . are at a relative distance of kpc, approach with radial velocity of up to km/s, have relative tangential velocities km/s, and have no other halos of comparable mass within Mpc. Details of the selection procedure are presented in Sawala et al. (2016); Fattahi et al. (2016).
All APOSTLE volumes were evolved with the same code developed for the EAGLE project (Schaye et al., 2015; Crain et al., 2015), which has been calibrated to reproduce some general properties of the galaxy population, such as the galaxy stellar mass function and the average size of galaxies as a function of mass. APOSTLE extends the EAGLE project to LG-like volumes and higher resolutions, using the fiducial choice of parameters (i.e., the ‘Ref’ model in the nomenclature of Schaye et al., 2015).
Each APOSTLE volume is resimulated at three different numerical resolutions, differing by successive factors of in particle mass (labelled L1 to L3 in order of decreasing particle count). Our analysis will focus on the most massive of the two main galaxies in volume AP-4-L2, using the nomenclature introduced by Fattahi et al. (2016). For that run, the high-resolution particle masses are (dark matter) and (gas). These particles interact with a gravitational potential softened with a Plummer-equivalent spline kernel length fixed at pc for and constant in comoving units at higher redshifts.
APOSTLE runs assume a WMAP-7 cosmology throughout, with parameters , , , and (Komatsu et al., 2011).
2.2 The simulated galaxy
Dark matter halos are identified using a friends-of-friends (FoF; Davis et al., 1985) algorithm run with linking length 0.2 times the mean interparticle separation. Individual self-bound structures within these halos are then identified recursively by the groupfinding algorithm SUBFIND (Springel et al., 2001; Dolag et al., 2009). The galaxies inhabiting the main subgroup of each FoF halo are referred to as ‘centrals’ and the rest as its ‘satellites’, if found within the virial radius of the central halo.
As stated above, we focus our analysis on a single galaxy, the more massive of the pair of central galaxies in APOSTLE volume AP-4-L2. This galaxy was selected because, at , it (i) has a prominent stellar disc; (ii) shows no obvious morphological peculiarities; (iii) has formed stars steadily throughout its history; (iv) has had a relatively quiet recent merging history; and (v) has a gas/stellar mass ratio not unlike that of the Milky Way. We emphasize that this galaxy is not a model of the Milky Way. The properties of this galaxy are representative of other disc-dominated galaxies in APOSTLE selected using similar criteria; however, since our analysis involves the detailed tracking of individual star particles it is easier to illustrate our main results by focussing on a single galaxy. Fig. 1 shows face-on and edge-on views of the galaxy, in a box kpc on a side.
3.1 The assembly of the galaxy
The top panel of Fig. 2 shows the mass assembly history of the galaxy, as measured by the total mass of each component within a radius of kpc (physical) from the centre of the progenitor halo. Most of the galaxy mass within that radius is assembled quickly: of the dark matter mass was already in place Gyr ago; the corresponding lookback time for the baryonic component is Gyr.
Most baryons are accreted as gas and then turn gradually into stars, as shown in the middle panel of Fig. 2. Gas is delivered to the centre of the galaxy during an early period of rapid merging that is essentially over Gyr ago. The last episode of early accretion is depicted in Fig. 3, where we see that the accretion of a substantial amount of gas along a filamentary structure accompanies a minor merger that happened roughly Gyr ago.
After that accretion event the evolution of the galaxy is largely quiescent, except for an episode of gas accretion that started roughly Gyr ago. This late accretion adds of gas, but not many stars. This addition increases the final baryonic mass of the galaxy by , but nearly doubles its gas mass at the time of accretion. The extra gas (which is of relatively low metallicity) serves to ramp up the star formation rate of the galaxy and to decrease its average metallicity at late times.
The middle panel of Fig. 2 also makes clear that the large majority of stars in this galaxy were born in situ444At very early times, when only a small fraction of the final galaxy has been assembled, the main progenitor is not well defined. This leads to temporary uncertainties in the assignment of in-situ vs accreted stars. These uncertainties have no discernible consequences on our main conclusions.: fewer than of all star particles at formed in other systems, and roughly of those were accreted more than Gyr ago.
The quiet late merging history of the galaxy, together with the fact that most baryons are accreted in gaseous form, favour the formation of a centrifugally supported disc, as may be seen in the bottom panel of Fig. 2. This panel shows the circularity parameter of star particles, , measured at , as a function of their formation time. The parameter is defined as the -component of the specific angular momentum of a star in units of that of a circular orbit of the same binding energy, and it therefore varies between for stars in the midplane. Some stars may have slightly exceeding unity if they are outside the disc midplane. The axis of the disc is defined at all times as the direction of the angular momentum vector of star particles younger than Gyr old.
Star particles older than Gyr formed before the early merging period of the galaxy was over: indeed, of them are accreted and all together they have no net sense of rotation. Over the next Gyr, however, a disc gradually forms and the average circularity climbs to unity. Stars formed during that epoch (i.e., ages between and Gyr) have a broad distribution of circularities, and the great majority (more than ) formed in-situ. Stars younger than Gyr, on the other hand, nearly all formed in-situ and are found in approximately circular orbits in a thin coplanar disc (i.e., ).
The various panels of Fig. 4 show a few snapshots555This figure was inspired by Fig. 1 of Ma et al. (2017). of the evolution of the galaxy. Here, from right to left, each column shows different orthogonal projections of star particles formed, within a narrow interval of time, , , , and Gyr ago, respectively. The first two rows show each coeval population of stars at . The third and fourth rows, on the other hand, show them at the time of their formation. The gas component of the galaxy is shown in the bottom two rows, for each time.
Fig. 4 shows that, as expected, star particles trace, at the time of formation, the densest regions of the gaseous disc. Interestingly, however, the spatial configuration of stars evolves little after formation, and remains similar at . The stars of the simulated disc, therefore, seem to provide today a relatively faithful tracer of the properties of the gaseous disc at the time of their formation. This is the central result of our study, and will be a recurrent theme in the analysis that follows.
3.2 Age and metallicity gradients
The evolution of the gaseous disc discussed in the previous subsection leaves discernible gradients in the age and metallicity of its descendent stars at . We show this in Fig. 5, where the heat maps show, as a function of the cylindrical coordinates, and , the average age and metallicity of star particles. The gradients shown are characteristic of the ‘inside-out/upside-down’ disc formation scenario described in earlier work (see, e.g., Bournaud et al., 2009; Brook et al., 2012; Bird et al., 2013; Stinson et al., 2013; Miranda et al., 2016; Ma et al., 2017).
In the disc midplane the average age and metallicity decrease monotonically with increasing . At fixed , the ages of star particles increase and their metallicities decrease with increasing . These gradients are steeper near the centre of the galaxy than in its outskirts, implying that the characteristic scaleheight of stars of fixed metallicity (or age) ‘flares’ outwards. This flare is the reason why the radial metallicity gradient at fixed away from the plane ‘inverts’: for example, at fixed kpc, the metallicity increases with , a trend opposite to that seen in the midplane.
3.3 The origin of radial and vertical gradients
The origin of the gradients highlighted in the previous subsection may be elucidated by contrasting the properties of stars at with those of their parent gaseous disc at the time of their formation. We do this in Fig. 6, selecting for analysis star particles formed between and Gyr ago. This period is simpler to analyze because it is characterized by a single major accretion episode, which, as discussed in Sec. 3.1, delivered a relatively large amount of gas to the forming galaxy. The accreted gas quickly assembles into a thick disc structure that slowly thins down as it forms stars (see; e.g., Fig. 4).
This evolution is depicted in the left column of Fig. 6, where, from top to bottom, we show the gas surface density profile, , its half-mass scaleheight profile, , and its average metallicity profile, [Fe/H], coloured as a function of time. These three panels show the evolution of a self-enriching, slowly-thinning, flared gas disc that gradually transforms most of its gas into stars.
Interestingly, as shown in the right column of Fig. 6, stars show, at , a very similar structure to the gas in the radial and vertical directions, when selecting their ages to correspond to the same times chosen for the gas in the left panels. The clear correspondence between gas and stars as a function of time/age indicates that the origin of the stellar gradients lies in the evolving structure of the gaseous disc. Stars inherit the gas properties at birth and, to first order, preserve them until .
Note that because the gaseous disc is flared at all times, and thins down as it enriches, stars of fixed metallicity form at different times at different radii and with different scaleheights. The flare of stars of fixed metallicity thus reflects the flaring of the gas disc, modulated by the gradual thinning that occurs at each radius.
The open circles and squares in Fig. 6 illustrate this. These symbols track two different metallicities, [Fe/H] and , respectively. In the bottom-left panel, the colour underneath each symbol indicates at what lookback time, as a function of radius, the average gas metallicity reached each of those values. In the bottom-right panel, on the other hand, the colours indicate, as a function of radius, the average age of star particles with each of those metallicities, at . The good agreement between gas and stars demonstrates that most stars have not migrated far from the radii where they were formed.
The lookback times (for the gas) and ages (for the stars) mentioned in the preceding paragraph are then used in the middle panels of Fig. 6 to indicate the half-mass scaleheight of the gas at each of those times (left panel) or the half-mass scaleheight of the stars of each of those ages at . Again, the close resemblance between the flare of gas and stars selected in this way indicates that the stellar gradients originate in the properties of the gaseous disc at the time of the formation of each coeval stellar population.
One consequence of this resemblance is that, at , star particles of fixed metallicity are flared (see Fig. 7). The steepness of the flare depends on a combination of several factors: how strong the gas disc flare is, how quickly the gas self-enriches, and how fast the disc thins down at various radii.
Had the disc thinned down much more rapidly, for example, so that at late times the scaleheight in the outskirts was similar to that of the inner regions at early times, the flare in stars of fixed metallicity would be much weaker, or might even disappear altogether. This interplay between thinning, enrichment, and the flare of the original gaseous disc explains also why the flaring depends on metallicity. Taking stars of [Fe/H] and repeating the exercise, we find that the flare is less pronounced than that of [Fe/H] at (see the right panel of Fig. 7).
The same exercise also explains the origin of vertical metallicity gradients at fixed . This is shown in Fig. 8, where the red curves indicate the half-mass scaleheight of stars as a function of metallicity, at the present time, and for two different radii, kpc (left) and kpc (right). The negative gradients (scaleheights decrease with increasing [Fe/H]) result from the fact that the parent gaseous disc was gradually thinning down as it enriched.
This is confirmed by the blue curves in Fig. 8, which indicate the half-mass scaleheight of the gas disc at the time when its average metallicity at that radius reached the value listed in the abscissa of Fig. 8. In other words, [Fe/H] is a proxy of time for the gas, and, in the absence of substantial accretion, tracks the enrichment process of the gaseous disc. The excellent agreement between the vertical gas evolution and the stellar gradients at indicate, again, that stars provide, to first order, a snapshot of the properties of the gas disc at the time of their formation.
3.4 Stellar migration after formation
The previous discussion suggests that the kinematics of star particles evolve little after their formation. We show this quantitatively in Fig. 9, where the top panel shows, as a function of initial cylindrical birth radius, , and for different stellar age bins, the average net change in distance from the galaxy rotation axis. The net changes are rather small ( over the whole radial extent of the disc), indicating that little net migration has happened since stars were born. More importantly, the rms radial change is also rather small (as seen in the bottom panel of Fig. 9), dropping from near the centre to in the outskirts of the disc, even for stars as old as Gyr.
The changes induced on radial gradients by migration are consequently rather small. We demonstrate this explicitly in Fig. 10, where we show the radial profiles of the vertical scaleheight (top) and of the metallicity (bottom) of stars formed at two different lookback times. Dashed lines indicate the properties of such stars at and solid lines the same at the time of their formation. This figure shows explicitly that there is little radial mixing of stars after their formation: the radial metallicity gradients of stars of given age are essentially the same at as at the time of formation (bottom panel). There is also little ‘heating’ in the vertical kinematics of stars after formation: the vertical gradients are again very similar at the time of formation as at (top panel).
Since the simulated disc lacks a prominent bar, massive spiral arms, or dynamical analogues of giant molecular clouds, it is perhaps not unexpected that secular evolutionary processes such as radial migration play a minor role in the disc’s dynamical evolution. Migration in galaxies with stronger spiral patterns (and, in particular those where the pattern recurs periodically; see, e.g., Vera-Ciro et al., 2014) is expected to be more important, and our results thus do not exclude that migration may have played an important role in the Milky Way. This, however, does not affect our main result: the radial and vertical trends in the simulated disc reflect the conditions of the gaseous disc at birth, with secular evolutionary processes playing a minor role.
3.5 Vertical thinning of the gas disc
The results presented above indicate that star particles of fixed age faithfully trace the properties of the gaseous disc at the time of their formation: their vertical and radial gradients are set at birth, and evolve little thereafter. This elicits two important questions: (i) what determines the vertical scaleheight of the gas, its outward flare, and the timescale of its thinning?; and (ii) since, unlike the gas, stars are not subject to hydrodynamical forces, why do stars trace the properties of the gas so closely?
The scaleheight of an equilibrium gaseous disc is given by the balance between pressure forces and the combined vertical compressive forces of the dark matter halo and the disc (see Benitez-Llambay et al., 2017, for a recent discussion). For an isothermal gas disc with sound speed , embedded in a dark matter halo with circular velocity profile , and surface density profile , it is straightforward to show that, in the thin disc approximation, the vertical compressive force is given by
where the first term of the right side indicates the contribution of the dark matter halo and the second that of the disc, which may also include a stellar component. ( is the surface density enclosed between .)
When the halo term dominates, the disc is usually termed “non-self-gravitating” (NSG) and its characteristic scaleheight is given by
On the other hand, when the disc term dominates the vertical force the disc is “self-gravitating” (SG) and its characteristic scaleheight is given by
where is total surface density integrated over all . Because the vertical density law differs when the disc is SG or NSG the characteristic scaleheight values given by Eqs. 2 and 3 are not direct measures of the half-mass scaleheights, , which is what we actually measure in the simulations. The proportionality factors, however, may be easily computed in each case and may be consulted in Benitez-Llambay et al. (2017).
Note that, in general, we expect isothermal exponential gas discs with flat circular velocity curves to ‘flare’ outwards, since, in that case, the thickness would be if non-self-gravitating, and if self-gravitating.
We compare in Fig. 11 the above expectation with the gas disc half-mass scaleheights measured at two different lookback times, and Gyr ago. The thick grey curve shows the result of the simulation, whereas the SG and NSG scaleheights are shown with dotted and dashed red curves, respectively. We use the midplane sound speed to account for the density-dependent EAGLE equation of state; as a result, the effective sound speed in Eq. 2 depends (weakly) on . The solid red line is the actual expected value when considering the disc and halo combined vertical forces. Clearly the simulated disc is much thicker than expected given these considerations, at essentially all radii.
Indeed, what sets the thickness of the disc is actually the balance between the vertical forces and the ‘pressure’ provided by random bulk motions in the gas. The blue dashed curves in Fig. 11 show again the NSG solution, but replacing the sound speed, , with the vertical velocity dispersion of the gas particles, , in Eq. 2. The agreement between this curve and the simulation results indicates that the gas disc is kept thick by the vertical random motions of the gas, and, in particular, of the star forming gas clouds that may be clearly seen in the bottom panels of Fig. 4, especially at early times.
These random motions should be quelled in a couple of dynamical timescales, but, as seen in the right panel of Fig. 11, they still dominate in the outskirts of the disc after Gyr of evolution, or more than circular orbit periods at kpc. (Orbital times as a function of radius, as well as circular velocity profiles of the simulated galaxy are provided in the Appendix.) The gas disc is actually kept thick by the feedback energy provided by evolving stars—this is deposited as thermal energy in the disc, where it is able to blow bubbles and to push gas outside the disc quite efficiently, especially when star formation rates are high.
We show this in Fig. 12, where we plot the half-mass scaleheight evolution of a gas disc formed in an idealized simulation which evolved a system tailored to match approximately the mass, size, and angular momentum of the gaseous disc at Gyr (see Fig. 6). In this simulation, gas is allowed to cool and collapse from a kpc-thick rotating ‘slab’ embedded in a spherical Navarro-Frenk-White (NFW, Navarro et al., 1996, 1997) dark matter halo potential with parameters chosen to match the circular velocity profile of the galaxy at that time (Fig. 16).
The solid lines in Fig. 12 indicate the evolution of the half-mass scaleheight at two different radii ( and kpc), when the system is evolved in the absence of star formation or feedback. As expected, the disc quickly collapses vertically to a pressure-supported disc roughly pc thick. The simulation is ended after just Myr because the disc breaks into massive self-bound clumps, making estimates of its thickness unreliable. When star formation is included, however, the gaseous disc settles into a much thicker structure which is supported by the gas bulk motions induced by the feedback energy released by evolving stars. (See dashed lines in Fig. 12.)
Interestingly, the idealized disc with feedback is still thinner than the APOSTLE galaxy (compare with the left panel of Fig. 11), suggesting that there may be an additional source of vertical support in the APOSTLE simulation, possibly related to bulk motions induced by the continuous gas accretion, which is not included in the idealized runs. In any case, Fig. 12 shows clearly that the thickness of our simulated disc is largely set and controlled by feedback-induced vertical motions in the gas, whose effective pressure far exceeds the thermal support.
Feedback is thus a critical ingredient for understanding the origin of vertical gradients in our simulations, connecting the star formation history at each radius and its consequent enrichment with the evolving thickness of the star forming disc. In this scenario, the disc thins down as the gaseous disc is depleted, star formation abates, and feedback (and, possibly, accretion) heating becomes less effective. A strong prediction is then that the star formation, enrichment, and thinning timescales should all be linked at all radii, as it may have already been observed in the inner Galaxy (see; e.g., Freudenburg et al., 2016).
The same scenario explains why the vertical structure of stars tracks so closely that of the gas at the time of their formation. This is because the star forming gas is best described as a collection of dense clouds with appreciable bulk velocities, rather than as a fluid in vertical hydrostatic equilibrium supported by thermal pressure. Indeed, had star particles been born of gas in thermal hydrostatic equilibrium the stellar disc would be much thinner than its parent gas disc due to the loss of vertical pressure (see Benitez-Llambay et al., 2017). In our simulations, however, star particles inherit the bulk motions of the gas clouds from which they are born. Since these motions are dominant over thermal pressure forces there is little difference between gas clouds and newly formed stars, explaining why stars trace faithfully the properties of the gas at the time of their formation.
We end this discussion with a caveat by noting that, when feedback is included, the gaseous disc thickness exceeds that of the thin disc of the Milky Way (see Fig. 12). This elicits questions as to whether a more realistic star formation/feedback model would yield radial and, in particular, vertical trends similar to those we report here. A definitive answer to this question will need to wait until our modeling improves to match the observed disc thickness of Miky Way-like spirals. However, supporting evidence for our conclusions may be gleaned from simulations that yield much thinner discs than ours but, qualitatively, the same trends we report here (see, e.g., Ma et al., 2017, and references therein). We are planning to revisit this issue in future work.
4 A model for the origin of metallicity gradients
The scenario discussed above, where gradients result from the intertwined evolution of the gaseous disc’s star formation, enrichment, and vertical structure, may be encapsulated in a simple model that can be used to interpret the origin of various chemical and kinematical trends in the Galaxy. The model assumes that stars inherit the properties of the gas at the time of formation and preserve them to the present. It also requires, at a minimum, a characterization of the disc surface density profile, , its thickness profile, , and its average metallicity profile, [Fe/H], as a function of time.
4.1 Gas disc evolution parametrisation
The evolution of the vertical structure of the APOSTLE gaseous disc may be approximated simply by
where is lookback time, is the central half-mass disc scaleheight at the present time, is a characteristic flaring radius, and we have assumed that the thinning timescale of the disc, , is independent of .
A similar description may be adopted for the evolution of the metallicity gradient. We assume that the average metallicity profile is given by [Fe/H], and that it may be approximated by
where is the central average gas metallicity at the present time, , and we have assumed that the enrichment timescale, , is independent of .
We compare these parametrisations with the evolution of the simulated disc in Fig. 13, for three different radii. The thinning and enrichment timescales that best approximate the simulated disc are Gyr and Gyr, respectively. This figure shows that the functional forms adopted above are adequate, at least for Gyr. At more recent lookback times the disc accretes a substantial amount of metal-poor gas (see middle panel of Fig. 2), lowering the average metallicities at all radii. This non-monotonic behaviour cannot be reproduced by our simple formulation.
4.2 Application to the Milky Way
Given the assumptions of the parametrisation adopted above, the radial and vertical gradients that result will be largely set by the flaring profile and metallicity gradient adopted for the gaseous disc at , as well as by the ratio between the thinning timescale, , and the enrichment timescale, . Some intuition may be gained from considering hypothetical cases when one timescale is much longer than the other.
For example, if enrichment proceeded much faster than thinning, then all populations, regardless of metallicity, would have the same flaring profile (that of the gaseous disc), and there would be no vertical metallicity gradient at fixed . If, on the other hand, thinning proceeded much faster than enrichment, then we would expect strong vertical gradients to develop at fixed . In this case, the flaring profile of stars of fixed metallicity would depend in detail on the gas disc flare, as well as on the ratio of thinning and enrichment timescales. Depending on these parameters, stars of given might show no flare even though the gaseous disc is always flared.
We illustrate this next with a simple application to the Milky Way. We emphasize that this is not meant to be quantitative model of the chemical evolution of the Milky Way disc(s) but rather an illustration of how the results discussed in the previous subsections may be used to interpret the data.
The model requires the present-day gas MW density/metallicity profiles. We show these in the left panels of Fig. 14: the top two show the vertical gas scaleheight profile taken from the fits to 21-cm data presented by Kalberla & Dedes (2008). The gas disc is clearly flared, and can be well approximated, at , by Eq. 4 ( kpc and kpc), as shown by the bue lines. The bottom-left panel shows the Cepheid metallicity profile from Genovali et al. (2014), which we shall take to represent the present-day radial metallicity profile of the Galaxy midplane, and which we approximate with and dex/kpc. This gradient can also be approximated by the parametrisation adopted in Eq. 5, as shown by the blue lines in that panel.
The only remaining parameters are and . We choose two cases, again mainly for illustration. Model 1 (M1) assumes ; i.e., Gyr, and Gyr. Model 2 (M2), on the other hand, assumes that thinning proceeds on a much longer timescale than enrichment; i.e., Gyr, and Gyr. Note that the two models differ only in the thinning timescale.
The evolution of both thickness and average metallicity is shown for three different radii; , , and kpc, in the right panels of Fig. 14. In M1 the disc thins down quickly from a much thicker earlier configuration to its final state, whereas it evolves much more gradually in M2. The metallicity evolution is identical in both models, which assume that the disc has enriched at all radii by more than one dex666We note that the latter result differs from the solar neighbourhood, where the age-metallicity relation has large scatter and is relatively flat in the age range - Gyr (Nordström et al., 2004; Haywood et al., 2013). We emphasize again that our simulated galaxy is not meant to be a model for the Milky Way, but rather a basic illustration of how our results may be applied to gain insight into the Milky Way evolution. in the past Gyr.
The effect of these two different evolutionary patterns on the flaring profile of stars of fixed metallicity is shown in the left panels of Fig. 15. In the case of Model 1 the enrichment and thinning timescales are comparable and, for the flaring profile of Kalberla & Dedes (2008), leads to populations of stars that show no flare at fixed metallicity.
Model 2, on the other hand, leads to well-defined flares in populations of fixed metallicity, as shown in the bottom-left panel of Fig. 15. Interestingly, both trends are seen in the Milky Way disc, where there is compelling evidence for the presence of two chemically-distinct populations (Bovy et al., 2016): an -rich disc (the traditional ‘thick disc’) with no flares (like M1), and an -poor (‘thin’) disc whose subpopulations clearly flare outwards (like M2). (See right panels of Fig. 15.)
Although the presence of these flares has been taken as evidence for the effects of radial migration, our results argue that such flares might also result from the gradual thinning of the gaseous disc. This should not be viewed as implying that radial migration does not occur in the Milky Way; only that the mere presence of a flare does not guarantee that it has been caused by radial migration, and that other explanations should also be considered (see; e.g., Schönrich & Binney, 2009; Loebman et al., 2011; Minchev et al., 2012; Kubryk et al., 2013; Roškar et al., 2013; Vera-Ciro et al., 2014; Grand et al., 2015; Vera-Ciro et al., 2016; Kawata et al., 2017, for a more comprehensive discussion of the role of radial migration on disc flaring).
5 Summary and Discussion
We use a cosmological simulation from the APOSTLE project to study the origin of radial and vertical gradients in the age and metallicity of disc stars. We focus our analysis on one particular galaxy selected because of its overall resemblance to the Milky Way: most stars at lookback time are in a well defined, coplanar, centrifugally-supported disc component; it has formed stars throughout its history; it has had no recent major mergers; has no major morphological peculiarities; has had a quiet recent merger history; and has formed most of its stars in-situ.
Some radial and vertical trends in the simulated galaxy resemble those of the Milky Way. In particular, average stellar ages and metallicities decrease with increasing radius in the disc midplane; at fixed , metallicity decreases and age increases with increasing ; and the disc exhibits a pronounced ‘flare’ (i.e., the scaleheight of stars of fixed age or metallicity increases monotonically outwards).
We trace the origin of these gradients to the properties of the parent gaseous disc, which first assembles into a thick, slowly-rotating, flared structure and then gradually thins down, settles, and cools kinematically as gas turns into stars and enriches itself in the process. The disc is denser near the centre, and therefore enrichment and thinning proceed more rapidly there than in the outskirts.
Stars inherit the properties of the gas disc at the time of their formation, and evolve little thereafter. In other words, the resulting gradients are predominantly imprinted at birth, and are not the result of secular evolutionary processes such as radial migration or disc instabilities.
The similarity in the vertical structure of stars at late times and that of the gas at the time of their formation results because the gaseous disc is not in thermal hydrostatic equilibrium. Rather, its vertical structure is largely set by the bulk motions of star forming gas clouds, which are in turn induced and sustained by the feedback energy of evolving stars and, possibly, energy injected by continuous gas accretion. As a result, the timescales of star formation, enrichment, and equilibration are all intertwined, and leave behind clear radial and vertical gradients in the gaseous disc and its descendent stars.
Our results suggest that some trends that are often ascribed to secular evolutionary processes may actually be the result of the gradual equilibration process of a star forming, self-enriching disc. This should not be taken to imply that secular evolution plays no role in the establishment of such gradients in the case of the Milky Way, where the bar and spiral patterns, which are lacking in our simulated disc) are expected to lead to some radial and vertical redistribution of stars after formation. Rather, our results simply demonstrate that such gradients, per se, should not be regarded as necessarily due to secular evolution.
Our simulation is the latest of a number of cosmological simulations of disc galaxy formation that argue that the conditions ‘at birth’ play a critical role in our understanding of the origin of Galactic gradients (see; e.g., Brook et al., 2012; Bird et al., 2013; Stinson et al., 2013; Grand et al., 2015; Miranda et al., 2016; Ma et al., 2017; Grand et al., 2017, and references therein). Although the simulations differ in their details, the qualitative scenario they recount is common to all: a feedback-thickened, star-forming, flared gaseous disc that gradually turns itself into stars as it equilibrates and settles down, results in a disc galaxy that resembles the Milky Way in a surprising number of ways. These simulations also show that internal evolutionary processes (such as migration and vertical thickening) do matter, but are neither responsible for the main structural properties of the galaxies, nor for the origin of their gradients.
This scenario offers a theoretical template that may be used to interpret observations, such as those that will result from ongoing and upcoming surveys of the Milky Way. Its crucial prediction is that locally the timescales of star formation, enrichment, equilibration, and thinning should all be connected through a simple physical model of feedback and accretion. The relations between these timescales and their observational consequences need to be spelled out in more detail by future work to yield falsifiable predictions that may be used to assess the validity and general applicability of this scenario.
The research was supported in part by the Science and Technology Facilities Council Consolidated Grant (ST/F001166/1), and the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 278594-GasAroundGalaxies. CSF acknowledges ERC Advanced Grant 267291 COSMIWAY. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC system was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. This research has made use of NASA’s Astrophysics Data System.
We present in this brief appendix, for the sake of completeness and reference, the circular velocity profile of the galaxy at various times (Fig. 16), as well as the vertical density profile of star particles at a cylindrical radius kpc (Fig. 17). The former shows that our galaxy has a reasonably flat circular velocity profile, as commonly observed for spiral galaxies. The latter shows that, as for the Milky Way (Gilmore & Reid, 1983), the vertical density profile of stars is not well approximated by a simple exponential law, but can be described by the superposition of a ‘thick’ and a ‘thin’ disc.
- Anders et al. (2014) Anders F., et al., 2014, A&A, 564, A115
- Benitez-Llambay et al. (2017) Benitez-Llambay A., Navarro J. F., Frenk C. S., Ludlow A. D., 2017, preprint, (arXiv:1707.08046)
- Bensby et al. (2005) Bensby T., Feltzing S., Lundström I., Ilyin I., 2005, A&A, 433, 185
- Bergemann et al. (2014) Bergemann M., et al., 2014, A&A, 565, A89
- Bird et al. (2013) Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., Callegari S., Mayer L., Madau P., 2013, ApJ, 773, 43
- Boeche et al. (2013) Boeche C., et al., 2013, A&A, 559, A59
- Bournaud et al. (2009) Bournaud F., Elmegreen B. G., Martig M., 2009, ApJ, 707, L1
- Bovy et al. (2016) Bovy J., Rix H.-W., Schlafly E. F., Nidever D. L., Holtzman J. A., Shetrone M., Beers T. C., 2016, ApJ, 823, 30
- Brook et al. (2012) Brook C. B., et al., 2012, MNRAS, 426, 690
- Casagrande et al. (2016) Casagrande L., et al., 2016, MNRAS, 455, 987
- Chiappini et al. (2001) Chiappini C., Matteucci F., Romano D., 2001, ApJ, 554, 1044
- Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
- Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
- Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
- Eggen et al. (1962) Eggen O. J., Lynden-Bell D., Sandage A. R., 1962, ApJ, 136, 748
- Fattahi et al. (2016) Fattahi A., et al., 2016, MNRAS, 457, 844
- Freeman & Bland-Hawthorn (2002) Freeman K., Bland-Hawthorn J., 2002, ARA&A, 40, 487
- Freudenburg et al. (2016) Freudenburg J. K. C., Weinberg D. H., Hayden M. R., Holtzman J. A., 2016, preprint, (arXiv:1608.06342)
- Genovali et al. (2014) Genovali K., et al., 2014, A&A, 566, A37
- Gilmore & Reid (1983) Gilmore G., Reid N., 1983, MNRAS, 202, 1025
- Grand et al. (2015) Grand R. J. J., Kawata D., Cropper M., 2015, MNRAS, 447, 4018
- Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
- Hayden et al. (2014) Hayden M. R., et al., 2014, AJ, 147, 116
- Haywood et al. (2013) Haywood M., Di Matteo P., Lehnert M. D., Katz D., Gómez A., 2013, A&A, 560, A109
- Holmberg et al. (2009) Holmberg J., Nordström B., Andersen J., 2009, A&A, 501, 941
- Kalberla & Dedes (2008) Kalberla P. M. W., Dedes L., 2008, A&A, 487, 951
- Kawata et al. (2017) Kawata D., Grand R. J. J., Gibson B. K., Casagrande L., Hunt J. A. S., Brook C. B., 2017, MNRAS, 464, 702
- Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
- Kubryk et al. (2013) Kubryk M., Prantzos N., Athanassoula E., 2013, MNRAS, 436, 1479
- Loebman et al. (2011) Loebman S. R., Roškar R., Debattista V. P., Ivezić Ž., Quinn T. R., Wadsley J., 2011, ApJ, 737, 8
- Ma et al. (2017) Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, MNRAS, 467, 2430
- Martig et al. (2016) Martig M., Minchev I., Ness M., Fouesneau M., Rix H.-W., 2016, ApJ, 831, 139
- Matteucci & Francois (1989) Matteucci F., Francois P., 1989, MNRAS, 239, 885
- Mikolaitis et al. (2014) Mikolaitis Š., et al., 2014, A&A, 572, A33
- Minchev et al. (2012) Minchev I., Famaey B., Quillen A. C., Dehnen W., Martig M., Siebert A., 2012, A&A, 548, A127
- Minchev et al. (2013) Minchev I., Chiappini C., Martig M., 2013, A&A, 558, A9
- Minchev et al. (2014) Minchev I., Chiappini C., Martig M., 2014, A&A, 572, A92
- Miranda et al. (2016) Miranda M. S., et al., 2016, A&A, 587, A10
- Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Navarro et al. (2011) Navarro J. F., Abadi M. G., Venn K. A., Freeman K. C., Anguiano B., 2011, MNRAS, 412, 1203
- Ness et al. (2016) Ness M., Hogg D. W., Rix H.-W., Martig M., Pinsonneault M. H., Ho A. Y. Q., 2016, ApJ, 823, 114
- Nordström et al. (2004) Nordström B., et al., 2004, A&A, 418, 989
- Quillen & Garnett (2001) Quillen A. C., Garnett D. R., 2001, in Funes J. G., Corsini E. M., eds, Astronomical Society of the Pacific Conference Series Vol. 230, Galaxy Disks and Disk Galaxies. pp 87–88
- Recio-Blanco et al. (2014) Recio-Blanco A., et al., 2014, A&A, 567, A5
- Rix & Bovy (2013) Rix H.-W., Bovy J., 2013, A&ARv, 21, 61
- Roškar et al. (2013) Roškar R., Debattista V. P., Loebman S. R., 2013, MNRAS, 433, 976
- Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
- Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
- Schönrich & Binney (2009) Schönrich R., Binney J., 2009, MNRAS, 399, 1145
- Schönrich & McMillan (2017) Schönrich R., McMillan P. J., 2017, MNRAS, 467, 1154
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Stinson et al. (2013) Stinson G. S., et al., 2013, MNRAS, 436, 625
- Venn et al. (2004) Venn K. A., Irwin M., Shetrone M. D., Tout C. A., Hill V., Tolstoy E., 2004, AJ, 128, 1177
- Vera-Ciro et al. (2014) Vera-Ciro C., D’Onghia E., Navarro J., Abadi M., 2014, ApJ, 794, 173
- Vera-Ciro et al. (2016) Vera-Ciro C., D’Onghia E., Navarro J. F., 2016, ApJ, 833, 42
- Wielen (1977) Wielen R., 1977, A&A, 60, 263