Metallicity Gradient of the Thick Disc Progenitor at High Redshift

TL;DR

This study uses a novel MCMC chemical painting technique to infer the initial metallicity gradients of the thick disc progenitor, revealing insights into its formation and evolution at high redshift.

Contribution

It introduces a new MCMC-based method to analyze metallicity gradients and applies it to APOGEE data, proposing a formation scenario for the thick disc.

Findings

Initial radial metallicity gradient was flat or positive.

Thick disc formed inside-out and upside-down.

Radial mixing influences metallicity gradient slopes.

Abstract

We have developed a novel Markov Chain Mote Carlo (MCMC) chemical "painting" technique to explore possible radial and vertical metallicity gradients for the thick disc progenitor. In our analysis we match an N-body simulation to the data from the Apache Point Observatory Galactic Evolution Experiment (APOGEE) survey. We assume that the thick disc has a constant scale-height and has completed its formation at an early epoch, after which time radial mixing of its stars has taken place. Under these assumptions, we find that the initial radial metallicity gradient of the thick disc progenitor should not be negative, but either flat or even positive, to explain the current negative vertical metallicity gradient of the thick disc. Our study suggests that the thick disc was built-up in an inside-out and upside-down fashion, and older, smaller and thicker populations are more metal poor. In this case, star forming discs at different epochs of the thick disc formation are allowed to have different radial metallicity gradients, including a negative one, which helps to explain a variety of slopes observed in high redshift disc galaxies. This scenario helps to explain the positive slope of the metallicity-rotation velocity relation observed for the Galactic thick disc. On the other hand, radial mixing flattens the slope of an existing gradient.

Tables

Table 1. Table 1: N-body model parameters.

Thick1					Thick2					Thin
$M_{\text{d,thick1}}$	$R_{\text{d,thick1}}$	$z_{\text{d,thick1}}$	N${}_{\text{p,thick1}}$		$M_{\text{d,thick2}}$	$R_{\text{d,thick2}}$	$z_{\text{d,thick2}}$	N${}_{\text{p,thick2}}$		$M_{\text{d,thin}}$	$R_{\text{d,thin}}$	$z_{\text{d,thin}}$	N${}_{\text{p,thin}}$
$M_{☉}$	kpc	kpc			$M_{☉}$	kpc	kpc			$M_{☉}$	kpc	kpc
$1.0\times {10}^{10}$	2.0	1.0	1,000,000		$5.0\times {10}^9$	3.0	0.5	500,000		$4.5\times {10}^{10}$	3.5	0.25	4,500,000

Table 2. Table 2: Results of the MCMC fitting

	Thick1+Thick2 or Thick1 only					Thick2
Case	[Fe/H]0	(d[Fe/H]/dR)0	(d[Fe/H]/dz)0	${\sigma }_{\left[\mathrm{Fe}/H\right],0}$		[Fe/H]0	(d[Fe/H]/dR)0	(d[Fe/H]/dz)0	${\sigma }_{\left[\mathrm{Fe}/H\right],0}$
	dex	dex kpc-1	dex kpc-1	dex		dex	dex kpc-1	dex kpc-1	dex
C1	$-0.289\pm 0.016$	$0.023\pm 0.002$	$-0.242\pm 0.008$	$0.192\pm 0.003$		$-$	$-$	$-$	$-$
C2	$-0.271\pm 0.022$	$0.012\pm 0.004$	$-0.206\pm 0.012$	$0.151\pm 0.006$		$0.451\pm 0.033$	$-0.048\pm 0.005$	$-0.726\pm 0.040$	$0.131\pm 0.014$

Table 3. Table 3: Radial and Vertical Metallicity Gradients.

Case	d[Fe/H]/dRa	d[Fe/H]/dzb
	dex kpc-1	dex kpc-1
APOGEE	$0.0005\pm 0.0014$	$-0.1760\pm 0.0068$
C1($t=8$)	$0.0144\pm 0.0001$	$-0.1173\pm 0.0015$
C2($t=8$)	$0.0020\pm 0.0001$	$-0.1477\pm 0.0015$
C2($t=0$)	$-0.0066\pm 0.0006$	$-0.3537\pm 0.0014$
C2 thick1($t=8$)	$0.0072\pm 0.0001$	$-0.0757\pm 0.0015$
C2 thick2($t=8$)	$-0.0204\pm 0.0002$	$-0.1717\pm 0.0037$

Equations

$${\rho }_{\text{dm}}=\frac{3H_0^2}{8\pi G}\frac{{\Omega }_0-{\Omega }_b}{{\Omega }_0}\frac{{\delta }_c}{cx(1+cx{)}^2},$$

$$c=\frac{r_{200}}{r_{\text{s}}},x=\frac{r}{r_{200}},$$

$$r_{200}=1.63\times 10^{-2}{\left(\frac{M_{200}}{h^{-1}{M}_{\odot }}\right)}^{\frac{1}{3}}{h}^{-1}\text{kpc},$$

$${\rho }_{\text{d}}=\frac{M_{\text{d}}}{4\pi z_{\text{d}}{R}_{\text{d}}^2}{\text{sech}}^2\left(\frac{z}{z_{\text{d}}}\right)\exp \left(-\frac{R}{R_{\text{d}}}\right),$$

$$\ln L$$

Full text

††pubyear: 2017††pagerange: Metallicity Gradient of the Thick Disc Progenitor at High Redshift–Metallicity Gradient of the Thick Disc Progenitor at High Redshift

1 Introduction

