The imprint of clump formation at high redshift. I. A disc alpha-abundance dichotomy

TL;DR

This paper proposes that the chemical dichotomy in the Milky Way's disc arises from early gas-rich disc fragmentation into clumps, which leads to rapid enrichment and the formation of a thick disc, explaining observed chemical bimodalities.

Contribution

The study introduces a model linking high-redshift clump formation to the Milky Way's chemical dichotomy, providing a natural explanation for the thick and thin disc properties.

Findings

Clump formation at high redshift explains the chemical dichotomy.

The model reproduces properties of the Milky Way's thick disc.

Chemical bimodalities are likely common in massive galaxies.

Abstract

The disc structure of the Milky Way is marked by a chemical dichotomy, with high-alpha and low-alpha abundance sequences, traditionally identified with the geometric thick and thin discs. This identification is aided by the old ages of the high-alpha stars, and lower average ages of the low-alpha ones. Recent large scale surveys such as APOGEE have provided a wealth of data on this chemical structure, including showing that an identification of chemical and geometric thick discs is not exact, but the origin of the chemical dichotomy has remained unclear. Here we demonstrate that a dichotomy arises naturally if the early gas-rich disc fragments, leading to some fraction of the star formation occurring in clumps of the type observed in high-redshift galaxies. These clumps have high star formation rate density. They, therefore, enrich rapidly, moving from the low-alpha to the high-alpha sequence, while more distributed star formation produces the low-alpha sequence. We demonstrate that this model produces a chemically-defined thick disc that has many of the properties of the Milky Way's thick disc. Because clump formation is common in high redshift galaxies, we predict that chemical bimodalities are common in massive galaxies.

Equations

$$\mbox{$\rm[O/Fe]$}<\begin{cases}0.3&\text{if \mbox{$\rm[Fe/H]$}\ $<-0.8$}\\ -0.3913\mbox{$\rm[Fe/H]$}-0.0130&\text{otherwise.}\end{cases}$$

$$C_{\lambda }^{gas}=C_{\lambda ,MW}^{gas}\frac{Z}{Z_{MW}}$$

Full text

1. Introduction

The Milky Way (MW) is important to our understanding of the structure and evolution of galactic discs. The geometical thick disc was originally discovered in the MW by decomposing the vertical star count profile into two components, the thin and thick discs (Yoshii, 1982; Gilmore & Reid, 1983). Spectroscopic studies of Solar neighborhood stars found that the thick disc is composed of mainly older ($\sim$12–7 Gyr) stars that are enhanced in their $\alpha$ abundances, while thin disc stars are younger ($<$7 Gyr) and have near-Solar $\alpha$ abundances (e.g., Fuhrmann, 1998; Haywood et al., 2013). Bovy et al. (2012) studied the spatial distribution of mono-abundance populations which showed smooth variations of scale-height with abundance and suggested that there was no distinct separation of “geometrical” discs (Bovy et al., 2012), nor, perhaps, in chemistry. Even so, the “chemical” separation of high and low-$\alpha$ populations persists in recent high-resolution abundance surveys such as APOGEE (e.g., Anders et al., 2014; Nidever et al., 2014; Hayden et al., 2015) and unbiased Solar neighborhood samples (Adibekyan et al., 2012).

It is now believed that the chemically and geometrically defined thick discs are different entities. The chemically defined thick disc has a shorter scale-length than the geometric thick disc (Robin et al., 1996; Ojha, 2001; Jurić et al., 2008; Bensby et al., 2011; Bovy et al., 2012; Martig et al., 2016), although there is still some uncertainty about the scale-lengths (Robin et al., 2014). Old, $\rm\alpha$-enhanced stars are concentrated to the inner MW Minchev et al. (2015). Stars in the chemical thick disc are generally old, whereas geometric thick disc stars can have median ages as young as 5\mbox{>{\rm Gyr}} in the outer disc (Martig et al., 2016).

Several explanations have been proposed for the $\alpha$-bimodality. The “two-infall” model (Chiappini et al., 1997; Chiappini, 2009) suggests that the early MW had a high star formation rate (short gas consumption timescale) that changed dramatically around 8 Gyr when there was a drop in the overall MW star formation rate and a large infall of pristine gas that diluated the overall metallicity of the Milky Way and gave rise to the thin disc, low-$\alpha$ population. Since this model requires a particular event to occur (the large infall of pristine gas at a particular time), it suggests that the bimodality might not be a generic feature of all galaxies.

On the other hand, a “superposition” model (which is similar to the scenario proposed in Schönrich & Binney 2009) postulates that the $\alpha$-bimodality is a generic consequence of the chemical and dynamical processes in a spiral disc galaxy. The chemical evolution in $\alpha$ abundance is quite rapid at early times but slows down at later times when most of the stars are formed. In addition, a radial gradient in outflow rate exists in the disc due to the decrease in the energy required to unbind the gas. The large outflow in the outer disc inhibits the gas from attaining a high metallicity while in the inner galaxy chemical evolution is able to advance fairly unimpeded. Finally, the dynamical process of radial migration then moves stars away from their birth radii and thereby smoothes out the separate chemical evolution tracks. The $\alpha$-bimodality is created by the fact that the chemical evolution proceeds most quickly as SNIa’s turn on and evolution is moved from high-$\rm\alpha$ to low-$\rm\alpha$ thus producing fewer stars in the $\alpha$ “valley”. Producing chemical bimodality via these processes requires fine-tuning; indeed simulations which include these processes have failed to produce any chemical bimodality (e.g. Loebman et al., 2011; Minchev et al., 2013), which instead only produce a thick band in the ([$\alpha$/Fe], [Fe/H]) space representing the overall chemical evolution of the galaxy with small variations with radius. Nidever et al. (2014) presents a detailed discussion and example figures of the chemical evolution of these two scenarios.

