Free Access
Issue
A&A
Volume 534, October 2011
Article Number A75
Number of page(s) 10
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/201117566
Published online 05 October 2011

© ESO, 2011

1. Introduction

Disk galaxies are highly nonlinear systems which are driven by external forcing (satellites, interaction with nearby galaxies), internal instabilities (related to the formation of internal structures such as spiral arms, rings, bars) or a combination of both. Chaos and complexity, which are two different aspects of nonlinear response, dominate the dynamics of galactic systems. Disk galaxies are “complex” in the sense that they are made of many components (stars, gas, dark matter in the bar-bulge, disk and halo components) whose interactions can give rise to spontaneous self-organization and to the emergence of coherent, collective phenomena. Examples of emergent behavior in disk galaxies are the formation of a central bar or the onset of spiral arms and warps, which can develop in dynamically cold disks (Revaz & Pfenniger 2004). Disk galaxies are not only complex systems, but even display chaotic behavior, since tiny differences in initial orbits of stars can exponentially blow up, as indicated by numerical simulations. Since the chaotic orbits are more sensitive to perturbations than regular (periodic) ones, external/internal forcing are more effective on this chaotic component. In this paper, we investigate chaotic and complex phenomena related to the formation of a central bar which give rise to the diffusion of the stellar component in the disk.

The bar has a time-dependent activity, with a pattern speed which typically decreases in isolated galaxies (Sellwood 1981). However, the system can be cooled or heated by energy dissipation or infall of gas, or by forming stars on low-velocity dispersion orbits, with the net effect of impacting the amplitude of spiral waves and the strength of the bar, or even destroying it. In this way bars (and spiral waves) can be seen as recurrent patterns which can be rebuilt during their long history until the present configuration at redshift z = 0 (Bournaud & Combes 2002). Under the action of these non-axisymmetric patterns, stars move in the disk and gradually increase their velocity dispersion, as suggested by observations in the solar neighborhood (see Holmberg et al. 2009, and references therein) and in external galaxies (Gerssen et al. 2000; Shapiro et al. 2003). The origin and the amount of disk heating are still open to debate.

First attempts to explain such a heating process in disk galaxies were made by empirically modeling the observed increase of the stellar velocity dispersion with age in the solar neighborhood. Wielen (1977) suggested a diffusion mechanism in velocity space, which gives rise to typical relaxation times for young disk stars of the order of the period of revolution and to a deviation of stellar positions of 1.5 kpc in 200 Myr. The result was obtained without making detailed assumptions on the underlying local acceleration process responsible for the diffusion of stellar orbits. Global acceleration processes, such as the gravitational field of stationary density waves or of central bars with constant pattern speed, were ruled out since their contribution to the velocity dispersion of old stars was found to be negligible and concentrated in particular resonance regions (Wielen 1977; Binney & Tremaine 2008, p. 693). In isolated galaxies, different local accelerating mechanisms have been investigated, such as the gravitational encounters between stars and giant molecular clouds (Spitzer & Schwarschild 1951; 1953; Lacey 1984), secular heating produced by transient spiral arms (Barbanis & Woltjer 1967; Carlberg & Sellwood 1985; Fuchs 2001) or the combination of the two processes (Binney & Lacey 1988; Jenkins & Binney 1990). Another heating mechanism was suggested by Minchev & Quillen (2006), who showed that the stellar velocity dispersion can increase with time due to the non-linear coupling between two spiral density waves.

Such local acceleration mechanisms suggest the existence of a significant component of the galactic gravitational field with a rather chaotic behavior. Pfenniger (1986) investigated the relation between diffusion and chaotic orbits. The latter typically react promptly to small perturbations. He pointed out that the effect of perturbations on regular orbits, such as epicyclic orbits, underestimates strongly the stellar diffusion rate when a stellar system becomes non-integrable, as for example in the presence of a central bar. As the central bar develops, reaches its maximal amplitude and then settles down to an almost steady state, its gravitational potential changes in time. In time-dependent potentials, the number of chaotic orbits typically decreases while the system secularly evolves toward a quasi-steady state through collective effects. The system is then ready again to respond (mainly through the remaining irregular orbits) to external perturbations, such as new infall of gas, and to recurrently restore a strong bar. The bar is thus able to perturb orbits of stars born or passing through its region, which can visit at later times the solar neighborhood (Raboud et al. 1998). Indeed, most of the observational signatures of radial mixing reported in the literature (Grenon 1972; 1999; Castro et al. 1997) point to stars coming from a region next to the bulge/bar intersection, suggesting the bar to be a key player in the radial migration process.

The subject of radial migration of stars was revived when a large scatter in the observed age-metallicity relation (AMR) was reported by Edvardsson et al. (1993), later confirmed by the larger Geneva-Copenhagen Survey sample (Nordström et al. 2004; Holmberg et al. 2007; 2009; see also Casagrande et al. 2011). It must be said that even though the AMR has been extensively studied in the solar vicinity, the results are still controversial due to the large uncertainties in stellar ages (see Pont & Eyer 2004). For instance, using a sample of stars for which it was possible to obtain chromospheric ages, Rocha-Pinto et al. (2000) have reported a much tighter AMR.