The Milky Way’s thick disc was originally identified from the analysis of star counts as a function of the vertical height (Yoshii, 1982; Gilmore & Reid, 1983). This traditional definition of the thick disc corresponds to a geometrically-defined thick disc (e.g. Robin & Creze, 1986; Yamagata & Yoshii, 1992; Reid & Majewski, 1993; Jurić et al., 2008). On the other hand, high-resolution spectroscopic studies of the solar neighbourhood stars show two sequences in the $\alpha$-element to the iron abundance ratio, [$\alpha$/Fe], as a function of the iron abundance, [Fe/H] (e.g. Fuhrmann, 1998; Prochaska et al., 2000; Feltzing et al., 2003). The higher [$\alpha$/Fe] sequences is usually considered as a chemically-defined thick disc. The chemically defined thick disc could be different from the geometrically-defined thick disc (Minchev et al., 2015) (see also Kawata & Chiappini, 2016, for a discussion about the various definitions of thick and thin discs). It is now confirmed that a similar high-[$\alpha$/Fe] sequence exists over a large radial and vertical range of the Galactic disc (e.g. Nidever et al., 2014; Hayden et al., 2015). From the Gaia-ESO survey, Mikolaitis et al. (2014) analysed the radial metallicity gradients of the chemically-defined thin and thick disc populations within the region at $6<R<12$ kpc and $|z|<0.6$ kpc. They found that the chemically-defined thick disc shows no radial metallicity gradient (see also Allende Prieto et al., 2006), but a shallow negative vertical metallicity gradient. This indicates that the chemical properties of the chemically-defined thick disc are well mixed radially, but maintain a negative vertical gradient (see also Bensby et al., 2011; Nidever et al., 2014; Hayden et al., 2015). Comparing the spectroscopically derived stellar parameters with theoretical isochrones, Haywood et al. (2013) analysed the stellar ages for chemically-defined thick and thin discs. In general, the chemically-defined thick disc stars are older than the chemically-defined thin disc stars (see also Bergemann et al., 2014; Masseron & Gilmore, 2015). Although some young high-[$\alpha$/Fe] disc stars are also found (Chiappini et al., 2015; Martig et al., 2015), some of them are likely to be blue stragglers (e.g. Jofré et al., 2016; Yong et al., 2016) and their origin is still debated. Therefore, in this paper, we consider that the chemically-defined thick disc stars are generally old, and define our thick disc as the chemically-defined old thick disc.

Advanced analysis of spectroscopic data starts to reveal the radial metallicity gradients of disc galaxies at $z\sim 1-3$, which is likely the formation epoch for the thick disc. Yuan et al. (2011) found a very steep radial metallicity gradient of $-0.16\pm 0.02$ dex kpc*-1* for a lensed disc galaxy at $z=1.49$. Jones et al. (2013) also reported two out of four lensed rotation-dominated galaxies observed at $z\sim 2$ showing radial metallicity gradients steeper than $-0.2$ dex kpc*-1*. If the thick disc of the Milky Way had this steep metallicity gradient, there must be some mechanism to flatten this gradient (without heating the thin disc component) in order to explain the current lack of a radial gradient in the Milky Way’s thick disc.

Interestingly, observations of more disc galaxies at high redshift show a variety of the gradients, including a significant fraction of galaxies with a flat metallicity gradient (Swinbank et al., 2012; Leethochawalit et al., 2016; Wuyts et al., 2016), and some positive ones (e.g. Cresci et al., 2010). Wuyts et al. (2016) analysed radial gradients of 180 KMOS3D near-infrared integral-field-unit (IFU) data at $z=0.6-2.7$, and found that most of their galaxies show flat radial metallicity gradient within uncertainties. Only 15 out of 180 galaxies show a radial metallicity gradient that is significantly negative (13 galaxies) or positive (2 galaxies).

The evolution of radial metallicity gradients should provide strong constraints on disc formation scenarios (e.g. Gibson et al., 2013), and therefore it is important to compare gradients of different age populations of the Milky Way with Milky Way-like disc galaxies at different redshifts. However, it is not straightforward to infer initial metallicity distribution for mono-age populations at the time they formed, from the current Galactic metallicity distribution as a function of age, because radial mixing can alter the distributions (e.g. Schönrich & Binney, 2009; Minchev et al., 2014; Grand et al., 2015; Kawata et al., 2017). Following Kawata et al. (2017), we refer to “radial mixing" to describe the overall radial re-distribution due to both “churning" and “blurring" (Schönrich & Binney, 2009; Schönrich & McMillan, 2017). Churning indicates the change of angular momentum of stars, and blurring describes the radial re-distribution of stars due to their epicyclic motion. Churning can be split into two mechanisms, “radial scattering" and “co-rotation radial migration". Radial scattering describes the change of angular momentum with kinematic heating, i.e. the increase in random energy of the orbit. This includes heating from mergers and satellite bombardment (e.g. Toth & Ostriker, 1992; Bird et al., 2012). On the other hand, co-rotation radial migration indicates a gain or loss of angular momentum of stars at the co-rotation resonance of the spiral arms (Sellwood & Binney, 2002) or bar, and does not involve significant kinematic heating.

As mentioned, the Milky Way’s thick disc has a flat radial metallicity gradient. The main question addressed in this paper is, what range of radial metallicity gradients is plausible for the progenitor of the Milky Way’s thick disc when it formed at high redshift. We use an idealised N-body simulation model for the thick and thin discs, where the thick disc formed at an early epoch, 8 Gyr ago, and experienced radial mixing due to bar and spiral arms developed in the thin disc. To limit the complexity of the simulation and highlight the effects of radial mixing, our numerical experiment does not include the growth of the thin disc. We also assume a constant scale-height for the thick disc, independent of radius, for simplicity. Within these assumptions, we introduce a novel Markov Chain Monte Carlo (MCMC) painting technique which assigns metallicity tags to N-body particles by exploring the possible radial and vertical metallicity gradients before they suffered from radial mixing, attempting to fit the current metallicity distribution of the Milky Way’s thick disc observed by the Apache Point Observatory Galactic Evolution Experiment (APOGEE) at different radii and vertical heights (Hayden et al., 2015).