Nevertheless, in recent years some simulations have begun to find hints of $\alpha$-bimodalities. The EAGLE simulations have found $\alpha$-bimodalities produced by an early centralized starbust (important for the inner galaxy) and disc contraction (important for the outer galaxy) that produce chemical bimodalities somewhat similar to the MW in one of six simulated galaxies (Grand et al., 2017). Others have found that a late infall of a gas-rich galaxy can produce some metal-poor, $\alpha$-poor stars that give rise to a bimodality (Snaith et al., 2016; Mackereth et al., 2018). However, these bimodalities are weaker than the one observed in the MW and the age distributions are not always consistent with the MW, suggesting they are not the main mechanism at work in the MW.

One option that has not been previously considered is the influence of star forming “clumps”. Clumps are commonly found in Hubble Space Telescope (HST) images of high redshift galaxies (e.g. Elmegreen & Elmegreen, 2005; Ravindranath et al., 2006; Elmegreen et al., 2007; Förster Schreiber et al., 2011; Genzel et al., 2011; Guo et al., 2012, 2015). Such clumps were first identified by eye in HST images (e.g. Cowie et al., 1995; van den Bergh et al., 1996) and are common in MW progenitors at high redshift (present in $\sim 55\%$ of such galaxies at a redshift of 3, Guo et al., 2015). They appear \sim 1\mbox{>{\rm kpc}} in size when observed at low resolution, but at higher resolution appear to have sizes 100\mbox{>{\rm pc}} to 500\mbox{>{\rm pc}} (Livermore et al., 2012, 2015; Fisher et al., 2017; Cava et al., 2018). Clumps can contribute up to $20\%$ of the integrated star-formation rate of a galaxy, but can be less massive (up to $7\%$ of the total mass, Wuyts et al., 2012). It is proposed that they form a bulge in less than a half a Gyr (Dekel et al., 2009). They have also been proposed for the origin of the geometrical thick disc (Bournaud et al., 2009).

In this paper, we demonstrate that discs evolving with significant clump formation produce a bimodality in the [$\alpha$/Fe]-[Fe/H] space by virtue of the very rapid star formation which occurs in clumps. We show that the resulting high-$\alpha$ track in chemical space is comprised of old stars and suggest that this might be the process forming the $\alpha$-bimodality in the MW. The layout of this paper is as follows. In Section 2 we describe the simulation we use in this paper. Section 3 then presents the results of the simulation. We present our conclusions in Section 4.

2. Simulation

We set up a spherical Navarro-Frenk-White (Navarro et al., 1997) halo with an equilibrium spherical gas distribution containing 10% of the total mass and following the same density distribution. The starting halo is set up as described in Roškar et al. (2008): it has a mass within the virial radius ($r_{200}\simeq 200$ kpc) of $10^{12}\>{\rm M_{\odot}}$. A temperature gradient ensures an initial gas pressure equilibrium for an adiabatic equation of state. Gas velocities are initialized to give a spin parameter of $\lambda=0.065$ (Bullock et al., 2001; Macciò et al., 2007), with specific angular momentum $j\propto R$, where $R$ is the cylindrical radius. Both the gas corona and the dark matter halo are comprised of $10^{6}$ particles; gas particles initially have masses $1.4\times 10^{5}\>{\rm M_{\odot}}$ and force softening 50 pc. Dark matter particles instead come in two mass flavours ($10^{6}\>{\rm M_{\odot}}$ and $3.5\times 10^{6}\>{\rm M_{\odot}}$ inside and outside 200 kpc, respectively), with a softening of 100 pc. The stars in the simulation form self-consistently out of cooling gas, inheriting the softening and chemistry of the parent gas particle.

We evolve the simulation for 10 Gyr with gasoline (Wadsley et al., 2004, 2017), the smooth particle hydrodynamics (SPH) extension of the $N$-body tree-code pkdgrav (Stadel, 2001). Gas cooling includes the metal-line cooling implementation of Shen et al. (2010); in order to prevent the cooling from dropping below the resolution of our simulation, we set a pressure floor on gas particles $p_{floor}=3G\epsilon^{2}\rho^{2}$, where $G$ is Newton’s gravitational constant, $\epsilon$ is the softening length and $\rho$ is the gas particle’s density (Agertz et al., 2009). Gas cools and settles into a disc; once the density exceeds $1\mathrm{cm}^{-3}$ and the temperature drops below 15,000 K, star formation and supernova feedback cycles are initiated as described in Stinson et al. (2006). A consequence of the metal-cooling is that gas can reach lower temperatures where it is prone to fragmentation and clump formation. Supernovae feedback couples 10% of the $10^{51}$ erg per supernova to the interstellar medium as thermal energy. We include the effects of turbulent diffusion (Shen et al., 2010) allowing the gas to mix, reducing the scatter in the age-metallicity relation (Pilkington et al., 2012). We use a base time step of $\Delta t=10$ Myr with timesteps refined such that $\delta t=\Delta t/2^{n}<\eta\sqrt{\epsilon/a_{g}}$, where $a_{g}$ is the acceleration at a particle’s position and the refinement parameter $\eta=0.175$. We set the opening angle of the tree-code gravity calculation to $\theta=0.7$. In addition the time step of gas particles satisfies the condition $\delta t_{\mathrm{gas}}=\eta_{\mathrm{courant}}h/[(1+\alpha)c+\beta\mu_{\mathrm{max}}]$, where $\eta_{\mathrm{courant}}=0.4$, $h$ is the SPH smoothing length set over the nearest 32 particles, $\alpha$ and $\beta$ are respectively the linear and quadratic viscosity coefficients and $\mu_{\mathrm{max}}$ is described in Wadsley et al. (2004). Star particles represent single stellar populations with a Miller-Scalo initial mass function. SNII and SNIa yields of Oxygen and Iron are taken from Raiteri et al. (1996). As in Raiteri et al. (1996), Padova stellar lifetimes are used to determine SNII rates, and SNIa rates are determined from those same lifetimes in a binary evolution model.

