The structural evolution of galaxies with both thin and thick discs

TL;DR

This study uses N-body simulations to explore how the growth of thin and thick galactic discs influences galaxy structure, bar formation, and age distribution, providing insights into Milky Way-like galaxy evolution.

Contribution

It introduces controlled simulation models of disc galaxy evolution with thick and thin components, revealing their impact on galaxy morphology and dynamics.

Findings

Double-exponential vertical profiles similar to the Milky Way.

Thick discs delay and influence bar formation.

Vertical age distribution shaped by bar buckling.

Abstract

We perform controlled N-body simulations of disc galaxies growing within live dark matter (DM) haloes to present-day galaxies that contain both thin and thick discs. We consider two types of models: a) thick-disc initial conditions to which stars on near-circular orbits are continuously added over ~10 Gyr and b) models in which the birth velocity dispersion of stars decreases continuously over the same time-scale. We show that both schemes produce double-exponential vertical profiles similar to that of the Milky Way (MW). We indicate how the spatial age structure of galaxies can be used to discriminate between scenarios. We show that the presence of a thick disc significantly alters and delays bar formation and thus makes possible models with a realistic bar and a high baryon-to-DM mass ratio in the central regions, as required by microlensing constraints. We examine how the radial mass distribution in stars and DM is affected by disc growth and non-axisymmetries. We discuss how bar buckling shapes the vertical age distribution of thin- and thick-disc stars in the bar region. The extent to which the combination of observationally motivated inside-out growth histories and cosmologically motivated dark halo properties leads to the spontaneous formation of non-axisymmetries that steer the models towards present-day MW-like galaxies is noteworthy.

Tables

Table 1. Table 1: An Overview over the different ICs and their parameters: 1st Column : IC name; 2nd Column : IC disc mass M disc , i subscript 𝑀 disc i M_{\rm disc,i} ; 3rd Column : IC bulge mass M bulge , i subscript 𝑀 bulge i M_{\rm bulge,i} ; 4th Column : number of baryonic particles N b , i subscript 𝑁 b i N_{\rm b,i} in the ICs; 5th Column : concentration parameter for IC DM halo, c halo subscript 𝑐 halo c_{\rm halo} ; 6th Column : IC DM halo scalelength a halo subscript 𝑎 halo a_{\rm halo} ; 7th Column : IC radial disc scalelength h R , disc subscript ℎ 𝑅 disc h_{R,{\rm disc}} ; 8th Column : IC vertical disc scaleheight z 0 , disc subscript 𝑧 0 disc z_{0,{\rm disc}} ; 9th Column : IC bulge scalelength a bulge subscript 𝑎 bulge a_{\rm bulge} ; 10th Column : IC bulge flattening s 𝑠 s ; 11th Column : IC disc velocity ellipsoid ratio σ R 2 / σ z 2 superscript subscript 𝜎 𝑅 2 superscript subscript 𝜎 𝑧 2 \sigma_{R}^{2}/\sigma_{z}^{2} ; 12th Column : IC bulge rotation.

1st	2nd	3rd	4th	5th	6th	7th	8th	9th	10th	11th	12th
Name	$M_{\mathrm{disc},i}$	$M_{\mathrm{bulge},i}$	$N_{b,i}$	$c_{\mathrm{halo}}$	$a_{\mathrm{halo}}$	$h_{R,\mathrm{disc}}$	$z_{0,\mathrm{disc}}$	$a_{\mathrm{bulge}}$	$s$	${\sigma }_R^2/{\sigma }_z^2$	rotation
	$\left[{10}^9{M}_{\odot }\right]$	$\left[{10}^9{M}_{\odot }\right]$			$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$		disc	bulge
Y	$5$	–	$\mathrm{500\hspace{0.17em}000}$	9	30.2	1.5	0.10	–	–	2.0	–
P	$15$	–	$\mathrm{1\hspace{0.17em}500\hspace{0.17em}000}$	9	30.2	2.5	1.75	–	–	1.0	–
Q	$25$	–	$\mathrm{2\hspace{0.17em}500\hspace{0.17em}000}$	6.5	37.9	2.5	1.75	–	–	1.0	–
R	$15$	–	$\mathrm{1\hspace{0.17em}500\hspace{0.17em}000}$	6.5	37.9	2.5	1.75	–	–	1.0	–
U	$25$	–	$\mathrm{2\hspace{0.17em}500\hspace{0.17em}000}$	7.5	34.4	2.0	1.70	–	–	1.0	–
T	$20$	–	$\mathrm{2\hspace{0.17em}000\hspace{0.17em}000}$	9	30.2	2.0	1.70	–	–	1.8	–
C	–	$5$	$\mathrm{500\hspace{0.17em}000}$	9	30.2	–	–	0.45	1.15	–	no
K	–	$15$	$\mathrm{1\hspace{0.17em}500\hspace{0.17em}000}$	9	30.2	–	–	1.51	2.0	–	yes
M	–	$5$	$\mathrm{500\hspace{0.17em}000}$	9	30.2	–	–	1.51	2.0	–	yes
O	–	$15$	$\mathrm{1\hspace{0.17em}500\hspace{0.17em}000}$	9	30.2	–	–	1.51	3.0	–	yes
V	–	$5$	$\mathrm{500\hspace{0.17em}000}$	6.5	37.9	–	–	1.51	2.0	–	yes
W	$20$	$5$	$\mathrm{2\hspace{0.17em}500\hspace{0.17em}000}$	7.5	34.4	2.5	1.70	0.46	1.15	1.0	no
X	$20$	$5$	$\mathrm{2\hspace{0.17em}500\hspace{0.17em}000}$	7.5	34.4	2.0	1.70	0.70	1.15	1.8	no

Table 2. Table 2: List of models analysed in this paper. 1st Column : model name; 2nd Column : initial conditions; 3rd Column : final total baryonic mass M f subscript 𝑀 f M_{\rm f} (including initial baryonic mass); 4th Column : final time t f subscript 𝑡 f t_{\rm f} ; 5th Column : initial disc scalelength h R , i subscript ℎ 𝑅 i h_{R,{\rm i}} ; 6th Column : final disc scalelength h R , f subscript ℎ 𝑅 f h_{R,{\rm f}} ; 7th Column : scalelength growth parameter ξ 𝜉 \xi ; 8th Column : type of SFH law; 9th Column : exponential decay time-scale t SFR subscript 𝑡 SFR t_{\rm SFR} for the star formation rate; 10th Column : radial-to-vertical dispersion ratio for inserted particles λ 𝜆 \lambda ; 11th Column : prescription for initial velocity dispersion for inserted stellar particles, σ 0 ​ ( t ) subscript 𝜎 0 𝑡 \sigma_{0}(t) ; 12th-16th Column : parameters σ 1 subscript 𝜎 1 \sigma_{1} , σ 2 subscript 𝜎 2 \sigma_{2} , t 1 subscript 𝑡 1 t_{1} , t 2 subscript 𝑡 2 t_{2} and ι 𝜄 \iota , which determine σ 0 ​ ( t ) subscript 𝜎 0 𝑡 \sigma_{0}(t) ; 17th Column : GMC star formation efficiency ζ 𝜁 \zeta .

1st	2nd	3rd	4th	5th	6th	7th	8th	9th	10th	11th	12th	13th	14th	15th	16th	17th
Name	ICs	$M_f$	$t_f$	$h_{R,i}$	$h_{R,f}$	$\xi$	SF type	$t_{\mathrm{SFR}}$	$\lambda$	${\sigma }_0$	${\sigma }_1$	${\sigma }_2$	$t_1$	$t_2$	$\iota$	$\zeta$
		$\left[{10}^{10}{M}_{\odot }\right]$	$\left[\mathrm{Gyr}\right]$	$\left[\mathrm{kpc}\right]$	$\left[\mathrm{kpc}\right]$			$\left[\mathrm{Gyr}\right]$		type	$\left[\mathrm{km}{s}^{-1}\right]$	$\left[\mathrm{km}{s}^{-1}\right]$	$\left[\mathrm{Gyr}\right]$	$\left[\mathrm{Gyr}\right]$
Y1	Y	5	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.08
Y1$\zeta$-	Y	5	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.04
Y2	Y	5	10	2.5	2.5	0.0	1	8.0	1.0	const	6	–	–	–	–	0.08
P1$\zeta$-	P	5	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.04
P1s6	P	5	10	1.5	4.3	0.5	0	–	1.0	const	6	–	–	–	–	0.08
P1$\sigma$	P	5	10	1.5	4.3	0.5	1	8.0	1.0	const	10	–	–	–	–	0.08
P2	P	5	10	2.5	2.5	0.0	1	8.0	1.0	const	6	–	–	–	–	0.08
Q1	Q	6	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.08
Q1$\zeta$-	Q	6	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.04
R2$\sigma$	R	5	10	2.5	2.5	0.0	1	8.0	1.0	const	10	–	–	–	–	0.08
U1	U	6	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.08
U1$\sigma \zeta$*	U	6	10	1.5	4.3	0.5	1	8.0	1.0	const	10	–	–	–	–	0.06
W1	W	6	10	1.5	4.3	0.5	1	8.0	1.0	const	6	–	–	–	–	0.08
X5$\zeta$*	X	6	10	1.5	3.5	0.5	1	8.0	1.0	const	6	–	–	–	–	0.06
T5$\zeta$*	T	5.5	10	1.5	3.5	0.5	1	8.0	1.0	const	6	–	–	–	–	0.06
K2	K	5	10	2.5	2.5	0.0	1	8.0	1.0	const	6	–	–	–	–	0.08
O2	O	5	10	2.5	2.5	0.0	1	8.0	1.0	const	6	–	–	–	–	0.08
C$\alpha$2	C	5	10	2.5	2.5	0.0	1	8.0	1.0	atan	20	2	2.0	–	–	0.08
C$\beta$2s5	C	5	12	2.5	2.5	0.0	1	8.0	1.0	plaw	51	15	1.57	2.7	0.56	0.08
M$\alpha$1	M	5	10	1.5	4.3	0.5	1	8.0	1.0	atan	16	8	2.0	–	–	0.08
M$\alpha$1$\zeta$*	M	5	10	1.5	4.3	0.5	1	8.0	1.0	atan	16	8	2.0	–	–	0.06
M$\beta$1s5	M	5	12	1.5	4.3	0.5	1	8.0	1.0	plaw	51	15	1.57	2.7	0.47	0.08
M$\alpha$8s7	M	5	12	1.0	4.3	0.6	2	12.0	1.0	atan	16	8	3.5	–	–	0.08
V$\alpha$1	V	5	10	1.5	4.3	0.5	1	8.0	1.0	atan	16	8	2.0	–	–	0.08
V$\alpha$5$\lambda$	V	5	10	1.5	3.5	0.5	1	8.0	1.25	atan	16	8	2.0	–	–	0.08
V$\alpha$8s5	V	6	12	1.0	4.3	0.6	1	8.0	1.0	atan	16	8	2.0	–	–	0.08
V$\alpha$8s7	V	6	12	1.0	4.3	0.6	2	12.0	1.0	atan	16	8	3.5	–	–	0.08
V$\beta$8s5	V	6	12	1.0	4.3	0.6	1	8.0	1.0	plaw	51	15	1.57	2.7	0.47	0.08
V$\alpha$9s7$\lambda \zeta$*	V	6	12	1.0	3.5	0.6	2	12.0	1.25	atan	16	8	3.5	–	–	0.06
V$\alpha$9s8$\lambda \zeta$*	V	6	12	1.0	3.5	0.6	1	6.0	1.25	atan	16	8	2.0	–	–	0.06

Equations

$${\rho }_{\mathrm{DM}}(r)=\frac{M_{\mathrm{DM}}}{2\pi }\frac{a}{r{\left(r+a\right)}^3}.$$

$${\rho }_{\mathrm{disc},\mathrm{i}}(R,z)=\frac{M_{\mathrm{disc},\mathrm{i}}}{4\pi z_{0,\mathrm{disc}}{h_{R,\mathrm{disc}}}^2}{\mathrm{sech}}^2\left(\frac{z}{z_{0,\mathrm{disc}}}\right)\exp \left(\frac{-R}{h_{R,\mathrm{disc}}}\right),$$

$${\rho }_{\mathrm{bulge}}(R,z)=s{\rho }_{\mathrm{Hernquist}}\left(\sqrt{R^2+s^2{z}^2}\right).$$

$${⟨v_{\varphi }⟩}^2=k^2\left(⟨v_{\varphi }^2⟩-{\sigma }_R^2\right).$$

$${\sigma }_0(t)=\left(6+30{\mathrm{e}}^{-t/1.5\mathrm{Gyr}}\right)\mathrm{km}{\mathrm{s}}^{-1},$$

$${\sigma }_0(t)={\sigma }_1\left[\arctan \left(\frac{t_1-t}{1\mathrm{Gyr}}\right)+\frac{\pi }{2}\right]+{\sigma }_2(\mathrm{type}\mathrm{atan})$$

$${\sigma }_0(t)={\sigma }_1{\left(\frac{t+t_1}{t_2}\right)}^{-\iota }-{\sigma }_2(\mathrm{type}\mathrm{plaw})$$

$$\mathrm{SFR}(t)={\mathrm{SFR}}_0\times \exp (-t/t_{\mathrm{SFR}})(\mathrm{type}1),$$

$$\mathrm{SFR}(t)={\mathrm{SFR}}_0\times \exp (-t/t_{\mathrm{SFR}}-0.5\mathrm{Gyr}/t)(\mathrm{type}2),$$

$$h_R(t)=h_{R,\mathrm{i}}+(h_{R,\mathrm{f}}-h_{R,\mathrm{i}})(t/t_{\mathrm{f}}{)}^{\xi }.$$

$${\Sigma }_{\mathrm{GMC}}(\varphi )\propto {\left[{\Sigma }_{\mathrm{ys}}(\varphi )\right]}^{\alpha },$$