The formation of the thick disc is expected to have had been completed at an early epoch, before the thin disc starts growing, as suggested by cosmological simulations (e.g. Brook et al., 2004, 2006; Brook et al., 2012; Bird et al., 2013; Stinson et al., 2013; Minchev et al., 2013; Minchev et al., 2015; Grand et al., 2016). These simulations also suggest that the thick disc formed in an inside-out and upside-down fashion, i.e. the early thick disc was smaller and thicker (e.g. Brook et al., 2006; Brook et al., 2012; Bird et al., 2013; Minchev et al., 2015). Using two components for the thick discs, a thicker, smaller one and a thinner, larger one, we also qualitatively study how the inside-out and upside-down formation of the thick disc affects the possible metallicity distribution of the thick disc progenitor. Section 2 describes briefly the numerical simulation code and numerical models. In Section 3, we present our results. A summary and discussion of this study is presented in Section 4.

2 Method and Models

We use our original Tree N-body code, GCD+ (Kawata & Gibson, 2003; Kawata et al., 2013) for the numerical experiments presented in this paper, and model the evolution of a barred disc galaxy similar in size to the Milky Way. We initially set up an isolated disc galaxy which consists of stellar discs, with no bulge component, in a static dark matter halo potential (Rahimi & Kawata, 2012; Grand et al., 2012; Kawata et al., 2017). A live dark matter halo can respond to the disc particles by exchanging angular momentum. However, if not properly modelled, a live dark matter halo may introduce some numerical scattering and heating (e.g. D’Onghia et al., 2013). In the interest of performing a more controlled experiment and reducing the computational burden, we use a static dark matter halo. We use the standard Navarro-Frenk-White (NFW) dark matter halo density profile (Navarro et al., 1997), assuming a $\Lambda$-dominated cold dark matter ($\Lambda$CDM) cosmological model with cosmological parameters of $\Omega_{0}=0.266=1-\Omega_{\Lambda}$, $\Omega_{\rm b}=0.044$, and $H_{0}=71{\rm kms^{-1}Mpc^{-1}}$:

[TABLE]

where

[TABLE]

and

[TABLE]

where $\delta_{c}$ is the characteristic density of the profile (Navarro et al., 1997), $r$ is the distance from the centre of the halo and $r_{s}$ is the scale radius. The total halo mass is set to be $M_{200}=1.2\times 10^{12}M_{\odot}$ and the concentration parameter is set at $c=14$. The halo mass is roughly consistent with recent estimates of the mass of the Milky Way (e.g. McMillan, 2011; Bland-Hawthorn & Gerhard, 2016).

The stellar disc is assumed to follow an exponential surface radial density profile and $\textup{sech}^{2}$ vertical density profile:

[TABLE]

where $M_{\text{d}}$ is the stellar disc mass, $R_{\text{d}}$ is the scale length and $z_{\text{d}}$ is the scale height. In this paper, we consider two thick disc components and one thin disk component. The model parameters are summarised in Table 1. The first thick disc component (thick1) is smaller and thicker than the 2nd thick disc component (thick2). The rationale behind these two thick disc components is explained later. In Table 1, N${}_{\text{p}}$ is the number of particles for each component. This model has 10,000 M*⊙* for each particle. We apply the spline softening suggested by Price & Monaghan (2007). We set a softening length of 341 pc (the equivalent Plummer softening length is about 113 pc). To relax the initial system, we ran a simulation for 2 Gyr, moving randomly selected particles azimuthally, which suppresses the development of non-axisymmetric structures (McMillan & Dehnen, 2007). We then run the simulations for 8 Gyr which roughly corresponds to the age of the thin disc of the Milky Way. Note that the time of $t=0$ in this paper indicates the starting time of this 8 Gyr simulation using the the end-product of the relaxing run as initial conditions. The N-body simulation is similar to the simulation in Kawata et al. (2017) where we discuss that numerical heating is negligible in this simulation setting.

The simulation of the 8 Gyr evolution of the system mimics the evolution history of the Galactic disc after the thick disc formed. For simplicity, we set a thin disc already fully grown from the start, ignoring its growth. This is too simplistic, but we are interested in the effect of radial mixing driven by the thin disc on the thick disc stars. To focus on radial mixing driven by non-axisymmetric structures in the thin disc, we do not include mergers or tidal interaction with the sub-halos or satellite galaxies. This assumption could be justified for the Milky Way, since there is ample observational evidences that the Milky Way has had a quiescent accretion history (e.g. Ruchti et al., 2014; Deason et al., 2017). If there was a merger which was violent enough to induce a strong phase mixing of disc stars after they formed, we expect that it would be difficult to maintain the observed negative metallicity gradient of the thick disc and mono-age populations of the Galactic disc (see also Xiang et al., 2015; Ciucă et al., 2017), because strong mixing would erase the gradient.

Two component thick discs allow us to explore the inside-out and upside-down formation scenario of the thick disc (e.g. Brook et al., 2004; Bird et al., 2013). We set the thick1 disc being heavier than the thick2 disc. This mimics a scenario that the thick1 disc has grown before the thick2 disc formed, and the thick2 disc represents the latest generation of the thick disc. We study what radial and vertical metallicity gradients these two thick disc components should "originally" have had by comparing with the APOGEE observations of metallicity distribution function (MDF) at different Galactocentric radii and height (Hayden et al., 2015).

3 Results