2.1. Global properties of the final system

Figure 1 presents the face-on and edge-on views of the model at the end of the simulation, showing it to be a reasonable example of a disc galaxy. Figure 2 presents the final rotation curve and velocity dispersions of the system. Whilst our simulation is not designed to match any specific galaxy, the rotation curve of the final system has a circular velocity at 8\mbox{>{\rm kpc}} of 242~{}\mbox{>{\rm km,s^{-1}}}, making it comparable to the Milky Way. Its structural and kinematic properties are somewhat different from the Milky Way’s, for instance the model is unbarred, but is otherwise similar enough that it can be used to understand the origin of trends observed in the Milky Way.

3. Results

3.1. Chemical bimodality

Figure 3 shows the final (10\mbox{>{\rm Gyr}}) distribution of stars in the $\rm[O/Fe]$-$\rm[Fe/H]$ (chemical) plane at different radii and heights. As in the APOGEE survey (Anders et al., 2014; Nidever et al., 2014; Hayden et al., 2015) and other datasets (e.g., Fuhrmann, 1998; Bensby et al., 2005; Reddy et al., 2006; Adibekyan et al., 2012), we find two distinct sequences, a high-$\rm\alpha$ and a low-$\rm\alpha$ one. The low-$\rm\alpha$ sequence dominates over the high-$\rm\alpha$ one close to the mid-plane, while the high-$\rm\alpha$ sequence is more prominent at larger heights. The two sequences are bridged by stripes which are not observed in the MW. If we convolve the $\rm[O/Fe]$ and $\rm[Fe/H]$ of stellar particles in the Solar neighbourhood (defined throughout this paper as $7<R/kpc<9$) with the APOGEE observational errors, \mbox{\sigma_{\mbox{$\rm[Fe/H]$}}}=0.1 dex and \mbox{\sigma_{\mbox{$\rm[O/Fe]$}}}=0.03 dex (Nidever et al., 2014) the stripes visible in Figure 3 are masked, without erasing the chemical bimodality. Figure 4 directly shows the bimodality in $\rm[O/Fe]$ at the Solar neighbourhood for the cut through the chemical space -0.5<\mbox{\rm[Fe/H]}<-0.4. The distribution is bimodal, with a very clearly defined minimum, even when the distribution is convolved with the observational errors, as it must be.

We define two lines by eye to separate the chemical space into low-$\rm\alpha$ and high-$\rm\alpha$ sequences:

[TABLE]

Throughout this paper, we refer to stars above this locus as the high-$\rm\alpha$ sequence and those below as the low-$\rm\alpha$ sequence. These lines are overplotted on Figure 3. Across the entire galaxy, the mass of the high-$\rm\alpha$ sequence is $1.7\times 10^{10}\>{\rm M_{\odot}}$ out of a total stellar mass of $6.23\times 10^{10}\>{\rm M_{\odot}}$, *i.e. * the high-$\rm\alpha$ sequence contains $28\%$ of the stellar mass, which decreases to below $20\%$ if only the mass outside 3\mbox{>{\rm kpc}} is considered. In comparison, Snaith et al. (2014) estimate that the entire chemically-defined thick disc of the MW constitutes $47\%$ of the stellar mass.

Figure 5 compares the chemical distribution in the MW as observed by APOGEE throughout a large extent of the galaxy (DR14 red giants, Holtzman et al., 2018) with that of the model convolved with the APOGEE observational errors. The lower sequence has the opposite curvature compared with the model but this may be sensitive to the star formation history and the specific elements used (e.g. Ramírez et al., 2007), as is the scale of the vertical axis. The $\rm[\alpha/Fe]$ defined by APOGEE includes oxygen but also other $\rm\alpha$-elements; when we consider just $\rm[O/Fe]$ in the APOGEE data, although the two sequences are noisier, the curvature of the low-$\rm\alpha$ track is similar to the simulation’s at the low $\rm[Fe/H]$ end, though still opposite at high $\rm[Fe/H]$. Nonetheless, an unambiguous and sustained separation into two tracks is readily apparent in the model, as in the MW. Additionally, the extent and curvature of the upper sequence is similar to that in the APOGEE data. The facts that the high-$\rm\alpha$ sequence covers nearly the same $\rm[Fe/H]$ range as the low-$\rm\alpha$ sequence, and that the two sequences approach each other at the high $\rm[Fe/H]$ end, are especially intriguing constraints on the formation of the bimodality because it implies that it did not form in a single discrete event such as infall of pristine gas (e.g. Mackereth et al., 2018). The extension of the high-$\rm\alpha$ sequence to lower metallicities than the low-$\rm\alpha$ also matches the MW. Still, the model and the MW do not match in detail; one of the more significant differences is the larger separation between the high-$\rm\alpha$ and the low-$\rm\alpha$ sequences. This may be a problem of either the yields used (also supported by the relative vertical offset of the model and the MW) or the star formation rates that produce the high-$\rm\alpha$ sequence.