The mechanisms driving radial diffusion and heating are still hotly debated, and in many cases the role of the bar is not taken into account. Assuming the whole scatter seen in Edvardsson et al. (1993) data was real, Sellwood & Binney (2002) pointed out that the radial excursion predicted by Wielen (1977) was not sufficient to explain the weakness of the AMR in the solar neighborhood. In order to explain both the large scatter in the AMR and the evidence that even old disk stars today have nearly circular orbits, Sellwood & Binney (2002) suggested a new mechanism based on the resonant scattering of stars under the effect of transient spiral waves. In this process, a star initially on a nearly circular orbit resonates with a rotating wave and changes its angular momentum. If the duration of the peak amplitude of the perturbing potential is less than the period of the “horseshoe” orbits, i.e. orbits of particles trapped at the corotation radius of the spiral wave, the star can escape from the potential well without changing its eccentricity. The net effect of this scattering mechanism is that stars migrate radially without heating the disk. In other words, the overall distribution of angular momentum is preserved, except near the corotation region of the transient spiral wave, where stars can have large changes of their angular momenta. Haywood (2008) estimated upper values for the migration rate from 1.5 to 3.7 kpc/Gyr, which agree with the values in Lèpine et al. (2003) for the radial wandering due to the scattering mechanism assumed by Sellwood & Binney (2002).

Radial diffusion of stars (and gas) could have important implications for the interpretation of key observational constraints for the formation of the Galaxy, such as the AMR, metallicity distributions, or the metallicity gradients, since old, probably more metal rich stars that formed at small galactocentric radii, as well as young metal-poor stars formed at large radii are enabled to appear in solar-neighborhood samples (e.g. Haywood 2008). Due to the lack of detailed information on the processes driving stellar radial migration, models of the Galactic chemical evolution have evaluated past history of the solar neighborhood and the formation and evolution of the abundance gradients assuming that the Galaxy can be divided into concentric wide (~1–2 kpc) cylindrical annuli, which evolve independently (van den Bergh 1962; Schmidt 1963; Pagel 1997; Chiappini et al. 1997; Chiappini et al. 2001). Schönrich & Binney (2009a) explored the consequences of mass exchanges between annuli by taking into account the effect of the resonant scattering of stars described before. This approach appears to be successful to replicate many properties of the thick disk in the solar neighborhood without requiring any merger or tidal event (Schönrich & Binney 2009b). High resolution cosmological simulations (Roškar et al. 2008; Loebman et al. 2010) give support to the view that such scattering mechanism determines a significant migration in the stellar disk. However, the strong mixing driven by bar resonances was not taken into account (see below), casting thus doubts on some of the conclusions in the papers quoted above.

The Milky Way (MW) is a barred galaxy and it is clear that the process above, not accounting for the existence of the bar, is probably just one of the processes at play among others. Indeed, the role of resonant couplings between bars and spirals (Tagger et al. 1987) in the distribution of energy and angular momentum in disk galaxies could also play a major role. Recently, Minchev and collaborators (Minchev & Famaey 2010; Minchev et al. 2011) have further analyzed this mixing mechanism finding that resonances between the bar and the spiral arms can act much more efficiently than transient spiral structures, dramatically reducing the predicted mixing time-scales. Moreover, while for the Sellwood & Binney (2002) mechanism to work short-lived transient spirals are required, in barred galaxies, such as the MW, spirals are most likely coupled with the bar as shown by Sparke & Sellwood (1987), and thus longer lived (Binney & Tremaine 2008; Quillen et al. 2011). As a consequence the radial migration process in the MW could have been different than currently predicted.

In order to include the effect of radial migration in chemical evolution models and to gain a global (chemical and kinematic) understanding of the processes at play in the galactic disks, many dynamical aspects need to be further investigated and in particular the role of the bar, that is the strongest non-axisymmetric component in disk galaxies. In this paper, we present N-body simulations of barred spiral galaxies, and study how disks with different degrees of stability, ranging from marginally stable disks with Safronov-Toomre parameter QT ~ 1 to hot disks with QT > 1, respond to the presence of bar patterns. Our aim is to estimate the time and length scales of stellar diffusion in the radial direction and to relate these quantities to the strength of the bar and to the number of hot particles in the disk, i.e. generally chaotic particles which are susceptible to cross the corotation barrier and to explore all space, being characterized by values of the Jacobi integral H larger than the value at the Lagrangian points, H > H(L1,2) (Sparke & Sellwood 1987; Pfenniger & Friedli 1991). We investigate how these characteristic scales evolve in time, and depend on the activity of the bar. We consider different scenarios of diffusion, and discuss their implications for chemical evolution constraints in our Galaxy.

The paper is organized as follows. In Sect. 2 we describe the simulations and the relevant parameters. In Sect. 3 we solve the diffusion equation in axisymmetric systems, we define the diffusion coefficient, the diffusion time-scale and the diffusion length-scale, and the methods used to estimate these quantities from the simulation results. In Sect. 4 we present our results. In Sect. 5 we discuss the implications for chemical evolution models of the MW and we summarise our findings.

2. N-body simulations

We have run self-consistent N-body simulations starting from a bar-unstable axisymmetric model. We have analyzed initial configurations with disk, bulge and dark halo components which differ on the initial value of the Safronov-Toomre parameter QT = σr   κ/(3.36   G   Σ) (Safronov 1960; Toomre 1964), where σr is the radial velocity dispersion of the disk component, G is the gravitational constant, Σ is the disk surface density, κ is the epicycle frequency defined by κ2 = R2/dR + 4Ω2, where Ω is the circular frequency related to the global gravitational potential Φ(R,z,t) in the disk plane z = 0 by Ω2 = (1/R)   Φ/∂R.