Fig. 1 shows the evolution of the simulated galaxy. The system fully develops a central bar by $t=3.0$ Gyr, because of the dominated stellar mass compared to the dark matter mass in the central region. Fig. 1 shows that both thick disc components also develop a bar, although it is rounder than the thin disc one. This is consistent with what is shown in Bekki & Tsujimoto (2011); Kawata et al. (2017). Fig. 2 presents the angular momentum change between $t=0$ and $t=8$ Gyr, i.e. churning which includes both radial scattering and co-rotation radial migration. As expected, churning is relatively stronger for the thin disc than the thick discs. However, there are significant churning in both the thick1 and thick2 discs. Because the thick2 disc is kinematically colder than the thick1 disc, more significant churning can be seen in the thick2 disc than in the thick1 disc (cf. Solway et al., 2012).

Fig. 3 shows the radial profile of the cumulative specific angular momentum, $L/M$, surface density, $\Sigma$, root mean square (RMS) of the height, $<\,z^{2}>^{1/2}$, and velocity dispersions of the radial, $\sigma_{R}$, azimuthal, $\sigma_{\phi}$, and vertical, $\sigma_{z}$, directions for the thin (blue), thick1 (black) and thick2 (red) disc components at $t=0$ (dotted lines) and 8 (solid lines) Gyr. Due to bar formation, all disc particles in the radial range plotted in Fig. 3 lost angular momentum so the surface density is more centrally concentrated at $t=8$ Gyr. There is a slight increase in surface density in the outer disk $R>8$ kpc especially for the thick1 disc. This is because thick1 disc was initially the most compact component, and more stars gained angular momentum compared to the other components as seen in the top panel. Hence, radial mixing driven by a larger disc component can result in radial expansion of the more compact disc component, though the effect is not significant (see also Aumer et al., 2016).

For all the components, the RMS height is smaller in the central regions, as a result of the deeper gravitational potential. In the outer region, $R>\sim 3$ kpc, the RMS height becomes larger compared to the initial state. This is mainly due to the combination of heating, which can be seen as an increase in $\sigma_{z}$ in the bottom panel, and enhanced vertical oscillations of the stars that migrated outward. The latter mechanism is caused by higher initial vertical action for the outward migrators than the non-migrated ones, because of the deeper gravitational potential in the inner region (e.g. Minchev et al., 2012). Interestingly, thick1 disc shows an obvious increase in $<\,z^{2}>^{1/2}$ with radius, i.e. flaring. The radial dependence of the scale-height for the thick disc is still under debate. Using APOGEE data, Bovy et al. (2016) showed that the chemically-defined thick disc has a constant scale-height, independent of radii, as inferred from the analysis of mono-abundance populations of the Galactic disc. On the other hand, Minchev et al. (2017) argued that the mono-abundance populations are not same as the mono-age populations, and the old chemically-defined thick disc population can have a flared disc. This idea is supported by Mackereth et al. (2017), who analysed mono-age components of mono-abundance populations in the APOGEE data, and found that the old mono-age thick disc component shows flaring (see also Robin et al., 2014). Although the structure of our disc components is beyond the scope of this study, it is interesting to note that, in our simulations, radial mixing leads to a flaring thick disc.

The thin disc of our simulation is too hot compared to the velocity dispersion of stars in the solar neighbourhood. Our thin disc stars in $7.5<R<8.5$ kpc and $-0.5<z<0.5$ kpc have velocity dispersions of $(\sigma_{R},\sigma_{\phi},\sigma_{z})=(48.2,36.5,16.3)$ km s*-1*, while the solar neighbourhood thin disc stars defined as [$\alpha$/Fe]$<0.1$ measured in Gaia DR1 combined with APOGEE data shows $(37\pm 2,23\pm 1,18\pm 1)$ km s*-1* (Allende Prieto et al., 2016). In this idealised numerical experiment, we do not include newly formed stars, and those would add kinematically colder stars to the thin disc component (e.g. Sellwood & Carlberg, 1984; Aumer et al., 2016). This simplification leads to a velocity dispersion of the thin disc that is too high at $t=8$ Gyr. In this paper, we focus on the thick disc, and the thin disc component is merely used as a driver of radial mixing. We do not analyse the properties of the thin disc component, and therefore, this discrepancy is unlikely to be an issue for our qualitative study of the effect of radial mixing on the metallicity distribution of the thick disc.

The relatively high velocity dispersion mentioned above implies that our N-body simulation likely overestimates radial mixing by radial scattering mainly due to the bar formation, compared to what is expected to happen in the Milky Way. On the other hand, a high velocity dispersion weakens the spiral arms and reduces co-rotation radial migration, which underestimates radial mixing. Because we do not know what is the degree of radial mixing in the Milky Way, it is difficult to assess whether or not our simulation has an appropriate level of radial mixing for the thick disc compared to the Milky Way. However, as seen in Fig. 2, both thick disc components experience significant radial mixing. Hence, our numerical experiment is suitable to provide a qualitative evaluation of how radial mixing affects the radial and vertical metallicity gradient, but we will avoid a quantitative discussion.

Fig. 4 shows the [$\alpha$/Fe] and [Fe/H] distribution of APOGEE stars in Sloan Digital Sky Survey (SDSS) data release 12 from Hayden et al. (2015). Although the data show two distinct sequences for the [$\alpha$/Fe] and [Fe/H] relation, there is no obvious definition for chemically identified thick and thin disc stars. In this paper, we define thick disc stars to be ${\rm[\alpha/Fe]}>0.17$ at [Fe/H]$<-0.6$ and ${\rm[\alpha/Fe]}>-0.0875(\text{[Fe/H]}-0.2)+0.1$ at [Fe/H]$>-0.6$, which shows a reasonable division of the sequences in Fig. 4. We consider that these chemically-defined thick disc stars represent the old thick disc population, and we refer to them as the APOGEE thick disc stars. We explore what radial and vertical metallicity gradients are allowed in the thick disc when is born, to explain the MDF of the APOGEE thick disc stars at different radii and vertical heights.

3.1 Single metallicity distribution case