3.2. Present day properties of the high-$\rm\alpha$ sequence

The spatial distribution of the high-$\rm\alpha$ stars is also evident in Figure 3. High-$\rm\alpha$ stars are more prominent at large heights at the present time. The high-$\rm\alpha$ distribution is also more centrally concentrated than the low-$\rm\alpha$ ones; in the MW Bovy et al. (2012) found shorter scale-lengths for low-$\rm\alpha$ stars while Hayden et al. (2015) found little evidence of the high-$\rm\alpha$ sequence outside the Solar neighbourhood. Figure 6 shows the present day root mean square vertical height of stars across the chemical space. The high-$\rm\alpha$ stars are distributed in a thicker disc than the low-$\rm\alpha$ stars outside the bulge region (e.g. Bovy et al., 2012). The variation of height is continuous across the chemical space, as found in the MW by Bovy et al. (2012). This important constraint suggests that the high-$\rm\alpha$ stars are not the product of a single violent event which had suddenly heated the disc vertically. In spite of this continuous height variation the vertical stellar density profile is very clearly double-exponential, as can be seen in Fig. 8 for the Solar neighbourhood (other regions in the disc look qualitatively similar). Fitting two exponentials to the Solar neighbourhood vertical profile of the model gives scale-heights h_{1}=240\pm 8\mbox{>{\rm pc}} and h_{2}=1083\pm 52\mbox{>{\rm pc}}; in comparison Jurić et al. (2008) measured h_{1}=300\pm 50\mbox{>{\rm pc}} and h_{2}=900\pm 180\mbox{>{\rm pc}} in the Solar neighbourhood. The geometric thin and thick discs have a transition at about 1\mbox{>{\rm kpc}}, which is comparable to the transition in the MW (e.g. Gilmore & Reid, 1983).

Nidever et al. (2014) and Hayden et al. (2015) showed that the metallicity distribution function (MDF) of $\rm\alpha$-enhanced stars is essentially independent of the location within the MW’s disc. Figure 7 presents the MDFs of stars with \mbox{\rm[O/Fe]}>0.2 (note this $\rm\alpha$-enhanced population is not identical to the high-$\rm\alpha$ sequence, but is chosen in the same spirit as in Nidever et al., 2014), in the simulation for various radial and height cuts. The overall shape and width of these MDFs are similar in all regions considered. Nidever et al. (2014) and Hayden et al. (2015) suggested that the similarity of the MDFs of $\rm\alpha$-enhanced stars indicates these stars formed from well-mixed gas. Loebman et al. (2016) used simulations to demonstrate that the similarity of the MDF of $\rm\alpha$-enhanced stars is due to them having formed in a relatively small region in the inner Galaxy over a short time interval (\lesssim 2\mbox{>{\rm Gyr}}) and were subsequently migrated throughout the disc. Figure 9 presents the average age of stars in the chemical space. The high-$\rm\alpha$ stars are old, having mean age >8\mbox{>{\rm Gyr}}, whilst the low-$\rm\alpha$ stars, particularly those at high $\rm[Fe/H]$, are predominantly young to intermediate age. Figure 10 presents the mean formation radius, $R_{\mathrm{f}orm}$, across the chemical space. The high-$\rm\alpha$ stars largely formed within the inner 2\mbox{>{\rm kpc}}, while the low-$\rm\alpha$ stars formed over a larger radial region. In agreement with Loebman et al. (2016), therefore, we find that the high-$\rm\alpha$ stars, and especially the $\rm\alpha$-enhanced stars, formed over a short time interval in a relatively restricted volume.

Figure 11 shows the star formation history (SFH) of the high-$\rm\alpha$ and low-$\rm\alpha$ stars. The two populations overlap at old ages, with the SFH of the high-$\rm\alpha$ stars peaking very early, persisting to \sim 4\mbox{>{\rm Gyr}}, but largely over by 2\mbox{>{\rm Gyr}}. The SFH of the low-$\rm\alpha$ stars is relatively flat after 1\mbox{>{\rm Gyr}}, becoming the main mode of star formation thereafter. The early formation of the majority of the high-$\rm\alpha$ stars ensures that they are more centrally concentrated then the subsequent low-$\rm\alpha$ disc stars, in good agreement with what is found in the MW (Bovy et al., 2012; Bovy et al., 2016).

Figure 12 compares the stellar ages in the MW from the APOGEE survey and the simulation. The left panel, showing the ages in the simulation, uses star particles in the simulation with |z|<1\mbox{>{\rm kpc}}, and across Galactic radii similar to the red clump (3<x<13\mbox{>{\rm kpc}} and |y|>5\mbox{>{\rm kpc}}). The right panel uses $\sim 20,000$ red clump stars in APOGEE with ages from Ness et al. (2016). While ages of the high-$\rm\alpha$ stars are similar between the model and the MW, the two are quite different in detail. The biggest difference is at low metallicity, where the model has young stars only at low $\rm\alpha$. The model’s chemical enrichment history therefore is different from that of the MW. This makes the similarities in the bimodality even more remarkable.