$$\rho (z,R=8\mathrm{kpc})={\rho }_0\left[\exp (-|z|/h_{\mathrm{thin}})+f\exp (-|z|/h_{\mathrm{thick}})\right].$$

$$A_2(R)\equiv \frac{1}{N(R)}\sum _{j=1}^{N(R)}{e}^{2\text{i}{\varphi }_j}$$

Full text

††pubyear: 2017††pagerange: The structural evolution of galaxies with both thin and thick discs–References

1 Introduction

The vertical density profile of stars in the Milky Way (MW) is fitted well by a sum of two exponentials (Gilmore & Reid, 1983). Jurić et al. (2008) find scaleheights of $\sim 900\,{\rm pc}$ for the geometrical thick disc and $\sim 300\,{\rm pc}$ for the thin disc. Similarly, the vertical surface brightness profiles of the majority of bright edge-on spiral galaxies show thin and thick components (Yoachim & Dalcanton, 2006). Studies of the chemical abundances of solar neighbourhood (Snhd) stars reveal that populations with hotter vertical kinematics and thus larger scaleheights have abundances of $\alpha$ elements relative to iron ($\left[{\rm\alpha/Fe}\right]$) that are larger than those of populations with small scaleheights and comparable iron abundances $\left[{\rm Fe/H}\right]$ (Fuhrmann, 1998; Bensby et al., 2003).

In the $\left[{\rm\alpha/Fe}\right]-\left[{\rm Fe/H}\right]$ plane, the low- and high-$\alpha$ populations are generally found to separate into fairly distinct sequences. Higher $\left[{\rm\alpha/Fe}\right]$ indicates shorter chemical enrichment time-scales and age determinations find systematically older ages for stars of the $\alpha$ enhanced chemical thick disc (Masseron & Gilmore, 2015; Martig et al., 2016a): they are generally found to be older than 8 Gyr. The distinct sequences have motivated models in which the two components formed in two temporally separated phases in very different conditions (e.g. Chiappini et al., 1997), but can also be explained as the result of continuous star formation and chemical enrichment (Schönrich & Binney, 2009).

It has now become clear that chemical and geometrical definitions of the thick disc yield different results. Whereas Jurić et al. (2008), who determined the density of all stars in the Snhd independent of chemistry, found that the geometrically thick disc has a longer radial scalelength than the thin disc, high-$\left[{\rm\alpha/Fe}\right]$ stars are found to form a thicker, but more centrally concentrated component than low-$\left[{\rm\alpha/Fe}\right]$ stars (Bovy et al., 2012; Hayden et al., 2015). A scenario in which the disc forms inside-out and each mono-age population flares, i.e. is thicker at outer than at inner radii, could potentially explain these observations (Minchev et al., 2015; Schönrich & McMillan, 2017). Flaring can be caused by radial migration of stars (Schönrich & Binney, 2012; Roškar et al., 2013) or by vertical heating of the outer disc through satellite interactions (Kazantzidis et al., 2008) or misaligned gas infall (Jiang & Binney, 1999).

Aumer et al. (2016a, b) (hereafter Papers 1 and 2) presented $\sim 100$ idealized $N$-body models of disc galaxies growing within live dark matter (DM) haloes over $\sim 10\,{\rm Gyr}$. These models covered a large variety of star formation and radial growth histories and most of them followed the assumption that all stars are born on near-circular orbits as in the MW today. Structural and kinematical properties of the MW’s thin disc, such as an exponential profile with scaleheight $\sim 300\,{\rm pc}$ or the local age-velocity dispersion relations could be reproduced if giant molecular clouds (GMCs) were included. Additionally, bars of comparable size and structure to that of the MW formed in these models. However, none of the models produced a realistic thick disc. The conclusion was that additional sources of heating were required early in the disc’s life, prior to the onset of thin-disc formation.

In this paper, we create models similar to those of Paper 1 that, in addition to a realistic thin disc, also contain an appropriate $\sim 10\,{\rm Gyr}$ old thick disc. Demanding the presence of an old thick component makes it much harder to steer a model to a configuration consistent with current data, in part because the number of observations that need to be explained simultaneously increases roughly twofold. Moreover, the thick disc must be steered into its present form by adjusting the conditions at the onset of disc formation, which will modify the subsequent formation of the thin-disc component by altering spiral and bar structures and their interaction with the dark halo. These changes to the thin disc and dark halo will themselves modify the appearance of the current thick disc. Moreover, the old chemical thick disc of the MW is centrally concentrated and should be important in the central $\sim 5\,{\rm kpc}$, a region dominated by the Galactic bar (Portail et al., 2017a), a structure that is supposed to have formed from a rather cold disc.

To create a thick disc we follow two approaches. In one scheme, we assume that the velocity dispersions of newborn stars decline with time. This is motivated by observations of redshift $z_{\rm rs}\sim 2$ galaxies, which show a high fraction of galaxies with H$\alpha$ kinematics consistent with ordered rotation, but significantly higher velocity dispersions than today’s disc galaxies (Förster Schreiber et al., 2009) and of galaxies at lower redshifts that indicate a continuous decline of gas velocity dispersion with decreasing redshift (Kassin et al., 2012; Wisnioski et al., 2015, but see Di Teodoro et al., 2016 for a different conclusion). Such a declining birth dispersion has also been found in hydrodynamical cosmological simulations of disc galaxy formation (Bird et al., 2013; Grand et al., 2016), and can be understood in terms of galaxies that gradually become less gas-rich and are characterized by gravitationally driven turbulence that decays (Forbes et al., 2012).

In an alternative scheme, we model thin+thick disc systems by creating thick initial conditions and growing thin discs within them. These models are thus representations of two-phase formation scenarios. We here do not model the formation of the thick disc prior to redshift $z_{\rm rs}\sim 2$, but several overlapping scenarios envisage the production of such an object: heating of an initially thinner disc by a merger (e.g. Quinn et al., 1993); formation in early gas-rich mergers (Brook et al., 2004); the formation in a clumpy, turbulent disc (Bournaud et al., 2009). It should be noted that the latter two scenarios could also be fitted into the picture of declining birth dispersions that constitutes our alternative modelling scheme.

In this paper, we concentrate on the setup of our models and on their structural evolution. We explore how models with double-exponential vertical profiles, circular speed curves like that of the MW and bars of appropriate structure can be constructed. We examine how the presence of a thick disc changes the preferred density of the dark halo and influences the evolution of the thin disc and the formation of a bar. We show that different scenarios for the formation of the thick disc leave signatures in the current distribution of age with the $(R,z)$ plane. A companion paper (Aumer et al., 2017, hereafter Paper 4) focuses on disc heating and radial migration in the models.

Our paper is structured as follows. Section 2 describes the setup and parameters of our simulations. Section 3 discusses the evolution of vertical density profiles and how this shapes the final age structure of the disc. Section 4 analyses the radial distribution of dark and baryonic mass components and Section 5 illustrates the evolution and structure of bars in the presence of thick discs. Section 6 discusses the successes and problems of our models. Section 7 concludes.

2 Simulations

The simulations analysed in this paper are similar to the models presented in Paper 1. These are simulations of growing disc galaxies within non-growing live DM haloes. They are run with the Tree code GADGET-3, last described in Springel (2005). We focus on standard-resolution models, which contain $N=5\times 10^{6}$ particles in the DM halo and a similar number of particles in the final stellar system. We here rely on collisionless simulations; models that contain an isothermal gas component were discussed in Papers 1 and 2 and shown to differ only mildly from collisionless models in terms of the structure and kinematics of the stellar component. In addition, and crucially, all simulations contain a population of short-lived, massive particles representing GMCs. Papers 1 and 2 demonstrated the importance of GMC heating in creating thin-disc components with realistic vertical structure as non-axisymmetric structure contributes little to vertical disc heating (see also Sellwood, 2013).

2.1 Initial conditions

Table 1 gives an overview of the initial conditions (ICs) of our models. The ICs were created using the GALIC code (Yurin & Springel, 2014). We refer the reader to Paper 1 for a full description of the IC creation process and here focus on how the ICs used in this paper differ from each other. All models discussed start with a spherical DM halo with a Hernquist (1990) profile

[TABLE]

The total mass of all ICs is $M_{\rm tot}=10^{12}\,M_{\odot}$, of which $M_{\rm b,i}=M_{\rm disc,i}+M_{\rm bulge,i}$ is in a stellar component and the rest $M_{\rm{DM}}$ is in DM. The inner density profile is adjusted so that it is similar to an NFW profile with concentration $c_{\rm halo}=6-9$ and thus scale radii are in the range $a=30-40\,{\rm kpc}$. The kinematics of the DM particles are initially isotropic (i.e. have equal velocity dispersions $\sigma_{r}=\sigma_{\phi}=\sigma_{\theta}$). Each DM halo studied here is resolved with $N_{\rm DM}=5\,000\,000$ particles.

ICs contain either a disc or a bulge component or both. IC discs have a mass in the range $M_{\rm{disc,i}}=5-25\times 10^{9}\,M_{\odot}$ and a density profile of the form

[TABLE]

where $h_{R,{\rm disc}}=1.5-2.5\,{\rm kpc}$ is the exponential disc scalelength. Radially constant isothermal vertical profiles with scaleheights $z_{0,{\rm disc}}=0.1-1.7\,{\rm kpc}$ are assumed. The ratios of the radial to vertical velocity dispersions are also radially constant in the IC discs and have values in the range ${\sigma_{R}^{2}}/{\sigma_{z}^{2}}=1.0-2.0$.

IC bulge components are set up with Hernquist density profiles with scalelengths $a_{\rm bulge}=0.4-1.5\,{\rm kpc}$ that are distorted to be oblate spheroids with axis ratios in the range $s=1-3$ as

[TABLE]

They are supposed to model either compact, non-rotating bulges or rotating spheroidal components. Following section 3.2 of Yurin & Springel (2014), the rotation of the axisymmetric bulge components is controlled via the Satoh (1980) parametrization

[TABLE]

Non-rotating bulges assume $k=0$, whereas rotating bulges are modelled as isotropic rotators with $k=1$.

Stellar particles in the ICs have a particle mass $m=1\times 10^{4}\,M_{\odot}$ and our ICs thus contain $N_{\rm b,i}=5-25\times 10^{5}$ stellar particles. As in Paper 1, the applied force softening lengths for the given resolutions are $\epsilon_{\rm b}=30\,{\rm pc}$ for baryonic particles (including GMCs) and $\epsilon_{\rm DM}=134\,{\rm pc}$ for DM particles.

The first letter of a model’s name specifies its ICs according to the scheme laid out in Table 1. In summary:

1. (i)

Y ICs have a compact and thin baryonic disc; 2. (ii)

P, Q, R, U and T ICs contain more extended and thicker discs, with varying size, thickness, mass and DM halo concentration; 3. (iii)

C ICs have a compact non-rotating bulge; 4. (iv)

W and X ICs combine a compact non-rotating bulge with a thick disc; 5. (v)

K, M, O and V ICs contain rotating oblate spheroids that crudely represent elliptical galaxies.

2.2 Growing the disc

To simulate the growth of galaxies over cosmological time-scales, stellar particles with a particle mass $m=1\times 10^{4}\,M_{\odot}$ are continuously added to the simulations. As we are interested in the coevolution of thin and thick discs, we pursue two different ideas, which were already touched upon in Paper 1. a) We assume that a thick and rotating stellar component was formed early on in a galaxy’s history (e.g. through a merger), represent this component with a thick disc or rotating spheroid in our ICs and over the timespan of the simulation add star particles to the system on near-circular orbits. b) We assume that the birth velocity dispersions of stars have been continuously declining over the history of a galaxy and add stellar populations with continuously decreasing dispersions.

2.2.1 Evolution of input velocity dispersion

In case a), the young stellar populations are assigned low birth velocity dispersions in the range $\sigma_{0}=\sigma_{1}=6-10\,{\rm km\,s^{-1}}$ in all three directions $R$, $\phi$ and $z$ as observed in the MW (see blue stars in Aumer & Binney, 2009). The mean rotation velocity $v_{\phi}(R)$ at radius $R$ is set to the circular velocity $v_{\rm circ}=\sqrt{{a_{R}(R)}R}$, where $a_{R}(R)$ is the azimuthal average of the radial gravitational acceleration, $\partial\Phi/\partial R$. As was shown in Papers 1 and 2, the in-plane dispersions quickly adjust to higher values and an appropriate ratio $\sigma_{\phi}/\sigma_{R}$ due to spiral and bar heating.

In case b) we choose a value of $\sigma_{0}$ that is declining with time. Paper 1 had applied

[TABLE]

which provided significantly too few hot stars. So, here we consider different functional forms:

[TABLE]

or

[TABLE]

In observations of high-redshift galaxies, it is common to study kinematics of the H$\alpha$ emission line and to assign a single characteristic velocity dispersion $\sigma_{\rm{H\alpha}}$ to each galaxy. Applying such a procedure, Wisnioski et al. (2015) find a dependence of observed H$\alpha$ dispersions $\sigma_{\rm{H\alpha}}\propto(1+z_{\rm rs})$, where $z_{\rm rs}$ is redshift (see also Kassin et al., 2012). The decline of our input velocity dispersion parameter $\sigma_{0}(t)$ for simulations of type ‘plaw’ (Equation 7) resulted from a rough approximation to this proportionality under the crude assumption that the kinematics of young stars follow H$\alpha$ kinematics. Note that this trend is observed for populations of galaxies at varying redshifts $z_{\rm rs}$. The spread in $\sigma_{\rm{H\alpha}}$ for each given redshift range in Wisnioski et al. (2015) is significant and the expected formation histories for galaxies are diverse, so that it is reasonable to assume that individual galaxies can show evolutions of velocity dispersions very different from $\sigma_{\rm{H\alpha}}\propto(1+z_{\rm rs})$. We thus also test decline histories of type ‘atan’ (Equation 6), which represents a scenario in which the transition from hot formation to cold formation is faster, yielding rather distinct hot and cold phases. Figure 1 visualizes the difference between assumed $\sigma_{0}$ histories.