First, we assume both the thick1 and thick2 discs to have the same initial metallicity distribution. We call this Case C1. Grey lines and error bars in Fig. 5 represent the MDF for the APOGEE thick disc stars at different radial and height regions. We divide the sample of stars within $3<R<13$ kpc with a 2 kpc bin size, and 3 different height bins of $0<|z|<0.5$, $0.5<|z|<1$ and $1<|z|<2$ kpc. The error bars correspond to Poissonian statistical errors. We run MCMC with 4 parameters, the initial $R=0$ and $|z|=0$ metallicity, [Fe/H]0, radial metallicity gradient, (d[Fe/H]/dR)0, vertical metallicity gradient, (d[Fe/H]/dz)0, and metallicity dispersion, $\sigma_{\rm[Fe/H],0}$, at $t=0$, to explore what parameter set would best match with the APOGEE thick disc MDFs. For each MCMC chain, we assign [Fe/H] to particles using the position at $t=0$ to follow the parameters of [Fe/H]0, (d[Fe/H]/dR)0, (d[Fe/H]/dz)0 and $\sigma_{\rm[Fe/H],0}$. We also add Gaussian random errors to [Fe/H] of $0.05$ dex to take into account the typical observational errors for APOGEE data (Hayden et al., 2015). Then, using their position at $t=8$ Gyr, we analyse the MDF in the same volumes as for the APOGEE thick disc data in Fig. 5, and calculate the following log likelihood function.

[TABLE]

where $dN({\rm[Fe/H]_{\it i}})$ are the normalised number of APOGEE thick disc stars or model particles which are in the $i$-th [Fe/H] bin, and $i$ indicates each [Fe/H]i bin in the MDFs for all different volumes described above, except at the location where $dN({\rm[Fe/H]}_{i,{\rm obs}})=0$. We define the observational errors, $\sigma_{\rm[Fe/H]_{{\it i},obs}}$, of the $dN({\rm[Fe/H]_{\it i}})$ as Poisson errors from the number of stars in each bin. We use flat priors in the ranges of $-2.0<{\rm[Fe/H]}_{0}<2.0$ dex, $-0.5<{\rm(d[Fe/H]/dR)_{0}}<0.5$ dex kpc*-1*, $-1.0<{\rm(d[Fe/H]/dz)_{0}}<1.0$ dex kpc*-1* and $0.0<\sigma_{\rm[Fe/H],0}<1.0$ dex. We use emcee (Foreman-Mackey et al., 2013) for the MCMC sampler, and use 128 walkers and 2,000 chains. We call this assigning [Fe/H] to particles with "MCMC chemical painting" or "MCMC painting". Note that chemically painting N-body particles has been used in Minchev et al. (2013, 2014); Kubryk et al. (2013, 2015) where a chemical evolution model is used to assign metallicities to the particles. Our method is simpler, but allows us to fit the observational data with MCMC in a reasonable computational time. Hence, these different methods are complementary to each other.

The best fit parameters found with the MCMC chemical painting analysis are shown in Table 2 and the MDFs of the best model are shown in Fig. 5 (black lines). Our MCMC painting analysis for Case C1 suggests that a positive radial metallicity gradient of (d[Fe/H]/dR)${}_{0}=0.023\pm 0.002$ dex kpc*-1* for the initial disc is preferred to reproduce the observed MDF. A positive metallicity gradient may come as a surprise, because this is rare in disc galaxies. However, we find that this is required to reproduce the flat radial metallicity gradient and the negative vertical metallicity gradient observed in the thick disc of the Milky Way at present.

We measured the radial metallicity gradient of the APOGEE thick disc stars in the region $3<R<13$ kpc and the vertical height within $|z|<2$ kpc, and found d[Fe/H]/dR$=0.0005\pm 0.0014$ dex kpc*-1* as shown in Table 3. We also measured the vertical metallicity gradient in the region of $7<R<9$ kpc, and found d[Fe/H]/dz$=-0.1760\pm 0.0068$ dex kpc*-1* at $|z|<2$ kpc. As mentioned in Section 1, the thick disc stars have a flat radial metallicity gradient and a clear negative vertical metallicity gradient. We also measured the radial metallicity gradients in the same regions with the same way for our best fit model at $t=8$ Gyr (Tab. 3). We found d[Fe/H]/dR$=0.0144\pm 0.0001$ dex kpc*-1*, which is reasonably similar to the flat radial metallicity gradient found from the APOGEE thick disc stars111 Interestingly, there is a hint of positive radial metallicity gradient of the old thick disc discussed in Fig. 18 of Casagrande et al. (2011).. Hence, radial mixing induced by the bar and spiral arms seems to have effectively flattened the radial metallicity gradient from the initially assumed positive metallicity gradient. An important mechanism is that the stars brought from the inner (outer) region due to the radial mixing tend to populate the high (low) vertical height region, as shown in Kawata et al. (2017). Hence, if there is initially a negative radial metallicity gradient, radial mixing will make the vertical metallicity gradient positive, which is inconsistent with the clear negative vertical metallicity gradient seen in the APOGEE data. This is why our MCMC result arrived at a positive radial metallicity gradient for the initial disc; to reproduce the negative vertical metallicity gradient after radial mixing brought stars from the inner disc and populated with them the regions far from the plane. This result suggests that the negative vertical metallicity gradient observed in the thick disc of the Milky Way provides a strong constraint on the overall radial metallicity gradient of the thick disc progenitor and calls for a positive radial metallicity gradient.

Fig. 5 shows the best model MDF misses higher metallicity stars at $|z|<0.5$ kpc, which explains the shallower vertical metallicity gradient of the best model of Case C1 compared to the APOGEE thick disc data, as seen in Table 3. Still, this is the best model found by our MCMC run with the one component model. This model is informative to understand the general requirement to explain the observed MDF of the APOGEE thick disc data. However, this inconsistency indicates that we need more components to explain the observed MDF. Hence, in the next section, we consider a two component model.