3.3. Two modes of star formation

We now investigate in greater detail the origin of the high-$\rm\alpha$ sequence. High-$\rm\alpha$ stars are generally thought to form at high star formation rate (SFR); we compute the local SFR density, $\Sigma_{\mathrm{SFR}}$, at birth of each star particle from the distribution of stars forming in each 10\mbox{>{\rm Myr}} interval. We estimate the density of these new stars by recovering their positions at the beginning of the time interval using the mean disc rotation velocity at their radius. We then calculate the density of the particles in the disc plane by binning the particles in Cartesian coordinates with a bin size of 100 pc $\times$ 100 pc. The grid is then smoothed with a Gaussian of FWHM=200\mbox{>{\rm pc}} and the values are converted to physical units of \>{\rm M_{\odot}}\mbox{>{\rm yr}}^{-1}\mbox{>{\rm kpc}}^{-2}. Finally we assign local $\Sigma_{\mathrm{SFR}}$ values to each star particle on the basis of the bin within which it falls.

Figure 13 presents snapshots of the local star formation rate computed this way over the course of the simulation. At early times a considerable amount of star formation occurs in clumps, which have a high SFR density. As in the two-infall model of Chiappini et al. (1997) this mode of star formation becomes negligble after a few Gyr, aside from at the bulge. At the same time as the clumpy star formation is occurring, stars are forming in the disc in a more distributed (“smooth”) configuration with much lower $\Sigma_{\mathrm{SFR}}$. This is the main mode of star formation after 3\mbox{>{\rm Gyr}}.

Figure 14 shows $\Sigma_{\mathrm{SFR}}$ of star particles at the time of formation across the chemical plane. The high-$\rm\alpha$ sequence has a significantly higher $\Sigma_{\mathrm{SFR}}$, which results in being $\rm\alpha$-rich, while the low-$\rm\alpha$ stars have much lower local $\Sigma_{\mathrm{SFR}}$ and chemically evolve more slowly. Unlike Figure 6, Figure 14 has an abrupt transition between the high-$\rm\alpha$ and the low-$\rm\alpha$ sequences. Figure 15 presents histograms of $\Sigma_{\mathrm{SFR}}$ at formation, showing that there is a bimodality in $\Sigma_{\mathrm{SFR}}$. The two modes of SF differ by a factor of roughly 100 in $\Sigma_{\mathrm{SFR}}$. The high-$\rm\alpha$ sequence forms via the high-SFR density, clumpy star formation mode, while the low-$\rm\alpha$ sequence forms via the low-SFR density, smooth star formation mode that is dominant at later times. As already shown in Figure 11, the low-$\rm\alpha$ mode of star formation is already present at early times at the same time as the high star formation rate, high-$\rm\alpha$ mode. This simultaneous formation of both $\rm\alpha$-sequences, via the dual-modes of star formation, appears to be in agreement with the recent observational results of Silva Aguirre et al. (2017) and Hayden et al. (2017).

Figure 16 shows the evolution of the radius at which stars are forming in the high-$\rm\alpha$ and low-$\rm\alpha$ sequences. The high-$\rm\alpha$ sequence is seen to derive from a decreasing number of discrete star forming clumps that sink towards the centre. The star formation in the low-$\rm\alpha$ sequence is more broadly distributed, although the effect of spirals is evident as gaps in the radial distribution at R_{f}<3\mbox{>{\rm kpc}} at certain times. These spirals are associated with the clumps, by being partly excited by them (D’Onghia et al., 2013) and also partly needing the density enhancement provided by spirals to assemble the gas needed to form clumps. Low-$\rm\alpha$ star formation is also associated with the clumps throughout a large fraction of their trajectories, but the bulk of low-$\rm\alpha$ stars form outside clumps, particularly after the first $\>{\rm Gyr}$. A similar plot to Figure 16 but separated into a low and high SFR density at \Sigma_{\mathrm{SFR}}=30\>{\rm M_{\odot}}\mbox{>{\rm kpc}}^{-2}\mbox{>{\rm yr}}^{-1} is virtually identical apart from having more star formation in the bulge, persisting for the full 10\mbox{>{\rm Gyr}}. About $43\%$ of stars form in the high SFR density component, which reduces to $31\%$ if the inner 0.5\mbox{>{\rm kpc}} is excluded. In comparison we earlier found that the high-$\rm\alpha$ sequence comprised $28\%$ of all stars. Comparing these two numbers we conclude that the majority of high-$\rm\alpha$ stars form at high $\Sigma_{\mathrm{SFR}}$ within clumps.