We assume that $\sigma_{0}$ is radially constant and always set $\sigma_{z}(R,t)=\sigma_{0}(t)$. Resolution effects prevent high-$z_{\rm rs}$ observations from constraining the radial dependence of $\sigma_{\rm{H\alpha}}$ reliably, but the following two observational findings motivate this simple modelling approach. a) Newman et al. (2013) find that, if one characteristic value is assigned to each galaxy, $\sigma_{\rm{H\alpha}}$ depends little on the size of the galaxy. b)Jones et al. (2013) show radial $\sigma_{\rm{H\alpha}}(R)$ profiles of lensed galaxies and these show various shapes and high dispersions at outer radii.

In the case of high input dispersions, it is not desirable to use $\sigma_{\phi}/\sigma_{R}=1$ for an input, as this creates velocity distributions out of equilibrium. When $\sigma_{0}<10\,{\rm km\,s^{-1}}$, heating by non-axisymmetries is very efficient for young stars (see Paper 2) and the dispersions will quickly adjust to an appropriate ratio $\sigma_{\phi}/\sigma_{R}$, but heating after birth will play a minor role for high $\sigma_{0}$. Since $\sigma_{\phi}/\sigma_{R}$ depends on the shape of the rotation curve and the radial profile $\sigma_{R}(R)$ (see e.g. Binney & Tremaine, 2008, section 4.4.3), we apply a simple approach: We assume $\sigma_{\phi}=\sqrt{0.5}\sigma_{R}$ and an asymmetric drift correction for $\left<v_{\phi}\right>$ derived from Equation 4.228 in Binney & Tremaine (2008). We find that $\sigma_{\phi}/\sigma_{R}$ for these assumptions only adjusts mildly after insertion.

The choice for the ratio $\sigma_{z}/\sigma_{R}$ is also unclear. For high-$z_{\rm rs}$ discs, Genzel et al. (2011) studied $\sigma_{\rm{H\alpha}}$ as a function of disc inclination. Edge-on galaxies would thus be dominated by in-plane dispersions, whereas face-on galaxies would be dominated by vertical dispersions. They found mild indications for $\sigma_{R}>\sigma_{z}$, but Wisnioski et al. (2015) could not confirm this trend with a larger sample of galaxies. From the perspective of the thick-disc stars in the MW, Piffl et al. (2014) find $\sigma_{z}\approx\sigma_{R}$. We test values $\sigma_{R}=\lambda\sigma_{0}$ with $\lambda=1-1.3$.

As we place all newly added star particles in the midplane $z=0$, the measured vertical dispersion $\sigma_{z}$ of a young stellar component is smaller than $\sigma_{0}$, as the midplane is at the bottom of the vertical potential well and particles quickly lose kinetic energy moving away from it.

2.2.2 Star formation history

The star formation rate (SFR) is either constant (type 0) or

[TABLE]

or

[TABLE]

with $t_{\rm SFR}=6-12\,{\rm Gyr}$. Our simulations run for a total time of $t_{\rm f}=10-12\,{\rm Gyr}$ and the constant ${\rm SFR}_{0}$ is adjusted to produce at $t_{\rm f}$ a target final baryonic mass $M_{\rm f}$ in the range $5-6\times 10^{10}\,M_{\odot}$, including the mass of the stars in the ICs. Paper 1 showed that for the range of DM density profiles adopted for our models, galaxy masses $M_{\rm f}$ in this range provide the right level of self-gravity to explain the age-velocity dispersion relation of the Snhd. McMillan (2017) favours similar Galaxy masses for his MW mass models.

2.2.3 Radial growth history

Particles are added randomly distributed in azimuth every five Myr with an exponential radial density profile $\Sigma_{\rm SF}(R)\propto\exp(-R/h_{R}(t))$. The scalelength $h_{R}(t)$ of the newly added particles grows in time as

[TABLE]

To avoid inserting particles in the bar region, where near-circular orbits do not exist, particles are not added inside the cutoff radius $R_{\rm cut}$, which is determined by the current bar length (‘adaptive cutoff’; see Paper 1 for details).

2.3 GMCs

GMCs are modelled as a population of massive collisionless particles drawn from a mass function of the form ${\rm d}N/{\rm d}M\propto M^{\gamma}$ with lower and upper mass limits $M_{\rm low}=10^{5}\,M_{\odot}$ and $M_{\rm up}=10^{7}\,M_{\odot}$ and an exponent $\gamma=-1.6$. Their radial density is proportional to the star formation surface density $\Sigma_{\rm SF}(R)$, and their azimuthal surface density is given by

[TABLE]

where $\Sigma_{\rm ys}(\phi,R)$ is the surface density of young stars with ages 200-400 Myr and $\alpha=1$. The mass in GMCs is determined by the SFR efficiency $\zeta$. Specifically, for each $\Delta m_{\rm stars}$ of stars formed, a total GMC mass $\Delta m_{\rm GMC}=\Delta m_{\rm stars}/\zeta$ is created. GMC particles live for only $50\,{\rm Myr}$: for $25\,{\rm Myr}$ their masses grow with time as $m\propto t^{2}$, and for the final $25\,{\rm Myr}$ of their lives their masses are constant, before they disappear instantaneously. GMCs are added on orbits with $\sigma_{0}=6\,{\rm km\,s^{-1}}$. See Paper 1 for more details.

2.4 Overview of Models and their Naming

Table 2 gives an overview of the models discussed. We use an extended version of the naming convention described in Paper 1. All model names start with a capital letter identifying their IC according to the scheme defined by Table 1.

A Greek letter $\alpha$ or $\beta$ following the initial capital letter indicates that the model has declining birth velocity dispersion $\sigma_{0}$ rather than a thick disc in its ICs. In $\alpha$ models the decline of $\sigma_{0}$ follows an ‘atan’ shape (Equation 6) while in $\beta$ models it follows a ‘plaw’ shape (Equation 7). For simplicity the model names do not reflect the specific choices for the parameters of Equations (6) and (7), as they are of minor importance for the analyses in this paper.

These capital and, if present, Greek letters are followed by a number between 1 and 9 describing the radial growth history of the model, determined by parameters $h_{R,{\rm i}}$, $h_{R,{\rm f}}$ and $\xi$. Growth histories ‘1’, ‘2’ and ‘5’ were already used in Paper 1, ‘8’ ($h_{R,{\rm i}}=1.0\,{\rm kpc}$, $h_{R,{\rm f}}=4.3\,{\rm kpc}$ and $\xi=0.6$) and ‘9’ ($h_{R,{\rm i}}=1.0\,{\rm kpc}$, $h_{R,{\rm f}}=3.5\,{\rm kpc}$ and $\xi=0.6$) are new and represent inside-out growth from a very compact disc into an extended disc.

The final baryonic masses $M_{\rm f}$ of all models lie in the rather narrow range $5-6\times 10^{10}\,M_{\odot}$ and we do not include the variations in the naming convention. For all other parameters we define standard values and additional digits added to the model name only when a model deviates in one or more parameters from the standard. The meanings of the additional digits are:

• •

The overall star formation history (SFH) of a model is described by the SFR type, the final time $t_{\rm f}$ and the SF time-scale $t_{\rm SFR}$. Our standard choice is a type 1 SFR with $t_{\rm f}=10\,{\rm Gyr}$ and $t_{\rm SFR}=8\,{\rm Gyr}$. We have applied four additional SFHs, which are labelled by ‘s5’,…,‘s8’.

• •

The standard input velocity dispersion in models without declining $\sigma_{0}$ is $\sigma_{0}=6\,{\rm km\,s^{-1}}$ . Models with $\sigma_{0}=10\,{\rm km\,s^{-1}}$ are labelled as ‘$\sigma$’.

• •

The GMC star formation efficiency $\zeta$ has a standard value of 0.08. Models with $\zeta=0.04$ are labelled as ‘$\zeta$-’, and models with $\zeta=0.06$ are labelled as ‘$\zeta$*’.

• •

The standard choice for the radial-to-vertical dispersion ratio for inserted particles is $\lambda=1$. Models with $\lambda=1.25$ are labelled as ‘$\lambda$’.

3 Vertical Profiles

Jurić et al. (2008) studied the vertical stellar density profile in the Snhd and presented a bias-corrected model fit to Sloan Digital Sky Survey data of the form

[TABLE]

They found $h_{\rm thin}=300\,{\rm pc}$ and $h_{\rm thick}=900\,{\rm pc}$ with 20 per cent uncertainty each and $f=0.12$ with 10 per cent uncertainty. $f$ here is the ratio of local densities, whereas $f_{\Sigma}=fh_{\rm thick}/h_{\rm thin}=0.36$ is the ratio of local surface densities. Bland-Hawthorn & Gerhard (2016) show that while most studies in the literature find similar scaleheights $h_{\rm thick}$ and $h_{\rm thin}$ to Jurić et al. (2008), this determination ranks at the upper end in terms of $f_{\Sigma}$ for photometric surveys, the average literature value of which is significantly lower, $f_{\Sigma}=0.12$. This strong variation in $f_{\Sigma}$ between different studies is likely caused by significant differences in the survey selection functions and degeneracies between $f_{\Sigma}$ and the scaleheights.

3.1 Vertical profile shape

In the models of Paper 1, GMC heating created remarkably exponential vertical profiles with scaleheights $h_{\rm thin}=200-350\,{\rm pc}$. Thus these models can reproduce the scaleheight of the MW’s thin disc but fail to create a realistic thick disc. In Figure 2 we examine the vertical profiles of the present models at $t=t_{\rm f}$ and $R=8\pm 0.5\,{\rm kpc}$ as black symbols. We overplot in red fits of Equation (12) to these profiles with the scaleheights given in the top-left corner of each panel. The values of $f$ and $f_{\Sigma}$ are given in the top-right corner of each panel. The figure illustrates that both approaches to producing thin+thick disc systems presented here produce double-exponential profiles similar to the one observed in the MW.

In the upper row of Figure 2 we show profiles for models with thick-disc ICs. We find that using disc-like ICs with $z_{0}\sim 1.7\,{\rm kpc}$ yields values of $h_{\rm thick}$ in the range $815-1133\,{\rm pc}$ at $t_{\rm f}$ similar to the one inferred for the MW. As was already shown in Paper 1, $z_{0}\sim 1\,{\rm kpc}$ yields final thick discs that are too thin. For elliptical ICs we find that setting the axis ratio $s=2$ in model K2 yields $h_{\rm thick}\approx 2.4\,{\rm kpc}$ and $s=3$ in model O2 yields $h_{\rm thick}\approx 1.4\,{\rm kpc}$, which are both too thick. We thus decided to focus on disc-like thick ICs.

In terms of thin-disc scaleheights, we find values $h_{\rm thin}=200-339\,{\rm pc}$. As was already discussed in Paper 1, for a given GMC mass function, lower values of the star formation efficiency $\zeta\sim 0.05$ and thus a higher total mass in GMCs per mass of formed stars are required to obtain $h_{\rm thin}\sim 300\,{\rm pc}$ as in the MW. Lowering the DM halo concentration reduces the vertical force contribution from the halo and thus also mildly increases $h_{\rm thin}$. The highest values for both $h_{\rm thick}$ and $h_{\rm thin}$ for disc IC models are found for model U1 that at $t_{\rm f}$ features an overly extended and thick bar as is discussed in Section 5. As was discussed in Paper 1, bars that extend beyond $R=5\,{\rm kpc}$ can significantly thicken vertical profiles at $R=8\,{\rm kpc}$.

For the density and surface density ratios in models with thick-disc ICs, we find ranges of $f=0.04-0.167$ and $f_{\Sigma}=0.14-0.53$, which include the Jurić et al. (2008) values and are in the upper half of the values of the Bland-Hawthorn & Gerhard (2016) literature compilation. Naturally, increasing the thick-disc mass $M_{\rm disc,i}$ at fixed thick-disc scalelength $h_{R,{\rm disc}}$ increases these ratios (Q versus P models) and decreasing $h_{R,{\rm disc}}$ at fixed $M_{\rm disc,i}$ lowers them (U versus Q models). As far as other model parameters are concerned there are no clear patterns apparent. This is likely connected to competing effects. For example, few GMCs produce less vertical heating, but lead to stronger bars, which, if long enough, can thicken the vertical profile.

The lower row of Figure 2 shows models with declining birth dispersions. We find values for the scaleheights in the ranges $h_{\rm thin}=214-284\,{\rm pc}$ and $h_{\rm thick}=863-1234\,{\rm pc}$, very similar to the ranges found for thick IC disc models. Model C$\alpha$2 has the thinnest thin disc as at late times it has a very low value of $\sigma_{0}<5\,{\rm km\,s^{-1}}$ and it also has a high star formation efficiency $\zeta=0.08$. V$\alpha$9s8$\lambda\zeta$* has $\zeta=0.06$ and $\sigma_{0}\sim 10\,{\rm km\,s^{-1}}$ at late times and thus the highest $h_{\rm thin}$ among these models. V$\alpha$8s5 and V$\beta$8s5 only differ in the shape of the declining $\sigma_{0}$ curve: V$\alpha$8s5 is of type ‘atan’, whereas V$\beta$8s5 is of type ‘plaw’ (see Equations 6 and 7). Their final vertical profiles are rather similar, so both types of decline are acceptable. The value of $h_{\rm thick}$ is also mildly influenced by the ICs, which for V and M models is a low-mass $s=2$ elliptical. As models V$\alpha$8s5/V$\beta$8s5 have a higher final mass than M$\alpha$1 and thus a lower fraction of IC stars, their thick discs appear mildly thinner.