The initial mass distribution in our simulations corresponds to a superposition of a pair of axisymmetric Miyamoto-Nagai disks of mass MB, MD, horizontal scales AB + B, AD + B, and identical scale-height B, (1)The first component represents the bulge (B), while the second the disk (D) (Pfenniger & Friedli 1991). The parameters have been set to AB = 0.07 kpc, AD = 1.5 kpc, B = 0.5 kpc, MB/MD = 3/17. The initial particle positions and velocities are found by a Monte-Carlo draw following the density law corresponding to Eq. (1), truncated to a spheroid of semi-axes R = 30 kpc, z = 10 kpc. The number of particles in the disk-bulge component is N = 4 × 106 and the total mass is MBD = 4.2 × 1010   M.

In order to progressively heat the disk, we have added to this bulge-disk component an oblate pseudo-isothermal halo with the following density distribution (except in models m1 and m2): (2)The number of particles in this halo component is NH = 2 × 106, which is a value in the range suggested in Dubinski et al. (2009) in order to obtain convergent behavior in studies of bar formation and evolution. The length-scales are Rh = 7.5 kpc and zh = 3.5 kpc. The density distribution has been truncated to R = 30 kpc, z = 15 kpc. We set the total mass in the dark halo MH to be four times the total mass in the bulge-disk component MBD, except in the model m3, where MH = 2MBD (see Table 1). The effect of adding the halo component is that the bar becomes progressively smaller and with higher pattern speed, the disk is hotter and less sensitive to the bar perturbations.

We then impose the equilibrium of the first and second moments of velocities by solving Jeans’ equations (see e.g., Binney & Tremaine 2008) with a constant QT. The resulting distribution is relaxed for a couple of rotations until ripples spreading through the disk from the center disappear. We use this as the initial condition for the N-body simulations performed by using the Gadget-2 free source code (Springel et al. 2001; Springel 2005).

The initial Gadget-2 configurations considered in this work differ on the values of the Safronov-Toomre parameter, which ranges from QT ~ 5 at two scale lengths from the center for hot disks to QT ~ 1 for marginally stable disks. These QT values are listed in Table 1, along with the initial values of the radial and vertical velocity dispersions at two scale lengths from the center. We have considered these initial values in order to investigate how the radial diffusion depends on the disk sensitivity to perturbations. Thus, we have considered two extreme cases: in one case the disk is marginally stable and spiral waves develop (models m1 and m2), with the main global effect of heating in the radial direction (the Araki parameter σz/σr ~ 0.5 at two scale lengths from the center at the end of the simulations and it decreases at larger radii, see the middle panel of Fig. 1). In the other case (model m6), the disk is heated and the velocity dispersions increase in both the radial and the vertical direction (the Araki parameter σz/σr ~ 0.8 at large radii at the end of the simulation, see the corresponding curve in the middle panel of Fig. 1). The other models considered have intermediate values of QT and σz/σr between these two extreme cases. The Safronov-Toomre and Araki parameters at the final times are shown in the left and middle panels, respectively, of Fig. 1.

thumbnail Fig. 1

Left: Safronov-Toomre-parameter at the final time. Middle: Araki parameter at the final time. Right: rotation curves.

Open with DEXTER

Table 1

Initial configuration of the N-body simulations.

In the right panel of Fig. 1 the rotation curves for all the N-body models listed in Table 1 are shown at the final time. When the halo component is included (in models m3, m4, m5, m6), the rotation curve is flatter and Vc ~ 220   km   s-1 at two scale lengths from the center, that is in the region between 6 and 10 kpc, depending on the model. The rotation curve of model m3, where the ratio between halo and disk-bulge mass is MH/MBD = 2, has intermediate values between that of models without the halo (m1, m2) and the others.

We classify stellar orbits into three dynamical categories (Sparke & Sellwood 1987; Pfenniger & Friedli 1991). The first two dynamical categories are the bar and the disk orbits with the Jacobi integral H = E − ΩpLz smaller than the value at the Lagrangian points L1,2, H < H(L1,2), where E is the total energy and Lz is the z-component of the angular momentum. The separation of particles in the bar or disk component can be easily done since bar orbits typically have smaller values of Lz and E than disk orbits. The third category includes hot orbits for which H ≥ H(L1,2). The models considered here differ in the number of hot particles after the formation of the bar (the ratio between the number of stars in the hot component and the total number of the visible stars is listed in the last column of Table 1).

The bar pattern speed , where θ is the azimuthal angle of the bar major axis (in the inertial frame) calculated by diagonalising the moment of inertia tensor of the bar particles, ranges from 35   km   s-1   kpc-1 for model m1 to 40   km   s-1   kpc-1 for m6, at final times. These values are comparable to those found by Fux (1997 – see his Fig. 5), while recent estimates of the MW bar pattern speed are of the order of 50–60   km   s-1   kpc-1 (Dehnen 2000; Minchev et al. 2007; 2010). The pattern speed is typically slowly decreasing in time with a rate of a few km s-1 kpc-1 Gyr-1 (Fux 1997; Bournaud & Combes 2002), so that the values at the beginning of the simulations are 60   km   s-1   kpc-1 and 80   km   s-1   kpc-1 for models m1 and m6, respectively. The corresponding corotation radius Rc, obtained by the intersection of ΩP with the circular frequency, ΩP(t) = Ω(Rc,t), increases in time (it typically ranges from Rc = 2 to 5 kpc at final times in our models).

