© ESO, 2012

1. Introduction

Most local galaxies are observed to host a thick disc component (e.g., Comerón et al. 2011), including our Milky Way (MW), as first proposed by Gilmore & Reid (1983). Compared to the thin disc, the thick disc is composed of older, metal-poor, alpha enriched stars. An important question of modern galactic astronomy is how thick discs came into existence.

One possibility is that they were born thick at high redshift from the internal gravitational instabilities in gas-rich, turbulent, clumpy discs (Bournaud et al. 2009) or in the turbulent phase associated with numerous gas-rich mergers (Brook et al. 2004, 2005). They could also have been created through accretion of galaxy satellites (Meza et al. 2005; Abadi et al. 2003), where thick disc stars then have an extragalactic origin.

Another possibility is that thick discs are created through the heating of preexisting thin discs. This can either happen fast, on a Gyr timescale, in an early violent epoch, or as a continuous process throughout the galaxy lifetime. In the first case, the thick disc would appear as a clear distinct component in chemistry and phase-space, while in the latter case it would rather be seen as a gradual, continuous transition. Observationally, whether the MW thick disc is consistent with the former or the latter is heavily debated (e.g, Schönrich & Binney 2009; Bovy et al. 2012).

Preexisting thin discs can be heated fast through the effect of multiple minor mergers (Quinn et al. 1993; Villalobos & Helmi 2008; Di Matteo et al. 2011), the rate of which decreases with decreasing redshift. Evidence for merger encounters can be seen in structure in phase-space of MW disc stars (e.g., Minchev et al. 2009; Gómez et al. 2012b; Purcell et al. 2011), which could last for as long as 4 Gyr (Gómez et al. 2012a). The recent work by Minchev et al. (2012c) presented numerical results indicating that the velocity dispersion of the oldest MW-disc stars cannot be obtained without substantial merger activity at high redshift.

The fact that galactic discs heat with time has now been known for over 60 years, established from both observational and theoretical work. To explain the observed correlation between the ages and velocity dispersions of solar neighborhood stars, Spitzer & Schwarzschild (1951, 1953) suggested that massive gas clouds (then undetected) scattered stars from initially circular, into more eccentric and inclined, orbits. Giant molecular clouds were thought to be the sole scattering agents (e.g., Mihalas & Binney 1981) until Lacey (1984) showed that the observed ratio of the dispersion in the direction perpendicular to the Galactic plane and that toward the Galactic center, σ_z/σ_R, was too low to be consistent with the predictions from this scattering process. This resulted in the development of models that incorporated the heating of the stellar disc from transient spiral structure (Barbanis & Woltjer 1967; Sellwood & Carlberg 1984) in addition to scattering from molecular clouds (Jenkins & Binney 1990; Jenkins 1992). Other proposed models for the heating of stars include scattering by halo black holes (Lacey & Ostriker 1985) or dark clusters (Carr & Lacey 1987), scattering by giant molecular clouds and halo black holes (Hänninen & Flynn 2002, 2004), and vertical disc heating through “popping” star clusters (Kroupa 2002, 2005) or through “levitation" caused by the 2:2 resonance between the horizontal and vertical oscillations of stars in an inside-out forming disc, converting their radial actions into vertical actions (Sridhar & Touma 1996b,a). An investigation of the origin of the age-velocity relation in the MW was presented recently by House et al. (2011), by comparing with cosmological simulations. Interestingly, a recent work has shown that the age-velocity relation could arise not because of external perturbers or the gradual heating of a thin disc, but because the velocity dispersion of newly born stars decreases with redshift (Forbes et al. 2012). Finally, a viable internal mechanism for disc heating is the interaction of multiple spiral density waves, as proposed by Minchev & Quillen (2006) or bar and spirals (Minchev & Famaey 2010). Multiple waves are now known to exist in both observation (Elmegreen et al. 1992; Rix & Zaritsky 1995; Meidt et al. 2009) and simulations (e.g., Masset & Tagger 1997; Quillen et al. 2011; Minchev et al. 2012b; Grand et al. 2012; Roškar et al. 2012), therefore this heating process may be ubiquitous in disc galaxies.

Interestingly, non-axisymmetric patterns also give rise to another important process regarding the secular evolution of galaxy discs: radial migration. This process is related to the redistribution of angular momentum of stars within the disc, generally resulting in an increase in disc scale-length. Several radial migration mechanisms have been described in the literature: (i) the effect of the corotation resonance of transient spiral density waves (Sellwood & Binney 2002; Roškar et al. 2012), (ii) the effect of the non-linear coupling between multiple spiral waves (Minchev & Quillen 2006) or bar and spirals (Minchev & Famaey 2010; Minchev et al. 2011; Brunetti et al. 2011), and (iii) perturbations due to minor mergers (Quillen et al. 2009; Bird et al. 2012). Note that in case (ii) above, spirals do not need to be short-lived transients for efficient migration; relatively long-lived waves are expected to exist in barred discs (Binney & Tremaine 2008; Quillen et al. 2011; Minchev et al. 2012b), but are also found in non-barred simulations (e.g. D’Onghia et al. 2012). Recently, Comparetta & Quillen (2012) showed that, even if patterns are long-lived, radial migration can result from short-lived density peaks arising from interference among density waves overlapping in radius. This “disc mixing” can give rise to a number of phenomena, such as flattening in radial metallicity gradients (e.g., Minchev et al. 2011; Pilkington et al. 2012) and extended stellar density profiles (e.g., Roškar et al. 2008b; Sánchez-Blázquez et al. 2009; Minchev et al. 2012b).

Recently there has been a growing conviction that radial migration can result in a thick disc formation by bringing out high-velocity-dispersion stellar populations from the inner disc and the bulge. Such a scenario was used, for example, in the analytical model of Schönrich & Binney (2009), where the authors managed to explain the Milky Way thick- and thin-disc characteristics (both chemical and kinematical) without the need of mergers or any discrete heating processes. Similarly, the increase of disc thickness with time found in the simulation by Roškar et al. (2008a) has been attributed to migration in the works by Sales et al. (2009) and Loebman et al. (2011). We have to point out, however, that how exactly radial migration affects disc thickening in dynamical models has never been demonstrated.

In this paper we would like to examine in detail the effect of radial migration on the increase of radial and vertical velocity dispersion and, thus, on its possible relation to the formation of thick discs. For this purpose we study three isolated galaxy simulations and three simulations in the cosmological context. We must note that, while the simulations we examine address some of the most widely used initial conditions for simulating galactic discs: preassembled N-body discs and discs growing in the cosmological context, in this paper we are not testing the validity of other models of formation, such as the one used by Sales et al. (2009) and Loebman et al. (2011).

Fig. 1 First column: m=2 Fourier amplitudes, A2, as a function of time (color bar) for our model galaxies. Columns 2–3: edge-on view for each disc and bulge (separately) for earlier and later times of their evolution, as indicated in each panel. Columns 4–5: scale-heights, h, at  ~2.5 disc scale-lengths. The dotted red line shows a single exponential fit. Note that in 3–4 Gyr of evolution with strong perturbation from bars and spirals, h only doubles despite the efficient stellar radial migration.

2. Description of the simulations

2.1. Tree-SPH N-body simulations