The corresponding thin-to-thick disc ratios are $f=0.025-0.055$ and $f_{\Sigma}=0.11-0.22$, which are lower than the Jurić et al. (2008) values, but comfortably within the range found in the literature. For models with declining $\sigma_{0}(t)$ the thick-disc mass fraction depends on the fraction of mass formed during early formation stages with high $\sigma_{0}$ and thus on the detailed forms of $\sigma_{0}(t)$ and the SFH. Moreover, the radial growth history determines how many stars are formed at a certain radius during this period. Dynamical heating and migration processes also influence the number of old stars found at $t_{\rm f}$, so the final value of $f$ is not easily predicted.

3.2 Vertical profile evolution with time

In Figure 3 we analyse the temporal evolution of the shapes of vertical profiles in two different types of simulations. In the models of Paper 1 that included GMC heating, thin-disc vertical profiles are at all times exponential and their scaleheights change very little with time. This finding reflects balance between mass growth, which continuously supplies cold particles and deepens the vertical potential well, and GMC heating which efficiently increases the vertical velocity dispersions of young stars.

The thick-disc IC model P1$\sigma$ shows a double-exponential profile from early on. The thin-disc part of its profile is very constant, just like the thin-disc-only models of Paper 1. What changes are the surface density ratio $f_{\Sigma}$, which by construction becomes more and more thin-disc-dominated, and the scaleheight of the thick disc, which becomes smaller with time, because the growth in the thin disc’s mass deepens the vertical potential well and the thick disc is not heated significantly.

M$\alpha$1 represents models with declining $\sigma_{0}(t)$. Its profile at early times is closer to a single- than a double-exponential, as there are no cold, thin-disc populations present. Only as $\sigma_{0}(t)$ falls below $\sim 20\,{\rm km\,s^{-1}}$ does a thin disc build up. The scaleheight of the thin disc $h_{\rm thin}$ becomes smaller with time as the decline in $\sigma_{0}(t)$ cools the thin-disc population as a whole despite the vertical heating due to the GMCs. The deepening of the vertical potential well adds to that effect and reduces $h_{\rm thick}$ in the same way as in P1$\sigma$.

3.3 Radial dependence of vertical profiles

The vertical profiles of observed disc galaxies are very constant radially (van der Kruit & Searle, 1982). Paper 1 showed that in thin-disc-only models that include GMC heating, vertical profiles are almost independent of radius $R$, unless there is a buckled bar, which thickens only the central region. Figure 4 shows the radial variations of the vertical profiles in five of our thin+thick models. The inner parts of models U1$\sigma\zeta$* and O2 show thickening by a bar. U1$\sigma\zeta$* otherwise shows a profile that is almost constant in radius.

Model P1$\sigma$ shows a thin disc that becomes mildly thicker and has a higher mass fraction towards larger $R$, whereas $h_{\rm thick}$ stays roughly constant. This is caused by inside-out formation, on account of which the young thin disc has a longer scalelength than the thick disc, and higher $\sigma_{0}=10\,{\rm km\,s^{-1}}$ than the models studied in Paper 1. Combined with a shallower potential well at outer radii, where heating is limited, this value of $\sigma_{0}$ yields thicker thin discs.

T5$\zeta$* shows a stronger fading and a mild thickening of the thick populations towards the outskirts. Compared to P ICs, the T ICs have a more compact and more massive thick disc and also a bulge component. The effective vertical profile of the two components thus varies with radius already in the ICs. O2 has elliptical ICs, which lead to a thick population that becomes thicker with increasing $R$ and also attains a higher mass fraction in the outskirts.

Model M$\alpha$1 with declining $\sigma_{0}$ shows a thin disc that thickens mildly with $R$, because at late times it has $\sigma_{0}\sim 10\,{\rm km\,s^{-1}}$ and thus behaves similarly to model P1$\sigma$ discussed above. The thick disc becomes thicker with $R$. This is characteristic for models of this type. It is a consequence of our assumption that $\sigma_{0}$ is constant with radius. This leads to a flaring of the hot component, which is hardly influenced by vertical heating.

Note that in contrast to the thick discs in models with declining $\sigma_{0}$, the thick-disc ICs of our alternative modelling scheme were set up with a radially constant scaleheight and thus a vertical velocity dispersion that declines with $R$.

3.4 Flaring of mono-age components

In Figure 5 we examine for a selection of models how the median distance from the midplane $|z|_{\rm med}$ varies with radius $R$ for populations of different ages. It has been suggested that the lack of significant changes with disc radius $R$ in the double-exponential vertical mass profiles of disc galaxies is a consequence of inside-out growth combined with $|z|_{\rm med}$ being an increasing function of $R$ for all mono-age populations (‘flaring’; Minchev et al., 2015).

Flaring can be caused by satellite interactions or misaligned infall of gas, which are not present in our models. However, Paper 1 showed that disc galaxies formed in isolation also have flaring mono-age components. The amount of flaring is determined by the radial mass profile of the disc and the radial profile of the vertical velocity dispersion of mono-age components $\sigma_{z}(R,\tau)$, which is determined by the birth dispersions of stars and the vertical heating mechanism(s) at work. Radial migration of stars can also influence $\sigma_{z}(R,\tau)$ (Schönrich & Binney, 2012).

The left most panel of Figure 5 depicts $|z|_{\rm med}(R)$ for the mono-age populations of model P1$\sigma$, which has thick-disc ICs. The thin-disc populations are at all radii substantially thinner than the old thick-disc stars. These mono-age populations all show flaring, and for the youngest populations $|z|_{\rm med}$ increases from $R=1$ to $15\,{\rm kpc}$ by up to a factor of $\eta\equiv|z|_{\rm med}(15\,{\rm kpc})/|z|_{\rm med}(1\,{\rm kpc})\sim 5$. In fact, in P1$\sigma$ the structure of the thin-disc populations is similar to that in thin-disc-only model Y1 examined in Paper 1. The relative increase in $|z|_{\rm med}$ with $R$ becomes smaller with increasing age, and we find $\eta\sim 2.5$ for the oldest inserted stars. In model U1$\sigma\zeta$* (second panel in Figure 5) the situation is altered by a vertically extended bar. On account of bar buckling, $|z|_{\rm med}$ for intermediate-age populations now peaks around $R\sim 3\,{\rm kpc}$, then declines slightly to $R\sim 8\,{\rm kpc}$, and gradually increases further outward. The youngest populations have not been affected by bar buckling and have $\eta\sim 3$.

In all our models, non-IC stars at any given time $t$ are inserted with a radially constant birth dispersion $\sigma_{0}(t)$. On account of the outward decrease in surface density, this results in flaring. GMC heating increases $\sigma_{z}$ more strongly in the centre than in the outskirts and thus flattens the increase in $|z|_{\rm med}(R)$. For models without buckled bars and with $\zeta=0.04$ and thus more GMCs per unit mass of inserted stars, $\eta$ can be as low as $\sim 1.5$ for old inserted stars. Bars can additionally heat the central regions of galaxies (see e.g. Grand et al., 2016) and thus cause even flatter runs of $|z|_{\rm med}(R)$.

The three panels on the right of Figure 5 show models with declining $\sigma_{0}$, and in these models the structure of $|z|_{\rm med}(R)$ is quite different from what it is in the models with thick IC discs. Now the low-mass elliptical ICs are unimportant because the thick disc is formed mainly by old added stars. Due to the continuous decline in $\sigma_{0}(t)$, the curves form a continuum rather than a bimodal grouping. The young, thin-disc populations behave very similarly to those in the models with thick-disc ICs, the local maxima in $|z|_{\rm med}$ at $R\sim 2\,{\rm kpc}$ in model V$\alpha$8s5 being caused by a buckled bar.

The thick-disc stars in models with declining $\sigma_{0}$ (red and orange curves) yield a similar value of $\eta$ to the youngest thin-disc stars in all models. Again this reflects our decision to make the declining birth dispersion $\sigma_{0}$ independent of $R$. Given that vertical GMC heating has little influence on the thick-disc populations, and that the depth of the vertical potential well declines with $R$, strong flaring is an inevitable consequence. The flattest curves are thus found for the intermediate-age populations, which were already born on relatively cold orbits and have been significantly affected by GMC and bar heating. The differences between models M$\alpha$1 and M$\beta$1s5 can be explained by the different shapes of $\sigma_{0}(t)$ applied (see Figure 1). Whereas the oldest component of each model flares in a similar way, the two next-oldest populations show stronger flaring in the atan model M$\alpha$1 than in the plaw model M$\beta$1s5, in which $\sigma_{0}$ declines more gradually. Also, in the atan model $|z|_{\rm med}$ declines with decreasing age faster than in the plaw model.

The flaring of all disc populations discussed so far can be qualitatively explained by birth-dispersion profiles and disc heating mechanisms. A comparison between the thick-disc IC stars in models P1$\sigma$ and U1$\sigma\zeta$* indicates that radial migration plays a role as well. The IC stars in P1$\sigma$ only show a strong outward increase in $|z|_{\rm med}$ in the centre and hardly any flaring at larger radii, whereas in U1$\sigma\zeta$, the flaring of the IC stars in the outer disc is stronger. The thick-disc ICs were created with a radially constant scaleheight and are hardly affected by vertical heating mechanisms, so their curve of $|z|_{\rm med}(R)$ at $t_{\rm f}$ is determined by the change in the vertical potential well together with extent to which their stars migrate radially. Bar formation funnels a lot of mass to the centre and consequently the disc’s thickness decreases there. The stronger flaring in U1$\sigma\zeta$ is likely explained by higher levels of radial migration.

When stars migrate radially, their vertical actions $J_{z}$ are conserved (Solway et al., 2012). Populations of stars that are born with a radially constant scaleheight in a MW-like disc have their mean $J_{z}$ decreasing with $R$. Consequently, outward migrators at a given $R$ have higher $J_{z}$ (and thus higher $\sigma_{z}$) than non- or inward migrators. If there are more outward than inward migrators, as expected for the outer disc regions, such populations are expected to flare (Schönrich & Binney, 2012). This principle is complicated by the finding that stars migrate less if they have high $J_{z}$ (Vera-Ciro et al., 2014), but Paper 4 shows that for our thin-disc populations and for stars from thick-disc ICs at $t=t_{\rm f}$, outward migrators are indeed more numerous in the outer disc and have higher $\sigma_{z}$ than inward migrators (see also Roškar et al., 2013). Model U1$\sigma\zeta$* has both a more compact IC disc and a lower density dark halo than model P1$\sigma$, which leads to stronger non-axisymmetries and a higher fraction of outward migrators in the outer disc (Paper 4), which in turn explains the stronger flaring of IC stars in U1$\sigma\zeta$*.

Compared to the oldest inserted stars in models M$\alpha$1 and M$\beta$1s5, in model V$\alpha$8s5 this population shows a much flatter curve $|z|_{\rm med}(R)$ at $R>5\,{\rm kpc}$. Model V$\alpha$8s5 has a lower density dark halo and a more compact disc at early times than the M models. Hence in this model the oldest inserted stars in the outer disc have a higher fraction of outward migrators. For models with declining $\sigma_{0}$, thick-disc stars are born with radially constant $\sigma_{z}$ and thus their mean $J_{z}$ increasing with $R$. As they are hardly affected by disc heating, in these models old outward migrators at a given $R$ have lower $J_{z}$ and thus lower $\sigma_{z}$ than inward migrators. This likely explains the flatter curve $|z|_{\rm med}(R)$ for the oldest inserted stars in model V$\alpha$8s5. So depending on the shape of $\sigma_{z}(R)$ at birth, radial migration can both strengthen and weaken the flaring of a mono-age population. A detailed analysis of radial migration in our models is presented in Paper 4.

3.5 Radial and vertical age structure

The combination of recent and ongoing astrometric and spectroscopic surveys of MW stars is about to increase vastly our knowledge of the age structure in the Galactic disc(s) (e.g. Martig et al., 2016a). In the top row of Figure 6 we therefore show maps of median age $\tau_{\rm med}$ as a function of $R$ and $|z|$ for various models. To do so we assume that the oldest stars in all models are $13\,{\rm Gyr}$ old and thus randomly assign ages in the range $[t_{f},13\,{\rm Gyr}]$ to IC star particles.

The leftmost panel shows the thin-disc-only model Y1 from Paper 1. The innermost region $R<4\,{\rm kpc}$ at all altitudes and high altitudes $|z|>1.5\,{\rm kpc}$ at all radii are dominated by old stars. The youngest $\tau_{\rm med}$ are found at high $R$ and low $|z|$. Intermediate $\tau_{\rm med}$ are confined to $|z|<200\,{\rm pc}$ at $R<6\,{\rm kpc}$, to $|z|<500\,{\rm pc}$ at $R<10\,{\rm kpc}$, but can populate regions up to $|z|\sim 1500\,{\rm pc}$ at $R\sim 15\,{\rm kpc}$. This characteristic $\tau_{\rm med}$ pattern is caused by a combination of inside-out formation and disc heating. The oldest population is more compact and thicker, whereas the younger population lives closer to the midplane and preferentially at larger radii, where it extends to higher $|z|$.

Model P1$\sigma$ has an old thick disc, which is more massive, much thicker and more extended than the oldest disc population in Y1. This changes the $\tau_{\rm med}$ map only marginally. The high-$|z|$ populations are by construction older than in Y1 and the young, outermost populations are also mildly older than their counterparts in Y1. The $\tau_{\rm med}$ pattern is however rather similar.

P2 differs from P1$\sigma$ in lacking inside-out formation, which makes the outer populations older. It also has lower $\sigma_{0}$, which reduces the flaring of the outer populations. Consequently, at all radii $R$, all altitudes $|z|>500\,{\rm pc}$ are dominated by old stars and the zone populated by intermediate $\tau_{\rm med}$ is confined to $|z|>300\,{\rm pc}$ at $R<10\,{\rm kpc}$ and to $|z|>500\,{\rm pc}$ at all $R$.