3.2 Double metallicity distribution case

We now consider the thick1 and thick2 discs separately, i.e. we allow different radial and vertical metallicity gradients to be assigned for these two discs. We run MCMC with 8 parameters, i.e. two different sets of [Fe/H]0, (d[Fe/H]/dR)0, (d[Fe/H]/dz)0 and $\sigma_{\rm[Fe/H],0}$ for thick1 and thick2 discs at $t=0$, and analyse the log likelihood in the same way as equation (5) in previous section. We set the same flat priors as mentioned in Section 3.1. We call this Case C2.

The best fit parameters for Case C2 are summarised in Table 2, and the MDFs of the best fit model are shown in Fig. 6. Fig. 6 shows that the two component disc case reproduces better the APOGEE thick disc data than Case C1. The higher metallicity tails of the MDFs at $|z|<0.5$ kpc and $R<9$ kpc are described better. Overestimates of high metallicity stars at $1<|z|<2$ kpc and $R>11$ kpc in Case C1 are mitigated in Case C2.

There are some regions where our best fit model of Case C2 still struggles to reproduce the observed MDF. However, it is not our aim to reproduce the observed MDF perfectly. We are using only one N-body simulation result with one evolution history. The aim of this paper is to qualitatively discuss the effect of radial mixing on the thick disc metallicity distribution, and not to find the best evolution history of the Galactic disc. Therefore, we consider that we have a reasonable success in qualitatively reproducing the APOGEE data, and discuss what we can learn from the fact that two component case is better than one component case.

Table 2 shows that to reproduce the APOGEE thick disc MDF, the thick1 disc should have a positive radial metallicity gradient of (d[Fe/H]/dR)${}_{0}=0.012\pm 0.004$ dex kpc*-1*, like Case C1, while the thick2 disc should have a negative radial metallicity gradient of (d[Fe/H]/dR)${}_{0}=-0.048\pm 0.005$ dex kpc*-1*. Also, the slope required for the thick1 disc is shallower than the slope required for Case C1. Note that the thick1 disc is more massive than the thick2 disc, and the thick1 disc is thicker in scale height and smaller in scale radius. Our best fit MCMC result shows that the thick2 disc should be more metal rich than the thick1 disc. Hence, if we combine the thick1 and thick2 discs, the inner, higher regions are dominated by the metal poor thick1 disc, and the outer region and lower height region is more dominated by the metal rich thick2 disc. As a result, this combined disc provides a similar metallicity distribution to a single positive radial metallicity gradient in Case C1. We measured the overall radial metallicity gradient at $t=0$ for the best model of Case C2 in the same way as Section 3.1, and found d[Fe/H]/dR$=-0.0066\pm 0.0006$, which is almost flat or slightly negative (labeled "C2 (t=0)" in Table 3). However, note that we measured the metallicity gradients in Table 3 by a linear regression of the mean [Fe/H] of all the particles of both thick1 and thick2 discs at radial and vertical bins. In other words, this does not take into account the density variations with radii. Hence, this result likely underestimates the contribution of the metal poor particles in the inner region from the thick1 disc. Fig. 7 shows the radial metallicity distribution for the whole disc and thick1 and thick2 disc components at $t=0$ and 8 Gyr. In fact, metal poor particles are dominant in the inner region, but they are not contributing to our linear regression, because our linear regression is done outside of this region, as seen in the top right panel of Fig. 7. We therefore consider that, overall, the best model in Case C2 has more metal poor stars in the inner region and more metal rich stars in the outer region. As discussed in Section 3.1, this is required to reproduce the negative vertical metallicity gradient in the APOGEE data.

We conclude that the thick disc of the Milky Way can be explained with two or multiple generations of the discs. Considering the age-metallicity relation (older thick disc stars having lower [Fe/H]) observed in the thick disc (Haywood et al., 2013; Bensby et al., 2014; Bergemann et al., 2014), we can suggest a scenario where a thicker and smaller thick1-like disc formed first, and then later a more metal rich thinner and larger thick2-like disc was built up. Alternatively and more naturally, the thick disc was built up by gradually increasing in radial size and metallicity, and decreasing in thickness, i.e. an inside-out and upside-down fashion (Brook et al., 2006; Bird et al., 2012). In the latter scenario, the earlier generations of the thick disc are more metal poor and more compact, which leads to an overall positive radial metallicity gradient, despite individual generations of the disc can have had a flatter or even negative radial metallicity gradient. This is first shown in Fig. 10 of Minchev et al. (2014) and also discussed in Schönrich & McMillan (2017). These results mean that the inside-out formation can allow the thick disc star-forming region to have a variety of radial metallicity gradients. This helps to explain the range of radial metallicity gradients observed in high redshift disc galaxies, as mentioned in Section 1. On the other hand, radial mixing limits the range of radial metallicity gradient possible for the thick disc progenitor, since the current thick disc has a clear negative vertical metallicity gradient.

It is worth noting that the preferred solution for the thick2 disc is an initial steep vertical metallicity gradient of d[Fe/H]/dz$=-0.726\pm 0.040$. This is required to reproduce the negative vertical metallicity gradient in APOGEE data, because the negative radial metallicity gradient and radial mixing make the vertical metallicity gradient more positive. It is difficult to imagine the star forming gas disc having such a steep vertical metallicity gradient, especially if turbulent mixing is strong in high redshift galaxies (Yang & Krumholz, 2012; Grand et al., 2014; Petit et al., 2015). If we assume that metal mixing is effective and there should be no vertical metallicity gradient in the star forming disc at a fixed time, the steep vertical metallicity gradient indicates that the coeval population of the thick disc progenitor forming at different epochs gradually became thinner with time, and later populations are more metal rich. Then, the overall negative vertical metallicity gradient can be established by the thick disc built up in an upside-down fashion, with younger populations having higher [Fe/H]. Hence, to describe the thick disc formation, we need to consider multiple-generations of discs with different radial sizes and thicknesses and different metallicities or smooth growth of the disc, rather than only a two component disc.