thumbnail Fig. 2

Evolution in time of the bar’s strength for models m1 and m6.

Open with DEXTER

In order to understand the role played by the central bar on the distribution of stars in the disk, we follow the evolution of the strength of the bar in time for each model. If Cm is the amplitude of the mode m in the density distribution, (3)the bar’s strength is defined as the mode C2 when the stars j are restricted to the bar component. We normalise C2 with respect to the number of stars in the bar component, C0. This quantity is shown in Fig. 2 for models m1 and m6. Since we do not include the gas component in our models, we are not able to follow their evolution for times longer than a few Gyrs. After this typical time-scale, the bar amplitude saturates and the systems reach a quasi-steady state. Since the inclusion of the gas component necessarily requires the introduction of other less controlled parameters, such as the cooling rate of star formation, we prefer in this first study to limit the integration time over a couple of Gyrs, which is already enough to study the role of the bar in the radial migration process.

The face-on and edge-on views of the density distribution of models m1 (upper panels) and m6 (bottom pannels) are shown in Fig. 3 at time t ~ 550 Myr, that is just after the maximum strength of the bar (cf. Fig. 2). In the first case, both bar and spiral arms develop since the disk is sufficiently cold, while in the second case the disk is hot and only a bar (with a smaller corotation radius than in the previous case) develops in the central region. In the external regions of models m1 and m6, other patterns can be observed, the dominant being the pattern with m = 1. This mode is related to asymmetric distributions of mass pushed by the system rotation and it can give rise to filaments of stars or ring-like structures on long time-scales.

thumbnail Fig. 3

Density maps at t ~ 550 Myr. Right: face-on views, left: edge-on views. Top panels: model m1, bottom panels: model m6.

Open with DEXTER

3. Diffusion equation in axisymmetric systems

In this paper we model the diffusion along the radial direction in the galactic plane by introducing a distribution function F(R,t) which satisfies the phenomenological spatial diffusion equation (from Fourier’s law) in cylindrical coordinates: (4)where D is the diffusion coefficient (in unit of area per time). This diffusion model has already been more succinctly presented in Brunetti et al. (2010). We suppose that for a limited time all the stochastic processes ongoing in a spiral galaxy can be described by such an equation characterized locally by a single positive diffusion coefficient D, to be empirically measured in the full dynamical simulations for a range of R and t. The parameter D conveys the notion that the diffusion process is a surface increase per unit time. For the sake of simplicity we choose D to be constant and independent of F, which makes Eq. (4) a linear partial differential equation for F. Another choice could have been to take C ≡  as a local parameter, where ρ is the mass density. The unit of C is mass per length per time which conveys then the notion of a mass flux gradient. However this choice would have led to a more difficult empirical determination of C since the equation involves then ρ and its gradient. Thus, the model considered here gives a lower limit on the diffusion happening in a barred galaxy, since other more complex diffusion processes are neglected in the present analysis.

The general solution of Eq. (4) with constant D, which is non-singular at R = 0 is given by: (5)where J0(x) = J0( − x) is the Bessel function of the first kind. The function A can be determined by taking the Hankel transform of F(R,0): (6)By inserting Eq. (6) into (5) and assuming that the particles are initially localized at a certain radius R0 at time t0 = 0, F(R,0) = F0R0δ(R − R0), we obtain: (7)Thus, the diffusion of this distribution can be expressed in terms of Bessel and elementary functions (Gradsteyn & Ryzhik 2007, formula 6.633.2). The time-evolution of an initial set of localized particles reads: (8)where I0(x) is the modified Bessel function of the first kind, which is finite at the origin, I0(0) = 1. In Fig. 4 two initial distributions with D = 1 and centered in R0 = 0.5 and 2 (solid lines) evolve in time, as described by Eq. (8) (dashed and dotted lines, respectively). At large radii, the distributions are essentially Gaussian, while at small radii they are strongly modified from the contributions of particles at the center of the cylinder. Equation (8) describes the distribution of the radial positions of stars in the disk at the initial time ti which diffuse toward position R0 at time t0, such that t0 − ti = Δt ≤ TD, where TD is the diffusion time-scale or, equivalently, the distribution of stars which initially are in R0 at t0 and then diffuse toward R with Δt ≤ TD.

thumbnail Fig. 4

Two distributions with D = 1 and centered in R0 = 0.5 and 2 (solid lines) evolve in time (dashed and dotted lines, respectively), as described by Eq. (8).

Open with DEXTER

When R goes to zero, Eq. (8) reduces to: (9)For large values of the argument the modified Bessel function I0(x) → (2πx) − 1/2   exp(x) and thus F(R,t) reduces to: (10)where is the Gaussian distribution with mean value μ =  ⟨ R ⟩  = R0 and standard deviation .

In order to obtain a simple model of the distribution of the stars in a galactic disk, one can consider to envelop Eq. (10) by an exponential surface density Σ(R) ∝ exp( − R/Rd) (see, for example, Sellwood & Binney 2002), where Rd is the disk scale length, thus obtaining: (11)where C and C′ are normalization constants and is the Gaussian distribution with mean value μ′ =  ⟨ R ⟩  = R0 − σ2/Rd and standard deviation . It is important to remember that the diffusion model used for obtaining Eq. (10) is only valid for times smaller than the diffusion time-scale.