Model U1 is an inside-out model, but compared to P1$\sigma$ has a lower concentration halo, a more massive and more compact IC thick disc and lower $\sigma_{0}$. It has a thicker and longer bar than the other depicted models. The buckled bar causes an area at $R<5\,{\rm kpc}$, which is populated by stars with old $\tau_{\rm med}$ and shows no noticeable vertical age gradient. The bar also heats the disc at $R\sim 5-10\,{\rm kpc}$ vertically and thus increases the altitudes at which younger stars are found, so the radial increase of the maximum $|z|$ at which intermediate $\tau_{\rm med}$ are found is flatter in this model.

Model M$\beta$1s5 is a model with declining $\sigma_{0}$. It shows a $\tau_{\rm med}$ pattern that differs from those of the other models shown in Figure 6. At $R<5\,{\rm kpc}$, its vertical age structure is similar to that in the P models, but at $R>5\,{\rm kpc}$ the intermediate-age stars reach higher altitudes, which leads to declining radial age gradients at all $|z|$. This is caused by our assumption of radially constant $\sigma_{0}(t)$, as was also discussed in relation to the flaring of mono-age populations shown in Figure 5. Moreover, vertical heating is inefficient for stars with high $\sigma_{0}$ and for stars at large radii.

Martig et al. (2016b) have recently presented measurements of radial gradients in $\tau_{\rm med}$ at various altitudes $|z|$ in the Snhd. At $|z|>500\,{\rm pc}$, they find significant declines in $\tau_{\rm med}$ with radius $R$ at all radii. However, the value of $\partial\tau_{\rm med}/\partial R$ is still very uncertain. Close to the plane near the solar radius $R_{0}$, $\tau_{\rm med}(R)$ is rather flat. In the middle row of Figure 6 we show $\tau_{\rm med}$ versus $R$ at four ranges in $|z|$.

As model Y1 lacks a thick disc, we will not discuss it in detail. The P and U models with thick IC discs clearly show no decline in $\tau_{\rm med}$ with $R$ at $|z|\sim 1.5\,{\rm kpc}$ and only P1$\sigma$ shows a negative ${\rm d}\tau_{\rm med}/{\rm d}R$ at $|z|\sim 1.0\,{\rm kpc}$ and $R>10\,{\rm kpc}$ due to inside-out formation and a stronger flaring in the young disc due to higher $\sigma_{0}$. By contrast, M$\beta$1s5 shows clear negative gradients throughout the whole disc at all latitudes and is thus more in agreement with the age determinations of Martig et al. (2016b) .

At $|z|\sim 0.5\,{\rm kpc}$, all models show a negative radial $\tau_{\rm med}$ gradient. This gradient is weaker in U1 due to the influence of the unrealistically long and thick bar. In the midplane, the inside-out models also show clear negative age gradients, whereas P2, which has a constant radial feeding scalelength, shows a rather flat age profile, as does U1, again strongly affected by the long bar. As the age gradients in M$\beta$1s5 and P1$\sigma$ are flatter closer to the midplane and the observations are still very uncertain, little can be deduced yet about inside-out growth.

Casagrande et al. (2016) have presented evidence for a vertical age gradient in the Snhd from asteroseismology. They find a decline by $\sim 4\pm 2\,{\rm Gyr}/\,{\rm kpc}$, but have little knowledge of the shape of the decline. In the lower row of Figure 6 we present $\tau_{\rm med}$ as a function of $|z|$ at three radii, the cyan line representing a solar-like radius and the dashed line showing a constant vertical gradient of $4\,{\rm Gyr}/\,{\rm kpc}$.

Due to the dominance of the thick disc at $|z|>1\,{\rm kpc}$ and the weak radial age gradient, at all radii in P2, age increases more strongly than at $4\,{\rm Gyr}/\,{\rm kpc}$ up to $|z|\sim 1\,{\rm kpc}$ and then flattens out. In P1$\sigma$ at $R=13\,{\rm kpc}$, the enhanced presence of younger stars away from the plane results in a flatter increase of $\tau_{\rm med}$ with $|z|$, but at $R=8\,{\rm kpc}$ the situation is similar to that in P2. In U1 the long bar causes an almost flat $\tau_{\rm med}$ versus $|z|$ plot in the bar region and a flatter gradient at $R=8\,{\rm kpc}$. In M$\beta$1s5, $\tau_{\rm med}(|z|)$ at $R=8\,{\rm kpc}$ flattens more gradually and at $R=13\,{\rm kpc}$ has an almost constant slope due to the stronger flaring of mono-age populations.

Averaged over the studied vertical extent of $2\,{\rm kpc}$, all models show a vertical gradient consistent with the findings of Casagrande et al. (2016). As they observe very few stars above $|z|>1.2\,{\rm kpc}$, model M$\beta$1s5 shows the best agreement with the still very uncertain data.

4 Radial mass distribution

In this section we investigate the radial distribution of baryonic and dark matter in our models and the circular speed curves $v_{\rm circ}(R)$ that result from them.

4.1 Dark matter density

As discussed in Paper 1 the parameters for our DM haloes as set up in the ICs are motivated by what $\Lambda$ cold dark matter ($\Lambda$CDM) predicts for haloes associated with MW mass galaxies. The DM profile $\rho_{\rm DM}(r)$ will be modified by growing a massive baryonic disc within the DM halo and by interaction with non-axisymmetric disc structures such as the bar and spirals. As the halo is always spherical in the ICs, but the disc mass fraction in the ICs varies strongly between ICs and the various galaxy models evolve differently, the final haloes differ even if two models share the same DM IC parameters.

Our best constraints on the DM content in the MW come from dynamical measurements of the total matter surface density in the Snhd. Subtracting the baryonic components, McKee et al. (2015) find a local DM density of $\rho_{\rm DM}=0.013\pm 0.003\,M_{\odot}\,{\rm pc}^{-3}$. In Figure 7 we plot in red the DM density profiles $\rho_{\rm DM}(R)$ as measured in the midplane of the galaxy at $t=t_{\rm f}$ and compare them to the Snhd constraints assuming that the solar Galactic radius is $R_{0}=8.3\pm 0.3\,{\rm kpc}$ (Schönrich, 2012).

We find that at $t_{\rm f}$ all models fall within the constraints of McKee et al. (2015). The four shown in Figure 7 are a representative selection. The main drivers for $\rho_{\rm DM}(R_{0})$ at $t_{\rm f}$ are, as expected, initial halo concentration $c$ and the added disc mass $M_{\rm add}=M_{\rm f}-M_{\rm disc,i}-M_{\rm bulge,i}$. Consequently R2$\sigma$, which has $c=6.5$ and $M_{\rm add}=3.5\times 10^{10}\,M_{\odot}$ has the lowest $\rho_{\rm DM}(R_{0})$, whereas M$\alpha$1 with $c=9$ and $M_{\rm add}=4.5\times 10^{10}\,M_{\odot}$ has the highest $\rho_{\rm DM}(R_{0})$ among the models shown. P2 with $c=9$ and $M_{\rm add}=3.5\times 10^{10}\,M_{\odot}$ and V$\alpha$8s5 with P2 with $c=6.5$ and $M_{\rm add}=5.5\times 10^{10}\,M_{\odot}$ show intermediate $\rho_{\rm DM}(R_{0})$. It is worth noting that the models with $c=4$ presented in Paper 1 and discarded because of overly strong bars indeed show too low $\rho_{\rm DM}(R_{0})$.

We also note that at $t_{\rm f}$ none of our models shows a cored DM profile in the centre, as was recently favoured by Cole & Binney (2017). Our IC DM profiles do not contain a core as is indicated by the pink dashed lines. The DM densities $\rho_{\rm DM}(R_{0})$ of the ICs are significantly lower than in the final models. Initial profiles with $c=9$ lie at the lower allowed limit for today’s Snhd and models with $c=6.5$ are clearly below this limit. During the simulations they are altered by compression due to the added mass in stars and by angular momentum transfer from stars to DM due to bars and spirals. As all models have declining ${\rm SFR}(t)$ and bars form in the later evolution stages as shown in Section 5, the increase in the DM density at $R<15\,{\rm kpc}$, where the disc grows, is strong up to $t=0.3t_{\rm f}$ as indicated by the green dashed lines and rather weak afterwards. In the four models shown, the relative increase in $\rho_{\rm DM}(R<15\,{\rm kpc})$ is strongest in models V$\alpha$8s5 as it has the largest baryonic mass fraction and the largest $M_{\rm add}$, and the increase is weakest in P2. Angular momentum transfer to the halo by spirals and the bar is not strong enough to create cores, as was also shown by Sellwood (2008).

4.2 Solar Neighbourhood surface density

Paper 1 showed that despite having control over the evolution of the input scalelength $h_{R}(t)$, there was little control over the final surface density profile of the models. The more compactly a disc was fed, the earlier it grew a bar, which redistributed matter and, as we avoid inserting particles into the bar region, shifted the inner cutoff radius outwards. In the end, the surface density profiles of a range of models with different radial growth histories were thus rather similar. Our surface density was thus decided by the total mass of the final model, which we justified from the reasonable agreement of our models with a) the vertical scaleheight of the thin disc, b) an appropriate local circular speed, c) an appropriate amount of radial migration to $R_{0}$ and d) appropriate vertical and radial velocity dispersions.

In Figure 8 we examine how well our models fulfil constraints on the Snhd baryonic surface density $\Sigma_{\rm b}(R_{0})$. Table 3 of McKee et al. (2015) gives an overview of determinations of $\Sigma_{\rm b}(R_{0})$ including gas of all phases, stars and stellar remnants. $\Sigma_{\rm b}$ is consistently found to be in the range $\Sigma_{\rm b}(R_{0})=40-60\,M_{\odot}\,{\rm pc}^{-2}$, where the given errors are included in the interval. Gas is found to contribute $25-30$ per cent. As in our simulations gas is only represented by GMCs, our gas fractions are much lower.

We choose $\Sigma_{\rm b}(R)$ over the stellar surface density $\Sigma_{\star}(R)$, as we are interested in the connection between disc structure and kinematics, and the strength of non-axisymmetries and thus the levels of radial disc heating and radial migration are determined by $\Sigma_{\rm b}(R)$. Moreover, the interplay between vertical profiles and vertical velocity dispersions depends on the total mass surface density and not only on $\Sigma_{\star}(R)$ and it is therefore appropriate to compare $\Sigma_{\rm b}$ in models and observations, although the division of mass between gas and stars is very different. It should be noted that the neutral hydrogen component, which is missing in our models, will have a smaller scaleheight than the stars and thus a model that has the right vertical profile, kinematics and DM halo is not expected to agree with the Snhd $\Sigma_{\rm b}(R_{0})$.

The local radial exponential scalelength $h_{R}(R_{0})$ of the MW is rather uncertain. Licquia & Newman (2016) recently compiled a variety of measurements in the optical and infrared, the vast majority of which fall in the range $2-4\,{\rm kpc}$. Their meta-analysis of 29 previous measurements yields an estimate of $h_{R}(R_{0})\sim 2.65\,{\rm kpc}$. Bovy et al. (2012) showed that populations of stars with different chemical abundances show widely varying scalelength, the most compact of which have $h_{R}\sim 1.5\,{\rm kpc}$ and the most extended of which are consistent with locally flat profiles. As we are plotting $\Sigma_{\rm b}(R_{0})$ that includes GMCs, the comparison is not exact, but because at $t=t_{\rm f}$ the GMC mass fractions are $2-3$ per cent as for molecular gas in the MW today, the correction is negligible for our purposes. A more relevant question is whether the missing neutral gas mass, which is a highly relevant mass component in the outer MW disc, is properly represented in our models.

In the upper row of Figure 8 we plot $\Sigma_{\rm b}(R)$ for various models and overplot a blue box indicating $\Sigma_{\rm b}(R_{0})=40-60\,M_{\odot}\,{\rm pc}^{-2}$ at $R_{0}=8.3\pm 0.3\,{\rm kpc}$ and a dashed line indicating an exponential with $h_{R}(R_{0})=2.65\,{\rm kpc}$. Due to the connection of $\Sigma_{\rm b}(R)$ to bars, we also plot the $m=2$ Fourier amplitude

[TABLE]

in the lower row of Figure 8. The dashed line marks $\ln{A_{2}}=-1.5$, which is used for determining the adaptive cutoff region, within which no particles are inserted in our models.

Model P2 has a constant feeding scalelength $h_{R}=2.5\,{\rm kpc}$, an IC disc scalelength $h_{R,{\rm disc}}=2.5\,{\rm kpc}$ and lives in a $c=9$ halo. Despite the constant input scalelength the final profile is very different from a simple exponential. At $R<5\,{\rm kpc}$, the profile is shaped by the bar, which at $t_{f}$ has a length of $\sim 5\,{\rm kpc}$, similar to that of the MW bar. The bar steepens the profile in the centre and flattens it at radii similar to those of the bar tips. At $R=5-10\,{\rm kpc}$ the surface density profile is mildly flatter than the dashed $h_{R}(R_{0})=2.65\,{\rm kpc}$ line, whereas at $R>10\,{\rm kpc}$ the profile is steeper. The Snhd surface density is close to the upper limit of the observed range $=40-60\,M_{\odot}\,{\rm pc}^{-2}$.

Model T5$\zeta$* has a more compact and more massive IC thick disc than P2, inside-out formation in the range $1.5-3.5\,{\rm kpc}$, a mildly higher final mass and a higher GMC mass fraction. Its thick disc also has a higher-than-average ratio of radial to vertical velocity dispersions $\sigma_{R}^{2}/\sigma_{z}^{2}=1.8$. The outcome is a model with a weaker bar and thus a $\Sigma_{\rm b}(R)$ profile that is well-fit by an exponential at $R=3-15\,{\rm kpc}$. $\Sigma_{\rm b}(R_{0})$ agrees well with the Snhd constraints. This is one of the models which comes closest to the inferred local profile of the MW.