Figure 17 presents the evolution of two typical clumps chosen by eye. Clump 1 forms at \sim 6\mbox{>{\rm kpc}}; initially on the low-$\rm\alpha$ sequence, it quickly shoots up to the high-$\rm\alpha$ sequence as $\Sigma_{\mathrm{SFR}}$ rises, becoming more $\rm[Fe/H]$-rich along this sequence as supernovae of type Ia commence. Its orbit decays and by the time it reaches 1\mbox{>{\rm kpc}} it starts falling back on to the low-$\rm\alpha$ sequence. At this point the clump is lost in the bulge. Clump 2 also ascends from the low-$\rm\alpha$ to the high-$\rm\alpha$ sequence. Its orbital radius is rather more constant, and it evolves in $\rm[Fe/H]$ less over its shorter lifetime (300\mbox{>{\rm Myr}} versus 500\mbox{>{\rm Myr}} in clump 1). Towards the end of its life it merges with another clump and quickly drops to the centre of the galaxy, dropping from the high-$\rm\alpha$ to the low-$\rm\alpha$ sequence. These two clumps illustrate an important point: contrary to the usual interpretation, the high-$\rm\alpha$ sequence is reached from the low-$\rm\alpha$ sequence, rather than forming before it. These examples also demonstrate that the stripes connecting the low-$\rm\alpha$ to the high-$\rm\alpha$ sequence are due to the evolution of individual clumps in chemical space.

The clump mass distribution is shown in Figure 18. Clumps are detected in 100 Myr snapshots using a 10$\sigma$ threshold in the 100 pc $\times$ 100 pc binned mass map with a 1.5 kpc smoothed background removed. The total clump mass is measured by summing up the mass in a 500\mbox{>{\rm pc}} radius around the clump centre and removing the background level determined in a 1.5-2.5\mbox{>{\rm kpc}} annulus. Only clumps with well-determined masses (error $\lesssim 10\%$; to remove noise fluctuations) and galactocentric radii greater than 100\mbox{>{\rm pc}} (to remove the inner bulge) are kept and these correspond well to the visually identifiable clumps. Note that because clumps survive for \sim 500\mbox{>{\rm Myr}} before falling into the galactic centre, they are detected multiple times in our technique. The range of clump masses is $3\times 10^{7}$ to $8\times 10^{9}\>{\rm M_{\odot}}$ with a mean of $5\times 10^{8}\>{\rm M_{\odot}}$, although our detection technique is not sensitive to clumps below $\sim 3\times 10^{7}$. This mean mass is of the same order, but higher than, the median clump mass found in high-redshift galaxies observed at high resolution through gravitational lensing, $10^{8}\>{\rm M_{\odot}}$, (e.g. Cava et al., 2018, their Figure 2).

3.4. A high redshift view

Because our mass resolution results in only massive clumps being identifiable in our simulation, we now explore whether our simulation is an excessively clumpy galaxy when observed at the \sim\mbox{>{\rm kpc}} resolution typical of galaxies at high redshift. We use the ray-tracing 3D dust radiative transfer code DART-Ray (Natale et al., 2014, 2015, 2017), which self-consistently calculates the effect of dust attenuation on the stellar radiation (direct stellar light), and the subsequent dust emission (dust re-radiated stellar light). We use DART-Ray to produce face-on dust attenuated images at rest UV, B and V bands of our simulated galaxy at t=1\mbox{>{\rm Gyr}}, shortly after the peak SFR in the high-$\rm\alpha$ sequence (see Figure 11), as though it were a galaxy at $z=2$. The direct light is assumed to be strictly due to star particles, whereas emission by gas is ignored. The stellar emissivity is computed at wavelengths ranging from $0.1\mu$m to $30\mu$m for the mass, position, age and metallicity of stars in the simulation at 1\mbox{>{\rm Gyr}}. We adopt spectral energy distributions from the starburst99 spectral synthesis code assuming a Kroupa (Kroupa, 2001) initial mass function and the Padova AGB tracks. For the optical properties of grains and grain sizes we use the trust benchmark dust model, which is based on the model of Zubko et al. (2004). We assume that the gas $\rm[Fe/H]$ traces the metal abundance $Z$ according to Z=Z_{MW}10^{\mbox{\rm[Fe/H]}} where $Z_{MW}$ is the metal abundance of the Milky Way, assumed to be equal to 0.018. We assumed that the dust cross-section $C^{gas}_{\lambda}$ is proportional to the metal abundance through a linear relation:

[TABLE]

where $C^{gas}_{\lambda,MW}$ is the dust cross-section per unit gas mass derived from the trust dust model. Then the radiation transfer is computed, reaching a highest spatial resolution of 82\mbox{>{\rm pc}}.

The images from DART-Ray are noise-free and have a spatial resolution comparable to that of the simulation. To enable a better comparison with real observations, we degrade the signal-to-noise and spatial resolution of the images to that typical of high-redshift surveys (see also Snyder et al., 2015). Therefore, we generate synthetic Hubble Space Telescope Wide Field Camera 3 (HST/WFC3) images matched to the observing conditions of the CANDELS (Grogin et al., 2011) and HUDF09 (Bouwens et al., 2010) surveys, adopting the typical depth and resolution of these surveys as listed in Table 1 of Guo et al. (2013). We create synthetic HST/WFC3 images for three WFC3 filters — F098M, F125W and F160W. At $z=2$, these trace 0.3 - 0.5 $\mu$m in the rest-frame, i.e., UV-, B- and V-band. To generate the synthetic images, we first rescale the surface brightness of each pixel as $I_{observed}=I_{emitted}\times(1+z)^{-4}$ to account for cosmological surface brightness dimming, ignoring k-corrections. The angular sizes of the pixels are then scaled with the redshift and the pixels are rebinned to 0.″6/pixel — the typical pixel scale of drizzled HST WFC3 observations. Next, we convolve the images with a 2D Gaussian kernel to simulate the WFC3 point spread function, FWHM 0.″13 for F098M and 0.″16 for F125W and F160W. Finally, we add shot noise to each pixel to match the typical depth of each survey — 5$\sigma$ limits of $\sim$ 27.5 and 28.3 mag/arcsec2 for the wide and deep portions of the CANDELS survey, respectively, and 29.7 mag/arcsec2 for the HUDF09.