We set R0 = R = 8 kpc and Rd ~ 3 kpc. We consider the distribution of stars in R ~ R at the initial time. The previous expression, Eq. (11), can be used to estimate the relative fraction of stars which remain in a diffusion time-scale within the local volume |R − R| ≤ d. This is given by: (12)The integral can be written as: (13)where . If d ≪ σ, we get x +  ~ x −  and thus the probability of staying in the local volume is nearly zero. If d ~ σ ≪ Rd, we have . If d ~ σ ~ Rd, we have . Thus, the fraction of stars which remain in the local volume in a diffusion time-scale strongly depends on the ratio d/σ and on the value of Rd.

The diffusion coefficient D in Eqs. (4)–(8) can be regarded as an instantaneous coefficient which depends on the position at which particles are initially localised and on the diffusion time, since, as already mentioned, modeling the stellar migration as a diffusion process is valid only for time intervals less than the diffusion time-scale, Δt ≤ TD. As we described in Brunetti et al. (2010), the diffusion time-scale can be estimated from the simulation results and it turns out to be of the same order of the rotation period, TD ~ Trot = 2π/Ω(R,t). The diffusion coefficient is calculated by applying the nonlinear least-square method which minimizes the difference between the numerical results and the general solution of the diffusion equation described by Eq. (8) for times less than TD. At each time and radial position in our N-body simulations, we estimate the instantaneous diffusion coefficient and related quantities, thus obtaining a description of stellar migration along the whole simulation.

thumbnail Fig. 5

Top row: contour maps of the diffusion coefficient D. Middle row: contour maps of the radial dispersion . Bottom row: contour maps of the diffusion velocity . Left: model m1, right: model m6.

Open with DEXTER

4. Results

The equations and the numerical methods described in the previous section and in Brunetti et al. (2010) allow us to calculate the diffusion coefficient D(R,t) which is shown in the contour maps of Fig. 5, top row, for the models m1 (left panel) and m6 (right panel). The others models m2, ..., m5 have intermediate values of the diffusion coefficient between these two extreme cases. The bar’s corotation radius is shown in Fig. 5 as a dashed white line. It can be seen that the diffusion coefficient is not constant in time nor in radius. If the disk is not too hot, D has the largest values D ~ 0.12   kpc2   Myr-1 outside the corotation radius of the bar (which increases in time in our simulations), where the density is strongly perturbed by a m = 2 pattern created by the bar and the transient spiral arms, and in the external regions R > 8 kpc, where the density is modulated by a m = 1 pattern. The stars respond collectively to these modulations and the process of migration corresponds to a diffusion in an axisymmetric system. In hot disks, the diffusion coefficient is large only in the external region R > 10 kpc, where the m = 1 mode appears, with values of the order of D ~ 0.08   kpc2   Myr-1. In this latter case, the disk is not sufficiently cold to respond to the m = 2 perturbation created by the central bar.

When the disk is marginally stable with initial Safronov-Toomre parameter QT ~ 1, the diffusion coefficient has the largest values just outside the corotation region, where two different families of orbits are present, as can be inferred by the total energy and angular momentum values of the particles in this region. We consider three different bins of particles located respectively in R = (1.5 ± 0.5) kpc, R = (3.0 ± 0.5) kpc and R = (8.0 ± 0.5) kpc at t = 2.2 Gyr in the model m1. In the first bin, particles are mainly inside the bar. Their energies E span small negative values and the z-component of the angular momentum is nearly centered in Lz ~ 0 (see the two panels of Fig. 6, blue lines, labeled as “bin 1”). Particles in the second bin centered in R = (3.0 ± 0.5) kpc at t = 2.2 Gyr belong to two different types of orbits: one family can migrate only inside the bar, the other can go outside the bar, in the disk. Large values of the diffusion coefficient D near the corotation region (see the left panel of Fig. 5, top row) are related to this superposition of two families of orbits. The corresponding values of E and Lz are shown in the panels of Fig. 6 (red lines, labeled as “bin 2”). The two peaks in the energy distribution correspond to the two families of orbits: bar particles have large negative energies, while particles which can go into the disk have small negative energies. Intermediate values correspond to the so called hot particles, which as already mentioned before have a Jacobi integral H larger than the value of H at the Lagrangian points L1,2, H ≥ H(L1,2). We will discuss the hot particles later on in this section. It is important to note that the two peaks are increasingly less evident in hotter disks and the number of hot particles decreases as QT increases (see the last column in Table 1). In the third bin, particles belong to the disk component: these disk particles have small negative energies and large values of Lz (black lines, labeled as “bin 3”).

From D we can estimate the radial dispersion σ in a rotation period (which is of the order of the diffusion time-scale), . In the middle row of Fig. 5 we show the contour maps of the radial dispersion for models m1 (left panel) and m6 (right panel). If the disk is sufficiently cold, the radial dispersion is high near the corotation region of the bar, where it can recurrently assume values of the order of σ ~ 6 kpc. This implies that internal stars can recurrently be forced by the activity of the bar to migrate in the external region of the disk. At an intermediate radius, such as R = 6 kpc, the radial dispersion is of the order of σ ~ 3 kpc and it increases in the external less dense regions where the pattern m = 1 dominates. In the external region, the diffusion of the stellar component is related to the presence of both patterns m = 1 and m = 2 and it can be enhanced in regions where the bar’s outer Lindblad resonance overlaps with the spiral arms resonances (Minchev et al. 2011).