Model M$\beta$1s5 has declining $\sigma_{0}$. It lives in a $c=9$ halo and has inside-out growth in the range $1.5-4.5\,{\rm kpc}$, which generates a final exponential mildly flatter than $h_{R}(R_{0})=2.65\,{\rm kpc}$. Its value of $\Sigma_{\rm b}(R_{0})$ is in agreement with the Snhd constraints. It has a rather weak bar, which influences the profile only at the inner radii.

Model U1 lives in a $c=7.5$ halo and has a massive and compact thick IC disc and inside-out formation in the range $1.5-4.5\,{\rm kpc}$. It has a stronger and longer bar compared to the two previous models with $c=9$ haloes due to a higher baryon fraction as discussed in Paper 1. Its profile shows a steep bar region out to $R\sim 5\,{\rm kpc}$, a flat region at $R\sim 5-8\,{\rm kpc}$ and a shallow exponential decline at $R>8\,{\rm kpc}$. Like for the previous models, $\Sigma_{\rm b}(R_{0})$ agrees well with Snhd constraints.

The declining $\sigma_{0}$ model V$\alpha$9s8$\lambda\zeta$* has a $c=6.5$ halo, a high final mass $M_{\rm f}=6\times 10^{10}\,M_{\odot}$ and grows inside out in the range $1.0-3.5\,{\rm kpc}$. Its $\Sigma_{\rm b}(R_{0})$ is too high for MW constraints and although its bar is weaker than in U1 and also P2, its $\Sigma_{\rm b}(R)$ profile is significantly flattened at $R\sim 4-7\,{\rm kpc}$. At $R\sim 7-12\,{\rm kpc}$ the profile agrees well with an exponential with $h_{R}(R_{0})=2.65\,{\rm kpc}$.

There is a clear tendency of $A_{2}$ being lower in models with declining $\sigma_{0}$ compared to models with thick IC discs. We will discuss this further in Section 5.

4.3 Radial profile evolution of mono-age components

Paper 1 and Section 4.2 have shown that the output scalelengths are somewhat independent of the input scalelengths, as bars and spirals redistribute matter. We know from observations of stars in the Snhd that the old thick populations are more compact than the young thin ones (Bovy et al., 2012). To test, how different populations of stars are affected by changes in the radial distribution, in Figure 9 we plot local scalelengths $h_{R}(R_{0})$ as a function of time for populations of different ages. The pink line is for IC stars and the other colours are for 10 equally spaced age bins of all stars inserted during the simulations. $h_{R}(R_{0})$ is determined by a single-exponential fit to the surface density profile $\Sigma(R)$ at $R=6-10\,{\rm kpc}$, irrespective of how good the fit is.

Model P1s6 has little bar activity at any stage of its evolution and is thus well suited to understand the plots. As it has inside-out growth from $1.5$ to $4.3\,{\rm kpc}$ paired with a thick IC disc with $h_{R,{\rm disc}}=2.5\,{\rm kpc}$, different age components are rather well separated in size. The IC component has a shorter scalelength than in the setup, because it is compressed by the disc’s gravitational field. Apart from a mild shrinking of all populations due to compression and a mild level of noise, which is likely caused by spiral activity, the output scalelengths are essentially set by the input scalelengths.

P2 has a constant input scalelength $h_{R}=2.5\,{\rm kpc}$ and the same ICs as P1s6 and thus the populations of different ages are only mildly separated in size due to continuous compression. At $t\sim 8\,{\rm Gyr}$ bar formation causes an increase in $h_{R}(R_{0})$, which is stronger for younger populations, so the measured scalelengths increase to $h_{R}(R_{0})\sim 2.5\,{\rm kpc}$ for the oldest and $\sim 4\,{\rm kpc}$ for the youngest components. As already shown in Figure 8, U1 is more strongly affected by a bar. U1 has inside-out growth as in P1s6 and a shorter IC disc scalelength $h_{R,{\rm disc}}=2.0\,{\rm kpc}$. Bar formation at $t\sim 4\,{\rm Gyr}$ causes a strong increase in $h_{R}(R_{0})$ for all age groups and at $t\sim 9\,{\rm Gyr}$ bar growth causes another increase, so at $t=t_{f}$, the oldest population has $h_{R}(R_{0})\sim 2.5\,{\rm kpc}$ and the youngest populations have an essentially flat profile.

The declining $\sigma_{0}$ models M$\alpha$1 and V$\alpha$9s8$\lambda\zeta$* have low-mass elliptical ICs, for which we find $h_{R}(R_{0})\sim 2.5\,{\rm kpc}$ fits at early times. They both grow inside out, M$\alpha$1 from $1.5$ to $4.3\,{\rm kpc}$ and V$\alpha$9s8$\lambda\zeta$* from $1.0$ to $3.5\,{\rm kpc}$. As the initial mass of the baryonic ICs is much lower than in P and U models, the amount of compression for the oldest components is stronger. This is especially true for V$\alpha$9s8$\lambda\zeta$, which has more compact feeding scalelengths at early times and also a shorter-than-average SFR time-scale $t_{\rm SFR}=6\,{\rm Gyr}$. The latter increases the mass in stars added at early times and thus also the mass in GMCs present at these formation stages. M$\alpha$1 has a bar from $t\sim 6\,{\rm Gyr}$, which causes a mild increase for all $h_{R}(R_{0})$, whereas bar activity is measurable from $t\sim 3\,{\rm Gyr}$ onwards in V$\alpha$9s8$\lambda\zeta$.

Irrespective of how long a bar is and how strongly it affects disc evolution and how high is the level of compression, the final ordering of $h_{R}(R_{0})$ always reflects the ordering of scalelengths at input. Additionally, in combination with the results of Section 5, it is clear that all models with a bar similar to that of the MW show an increase for $h_{R}(R_{0})$ of all age components with time due to bar formation and growth.

4.4 Circular speed curves

Recently, various surveys of bulge/bar stars and microlensing data have enabled more detailed mass models of the centre of the MW (Wegg et al., 2016; Portail et al., 2017a; Cole & Binney, 2017). These models agree in the following points: 1) The centre of the MW is baryon dominated; 2) The baryonic contribution to the rotation curve at $R\sim 3\,{\rm kpc}$ is $v_{\rm circ,b}\sim 185\,{\rm km\,s^{-1}}$ [although Wegg et al., 2016 find an uncertainty $\sim\pm 25\,{\rm km\,s^{-1}}$]; 3) The contributions of DM and baryons to $v_{\rm circ}$ are roughly equal at $R_{0}$. Further constraints on $v_{\rm circ}(R)$ come from the motion of stars in the Snhd: Schönrich (2012) finds $v_{\rm circ}(R_{0})=238\pm 9\,{\rm km\,s^{-1}}$ and $R_{0}=8.3\pm 0.3\,{\rm kpc}$.

Aumer & Schönrich (2015) (hereafter AS15) presented an inside-out growing model in a $c=9$ halo. Its rotation curve fulfilled the Schönrich (2012) constraints, but has too few baryons in the centre to match any of the constraints from microlensing. In Figure 10 we present circular speed curves for a selection of our models: black is total $v_{\rm circ}(R)$, red is the baryonic contribution $v_{\rm circ,b}(R)$ and blue is the DM contribution $v_{\rm circ,DM}(R)$. The pink boxes mark the Schönrich (2012) constraints. $v_{\rm circ}(R)$ is measured in the midplane of the disc and averaged azimuthally.

Figure 10 shows two $c=9$ models: P2 and M$\beta$1s5. P2 has $h_{R}=2.5\,{\rm kpc}$ both in the ICs and at all times through the simulation. Its value of $v_{\rm circ}(R_{0})$ is at the upper end of allowed values. Its central baryonic contribution is higher than that of AS15 but still too low for the microlensing constraints. M$\beta$1s5 grows inside out in the range $1.5-4.3\,{\rm kpc}$ and has a higher added mass $M_{\rm add}$. The former leads to a weaker central baryonic contribution than in P2 and the latter causes a stronger compression of the halo and thus an unacceptably high $v_{\rm circ}(R_{0})$.

The microlensing constraints suggest shifting mass from the halo to the discs. As was shown above, the Snhd DM density allows concentrations as low as $c=6.5$ given a constant IC halo mass of $M_{\rm tot}=10^{12}\,M_{\odot}$. Models U1 and X5$\zeta$* have haloes with $c=7.5$ and V$\alpha$9s8$\lambda\zeta$* has a $c=6.5$ halo. U1 has a massive thick IC disc with $M_{\rm disc,i}=2.5\times 10^{10}\,M_{\odot}$ and $h_{R,{\rm disc}}=2.0\,{\rm kpc}$. It grows inside out from $1.5$ to $4.3\,{\rm kpc}$ reaching a final mass of $M_{\rm f}=6\times 10^{10}\,M_{\odot}$. Its baryonic contribution to the rotation curve peaks at $v_{\rm circ,b}\sim 190\,{\rm km\,s^{-1}}$ and falls below the DM contribution at $R\sim 8\,{\rm kpc}$ and thus fulfils all microlensing constraints. Its Snhd $v_{\rm circ}(R_{0})$ is lower than that found by Schönrich (2012).

Model X5$\zeta$* has a thick IC disc with $M_{\rm disc,i}=2.0\times 10^{10}\,M_{\odot}$ and $h_{R,{\rm disc}}=2.0\,{\rm kpc}$ and in addition an IC bulge with $a_{\rm bulge}=0.7\,{\rm kpc}$ and $M_{\rm bulge,i}=0.5\times 10^{10}\,M_{\odot}$. It has inside-out growth in the range $1.5-3.5\,{\rm kpc}$ and the same final mass as U1. Consequently, its peak $v_{\rm circ,b}$ is higher at $\sim 200\,{\rm km\,s^{-1}}$, but it still fulfils all constraints from microlensing, as well as the Snhd $v_{\rm circ}(R_{0})$ constraints. The declining $\sigma_{0}$ model V$\alpha$9s8$\lambda\zeta$* also fulfils all constraints but in a different way. It starts from a low-mass elliptical IC and grows a disc with $M_{\rm f}=6\times 10^{10}\,M_{\odot}$ like those of two previous models. Its inside-out growth is from $1.0$ to $3.5\,{\rm kpc}$ and its final $v_{\rm circ,b}(R)$ is rather constant at $170-180\,{\rm km\,s^{-1}}$ in the range $R=2-10\,{\rm kpc}$. Unlike X5$\zeta$, which has a flat total circular speed with $240\,{\rm km\,s^{-1}}$ for $R=3-10\,{\rm kpc}$, V$\alpha$9s8$\lambda\zeta$ has $v_{\rm circ}(R)$ increasing in this radial range from $220\,{\rm km\,s^{-1}}$ to $250\,{\rm km\,s^{-1}}$.

5 Bar formation and evolution

Paper 1 showed that models starting with a thin-disc IC and having no GMC heating undergo strong bar activity from early times. GMC heating can delay and weaken bar formation and evolution and in extreme cases prevent the formation of a strong bar over cosmological time-scales. In this section, we examine how this picture is modified by an old thick-disc component.

5.1 Bar strengths

Figure 11 displays for 10 models the evolution of the $m=2$ Fourier amplitude $A_{2}$ (see Equation 13) for all the stars within $R=3\,{\rm kpc}$. Model Y2 represents thin-disc-only models and we see that $\ln(A_{2})$ instantaneously increases to $-2.5$ as the addition of mass to the thin and compact IC disc makes the system develop non-axisymmetries. As discussed in Paper 2, the radial heating by GMCs at low disc mass and high SFR is important and delays the formation of a strong bar during the first $3\,{\rm Gyr}$ of evolution in Y2. This effect is enhanced in model Y1$\zeta$-, which has $\zeta=0.04$ and thus twice as many GMCs per unit mass of added stars: $\ln(A_{2})$ is kept at $\sim-3$ until $t\sim 7\,{\rm Gyr}$, when it increases to $\sim-1.7$, indicating a rather weak bar.

The curve for model P2 is very different. P2 shares with Y2, the constant radial growth history $h_{R}(t)=2.5\,{\rm kpc}$, the final mass $M_{\rm f}=5\times 10^{10}\,M_{\odot}$ and the shape of the SFH. Its IC disc is thicker, more extended and more massive than the one in Y2 and the normalization of its SFH is thus lower. The existence of a hot disc, which is stable against bar formation prevents the growth of $A_{2}$, although stars are continuously added on cold orbits throughout the simulation. Only at $t\sim 7\,{\rm Gyr}$ has enough thin disc been accumulated to make the composite system unstable to bar formation.

Paper 1 showed that, in the absence of a thick disc, lowering the halo concentration from $c=9$ is problematic, as the system becomes more self-gravitating and bar unstable. U1 has a more massive and more compact IC thick disc than P2, lives in a $c=7.5$ halo and has an inside-out growth history. Despite the higher thick-disc central surface density and the lower halo concentration than P2, the thick disc still suppresses bar formation for $4\,{\rm Gyr}$.

Another ingredient here is the radial-to-vertical velocity dispersion ratio in the IC thick disc. Model T5$\zeta$* has a more massive and more compact IC thick disc and a slightly higher final mass than P2 and inside-out formation from $1.5$ to $3.5\,{\rm kpc}$. Despite the higher surface densities, $\sigma_{R}^{2}/\sigma_{z}^{2}=1.8$ for the ICs in T5$\zeta$* compared to a value of 1 in P2 makes the system more stable against bar formation and thus weakens bar formation more strongly than in P2.

To understand bar evolution in models with declining $\sigma_{0}$, which start with low-mass elliptical ICs, we first examine model K2 that grows a thin disc inside a higher mass elliptical IC. K2 and P2 differ only in that P2 has a thick-disc IC of the same mass. During the early evolution phases of K2, $A_{2}$ is suppressed due to the elliptical ICs. However, in K2 $A_{2}$ increases at an earlier time than in P2 and from $t\sim 3\,{\rm Gyr}$ on shows a significant bar. By construction, the in-plane velocity dispersions of the IC stars at $R\la 2\,{\rm kpc}$ are similar in these models. However, the surface densities are higher and thus the rotation velocities of the IC stars are faster in K2. Consequently, K2 is more unstable to bar formation than P2.