Figure 19 shows the resulting images, at the original resolution, and convolved with HST resolution, with surface brightness dimming included and with noise added as in the CANDELS survey (Grogin et al., 2011) and the Hubble Ultra Deep Field (HUDF, Beckwith et al., 2006). We use these images to identify clumps. We first fit a single-component Sérsic (Sérsic, 1968) model to the F098M image (0.3 $\mu$m rest wavelength) and search in the residual image for off-centre clumps. Following the definition of Guo et al. (2018), we identify one clump containing $\sim 14\%$ of the total UV flux in the CANDELS-wide image. This clump, visible in the top right of the galaxy, is comprised of 4 sub-clumps which merge together at this resolution. It is at 4\mbox{>{\rm kpc}} from the galaxy centre at this time and has a F098M-F160W ($U-V$ rest wavelength) colour of $0.34\pm 0.15$, which is significantly bluer than the rest of the galaxy, but within the range of observed clumps (e.g. Guo et al., 2018). Our analysis of the mock deep and ultra-deep images identifies three additional mini-clumps, having less than $1\%$ of the total UV flux and would not be considered clumps in the usual definition. Their $U-V$ colours range from $0.2$ to $0.8$.

We have repeated this analysis at t=2\mbox{>{\rm Gyr}} and at t=3\mbox{>{\rm Gyr}} (still pretending that the galaxy is at $z=2$). At t=2\mbox{>{\rm Gyr}} we find 2 clumps, each contributing $\sim 10\%$ of the UV flux, while we find only 1 clump and only at the UDF depth for t=3\mbox{>{\rm Gyr}}, contributing $\sim 9\%$ of the UV flux.

This analysis shows that our simulation is not plagued by an excess clump formation and would not be viewed as unusual if found amongst galaxies at $z\sim 2$.

4. Discussion and Conclusions

We have shown that the early, gas-rich phase of galaxy formation, which supports two modes of star formation, distributed and clumpy, as generally observed by HST at high redshift, naturally gives rise to a bimodal distribution in chemical space, similar to what is seen in the Milky Way. The clumpy mode of star formation has a star formation rate density about 100 times higher than the distributed mode and leads to rapid self-enrichment with high $\rm\alpha$-abundance, whereas the distributed star formation is responsible for less $\rm\alpha$-enchanced stars. As the stellar mass of the galaxy builds up, the gas fraction drops, leading to less star formation taking place in clumps, which naturally shuts off most of the formation of high-$\rm\alpha$ stars (with star formation in the bulge being the exception).

While the vertical distribution of stars in Figure 6 and Figure 8 is reminiscent of the thick disk in the Milky Way, the density distribution of high-$\rm\alpha$ stars flares, which means there is not a direct one-to-one match between our model and the Milky Way. The high-$\rm\alpha$ vertical scale-height is consistent with Milky Way observations at larger radii, but it is too small in the inner galaxy to match the MW. This inconsistency could be a result of the idealised nature of our simulation, which is isolated and thus excludes the effects of galaxy accretion that would dynamically heat the disk (and increase the scale-height) further and is most important at early times, when most of the high-$\rm\alpha$ stars are forming. The inclusion of these effects would bring our model into closer agreement with the observations. Moreover, the age of stars in chemical space are qualitatively different from those in the Milky Way, although this affects the low-$\rm\alpha$ sequence (the “thin” disc) more than the high-$\rm\alpha$ one. The main result we highlight here is the fact that the current simulation is naturally able to produce the chemical bimodality that is seen in the surveys but which, to date, has been challenging to produce in simulations.

Clumps are thought to form via gravitational instability in gas rich discs (Noguchi, 1999; Immeli et al., 2004; Elmegreen et al., 2008; Inoue et al., 2016), although some clumps may have an “ex-situ” (merger) origin (e.g. Puech et al., 2009; Hopkins et al., 2013; Mandelker et al., 2014). The fate of clumps in simulations has been a subject of debate, with either their destruction by supernova and/or radiative feedback (Hopkins et al., 2012; Genel et al., 2012; Buck et al., 2017; Oklopčić et al., 2017) or sinking to the centre of the galaxy to contribute to the bulge (Noguchi, 1999; Immeli et al., 2004; Bournaud et al., 2007; Genzel et al., 2008; Elmegreen et al., 2008; Dekel et al., 2009; Ceverino et al., 2010; Mandelker et al., 2017) suggested. We have found that our clumps sink to the centre of the galaxy on timescales of $200-500$ Myr. In a companion paper, we explore the consequences of clumps falling to the centre for the chemistry of the bulge, which further helps constrain this model.

We have verified that if we increase the supernova feedback in our simulations that clump formation is substantially inhibited and no chemical bimodality results. Although we use a lower feedback efficiency than is often employed in cosmological simulations, we have demonstrated that the clumps in our simulation are reasonable compared to observations, such as those in the CANDELS-wide (Grogin et al., 2011) and the HUDF (Beckwith et al., 2006), and that our model is not excessively clumpy. At present there is some uncertainty regarding the overall importance of clumps; our results suggest that even with modest clumps the chemical consequences for the MW could be significant.