The radial migration driven by the bar seems to be efficient in cold disks. According to our results for model m1, the number of stars which stay always in a local volume of 100 pc around R = 8 kpc is low, since d ~ 100   pc ≪ σ ~ 5 kpc (see discussion after Eq. (13)). If the disk is hot, the values of σ are lower than the corresponding values for cold disks (at R = 6 kpc the radial dispersion is of the order of σ ~ 1 kpc for model m6, see Fig. 5, middle row, right panel) and consequently radial migration is less effective than before.

thumbnail Fig. 6

Energy values (left panel) and Lz-values (right panel) of the particles at time 2.2 Gyr within the radial ranges R = (1.5 ± 0.5) kpc (bin 1, blue lines), R = (3.0 ± 0.5) kpc (bin 2, red lines) and R = (8.0 ± 0.5) kpc (bin 3, black lines) for model m1.

Open with DEXTER

Another interesting quantity is the diffusion velocity in a rotation period, defined by . We show the contour maps of this quantity in the bottom panels of Fig. 5 for models m1 (left) and m6 (right). The diffusion velocity can reach values of vD ~ 40 km s-1 near the corotation region in the model m1, while it is always less than vD ~ 20 km s-1 in the model m6.

The fact that the corotation region plays a crucial role in stellar diffusion can be also seen from Fig. 7. Here it is shown how the radial position of stars at some final time, Rnow at tnow = 2.2 Gyr (on the vertical axis), depends on the corresponding radial distribution Rpast at time tpast = 300 Myr (on the horizontal axis). The different colors in the map are related to the ratio of the number of stars at the two times, N(Rpast,tpast)/N(Rnow,tnow). It can be seen that in the case of model m1 (left panel) particles near the corotation radius, Rnow ~ Rc = 4 kpc (which is the value of the corotation radius at tnow = 2.2 Gyr in the considered model) came from inside and outside the corotation region, since they belong to two different families of orbits, as discussed before (see Figs. 6). Particles outside the corotation at the final time were spread over the disk in the past, with Rc < Rpast ≤ 10 kpc, while particles well inside the corotation radius, Rnow < 2 kpc, were confined into the bar also in the past. On the contrary, in the hot model m6 (see the right panel of Fig. 7) particles were essentially located at the same positions in the past, with Rpast ~ Rnow ± ΔR, where ΔR ~ 1 kpc and it slightly increases with Rnow. In this case, no particular activity is observed in the corotation region, Rc ~ 2 kpc.

thumbnail Fig. 7

Relative fraction of stars which at time tnow = 2.2 Gyr are at radial position Rnow (vertical axis) and at time tpast = 300 Myr were at Rpast (horizontal axis): (left) m1, (right) m6.

Open with DEXTER

thumbnail Fig. 8

Evolution history of the radial position of stars which at time tnow = 2.2 Gyr are in Rnow = (8.0 ± 0.1) kpc: (left) m1, (middle) m6, (right) m1 without hot particles.

Open with DEXTER

As a further example, let us consider a bin of stars located in Rnow = (8.0 ± 0.1) kpc at time tnow = 2.2 Gyr. Their evolution history is shown in Fig. 8, for the two models, m1 (left panel) and m6 (middle panel). When the bar’s strength is maximum (at t ~ 350 Myr, cf. Fig. 2), stars localized near the corotation radius are forced to move toward larger radii, since there the diffusion coefficient is large (cf. Fig. 5, top row). Most of the stars bounce then back and forth in the region Rc < R < 10 kpc, to reach the final bin position. The evolution history in heated disks such as that of model m6 (see Fig. 8, middle panel) is very different, since stars are always localized in the region R = Rnow ± σ, with Rnow = 8 kpc and σ ~ 2 kpc, in agreement with the corresponding value of the radial dispersion which can be inferred from the right panel of Fig. 5 (middle row).

thumbnail Fig. 9

Histograms of radial position of stars in the bins R = (8 ± 1) kpc (left panel) and R = (3 ± 1) kpc (right panel) at the final time (black line) and 2.2 Gyr before (red for model m1 and blue for model m6).

Open with DEXTER

Hot particles, characterized by a Jacobi integral H > H(L1,2), have a distribution which is maximum in the corotation region. In order to investigate the role of such hot particles in the diffusion process, we have removed them from the bin of stars localized in Rnow = (8.0 ± 0.1) kpc at time tnow = 2.2 Gyr which we have considered just before. In the right panel of Fig. 8, the evolution history of the bin where the hot particles have been eliminated is shown. It must be compared with the evolution of the bin which includes them, shown in the left panel of the same figure. It can be seen that the main difference between the two panels is in the relative number of stars which are able to reach the corotation region. The hot particles migrate radially much more than the other particles in the disk. The number of hot particles in cold disks (such as m1) is nearly three times larger than the number in hot disks (see last column in Table 1).

5. Discussion and conclusions