Compared to the K ICs, the M ICs contain an elliptical, which is three times less massive. So a cold disc model in M would have high $A_{2}$ at an earlier time than K2. The hot input dispersions for the old populations in M$\alpha$1, however, act in the same way as the thick IC disc in the P models and delay bar formation until $t\sim 6\,{\rm Gyr}$. We find that the specific shape of declining $\sigma_{0}(t)$ does not significantly influence bar formation.

V$\alpha$1 is the equivalent model to M$\alpha$1, but it lives in a lower concentration $c=6.5$ halo. Still, bar formation is delayed until $t\sim 4\,{\rm Gyr}$. V$\alpha$9s8$\lambda\zeta$* has a higher mass, a more compact feeding history and a shorter SFR time-scale $t_{\rm SFR}$ than V$\alpha$1. All three factors lead to a much faster increase in surface density at early times, which outweighs the fact that at feeding the radial-to-vertical input dispersion ratio $\sigma_{R}/\sigma_{z}=1.25$ is higher than in V$\alpha$1. Consequently, $A_{2}$ is higher at early times.

The specifics of the SFH, the radial growth history, the dispersions of the old components and GMC heating thus determine the bar formation history of an individual model. Figures 8 and 11 however show that halo concentrations $c=6-7$ allow models with reasonable final bars in the presence of hot disc components. The fact that in the lower row of Figure 8, the $m=2$ amplitudes at final times and radii $R>5\,{\rm kpc}$ are lower in models with declining $\sigma_{0}$ than in the thick IC disc models is connected to the gradients of $\sigma_{R}$. In models with thick-disc ICs, the oldest stars have radially constant scaleheights and $\sigma_{R}/\sigma_{z}$ and thus declining $\sigma_{R}(R)$, whereas models with declining $\sigma_{0}$ assume a radially constant input dispersion $\sigma_{R}$. Thus in the end, the outer thick components are radially hotter in models with declining $\sigma_{0}$ and thus less unstable to $m=2$ modes.

5.2 Bar morphology

The central region of our Galaxy is dominated by a bar, the inner part of which consists of a boxy/peanut-shaped bulge at $R<2\,{\rm kpc}$ with an X-shape at $|z|>500\,{\rm pc}$ (Wegg & Gerhard, 2013) surrounded by a vertically thin part, the long bar, extending to $R\sim 4-5\,{\rm kpc}$ (Wegg et al., 2015). AS15 demonstrated that an inside-out growing model without GMCs in a $c=9$ halo produces a bar with X-shaped structure with the tips of this structure at $(x,z)\sim(2,1.3)\,{\rm kpc}$, very similar to the structure of the MW bulge/bar region inferred by Wegg & Gerhard (2013). Paper 1 showed that some thin-disc-only models with GMCs also displayed bars very similar to the one in the MW, but also noted that not all of these bars are buckled and that lower concentration haloes favoured unrealistically long bars.

Here we test how well our thick-disc models can reproduce the MW bar. This is interesting as the chemically defined thick disc is concentrated and should thus have a high mass fraction in the bar region, but can only form a bar if the thin-disc fraction is high enough, as discussed above. Figure 12 shows edge-on and face-on surface density maps of several galaxies at $t_{\rm f}$.

Like the standard Y models in Paper 1 and the model in AS15, P2 lives in $c=9$ halo and shows a bar that is $\sim 5\,{\rm kpc}$ long. Model P1s6, which has a more radially extended feeding history shows only a small, weak bar in the central $R<2\,{\rm kpc}$. The vertical structure of P1s6 is indistinguishable from a pure disc galaxy, whereas P2 has a boxy shape with a lateral extent of $\pm\sim 2\,{\rm kpc}$ and a vertical extent of $\pm\sim 1\,{\rm kpc}$ and a mild hint of an X structure. As in the MW bar, the outer regions are thinner. T5$\zeta$* has a more compact and more massive, but radially hotter thick disc and does not show a bar, just mildly elliptical surface density contours.

The galaxy in model O2 evolves in a $c=9$ halo from an elliptical IC. At $t_{\rm f}$, its bar is $\sim 4\,{\rm kpc}$ long and is currently buckling as indicated by the broken mirror symmetry relative to the $x$-axis as first observed in a simulation by Raha et al. (1991). This event will eventually produce an edge-on peanut bulge with a characteristic X-shape (see Combes & Sanders, 1981). U1 has a compact, massive IC thick disc, a higher-than-average final mass and lives in a $c=7.5$ halo. Its bar grows to a length $\sim 6\,{\rm kpc}$ and is thus only mildly longer than the P2 bar. However, bar buckling has created a significantly more extended X-shape with a lateral extent of $\pm\sim 4\,{\rm kpc}$ and a vertical extent of $\pm\sim 2\,{\rm kpc}$.

The lower row displays bars in models with declining $\sigma_{0}$. Model C$\beta$2s5 grows a disc around a compact bulge IC in a $c=9$ halo. At $t_{\rm f}$, it clearly shows an X-shaped edge-on structure in the central $\sim 2\,{\rm kpc}$, stronger than the similarly sized one in P2. However, its face-on image reveals an almost axisymmetric image. Going back in time, we find that at $t=8.5\,{\rm Gyr}$ C$\beta$2s5 exhibited a strong bar that had not yet buckled. From the evolution of $A_{2}(R<3\,{\rm kpc})$ depicted in Figure 11, we learn that the bar formed around $t=6\,{\rm Gyr}$. Between $t=8.5\,{\rm Gyr}$ and $t_{\rm f}$ C$\beta$2s5 undergoes buckling, but in the final $\sim 2\,{\rm Gyr}$ of the simulation its bar becomes continuously weaker.

Models M$\alpha$1$\zeta$* ($c=9$ halo) and V$\beta$8s5 and V$\alpha$9s8$\lambda\zeta$* ($c=6.5$ halo) are models with declining $\sigma_{0}$ starting from low-mass elliptical ICs. They display bars of reasonable sizes ($3.5-5.5\,{\rm kpc}$) and varying strengths. Their vertical profiles all show boxy edge-on shapes with vertical extents $\pm\sim 1.0-1.5\,{\rm kpc}$ and lateral extents that are wider for longer bars.

In summary, it is possible to create models that show reasonable agreement with the vertical profile of the MW, its circular speed curve and at the same time contain a bar, which in length, strength and vertical extent agrees reasonably with that of the MW. The problem is that the details depend on mass and size growth history, DM halo density and GMC heating and that, additionally, the evolution of bar length and strength is to some degree stochastic (Sellwood & Debattista, 2009), so that, not even within the limits of our methods, it is possible to determine which model best represents our Galaxy.

5.3 Bar age structure

In Figure 13, we examine the age structure of the bar. For several models we plot a map of median age $\tau_{\rm med}$ in the $x$-$z$ plane, where $x$ is measured along the major axis of the bar and $z$ is perpendicular to the disc. All stars with $|y|<1.5\,{\rm kpc}$ are considered for the map. If $\tau_{\rm med}$ at a certain position corresponds to a star particle from the ICs, we apply a distinct yellow colour, whereas inserted star ages range from black (young) to orange (old). We overplot in white edge-on density contours.

We start by analysing model T5$\zeta$*, which does not show a noticeable bar, in the left-hand panel of the middle row. Clearly, above $|z|\sim 600\,{\rm pc}$ the IC thick disc dominates at all $x$. The thin disc is youngest in the plane and has a vertical age gradient, as an effect of stars being born cold and being vertically heated by GMCs. It also has a radial age gradient, due to inside-out formation.

In contrast to model T5$\zeta$*, the three models displayed in the the top row of Figure 13 at $t_{\rm f}$ show X-shaped edge-on structures. These three models all grow thin discs within thick IC discs. The X-shape is beautifully visible in the $\tau_{\rm med}$ maps. The cones of relatively lower surface density above and below the galactic centres are dominated by old, thick-disc stars, whereas younger stars fill diamonds that lie at both sides of the cones in the $x$-direction. These diamond-shaped regions have vertical age gradients of different strengths. This is connected to the fact that we do not feed new stars into regions of strong bars. The bar in U1 forms already at $t=4\,{\rm Gyr}$, whereas the bar in P2 only forms at $t=7.5\,{\rm Gyr}$. As shown in AS15, stars can be captured by the bar, but the rates of this process are low, so the age structure is only mildly affected. Consequently, the number of relatively young stars present in P2 is high and vertical and radial age gradients are present in the bar region. In U1, the region at $|x|<4\,{\rm kpc}$, which is part of the diamonds, is well mixed in $\tau_{\rm med}$. The bar in W1 forms at $t=6\,{\rm Gyr}$ and the model is thus intermediate to P2 and U1 in terms of age gradients.

The vertical mixing of stars of various ages happens during bar buckling, as is illustrated by model O2, which is undergoing the process at the depicted moment. We see how the younger stars are spread vertically at around $|x|\sim 2\,{\rm kpc}$ and the model thus transitions from a disc-like to an X-shaped age structure. Note that, as bars can vary in length and strength with time, cold, young stars can be added in the plane after buckling events and there can be multiple buckling events, which complicates the age structure.

The remaining four models in Figure 13 feature declining $\sigma_{0}$. As the thick discs of these models comprise fed-in stars, which have radial and vertical density distributions that vary continuously with age, the age structure is more complicated to interpret, as in real galaxies. At $t_{\rm f}$, V$\alpha$1, like T5$\zeta$*, has no X-shaped edge-on structure. The thin discs of these galaxies exhibit qualitatively similar vertical and radial age gradients. The thick disc of V$\alpha$1 shows a continuation of these gradients at higher $|z|$; the gradients in these regions are, however, much shallower.

M$\alpha$1 and V$\alpha$9s8$\lambda\zeta$* have edge-on peanut-shaped density structures similar to that of P2 and, indeed, the age structures of the three bars are qualitatively similar. As the thick disc background, however, is, again, much less distinct, the structure combining old cones plus younger diamonds is much less evident and might be hard to measure in a real galaxy.

Model C$\beta$2s5 differs from the other models with declining $\sigma_{0}$ in having compact bulge ICs and thus the density of IC stars near the galactic centre is higher and the $\tau_{\rm med}$ gradients away from the centre are steeper than in models with low-mass elliptical ICs. A relatively recent and strong buckling event before the disappearance of the bar and a radial growth history with constant $h_{R}=2.5\,{\rm kpc}$ has led to relatively young areas at $|x|\sim 2\,{\rm kpc}$.

6 Discussion

6.1 Vertical profiles and age gradients

In all final models the vertical profile of the Snhd is fitted well by the sum of two exponentials similar to that of the MW. When the IC contains a thick disc, the vertical profile is double-exponential from the outset, with a growing thin-to-thick density ratio, whereas when $\sigma_{0}$ declines, the vertical profile gradually evolves a double-exponential structure as a sufficiently massive, cold and thin population forms, and the scaleheights of both components change continuously.

We cannot directly measure the evolution of scaleheights over cosmic time, but in nearby galaxies we can probe this evolution through observations of the radial variation of the vertical profile and age structure of the disc. The radial variation of the vertical profile is dominated by the evolution of the thick disc because our thin discs have scaleheights that are almost independent of both $R$ and $t$. Our thick-disc ICs are set up with declining $\sigma_{z}(R)$ and radially constant scaleheights, and the latter property is roughly conserved to the present epoch, as radial migration causes only mild levels of flaring. Models with declining $\sigma_{0}$ have, by construction, radially constant $\sigma_{z}$ at birth. As the oldest and hottest stars are barely affected by vertical heating, their thick discs flare strongly and their scaleheights increase with radius. In these models, radial migration weakens flaring, but the effect is not strong enough to balance the outward increase in scaleheight imprinted at birth.

This has important consequences on the age structure of the disc. At radii $R=5-12\,{\rm kpc}$ and altitudes $|z|>1\,{\rm kpc}$ models with thick-disc ICs are dominated by old IC disc stars; they thus show no radial variation in median age $\tau_{\rm med}$. At lower altitudes these discs become younger with increasing $R$. The age structure is markedly different in models with declining $\sigma_{0}$. On account of the strongly flaring old and intermediate-age components, they show negative $\tau_{\rm med}$ gradients at all $|z|$. Models with declining $\sigma_{0}$ thus agree better with recent measurements of the radial age structure at various altitudes in the MW by Martig et al. (2016b). They also show somewhat flatter vertical age gradients at $R_{0}$, in rough agreement with measurements by Casagrande et al. (2016). For thick IC disc models, strong inside-out growth improves the agreement with these measurements.

Our discs evolve in isolation, whereas at least the low-density outskirts of discs are likely to be affected by processes capable of significant vertical thickening, such as disc-satellite interactions (Kazantzidis et al., 2008) or infall of gas with misaligned angular momentum (Jiang & Binney, 1999). Including minor mergers or adding stars in tilted outer discs would thus be a valuable extension of our models.

6.2 The interplay between the dark halo, the thick disc and the bar

The evolution of a disc depends on the local density of the dark halo because increasing the latter reduces the extent to which the disc controls the gravitational field in which it moves. In particular, decreasing the local DM density increases the amplitude of non-axisymmetric structure. We varied the local DM density by varying the initial concentration parameter $c$ at a fixed halo mass, $M_{\rm tot}=10^{12}\,M_{\odot}$. Concentrations in the range $c=6-9$ work well. Indeed, after the DM has been compressed by the disc, the final DM density in the Snhd then agrees with observational constraints, and in many models a bar similar to that of the MW emerges before the current epoch. Although non-axisymmetric structures in the disc transfer angular momentum to the DM, the dark halo does not acquire a core like that favoured in the MW (Cole & Binney, 2017). If such a core exists in a $\Lambda$CDM context, it thus probably formed in the very early evolution stages of the Galaxy that are not modelled here.