3.3 [Fe/H] and the Rotation Velocity Relationship

The thick disc stars defined by their [$\alpha$/Fe] abundance have a distinctive feature in the [Fe/H]-$V_{\phi}$ relation, compared to the chemically-defined thin disc in the solar neighbourhood stars. Observational studies demonstrate that thick disc stars have a positive slope of d$V_{\phi}$/d[Fe/H], while the thin disc stars show a negative slope (e.g. Spagna et al., 2010; Lee et al., 2011; Adibekyan et al., 2013; Recio-Blanco et al., 2014; Allende Prieto et al., 2016). There are several theoretical possibilities to explain this trend (e.g. Loebman et al., 2011; Curir et al., 2012; Allende Prieto et al., 2016; Schönrich & McMillan, 2017). In this section, we discuss the [Fe/H]-$V_{\phi}$ relation for the thick disc only, given it is the focus of this paper. Note that the negative gradient d$V_{\phi}$/d[Fe/H] for the thin disc can be easily explained by the observed radial metallicity gradient and stellar epicycle motion, i.e. blurring (Vera-Ciro et al., 2014; Allende Prieto et al., 2016).

Fig. 8 shows the rotation velocity, $V_{\phi}$, as a function of [Fe/H] for the thick disc particles for the best model for Cases C1 and C2 and thick1 and thick2 discs in Case C2. As discussed in Curir et al. (2012); Allende Prieto et al. (2016); Schönrich & McMillan (2017), a positive slope of d$V_{\phi}$/d[Fe/H] can be explained with an overall positive radial metallicity gradient, which explains the clear positive slope in Case C1 and thick1 disc in Case C2. In these cases, the stars whose guiding radius is smaller (larger) are more metal poor (rich), and tend to be rotating slower (faster), because they tend to be near apocentre (pericentre) phase. This is due to the blurring effect, not due to churning, i.e. the change in angular momentum. The steeper slope in the inner region is due to the increase in the velocity dispersion in the inner disc (see Fig. 3), which leads to a stronger change in $V_{\phi}$ at a given radius.

It is interesting that in Case C1 and thick1 disc in Case C2 the slope of d$V_{\phi}$/d[Fe/H] is minimum at $R=8$ kpc, and becomes steeper at $R>8$ kpc. We found that this is because around $R=8$ kpc churning is most effective. The top panels of Fig. 9 show the initial radii, $R_{\rm ini}$, of particles at $R=4$, 8 and 12 kpc at $t=8$ Gyr. In the inner (outer) region, i.e. at $R=4$ (12) kpc, more particles came from the outer (inner) disc, i.e. $R_{\rm ini}$ is larger (smaller) than the current radius, because of loss (gain) of angular momentum. On the other hand, at $R=8$ kpc the particles came from both inner and outer region. We calculated the guiding radius, $R_{\rm g}$, at $t=8$ Gyr, by comparing angular momentum of the particles and the circular velocity at the disc plane, which are shown in the middle panels of Fig. 9. It is interesting to see that the guiding radius tends to be smaller (higher) than the current radius for the particles which lost (gained) their angular momentum. For example, for the particles at $R=8$ kpc, those with negative (positive) $\Delta L$ tend to have smaller (larger) $R_{\rm g}$ than the red line of $R=8$ kpc in the middle panel of Fig. 9. As a result, the particles that lost (gained) angular momentum tend to be detected at apocentre (pericentre) phase, and therefore are rotating more slowly (faster). The particles from larger (smaller) $R_{\rm ini}$ are more metal rich (poor), because of the initial positive radial metallicity gradient in Case C1. Consequently, the particles which changed their guiding centre from the outer (inner) region to the inner (outer) region contributes more to the metal rich (poor) and slowly (fast) rotating region in the [Fe/H]-$V_{\phi}$ relation, flattening a positive slope of d$V_{\phi}/$d[Fe/H]. This result implies that churning can flatten a positive d$V_{\phi}$/d[Fe/H] slope. We also confirmed that if we picked up particles with small $|\Delta L|$, the d$V_{\phi}$/d[Fe/H] slope is more positive in Case C1. Because at $R=8$ kpc there are particles that came from both sides of the disc, d$V_{\phi}$/d[Fe/H] becomes more flatten by churning. Some indication of this kind of trend is observed in Kordopatis et al. (2017).

It is interesting that the thick2 disc shows negative d$V_{\phi}$/d[Fe/H] slopes at all radii. It is basically the same mechanism as the thin disc case as discussed above, i.e. blurring and negative radial metallicity gradient drives a negative d$V_{\phi}$/d[Fe/H] slope. Because the outer region is dominated by the thick2 disc in Case C2, the slope of d$V_{\phi}/$d[Fe/H] becomes shallower, or even negative, in the outer disc (2nd top panels in Fig. 8). We expect that if the thick disc was built up in an inside-out fashion and the star forming disc at different epochs had a negative radial metallicity gradient, an mono-age thick disc population at $z=0$ could have a negative d$V_{\phi}$/d[Fe/H] like the thick2 disc. Then, the overall d$V_{\phi}/$d[Fe/H] as a function of radius for the thick disc depend on the inside-out formation history and the radial metallicity gradient for each mono-age population.