Recently Mackereth et al. (2018) suggested that a bimodality in chemical space can be produced through gas accretion. Their cosmological simulations produce chemical bimodality less than half the time. In contrast, the scenario presented here will work for a large fraction of galaxies with mass comparable to the Milky Way; Wuyts et al. (2012) find that $\sim 40\%$ of galaxies of Milky Way mass are clumpy at $1.5<z<2.5$ when considering their mass distribution (rising to $60\%$ in the UV) in which case most such galaxies will have chemical bimodalities. At the moment this prediction cannot be tested but future facilities, such as the Extremely Large Telescope, Giant Magellan Telescope and the Thirty Meter Telescope, may be able to map the chemical space of stars in nearby galaxies.

4.1. Summary

The results and implications of this study are as follows:

We produce a bimodal distribution of stars in the $\rm[O/Fe]$-$\rm[Fe/H]$ space via clumpy$+$distributed star formation. The clumps self-enrich in $\alpha$-elements due to their high star formation rate density and produce the high-$\alpha$ sequence while the low-$\alpha$ sequence is produced by the distributed star formation. The clumpy mode becomes less efficient as the gas mass fraction drops, leaving an $\rm\alpha$-enhanced population that is old at the present day. The clumps produced in our simulation are in reasonable agreement with clumps observed at $z\sim 2$ (e.g. Guo et al., 2018). 2. 2.

Clumps form on the low-$\rm\alpha$ sequence, and ascend to the high-$\rm\alpha$ sequence and their star formation rate increases. Towards the ends of their lives, they tend to drop back to the low-$\rm\alpha$ sequence, although this typically happens close to the bulge. 3. 3.

The two $\rm\alpha$-sequences form simultaneously early on in the evolution of the Milky Way, resulting in overlapping ages. Recent results by Silva Aguirre et al. (2017) and Hayden et al. (2017) indicate that temporal overlap of the two sequences is seen in the observations as well. For the foreseeable future this age overlap is the most promising way to test this hypothesis further. 4. 4.

The model predicts that bimodal chemical thick discs should be common at the mass of the Milky Way, since at least $\sim 40\%$ of galaxies of Milky Way mass are clumpy at $z\sim 2$ (Wuyts et al., 2012). An important test of this scenario therefore is that chemical bimodality should be common in Milky Way mass galaxies.

The importance of clumps to the evolution of disc galaxies continues to be hotly debated (e.g. Ceverino et al., 2010; Hopkins et al., 2012; Moody et al., 2014; Tamburello et al., 2015; Mayer et al., 2016; Mandelker et al., 2017; Buck et al., 2017; Tamburello et al., 2017; Benincasa et al., 2018). The lifetimes and masses of clumps are therefore somewhat uncertain, and in simulations is dependent on resolution and sub-grid physics (e.g. Tamburello et al., 2015), and possibly also on the mode of gas accretion (Ceverino et al., 2010). However there is a growing view from high redshift observations that clumps cannot be too massive given the high degree of rotation in discs (e.g. Wuyts et al., 2012; Wisnioski et al., 2015, 2018). For instance, by seeding their models with a spectrum of velocity perturbations matching the turbulent cascade, Benincasa et al. (2018) found that likely bound masses cannot be larger than $10^{9}\>{\rm M_{\odot}}$ and typical masses $\sim 10^{6}\>{\rm M_{\odot}}$ to $\sim 10^{8}\>{\rm M_{\odot}}$. At lower mass scales, simulations of giant molecular clouds have shown that photoionizing radiation, stellar-winds, and supernova feedback reduce star formation by only a modest amount (e.g. Rogers & Pittard, 2013; Dale, 2017; Howard et al., 2017), indicating that understanding the role of feedback in controlling clumps is still an open question.

The high resolution, the inclusion of metal-line cooling and the absence of any ab initio stars make our simulation susceptible to clump formation. Our results in this and the companion papers show that clumps may solve some important problems in galaxy formation that have not been explained fully satisfactorily otherwise. Nonetheless, the model in this paper is quite simple and there is much room for improvement. In particular we hope that future works will provide even more realistic models of clumps which will permit more detailed comparisons with data, rather than the proof-of-concept approach we have taken here.

V.P.D. is supported by STFC Consolidated grant # ST/R000786/1 and acknowledges the personal support of George Lake during part of this project, as well as the Pauli Center for Theoretical Studies, which is supported by the Swiss National Science Foundation (SNF), the University of Zürich, and ETH Zürich. S.R.L. acknowledges support from the Michigan Society of Fellows. S.R.L. was also supported by NASA through Hubble Fellowship grant HST-HF2-51395.001-A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. DBF acknowledges support from ARC Future Fellowship FT170100376. V.P.D. anknowledges the importance of ideas developed during the workshops “Disk instabilities across cosmic scales” (Sexten, July, 2017) and “Thin, thick and dark disks” (Ascona, July, 2017). V.P.D. thanks Lucio Mayer for his support in attending the former, and the Congressi Stefano Franscini for their financial support in organising the latter. The authors thank Giovanni Natale for his help with the use of DART-Ray. The simulation in this paper was run at the DiRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/J005673/1, STFC capital grant ST/H008586/1 and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure. Images with DART-Ray were built using the High Performance Computing Platform of Peking University. This work was completed at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This work was partially supported by a grant from the Simons Foundation. We thank the anonymous referee for their useful comments that helped improve this paper.