Observational constraints on the circular speed, $v_{\rm circ}(R_{0})$, near the Sun (Schönrich, 2012) and on the microlensing optical depth towards the MW bar/bulge (Wegg et al., 2016; Cole & Binney, 2017) indicate a baryon-dominated central MW and roughly equal contributions of DM and baryons to $v_{\rm circ}(R_{0})$. To achieve this in our models, initial concentrations $c\sim 6-7$ are favoured as higher values do not allow for enough baryonic mass in the central regions.

Paper 1 favoured $c\sim 9$ because thin-disc-only models with lower values of $c$ showed unrealistically long and strong bars. The presence of old thick, and thus kinematically hot, disc components alleviates this problem. By shifting mass from the cold thin disc to a thick, radially hot component, the formation of a bar can be delayed by several Gyr. It is immaterial whether the hot component is included in the ICs or arises from declining $\sigma_{0}$, and the ratio of radial to vertical dispersions can be $\sigma_{R}/\sigma_{z}\sim 1$, so smaller than the values $\sigma_{R}/\sigma_{z}\sim 2$ in thin discs. As a consequence, the time available for the bar to grow in strength and length is limited. Indeed, models with $c=9$ haloes tend to show an unrealistically weak final bar. Models with $c\sim 6-7$ generally have bars with lengths similar to that of the MW’s bar. The delaying of bar formation by old thick discs also explains the observation that the fraction of barred disc galaxies decreases with increasing redshift (Sheth et al., 2008).

6.3 The edge-on structure of bars

Several models have bars that are morphologically similar to that of the MW: a boxy bulge in the central $R\la 2\,{\rm kpc}$ with an X-shape up to $|z|\sim 1-1.5\,{\rm kpc}$ at the centre of a thinner outer bar that extends to $R\sim 5\,{\rm kpc}$ (Wegg & Gerhard, 2013; Wegg et al., 2015). In one of these models the bar dissolved at the very end of the simulation. However, notwithstanding the face-on surface density being almost axisymmetric, the edge-on peanut shape survived. Thus not all observed boxy edge-on bulges need be bars.

In edge-on density projections of these models, the X-shapes are not as striking as in the model of AS15, which lacks a thick disc and GMCs, but features an isothermal gas component. This is not surprising as the MW bulge does not show an X-shape in all stellar components. Dékány et al. (2013) found that old and metal-poor RR Lyrae stars appear to have a more spheroidal shape and Portail et al. (2017b) found that low-metallicity stars ($\left[{\rm Fe/H}\right]<-0.5$) in the bar/bulge contain a much lower fraction of stars on bar-supporting orbits than stars with higher metallicities. Moreover, the X is possibly absent in younger populations as well (López-Corredoira, 2016). Analysis on the edge-on age structure of our bars offers insight into why the shape should vary with stellar population.

The X structures are particularly pronounced in models with thick-disc ICs that have final boxy bulge/bar regions that are more extended than in the MW. Here bar buckling spreads stars from the thin disc vertically, and because they mainly populate the 2:2:1 resonant orbit family (Pfenniger & Friedli, 1991), these stars form a structure that, seen edge-on, resembles two diamonds overlapping at the galactic centre. The cones above and below the centre, which are not populated by the 2:2:1 orbits, are dominated by thick-disc stars. Young stars that were captured by the bar after the buckling event will be found in the plane. Consequently, a distinct age pattern should be observed in buckled edge-on bars.

A similar separation in edge-on morphology between the oldest, the intermediate-age and the youngest stars in a barred galaxy has also been found in the simulations of Athanassoula et al. (2016). Debattista et al. (2016) recently studied bar formation in galaxies that contain disc populations with differing random motions. They showed that radially cooler populations form stronger bars, the edge-on profiles of which are vertically thinner and peanut-shaped, whereas the hotter populations form a weaker bar with a vertically thicker edge-on box shape (see also Fragkoudi et al., 2017).

The thick and thin discs of models with declining $\sigma_{0}$ are not strictly separated in age, as they are in models with thick-disc ICs, and in consequence their characteristic age structure is less clear in an edge-on age map. Thus while observations of the edge-on age structures of bars have the potential to betray the formation history of the bulge, central thick and thin discs and the timing of the bar buckling event, the constraints will be less tight if the scenario with declining $\sigma_{0}$ is more appropriate than that in which the thick disc is included in the ICs.

6.4 Radial redistribution and inside-out growth

Bars and spiral structure make it hard to steer the disc’s radial scalelength $h_{R}$ at the solar radius $R_{0}$ to a preferred value, $h_{R}(R_{0})\sim 2.6\,{\rm kpc}$. As was already discussed in Paper 1, making the disc more compact results in stronger non-axisymmetries, which in turn leads to more mass redistribution and larger $h_{R}(R_{0})$. We have shown that this is the case for all age components and for all models with appropriate bars, so in the past $h_{R}(R_{0})$ would likely have been smaller than it is today. Assigning higher values of $\sigma_{R}/\sigma_{z}$ to the old thick disc can yield steeper-than-average profiles, but doing so weakens bars inappropriately.

Still, for models that grow inside out, $h_{R}(R_{0})$ always increases with decreasing age, just as observations of the Snhd suggest (Bovy et al., 2012). As these observations show scalelengths $h_{R}(R_{0})\sim 1.5-2.0\,{\rm kpc}$ for the most compact and oldest mono-abundance populations (see also Cheng et al., 2012), $h_{R}=2\,{\rm kpc}$ can be regarded as an upper limit on the birth scalelength at the earliest times, but model V$\alpha$9s8$\lambda\zeta$* demonstrates that the scalelength at birth could have been as small as $h_{R}=1\,{\rm kpc}$, a conclusion similar to that of Schönrich & McMillan (2017). The flat age-metallicity relation of the Snhd and the radial metallicity gradient in the MW make it hard to infer the scalelength of current star formation by studying mono-abundance populations. In our models input scalelengths at late times in the range $h_{R}\sim 3-4\,{\rm kpc}$ give reasonable results.

6.5 Thick-disc formation scenarios

The age structure of the MW disc points towards a model with declining $\sigma_{0}$. Most hydrodynamical cosmological simulations of disc galaxies support a picture in which birth dispersions and gas fractions decline continuously with time (e.g. Bird et al., 2013; Stinson et al., 2013; Ma et al., 2017). Such a scenario has also been inferred from observations of H$\alpha$ kinematics (e.g. Kassin et al., 2012; Wisnioski et al., 2015). However, Di Teodoro et al. (2016) argue that there is no substantial difference between the gas kinematics of galaxies at redshifts $z_{\rm rs}\sim 1$ and today. Moreover, Martig et al. (2014) find that models that have kinematics in line with those found in the Snhd favour a two-phase formation scenario, in which the thick-disc stars are born in a turbulent, merger-dominated phase and the thin-disc stars are born cold and heated subsequently.

Such a two-phase scenario motivated our models with thick-disc ICs, but in setting up an equilibrium thick stellar disc with a radially constant scaleheight we have ignored correlations between stellar ages, kinematics and density profiles that would naturally arise during formation of the proto-thick disc. Declining $\sigma_{0}$ inherently produces such correlations. However, an inappropriate thick disc is still liable to emerge through poor choices for $\sigma_{0}(t)$ or $h_{R}(t)$, or the choice of a radially constant $\sigma_{0}$.

The main limiting factor of our models is thus the lack of a self-consistent heating mechanism for thick-disc stars: in both scenarios thick-disc stars are created ad hoc. The heating mechanism will affect non-axisymmetries and the radial distribution of matter, which are crucial for disc evolution. As measurements of gas fractions at redshifts $z_{\rm rs}\sim 2$, a time consistent with the formation of the chemically defined thick disc of the MW, indicate that molecular gas makes up 50 per cent or more of the baryonic masses of galaxies (e.g. Genzel et al., 2015), the lack of a realistic gas component is a connected problem. Although at early times models such as V$\alpha$9s8$\lambda\zeta$* have as much as 45 per cent of their baryonic mass in GMCs, the GMCs have the same mass function as a present-day spiral galaxy, in which the gas fraction is lower and stars form cold. In a picture in which turbulence driven by gravitational disc instabilities causes stars to form with large dispersions (Forbes et al., 2012), molecular complexes would be expected to be more massive. Scattering of stars by massive clumps could contribute to thick-disc formation (Bournaud et al., 2009).

7 Conclusions

We have presented a new set of idealized $N$-body simulations of disc galaxies with both thin and thick discs within live dark haloes. These models are grown over $10-12\,{\rm Gyr}$ by continuously adding new stellar particles with specified age-dependent velocity dispersions. Short-lived massive particles represent GMCs. Thin discs grow by the addition of stars on near-circular orbits, whereas for thick-disc components we rely on two different concepts: a) create an appropriate thick disc in the ICs and only add thin-disc stars during the simulation, or b) start with low-mass, diffuse elliptical or compact bulge ICs and add stars with continuously declining input velocity dispersion $\sigma_{0}(t)$. Hence in scenario b) we form kinematically hot thick-disc stars at early times and cold thin populations at late times, whereas in scenario a) the simulation starts after the structure that will morph into the thick disc is fully formed. To understand the evolution of our models, we simulate a variety of histories of star formation, dark halo densities and thick disc properties.

Both types of models can produce at final time $t_{\rm f}$ models that are similar in structure to the MW. We find:

• •

Both scenarios create double-exponential vertical profiles. The scaleheight of the thin disc is governed by GMC heating. To achieve a MW-like thick-disc exponential scaleheight $h_{\rm thick}\sim 1\,{\rm kpc}$ at $t_{\rm f}$, thick-disc ICs with isothermal vertical scaleheights $z_{0}\sim 1.7\,{\rm kpc}$ are suitable. For declining $\sigma_{0}$ models, the input velocity dispersions should be $\sigma_{0}\sim 40-50\,{\rm km\,s^{-1}}$ at the earliest formation stages. Thick-disc scaleheights are not affected by GMC heating.

• •

Models need to undergo inside-out growth to reproduce the observed dependence of radial scalelength $h_{R}$ on chemical composition of disc stars. We find that models that grow from $h_{R}\sim 1-2\,{\rm kpc}$ at early times to $h_{R}\sim 3-4\,{\rm kpc}$ today are suitable.

• •

To explain the baryon dominance of the Galactic Centre, the circular speed curve of the MW and the structure of the bar in the presence of a thick disc, we favour DM haloes that at mass $M_{\rm{DM}}=10^{12}\,M_{\odot}$ have an initial concentration parameter $c\sim 7$.

It is essential that thick-disc stars are already hot when the thin disc starts forming because Paper 1 showed that heating by GMCs and non-axisymmetries is incapable of producing the thick disc, although it explains the properties of the thin disc. The presence of the thick disc modifies the evolution of the thin disc, but the final properties of the thin discs in our thin+thick disc models are similar to those of thin-disc-only models in slightly more concentrated dark haloes. Crucially, this change in halo density and the presence of a hot and thick disc make it possible to bring models with appropriate bars into agreement with the baryon fractions inferred for the central MW.

Regarding the non-axisymmetric structures of the disc, we find:

• •

Bars with a structure similar to that of the MW bar, i.e. a boxy/peanut-shaped bulge at $R<2\,{\rm kpc}$ with an X-shape surrounded by a vertically thin part extending to $R\sim 4-5\,{\rm kpc}$, can be found in some of our viable models. Stochasticity in the evolution of bar lengths and strengths complicates the comparison.

• •

The presence of a hot, thick-disc stellar population at the start of thin-disc formation suppresses non-axisymmetries and delays the formation of the bar.

• •

In models with an appropriate bar, the local exponential scalelengths $h_{R}(R_{0})$ of all mono-age populations are increased by the radial redistribution of matter that the bar and spirals generate. Populations measured in the Snhd today have thus likely had lower $h_{R}(R_{0})$ in the past.

• •

The dark halo’s density profile is modified by the growth of the disc and the non-axisymmetric structures that form in the disc, but the profile does not develop a central core as is currently favoured for the centre of the MW.

To distinguish between formation scenarios, it is helpful to study the radial and vertical age structure of disc galaxies. We find:

• •

Our two types of models for the creation of the thick disc differ significantly in the predicted age maps of the discs. The observed Snhd radial age gradient at $|z|>1\,{\rm kpc}$ and the vertical age gradient both favour models with declining $\sigma_{0}$. However, the measurements are still rather uncertain and models with thick IC discs by construction ignore any internal structure of the thick disc. Such structure could be added to these models.

• •

Bar buckling in thin+thick disc systems creates characteristic age patterns in edge-on views of the bar region. Buckling predominantly affects thin-disc stars and causes them to form a structure in the $(R,z)$ plane resembling two diamonds overlapping at the galactic centre. The cones above and below the centre are dominated by thick-disc stars.

Considering the wealth of data on the structure of the MW that will soon become available from surveys such as Gaia (Gaia Collaboration, 2016), evolutionary models of disc galaxies that grow over cosmological time-scales and contain both thick and thin discs will be essential to connect the data to the formation history of the MW. Our models allow for a relatively controlled and flexible setup, can be produced in large numbers and capture a wealth of important dynamical processes.

We have demonstrated that our models can reasonably reproduce a variety of observations of the structure of the MW. No model sticks out as particularly similar to the MW in all aspects, but this is to be expected given the remaining shortcomings in modelling. In a companion paper (Paper 4) we examine the models presented here in light of Snhd kinematics and constraints on radial migration.

Acknowledgements

We thank the referee for comments that helped improve the paper. It is a pleasure to thank Ralph Schönrich for valuable discussions and comments on the manuscript.

This work was supported by the UK Science and Technology Facilities Council (STFC) through grant ST/K00106X/1 and by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 321067. This work used the following compute clusters of the STFC DiRAC HPC Facility (www.dirac.ac.uk): i) The COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology. This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. ii) The DiRAC Complexity system, operated by the University of Leicester IT Services. This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. iii) The Oxford University Berg Cluster jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Oxford. DiRAC is part of the National E-Infrastructure.