At $R=8$ kpc, our models show too shallow slope of d$V_{\phi}$/d[Fe/H], compared to the observed one for the Galactic thick disc, e.g. one of the shallowest slope of d$V_{\phi}$/d[Fe/H]$=23\pm 10$ km s*-1* is reported in (Allende Prieto et al., 2016). If the thick disc consists of many generations of discs with different thickness, scale-length and [Fe/H], we can change this slope by adjusting it to reproduce the observed slopes, as Schönrich & McMillan (2017) demonstrated. However, such work is beyond the scope of this paper. Our qualitative study demonstrates that an overall positive d$V_{\phi}$/d[Fe/H] slope can be driven by an overall positive radial metallicity gradient, and radial mixing of churning does flatten this slope. Hence, the measurement of d$V_{\phi}$/d[Fe/H] as a function of radius for different mono-age populations provides strong constraints on the formation history of the thick disc.

4 Summary and Discussion

Combining an idealised N-body model of a disc galaxy similar in size to the Milky Way with a new MCMC chemical painting technique, we explore what metallicity distribution for the thick disc progenitor is needed, to explain the current thick disc metallicity distribution observed by the APOGEE survey. The main assumptions in our study are that thick disc formation was completed a long time (8 Gyr) ago and the thick disc stars experienced effective radial mixing for a long time since then. In addition, we assumed that the vertical scale height of the overall thick disc component is constant at different radii. Our key findings are the followings.

(i) The current flat radial metallicity gradient and negative vertical metallicity gradient of the thick disc stars observed in the APOGEE data provide strong constraints on the radial metallicity gradient for the thick disc progenitor; it should be flat or even positive. A negative radial metallicity gradient is not likely. This is because radial mixing with a negative radial metallicity gradient leads to a positive vertical metallicity gradient as demonstrated in Kawata et al. (2017), which is in conflict with observations.

(ii) Two or more components for the thick disc can explain better the MDF of the thick disc stars. Our result supports a scenario in which the thick disc formed in an inside-out and upside-down fashion, and an older, thicker and smaller disc is more metal poor. In this case, each generation of thick disc could have an initial negative radial metallicity gradient. Then, inside-out formation and age-metallicity relation produced an overall flat or positive metallicity gradient when their formation was completed.

(iii) Hence, the inside-out formation can allow the thick disc star-forming region to have a variety of (positive, flat or negative) radial metallicity gradients. This helps to explain the variety of radial metallicity gradients observed at high redshift disc galaxies. On the other hand, radial mixing rather limits the range of allowed radial metallicity gradient as discussed in (i).

(iv) The overall positive radial metallicity gradient, built-up by inside-out formation and the age-metallicity relation, helps to explain a positive d$V_{\phi}$/d[Fe/H] slope observed for the thick disc stars, as suggested by Schönrich & McMillan (2017). In this case, radial mixing flattens the positive d$V_{\phi}$/d[Fe/H] slope, especially around the solar radius.

To our knowledge, (i) and (iii) are new finding, although (iii) can be inferred from the results of Schönrich & McMillan (2017). (ii) is consistent with the scenario discussed in the semi-analytic chemo-dynamical model of Schönrich & McMillan (2017) and in the cosmologically motivated chemo-dynamical model of Minchev et al. (2014). The first part of (iv) is an independent confirmation of the scenario suggested by Schönrich & McMillan (2017), but the radial dependence of the d$V_{\phi}$/d[Fe/H] slope is discussed here for the first time. These results indicate that if the thick disc formed in an inside-out and upside-down fashion with a clear age-metallicity relation, different generations of the thick disc progenitor should have associated different chemo-dynamical structures, i.e. different sizes and radial metallicity gradients and different d$V_{\phi}$/d[Fe/H] slopes at different radii. If we could measure the age of the old stars accurately enough, we would expect that the older thick disc would be smaller and more metal poor, and each generation of the thick disc would exhibit a more negative radial metallicity gradient than the overall gradient of the thick disc. Also, d$V_{\phi}$/d[Fe/H] for a mono-age thick disc population would be flatter or even negative (if their radial metallicity gradient is also negative).

Unfortunately, the thick disc formation time-scale is relatively short (a few Gyr, e.g. Haywood et al., 2013). It may be challenging to divide the sample of thick disc stars into age bins fine enough to test this scenario. Hence, this would be a key challenge for future work on Galactic archaeology. The European Space Agency (ESA)’s Gaia mission will soon provide us accurate distance and photometric information supplemented with metallicity measurements from high-resolution spectroscopic surveys. The data will bring age estimates for turn-off stars from well-calibrated isochrones. In addition, for both dwarf and giants, asteroseismology will provide age estimates from the high-quality light-curve data from the on-going K2 mission with Kepler (e.g. Stello et al., 2015) and up-coming TESS and Plato missions (e.g. Miglio et al., 2017). Our study demonstrates that fitting such observational data with N-body simulations with MCMC painting would be a valuable tool to understand the building history of the Galactic disc. In this paper, we only used one simplified simulation model without any growth of the disc mass. Encouraged by the success of this study, we will further develop N-body disc models including the growth of the disc, as in Aumer et al. (2016), and run many different models with various formation scenarios. Combining N-body simulation models and MCMC painting against future observational data must help us to decipher the formation and evolution of the Milky Way disc.

Acknowledgments

We thank an anonymous referee for their constructive comments and helpful suggestions which have improved the manuscript. DK also thanks Ivan Minchev and Ralph Schönrich for providing insightful comments. DK and IC acknowledge the support of the UK’s Science & Technology Facilities Council (STFC Grants ST/K000977/1 and ST/N000811/1). CAP is thankful to the Spanish MINECO for funding through grant AYA2014-56359-P. LC gratefully acknowledges support from the Australian Research Council (grants DP150100250, FT160100402). RJJG acknowledges support by the DFG Research Centre SFB-881 ’The Milky Way System’, through project A1. JH is supported by a Dunlap Fellowship at the Dunlap Institute for Astronomy & Astrophysics, funded through an endowment established by the Dunlap family and the University of Toronto. This work used the UCL facility Grace and the DiRAC Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.