For the majority of this paper we study three main runs of isolated disc galaxies from the GalMer database (Di Matteo et al. 2007): the giant S0, Sa, and Sbc runs (hereafter gS0, gSa, and gSb). In GalMer, for each galaxy type, the initial halo and the optional bulge are modeled as Plummer spheres, with characteristic masses M_H and M_B, and characteristic radii r_H and r_B, respectively. Their densities are given by ${\mathit{\rho }}_{\mathrm{H}\mathit{,}\mathrm{B}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\left(\frac{\mathrm{3}{\mathit{M}}_{\mathrm{H}\mathit{,}\mathrm{B}}}{\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{H}\mathit{,}\mathrm{B}}^{\mathrm{3}}}\right){\left(\mathrm{1}\mathrm{+}\frac{{\mathit{r}}^{\mathrm{2}}}{{\mathit{r}}_{\mathrm{H}\mathit{,}\mathrm{B}}^{\mathrm{2}}}\right)}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{·}$(1)On the other hand, the initial gaseous and stellar discs follow Miyamoto-Nagai density profiles with masses M_g and M_s, and vertical and radial scale-lengths given by b_g and a_g, and b_s and a_s, respectively: $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{g}\mathit{,}\mathrm{s}}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}& \mathrm{=}& \left(\frac{{\mathit{b}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{g}\mathit{,}\mathrm{s}}}{\mathrm{4}\mathit{\pi }}\right)\\ \mathrm{\times }& & \frac{{\mathit{a}}_{\mathrm{g}\mathit{,}\mathrm{s}}{\mathit{R}}^{\mathrm{2}}\mathrm{+}\left({\mathit{a}}_{\mathrm{g}\mathit{,}\mathrm{s}}\mathrm{+}\mathrm{3}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}}\right){\left({\mathit{a}}_{\mathrm{g}\mathit{,}\mathrm{s}}\mathrm{+}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}}\right)}^{\mathrm{2}}}{{\left[{\mathit{a}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}\mathrm{+}{\left({\mathit{a}}_{\mathrm{g}\mathit{,}\mathrm{s}}\mathrm{+}\sqrt{{\mathit{z}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right]}^{\mathrm{5}\mathit{/}\mathrm{2}}{\mathrm{(}{\mathit{z}}^{\mathrm{2}}\mathrm{+}{{\mathit{b}}_{\mathrm{g}\mathit{,}\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$(2)The parameters for the initial conditions of the three runs, including the number of particles used for the gaseous disc, N_g, the stellar disc and bulge, N_s, and the dark matter halo, N_DM, are given in Table 1. The initial Toomre parameter of both stars and gas is taken to be Q = 1.2 as the initial condition of the Tree-SPH simulations and particle velocities are initialized with the method described in Hernquist (1993). Full details of the simulations are given in Di Matteo et al. (2007) and Chilingarian et al. (2010).

These galaxies have been studied extensively by Minchev et al. (2011, 2012b). In these Tree-SPH N-body simulations the gSb model develops a bar similar to the one seen in the Milky Way (MW) (Minchev et al. 2007, 2010), while both the gS0 and gSa bars are substantially larger than that (see Fig. 1). We have shown (Minchev et al. 2011) that this has strong effects on the migration efficiency found in these discs, with larger bars giving rise to stronger mixing (also consistent with studies of diffusion coefficients, e.g., Brunetti et al. 2011; Shevchenko 2011). Most recently, we (Minchev et al. 2012b) demonstrated that, in these same models, the strongest changes of angular momentum occur near the bar’s corotation (CR) region, in contradiction to the expectation that patterns need to be transient for efficient migration (Sellwood & Binney 2002).

Fig. 2 First row: time evolution of the radial vertical dispersion profiles, σz(r), for the gSb, gSa, and gS0 models. Different colors and line styles correspond to the times indicated in the first row, rightmost panel. The bars’ corotation (CR) and outer 2:1 Lindblad resonance (OLR) are shown by the dotted and solid vertical lines, respectively, with colors corresponding to the different times. Second row: same but for the radial velocity dispersion, σr(r). Third row: same but for the ratio σz/σr. The slope in σz/σr outside the bar shows that the radial profile flattens more with time than the vertical one, while its overall decrease means that the disc heats more radially than vertically.

2.2. Cosmological re-simulations

We concentrate on the controlled models just described for most of the analyses we are about to do, since these discs are easier to analyze given their constrained nature, lack of gas accretion and mergers. However, we also will use much more realistic simulations, following the self-consistent assembly of galactic discs in the cosmological context, in order to test the validity of our results. These numerical experiments were first presented by Martig et al. (2012), where the authors studied the evolution of 33 simulated galaxies from z = 5 to z = 0 using the zoom-in technique described in Martig et al. (2009). This technique consists of extracting merger and accretion histories (and geometry) for a given halo in a Λ-CDM cosmological simulation, and then re-simulating these histories at much higher resolution (150 pc spatial, and 10⁴⁻⁵   M_⊙ mass resolution), replacing each halo by a realistic galaxy containing gas, stars and dark matter. We have chosen to examine three simulations, which we refer to as the C1, C2, and C3 models, with approximately flat rotation curves and circular velocities V_c ~ 200, 180, and 200 km s^-1, respectively. The bulge/bar/disc decomposition and disc mass growth for the three galaxies are presented in the third row of Fig. 25 (C1), the first row of Fig. 29 (C2), and the bottom row of Fig. 27 (C3) by Martig et al. (2012).

These simulations differ in several ways from the gSb, gSa, and gS0 models described above: (i) the discs grow self-consistently as the result of cosmological gas accretion from filaments and gas-rich mergers, as well as merger debris; (ii) the simulation techniques are different both in the way gravitational forces are computed (Particle Mesh vs hierarchical tree) and in the way gas dynamics is treated (sticky-particles vs. SPH); and (iii) the spiral and bar instabilities are present for a cosmological time, thus the effects we see are not the result of recent bar formation, for example.

Fig. 3 First row: initial (dashed black), final (solid blue), and net (dash-dotted red) vertical velocity dispersions, σz,i, σz,f, and Δσz, respectively, as functions of the initial radius, r0, for particles ending up in the green bin at the final time indicated in each panel. Two time periods for each of the three models are shown, as indicated. Migrating stellar samples decrease or increase their velocity dispersion is such a way as to match the value at the destination. Second row: same as the first row, but for the radial velocity dispersion. Third row: the corresponding initial radial distributions of stars ending up in the green bin at the final time indicated in each panel. The bar’s CR and OLR, measured at the final times, are shown by the vertical dotted and solid black lines. Most migrators come from the inner disc, with a strong contribution from particles originating near the bar’s CR.

3. Disc thickening and velocity dispersions

Up until Sect. 8 we only analyze our constrained, Tree-SPH N-body simulations described in Sect. 2.1. As we have shown in Minchev et al. (2012b), all of these models exhibit multiple patterns, giving rise to strong radial migration. To show the strength of the discs’ asymmetric components, in the first column of Fig. 1 we plot the m = 2 Fourier amplitudes, A₂/A₀, as a function of radius for each galaxy, where A₀ is the axisymmetric component and m is the multiplicity of the pattern. Different color curves show the time evolution during the first Gyr. A₂ indicate the bi-symmetric structure in the disc. For example, for the gS0 model the bar is identifiable by the smooth curve in the inner disc, which at later times has a maximum at  ~3 kpc and drops almost to zero at  ~8 kpc. Deviations from zero seen beyond that radius are due to the spiral structure. Note that the spirals are strongest for the gSa model, related to the initial 10% gas fraction in its disc. The gSb model is the one with parameters closest to the MW: bar size  ~3−4 kpc, v_c ~ 220 km s^-1, and a small bulge. The circular velocities and bar pattern speeds for all three models can be found in Fig. 2 by Minchev et al. (2012b).

To see how much discs thicken, in Cols. 2, 3 of Fig. 1 we plot the edge-on view for each disc and bulge (separately) for earlier (0.25 Gyr) and later (4 Gyr for gSb and gS0, and 3 Gyr for gSa) times of their evolution, as indicated in each panel. Significant thickening is seen for all models, especially for gSa and gS0 due to their larger bars. The bulges are also found to expand, not significantly so for the gSb. Although from these plots it is evident that the discs thicken, it is important to see how this is reflected in their scale-heights. These are shown in Cols. 4, 5 at a radius 12 kpc (~2.5 disc scale-lengths), which, when scaled, is comparable to the solar distance from the Galactic center. The dotted red line shows single exponential fits.

Surprisingly, we find that, despite the strong radial migration, the initial scale-heights, h ≈ 0.3 kpc, hardly double in the 3 − 4 Gyr of evolution for gSa and gS0, respectively, and increase by only  ~50% in 4 Gyr for gSb.

The disc thickening we see in Fig. 1 must be related to an increase in the vertical velocity dispersion. To see how our discs heat, in the first two rows of Fig. 2 we show the time evolution of the vertical, σ_z(r), and radial, σ_r(r), velocity dispersion profiles for the stellar populations of the gSb, gSa, and gS0 models in the first, second and third columns, respectively. Initial bulge components are not considered. In each panel the different color curves show different times, as indicated in the top left panel. The dotted and solid vertical lines indicate the radial positions of the bar’s corotation (CR) and 2:1 outer Lindblad resonance (OLR) with colors matching the different times.

An interesting observation is that the σ_r profiles flatten at later times, especially for the gS0 and gSa models, i.e., the galaxies with the larger bars and stronger migration histories. This is reflected in the ratio σ_z/σ_r, shown in the bottom row, where the later curves (blue and black) develop a negative slope. Moreover, σ_r increases more in the outer parts of the disc, compared to the region near the bar’s CR. In fact, we observe radial cooling near the CR for both the gS0 and gSa models: the solid-black curve falls below the dashed-green and dash-dotted blue ones. The reason for this is revealed in the next Sections.

4. Cooling/heating of stellar samples migrating outwards/inwards

Stellar samples in the inner galactic discs possess high velocity dispersions in all three components. It has been expected that by migrating radially outwards, stars retain their random energies and thus thicken the disc. If this were true, then how do we explain the lack of sufficient disc thickening even for the strong large-scale disc mixing seen in Fig. 1? To answer this question, we consider the samples of stars used to estimate the scale-heights in Fig. 1. In the first row of Fig. 3 we plot the initial (dashed black), final (solid blue), and net (dot-dashed red) vertical velocity dispersions, σ_z,i, σ_z,f, and Δσ_z, respectively, as functions of the initial radius, r₀, for particles ending up in the green bin (11 < r < 13 kpc) at the end of the given time interval. This is done for each galaxy in the time intervals 0.5 < t < 1.5 and 0.5 < t < 3 Gyr (0.5 < t < 4 Gyr for gS0 and gSb), as indicated in the figure. The second row shows similar plots but for the radial velocity dispersion. The distribution of initial radii of stars ending up in the green bins at the corresponding times are shown in the third row.

From the first two rows it is immediately apparent that stellar samples arriving from the inner discs decrease their velocity dispersions by as much as 50% when coming from the smallest radii, with the effect diminishing the smaller the radial distance traveled. Conversely, stars coming from radii exterior to the green bin increase their velocity dispersion. Therefore, on the average, newcomers change their random energies in such a way as to approximately match the non-migrating population at the radius at which they arrive. Thus, even though at later times (second, fourth and sixth columns) a larger fraction of migrators enters the annulus under consideration (as seen in the bottom row), this has no effect on the velocity dispersion increase.

Fig. 4 First row: contour density plots of the final versus initial radii for all particles in different time spans for the gSb or gSa models, as indicated. The orange dashed line shows the locus of non-migrating, cold particles. Nonuniform deviations from this line indicate that there exist preferential radii where migration efficiency is stronger, however, these are not clearly seen. Second row: non-migrating stars, extracted from the distributions above as described in the text. Third row: the migrators only, obtained by subtracting the distributions in the second from those in the first row. The dotted-red and solid-blue vertical/horizontal lines in all panels indicate the positions of the bar’s CR and OLR, respectively, at the initial/final times. Strong clumps can now be seen across the bar’s CR for example. Fourth row: number density plots of the changes in angular momentum, ΔL, versus the initial angular momentum (or equivalently initial radius), L, divided by the rotational velocity at each radius. The broken red and blue vertical strips indicate the location of the bar’s CR and OLR, where their widths corresponding to the change during each time period. The pink lines connect some corresponding groups of migrated stars. The bar’s CR is active throughout the entire simulations.

The effect we just described is more dramatic for the radial velocity dispersion: for the gS0 (dissipationless) model the non-migrating stars (blue line portion inside the green bin) are hotter at the final time than the migrators, as evident by the positive slope of the solid-blue curve. Since the majority of stars arrive from the inner disc, the final annulus is being cooled, rather than heated by the migrators! We can see that, for all three models, stars arriving from inside the bar’s CR are radially cooler than the non-migrators: note the positive slope in the blue curve inside the CR region (i.e., to the left of the dotted vertical line). Because for gS0 the radial bin considered happens to be between the CR and the OLR, we see the most cooling for this model. We find similar results when the in is placed at a similar location with respect to the bar’s resonances in the case of the gSa and gSb models (see Fig. 6).

The velocity distribution of migrating stars before and after they have migrated are smooth and appear approximately Gaussian, as can be seen in Fig. A.1 in the Appendix.

5. Separating migrators from non-migrators

To assess the contribution of migrating stars to the intrinsic velocity dispersion of the disc we devise a simple procedure for separating migrators and non-migrators.

• (1)

  We consider a disc annulus of a certain width at some initial andthen final time.

• (2)

  We separate all stars into “non-migrators” (those that are found in the selected radial bin at both the initial and final times) and “migrators” (those that were not present in the bin initially but are there at the final time).

• (3)

  We do this for radial bins over the entire extent of the disc, overlapping every 0.3 kpc, for example, in order to get information at every radius.

In step (1) above we need to adopt a certain width for each annulus considered. We chose this as follows. Epicycle sizes are on the order of a ~ σ_r/κ, where κ is the epicyclic frequency. In order to contain non-migrating stars, we allow for particles to oscillate around their guiding radii by  ± a, therefore, we select our bin size to be 2a. Since both σ_r and κ are functions of radius, so is the value of a. In Fig. A.2 we plot all three functions versus galactic radius. Note that the bin sizes increase with radius, where the effect is stronger for the hotter gSa and gS0.

We apply the above procedure to our gSb and gSa models and display the results in Fig. 4. The first row shows contour density plots of the final versus initial radii for all particles in different time spans for the gSb or gSa models, as indicated. We exclude a small portion of the inner disc (white rectangle in the lower left corner) in order to display better the outer contours. The orange dashed line shows the locus of non-migrating, cold particles. Nonuniform deviations from this line indicate that there exist preferential radii where migration efficiency is stronger, however, these are not clearly seen. We note that the particle distributions presented here are not the true density, but a combination from all migrators and non-migrators resulting from the overlap of the radial bins we consider. However, deviation from the true density are small, as shown in Fig. A.3 and discussed in the Appendix.

The second row of Fig. 4 presents the non-migrating stars extracted from the distributions shown in the first row as outlined in (1), (2), and (3) above. Deviation from the orange line are about 1 − 2 kpc, increasing outwards, in agreement with the bins-size variation shown in Fig. A.2.

The third row of Fig. 4 shows the migrators only, obtained by subtracting the distributions shown in the second row from those in the first row. The dotted-red and solid-blue vertical and horizontal lines indicate the location of the bar’s CR and OLR at the initial (horizontal) and final (vertical) times. We can now see clearly overdensities on each side of the orange dashed line. We show how our results change by considering different bin size values in the appendix (see Fig. A.5).

In the fourth row of Fig. 4 we plot number density contours of the changes in the specific angular momentum, ΔL, versus the initial specific angular momentum, L, estimated during each time periods. Both axes are divided by the rotational velocity at each radius, therefore L is approximately equal to the initial radius and ΔL gives the distance by which guiding radii change. The broken-red and solid-blue vertical strips indicate the location of the bar’s CR and OLR, where their widths corresponding to the change during each time period related to the bar’s slowing down. For all time periods considered we find a nice correlation with structure across the bar’s CR, where stars inside the CR migrate outward and those outside it migrate inward. The pink lines connect some corresponding groups of migrated stars in the r_final − r_initial and the L − ΔL planes. At the earlier times shown (0.5 < t < 1.5 Gyr, left column) the changes in angular momentum are the strongest due to the potent spirals and their interaction with the bar. At the intermediate times, as the spirals weaken, the typical feature of negative slope across the bar’s CR is much better seen, in addition to some outer structure related to the spirals. Finally, at the later times (3 < t < 4 Gyr, third column) the bar’s CR is clearly dominating. In general, although here we only show the entire simulation time for the gSb model and only the beginning of gSa, clumps across the bars’ CR are invariably seen for all time periods and all three models. This continuous bar activity was shown also by Minchev et al. (2012b) (see their Fig. 7).

Contours are denser above the orange line in the r_final − r_initial plane. This is most extreme during the early time period shown for the gSa model (rightmost column of Fig. 4). The reason for this behavior is that, in the outer disc, stars are preferentially shifted outwards due to the decreasing stellar density. The effect is much stronger for gSa because the outward transfer of angular momentum for this simulation is much stronger than for the gSb disc, which Minchev et al. (2012b) related to the very efficient non-linear coupling among the bar and spiral waves of different multiplicity. Note that, since the total angular momentum of the disc must be conserved (with a small contribution to the halo and bulge), the outward transfer of angular moment needs to be balanced. This balance is achieved by transferring large amounts of mass inward of the bar’s CR (with small negative changes in stellar guiding radii), in contrast to the exponentially decreasing density in the outer disc (but much larger increase in guiding radii). An indication of this process may be found in the strong clump of inward migrators just outside the bars’ CR, seen in all plots but the latest time period for gSb; the reason for the latter is due to the saturation of migration efficiency in the outer disc in the absence of gas accretion (and thus spiral strength).

Fig. 5 Migrators’ contribution to the disc heating. First row: radial profiles of the vertical velocity dispersion, σz(r), for non-migrators (dashed black curve), outward going migrators (solid brown curve), and inward migrators (dash-dotted green curve) for the gSb model, estimated in the time periods indicated in the top row of each column. Second row: σz(r) estimated for non-migrators (dashed black curve), all migrators (i.e., inward and outward going, solid green curve), and the total population (dashed-dotted blue curve). Third and Fourth rows: same as the first and second rows, respectively, but for the radial velocity dispersion profile σr(r). An unexpected result is that outward migrators cool the disc inside the bar’s OLR (and inward ones heat it), unlike in the case of σz. Some increase in σz and σr from migrators is seen beyond the OLR, while the inner disc is cooled both radially and vertically. Figure 6 quantifies better the migrators’ contribution to σz and σr.

6. Migrators’ contribution to the disc velocity dispersion

By separating stars in the entire disc into “migrators” and “non-migrators”, we are now in a position to investigate their individual effects on the disc velocity dispersion increase. We consider the time intervals 0.5 < t < 1.5 and 1.5 < t < 3 Gyr, corresponding to the first two columns in Fig. 4. Because the heating in both radial and vertical directions mostly saturates around 1−2 Gyr (see Fig. 2), the effect at later times is similar to what we find during the period 1.5 < t < 3 Gyr.

In the first row of Fig. 5 we show the radial profile of the vertical velocity dispersion, σ_z(r) for non-migrators (dashed black curve), outward going migrators (solid brown curve), and inward migrators (dash-dotted green curve) for the gSb model. The dotted-red and solid-blue vertical lines denote the bar’s CR and OLR, respectively, at the final time of the time period considered.

We find that inward migrators cool the disc and the outward ones heat it. This is consistent with the information contained in Fig. 3, where we found that particles migrating to a final annulus at r = 12 kpc arrive with σ_z values slightly lower/higher than the non-migrating population (evident by the blue line’s slightly negative slope in Fig. 3).

Fig. 6 Fractional changes in σz(r) and σr(r) resulting from the migrated stars. The dashed black curve plots the quantity $\mathrm{\Delta }{\mathit{\sigma }}_{\mathit{z}}\mathrm{=}\mathrm{(}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}\mathrm{-}{\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}\mathrm{)}\mathit{/}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}$, where σz,all and ${\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}$ are the vertical velocity dispersions obtained from the total population and the non-migrators, respectively. A similar expression estimates Δσr, shown by the solid orange curve. The dotted-red and solid-blue vertical lines denote the bar’s CR and OLR, corresponding to the final time shown in each time period. Note the inverse correlation between Δσz and Δσr, indicating coupling between the radial and vertical motions of stars.

From the plots we just described (first row of Fig. 5) it is not clear how inward and outward migrators contribute to the velocity dispersion profile of the total population since we did not normalize to the number of particles in each group. How much σ_z(r) of the total population is affected depends on the fraction of each kind of migrators at each radial bin. To demonstrate the true contribution of migrators, in the second row of Fig. 5 we plot the vertical velocity dispersion for non-migrators (dashed black curve), all migrators (i.e., inward and outward going, solid green curve), and the total population (dashed-dotted blue curve). Interestingly, we find that the velocity dispersion from both types of migrators is very similar to that of the non-migrating particles, except for the inner parts of the disc for the earlier times shown. It is surprising that not only do migrators not heat most of the disc, but they even cool it inside the CR. Some increase in σ_z can be seen outside the OLR.

We now examine the effect of migrators on the radial velocity dispersion profile, σ_r(r). The third and fourth rows of Fig. 5 are similar to the first two, except the plotted quantity is σ_r(r). A distinct difference compared to the vertical velocity dispersion is the oscillatory behavior of σ_r(r) found in all three curves in the third row. Comparing the third and fourth rows for the earlier time period, it becomes clear that the dominant effect inside the bar’s OLR comes from the outward migrators, which here cool, rather than heat the disc as in the case of σ_z(r). Just inside and outside the OLR, the effect of migrators is to reduce the wiggles seen in the non-migrators’ σ_r, consequently resulting in a smooth radial velocity dispersion profile.

To better quantify the contribution of migrators, we now calculate the fractional changes in σ_z(r) and σ_r(r) resulting from the migrated stars for the same time periods shown in Figs. 5. For the vertical velocity dispersion, we estimate these as $\mathrm{\Delta }{\mathit{\sigma }}_{\mathit{z}}\mathrm{=}\mathrm{(}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}\mathrm{-}{\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}\mathrm{)}\mathit{/}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}\mathit{,}$(3)where σ_z,all and ${\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}$ are the vertical velocity dispersions for the total population and the non-migrators, respectively. An equivalent expression estimates the contribution to the radial velocity dispersion. We plot the results in Fig. 6. The dotted-red and solid-blue vertical lines denote the bar’s CR and OLR, corresponding to the final time considered in each panel. It now becomes much clearer how exactly migrators affect the discs.

For all models shown Δσ_z (dashed black curve) shows oscillatory behavior, changing abruptly from negative inside the CR, to positive outside it. The maximum deviations from zero for gSb, for example, are  ~ 7% cooling in the inner disc and  ~ 4% heating at larger radii, and appear consistent with Fig. 5. Clearly, such contribution to the vertical velocity dispersion would results in disc flaring, as we show in the next Section.

Now we look at the migrators’ contribution to the radial velocity dispersion, Δσ_r (solid orange curve). Remarkably, there appears to be an inverse correlation between Δσ_z(r) and Δσ_r(r), where the two functions cross near the bars’ CR and OLR. This oscillatory behavior with zeros near the CR and the OLR is present for all three models during the first time period shown, which is when the discs heat the most. Note that the crossing point near the OLR occurs above zero as the OLR is shifted more outward in the disc (top to bottom), due to the decreasing number of non-migrators and prevailing influence of outward migrators. Outside the OLR both Δσ_z and Δσ_r are positive, with a contribution of less then 5% (gSb) or less than 10% (gSa). At the later times the general trend of an inverse correlation for the radial and vertical velocity dispersions is still seen, although the contribution to the Δσ_r is larger.

Fig. 7 Scale-height radial profiles for migrators, non-migrators and all stars in the time intervals 0.5 < t < 1.5 Gyr (top) and 1.5 < t < 3 Gyr (bottom). The green solid curves indicate the fraction of migrators as functions of radius during each time period. The scale-height for the non-migrating stars is approximately constant with radius for all models (blue dashed curves, top). Due to the inner disc cooling and outer disc heating in the vertical direction from migrators (dotted-red curve, see also Fig. 6), all discs flare (solid black curve).

The inverse functional behavior just described indicates that the vertical and horizontal motions of migrators are coupled, as expected for a system governed by non-linear dynamics.

7. The effect on the discs’ scale-height

We showed in Fig. 6 that the migrators’ contribution to the vertical velocity dispersion increase is small and mostly in the outer disc. This suggests that our discs flare. To find out if this is really the case, we now measure the scale-height at different radial bins in the gS0, gSa and gSb discs.

In Fig. 7 we show radial profiles of the discs’ scale-height for migrators, non-migrators and all stars in the time intervals 0.5 < t < 1.5 Gyr (top) and 1.5 < t < 3 Gyr (bottom). As in the figures so far, we take the initial time at 0.5 Gyr (or after the bars have formed) in order to avoid the effect of the initial formation of asymmetric structure. The green solid curves indicate the fraction of migrators as functions of radius during each time period. Note that during the later time period (1.5 < t < 3 Gyr) the stars we refer to as “non-migrators” may have experienced some migration previously. This is also true for the earlier times, i.e., some migration takes place before 0.5 Gyr, but the strongest effect occurs in the range 0.5 < t < 1.5 Gyr. All our procedure ensures is that they stayed at the same radii during the later time interval.

Figure 7 shows that the scale-height for the non-migrating stars is approximately constant with radius for all models (blue dashed curves, top) during the earlier time period. In fact, for gS0 the scale-height for non-migrators shows a decrease with radius. However, due to the inner disc cooling and outer disc heating in the vertical direction from migrators (dotted-red curve, see also Fig. 6), all discs flare (solid black curve). This is entirely consistent with Fig. 6, where we found that the contribution of migrators to the vertical velocity dispersion is negative inside the bar’s CR and positive in the outer disc.

By comparing the variation of the migration fraction with radius (green solid curve) to the scale-height of the total population (black solid curve) we note an inverse correlation inside the bars’ CR: the more migrators present at a given radius, the smaller the scale-height. This is especially obvious at the earlier time period and for the more efficient gSa and gS0. At larger radii, where all stars are in fact migrators (although fraction is smaller that 100%, these stars have migrated out before t = 0.5 Gyr), scale-height increases with migration fraction, which is expected.

Fig. 8 Application to three high-resolution cosmological re-simulations. First and second rows: face-on and edge-on density contours for the C1, C2, and C3 models. Third row: change in angular momentum as a function of disc radius: effective migration is seen for all models. Fourth row: vertical velocity dispersion profiles for non-migrators, outward and inward migrators. Fifth row: same but for the radial velocity dispersion. Bottom row: fractional contribution of migrators to the total vertical and radial velocity dispersions (as Fig. 6). Results are very similar to the simpler models presented earlier.

8. Application to cosmological re-simulations

We now apply our newly developed technique for separating migrators and non-migrators to our three, much more realistic simulations in the cosmological context, described in Sect. 2. To assess the effect of the internal perturbers and be able to compare to our previous results, here we only consider the last Gyr of evolution before z = 0, thus, avoiding strong merger activity. We defer a more detailed study of these models to a future work.

During the  ~1 Gyr of evolution we study here, disc penetration by satellites occurs only at (r,t) ≈ (18 kpc, 13.15 Gyr) for C1 and (r,t) ≈ (10 kpc, 13.4 Gyr) for C3. There are no disc encounters within r < 30 kpc for the C2 model. All of the aforementioned companions have mass ratios less than 1:100. Face-on and edge-on stellar density contours are plotted in the first and second rows of Fig. 8, respectively, depicting the time t = 13.4 Gyr. This is roughly in the middle of the timespan considered for the rest of the plots in the figure. Some small satellites are still seen in all three simulations. For C1 and C2 these are orbiting outside 20 kpc radius, but for C3 there is a small companion at r ~ 10 kpc, located close to the plane at this time.

8.1. First results for migration in the cosmological context

To see how much migration has taken place in the last Gyr of evolution, in the third row of Fig. 8, for each model, we plot the changes in the specific angular momentum of stellar particles in the range 4 < r < 23 kpc and  | z |  ≲ 4 kpc. The horizontal axis plots the initial angular momentum, L_i (estimated at t = 12.73 Gyr); this is divided by the circular velocity and, thus, approximately equals galactic radius. The vertical axis plots ΔL = L_f − L_i, where, L_f is the angular momentum estimated as t = 13.7 Gyr. Stars at ΔL ~ 4, for example, gain angular momentum which brings them  ~4 kpc outward in the disc and conversely for negative ΔL. The dotted-red and solid-blue vertical lines show the radial location of each bar’s CR and OLR. The high peak near the OLR for the C3 model is caused by the satellite closest to the galactic center seen in the disc face-on view above.

We note that this is the first time such plots are presented for simulations in the cosmological context. It is interesting that the strongest changes in angular momentum happen near the bar’s CR, just as noted recently by Minchev et al. (2012b) (using the gSb, gSa, and gS0 models), despite the fact that here the bars have been formed for a long time (more than 6 Gyr) and are not transient. The strength of the bar correlates with the effect on the r − ΔL plane; the bar length is related to the pattern speed, as expected, with longer bars rotating slower, since bars exist within their CR radii.

8.2. Contribution of migrators to the increase in velocity dispersion

Now that we have established that migration is efficient in these runs (in fact, very similar to our controlled simulations, see Minchev et al. 2012b, Figs. 7 and 4 in this paper), we can look for the effect of the migrating and non-migrating stars (during the time period considered) on the vertical and radial velocity dispersion profiles, just as we did in Figs. 5 and 6. In the fourth row of Fig. 8 we plot the vertical velocity dispersion of inward (moving in) and outward (moving out) migrators, and the non-migrators, similarly to Fig. 5. Cooling from inward migrators and heating from outward ones is seen, mostly inside the bar’s CR (dotted-red vertical). This is in agreement with what we found for our preassembled models, with the difference that the effect here is weaker in the outer disc, except for the C3 model, where a significant peak is seen between the CR and the OLR; the latter is most likely due to the small satellite mentioned earlier, taking place in the same radial range.

The fifth row of Fig. 8 is the equivalent of the top one, but for the radial velocity dispersion. As in our controlled models, cooling from outward migrators and heating from inward ones (the reverse of what happens for σ_z) occurs between the bar’s CR and its OLR. Even in the case of C3, which is perturbed by the satellite at that radius, we see the same behavior.

Finally, the bottom row of Fig. 8 plots the fractional contribution to the vertical and radial velocity dispersions from all migrators estimated from Eq. (3). Similarly to Fig. 6, migrators contribute to Δσ_z(r) by  ~5%, cooling the disc inside the CR for C1, just as we found for the controlled models. For C2 and C3, the radial variation of Δσ_z is more erratic, yet contribution is mostly less than 3%. The C1 model’s dynamics is dominated by the presence of its large bar, which could explain its closer similarity to the controlled simulations.

While outward migrators for all three models are radially cooler than the non-migrators in the region between the CR and the OLR (fifth row of Fig. 8) just as the case of our constrained simulations, Δσ_r(r) is only slightly below zero for C1 and C2 and slightly positive for C3. These differences are most likely due to tidal effects associated with the orbiting satellites seen in the face-on plots in the first row. An inverse relation between Δσ_z(r) and Δσ_r(r) for our cosmological simulations can still be seen.

It is remarkable that both our isolated discs and the just described much more sophisticated simulations (which also use a completely different simulation technique) exhibit very similar behavior. This makes our result, that migrators do not contribute significantly to the increase of the velocity dispersion of the disc, much stronger. We, therefore, suggest that our finding is a fundamental property of the secular evolution of galactic discs. It is possible that when stronger satellite perturbations are included, as is the case at the early stages of galaxy formation, the situation may be different.

9. Conservation of vertical and radial actions?

Our findings so far suggest the existence of a conserved dynamical quantity for stars moving around the galactic disc, resulting in the conspiracy to approximately match (on the average) the kinematics expected at their destination. Since from Figs. 3, 5, and 6 it is clear that this is not the energy, the obvious candidates are the radial and vertical actions.

9.1. The vertical action

Binney & McMillan (2011) have shown that, to approximate how the vertical motion of disc stars is affected by their radial motion in an axisymmetric potential, the adiabatic approximation is valid, i.e. that the vertical action is adiabatically conserved in the epicyclic sense¹. As first proposed by Minchev et al. (2012a) and recently confirmed by Solway et al. (2012), the vertical actions of migrating stars are closely conserved, as well.

Note that while in an axisymmetric disc the energy and angular momentum are conserved quantities (thus, conservation of actions is expected), since stars migrate due to resonant perturbation associated with disc asymmetries (such as the bar and spirals waves), the above mentioned classical integrals of motion are not necessarily applicable here. Nevertheless, following Minchev et al. (2012a), we now show that even in the case of a simple galactic disc model, assuming conservation of vertical action is in good agreement with our results.

We assume that for stars at a given radius, the vertical action is given by J_z = E_z/ν, where E_z and ν are the vertical energy and frequency, respectively. We also assume that this quantity is conceived (on average) as stellar samples shift guiding radii (i.e., migrate). This is a good approximation in the limit that (i) the vertical motion decouples from the planar motion (a good assumption when z_max ≪ r, which holds for most disc stars, but note that we did find evidence to the opposite in Sect. 6), and (ii) the stars move outwards much more slowly than their vertical oscillations (not really the case). From Gauss’ law and Poisson’s equations, $\mathit{\nu }\mathrm{~}\sqrt{\mathrm{2}\mathit{\pi G}\mathrm{\Sigma }}$, therefore, as stars migrate radially, their vertical frequency changes as ν_mig(r) ~ exp( − r/2r_d), given that the stellar density varies as Σ ~ exp( − r/r_d). To conserve the mean vertical action of a set of stars migrating outwards, we then require that their mean vertical energy, $\mathrm{⟨}{\mathit{E}}_{\mathit{z,}\mathrm{mig}}\mathrm{⟩}\mathrm{~}{\mathit{\sigma }}_{\mathit{z,}\mathrm{mig}}^{\mathrm{2}}$, decreases as ν. Therefore, the vertical velocity dispersion of a stellar sample migrating radially outwards will decrease as σ_z,mig(r) ~ exp( − r/4r_d), where r is the radius at witch it ends up.

On the other hand, the vertical velocity dispersion in the ambient (non-migrating) disc component changes with radius as ${\mathit{\sigma }}_{\mathit{z}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{~}\sqrt{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{h}}$ (e.g., Eq. (7) by van der Kruit and Freeman 2011). Assuming the scale-height, h, remains constant, we find that σ_z(r) ~ exp( − r/2r_d). Thus, while a star’s vertical energy does decay as it migrates outwards, it does not decay as fast as the underlying typical vertical energy of disc stars and, therefore, still increases the vertical temperature of the disc where it arrives. Similarly, stars migrating inwards would cool the disc. This is exactly what we see in Fig. 5 at radii outside the bars’ CR: the disc is heated/cooled by outward/inward migrators. We found that inside the CR the discs are, in fact, cooled vertically, since contribution is made only by inward migrators (see Figs. 5 and 6). Similarly, in the disc outskirts outward migrators dominate, thus the vertical temperature is increased.

9.2. The radial action

For the simple galaxy disc model considered above, the same argument which lead us to conclude that the vertical action is conserved can also be applied in the plane. We assume a conservation of radial epicyclic action J_p = E_p/κ where E_p is the radial epicyclic energy and and κ ~ 1/r is the radial epicyclic frequency, valid for a flat rotation curve. We, therefore, expect that orbits conserving J_p lose energy E_p at the rate of 1/r.

There are some differences when compared to the vertical motion analysis. We found that E_z of migrators needs to change as ν_mig(r) ~ exp(−r/2r_d), which is the same as the expected drop in velocity dispersion. The radial epicyclic frequency, on the other hand, decreases as 1/r. We may, therefore, expect that the changes in E_z and E_p for the same migrating sample are different, with migrators contributing more to the radial disc heating.

At this point a deficiency in the simple model we consider is found. As revealed by Fig. 5, third row, in our simulations outward migrators radially cool the disc inside the bar’s OLR instead of heating it, i.e., the opposite to the behavior of their vertical velocity dispersion. This effect is most prominent in the region between the CR and the OLR for all the models we study, while at r ≳ r_OLRσ_r(r) and σ_z(r) behave similarly.

9.3. The need for a more complex conserved quantity: a function of both J_z and J_p?

We found that the conservation of epicyclic vertical action for migrators is consistent with our numerical results. In the case of the horizontal motion, this simple model is less accurate, as we just discussed. This result is not really surprising, as (i) the adiabatic approximation is less adequate for horizontal motions than vertical ones even in the axisymmetric case ², and (ii) it is also well known that the classical epicycle theory greatly underestimates the radial action even in the case of pure planar motion (see Dehnen 1999 for a better approximations of J_p).

Additionally, the adiabatic invariance of J_p is not necessarily a good approximation for a time evolving system perturbed by multiple asymmetric patterns with stars migrating faster than their epicyclic oscillations. We note that the inverse correlation between Δσ_z(r) and Δσ_r(r) we found in Sect. 6 is suggestive of an exchange between the vertical and radial motions. An exchange between J_z and J_p can be achieved when vertical and horizontal resonances coincide in radius. Such locations may be found in the region inside the bar’s CR and between the CR and the OLR, i.e., where we see the largest differences between Δσ_z and Δσ_r. For example, Combes et al. (1990) have shown that the approximate coincidence in radius of the bar’s vertical and horizontal 2:1 ILRs (inside CR) are responsible for the bar’s peanut shape. We, therefore, speculate that a conserved quantity for migrating stars could be a function of both J_z and J_p, perhaps in a similar fashion to the conservation of the sum J_z + J_p proposed by Sridhar & Touma (1996b,a) in the context of an adiabatically growing galactic disc, which involves a time dependent perturbation as is the case of our discs perturbed by multiple patterns³. More work is needed to find a possibly similar invariant for the more complicated case of migrating stars due to disc asymmetries.

10. Discussion and conclusions

In this work we have used Tree-SPH N-body (controlled) simulations and unconstrained, cosmological re-simulations, to study the effect of internal disc evolution on disc thickening. All analyzed models have been previously found, or here shown, to develop strong multiple patterns, resulting in efficient radial disc mixing.

We devised a procedure for separating migrating from non-migrating stars in simulations. We found an excellent correspondence between the changes in angular momentum of stars as a function of initial radius and structure in the r_initial − r_final space, thus, verifying that our separating technique is correct (see Fig. 4). Using this simple procedure, we were able, for the first time, to study and contrast the kinematics and spatial distribution of migrators and non-migrators in galactic discs.

A number of recent works have assumed that stars migrating from the inner disc retain their high energies, thus populating a thick disc component. To test this expectation, we first compared the disc velocity dispersion profiles in the vertical and radial directions, σ_z(r) and σ_r(r), respectively, for the migrating and non-migrating populations. Surprisingly, we found that the migrators’ contribution to the vertical velocity dispersion, Δσ_z, is negative inside the bar’s corotation (CR) radius, r_cr, and less than 5% in the range r_cr < r < 3h_d, where h_d is the disc’s scale-length (see Fig. 6). This is related to the fact that migrating samples change their energy in such a was as to match closely the kinematical properties of the non-migrating population at the radius of arrival. We thus have shown that, not only do migrators not heat significantly, but they cool a large portion of the disc.

Given the migrators’ contribution to vertical heating in the outer discs and cooling at small radii, disc flaring could be expected. We showed that this is indeed the case, by plotting the disc scale-height as a function of radius for migrators and non-migrators for each simulation. As anticipated, non-migrators exhibit an approximately constant scale-height in contrast to the flared disc composed of migrators (see Fig. 7). We have to note that this flaring is more significant outside the disc’s initial radial distribution (~12 kpc), i.e., where all particles are migrators.

Using a simple galactic disc model, where both the radial, J_p = E_p/κ, and vertical, J_z = E_z/ν, actions, are conserved, we explained the decrease of velocity dispersion of outward migrators as disc cooling. In order to conserve the actions of a migrating sample of stars, their energies, and thus the velocity dispersions, need to change with radius similarly to κ and ν. While stellar samples migrating outwards lose vertical energy, $\mathrm{⟨}{\mathit{E}}_{\mathit{z}}\mathrm{⟩}\mathrm{~}{\mathit{\sigma }}_{\mathit{z}}^{\mathrm{2}}$, their reduced σ_z is still higher than what is expected from extrapolating σ_z ~ Σ(r). This is consistent with Fig. 5, where we found that the outward migrators slightly heat the disc vertically and inward ones cool it. However, instead of a larger contribution from migrators to the radial velocity dispersion, as predicted by this simple model, we found that outward migrators cool the disc instead of heating it inside about three scale-lengths. This means that, in the approximation we considered, the radial energy decreases with radius faster than the epicyclic frequency, i.e., the cooling is even more drastic than predicted by the simple model from Sect. 9.1. We suggested that a better conserved quantity is a function of both the vertical and radial actions.

We estimated that for the gSa and gS0 models, bulge contribution to the thick disc beyond  ~ 2.5 disc scale-lengths is only about 1%. For a MW-like bulge (the gSb model) no stars reach this radius. This is related to the bar’s CR lying mostly outside such a small bulge (see Fig. 1). It should be kept in mind that here the initial distributions of the bulges were prescribed as Plummer spheres. As evident form the third row of Fig. 4, a strong contribution comes from disc stars inside the bars’ CR.

Our controlled simulations did not consider gas accretion. In a more realistic case, the mass increase due to a continuous thin-disc formation would contract the thickened discs and further decrease the scale-heights (Bournaud et al. 2009), particularly in the outer discs, where for the gSb, gSa, and gS0 models all stars are migrators. To check this, we applied our technique to much more realistic, high-resolution, unconstrained simulations in the cosmological context, where discs build up self-consistently. Interestingly, we found very similar results to those obtained with our simpler models, where the migrators’ contribution to the total increase in random motions is typically less then 3%.

We do not expect insufficient resolution to be affecting our conclusions for the reason that low resolution should result in more heating, not less, as we find here. Hence, increasing the (intermediate) resolution of our controlled simulations would give rise to even less disc thickening⁴. Moreover, the recent work by Roškar et al. (2012) has shown that increasing the softening length from 50 to 500 pc does not result in any significant change in the disc dynamics, therefore, we expect that the migration efficiency is not affected. Finally, our results do not change significantly when we use the higher resolution cosmological re-simulations studied in Sect. 8. In any case, we are not interested here so much in the overall heating, but in the difference between the increase of velocity dispersion of migrators and non-migrators.

It is remarkable that both our isolated discs and the much more sophisticated simulations in the cosmological context (which also use a completely different simulation technique) exhibit very similar behavior. We believe the inability of migration to thicken discs is a fundamental property of secular disc evolution, irrespective of the migration mechanism at work. We note, however, that when mergers are involved (i.e., when the evolution becomes non-secular), the approximate conservation of the vertical and radial actions we find here, may not be valid anymore. We defer a study of this important topic to a future work.

10.1. Expected signatures of migration in galactic discs

Considering the results of this work, we now outline three possible observables for the effects of radial migration.

• (i)

  Old disc flaring: the effect of migration, driven by internaldisc evolution, on the disc heating is mostly seen in the outerdiscs. The reason for this is the following. As we found inSect. 4, as stars migrate outwards their veloc-ity dispersion decreases as to approximately match the energyof stars at the final radius (which did not migrate). However, asmall amount of energy may be retained, with the effect beingstronger, the farther the stars come from. In other words, whilestars do arrive with slightly larger energy at an outer radius, theeffect is negligible, unless the change in guiding radius is verylarge. It thus appears that migration would always work as toincrease/decrease the scale-height in the outer/inner discs,giving rise to flaring. We showed in Sect. 8 thatthe situation may be different when self-consistently growingdiscs are considered (see Fig. 8). In any case, we doexpect that older stellar populations, being exposed to migrationmechanisms the longest, would exhibit flared vertical profiles,even if the total mass in the disc would not flare. Since this effect isrelated to the conservation of the average vertical action of stars, inthe case of an active merger history, this signature may have beenerased.

• (ii)

  Flattening of the average vertical action radial profile: as pointed out by Minchev et al. (2012a), an important consequence of the conservation of actions for migrating stars is that radial migration also mixes the radial distributions of the actions  ⟨ J_p ⟩  and  ⟨ J_z ⟩ : the more the disc mixes, the weaker their variation with radius. This may offer a way to constrain the migration history in the Milky Way. Estimates of the average (better conserved) vertical action in observations for different populations of stars should reveal different variation with radius (flattening) for older groups of stars.

• (iii)

  Structure in vertical chemical gradients: while we here showed that a kinematically thick disc is not expected to result from radial migration, signatures of migration may be found in the current chemical distribution of stars, since outward migrators are preferentially deposited at high altitudes, the reverse being true for inward ones (see Fig. 5). Thus, for a given chemical evolution model for the Galaxy, predictions can be made for vertical chemical gradients as a function of radius as observed today.

Acknowledgments

We would like to thank P. Di Matteo, F. Combes, D. Pfenniger, L. Athanassoula, C. Chiappini, R. de Jong, and M. Steinmetz for helpful discussions and support with this work. We also thank the anonymous referee for valuable suggestions which have greatly improved the manuscript.

References

Appendix A

Fig. A.1 Distributions of migrators and non-migrators ending up in a radial bin in the range 10 < r < 13 kpc, estimated in the time period 0.5 < t < 1.5 Gyr for the gSa model. First row: the radial (left) and vertical (right) velocity distributions for stars which stay in the radial bin during the entire time. An increase of velocity dispersion is seen with time. Second row: stars which came from r < 6 kpc. Note that their velocity dispersions drop by about 40%. Third row: stars originting from r > 13 kpc. These are found to heat up upon arriving at the final radius. In all cases, final distributions (orange histogram) roughly match for each velocity component.

Fig. A.2 Variation of radial bins size with radius estimated as 2σr/κ (black curve in units of 10 km s-1), where σr is the radial velocity dispersion and κ is the epicyclic frequency. The dotted red and solid blue curves show σr and κ, respectively. Note that the bin sizes increase with radius.

Fig. A.3 Contrasting the true density of the gSb model (panels 1 − 3), estimated directly from the initial and final radii of stars in the disc for the times shown, and the density obtained in the method described in Sect. 5 (panels 4 − 6, same as top row of Fig. 4).

A.1. Velocity distributions of migrators and non-migrators

In Fig. 3 we showed that as stellar samples migrate, they decrease/increase their velocity dispersions both in the radial and vertical directions when coming from the inner/outer disc, respectively, so as to approximately match the kinematics of non-migrating stars in the final radial bin. Figure A.1 offers a different view of this process. We consider the gSa model in the time interval 0.5 < t < 1.5 Gyr (same as the third column of Fig. 3). The first row of Fig. A.1 shows the radial (left) and vertical (right) velocity distributions for stars which stay in the radial bin during the entire time. An increase of velocity dispersion is seen with time. The second row shows stars which came from r < 6 kpc. Note that their velocity dispersions drop by about 40%. The third row plots stars originating from r > 13 kpc. These are found to heat up upon arriving at the final radius. In all cases, final distributions (orange histogram) roughly match for each velocity component.

Fig. A.4 The effect of different bin overlapping used in the selection procedure. The top, right pair of panels is the same as the top row of Fig. 6 (bin overlap of 0.3 kpc), showing the fractional changes in the vertical and radial velocity dispersions due to migrators, Δσz(r) and Δσr(r), respectively. We contrast four possibilities: bin overlap of 0.1, 0.3, 0.5 and 1 kpc. The overall shapes of the functions plotted remain the same, while, naturally, structure is lost with an increase in the distance between overlapping bins, i.e., for 0.1 to 1 kpc.

Fig. A.5 Changes in Fig. 6 due to the use of different bin sizes in the selection procedure. In addition to our standard choice of 2a (middle two columns here), in the first two rows we also consider a bin size of a (columns 1 − 2) and 3a (columns 5 − 6) for the gSb and gSa models and earlier and later times, as indicated. The bottom row shows the effect of a constant bin size of 1, 2, and 3 kpc for the gSb model. The overall behavior is preserved in all cases.

A.2. Particle density from selection procedure

Here we discuss in more detail our procedure for separating migrators and non-migrators. Figure A.3 contrasts the true density (panels 1−3), estimated directly from the initial and final radii of stars in the disc for the times shown, and the density obtained in the method described in Sect. 5 (panels 4−6, same as top row of Fig. 4). Note that the overall shape of contours remains the same despite the overlap of particles resulting from our selection method. Considering that the contour levels are the same, we see that the true density drops faster with radius. This is due to the fact that the overlap of stars involves more stars coming from the inner disc, due to the asymmetric drift. It is clear that this overlap will not affect velocity measurements since contribution from multiple stars would cancel out. In any case, we remove multiples prior to estimating the velocity dispersions and disc scale-height.

A.3. Radial bin overlapping

To achieve smooth variation with radius when estimating velocity dispersion profiles, etc., for all figures presented in the paper we overlapped bins every 0.3 kpc. Note that a typical radial bin width is about 2−3 kpc, increasing with radius (see Fig. A.2). In Fig. A.5 we examine the effect of changing the distance between the radial bins. The top, right pair of panels is the same as the top row of Fig. 6 (bin overlap of 0.3 kpc), showing the fractional changes in the vertical and radial velocity dispersions due to migrators, Δσ_z(r) and Δσ_r(r), respectively. We contrast four possibilities: bin overlap of 0.1, 0.3, 0.5 and 1 kpc. The overall shapes of the functions plotted remain the same, while, naturally, structure is lost with an increase in the distance between overlapping bins, i.e., for 0.1 to 1 kpc.

A.4. The effect of different bin sizes

To construct the figures in this paper we have estimated the bin size as 2a, where a ~ σ_r/κ is the approximate epicyclic amplitude at a given radius, with σ_r and κ the radial velocity dispersion and epicyclic frequency, respectively. In order to contain non-migrating stars, we allow for particles to oscillate around their guiding radii by  ± a, therefore, we select our bin size to be 2a. Since both σ_r and κ are functions of radius, so is the value of a.

In Fig. A.5 we show how the migrators’ contribution to the disc velocity dispersion (as shown in Fig. 6) changes with bin size. In addition to the choice of 2a used in Fig. 6 (middle two columns here), in the first two rows we also consider a bin size of a (Cols. 1−2) and 3a (Cols. 5−6) for the gSb and gSa models and earlier and later times, as indicated. The bottom row shows the effect of a constant bin size of 1, 2, and 3 kpc for the gSb model.

In all cases, the innermost disc appears to be affected the most by the change in bin size. Concentrating on the earlier times (when most of the velocity dispersion increase takes place), we note that the oscillatory inverse functional behavior of Δσ_z(r) and Δσ_r(r) is seen for all bin sizes.

All Tables

All Figures

Fig. 1 First column: m=2 Fourier amplitudes, A2, as a function of time (color bar) for our model galaxies. Columns 2–3: edge-on view for each disc and bulge (separately) for earlier and later times of their evolution, as indicated in each panel. Columns 4–5: scale-heights, h, at  ~2.5 disc scale-lengths. The dotted red line shows a single exponential fit. Note that in 3–4 Gyr of evolution with strong perturbation from bars and spirals, h only doubles despite the efficient stellar radial migration.

In the text

Fig. 2 First row: time evolution of the radial vertical dispersion profiles, σz(r), for the gSb, gSa, and gS0 models. Different colors and line styles correspond to the times indicated in the first row, rightmost panel. The bars’ corotation (CR) and outer 2:1 Lindblad resonance (OLR) are shown by the dotted and solid vertical lines, respectively, with colors corresponding to the different times. Second row: same but for the radial velocity dispersion, σr(r). Third row: same but for the ratio σz/σr. The slope in σz/σr outside the bar shows that the radial profile flattens more with time than the vertical one, while its overall decrease means that the disc heats more radially than vertically.

In the text

Fig. 3 First row: initial (dashed black), final (solid blue), and net (dash-dotted red) vertical velocity dispersions, σz,i, σz,f, and Δσz, respectively, as functions of the initial radius, r0, for particles ending up in the green bin at the final time indicated in each panel. Two time periods for each of the three models are shown, as indicated. Migrating stellar samples decrease or increase their velocity dispersion is such a way as to match the value at the destination. Second row: same as the first row, but for the radial velocity dispersion. Third row: the corresponding initial radial distributions of stars ending up in the green bin at the final time indicated in each panel. The bar’s CR and OLR, measured at the final times, are shown by the vertical dotted and solid black lines. Most migrators come from the inner disc, with a strong contribution from particles originating near the bar’s CR.

In the text

Fig. 4 First row: contour density plots of the final versus initial radii for all particles in different time spans for the gSb or gSa models, as indicated. The orange dashed line shows the locus of non-migrating, cold particles. Nonuniform deviations from this line indicate that there exist preferential radii where migration efficiency is stronger, however, these are not clearly seen. Second row: non-migrating stars, extracted from the distributions above as described in the text. Third row: the migrators only, obtained by subtracting the distributions in the second from those in the first row. The dotted-red and solid-blue vertical/horizontal lines in all panels indicate the positions of the bar’s CR and OLR, respectively, at the initial/final times. Strong clumps can now be seen across the bar’s CR for example. Fourth row: number density plots of the changes in angular momentum, ΔL, versus the initial angular momentum (or equivalently initial radius), L, divided by the rotational velocity at each radius. The broken red and blue vertical strips indicate the location of the bar’s CR and OLR, where their widths corresponding to the change during each time period. The pink lines connect some corresponding groups of migrated stars. The bar’s CR is active throughout the entire simulations.

In the text

Fig. 5 Migrators’ contribution to the disc heating. First row: radial profiles of the vertical velocity dispersion, σz(r), for non-migrators (dashed black curve), outward going migrators (solid brown curve), and inward migrators (dash-dotted green curve) for the gSb model, estimated in the time periods indicated in the top row of each column. Second row: σz(r) estimated for non-migrators (dashed black curve), all migrators (i.e., inward and outward going, solid green curve), and the total population (dashed-dotted blue curve). Third and Fourth rows: same as the first and second rows, respectively, but for the radial velocity dispersion profile σr(r). An unexpected result is that outward migrators cool the disc inside the bar’s OLR (and inward ones heat it), unlike in the case of σz. Some increase in σz and σr from migrators is seen beyond the OLR, while the inner disc is cooled both radially and vertically. Figure 6 quantifies better the migrators’ contribution to σz and σr.

In the text

Fig. 6 Fractional changes in σz(r) and σr(r) resulting from the migrated stars. The dashed black curve plots the quantity $\mathrm{\Delta }{\mathit{\sigma }}_{\mathit{z}}\mathrm{=}\mathrm{(}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}\mathrm{-}{\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}\mathrm{)}\mathit{/}{\mathit{\sigma }}_{\mathit{z,}\mathrm{all}}$, where σz,all and ${\mathit{\sigma }}_{\mathit{z,}\mathrm{non}\begin{array}{}\end{array}\mathrm{mig}}$ are the vertical velocity dispersions obtained from the total population and the non-migrators, respectively. A similar expression estimates Δσr, shown by the solid orange curve. The dotted-red and solid-blue vertical lines denote the bar’s CR and OLR, corresponding to the final time shown in each time period. Note the inverse correlation between Δσz and Δσr, indicating coupling between the radial and vertical motions of stars.

In the text

Fig. 7 Scale-height radial profiles for migrators, non-migrators and all stars in the time intervals 0.5 < t < 1.5 Gyr (top) and 1.5 < t < 3 Gyr (bottom). The green solid curves indicate the fraction of migrators as functions of radius during each time period. The scale-height for the non-migrating stars is approximately constant with radius for all models (blue dashed curves, top). Due to the inner disc cooling and outer disc heating in the vertical direction from migrators (dotted-red curve, see also Fig. 6), all discs flare (solid black curve).

In the text

Fig. 8 Application to three high-resolution cosmological re-simulations. First and second rows: face-on and edge-on density contours for the C1, C2, and C3 models. Third row: change in angular momentum as a function of disc radius: effective migration is seen for all models. Fourth row: vertical velocity dispersion profiles for non-migrators, outward and inward migrators. Fifth row: same but for the radial velocity dispersion. Bottom row: fractional contribution of migrators to the total vertical and radial velocity dispersions (as Fig. 6). Results are very similar to the simpler models presented earlier.

In the text

Fig. A.1 Distributions of migrators and non-migrators ending up in a radial bin in the range 10 < r < 13 kpc, estimated in the time period 0.5 < t < 1.5 Gyr for the gSa model. First row: the radial (left) and vertical (right) velocity distributions for stars which stay in the radial bin during the entire time. An increase of velocity dispersion is seen with time. Second row: stars which came from r < 6 kpc. Note that their velocity dispersions drop by about 40%. Third row: stars originting from r > 13 kpc. These are found to heat up upon arriving at the final radius. In all cases, final distributions (orange histogram) roughly match for each velocity component.

In the text

	Fig. A.2 Variation of radial bins size with radius estimated as 2σr/κ (black curve in units of 10 km s-1), where σr is the radial velocity dispersion and κ is the epicyclic frequency. The dotted red and solid blue curves show σr and κ, respectively. Note that the bin sizes increase with radius.
In the text

	Fig. A.3 Contrasting the true density of the gSb model (panels 1 − 3), estimated directly from the initial and final radii of stars in the disc for the times shown, and the density obtained in the method described in Sect. 5 (panels 4 − 6, same as top row of Fig. 4).
In the text

Fig. A.4 The effect of different bin overlapping used in the selection procedure. The top, right pair of panels is the same as the top row of Fig. 6 (bin overlap of 0.3 kpc), showing the fractional changes in the vertical and radial velocity dispersions due to migrators, Δσz(r) and Δσr(r), respectively. We contrast four possibilities: bin overlap of 0.1, 0.3, 0.5 and 1 kpc. The overall shapes of the functions plotted remain the same, while, naturally, structure is lost with an increase in the distance between overlapping bins, i.e., for 0.1 to 1 kpc.

In the text

Fig. A.5 Changes in Fig. 6 due to the use of different bin sizes in the selection procedure. In addition to our standard choice of 2a (middle two columns here), in the first two rows we also consider a bin size of a (columns 1 − 2) and 3a (columns 5 − 6) for the gSb and gSa models and earlier and later times, as indicated. The bottom row shows the effect of a constant bin size of 1, 2, and 3 kpc for the gSb model. The overall behavior is preserved in all cases.

In the text