From the previous sections it is clear that the amount of radial migration in disk galaxies is strongly dependent on the bar and spiral strengths. As we analyzed N-body simulations without gas, in order to focus our study on the pure effect of the bar, we can only trace radial migration for  ~2–3 Gyr. It is out of the scope of the present study to compare our simulations with the observations of the solar neighborhood, which are the result of 10 Gyr of evolution. Indeed, our models were not intended to reproduce the conditions in our Galaxy, thus the precise values of, for example, the diffusion time-scale TD and the radial dispersion σ due to radial migration obtained from these models do not correspond to those valid for the Milky Way. However, the physical processes at play in cold disks, like our model m1, should be similar to those happening in our Galaxy, where the Safronov-Toomre-parameter in the solar vicinity is QT ~ 2. Thus, we expect that the order of magnitude of the relevant quantities obtained from our numerical results are reasonably valid also for the earliest phases of the thin disk of our Galaxy (within  ~2–3 Gyr from the bar formation).

5.1. Implications for chemical evolution models of the Milky Way

Chemical evolution models of our Galaxy traditionally assume that the majority of stars do not migrate over large distances, and model the Galaxy by introducing independent radial annuli which are wide enough (around 1–2 kpc wide) so that this approximation would be a valid one (van den Bergh 1962; Schmidt 1963; Pagel 1997; Chiappini et al. 1997, 2001). The expectation is that intruders from other galactocentric distances would not represent more than a few percent of the stars in the local samples. However, as discussed in Sect. 1, there are recent claims that radial migration was more efficient than previously assumed. This, in turn, is driven by the large scatter in the AMR of the Geneva-Copenhagen sample.

In the particular model of Chiappini et al. (2001), each annulus is 2 kpc wide (i.e. d = 2 kpc). Thus, in this case in the solar vicinity d is nearly half the value of the radial dispersion obtained in model m1 at intermediate radii 6–8 kpc, d ~ 2σ. Thus, the relation between the different length-scales which is approximately valid in the solar vicinity is d ~ σ ~ Rd (see discussion after Eq. (13)) and the percentage of stars which stayed in a volume |R − R| ≤ d = 2 kpc turns out to be of the order of 50% in a diffusion time-scale, TD ~ Trot = 2πR/Vc ~ 223 Myr near the Sun. The region from which the rest of stars comes from depends on the activity of the bar which is not a constant pattern, as assumed in the past (Wielen 1977). If we consider for example the recurrent bar scenario described in Bournaud & Combes (2002), the bar can be rebuilt several times in a Hubble time in galaxies with significant gas accretion. When the bar strength is high, such as at the beginning of our simulations and in general when the bar is rebuilt by episodes of dissipative infall of gas, stars come mainly from the corotation region (see Fig. 8, left panel, t < 500 Myr), which is closer to the center since the corotation radius typically becomes smaller after each reformation episode (Bournaud & Combes 2002). When the bar strength saturates to a constant value during a quiescent phase, stars can span the region between the corotation radius and  ~10–11 kpc from the galactic center (see Fig. 8, left panel, t > 500 Myr).

In Fig. 9 we show two examples of stellar diffusion. Stars localized in R = (8 ± 1) kpc and R = (3 ± 1) kpc at the end of the simulations are in the black bins, in the left and right panel, respectively. The red and blue distributions correspond to the radial positions of the stars 2.2 Gyr before, for model m1 and m6, respectively. As can be seen from the left panel in Fig. 9, when the disk is sufficiently cold (model m1, red distribution), the radial dispersion is higher than the corresponding values for the hot disk (model m6, blue distribution). In 2.2 Gyr (which is much larger than the diffusion time-scale), only 25% of the stars remain in the black bin, all the others come from outside. When the disk is hot, the percentage is nearly twice. In the right panel we show the case of stars near the corotation region (remember that the corotation radius is Rc = 2 kpc for model m1 and 1 kpc for model m6). Stars come mainly from the corotation radius when the disk is cold (red distribution), while more than 60% stay at the same position if the disk is hot.

We can conclude that the dynamical effects of stellar migration should be included in chemical evolution models of our Galaxy insofar as the distance d between each annulus is smaller than the radial dispersion σ, which is related to the diffusion coefficient D and to the diffusion time-scale TD by . The radial dispersion σ depends on the degree of marginality of the disk.

5.2. Conclusions

Disk galaxies are complex systems where collective phenomena give rise to the emergence of nonlinear structures, such as central bars and spiral arms. We have investigated the role of bars in marginally stable disks and overheated disks, and their forcing effect on the hot (chaotic) particle component, which is more sensitive to external/internal perturbations.

Modeling the migration of stars in marginally stable disks as a diffusion process in the radial direction is a powerful tool which allows us to estimate quantitatively two crucial parameters, the diffusion coefficient and the diffusion time-scale. With these quantities we are able to compute two other fundamental quantities which are the radial dispersion and the diffusion velocity at different radii and at different times. It is important to note that such a diffusion model makes sense only if there is a stochastic microscopic component at the origin of the diffusion. Ideally, the diffusion time-scale should not be much shorter than the microscopic e-folding time due to chaos, which in N-body systems is typically of the order of the dynamical time (Miller 1964). Over longer time-scales, the diffusion model becomes influenced by the mere global dynamical evolution of the disk, so its behavior may depart from a simple linear diffusion equation with constant coefficient. Thus, the present diffusion calculation is useful for diffusion time-scales in a range around the rotational period. At each time and radial position in N-body simulations, we are able to estimate the instantaneous diffusion coefficient and the related quantities (i.e., diffusion time-scale, radial dispersion and diffusion velocity), thus obtaining a description of the stellar diffusion on the whole simulation.

We have found that the diffusion time-scale is of the order of one rotation period and that the diffusion coefficient D depends on the evolution history of the disk and on the radial position. Larger values D are found in cold disks near the corotation region, which evolves in time, and in the external region, where asymmetric patterns develop. Marginally stable disks, with QT ~ 1, have two different families of bar orbits with different values of angular momentum Lz and energy E, which determine a large diffusion in the corotation region. In hot disks, QT > 1, stellar diffusion is much more reduced than in the case of marginal disks.

The calculations of both the diffusion coefficient and the diffusion time-scale give us a quantitative measure of the migration process in the disk. Another advantage of studying the diffusion of stars in real space, rather than in velocity space, is that it can be more easily related to the evolution of chemical elements, which can be modeled as tracers which follow the evolution of the stellar component. The diffusion process of stars and tracers can be directly implemented in chemical evolution codes, which we plan to do in the near future.

It is interesting to compare our results with those obtained in a recent analysis of Shevchenko (2011), where the Lyapunov and the diffusion times are estimated for the Quillen’s model (Quillen 2003) which describes the Hamiltonian motion in the solar neighborhood due to the interactions of bar and spiral arms resonances. He found that the Lyapunov time, of the order of 10 Galactic years, depends weakly on the model parameters, which can radically change the extent of the chaotic domain (as we obtain by varying the Safronov-Toomre parameter QT). The diffusion time, which characterizes the transport in the chaotic domain of the phase space, is calculated as the inverse of the diffusion rate in the energy variable (thus differing from our definition), with upper bounds of the order of 10 Gyr. It strongly depends on the radial position in the Galaxy, in agreement with our findings.

We call attention to the fact that although there are good reasons (both theoretical and observational) to expect that radial process took place during the evolution of our Galaxy, it could have been much weaker than what has been proposed so far, as implied by the existence of radial abundance gradients, and its mild time-variation. Unfortunately, the uncertainties in the observed abundance gradient evolution are still large and we need to wait until better ages and distances will be available, which will happen in the near future thanks to GAIA and asteroseismology.

Meanwhile, the theoretical work should focus on the relative importance of the main drivers of radial migration, and in the case of the Milky Way, on the role of the bar. Here we have shown that the radial migration process is not only time-dependent but also changes with galactocentric distance, in connection to the bar, plus spiral arms. We plan to analyze simulations with bar but including also gas accretion, where the radial migration process can be traced for several Gyrs. In this way we will be able to answer if radial migration could have had repeated peaks along the MW history or if it faded away after the first 2–3 Gyr. These models coupled with the chemical information will be essential to interpret the radial mixing effects on the local age-metallicity relation, the metallicity distributions of the local thick and thin disks, and on the evolution of the abundance gradients in both disks.

Acknowledgments

Simulations have been run on the REGOR cluster at Geneva Observatory. We thank Michel Grenon and Ralph Schönrich for useful discussions. This work has been supported by the Swiss National Science Foundation.

References

All Tables

Table 1

Initial configuration of the N-body simulations.

All Figures

thumbnail Fig. 1

Left: Safronov-Toomre-parameter at the final time. Middle: Araki parameter at the final time. Right: rotation curves.

Open with DEXTER
In the text
thumbnail Fig. 2

Evolution in time of the bar’s strength for models m1 and m6.

Open with DEXTER
In the text
thumbnail Fig. 3

Density maps at t ~ 550 Myr. Right: face-on views, left: edge-on views. Top panels: model m1, bottom panels: model m6.

Open with DEXTER
In the text
thumbnail Fig. 4

Two distributions with D = 1 and centered in R0 = 0.5 and 2 (solid lines) evolve in time (dashed and dotted lines, respectively), as described by Eq. (8).

Open with DEXTER
In the text
thumbnail Fig. 5

Top row: contour maps of the diffusion coefficient D. Middle row: contour maps of the radial dispersion . Bottom row: contour maps of the diffusion velocity . Left: model m1, right: model m6.

Open with DEXTER
In the text
thumbnail Fig. 6

Energy values (left panel) and Lz-values (right panel) of the particles at time 2.2 Gyr within the radial ranges R = (1.5 ± 0.5) kpc (bin 1, blue lines), R = (3.0 ± 0.5) kpc (bin 2, red lines) and R = (8.0 ± 0.5) kpc (bin 3, black lines) for model m1.

Open with DEXTER
In the text
thumbnail Fig. 7

Relative fraction of stars which at time tnow = 2.2 Gyr are at radial position Rnow (vertical axis) and at time tpast = 300 Myr were at Rpast (horizontal axis): (left) m1, (right) m6.

Open with DEXTER
In the text
thumbnail Fig. 8

Evolution history of the radial position of stars which at time tnow = 2.2 Gyr are in Rnow = (8.0 ± 0.1) kpc: (left) m1, (middle) m6, (right) m1 without hot particles.

Open with DEXTER
In the text
thumbnail Fig. 9

Histograms of radial position of stars in the bins R = (8 ± 1) kpc (left panel) and R = (3 ± 1) kpc (right panel) at the final time (black line) and 2.2 Gyr before (red for model m1 and blue for model m6).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.